lectures/lecture3/lecture3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Лекция 3\n",
    "## Часть 1. Умножение матрицы на вектор и умножение матриц, иерархия памяти и алгоритм Штрассена"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## На прошлой лекции\n",
    "\n",
    "- Ранг матрицы\n",
    "- Малоранговая аппроксимация\n",
    "- SVD и теорема Eckart'a-Young'a\n",
    "- Линейные системы\n",
    "- Ряд Неймана и число обусловленности\n",
    "- Линейная задача наименьших квадратов"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Примеры пиковой производительности\n",
    "\n",
    "**Flops** –– число операций над числами с плавающей точкой в секунду (floating point operations per second).\n",
    "\n",
    "Гига = $2^{30} \\approx 10^9$,  \n",
    "Тера = $2^{40} \\approx 10^{12}$,      \n",
    "Пета = $2^{50} \\approx 10^{15}$,       \n",
    "Экза = $2^{60} \\approx 10^{18}$ \n",
    "\n",
    "Какова **пиковая производительность**:\n",
    "\n",
    "1. Cовременных CPU?\n",
    "2. Cовремененных GPU?\n",
    "3. Cамого большого суперкомпьютера? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Тактовая частота CPU vs. производительность в flops\n",
    "\n",
    "FLOPS = sockets * (cores per socket) * (number of clock cycles per second) * (number of floating point operations per cycle).\n",
    "\n",
    "- Обычно sockets = 1\n",
    "- Число ядер обычно 2 или 4\n",
    "- Число тактов в секунду - это привычная частота\n",
    "- Число операций за один такт зависит от конкретного процессора"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "1. Современных CPU (Intel Core i7)   –– 400 Gflops\n",
    "2. Современных GPU (NVidia Titan X) –– 12 Tflops single precision \n",
    "3. [Самый большой суперкомпьютер](https://www.top500.org/lists/2019/06/) –– 200.79 Pflops "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Умножение матрицы на вектор (matvec)\n",
    "\n",
    "Умножение матрицы $A$ размера $n\\times n$ на вектор $x$ размера $n\\times 1$ ($y=Ax$):\n",
    "\n",
    "$$\n",
    "y_{i} = \\sum_{i=1}^n a_{ij} x_j\n",
    "$$\n",
    "\n",
    "требует $n^2$ умножений и $n(n-1)$ сложений. Таким образом, общая сложность равна $2n^2 - n =$ <font color='red'> $\\mathcal{O}(n^2)$ </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Насколько плоха сложность $\\mathcal{O}(n^2)$?\n",
    "\n",
    "- Пусть $A$ – матрица попарных гравитационных взаимодействий между планетами в галактике.\n",
    "\n",
    "- Число планет в средней галактике $10^{11}$, поэтому размер этой матрицы $10^{11} \\times 10^{11}$.\n",
    "\n",
    "- Для моделирования эволюции во времени нам нужно умножать эту матрицу на вектор для получения состояния в каждый момент времени.\n",
    "\n",
    "- Лучшие суперкомпьютеры делают примерно $10^{16}$ операций в секунду (flops), поэтому для умножения матрицы на вектор нам потребуется примерно\n",
    "\n",
    "$$\n",
    "\\frac{(10^{11})^2 \\text{ операций}}{10^{16} \\text{ flops}} = 10^6 \\text{ сек.} \\approx 11.5 \\text{ дней} \n",
    "$$\n",
    "\n",
    "для одного шага по времени. Если бы мы могли умножать матрицу на вектор за $\\mathcal{O}(n)$ операций, мы бы получили\n",
    "\n",
    "$$\n",
    "\\frac{10^{11} \\text{ операций}}{10^{16} \\text{ flops}} = 10^{-5} \\text{ сек}.\n",
    "$$\n",
    "\n",
    "Далее посмотрим видео, моделирующее столкновение двух галактик, которое было получено с использованием алгоритма умножения матрицы на вектор со сложностью $\\mathcal{O}(n \\log n)$:"
   ]
  },
  {
   "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": [
    "## Можно ли умножить матрицу на вектор быстрее $\\mathcal{O}(n^2)$?\n",
    "\n",
    "- Вообще говоря **НЕТ**. \n",
    "- Это связано с тем, что мы имеем $\\mathcal{O}(n^2)$ входных данных, поэтому для произвольной матрицы быстрее выполнить умножение нельзя.\n",
    "- К счастью, мы можем  умножать матрицу на вектор быстрее <font color='red'>для определённого класса матриц</font>.\n",
    "Например:\n",
    "\n",
    "    * <font color='red'>Малоранговые</font> матрицы (или матрицы с малоранговыми блоками)\n",
    "\n",
    "    * <font color='red'>Разреженные</font> матрицы с $\\mathcal{O}(n)$ ненулевыми элементами\n",
    "\n",
    "    * <font color='red'>Структурированные</font> матрицы:\n",
    "      * Фурье\n",
    "      * Циркулянтные\n",
    "      * Тёплицевы \n",
    "      * Ганкелевы"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Произведение матриц\n",
    "\n",
    "Рассмотрим композицию линейных операторов:\n",
    "\n",
    "1. $y = Bx$\n",
    "2. $z = Ay$\n",
    "\n",
    "Тогда, $z = Ay =  A B x = C x$, где $C$ – **произведение двух матриц**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Классическое определение\n",
    "\n",
    "**Определение**. Произведение матрицы $A$ размера $n \\times k$ и матрицы $B$ размера $k \\times m$ – это матрица $C$ размера $n \\times m$ такая что её элементы записываются как \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",
    "Для $m=k=n$ сложность наивного алгоритма составляет $2n^3 - n^2 =$ <font color='red'>$\\mathcal{O}(n^3)$</font>. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Обсуждение\n",
    "\n",
    "- Умножение матриц – это основная операция, необходимая для почти всех эффективных алгоритмов в вычислительной линейной алгебре.\n",
    "\n",
    "- Говоря нестрого, все алгоритмы для решения задач вычислительной линейной алгебры с плотными матрицами сводятся к последовательности умножений матриц.\n",
    "\n",
    "- Повышение эффективности умножения матриц, пропорционально увеличивает эффективность алгоритмов, которые используют эту операцию (то есть почти всех).\n",
    "\n",
    "- Однако реализация быстрых методов для умножения матриц – непростая задача!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Эффективная реализация алгоритмов вычисления произведения матриц\n",
    "\n",
    "**Вопрос**: Легко ли вычислить произведение матриц наиболее эфективным образом?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Ответ: нет, это не легко!\n",
    "\n",
    "Если вы хотите вычислить произведение как можно быстрее, используя обычный компьютер."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Демо\n",
    "\n",
    "Сравним время вычисления произведения матриц для следующих реализаций\n",
    "\n",
    "- стандартная функция `np.dot()`\n",
    "- рукописная реализация алгоритма \n",
    "- Numba версия рукописной реализации"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "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": 3,
   "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": [
    "Теперь давайте сравним время вычислений.\n",
    "\n",
    "Попробуйте угадать ответ :)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "536 ms ± 16 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "1.02 ms ± 12.7 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "18.4 µs ± 852 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
     ]
    }
   ],
   "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": [
    "Для всех ли размерностей выполнены такие соотношения?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dimension = 10\n",
      "Dimension = 20\n",
      "Dimension = 30\n",
      "Dimension = 40\n",
      "Dimension = 50\n",
      "Dimension = 60\n",
      "Dimension = 70\n",
      "Dimension = 80\n",
      "Dimension = 90\n",
      "Dimension = 100\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "dim_range = [10*i for i in range(1, 11)]\n",
    "time_range_matmul = []\n",
    "time_range_numba_matmul = []\n",
    "time_range_np = []\n",
    "for n in dim_range:\n",
    "    print(\"Dimension = {}\".format(n))\n",
    "    a = np.random.randn(n, n)\n",
    "    b = np.random.randn(n, n)\n",
    "\n",
    "    t = %timeit -o -q matmul(a, b)\n",
    "    time_range_matmul.append(t.average)\n",
    "    t = %timeit -o -q numba_matmul(a, b)\n",
    "    time_range_numba_matmul.append(t.average)\n",
    "    t = %timeit -o -q np.dot(a, b)\n",
    "    time_range_np.append(t.average)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAESCAYAAADe2fNYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xl81NW9//HXySSTdbKvJGRjExCBEIFWUWRRa1HQ1opKUVARrRarLVeprXqtAr0uqLgUWVKs4vWnUtFrXQGtrVUQUARUICQQAiQhIfsyy/n98Z1MErKQwCQzyXyej0ceyXznzHdOJpO8c9av0lojhBBCuIufpysghBCib5FgEUII4VYSLEIIIdxKgkUIIYRbSbAIIYRwKwkWIYQQbiXBIoQQwq0kWIQQQriVv6crcKaUUqHAc0ADsFlr/bKHqySEED5NeePKe6XUamAaUKS1PrvZ8UuBpwATsFJrvUQp9UvghNb6baXU/2qtrznV+WNjY3V6eno31V4IIfqmr776qkRrHXeqct7aYskBlgNrGw8opUzAs8BUoADYopTaAKQAO53F7J05eXp6Olu3bnVnfYUQos9TSuV3ppxXjrForT8FSk86PBbYp7XO1Vo3AK8C0zFCJsVZpt3vRyk1Tym1VSm1tbi4uDuqLYQQAi8NlnYkA4ea3S5wHnsT+JlS6nng7fYerLVeobXO1lpnx8WdsiUnhBDiNHlrV1hbVBvHtNa6GpjTqRModTlw+cCBA91aMSGEEE16U4ulAOjf7HYKUOihugghhGhHbwqWLcAgpVSGUsoMzAQ2dOUEWuu3tdbzIiIiuqWCQgghvDRYlFLrgM+BIUqpAqXUTVprG3AH8D6wB3hNa72ri+e9XCm1ory83P2VFkIIAXjpOpbulp2drWW6sRBCdI1S6iutdfapyvWmwfseV19fT2lpKZWVldjtnVoiI4RPMJlMWCwWoqOjCQwM9HR1RAfKa63kFleRW1xNbkkVs3+UTkJ4ULc+p08FS1dmhdXX13Pw4EGioqJIT08nICAApdqamCaEb9FaY7Vaqaio4ODBg6Smpkq4eJjV7uBQaY0rPHKLq11fl1Q1uMr5+ynOGxjb7cEiXWHtOHLkCAEBAcTGxvZQrYTofUpKSrBarSQlJXm6Kn2e1prS6gZyS6pdLZD9zvA4eLwGm6Ppb3lMqJnMuFAyY8OMz3HG59ToEAJMpz+0Ll1hbehKi6WyshLZT0yIjoWHh5OXlyfB4kb1Njv5x2vILa4ygqNZK6S81uoqZzb5kR4bwuB4C5cOT3SFx4DYMCJCAjz4HfhYsGit3wbezs7OvuVUZe12OwEBnv3hCOHtAgICZPzxNGitKaqsZ39xy26r3OJqCspqaNb4ICE8kMzYMKadk9QiPJKjgjH5eWf3vE8FS1fJmIoQHZPfkY7V2+zkldSwt6iS/UVN4XGgpJqqepurXHCAiYzYUM5JiWDG6GQGOLuxMuJCCQvsfX+me1+NhRDCy9RZ7eQWV7O3qJK9x6qMz0VV5B+vwe5sfigF/SKCyYwL5edjUlqMgSSGB+Hnpa2P0+FTwSJ7hQkhzkRNg439RdWu4Nh7rIp9RZUcLG3qvjL5KdJjjLGPn45IYmB8GAPjwxgQF0ZQgMmz30AP8alg6coYixDCd1XV29hXVMXeY5XsK6rih2NGkBSU1brKBJgUGbGhDO8XwfRRyQxKCGNwgoX0mFDM/l65qUmP8algEaK5iRMnkpeXR15enqerIjykvNbKPlf3lfGx71glheV1rjJmkx+ZcaGMTo3imuz+DEoIY2C8hbSYM5u625dJsAg2b97MRRddBMCvfvUrli9f3qpMUVERKSkpWK1WLrzwQjZv3tzl58nJyeHEiRPcddddZ1plIbqkrLrBGRyVzu4r4+tjFfWuMkEBfgyMD2NcZgwD48MYFB/GoAQL/aOC8ZcA6RKfChYZY+lYUFAQr7zyCo8//nirldQvvfQSWmv8/U//LZOTk0NeXp4Ei+g2ZdUN/HCskh+c3Vg/OLuymq8+DzGbGBQfxvkD4xiU4AyQeAspUcF9agDdk3wqWGSMpWNXXnkl69at46233uIXv/hFi/vWrFnDZZddxscff+yh2gnRpLzW6gwOY/zjB+fXJVVNLZCwQH8Gxocx6ax4BsVbjBBJsJDUx2ZgeSOfChbRsaysLHbv3s2aNWtaBMuXX37Jrl27+NOf/tQqWD744ANWrVrFli1bOHLkCIGBgYwdO5bf//73XHjhha5y6enp5OfnAy3XPmzatImJEye6xjs2b97Mb37zGzZt2oRSiunTp7N8+XJCQkJYsmQJL774IkeOHGHYsGE888wznHfeea5z5eTkMGfOHNc5m5PxlN6pos5qjH84g2NvkREizbuwGlsgE4fEMdgZHoMTLPSLCJJ1Nh4iwdJFD729i92FFZ6uRgvD+oXzwOXD3XKuOXPmcPfdd1NQUEBKSgoAq1evJj4+nmnTprUqn5OTQ2lpKbNnzyYlJYXDhw+zcuVKJk+ezKZNm5gwYQIAy5Yt47777qOkpIQnn3zS9fihQ4e6vq6urmbSpElccMEFLFmyhC1btrB69Wrq6uqIiYnhiy++4M4778RqtfLYY49x+eWXk5+fj8Viccv3Ljynqt7G3mPG+EfzrqwjzQbRgwL8GBRv4byBsQxOsBghEm8hOVK6sLyNBItoYdasWSxcuJC1a9eyaNEiamtrefXVV7n55pvbHF958cUXCQ0NbXFs/vz5DB8+nMWLF7uCZcaMGSxbtoza2lpmzZrV5nOXlJSwcOFCfve737nOU1ZWxmuvvUZWVhaff/65a5udoUOHMn36dF555RVuvfVWd74EohvVNNic03ebxkB+OFbF4RNN03jN/n4MjAtjXEa0q/UxJEHGQHoTCZYuclfLwFvFxMRwxRVXkJOTw6JFi3jzzTcpLy9n7ty5bZZvHipVVVXU19djMpkYN24c//nPf7r03CaTiTvvvLPFsQkTJrB+/Xrmz5/fYu+2xsDau3dvl55D9Iw6q921/sMVIkWVFJTV0riheuM03jFpUVw7tr8rRFKjQ7x2DyzROT4VLDIrrHPmzJnDT3/6Uz777DNWr17N2LFjGTZsWJtl9+/fz+9//3vef/99Tpw40eK+rvZvJyUlERTU8joRUVFRAGRkZLR5/Pjx4116DuF+Dq35Kr+Ubw9XsPNwOd8eLmdvUZVrK5PGhYTnpERy9Zj+rnGQtOgQmcbbR/lUsMissM655JJLSE5O5qGHHmLTpk08//zzbZarqqriggsuoLq6mrvuuosRI0ZgsVjw8/Nj8eLFbNy4sUvPazK1v91Fe/c1v55QR0Fms9navU90nt3hoNbqoLbBTp3VTm2DnSMn6rh57eeAcR2Qs5MjmDw0nmFJEQxOCCM9NlQWEvoYnwoW0Tkmk4nZs2ezePFigoODmTlzZpvlPv74YwoLC1m9ejVz5sxpcd/999/fqnx3z9CJjo4GoLS0tNV9Bw4ckMsgdJHN4aCuwU6t1U5tg4Naq516W9MW+QEmP4IDTFiC/HlxdjZnJ4eTGC4zsYQEi2jH/PnzMZvNZGZmEhER0WaZxlbEyVch/eCDD/jiiy9alQ8LC6OsrAytdbf88Rk8eDAAH330EVdddZXr+Lp16ygsLCQtLc3tz9lX2OxGcNQ6WyG1VjsNNofrfrPJj6AAE5EhAQQHmAg2m1ytkNriAMYNTfBU1YUXkmARbUpNTeXBBx/ssMz5559PYmIi99xzD3l5eaSkpLBjxw5eeuklRowYwc6dO1uUHz9+PO+88w533HEHP/7xjzGZTEyaNIn4+Hi31HnIkCFMmTKFv/zlL2itGTVqFDt27GD9+vUMHDgQq9V66pP4AKvd4QqPxi6tBnuzEPE3WiLRoWYjRAJMMhYiukSCRZy2yMhI3n//fRYuXMgzzzyDzWZjzJgxvPvuu6xatapVsNx1113k5uby+uuv88ILL+BwONi0aZPbggWMrWfuvPNOXn75ZV566SUmTJjApk2buO2223xucaTWGptdt2qJWJuFSKC/iRCzPzFmI0yCJESEG6iTuzF8QXZ2tt66dWuHZfbs2dNi8Z4Q3s7mcFBTb6emWWvE5jBCRGGESLDZCI9gs4ngAD9MfmceIvK74juUUl9prbNPVa7Xt1iUUpnA74EIrfXPPV0fIXqKze6gusFGdb2d6nobtVZjYF2hCAzwwxLk7wwQI0xkbYjoKR4NFqXUamAaUKS1PrvZ8UuBpwATsFJrvaS9c2itc4GblFKvd3d9hfAkq80Ikqp6I0waZ2j5KUWI2URCeBChZn9CzCZZoS48ytMtlhxgObC28YBSygQ8C0wFCoAtSqkNGCGz+KTHz9VaF/VMVYXoOVprGuwOV2ukusHmmqVlUoqQQH+iQgMINRutEj+Z4iu8iEeDRWv9qVIq/aTDY4F9zpYISqlXgela68UYrZvTopSaB8wDY8aTEN5Ea029zWGESL2d6gaba5Dd308RGuhPTGggYYFGt5asFRHezNMtlrYkA4ea3S4AxrVXWCkVAzwCjFZK3ecMoFa01iuAFWAM3ruvukJ0ndaaOqvdFSLV9U0D7QEmP0LN/oQGmggN9CfQ30+CRPQq3hgsbf0GtRsEWuvjwPxOnVj2ChMe4tCa2oamEKmpt2F3zsg0+xsD7aGB/oSaTZglSEQv543BUgD0b3Y7BSj0UF2EOC0Oh6bG6hwfqbdR02DH4QySQH8TESEBziDxx+wv60ZE3+KNwbIFGKSUygAOAzOB69xxYtmEUnQXh0M7WyPOFonV7trqJsi5ij3UbCIk0F82ZBR9nqenG68DJgKxSqkC4AGt9Sql1B3A+xgzwVZrrXe56fmkK0y4ReNge2Wdjco6K9UNRpAoFMFmE7FhZtfUX1nJLnyNp2eFXdvO8XeBd7vh+aTFIk6bze6gqt5GVZ2NyvqmWVtB/iZiQs1YgvwJMfvLQkTh87yxK0wIr6CdA+6V9TYq62zUNtjQgMlPERbojyUokLDAABkjEeIkPhUs0hUmTsVqN7q3quqsVNbbXFdBDDGbiLMEOVslso5EiI741L9aWuu3tdbz2ru+iPA9EydOJC0tnSPltfxwrJI9RyooKKuhqt5OeFAAqdEhDEsKZ2C8hcSIIEID/b0+VG688Uavr6Po23wqWJRSlyulVpSXl3u6Kl5l8+bNKKVQSnHHHXe0WaaoqAiz2YxSiokTJ572c+Xk5LBs2bLTfrw71FvtlFTVk1dSTXW9DZvDQUlVA/5+isSIIAbFhzE0yUL/6BAiQ8yuwff09HSUUpx//vltnrfxD3pJSUlPfjtCeB2fChZpsXQsKCiIV155hfr6+lb3vfTSS2it8fc/s95TTwSL3aGpqLVyuKyW745W8P2xSgpP1FJns+Nv8sPfz49hSeFkxoURbwki2Nxxq+Rf//oXb731Vg9+B0L0Lj4VLKJjV155JWVlZW3+0VyzZg2XXXYZgYGBHqhZ1zQOuhdV1pFbXMXuIxXkHa+mrKaBIH8T/SKDGZJg4azEcOd2KXR6JldaWhpxcXHcd9992O32Uz9ACB8kwSJcsrKyGDlyJGvWrGlx/Msvv2TXrl3MmTOnzcd98MEHXHPNNWRmZhIcHExkZCQXX3wxn3zySYty6enpfPLJJ+Tn57u63pRSbN68GTDGO9LT08nLy+PKK68kMjKSqKgobrzxRqqqqnA4HDz66KNkZGQQFBREVlYW//rXvwBjKvCJmgaeWP4X/Pz8eOWtf3C0vA6bQxMbZiYzNpQ7rp/OxOzhxIYFEhhgOq3XKCwsjPvvv589e/aQk5NzyvIdjXcopbjxxhtdt/Py8lBK8eCDD/Laa68xatQogoODGThwoOtncvDgQX7+858THR2NxWJh1qxZVFZWtnn+4uJiZs+eTUxMDKGhoUyePJnt27e3Kvfcc89x8cUXk5ycjNlsJikpiVmzZvncFTeF+8issK76x71wdOepy/WkxBHwk3YvWdMlc+bM4e6776agoICUlBQAVq9eTXx8PNOmtb25dE5ODqWlpcyePZuUlBQOHz7MypUrmTx5Mps2bWLChAkALFu2jPvuu4+SkhKefPJJ1+ObX32wurqaSZMmccEFF7BkyRK2bNnC6tWrqaurIyYmhi+++II777wTq9XKY489xrRpl7Ppq134mYMBXBe7igsLZGhSeLescp8/fz5PPfUUDzzwANdddx3BwcFuPf8777zDCy+8wO233050dDSrVq1i7ty5mM1mFi1axKRJk3j00Uddr01QUBArV65sdZ5LL72U6OhoHnzwQY4ePcry5cu54IIL+Pzzzzn7bNflj3jssccYP348v/71r4mOjubbb79l5cqVbNy4kZ07dxITE+PW70/0fT4VLLJA8tRmzZrFwoULWbt2LYsWLaK2tpZXX32Vm2++ud3xlRdffJHQ0NAWx+bPn8/w4cNZvHixK1hmzJjBsmXLqK2tZdasWW2eq6SkhIULF/K73/3OdZ6ysjJee+01srKy+Oxf/6bWBmU1DQTH9WfBTdfz9huvMe/WW7EE+tMv0vgjHx4c0G1bp5jNZh5++GGuv/56nnrqKe699163nn/Pnj3s3r2btLQ0AK655hr69+/PL3/5Sx577DHuvvtuoOm1Wbt2LcuWLSMsLKzFedLS0njjjTdcLaarrrqKc889l9/+9re89957rnI7d+5s9fO74oormDJlCqtWrWLhwoVu/f5E3+dTweIWbmoZeKuYmBiuuOIKcnJyWLRoEW+++Sbl5eXMnTu33cc0/6NUVVVFfX09JpOJcePG8Z///KdLz28ymbjzzjtbHDv//PNZv349V8+ay76SWhxaYzb5cfHkiwAoP3aIhPAggB6bZnvttdfy+OOPs3TpUubNm0d0dLTbzj1jxgxXqADExcUxZMgQdu3axa9+9asWZSdMmMD69evJy8tr0QoBWLhwYYvXY8yYMUydOpWPPvqIqqoqVxA1/vwcDgeVlZVYrVZGjhxJREQEX3zxhdu+L+E7ZIxFtDJnzhz27t3LZ599xurVqxk7dizDhg1rt/z+/fuZOXMmUVFRWCwWYmNjiYuL491336WsrKxLz52UlERQUBBaa6rrbRw+UUutMlohsUkpRIUEMCAujCGJFs5KSwLg+PHjp//NnialFEuWLOHEiRM88sgjbj13ZmZmq2NRUVEkJSW1mjwRFRUFtP0aNO9ibDRs2DDsdjv5+fmuYxs3bmTixImEhoYSGRlJXFwccXFxlJeXd/nnJwT4WItFVt53ziWXXEJycjIPPfQQmzZt4vnnn2+3bFVVFRdccAHV1dXcddddjBgxAovFgp+fH4sXL2bjxo1dem4/k4mjFXWcqGmgwebATykCA4z/fzLiLSRHhbR6TOMuwtBxi8Vms3WpLqcydepUpkyZwrPPPsuCBQvaLNNefTqqi8nU9sSC9o5Dy9egIyeX27JlCxdffDEDBw5kyZIlZGRkEBwcjFKKmTNn4nBefEyIrvCpYJExls4xmUzMnj2bxYsXExwczMyZM9st+/HHH1NYWMjq1atbzRq7//77W5Vv6w+t1e6gvMZKrdWOze6gqKKOsEB/4i1BRAT7szXM+C+9M9d1b+ySKi0tbXXfgQMHCAgIOOU5umLp0qVkZ2fzhz/8oc3vrXl9mneX5ebmurUebdmzZw/jx49vdcxkMrm62l555RXsdjv/+Mc/yMjIcJWrrq6W1oo4bdIVJto0f/58HnjgAV544QU6WlDa+F/0yf8Jf/DBB232z4eFhVFWVobN7qCspoEDJdV8d6SSwvJa0MZ6kqGJxmLF6FAzJr+uvUUHDx4MwEcffdTi+Lp16ygsdP/14rKyspg5cyZ/+9vf2Lmz9WzB9urz+OOPu70uJ/vzn//c4ueybds2PvroIyZPnuwaX2nv5/foo49Ka0WcNp9qsYjOS01N5cEHHzxlufPPP5/ExETuuece8vLySElJYceOHbz00kuMGDGixR9brTWjss/lnXfeYdZNtzJyzLmY/f25ZOpkBqWlEGw24acUAWewW/CQIUOYMmUKf/nLX4znGzWKHTt2sH79egYOHIjVaj3tc7fnkUce4Y033mDbtm2t7rv22mtZtGgR8+bN47vvviMmJoZ//OMfPbLtS35+PpdccglXXHEFR44cYfny5QQHB/M///M/rjJXXnklTz75JJdddhnz5s3DbDbz4Ycf8s033xAbG9vtdRR9k7RYxBmJjIzk/fffZ9y4cTzzzDPcc8897N69m3fffZesrCwAahpsFJ6oZc+RSn567U1cOXMWH7/7FvffdRu/vf0mSgoOEHSaCxbb8tJLL3HVVVfx8ssvuwJv06ZNJCcnu+05msvIyGD+/Plt3hceHs67777L8OHDefTRR3nggQfo169fi+m+3eW9994jISGBBx54gCeffJIxY8bwySefcM4557jKnHfeebzxxhuEhobyhz/8gQcffJDg4GA++eSTVlOQhegs1dlBv74kOztbb926tcMye/bsaXNWjeicequdE7VWTtQ0UG9zoJQiPMifyBDjglidGS8RvYP8rvgOpdRXWuvsU5Xzqa4wmRXWvax2B+W1Vk7UWKlpMGY9hQX6E+cchO/qeIkQonfyqWCRWWHuZ3doKuuslNVYqaqzodEEBZhIiggiItgsV1cUwgf5VLAI92mw2SmqqOdErRWH1gSY/Ii1mIkKMbt1vEQI0ftIsIgusdkdFFXWc7y6AQVEhgQQGWImVC7XK4RwkmARnWJ3aEqq6imurEdriAoNICE8qNs2ehRC9F4SLKJDDq0prW6gqKIem8NBRLARKNLdJYRojwSLaJPWmvJaK0cr6miwOQgN9CctPITQQHnLCCE61if+SiilZgA/BeKBZ7XWH3i4Sr1aZZ2Vo+V11FrtBAWYSI8NxRLY8XXghRCikcc7yJVSq5VSRUqpb086fqlS6nul1D6lVIdXUtJa/11rfQtwI3BNN1a3T6tpsJFbXMWBkmrsDk3/6BAGxYcRHhQgoSKE6DRvaLHkAMuBtY0HlFIm4FlgKlAAbFFKbQBMwOKTHj9Xa13k/Pp+5+NEF9Rb7RyrqONErRV/P0VSRDAxYWZZHS+EOC0eDxat9adKqfSTDo8F9mmtcwGUUq8C07XWi4FWF15Xxr/TS4B/aK1b7wRolJkHzANjg0VhrJQvqqintLoBpSDeEkScpes7CgshRHMeD5Z2JAOHmt0uAMZ1UP5OYAoQoZQaqLV+4eQCWusVwAow9gpzY117HbvDQXFlAyVVxtTh6NAA4mXqsBDCTbw1WNrqg2k3DLTWTwNPn/KkPr5XmENrSqsaKKpsmjqcGB5EoEwdFkK4kbf+i1oA9G92OwVw/1WafITWmrKaBn44alxQKyjAj4HxYaTFhBIYYGLz5s0opVBKsXLlyjbPoZRi2rRWvZBCCNGKtwbLFmCQUipDKWUGZgIbzvSkWuu3tdbzOroiYl+itaaizsreoioOldZg8lNkxIaSERtKiLntxuoDDzxAbW1tD9dUCNGXeDxYlFLrgM+BIUqpAqXUTVprG3AH8D6wB3hNa73LDc91uVJqRXl5+ZmeyuvV1NvILakmr6Qah9akRocwMD4MSwdTh7OzsyksLGTZsmU9XFshRF/i8WDRWl+rtU7SWgdorVO01qucx9/VWg/WWg/QWj/ipufq8y2WOqud/OPV7Cuuot7qoF9kMIMTLESGmE+5FuUXv/gFY8aMYenSpRw/frzDskopbrzxxlbHc3JyUEqxefNm17EHH3wQpRS7d+/mrrvuIikpidDQUCZPnsz3338PwJtvvklWVhbBwcGkp6ezYsWKdp/zo48+Yvz48YSEhJCYmMiCBQuorq52lXviiSdQSrW6zjxAfX090dHRTJ48ucPvTwhx+rx18N5rLf1yKd+VfufparRwVvRZ3D3mdxyrqKOs2opSkBAeRGxYICa/zq9FUUqxdOlSpkyZwiOPPMITTzzh1nrecMMNhIWFsWjRIoqLi3n88ce55JJLePjhh1m4cCG33XYbc+fOZdWqVdx6660MGzaM888/v8U5tm3bxuuvv84tt9zC7Nmz2bRpE08//TTffvstH374IX5+ftxwww0sWrSIVatWMWXKlBaPX79+PWVlZdx0001u/d6EEE18Klj64qwwDVTX2/j+aCUaiAkzE2cJPO2pw5MnT2bq1Kk899xzLFiwgLS0NLfVNTExkQ0bNrhaTrGxsSxYsIDbb7+dXbt2udYXXXPNNfTv359nn322VbDs3LmT9evXM2PGDABuv/12FixYwNNPP81rr73GzJkziYmJ4aqrruLNN9+ktLSU6Oho1+NXrVpFVFQUV111ldu+LyFESz4VLO64guR/jf0vN9bo9DkcmuPVDRRV1mF3aMKDA0gIDyTQ/8ynDi9dupQxY8bwhz/8gbVr1576AZ3061//ukV33IQJEwCYPn16i0WrcXFxDBkyhL1797Y6x5AhQ1yh0ujee+/l6aefZv369cycOROAefPmsW7dOl5++WXuvPNOAPLy8vj444/51a9+RVBQkNu+LyFESx4fY+lJfWXwvrymge+PVXKkvJYQsz+D4sNIjQ5xS6gAjB49mmuvvZaXX36Zb775xi3nBMjMzGxxOyoqCoCMjIxWZaOiotoc5xk6dGirY0lJSURGRpKbm+s6NnHiRAYPHsyqVatcx9asWYPWmptvvvm0vwchxKn5VLD09sF7rTVFlXXkl9bg76fIdE4dDm5n6vCZ+NOf/oS/vz//9V9da6HZbLZ27zOZ2g6+9o5r3XpNbHsTENoqe8stt/D111/z1Vdf4XA4yMnJITs7m5EjR7ZbRyHEmfOpYOnNtNYcKa/jaHkdkcFmBsSHERYU0G3Pl5GRwW233cZ7773Hpk2bWt0fHR1NaWlpq+PNWw3dYffu3a2OHTlyhPLy8lYtohtvvBGz2cyqVav48MMPOXjwoAzaC9EDJFh6AYfWHCqtoaSqntiwQPpHB/fIzsP3338/4eHhbbZaBg8ezOeff05NTY3rWFlZGWvWrOnWOn3//ff8/e9/b3Fs6dKlAK3GXmJjY5kxYwavvPIKy5cvJyQkhOuuu65b6yeE8LFg6Y3LEPxOAAAgAElEQVRjLHaHg7ySak7UWkmKCCIpIqjHro0SGxvL7373O7Zs2dLqvjvuuIOSkhImTZrEs88+yyOPPMLo0aPdOousLSNGjGDWrFn85je/4bnnnuPqq6/mqaee4sILL+Saa1pfimfevHmUl5fzzjvvcPXVVxMeHt6t9RNC+Fiw9LYxFqvdQW5xNdX1dvpHhRBn6blQaXT33XeTlJTU6vj111/Pn//8Z44cOcLdd9/N3/72N/74xz8yf/78bq1PVlYWf//73/n3v//NPffcw6effsodd9zB22+/jV8b2/1PmjSJxunl0g0mRM9QbQ169nXZ2dl669atHZbZs2dPmzOQekq91c6B49XY7JrUmBDCu3E8pbdQSnHDDTeQk5PTpccNHz4cu93Od99518LWvsLTvyui5yilvtJaZ5+qnE+1WHqLmgYb+4urcTg0mbGhEipnYOPGjezevZt58+Z5uipC+AyfWiDZG1beV9ZZyT9uTCdOjw0jSK6Vclo2btzI/v37Wbx4MXFxcdxyy2mviRVCdJFPtVi8fYzlRE0DecdrMPv7MSBeQuVM/Pd//ze33XYbYWFhvPHGG1gsFk9XSQif0aUWi1LKAvwGuBhIAGZrrT9XSsUCt2Nsby8d2aehuLKeI+W1hAb6kx4TItedb0NXxgOb764shOhZnQ4WpVQc8BmQCexzfg4G0FqXKKVuACKBu7uhnn2W1pqjFXUUV9YTERxA/6gQ/LqwI7EQQnibrrRY/gQkAuOAg0DRSfe/BfSpi1xorbt1eq9Daw6X1VJW00BMaCD9Int+OrEQZ8IXZ5WKU+tKf8s04Dmt9TaM3dpPlkvL69T3aiaTCavV2m3ntzs0+cdrKKtpICE8SEJF9EpWq7Xdvd6E7+pKsMRidIG1xwF49V7kXVl5b7FYqKio6JZ62OwODpRUUVVnJSUqmIRwCRXRO1VUVMjECNFKV4LlKDCgg/tHY3SRea2uzAqLjo6mrKyMkpISGhoa3Nbkb7DZ2V9cTZ3VQVpMKNGhgW45rxA9RWtNQ0MDJSUllJWVtbiQmhDQtTGWd4GblFLPAA3N71BKjQNmA8vcWDePCgwMJDU1ldLSUvLy8rDb7Wd8TqvdQUmVEVIxYYEcrvDjsBvqKkRPM5lMWCwWUlNTCQyUf45ES10JloeAK4DtwAaMcZYblFK3AFcBhcBSt9fQgwIDA0lKSmpzr6yu+nz/ceat3UpYkD9/nTuWwQnSfSCE6Js63RWmtT4KjAe+AOYCCvgl8AvgA2CC1rr1BToE7+48wg2rvyQxIog3bvuxhIoQok/r0gJJrfUhYLpSKhwYghEu+yRQ2vfS53n8ccMuslKjWHVDNpEhZk9XSQghutVp7RWmta4AWl+kQ7horXn8gx9YvmkfU4YmsPy60bJFixDCJ3Q5WJRSIUA6EIPRYmlBa/3pmVerS/UZCizAmA79sdb6+Z58/rbY7A5+v/5b/nfrIWae258/zTgbf5Ns0SKE8A1d2dIlFHgSY/ZXW/u4K4wB/U7/W66UWo2x8LJIa312s+OXAk85z7VSa72kvXNorfcA85VSfsCLnX3u7lLbYOfOddv5aM8x7pw0kLunDpY1KkIIz2uogdJciM4Ec0i3PlVXWiwvANcD64F/AmVueP4cYDmwtvGAUsoEPAtMBQqALUqpDRghs/ikx8/VWhcppa4A7nWey2NO1DRw01+3su1gGQ9PH84vf5TuyeoIIXyNrQHK8qB0PxzfB8edn0tzocK5uGHOPyDtx91aja4Ey3RgldbabRe20Fp/qpRKP+nwWIwJAbkASqlXgela68UYrZu2zrMB2KCU+j/glbbKKKXmAfMAUlNT3VL/5gpP1DJ79ZccPF7Ds9dlcdmIM5+iLIQQrTjscOKgMzz2NwuP/cZx7WgqGxwNMQMg4wKIHmB8HTu426vYlWCx0jMD9snAoWa3CzA2vmyTUmoixjqaQIxFnG3SWq8AVoBxaWJ3VLTRD8cquWH1l1TV2fjr3LH8aECMO08vhPA1WkPlkdatjuP7jBaJvdkadXOYERj9smDEL4yvYwYaXV4hntkVoSvBshHjD/yKbqpLo7YGJNoNAq31ZmBzp07cDVeQ3JpXytycLQQGmPjfW3/EsH7hbju3EKIP0xpqjrdscRzfB8dzja+tNU1lTYFGUMQOhiE/cQaHswUSlgBeNo7blWC5B/hUKbUAY5fj7tr6t4CWuySnYKzq9zof7j7GHa9sIzkymL/OHUv/6O4dEBNC9EINNVDyfVOANA+SumYb4ioTRKUZoZExwQiSmIFGeISnQC+6+J/qyuaKSqnrMAbaHcAR4OQNtLTWuqONKts6ZzrwTuOsMKWUP/ADxrVdDmN0v12ntd7VlfN2JDs7W2/duvWMzvHqlwdZtH4nI1IiWX1DNjFhsl+SED7NbjO6q4p2Gx/HdkHRHuNY806XiP4tQ6Ox9RGVBqa2Jtx6D6XUV1rr7FOV68p04xuBVRgbUH6PG2aFKaXWAROBWKVUAfCA1nqVUuoO4H2MmWCr3RUq7ugK01rzzMZ9PPHhD1w4OI7nZ2URYj6tdaZCiN5Ia6gobBYgu6FoFxT/APZ6o4zyM8IjYRiMuBriz4LYIRCdAQHBnq1/D+h0i0UptReoAC7RWpd0a6262em2WOwOzYMbdvHSf/K5KiuZpT87hwBZ+ChE31VbZrQ6ju1yBske43PzLixLEsQPg/ihkDDc+Bx3Vp8MELe3WDBmaz3Z20PlTNz7xjf8v68KuPXCTO699CxZ+ChEX2GtheLvncGxy9kK2QOVzYZ3AyOM0Dj7Z84gcYaJh2ZeebOuBMv3QK9+Bc+0K+yqrBSGJoUz9/wM91ZMCNEzHHYoPWCER/OWSGlu0/oPUyDEDTbWfjRvhYQne93sK2/Vla6wqzG2WRmrtS7o1lp1M3cM3gshvFxNKRzZAUe/bWqJFH8PtjpnAWWMgzQPj/jhxjGTjJu2pTu6woZizNLao5RaDxyg7VlhD3fhnD2qO9axCCG8QH0VHPkaCrdB4XY4vA3KDjTdH5ZoDKSfe7PRhZUwzBhM7+Y9s3xVV1osjlOXQmutvX5veGmxCNGLWevg2LdNAVK4zWiJNE7pjegP/UYbH8lZkHiOjIO4SXe0WGRgQQjRs+w2KN7jDJDtRogc2w0O5/rs0DhjK5PhVxqf+42GsDjP1ll0Pli01vndWREhhI9zOIzV6I2tkMPb4OhOsNUa9wdFGMHx4zuMEEnOkgF1L+VTI1QyxiKEl9Da2Im3sRVyeJsxRlJfYdwfEAJJIyF7blOXVlRGr9rWxJe1O8ailPojRqflI1prh/P2qXj14H0jGWMRoodVHmsKkMYwqTlu3OcXAIlnN7VC+mUZmy3KzCyv09kxlo6CxYERLMFa6wYZvBdCdEr1cWOab+H2po/Gi0wpP4gbCsnOwfV+WcZUX3/Za683cMfgfQaA1rqh+W0hhHCpPGp0YTV+FO6AimbL3KIHGFcrbGyNJI4Ac6jn6it6RLvBorXOV0qlKqWCtda1fWHwXsZYhDhNWhutjsIdLYOk6qizgDJ26U0dD/1GGeMjiedAcKRHqy08o8N1LEopO/BLrXWbl/vtraQrTIgOaG1cpfDI10aXVmOINI6JKD9jk8Wkkc6PUcYYSaDFo9UW3c9d61hkHp8QfVnjFN/GECncAUe/adq918/f2OpkyGVNIZIwXFasiw7JtAshfIXdBiU/NOvK2mGsE2moMu43BRqhMfyqpu6s+GEysC66TIJFiL7I1mCsWG8+HnL026bFhgEhxkD6qOuaWiJxQ7z+Coaid+hMsExwXi64U7TWa8+gPkKI01F9HPL/BXmfwaEvjK3g7c4JnWZL02LDxnGR2EHg5/UrA0QvdarB+8a1LJ06F16+jqXZrLBb9u7d6+nqCHH6qkuagiTvMyNIwGiJpGQba0QaWyKyYl24iTs3oVwB/OfMq+R5Wuu3gbezs7Nv8XRdhOiSjoKk/zjjqobpE4xA8Td7tq7C53UmWP7Z16YbC+H1OgqS1PEw4udGkCSNkiARXkcG74XwBtUlTSGS95kx8A4SJKJXkmARwhOqilu2SE4OknOuburakplaopeRYBGiJ7QbJKHOIPmFM0hGSZCIXq/DYNFay1QSIU5HVTHkN+/a+s44LkEifECfaLEopUKBT4EHtNbveLo+wgdVH4e8TzsIkmskSITP8GiwKKVWA9OAIq312c2OXwo8BZiAlVrrJac41X8Br3VbRYU4ma0eDv4H9m+E3E1w5BtAG0GS9iMYOdM52D5SgkT4HE+3WHKA5YBrtb5SygQ8C0wFCoAtSqkNGCGz+KTHzwXOAXYDQT1QX+GrtDam/O7fCPs3Qf6/je1R/Pwh5Vy4aBFkTpTBdiHwcLBorT9VSqWfdHgssE9rnQuglHoVmK61XozRumlBKXUREAoMA2qVUu9qrVtd7VIpNQ+YB5CamurOb0P0VZVHjRDJ3QS5m6HqmHE8djBkzYYBF0H6+bJdvBAn8XSLpS3JwKFmtwuAce0V1lr/HkApdSNQ0laoOMutwNhFgOzs7M5uUyN8SUO10RJpDJPGRYkhMUZrJPMiI0wiUjxZSyG8njcGS1vXgDllEGitc055YrmCpGjOYTd2/c3dZITJoS+MjRtNgcaA+5QHYcAkSBghe20J0QXeGCwFQP9mt1OAQg/VRfQ1Jw4aIbJ/Ixz4BGrLjOMJZ8O4W41WSeqP5EJWQpwBbwyWLcAgpVQGcBiYCVznjhPLJpQ+qK4cDvyzqVVSut84bkmCwT8xurYyJ0JYvCdrKUSf4unpxuuAiUCsUqoAYx3KKqXUHcD7GDPBVmutd7np+aQrrK+zW+HwV03jJAVbQduNrVLSz4dzbzbCJO4sUHLlbSG6Q4fXY+mrsrOz9datWz1dDeEux/c3TQM+8Ck0VALKmPo74CJjnCRlrGzeKMQZcuf1WITwPtZa2P0WbF0Dh5yXC4pMhRE/M8ZJMi6AkGjP1lEIH+VTwSJdYX1AyV4jTL5+xRh4j86Eqf8NZ00zvpbuLSE8zqeCRQbveylbPex5G77Kgbx/Gqvdz5oG2XMg/QKZCiyEl/GpYJEWSy9TmmuEyfaXoabE6Oqa/EcYNQssCZ6unRCiHT4VLNJi6QXsVvj+XaO7K3cTKBMM+QmMmWMMwkvrRAiv51PBIrzYiYPw1V9h+0vGnlzhyTBxEWT9EsL7ebp2QogukGARnmO3wd4P4Ks1sPdD49igi42xk4FTwSRvTyF6I5/6zZUxFi9RUQjb1hofFYchLBEu+K2xY3Ck7DwtRG8nCyRFz3DYjUWMW9fAD+8Zq+EHTILsuTD4UrmGiRC9gCyQFN6h8pgxbrLtr8Y4SkgsnPdryLoBojM8XTshRDeQYBHu53AY13/fuhq++z9w2IzL9E55yFh/IlurCNGn+VSwyBhLN6sugR0vG2tPSnMhOArGzTemCsfKay6Er5AxFnFmtIb8fxljJ3s2GBfKSv2xMbNr6BUQEOTpGgoh3ETGWET3qq+C7X+Draug5AcIjDAG4sfcCPFDPV07IYQHSbCIrqkugS/+Al+ugLoTkJwN05+D4VfKVReFEIAEi+issnz4fDlsewlstcYg/Hl3Qf9zPV0zIUQzNoeNyoZKKhoqqKivMD43+/ryAZeTGJrYrXWQYBEdO/ot/Osp+PYNUH5wzjXGdOG4IZ6umRB9ltVhNcKhjWBo9fVJt6ut1R2ee2TcSAkWd5JZYZ2kNeT/Gz57EvZ9COYwGH8bjL8dIpI9XTshepVqazV55XmU1pV2OiRqbDUdnjPYPxiL2UK4OZxwczhJYUkMMQ8xbgeGu45HBEa4vg4PDMdithBoCuz271lmhYkmDoexs/C/lkHBFmMx4/j5xnXig6M8XTshvFpZXRm55bnsP7GfA+UHXF8fqznWZvlg/+BWQdDqdjtfm02eWQsms8JE59kaYOdrRpdXyQ8QmQaXPQajZ0FAsKdrJ4TX0FpTVFPE/nIjPPaf2E9ueS4Hyg9QWlfqKhfsH0xGRAZjE8eSGZlJRngGsSGxLQIkoA9vYyTB4svqK43FjJ8/B5WFkDACfrYKhs2QnYWFT7M77BRWFRqtjvL95J7IJbfc+Gg+hhFuDmdA5AAu6n8RmRGZZEZmkhmRSWJoIn7Kd68dJH89fFFVMXzxAmx5EerKje1Wpj8DAybLNeOFT7HareRX5LtCozFA8iryqLfXu8rFBceRGZHJFQOuIDMikwGRA8iIyCAmKAYlvzOtSLD4ktIDxpTh7X8zriM/dBqc9xtIGePpmgnRrWqsNeRV5LnGPxq7sA5VHsKu7a5yyWHJZEZkMj5pvCs8MiMzCTeHe7D2vY8Eiy848o0xIL9rvXGp35Ez4bwFEDvI0zUT4oxZHVaO1x6nqKaI4ppiimqNz8dqjlFUU8TBioMUVhe6yvsrf/qH92dA5ACmpk0lMzKTAREDSI9IJ9hfxhTdodcHi1JqIvAwsAt4VWu92aMV8hZaQ94/4bNlsP9jY8rwj35lTBmWS/2KXsChHZTVlVFcW9wUGjVFruAoqimiqKaI0rpSNC1nt5qUidjgWOJD4hkZP5IrI65kQOQAMiMySbWk9umBc2/g0WBRSq0GpgFFWuuzmx2/FHgKMAErtdZLOjiNBqqAIKCgG6vbOzjsxlb1nz0JhdsgNA4m/xGyb4LgSE/XTgi01lRZq1ytivaCo7i2GJvD1urx0UHRxIfEExccx7CYYcbXIXHEB8e7vo4OivbpwXNP83SLJQdYDqxtPKCUMgHPAlMxgmKLUmoDRsgsPunxc4F/aq0/UUolAE8A1/dAvb2PrR6+fhX+/TQc3wdR6fDTJ2DUdTJlWPQYrTUn6k9wsPIgR6qOuFoVzVsZxbXF1NpqWz3WEmAxAiIknnMTzyUuOI64kDgSQhJcwREbHCutjV7Ao8Gitf5UKZV+0uGxwD6tdS6AUupVYLrWejFG66Y9ZUC7S0qVUvOAeQCpqX3ouup1FfDVGmPKcNVRSDwHfr4ahk6XKcOi25TXl3Ow4iD5lfnG54p81+3KhsoWZQNNga1aGI23G4MjNjiWkADZxLSv8Ma/PMnAoWa3C4Bx7RVWSl0FXAJEYrR+2qS1XgGsAGPlvVtq6klVRfCf52HLKqgvh4wL4MrnIfMimTIs3KKqoarN4DhYcZAT9Sdc5RSKpNAkUsNTuSzjMlItqaSFp9EvrB/xIfGEm8NlSq6P8cZgaesd2G4QaK3fBN7s1In7wl5h1jrY+DB8+aJxUa1hVxgzvJJlyrDouhprDQcrmwVHRb7rdvOV5AAJIQmkhacxJW0KaZY0UsONAEmxpPTI/lOi9/DGYCkA+je7nQIUtlPWtxTtgddvgqJdMGoWnP8bueSvOKVaWy2HKg+1Co6DFQcpri1uUTYuOI7U8FQm9p/oanmkhqfS39JfpuKKTvPGYNkCDFJKZQCHgZnAde44sdb6beDt7OzsW9xxvh6jtXGlxvd/D4EWuP51GDTV07USHqS1pt5eT0VDRatrbxyvPd6iC+vkTRCjg6JJC0/jx/1+7AqOtPA0Ui2pMs4h3MLT043XAROBWKVUAfCA1nqVUuoO4H2MmWCrtda73PR8va8rrPo4bLgTvv8/Y8uVK1+AsHhP10q4gUM7XKHQ4nN9s7Bo9nHyfVaHtd1zRwVGkRqeytjEsU3BEZ5KmiWNMHNYD36XwhfJtvneLPcTWH+rcTngqQ/BuNvAT+bmexutNSW1JRyuOtwqBFqFRrP7qqxVrRb2NWdSJtc1N1zX3ggMb33M+WExWwgPDCcqKEq2IBHdQrbN783sVtj0iLFqPmYgXPe/kDTS07XyaVprjtcd51DloVYD3QcrDrZ7YaaTL8iUGJLIoMhB7QaExWwhIjACi9lCiH+IzKYSvZJPBUuv6AorzTUG6Au3QdZsuHQJmEM9XSufoLWmrL6MgxUHW82UOlR5iCprlausSZlIDksmNTyVMQljSLWkkmJJITIwskVYyGI+4YukK8ybfP0q/N894GeCy5+G4TM8XaM+qby+nPyK/FYzpA5WHKTS2rS4z0/50S+0X6sB7rTwNJLCkgjwk9AQvkW6wtrgtS2WugojUHa+BmnnwVUrICLF07Xq1SoaKtpc2Hew8iDl9eWucgpFv7B+pFpSuSzzMtLC00gLT6O/pT8pYSnS4hDiNEiLxdMObYE3boLyAph4L0xwtljEKWmtOVp9lK9LviavPM8VIIcqDlFWX+Yqp1Akhia6ZkWlhqe6Wh4plhSPXT9ciN5GWizezmGHz56ATYshIhnm/ANS2925RgA2h43vy75nR9EOthdtZ3vRdopqilz3N64Mn5Q6qcX02v7h/WVluBA9SILFE8oL4M1bIf8zOPtnMO1JCIrwdK28TmVDJV8Xf832ou3sKNrBzpKdrl1xE0MTyYrPYlT8KEbFjyIzIlNWhgvhJXwqWLxijGX3BmPBo90KM56HkdfKppEY3VoFVQUtWiP7T+xHo/FTfgyJGsKMgTMYHT+a0fGjSQxN9HSVhRDtkDGWntJQA+/fB1/lQL/R8LNVEDOgZ+vgRax2K3tK97haIzuKd1BSWwJAWEAYI+NGMjJ+JKPjR3NO7Dmy1YgQXkDGWLzJ0Z3w+lwo+cHYifii+8HftwaMy+vLXa2RHcU7+LbkW+rt9QAkhyUzLmkco+NGMyp+FAMjB2KSCQxC9FoSLN1Ja/jiBfjwjxAcDb/8Owy4yNO16nZaa/Ir8tletN01RpJbnguAv/JnaMxQrh58NaPjjSCJD5G9z4ToS3wqWHp0jKWqGP5+G+z7EAb/BKY/C6Ex3f+8HtBgb2DX8V2uFsnXxV+7ruVhMVsYFTeKaZnTGBU/irNjz5ZBdiH6OBlj6Q77PoL1t0FdOVzyCJx7c58aoC+tKzXGRZxBsuv4LtdOu2nhaYyMG+kaZM+IyMBPycaZQvQFMsbiCbZ6+Pi/4fPlEDcUZr8FCcM8XaszorXmQMWBpvGRoh3kVeQBEOAXwLCYYVw/9Hpj2m/cKGKC+2arTAjReRIs7lKy1xigP/oNnHsLXPwwBPS+Lp96ez27Sna1mK3VeH3zyMBIRsWNck37HR47XBYeCiFakWA5U1rDtrXw3r3gHwQz18FZl3m6Vp3WUbdWeng6E/tPdA2yZ4RnyDbuQohTkmA5E7Vl8PYC2P0WZFwAV/4Fwvt5ulbt6ky31qyhs1yr2aODoj1bYSFEr+RTweLWWWH5/4Y3boGqozDlQfjxAq+7uqN0awkhPEFmhXWV3Qaf/hk+/R+ISoefrYTkMW6t3+k6VbfWqPhR0q0lhDhtMiusO5Tlw5u3wKEvYOR1cNmfIdDikapIt5YQwltJsHTFB/dD0R5jn68RP+/Rp7barew6bnRrbSvaxo4i6dYSQngnCZauuOx/wFZndIF1s8Yt47cd28b2ou3sLNnp2lsrLTyNC1MuJCshS7q1hBBeR4KlKyzdt1X70eqjRmvEGSQ/lP2ARmNSJoZGD+UXQ37huv5IbHBst9VDCCHOlASLBzi0g9wTuWwr2sa2om1sP7adwupCAEL8QxgZN5LbRt1GVnwWI2JHyJbxQohepdcHi1LKD3gYCAe2aq3/6uEqtdK4SWNja2R70XYqGioAiA2OZXT8aH457JdkJWQxOGow/n69/scihPBhHv0LppRaDUwDirTWZzc7finwFGACVmqtl3RwmulAMlAKFHRjdTutoqHCNVtr27FtfFvyLQ2OBgAyIjKYmjaV0fGjyYrPIsWSIuMjQog+xdP/GucAy4G1jQeUUibgWWAqRlBsUUptwAiZxSc9fi4wBPhca/0XpdTrwMc9UO8WjlQdMbq0nDO29pXtQ6PxV/4Mix3GdUOvc60fkWm/Qoi+zqPBorX+VCmVftLhscA+rXUugFLqVWC61noxRuumBaVUAdDgvGlv77mUUvOAeQCpqamnXWeHdrDvxD62HdvmCpOj1UcBCA0IZVTcKC5Ju4SshCy59ogQwid5usXSlmTgULPbBcC4Dsq/CTyjlJoAfNpeIa31CmAFGCvvT6diT297mle/e5VKayUA8cHxZCVkMXr4aLISshgUOUguqSuE8HneGCxtDTi0GwRa6xrgpk6d+Az3CksISeCSjEvIis8iKyGLfqH9ZHxECCFO4o3BUgD0b3Y7BSj0UF1auOasazxdBSGE8HretR2vYQswSCmVoZQyAzOBDe44sdb6ba31vIiICHecTgghRBs8GixKqXXA58AQpVSBUuomrbUNuAN4H9gDvKa13uWm57tcKbWivLzcHacTQgjRBtk2XwghRKd0dtt8b+wKE0II0Yv5VLBIV5gQQnQ/nwoWGbwXQoju51PBIi0WIYTofj4VLNJiEUKI7ueTs8KUUsVAvqfrcYZigRJPV8KLyOvRRF6LluT1aHKmr0Wa1jruVIV8Mlj6AqXU1s5M+/MV8no0kdeiJXk9mvTUa+FTXWFCCCG6nwSLEEIIt5Jg6b1WeLoCXkZejybyWrQkr0eTHnktZIxFCCGEW0mLRQghhFtJsAghhHArCRYvp5Tqr5TapJTao5TapZRa4DwerZT6UCm11/k5ytN17UlKKZNSartS6h3n7Qyl1BfO1+N/ndfy8QlKqUil1OtKqe+c75Mf+er7Qyn1G+fvybdKqXVKqSBfem8opVYrpYqUUt82O9bme0EZnlZK7VNKfaOUynJXPSRYvJ8NuEdrPRQYD/xKKTUMuBf4WGs9CPjYeduXLMC4Xk+jpcCTztejjE5errqPeAp4T2t9FjAS43XxufeHUioZ+DWQrbU+GzBhXCjQl94bOcClJx1r773wE2CQ82Me8Ly7KiHB4uW01ke01tucX1di/NFIBqYDf3UW+yswwzM17HlKqRTgp8BK520FTAJedxbxmddDKRUOXACsAtBaN2itT+C77w9/IFgp5Q+EAEfwofeG1vpToPSkw8KJT24AAAbzSURBVO29F6YDa7XhP0CkUirJHfWQYOlFlFLpwGjgCyBBa30EjPAB4j1Xsx63DFgIOJy3Y4ATzquPAhRghK8vyASKgTXOrsGVSqlQfPD9obU+DDwGHMQIlHLgK3z3vdGovfdCMnCoWTm3vTYSLL2EUioMeAO4S2td4en6eIpSahpQpLX+qvnhNor6yjx6fyALeF5rPRqoxge6vdriHDuYDmQA/YBQjO6ek/nKe+NUuu33RoKlF1BKBWCEysta6zedh481Nludn4s8Vb8edh5whVIqD3gVo5tjGUYz3t9ZJgUo9Ez1elwBUKC1/sJ5+3WMoPHF98cU4IDWulhrbQXeBH6M7743GrX3XigA+jcr57bXRoLFyznHD1YBe7TWTzS7awNwg/PrG4C3erpunqC1vk9rnaK1TscYmN2otb4e2AT83FnMl16Po8AhpdQQ56HJwG588/1xEBivlApx/t40vhY++d5opr33wgZgtnN22HigvLHL7EzJynsvp5Q6H/gnsJOmMYVFGOMsrwGpGL9QV2utTx6069OUUhOB32qtpymlMjFaMNHAdmCW1rrek/XrKUqpURgTGcxALjAH459Gn3t/KKUeAq7BmE25HbgZY9zAJ94bSql1wESM7fGPAQ8Af6eN94IzfJdjzCKrAeZorbe6pR4SLEIIIdxJusKEEEK4lQSLEEIIt5JgEUII4VYSLEIIIdxKgkUIIYRbSbAI4aSUulEppZ3TmPsUpVSeUmqzp+shfIMEi+iTlFITnSHR+GFXSpU5t1P/q1LqUuc8fiGEm8k6FtEnOVsdm4B1wLsY+yJZgCEYu7umAh9hLBY74XyMCQgAGrTWjjZO22sppQKB/9/euYZYVUVx/PcnESMdpiiM+tCDkCmjD5LQh8DI7GFQ2oegkBCEwkAoawZJelmQ9KIM61MWNRQU4WMwmoyY0kErkMqChprJIJyIahQde9C4+rD3jdPh3pk7wx21uf8fbPY5666z9zp34K7ZZ62zV0TEXyfaFjP1mTa2ijH/a/ZGRGdRIGk18CSwmuR4bgCIiBFg5LhbeByYqm+am5MTPwozTUdEjETEfcAu4Pq8bU7VGEtBtlDSQ5J+kPR7rkh4RdZZIGmXpGFJg5IerDavpMslbZb0i6Q/JfVJWlvYILGi15NjIufkKohDeexuSXNKujMkPZLHOirpoKR9kp4q6VWNsUhaIqlX0pHceiXdXEVvf7arTdJ2SYclHVKqXHl2vd+9aQ7sWEwz83Lub6xDdz3pEdrzwKOkOijdkpaQdtHdCdwPfAOsk7SseLGkxUAvMAd4hlTpcDewjrRqKnMa8DFpBfUAsJG0B9TW/MiuwkbSflB7SCuwtaQqgVePdUOS7gY2k/bQehx4LB9vkXRnlUvOBXpI+021A28AtwCvjTWXaTIiws1tyjXSj3CQNqmspTMv67yTz5fn86sKOhXZXmB6QX5Tlv8NzC/Ip5OKTO0uyGYAP5EcxbSSDfdWmbMnyzpKuu1Zfl1B9hvwbh3fx36gp3B+OnAE+A5oKchbgH7gMNBauj6AW0vjbszythP9N3c7eZpXLKaZqRRMa6lD96X4b+B7Z+73RMRnFWHW+ZRUR7zCImA28AqpNsiZlUZKLAC4tjTfMWBDSfZh7otjHwLmSrq0jnsosoi0KtoQhcJx+fgFYCapvkmRAxHxVg2bLhrn/GYK4+C9aWYqDqWeipwDxZOIGMrZyt9X0R0ilUuucHHuN40y/uzS+YGI+KMk+zX3xbHvAV4H9kkaIGXCdQFdMXpm2wW5/7rKZ1/l/sKSfKCsWMMm0+TYsZhm5rLc99WhWytbrJ4sssr7Mu3A5zV0ypX7Rhv33/dvImKrpPOBxcAC0ipjBbBT0jVRO714Iu/w1GWTMXYspplZkfvtkzzPt7kfjogPGj14pAJenUBnfulzPdBBqv/+do3L+nM/lxTsL3JJ7qutUIwZE8dYTNMh6RRJTwNXkgLfvZM8ZTepzvgaSWdUsedUSbPGO2i+j9aiLCKCVCURUoZXLXYAw8Cq4tz5eBUpsL9jvDYZA16xmKnPvELqb/HN+/OA94HbJ9uAiBiWdAepRGyfpE2kbKxWoI2UsruUlA02HmYBg5K2kZzJz6TYyUpSnKdrFJsOSuogZXV9IunV/NFyUiD+rog4NE57jAHsWMzU57bcjpH+C/8R+Ah4MyLeO15GRES3pPnAGmAZcBbpx78feBb4cgLDHgWeAxaSYiszSanO24AnIqIctynb9KKkQVLs5+Es/gJYGhFbJmCPMYD3CjPGGNNgHGMxxhjTUOxYjDHGNBQ7FmOMMQ3FjsUYY0xDsWMxxhjTUOxYjDHGNBQ7FmOMMQ3FjsUYY0xDsWMxxhjTUP4BsYURiWfI4nwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(dim_range, time_range_matmul, label=\"Matmul\")\n",
    "plt.plot(dim_range, time_range_numba_matmul, label=\"Matmul Numba\")\n",
    "plt.plot(dim_range, time_range_np, label=\"Numpy\")\n",
    "plt.legend(fontsize=18)\n",
    "plt.xlabel(\"Dimension\", fontsize=18)\n",
    "plt.ylabel(\"Time\", fontsize=18)\n",
    "plt.yscale(\"log\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Почему наивная реализация такая медленная?\n",
    "\n",
    "Это так по двум причинам:\n",
    "\n",
    "- Она не использует преимуществ быстрой памяти (кэшей) и в целом архитектуру памяти\n",
    "- Она не использует параллелизм (особенно важно для GPU)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Архитектура памяти  \n",
    "<img width=80% src=\"./Memory-Hierarchy.jpg\">\n",
    "\n",
    "- \"Быстрой\" памяти мало\n",
    "- \"Медленной\" памяти много"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Реализация алгоритмов наиболее вычислительно интенсивным образом\n",
    "\n",
    "<font color='red'>**Эффективные реализации для задач NLA**</font> используют блочные версии алгоритмов. <br>\n",
    "\n",
    "Эта идея лежит в основе **[BLAS (Basic Linear Algebra Subroutines)](http://www.netlib.org/blas/)**, написанной на Фортране много лет назад, однако до сих пор не потерявшей актуальность. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Разделим матрицу на блоки! Для примера рассмотрим разбиение на блочную матрицу размера $2 \\times 2$:\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",
    "Тогда  \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": "slide"
    }
   },
   "source": [
    "Если $A_{11}, B_{11}$ и их произведение помещается в кэш-память (которая достигает 12 Мб (L3) в [последнем CPU от Intel](https://en.wikipedia.org/wiki/List_of_Intel_microprocessors#Desktop)), тогда мы только один раз загружаем их в память.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## BLAS\n",
    "\n",
    "BLAS имеет три уровня:\n",
    "1. BLAS-1,  операции типа $c = a + b$\n",
    "2. BLAS-2, операции типа умножения матрицы на вектор\n",
    "3. BLAS-3, произведение матриц\n",
    "\n",
    "**Q:** Какова принципиальная разница между ними?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**A:** основная разница в соотношении между количеством операций и количеством данным!\n",
    "\n",
    "1. BLAS-1: $\\mathcal{O}(n)$ данных, $\\mathcal{O}(n)$ операций\n",
    "2. BLAS-2: $\\mathcal{O}(n^2)$ данных, $\\mathcal{O}(n^2)$ операций\n",
    "3. BLAS-3: $\\mathcal{O}(n^2)$ данных, $\\mathcal{O}(n^3)$ операций"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Почему BLAS так важна и актуальна?\n",
    "\n",
    "1. Стандартные реализации базовых операций линейной алгебры\n",
    "2. Сохраняются стандартные имена для всех операций, сохраняющиеся в каждой новой реализации (например [ATLAS](https://www.netlib.org/atlas/), [OpenBLAS](https://www.openblas.net/), [MKL](https://software.intel.com/en-us/mkl)). Вы можете вызвать функцию умножения матриц (GEMM) в вашем коде, слинковать его с любой реализацией BLAS, и всё будет корректнро работать\n",
    "3. Позволяет формулировать новые алгоритмы в терминах BLAS операций – гарантии скорости\n",
    "4. Есть обёртки для большинства популярных языков"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Реализации BLAS и их использование\n",
    "\n",
    "1. [ATLAS](http://math-atlas.sourceforge.net) - Automatic Tuned Linear Algebra Software. Она автоматически адаптируется под архитектуру конкретного компьютера.\n",
    "2. [LAPACK](http://www.netlib.org/lapack/) - Linear Algebra Package. Она выполняет высокоуровневые операции (например матричные разложения), которые опираются на вызовы функций из BLAS.\n",
    "3. [Intel MKL](https://software.intel.com/en-us/intel-mkl) - Math Kernel Library. Эта реализация является оптимизированной версией BLAS и LAPACK для процессоров Intel. Доступна в Anaconda: \n",
    "```\n",
    "conda install mkl\n",
    "```\n",
    "MATLAB использует Intel MKL по умолчанию.\n",
    "\n",
    "4. OpenBLAS – это оптимизированная версия BLAS, основанная на [GotoBLAS](https://en.wikipedia.org/wiki/GotoBLAS). \n",
    "\n",
    "5. PyTorch [поддерживает](https://pytorch.org/docs/stable/torch.html#blas-and-lapack-operations) некоторые вызовы из BLAS и LAPACK\n",
    "\n",
    "6. Для GPU реализован [cuBLAS](https://docs.nvidia.com/cuda/cublas/index.html).\n",
    "\n",
    "Для сравнения производительности OpenBLAS и Intel MKL, см [этот обзор](https://software.intel.com/en-us/articles/performance-comparison-of-openblas-and-intel-math-kernel-library-in-r)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Более быстрые алгоритмы для умножения матриц\n",
    "\n",
    "Напомним, что сложность алгоритма матричного умножения $\\mathcal{O}(n^3)$ операций.\n",
    "Однако, данных требуется только $\\mathcal{O}(n^2)$.\n",
    "\n",
    "**Q:** возможно ли уменьшить число операций необходимое для умножения матриц до $\\mathcal{O}(n^2)$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**A**: алгоритм умножения матриц со сложностью $\\mathcal{O}(n^2)$ до сих пор не найден."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "* Метод Штрассена – $\\mathcal{O}(n^{2.807\\dots})$, иногда используется на практике\n",
    "\n",
    "* [Текущий мировой рекорд](http://arxiv.org/pdf/1401.7714v1.pdf) $\\mathcal{O}(n^{2.37\\dots})$ – большая константа, не практичный, основан на [алгоритме Coppersmith-Winograd'a](https://en.wikipedia.org/wiki/Coppersmith%E2%80%93Winograd_algorithm).\n",
    "\n",
    "Рассмотрим метод Штрассена более детально."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Наивный алгоритм\n",
    "\n",
    "Пусть $A$ и $B$ две матрицы $2\\times 2$. Наивное умножение $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",
    "состоит из $8$ умножений и $4$ сложений."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Алгоритм Штрассена\n",
    "\n",
    "В работе [Gaussian elimination is not optimal](http://link.springer.com/article/10.1007%2FBF02165411?LI=true) (1969 г.) Штрассен обнаружил, что можно вычислить матрицу $C$, используя 18 сложений и только 7 умножений:\n",
    "\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",
    "где\n",
    "\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",
    "К счастью, эти формулы работают даже если $a_{ij}$ и $b_{ij}$, $i,j=1,2$ – блочные матрицы.\n",
    "\n",
    "Таким образом, схема алгоритма Штрассена такова \n",
    "- Сначала мы  <font color='red'>разбиваем</font> матрицы $A$ и $B$ размера $n\\times n$, $n=2^d$ <font color='red'> на 4 блока</font> размера $\\frac{n}{2}\\times \\frac{n}{2}$\n",
    "- После чего <font color='red'>вычисляем произведения</font> по указанным выше формулам <font color='red'>рекурсивно</font>\n",
    "\n",
    "Это приводит нас к идее **разделяй и властвуй**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Сложность метода Штрассена\n",
    "\n",
    "#### Число умножений\n",
    "\n",
    "Обозначим $M(n)$ число умножений необходимых для вычисления произведения двух матриц $n\\times n$ с помощью стратегии \"разделяй и властвуй\".\n",
    "Тогда для наивного алгоритма получаем следующее число умножений\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",
    "Итак, даже используя стратегию \"разделяй и властвуй\" мы не улучшаем сложность в $\\mathcal{O}(n^3)$ операций.\n",
    "\n",
    "Вычислим число умножений в методе Штрассена:\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": [
    "#### Число сложений\n",
    "\n",
    "Для алгоритма Штрассена:\n",
    "\n",
    "$$\n",
    "A_\\text{strassen}(n) = 7 A_\\text{strassen}\\left( \\frac{n}{2} \\right) + 18 \\left( \\frac{n}{2} \\right)^2\n",
    "$$\n",
    "поскольку на первом уровне нам нужно сложить матрицы размера $\\frac{n}{2}\\times \\frac{n}{2}$ 18 раз, а потом идти на следующий уровень рекурсии для каждого из 7 умножений. \n",
    "Таким образом,\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",
    "(так как $4^d = n^2$ и $7^d = n^{\\log_2 7}$).\n",
    "\n",
    "\n",
    "Асимптотику для $A(n)$ также можно найти с помощью следующей [теоремы](https://en.wikipedia.org/wiki/Master_theorem)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Итоговая сложность метода Штрассена\n",
    "\n",
    "Итоговая сложность –  $M_\\text{strassen}(n) + A_\\text{strassen}(n)=$ <font color='red'>$7 n^{\\log_2 7}$</font>. \n",
    "Метод Штрассена становится быстрее наивного алгоритма, если\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "2n^3 &> 7 n^{\\log_2 7}, \\\\\n",
    "n &> 667,\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "поэтому не самая хорошая идея спускаться на самый низкий уровень рекурсии."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Является ли метод Штрассена практически важным?\n",
    "\n",
    "- Долгое время считалось, что метод Штрассена имеет только теоретический интерес и плохо применим на практике\n",
    "- Статья [Strassen algorithm reloaded](http://jianyuhuang.com/papers/sc16.pdf)\n",
    "опровергает утверждение, что метод Штрассена не практичен, и предлагает реализацию, которая работает не только на больших матрицах."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Алгоритм Штрассена и тензорный ранг\n",
    "\n",
    "- Не очень понятно, как Штрасен вывел формулы для вычисления произведения матриц. \n",
    "- Однако сейчас мы увидим, что они не такие искусственные как может показаться.\n",
    "- Чтобы это показать воспользуемся общим подходом под названием **тензорные разложения**.\n",
    "- Под тензором мы понимаем многомерный массив – обобщение матрицы на большие размерности.\n",
    "\n",
    "Давайте пронумеруем элементы в матрице $2\\times 2$ следующим образом \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",
    "Это можно записать как\n",
    "\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}$ – трёхмерный массив, состоящий из нулей и единиц:\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",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Трилинейное разложение\n",
    "\n",
    "Для получения метода Штрассена сделаем следующий трюк – разложим $x_{ijk}$\n",
    "\n",
    "$$\n",
    "x_{ijk} = \\sum_{\\alpha=1}^r u_{i\\alpha} v_{j\\alpha} w_{k\\alpha}.\n",
    "$$\n",
    "\n",
    "Это разложение называется **трилинейным тензорным разложением** и означает разделение переменных: у нас есть $r$ слагаемых с разделёнными индексами $i$, $j$ и $k$. Число $r$ называется **тензорным рангом**. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Метод Штрассена как трилинейное разложение\n",
    "\n",
    "Теперь у нас есть\n",
    "\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",
    "\n",
    "Умножение на $u_{i\\alpha}$ или $v_{j\\alpha}$ или $w_{k\\alpha}$ не требует рекурсии, поскольку $u, v$ и $w$ известные предвычисленные матрицы. \n",
    "Поэтому у нас есть только $r$ умножений $\\left(\\sum_{i=1}^4  u_{i\\alpha} a_i \\right)$ $\\left( \\sum_{j=1}^4 v_{j\\alpha} b_j\\right)$, где оба множителя зависят от входных данных.\n",
    " \n",
    "Как вы можете догадаться ранг тензора $x_{ijk}$ равен 7, что приводит нас к методу Штрассена с 7 умножениями!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Резюме первой части лекции\n",
    "\n",
    "- Умножение матриц – одна из основных операций в вычислительной линейной алгебре. Если необходима высокая эффективность, необходимо использовать блочные алгоритмы\n",
    "- Это связано с оптимальным использованием иерархии памяти\n",
    "- Концепция блочных алгоритмов\n",
    "- Алгоритм Штрассена и трилинейные формы"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Часть 2. Матричные разложения и методы их вычисления"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## О каких разложениях пойдёт речь далее\n",
    "\n",
    "- LU разложение (краткое напоминание)\n",
    "- QR разложение\n",
    "- Разложение Шура\n",
    "- Несколько слов про вычисление SVD – содержательно на следующей лекции"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Общая концепция матричных разложений\n",
    "\n",
    "Вычислительная линейная алгебра занимается решением следующих задач:\n",
    "\n",
    "- Решение линейных систем $Ax = f$\n",
    "- Вычисление собственных значений и векторов\n",
    "- Вычисление сингулярных значений и векторов \n",
    "- Вычисление обратных матриц и иногда детерминантов\n",
    "- Вычисление **матричных функций** таких как $\\exp(A), \\cos(A)$ (это не поэлементные функции!)\n",
    "\n",
    "Для решения таких задач мы представляемм матрицу в виде суммы и/или произведения матриц **более простой структуры**, таких что мы можем решить эти задачи быстрее и/или более устойчивым образом.\n",
    "\n",
    "Что такое **более простая структура**?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Матрицы простой структуры\n",
    "\n",
    "Мы уже упоминали некоторые классы структурированных матриц. \n",
    "\n",
    "Для плотных матриц это следующие матрицы\n",
    "\n",
    "- **унитарные матрицы**\n",
    "- **нижне- верхнетреугольные матрицы** \n",
    "- **диагональные матрицы**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## PLU разложение\n",
    "\n",
    "Любая невырожденная матрица моджет быть разложена как\n",
    "\n",
    "$$\n",
    "A = P L U,\n",
    "$$\n",
    "\n",
    "где $P$ – матрица перестановки, $L$ – **нижнетреугольная матрица**, $U$ – **верхнетреугольная матрица**.\n",
    "\n",
    "**Основная цель** вычисления LU разложения – это решение систем линейных уравнений, потому что\n",
    "\n",
    "$$\n",
    "    A^{-1} f = (L U)^{-1} f = U^{-1} L^{-1} f, \n",
    "$$\n",
    "\n",
    "что сводится к решению двух линейных систем \n",
    "\n",
    "$$\n",
    "     L y = f,  \\quad U x = y\n",
    "$$\n",
    "\n",
    "с нижне- и верхнетреугольной матрицей, соответственно. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Положительно определённые матрицы и разложение Холецкого\n",
    "\n",
    "Если матрица эрмитова и положительно определена, то есть\n",
    "\n",
    "$$\n",
    "A = A^*, \\quad (Ax, x) > 0, \\quad x \\ne 0,\n",
    "$$\n",
    "\n",
    "тогда она может быть представлена в виде\n",
    "\n",
    "$$\n",
    "A = RR^*,\n",
    "$$\n",
    "\n",
    "где $R$ – нижнетреугольная матрица.\n",
    "\n",
    "Это нам понадобится для QR разложения."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR разложение\n",
    "\n",
    "Как следует из названия, это представление матрицы в виде произведения\n",
    "\n",
    "$$\n",
    "    A = Q R, \n",
    "$$\n",
    "\n",
    "где $Q$ – матрица с **ортогональными столбцами** и $R$ – **верхнетреугольная**.  \n",
    "\n",
    "Размеры матриц: $Q$ – $n \\times m$, $R$ – $m \\times m$, если $n\\geq m$.\n",
    "\n",
    "QR разложение определено для любой **прямоугольной матрицы**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Приложения QR разложения\n",
    "\n",
    "Это разложение играет ключевую роль при решении многих задач, например:\n",
    "- Получение ортогонального базиса в линейном пространстве\n",
    "- Используется для препроцессинга при вычислении SVD\n",
    "- QR алгоритм для вычисления собственных векторов и собственных значений ([один из 10 самых важных алгоритмов ХХ века](https://archive.siam.org/pdf/news/637.pdf)) основан на вычислении QR разложения\n",
    "- Решение переопределённых систем линейных уравнений"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR разложение и линейная задача наименьших квадратов: напоминание\n",
    "\n",
    "Пусть нам надо решить переопределённую систему $Ax = f$ в смысле наименьших квадратов невязок  \n",
    "\n",
    "$$\n",
    "\\Vert A x - f \\Vert_2 \\rightarrow \\min_x,\n",
    "$$\n",
    "\n",
    "где $A$ – $n \\times m$ и $n \\geq m$.\n",
    "\n",
    "Тогда мы разложим её \n",
    "\n",
    "$$\n",
    "A = Q R,\n",
    "$$\n",
    "\n",
    "и используем выражение для псевдообратной матрицы в случае матрицы полного ранга:\n",
    "\n",
    "$$\n",
    "x = A^{\\dagger}b = (A^*A)^{-1}A^*b = ((QR)^*(QR))^{-1}(QR)^*b = (R^*Q^*QR)^{-1}R^*Q^*b = R^{-1}Q^*b. \n",
    "$$\n",
    "\n",
    "тогда $x$ может быть найдён из решения следующей системы\n",
    "\n",
    "$$R x = Q^*b$$\n",
    "\n",
    "Заметим, что это квадратная система линейных уравнений с верхнетреугольной матрицей. Какова сложность решения такой системы?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Существование QR разложения\n",
    "\n",
    "**Теорема.**\n",
    "\n",
    "Каждая матрица $n \\times m$ может быть представлена в виде QR разложения. \n",
    "\n",
    "\n",
    "Существует несколько способов доказательства и вычисления:\n",
    "\n",
    "- Теоретический: используя матрицу Грама и разложение Холецкого\n",
    "- Геометрический: используя ортогонализацию Грама-Шмидта\n",
    "- Практический: используя преобразования Гивенса/Хаусхолдера"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Доказательство с использованием разложения Холецкого\n",
    "\n",
    "Если у нас есть разложение\n",
    "\n",
    "$$A = QR,$$\n",
    "\n",
    "тогда $A^* A = ( Q R)^* (QR)  = R^* (Q^* Q) R = R^* R$, матрица $A^* A$ называется **матрицей Грама**, и её элементы – скалярные произведения столбцов матрицы $A$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Случай матрицы полного ранга\n",
    "\n",
    "Пусть $A$ имеет **полный столбцовый ранг**. Тогда легко показать, что $A^* A$ положительно определена:\n",
    "\n",
    "$$\n",
    "   (A^* A y, y) = (Ay, Ay) = \\Vert Ay \\Vert^2  > 0, \\quad y\\not = 0.\n",
    "$$\n",
    "\n",
    "Поэтому, $A^* A = R^* R$ всегда существует.\n",
    "\n",
    "Тогда матрица $A R^{-1}$ унитарна:  \n",
    "\n",
    "$$\n",
    "   (A R^{-1})^* (AR^{-1})= R^{-*} A^* A R^{-1} = R^{-*} R^* R R^{-1} = I.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Случай матрицы неполного ранга\n",
    "\n",
    "- QR разложение по-прежнему существует.\n",
    "\n",
    "- Для любой матрицы неполного ранга существует последовательность матриц полного ранга $A_k$ такая что $A_k \\rightarrow A$ (почему?).\n",
    "\n",
    "- Каждая $A_k$ может быть разложена $A_k = Q_k R_k$. \n",
    "\n",
    "- Множество унитарных матриц образует компакт, таким образом найдётся сходящаяся подпоследовательность $Q_{n_k} \\rightarrow Q$ (почему?), и $Q^* A_k \\rightarrow Q^* A = R$, которая верхнетреугольная."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Устойчивость вычисления QR разложения с помощью разложения Холецкого \n",
    "\n",
    "Итак, простейший способ вычисления QR разложения следующий\n",
    "\n",
    "$$A^* A = R^* R,$$\n",
    "\n",
    "и\n",
    "\n",
    "$$Q = A R^{-1}.$$\n",
    "\n",
    "Это **плохая идея** с точки зрения устойчивости. Покажем, почему это так."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Built-in QR orth 5.061169633223855e-16\n",
      "Via Cholesky: 3.2416882870748012e-12\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "n = 50\n",
    "r = 4\n",
    "a = [[1.0 / (i + j + 0.5) for i in range(r)] for j in range(n)]\n",
    "a = np.array(a)\n",
    "q, Rmat = np.linalg.qr(a)\n",
    "e = np.eye(r)\n",
    "print('Built-in QR orth', np.linalg.norm(np.dot(q.T, q) - e))\n",
    "gram_matrix = a.T.dot(a)\n",
    "Rmat1 = np.linalg.cholesky(gram_matrix)\n",
    "q1 = np.dot(a, np.linalg.inv(Rmat1.T))\n",
    "print('Via Cholesky:', np.linalg.norm(np.dot(q1.T, q1) - e))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Второй способ: ортогонализация Грама-Шмидта\n",
    "\n",
    "- QR разложение – это способ записи процесса ортогонализации Грама-Шмидта\n",
    "- Для данного набора векторов $a_1, \\ldots, a_m$ мы хотим найти ортонормированный базис $q_1, \\ldots, q_m$ такой чтобы каждый вектор $a_i$ представлялся как линейная комбинация векторов из базиса.  \n",
    "\n",
    "**Метод Грама-Шмидта:**\n",
    "1. $q_1 := a_1/\\Vert a_1 \\Vert$\n",
    "2. $q_2 := a_2 - (a_2, q_1) q_1, \\quad q_2 := q_2/\\Vert q_2 \\Vert$\n",
    "3. $q_3 := a_3 - (a_3, q_1) q_1 - (a_3, q_2) q_2, \\quad q_3 := q_3/\\Vert q_3 \\Vert$\n",
    "4. И так далее  \n",
    "\n",
    "Заметим, что преобразование из $A$ в $Q$ имеет треугольную структуру, поскольку мы вычитаем из $k$-го вектора только предыдущие. Это следует из того, что произведение треугольных матриц – это треугольная матрица."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Модифицированный метод Грама-Шмидта\n",
    "\n",
    "- Метод Грама-Шмидта может быть очень неустойчивым (то есть генерировать векторы, которая не являются ортогональными, особенно если $q_k$ малой нормы).  \n",
    "Это называется **потеря ортогональности**.  \n",
    "\n",
    "- Этот недостаток метода Грама-Шмидта исправляется с помощью **модифицированного метода Грама-Шмидта**. Вместо вычисления\n",
    "\n",
    "$$q_k := a_k - (a_k, q_1) q_1 - \\ldots - (a_k, q_{k-1}) q_{k-1}$$\n",
    "\n",
    "мы будем вычислять это выражение шаг за шагом. Сначала присвоим $q_k := a_k$, затем последовательно ортогонализуем:\n",
    "\n",
    "$$\n",
    "   q_k := q_k - (q_k, q_1)q_1, \\quad q_k := q_{k} - (q_k,q_2)q_2, \\ldots\n",
    "$$\n",
    "\n",
    "- В точной арифметике, это одинаковые алгоритмы. В неточной арифметике они абсолютно разные!\n",
    "\n",
    "- Заметим, что сложность по-прежнему $\\mathcal{O}(nm^2)$ операций"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR разложение: почти практический способ\n",
    "\n",
    "Если $A = QR$, тогда  \n",
    "\n",
    "$$\n",
    "R = Q^* A,\n",
    "$$\n",
    "\n",
    "и нам нужно найти такую ортогональную матрицу $Q$, которая преобразует данную матрицу $A$ в верхнетреугольную.  \n",
    "Для простоты мы будем смотреть на матрицы $n \\times n$ такие что\n",
    "\n",
    "$$ Q^* A = \\begin{bmatrix} * & * & *  \\\\ 0 & * & * \\\\ 0 & 0 & * \\\\ & 0_{(n-m) \\times m} \\end{bmatrix} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Будем приводить матрицу к такому виду столбец за столбцом.\n",
    "\n",
    "Сначала найдём такую матрицу Хаусхолдера $H_1 = (I - 2 uu^{\\top})$ что \n",
    "\n",
    "$$ H_1 A = \\begin{bmatrix} * & * & * \\\\ 0 & * & * \\\\ 0 & * & * \\\\ 0 & * & * \\end{bmatrix} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Затем\n",
    "\n",
    "$$ H_2 H_1 A = \\begin{bmatrix} * & * & * \\\\ 0 & * & * \\\\ 0 & 0 & * \\\\ 0 & 0 & * \\end{bmatrix}, $$\n",
    "\n",
    "где\n",
    "\n",
    "$$ H_2 = \\begin{bmatrix} 1 & 0 \\\\ 0 & H'_2, \\end{bmatrix} $$\n",
    "\n",
    "и $H'_2$ матрица Хаусхолдера $3 \\times 3$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "И наконец, \n",
    "\n",
    "$$ H_3 H_2 H_1 A = \\begin{bmatrix} * & * & * \\\\ 0 & * & * \\\\ 0 & 0 & * \\\\ 0 & 0 & 0 \\end{bmatrix}, $$\n",
    "\n",
    "где $H_3=\\begin{bmatrix}I_2 & \\\\ & {\\widetilde H}_3 \\end{bmatrix}$ такая что\n",
    "\n",
    "$$ \n",
    "{\\widetilde H}_3 \\begin{bmatrix} \\boldsymbol{\\times}\\\\ \\boldsymbol{\\times} \\\\ \\boldsymbol{\\times}  \\end{bmatrix} = \n",
    "\\begin{bmatrix} \\times \\\\ 0 \\\\ 0  \\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "Попробуйте самостоятельно реализовать такой алгоритм, это просто!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Получение QR разложения\n",
    "\n",
    "Так как \n",
    "\n",
    "$$ H_3H_2H_1A = HA = R,$$\n",
    "\n",
    "где $H$ – унитарная матрица, то\n",
    "\n",
    "$$ A = H^*R. $$\n",
    "\n",
    "Таким образом $Q = H^*$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR разложение: практический способ\n",
    "\n",
    "- Поскольку мы работаем с плотной матрицей, то на практике нам нужен алгоритм, оперирующий блоками (почему?).  \n",
    "\n",
    "- Вместо использования преобразования Хаусхолдера, мы будем использовать **блочное преобразование Хаусхолдера** вида \n",
    "\n",
    "$$H = (I - 2UU^*), $$\n",
    "\n",
    "где $U^* U = I$.\n",
    "\n",
    "- Это позволяет использовать BLAS-3 операции."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR разложение: практический способ - 2\n",
    "\n",
    "Аналогично преобразованию Хаусхолдера для вычисления QR разложения можно использовать преобразования Гивенса\n",
    "\n",
    "$$\\begin{bmatrix} \\times & \\times & \\times \\\\ \\bf{*} & \\times & \\times \\\\ \\bf{*} & \\times & \\times \\end{bmatrix} \\to \\begin{bmatrix} * & \\times & \\times \\\\ * & \\times & \\times \\\\ 0 & \\times & \\times \\end{bmatrix} \\to \\begin{bmatrix} \\times & \\times & \\times \\\\ 0 & * & \\times \\\\ 0 & * & \\times \\end{bmatrix} \\to \\begin{bmatrix} \\times & \\times & \\times \\\\ 0 & \\times & \\times \\\\ 0 & 0 & \\times \\end{bmatrix} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Преобразование Гивенса vs. преобразование Хаусхолдера\n",
    "\n",
    "- Матрицы Хаусхолдера полезны для плотных матриц (сложность примерно в два раза меньше), в которых необходимо занулить большое число элементов.\n",
    "- Вращения Гивенса больше подходят для работы с разреженными матрицами или параллельными вычислениями, поскольку они зануляют только один элемент каждым действием."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Анонс следующей лекции (14 октября)\n",
    "\n",
    "- Собственные векторы и собственные значения\n",
    "- Какие задачи сводятся к их вычислению\n",
    "- Степенной метод\n",
    "- Теорема Гершгорина\n",
    "- Разложение Шура и алгоритм его вычисления"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## На следующей неделе будет два семинара!"
   ]
  },
  {
   "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": {
  "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}