lectures/lecture-3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "code_folding": [
     1
    ],
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 3: Matvecs and matmuls, memory hierarchy, Strassen algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Recap of the previous lecture\n",
    "\n",
    "- Most important matrix norms: spectral and Frobenius\n",
    "\n",
    "- Unitary matrices preserve these norms\n",
    "\n",
    "- There are two \"basic\" classes of unitary matrices: Householder and Givens matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Examples of peak performance\n",
    "\n",
    "**Flops** –– floating point operations per second.\n",
    "\n",
    "Giga = $2^{30} \\approx 10^9$,  \n",
    "Tera = $2^{40} \\approx 10^{12}$,      \n",
    "Peta = $2^{50} \\approx 10^{15}$,       \n",
    "Exa = $2^{60} \\approx 10^{18}$ \n",
    "\n",
    "What is the **peak perfomance** of:\n",
    "\n",
    "1. Modern CPU (Intel Core i7)?\n",
    "2. Modern GPU (NVidia Titan)?\n",
    "3. Largest supercomputer of the world? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "1. Modern CPU (Intel Core i7)   –– 400 Gflops\n",
    "2. Modern GPU (NVidia Titan X) –– 12 Tflops single precision \n",
    "4. [Largest supercomputer in the world](https://www.top500.org/lists/2018/06/) –– 122.3 Pflops \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Matrix-by-vector multiplication (matvec)\n",
    "\n",
    "Multiplication of an $n\\times n$ matrix $A$ by a vector $x$ of size $n\\times 1$ ($y=Ax$):\n",
    "\n",
    "$$\n",
    "y_{i} = \\sum_{i=1}^n a_{ij} x_j\n",
    "$$\n",
    "\n",
    "requires $n^2$ mutliplications and $n(n-1)$ additions. Thus, the overall complexity is $2n^2 - n =$ <font color='red'> $\\mathcal{O}(n^2)$ </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## How bad is $\\mathcal{O}(n^2)$?\n",
    "\n",
    "- Let $A$ be the matrix of pairwise gravitational interaction between planets in a galaxy.\n",
    "\n",
    "- The number of planets in an average galaxy is $10^{11}$, so the size of this matrix is $10^{11} \\times 10^{11}$.\n",
    "\n",
    "- To model evolution in time we have to multiply this matrix by vector at each time step.\n",
    "\n",
    "- Top supercomputers do around $10^{16}$ floating point operations per second (flops), so the time required to multiply the matrix $A$ by a vector is approximately\n",
    "\n",
    "$$\n",
    "\\frac{(10^{11})^2 \\text{ operations}}{10^{16} \\text{ flops}} = 10^6 \\text{ sec} \\approx 11.5 \\text{ days} \n",
    "$$\n",
    "\n",
    "for one time step. If we could multiply it with $\\mathcal{O}(n)$ complexity, we would get\n",
    "\n",
    "$$\n",
    "\\frac{10^{11} \\text{ operations}}{10^{16} \\text{ flops}} = 10^{-5} \\text{ sec}.\n",
    "$$\n",
    "\n",
    "Here is the YouTube video that illustrates collision of two galaxisies which was modelled by $\\mathcal{O}(n \\log n)$ algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from IPython.display import YouTubeVideo\n",
    "YouTubeVideo(\"7HF5Oy8IMoM\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Can we beat $\\mathcal{O}(n^2)$?\n",
    "\n",
    "- Generally speaking **NO**. \n",
    "- The point is that we have $\\mathcal{O}(n^2)$ input data, so there is no way to be faster for a general matrix.\n",
    "- Fortunately, we can be faster <font color='red'>for certain types of matrices</font>.\n",
    "Here are some examples:\n",
    "\n",
    "    * The simplest example may be a matrix of all ones, which can be easily multiplied with only $n-1$ additions. This matrix is of rank one. More generally we can multiply fast by <font color='red'>low-rank </font> matrices (or by matrices that have low-rank blocks)\n",
    "\n",
    "    * <font color='red'>Sparse</font> matrices (contain $\\mathcal{O}(n)$ nonzero elements)\n",
    "\n",
    "    * <font color='red'>Structured</font> matrices:\n",
    "      * Fourier\n",
    "      * Circulant\n",
    "      * Toeplitz \n",
    "      * Hankel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Matrix-by-matrix product\n",
    "\n",
    "Consider composition of two linear operators:\n",
    "\n",
    "1. $y = Bx$\n",
    "2. $z = Ay$\n",
    "\n",
    "Then, $z = Ay =  A B x = C x$, where $C$ is the **matrix-by-matrix product**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Matrix-by-matrix product (MM): classics\n",
    "\n",
    "**Definition**. A product of an $n \\times k$ matrix $A$ and a $k \\times m$ matrix $B$ is a $n \\times m$ matrix $C$ with the elements  \n",
    "$$\n",
    "   c_{ij} = \\sum_{s=1}^k a_{is} b_{sj}, \\quad i = 1, \\ldots, n, \\quad j = 1, \\ldots, m \n",
    "$$\n",
    "\n",
    "For $m=k=n$ complexity of a na&iuml;ve algorithm is $2n^3 - n^2 =$ <font color='red'>$\\mathcal{O}(n^3)$</font>.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Discussion of MM\n",
    "\n",
    "- Matrix-by-matrix product is the **core** for almost all efficient algorithms in linear algebra.  \n",
    "\n",
    "- Basically, all the dense NLA algorithms are reduced to a sequence of matrix-by-matrix products.\n",
    "\n",
    "- Efficient implementation of MM reduces the complexity of numerical algorithms by the same factor.  \n",
    "\n",
    "- However, implementing MM is not easy at all!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Efficient implementation for MM\n",
    "**Q1**: Is it easy to multiply a matrix by a matrix in the most efficient way?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Answer: no, it is not easy\n",
    "\n",
    "if you want it as fast as possible, using the computers that are at hand."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Demo\n",
    "Let us do a short demo and compare a `np.dot()` procedure which in my case uses MKL with a hand-written matrix-by-matrix routine in Python and also its numba version."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [],
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "def matmul(a, b):\n",
    "    n = a.shape[0]\n",
    "    k = a.shape[1]\n",
    "    m = b.shape[1]  \n",
    "    c = np.zeros((n, m))\n",
    "    for i in range(n):\n",
    "        for j in range(m):\n",
    "            for s in range(k):\n",
    "                c[i, j] += a[i, s] * b[s, j]\n",
    "                \n",
    "    return c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numba import jit # Just-in-time compiler for Python, see http://numba.pydata.org \n",
    "\n",
    "@jit(nopython=True)\n",
    "def numba_matmul(a, b):\n",
    "    n = a.shape[0]\n",
    "    k = a.shape[1]\n",
    "    m = b.shape[1]\n",
    "    c = np.zeros((n, m))\n",
    "    for i in range(n):\n",
    "        for j in range(m):\n",
    "            for s in range(k):\n",
    "                c[i, j] += a[i, s] * b[s, j]\n",
    "    return c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Then we just compare computational times.\n",
    "\n",
    "Guess the answer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "n = 100\n",
    "a = np.random.randn(n, n)\n",
    "b = np.random.randn(n, n)\n",
    "\n",
    "%timeit matmul(a, b)\n",
    "%timeit numba_matmul(a, b)\n",
    "%timeit np.dot(a, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Why is na&iuml;ve implementation slow?\n",
    "It is slow due to two issues:\n",
    "\n",
    "- It does not use the benefits of fast memory (cache) and in general memory architecture\n",
    "- It does not use available parallelization ability (especially important for GPU) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Memory architecture\n",
    "<img width=80% src=\"pics/Memory-Hierarchy.jpg\">\n",
    "\n",
    "- Fast memory is small\n",
    "- Bigger memory is slow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Making algorithms more computationally intensive \n",
    "\n",
    "<font color='red'>**Implementation in NLA**</font>: use block version of algorithms. <br>\n",
    "\n",
    "This approach is a core of **[BLAS (Basic Linear Algebra Subroutines)](http://www.netlib.org/blas/)**, written in Fortran many years ago, and still rules the computational world."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Split the matrix into blocks! For illustration consider splitting in $2 \\times 2$ block matrix:\n",
    "\n",
    "$$\n",
    "   A = \\begin{bmatrix}\n",
    "         A_{11} & A_{12} \\\\\n",
    "         A_{21} & A_{22}\n",
    "        \\end{bmatrix}, \\quad B = \\begin{bmatrix}\n",
    "         B_{11} & B_{12} \\\\\n",
    "         B_{21} & B_{22}\n",
    "        \\end{bmatrix}$$\n",
    "\n",
    "Then,  \n",
    "\n",
    "$$AB = \\begin{bmatrix}A_{11} B_{11} + A_{12} B_{21} & A_{11} B_{12} + A_{12} B_{22} \\\\\n",
    "            A_{21} B_{11} + A_{22} B_{21} & A_{21} B_{12} + A_{22} B_{22}\\end{bmatrix}.$$  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "If $A_{11}, B_{11}$ and their product fit into the cache memory (which is 1024 Kb for the [recent Intel Chip](http://en.wikipedia.org/wiki/List_of_Intel_Core_i7_microprocessors#.22Haswell-H.22_.28MCP.2C_quad-core.2C_22_nm.29)), then we load them only once into the memory.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## BLAS\n",
    "BLAS has three levels:\n",
    "1. BLAS-1, operations like $c = a + b$\n",
    "2. BLAS-2, operations like matrix-by-vector product\n",
    "3. BLAS-3, matrix-by-matrix product\n",
    "\n",
    "What is the principal differences between them?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The main difference is the number of operations vs. the number of input data!\n",
    "\n",
    "1. BLAS-1: $\\mathcal{O}(n)$ data, $\\mathcal{O}(n)$ operations\n",
    "2. BLAS-2: $\\mathcal{O}(n^2)$ data, $\\mathcal{O}(n^2)$ operations\n",
    "3. BLAS-3: $\\mathcal{O}(n^2)$ data, $\\mathcal{O}(n^3)$ operations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Why BLAS is so important and actual?\n",
    "\n",
    "1. The state-of-the-art implementation of the basic linear algebra operations\n",
    "2. Provides standard names for operations in any new implementations (e.g. [ATLAS](https://www.netlib.org/atlas/), [OpenBLAS](https://www.openblas.net/), [MKL](https://software.intel.com/en-us/mkl)). You can call matrix-by-matrix multiplication function (GEMM), link your code with any BLAS implementation and it will work correctly\n",
    "3. Formulate new algorithms in terms of BLAS operations\n",
    "4. There are wrappers for the most popular languages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Packages related to BLAS\n",
    "\n",
    "1. [ATLAS](http://math-atlas.sourceforge.net) - Automatic Tuned Linear Algebra Software. It automatically adapts to a particular system architechture. \n",
    "2. [LAPACK](http://www.netlib.org/lapack/) - Linear Algebra Package. It provides high-level linear algebra operations (e.g. matrix factorizations), which are based on calls of BLAS subroutines.\n",
    "3. [Intel MKL](https://software.intel.com/en-us/intel-mkl) - Math Kernel Library. It provides re-implementation of BLAS and LAPACK, optimized for Intel processors. Available in Anaconda Python distribution: \n",
    "```\n",
    "conda install mkl\n",
    "```\n",
    "MATLAB uses Intel MKL by default.\n",
    "\n",
    "4. OpenBLAS is an optimized BLAS library based on [GotoBLAS](https://en.wikipedia.org/wiki/GotoBLAS). \n",
    "\n",
    "For comparison of OpenBLAS and Intel MKL, see this [review](https://software.intel.com/en-us/articles/performance-comparison-of-openblas-and-intel-math-kernel-library-in-r)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Faster algorithms for matrix multiplication\n",
    "\n",
    "Recall that matrix-matrix multiplication costs $\\mathcal{O}(n^3)$ operations.\n",
    "However, storage is $\\mathcal{O}(n^2)$.\n",
    "\n",
    "**Question:** is it possible to reduce number operations down to $\\mathcal{O}(n^2)$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Answer**: a quest for $\\mathcal{O}(n^2)$ matrix-by-matrix multiplication algorithm is not yet done."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "* Strassen gives $\\mathcal{O}(n^{2.807\\dots})$ –– sometimes used in practice \n",
    "\n",
    "* [Current world record](http://arxiv.org/pdf/1401.7714v1.pdf) $\\mathcal{O}(n^{2.37\\dots})$ –– big constant, not practical. Based on [Coppersmith-Winograd_algorithm](https://en.wikipedia.org/wiki/Coppersmith%E2%80%93Winograd_algorithm).\n",
    "\n",
    "Consider Strassen in more details. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Na&iuml;ve multiplication\n",
    "\n",
    "Let $A$ and $B$ be two $2\\times 2$ matrices. Na&iuml;ve multiplication $C = AB$\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix} c_{11} & c_{12} \\\\ c_{21} & c_{22}  \\end{bmatrix}  =\n",
    "\\begin{bmatrix} a_{11} & a_{12} \\\\ a_{21} & a_{22}  \\end{bmatrix}\n",
    "\\begin{bmatrix} b_{11} & b_{12} \\\\ b_{21} & b_{22}  \\end{bmatrix} =\n",
    "\\begin{bmatrix} \n",
    "a_{11}b_{11} + a_{12}b_{21} & a_{11}b_{21} + a_{12}b_{22} \\\\ \n",
    "a_{21}b_{11} + a_{22}b_{21} & a_{21}b_{21} + a_{22}b_{22} \n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "contains $8$ multiplications and $4$ additions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Strassen algorithm\n",
    "\n",
    "In the work [Gaussian elimination is not optimal](http://link.springer.com/article/10.1007%2FBF02165411?LI=true) (1969) Strassen found that one can calculate $C$ using 18 additions and only 7 multiplications:\n",
    "$$\n",
    "\\begin{split}\n",
    "c_{11} &= f_1 + f_4 - f_5 + f_7, \\\\\n",
    "c_{12} &= f_3 + f_5, \\\\\n",
    "c_{21} &= f_2 + f_4, \\\\\n",
    "c_{22} &= f_1 - f_2 + f_3 + f_6,\n",
    "\\end{split}\n",
    "$$\n",
    "where\n",
    "$$\n",
    "\\begin{split}\n",
    "f_1 &= (a_{11} + a_{22}) (b_{11} + b_{22}), \\\\\n",
    "f_2 &= (a_{21} + a_{22}) b_{11}, \\\\\n",
    "f_3 &= a_{11} (b_{12} - b_{22}), \\\\\n",
    "f_4 &= a_{22} (b_{21} - b_{11}), \\\\\n",
    "f_5 &= (a_{11} + a_{12}) b_{22}, \\\\\n",
    "f_6 &= (a_{21} - a_{11}) (b_{11} + b_{12}), \\\\\n",
    "f_7 &= (a_{12} - a_{22}) (b_{21} + b_{22}).\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "Fortunately, these formulas hold even if $a_{ij}$ and $b_{ij}$, $i,j=1,2$ are block matrices.\n",
    "\n",
    "Thus, Strassen algorithm looks as follows. \n",
    "- First of all we <font color='red'>split</font> matrices $A$ and $B$ of sizes $n\\times n$, $n=2^d$ <font color='red'> into 4 blocks</font> of size $\\frac{n}{2}\\times \\frac{n}{2}$\n",
    "- Then we <font color='red'>calculate multiplications</font> in the described formulas <font color='red'>recursively</font>\n",
    "\n",
    "This leads us again to the **divide and conquer** idea."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Complexity of the Strassen algorithm\n",
    "\n",
    "#### Number of multiplications\n",
    "\n",
    "Calculation of number of multiplications is a trivial task. Let us denote by $M(n)$ number of multiplications used to multiply 2 matrices of sizes $n\\times n$ using the divide and conquer concept.\n",
    "Then for na&iuml;ve algorithm we have number of multiplications\n",
    "\n",
    "$$\n",
    "M_\\text{naive}(n) = 8 M_\\text{naive}\\left(\\frac{n}{2} \\right) = 8^2 M_\\text{naive}\\left(\\frac{n}{4} \\right)\n",
    "= \\dots = 8^{d-1} M(1) = 8^{d} = 8^{\\log_2 n} = n^{\\log_2 8} = n^3\n",
    "$$\n",
    "\n",
    "So, even when using divide and coquer idea we can not be better than $n^3$.\n",
    "\n",
    "Let us calculate number of multiplications for the Strassen algorithm:\n",
    "\n",
    "$$\n",
    "M_\\text{strassen}(n) = 7 M_\\text{strassen}\\left(\\frac{n}{2} \\right) = 7^2 M_\\text{strassen}\\left(\\frac{n}{4} \\right)\n",
    "= \\dots = 7^{d-1} M(1) = 7^{d} = 7^{\\log_2 n} = n^{\\log_2 7}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Number of additions\n",
    "\n",
    "There is no point to estimate number of addtitions $A(n)$ for naive algorithm, as we already got $n^3$ multiplications.  \n",
    "For the Strassen algorithm we have:\n",
    "$$\n",
    "A_\\text{strassen}(n) = 7 A_\\text{strassen}\\left( \\frac{n}{2} \\right) + 18 \\left( \\frac{n}{2} \\right)^2\n",
    "$$\n",
    "since on the first level we have to add $\\frac{n}{2}\\times \\frac{n}{2}$ matrices 18 times and then go deeper for each of the 7 multiplications. Thus,\n",
    "\n",
    "<font size=2.0>\n",
    "$$\n",
    "\\begin{split}\n",
    "A_\\text{strassen}(n) =& 7 A_\\text{strassen}\\left( \\frac{n}{2} \\right) + 18 \\left( \\frac{n}{2} \\right)^2 = 7 \\left(7 A_\\text{strassen}\\left( \\frac{n}{4} \\right) + 18 \\left( \\frac{n}{4} \\right)^2 \\right) + 18 \\left( \\frac{n}{2} \\right)^2 =\n",
    "7^2 A_\\text{strassen}\\left( \\frac{n}{4} \\right) + 7\\cdot 18 \\left( \\frac{n}{4} \\right)^2 +  18 \\left( \\frac{n}{2} \\right)^2 = \\\\\n",
    "=& \\dots = 18 \\sum_{k=1}^d 7^{k-1} \\left( \\frac{n}{2^k} \\right)^2 = \\frac{18}{4} n^2 \\sum_{k=1}^d \\left(\\frac{7}{4} \\right)^{k-1} = \\frac{18}{4} n^2 \\frac{\\left(\\frac{7}{4} \\right)^d - 1}{\\frac{7}{4} - 1} = 6 n^2 \\left( \\left(\\frac{7}{4} \\right)^d - 1\\right) \\leqslant 6 n^2 \\left(\\frac{7}{4} \\right)^d = 6 n^{\\log_2 7}\n",
    "\\end{split}\n",
    "$$\n",
    "</font>\n",
    "\n",
    "(since $4^d = n^2$ and $7^d = n^{\\log_2 7}$).\n",
    "\n",
    "\n",
    "Asymptotic behavior of $A(n)$ could be also found from the [master theorem](https://en.wikipedia.org/wiki/Master_theorem)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Total complexity\n",
    "\n",
    "Total complexity is $M_\\text{strassen}(n) + A_\\text{strassen}(n)=$ <font color='red'>$7 n^{\\log_2 7}$</font>. Strassen algorithm becomes faster\n",
    "when\n",
    "$$\n",
    "\\begin{split}\n",
    "2n^3 &> 7 n^{\\log_2 7}, \\\\\n",
    "n &> 667,\n",
    "\\end{split}\n",
    "$$\n",
    "so it is not a good idea to get to the bottom level of recursion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Strassen algorithm reloaded\n",
    "\n",
    "- Recent paper [Strassen algorithm reloaded](http://jianyuhuang.com/papers/sc16.pdf)\n",
    "claim to **break conventional wisdom** that Strassen algorithm is not very practical."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- Conventional wisdom: it is only\n",
    " practical for very large matrices. The proposed implementation is practical\n",
    " for small matrices. \n",
    "- Conventional wisdom: the matrices being\n",
    " multiplied should be relatively square. The proposed implementation is\n",
    " practical for rank-$k$ updates, where $k$ is relatively small (a shape\n",
    " of importance for libraries like LAPACK). \n",
    "- Conventional wisdom:\n",
    "  it inherently requires substantial workspace. The proposed implementation\n",
    "  requires no workspace beyond buffers already incorporated\n",
    "  into conventional high-performance DGEMM implementations.\n",
    "- Conventional wisdom: a Strassen DGEMM interface must pass\n",
    " in workspace. The proposed implementation requires no such workspace\n",
    " and can be plug-compatible with the standard DGEMM interface.\n",
    "- Conventional wisdom: it is hard to demonstrate speedup\n",
    " on multi-core architectures. The proposed implementation demonstrates\n",
    " speedup over conventional DGEMM even on an IntelR Xeon\n",
    " PhiTM coprocessor utilizing 240 threads. It is shown how a distributed\n",
    " memory matrix-matrix multiplication also benefits from\n",
    " these advances."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Strassen algorithm and tensor rank (advanced topic)\n",
    "\n",
    "- It is not clear how Strassen found these formulas. \n",
    "- However, now we can see that they are not artificial.\n",
    "- There is a general approach based on the so-called tensor decomposition technique. \n",
    "- Here by tensor we imply a multidimensional array - generalization of the matrix concept to many dimensions.\n",
    "\n",
    "Let us enumerate elements in the $2\\times 2$ matrices as follows\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix} c_{1} & c_{3} \\\\ c_{2} & c_{4}  \\end{bmatrix} =\n",
    "\\begin{bmatrix} a_{1} & a_{3} \\\\ a_{2} & a_{4}  \\end{bmatrix}\n",
    "\\begin{bmatrix} b_{1} & b_{3} \\\\ b_{2} & b_{4}  \\end{bmatrix}=\n",
    "\\begin{bmatrix} \n",
    "a_{1}b_{1} + a_{3}b_{2} & a_{1}b_{3} + a_{3}b_{4} \\\\ \n",
    "a_{2}b_{1} + a_{4}b_{2} & a_{2}b_{3} + a_{4}b_{4} \n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "This can be written as\n",
    "$$\n",
    "c_k = \\sum_{i=1}^4 \\sum_{j=1}^4 x_{ijk} a_i b_j, \\quad k=1,2,3,4\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "$x_{ijk}$ is a 3-dimensional array, that consists of zeros and ones:\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "x_{\\ :,\\ :,\\ 1} = \n",
    "\\begin{pmatrix}\n",
    "1 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 1 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "\\end{pmatrix}\n",
    "\\quad\n",
    "x_{\\ :,\\ :,\\ 2} = \n",
    "\\begin{pmatrix}\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "1 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 1 & 0 & 0 \\\\\n",
    "\\end{pmatrix} \\\\\n",
    "x_{\\ :,\\ :,\\ 3} = \n",
    "\\begin{pmatrix}\n",
    "0 & 0 & 1 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 1 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "\\end{pmatrix}\n",
    "\\quad\n",
    "x_{\\ :,\\ :,\\ 4} = \n",
    "\\begin{pmatrix}\n",
    "0 & 0 & 1 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 1 \\\\\n",
    "\\end{pmatrix}\n",
    "\\end{split}\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Trilinear decomposition\n",
    "\n",
    "To get Strassen algorithm we should do the following trick –– decompose $x_{ijk}$ in the following way\n",
    "\n",
    "$$\n",
    "x_{ijk} = \\sum_{\\alpha=1}^r u_{i\\alpha} v_{j\\alpha} w_{k\\alpha}.\n",
    "$$\n",
    "\n",
    "This decomposition is called **trilinear tensor decomposition** and has a meaning of separation of variables: we have a sum of $r$ (called rank) summands with separated $i$, $j$ and $k$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Strassen via trilinear\n",
    "\n",
    "Now we have\n",
    "$$\n",
    "c_k = \\sum_{\\alpha=1}^r w_{k\\alpha} \\left(\\sum_{i=1}^4  u_{i\\alpha} a_i \\right) \\left( \\sum_{j=1}^4 v_{j\\alpha} b_j\\right), \\quad k=1,2,3,4.\n",
    "$$\n",
    "Multiplications by $u_{i\\alpha}$ or $v_{j\\alpha}$ or $w_{k\\alpha}$ do not require recursion since $u, v$ and $w$ are known precomputed matrices. Therefore, we have only $r$ multiplications of $\\left(\\sum_{i=1}^4  u_{i\\alpha} a_i \\right)$ $\\left( \\sum_{j=1}^4 v_{j\\alpha} b_j\\right)$ where both factors depend on the input data.\n",
    " \n",
    "As you might guess array $x_{ijk}$ has rank $r=7$, which leads us to $7$ multiplications and to the Strassen algorithm!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary of MM part\n",
    "- MM is the core of NLA. You have to think in block terms, if you want high efficiency\n",
    "- This is all about computer memory hierarchy\n",
    "- Concept of block algorithms\n",
    "- (Advanced topic) Strassen and trilinear form"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}