seminars/seminar6/seminar6.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Семинар 6.\n",
    "## Часть 1. Повторение пройденного"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Что мы узнали по прошествии половины курса?\n",
    "\n",
    "- Какие форматы представления чисел в компьютере существуют и как устроена память в компьютере\n",
    "- Основные операции с матрицами и векторами, сложность их вычисления и частные случаи\n",
    "- Матричные разложения: SVD, скелетное разложение, LU/Холецкий, QR, разложение Шура\n",
    "- Приложения SVD\n",
    "- Вычисление собственных значений и собственных векторов: QR алгоритм и степенной метод\n",
    "- Алгоритмы вычисления SVD (узнаем в следующий раз)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Singular value decomposition (SVD)\n",
    "\n",
    "- Представление матрицы в виде\n",
    "\n",
    "$$ A = U\\Sigma V^*,$$\n",
    "\n",
    "- $U$ и $V$ – унитарные\n",
    "- $\\Sigma$ – диагональная\n",
    "\n",
    "- Существует для **любой** матрицы\n",
    "- Позволяет получить наилучшее малоранговое приближение матрицы в спектральной и Фробениусовой нормах (теорема Эккарта-Янга)\n",
    "- Как оценивается ошибка аппроксимации для этих норм?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Задача 1: решение системы линейных уравнений\n",
    "\n",
    "$$ Ax = b $$\n",
    "\n",
    "- LU разложение/разложение Холецкого\n",
    "- Задача наименьших квадратов и compressed sensing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "###  Случай 1: матрица квадратная\n",
    "\n",
    "Решение системы линейных уравнений с квадратной матрицей $A$\n",
    "\n",
    "$$A x = b$$\n",
    "\n",
    "существует тогда и только тогда, когда матрица невырождена, то есть\n",
    "* $\\det A \\ne 0$\n",
    "* матрица $A$ имеет полный ранг\n",
    "* существует $A^{-1}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Как решать?\n",
    "\n",
    "- Метод Гаусса или LU разложение\n",
    "\n",
    "**Определение**: LU разложение матрицы $A$ – это представление\n",
    "\n",
    "$$A =  LU,$$\n",
    "\n",
    "где $L$ – **нижнетреугольная** и $U$ – **верхнетреугольная** матрица.\n",
    "\n",
    "- Проход вперёд и проход назад в методе Гаусса\n",
    "    - Проход вперёд эквивалентен\n",
    "\n",
    "    $$ L y = f.$$\n",
    "\n",
    "    - Аналогично для прохода назад\n",
    "    \n",
    "    $$ U x = y.$$\n",
    "\n",
    "Всегда ли существует $LU$ разложение?\n",
    "- Сложность $\\mathcal{O}(n^3)$\n",
    "- **Q:** можно ли уменьшить степень?\n",
    "- Проблемы с устойчивостью и строго регулярные матрицы"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Выбор ведущего элемента в строке (pivoting)\n",
    "\n",
    "- Уже для любой матрицы существует PLU разложение\n",
    "\n",
    "$$ A = PLU, $$\n",
    "\n",
    "где $P$ – матрица перестановки.\n",
    "\n",
    "- Какая сложность?\n",
    "- Как с помощью такого разложения решить линейную систему?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Если матрица эрмитова и положительно определённая\n",
    "\n",
    "Эрмитова положительно определённая матрица $A$ строго регулярна и имеет **разложение Холецкого** вида\n",
    "\n",
    "$$A = RR^*,$$\n",
    "\n",
    "где $R$ нижнетреугольная матрица.\n",
    "\n",
    "Часто матрица $R$ называется \"квадратным корнем\" матрицы $A$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Устойчивость решения линейной системы\n",
    "\n",
    "- Поведение алгоритмов в точной и машинной арифметике **принципиально различаются**!\n",
    "- Для решения линейных систем введено понятие **числа обусловленности** матрицы\n",
    "\n",
    "$$\\mathrm{cond}(A) = \\Vert A \\Vert \\Vert A^{-1} \\Vert.$$\n",
    "\n",
    "- С его помощью можно оценить\n",
    "\n",
    "$$ \\frac{\\Vert \\widehat{x} - x \\Vert}{\\Vert x \\Vert} \\leq \\mathrm{cond}(A) \\Big(\\frac{\\Vert\\Delta A\\Vert}{\\Vert A \\Vert} + \\frac{\\Vert \\Delta f \\Vert}{ \\Vert f \\Vert}\\Big),$$\n",
    "\n",
    "где $(A + \\Delta A) \\widehat{x} = f + \\Delta f.$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Концепция предобусловливателя\n",
    "\n",
    "- Вместо решения системы $Ax = b$ будем решать систему $PAx = Pb$\n",
    "- Но число обусловленности $PA$ сильно меньше, чем у матрицы $A$!\n",
    "- **Q:** какая матрица $P$ идеальная с точки зрения минимизации числа обусловленности?\n",
    "- Далее в курсе подробно обсудим конкретные примеры и следствия использования такого подхода"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Случай 2: матрица неквадратная\n",
    "\n",
    "- строк больше, чем столбцов – переопределённая система\n",
    "    - решение может не существовать\n",
    "- строк меньше, чем столбцов – недоопределённая система\n",
    "    - решение неединственное"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Задача наименьших квадратов\n",
    "\n",
    "$$\\Vert A x - b \\Vert^2_2 \\rightarrow \\min_x$$\n",
    "\n",
    "- Геометрическая интерпретация\n",
    "- Стандартная задача: нужно вызвать функцию из библиотеки"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Матрица Грама и нормальное уравнение\n",
    "\n",
    "- Из условия оптимальности первого порядка получим нормальное уравнение\n",
    "\n",
    "$$\n",
    "A^* A x = A^* b\n",
    "$$\n",
    "\n",
    "Матрица $A^* A$ называется **матрицей Грама**.\n",
    "\n",
    "Свойства нормального уравнения:\n",
    "\n",
    "- Число обусловленности матрицы $A^* A$ равно квадрату числа обусловленности матрицы $A$ (проверьте!).\n",
    "- Поэтому решать нормальное уравнение в таком виде – не самая хорошая идея!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Псевдообратная матрица\n",
    "\n",
    "- Если матрица $A$ неполного ранга, то матрица Грама – вырождена\n",
    "- Поэтому введём понятие псевдообратной матрицы  $A^{\\dagger}$\n",
    "\n",
    "$$x = A^{\\dagger} b.$$\n",
    "\n",
    "- Матрица \n",
    "\n",
    "$$A^{\\dagger} = \\lim_{\\alpha \\rightarrow 0}(\\alpha I + A^* A)^{-1} A^*$$ \n",
    "\n",
    "называется псевдообратной матрицей Мура-Пенроуза для матрицы $A$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Свойства псевдообратной матрицы\n",
    "\n",
    "* Если матрица $A$ имеет полный ранг, тогда $A^* A$ невырождена, и мы получим \n",
    "\n",
    "$$A^{\\dagger} = \\lim_{\\alpha \\rightarrow 0}(\\alpha I + A^* A)^{-1} A^* = (A^* A)^{-1} A^*.$$\n",
    "\n",
    "* Если матрица $A$ квадратная и невырожденная, мы получим \n",
    "\n",
    "$$ A^{\\dagger} = \\lim_{\\alpha \\rightarrow 0}(\\alpha I + A^* A)^{-1} A^* = (A^* A)^{-1} A^* = A^{-1} A^{-*} A^* = A^{-1}$$\n",
    "\n",
    "Получили обычную обратную матрицу для $A$\n",
    "\n",
    "* Если $A$ имеет линейно зависимые столбцы, тогда $A^\\dagger b$ даёт решение минимальной евклидовой нормы – вспоминаем геометрию задачи! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Вычисление псевдообратной матрицы с помощью SVD\n",
    "\n",
    "Пусть $A = U \\Sigma V^*$ SVD для матрицы $A$. Тогда,\n",
    "\n",
    "$$A^{\\dagger} = V \\Sigma^{\\dagger} U^*,$$\n",
    "\n",
    "где $\\Sigma^{\\dagger}$ состоит из обращённых ненулевых сингулярных чисел матрицы $A$. Действительно,\n",
    "\n",
    "\\begin{align*}\n",
    "A^{\\dagger} = \\lim_{\\alpha \\rightarrow 0}(\\alpha I + A^* A)^{-1} A^* = \\lim_{\\alpha \\rightarrow 0}( \\alpha VV^* + V \\Sigma^2 V^*)^{-1} V \\Sigma U^* = \\lim_{\\alpha \\rightarrow 0}( V(\\alpha I + \\Sigma^2) V^*)^{-1} V \\Sigma U^* = V \\lim_{\\alpha \\rightarrow 0}(\\alpha I + \\Sigma^2)^{-1} \\Sigma U^* = V \\Sigma^{\\dagger} U^*,\n",
    "\\end{align*}\n",
    "\n",
    "* Вы можете проверить, что $\\Sigma^{\\dagger}$ состоит из обращённых ненулевых сингулярных чисел <br>\n",
    "* Если сингулярные числа малы, их можно не обращать. Это даст решение менее чувствительное к шуму в правой части\n",
    "\n",
    "**Q:** что произошло с числом обусловленности?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Для решения задачи используем QR разложение \n",
    "\n",
    "Любая матрица может быть представлена в виде \n",
    "\n",
    "$$\n",
    "A = Q R,\n",
    "$$\n",
    "где $Q$ – унитарная матрица, и $R$ – верхнетреугольная.\n",
    "\n",
    "- Ортогонализация Грама-Шмидта неустойчива\n",
    "- Модифицрованный алгоритм устойчив"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- Вычисление с помощью отражений Хаусхолдера. Это преобразование зануляет подвектор\n",
    "\n",
    "$$\n",
    "    H \\begin{bmatrix} \\times \\\\ \\times \\\\ \\times \\\\ \\times  \\end{bmatrix} = \n",
    "      \\begin{bmatrix} \\times \\\\ 0 \\\\ 0 \\\\ 0  \\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "- Матрица Хаусхолдера \n",
    "\n",
    "$$ H = I - 2 vv^*. $$\n",
    "\n",
    "Сложность умножения на вектор $\\mathcal{O}(n)$!\n",
    "\n",
    "- Сложность вычисления QR разложения – $\\mathcal{O}(n^3)$\n",
    "- Какая сложность для матрицы $m \\times n,$ где $m > n$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- Вычисление с помощью вращения Гивенса\n",
    "- Матрица Гивенса в 2D имеет вид\n",
    "\n",
    "$$\n",
    "    G = \\begin{bmatrix}\n",
    "          \\cos \\alpha & -\\sin \\alpha \\\\\n",
    "          \\sin \\alpha & \\cos \\alpha\n",
    "        \\end{bmatrix},\n",
    "$$\n",
    "- Последовательное зануление элементов\n",
    "- Подходит для разреженных матриц и параллельных архитектур"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Возвращаемся к задаче наименьших квадратов\n",
    "\n",
    "Тогда если $A$ имеет полный ранг, то\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",
    "Таким образом, необходимо решить следующую квадратную систему \n",
    "\n",
    "$$\n",
    "Rx = Q^* b.\n",
    "$$\n",
    "\n",
    "- $R$ верхнетреугольная\n",
    "- Решение требует $\\mathcal{O}(n^2)$ операций\n",
    "- Более устойчивый способ, чем использование псевдообратной матрицы напрямую"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Альтернатива для недоопределённой системы: compressed sensing\n",
    "\n",
    "- Альтернативным способов определения единственности решения является задача минимизации его первой нормы\n",
    "\n",
    "\\begin{align*}\n",
    "& \\Vert x \\Vert_1 \\rightarrow \\min_x \\\\\n",
    "\\mbox{s.t.  } & Ax = f\n",
    "\\end{align*}\n",
    "\n",
    "- Решается специальными методами оптимизации\n",
    "- Ожидается получение разреженного решения"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Всегда ли решение будет разреженным?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 7.07459892e-01 -7.06752786e-01  1.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]\n",
      " [-7.06752786e-01  7.07459892e-01  0.00000000e+00  1.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]\n",
      " [ 3.53553170e-04  3.53553170e-04  0.00000000e+00  0.00000000e+00\n",
      "   1.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]\n",
      " [ 3.53553170e-04  3.53553170e-04  0.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00  1.00000000e+00  0.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]\n",
      " [ 3.53553170e-04  3.53553170e-04  0.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00  1.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]\n",
      " [ 3.53553170e-04  3.53553170e-04  0.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00  0.00000000e+00  1.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]\n",
      " [ 3.53553170e-04  3.53553170e-04  0.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00\n",
      "   1.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]\n",
      " [ 3.53553170e-04  3.53553170e-04  0.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00  1.00000000e+00  0.00000000e+00  0.00000000e+00]\n",
      " [ 3.53553170e-04  3.53553170e-04  0.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00  1.00000000e+00  0.00000000e+00]\n",
      " [ 3.53553170e-04  3.53553170e-04  0.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00  0.00000000e+00  1.00000000e+00]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "eps = 1e-3\n",
    "n = 10\n",
    "b = np.ones(n)\n",
    "A = np.zeros((n, n+2))\n",
    "A[:, 2:] = np.eye(n)\n",
    "A[:, :2] = eps / 2.\n",
    "A[0, 0] += 1\n",
    "A[1, 1] += 1\n",
    "A[0, 1] -= 1\n",
    "A[1, 0] -= 1\n",
    "c = 1. / np.linalg.norm(A[:, 0])\n",
    "A[:, 0] = A[:, 0] / np.linalg.norm(A[:, 0])\n",
    "A[:, 1] = A[:, 1] / np.linalg.norm(A[:, 1])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 10.0\n",
      "10\n"
     ]
    }
   ],
   "source": [
    "x1 = np.zeros(n+2)\n",
    "x1[2:] = 1\n",
    "print(np.linalg.norm(A @ x1 - b), np.linalg.norm(x1, 1))\n",
    "print(np.sum(x1 != 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 2828.4288925125907\n",
      "2\n"
     ]
    }
   ],
   "source": [
    "x2 = np.zeros(n+2)\n",
    "x2[:2] = 1\n",
    "x2 = x2 / (eps * c)\n",
    "print(np.linalg.norm(A @ x2 - b), np.linalg.norm(x2, 1))\n",
    "print(np.sum(x2 != 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Почему так получилось?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Задача 2: вычисление собственных векторов и собственных значений\n",
    "\n",
    "$$ Ax = \\lambda x, \\qquad A = S \\Lambda S^{-1} $$\n",
    "\n",
    "- Степенной метод\n",
    "- QR алгоритм"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Теоретический подход\n",
    "\n",
    "- Характеристическое уравнение\n",
    "\n",
    "$$ \\det(A - \\lambda I) = 0 $$\n",
    "\n",
    "- Решение через вычисление коэффицентов полинома и нахождения его нулей\n",
    "- Очень неустойчивая операция! \n",
    "- Почему?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Практический подход: частичная задача на собственные векторы\n",
    "\n",
    "- Степенной метод и блочный степенной метод\n",
    "- Одна итерация\n",
    "\n",
    "$$ x_{k+1} = A x_k, \\quad x_{k+1} := \\frac{x_{k+1}}{\\Vert x_{k+1} \\Vert_2}.$$\n",
    "\n",
    "- Сходится к максимальному по модулю собственному значению\n",
    "- Для выполнения итераций достаточно иметь функцию умножения матрицы на вектор!\n",
    "\n",
    "**Q:** какая сложность?\n",
    "\n",
    "- Скорость сходимости – линейная (геометрическая прогрессия), константа $q = \\left|\\frac{\\lambda_{2}}{\\lambda_{1}}\\right| < 1$, где $\\lambda_1>\\lambda_2\\geq\\dots\\geq \\lambda_n$.\n",
    "- Не сходится, если $\\lambda_1 = \\bar{\\lambda}_2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Крыловское пространство\n",
    "\n",
    "- Собственный вектор, полученный степенным методом, лежит в **Крыловском подпространстве** $\\{x_0, Ax_0,\\dots,A^{k}x_0\\}$ и имеет вид $\\mu A^k x_0$, где $\\mu$ нормировочная постоянная. \n",
    "- Далее Крыловские пространства будут играть **ключевую роль** для задач большой размерности!\n",
    "- Подробности через лекцию..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Практический подход: полная задача на собственные векторы\n",
    "\n",
    "- Ищем весь спектр!\n",
    "- Разложение Шура\n",
    "\n",
    "$$A = UTU^*,$$ \n",
    "\n",
    "где $U$ унитарная, а $T$ верхнетреугольная.\n",
    "\n",
    "- Спектр $T$ и $A$ совпадает в силу унитарности $U$\n",
    "- Как вычислить разложение Шура?  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "###  QR алгоритм\n",
    "1. Инициализируем $A_0 = A$.\n",
    "2. Вычислим QR разложение матрицы $A_k$: $A_k = Q_k R_k$.\n",
    "3. Обновим аппроксимацию $A_{k+1} = R_k Q_k$.\n",
    "\n",
    "\n",
    "- Как проверить сходимиость?\n",
    "- Какая сложность?\n",
    "- Как найти матрицы $U$ и $T$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Сходимость и ускорение\n",
    "\n",
    "- Ускорим QR алгоритм до $\\mathcal{O}(n^3)$\n",
    "- Используем свойства матриц в **верхне-гессенберговой форме**\n",
    "- Матрица $A$ представлена в верхне-гессенберговой форме, если\n",
    "\n",
    "$$a_{ij} = 0, \\quad \\mbox{if } i \\geq j+2.$$\n",
    "\n",
    "или \n",
    "\n",
    "$$A = \\begin{bmatrix} * & * & * & * & * \\\\ * & * & * & * & * \\\\ 0 & * & * & * & *\\\\ 0 & 0 & * & * & *\\\\ 0 & 0 & 0 & * & * \\\\ \\end{bmatrix}.$$\n",
    "\n",
    "- Подробности на следующей лекции..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Задача 3: вычисление матричных функций\n",
    "\n",
    "$$ B = f(A) = \\sum_{k=0}^{\\infty} c_k A^k$$\n",
    "\n",
    "- Это не поэлементная функция\n",
    "- Нужно уметь вычислять ряд\n",
    "- Устойчивость"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Вычисление степени матрицы\n",
    "\n",
    "- Дана матрица $A$\n",
    "- Нужно вычислить $A^k$\n",
    "- Это необходимо при вычисление матричных функций, которые представляются в виде матричного ряда"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Используем спектральное разложение\n",
    "\n",
    "- $A = U\\Lambda U^{-1}$\n",
    "- Тогда $A^k = U \\Lambda^k U^{-1}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Матричная экспонента: знакомство или напоминание\n",
    "\n",
    "$$\\exp(A) = \\sum_{k=0}^{\\infty} \\frac{A^k}{k!} $$\n",
    "\n",
    "- Про матричные функции будет отдельное занятие!\n",
    "- Важно! Чаще нужно не само значение матричной экспоненты, а результат её умножения на некоторый вектор $v$: \n",
    "\n",
    "$$ \\exp(A)v $$\n",
    "\n",
    "- Это позволяет разработать быстрые методы для вычисления произведения, не считая саму матричную экспоненту!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10.14557721316647\n",
      "[[  88.1374  -20.1911   -4.2125    3.0542    6.2766   -1.0849  -28.8463]\n",
      " [ -20.1911   18.0428   -5.4419   41.6549  -13.4495  -18.449    24.4322]\n",
      " [  -4.2125   -5.4419   21.6938  -53.9761   40.6681    2.6968  -33.6354]\n",
      " [   3.0542   41.6549  -53.9761  204.2881 -105.3993  -51.6985  108.6377]\n",
      " [   6.2766  -13.4495   40.6681 -105.3993   80.2581    6.1887  -70.7844]\n",
      " [  -1.0849  -18.449     2.6968  -51.6985    6.1887   30.0549  -17.7194]\n",
      " [ -28.8463   24.4322  -33.6354  108.6377  -70.7844  -17.7194   74.2981]]\n"
     ]
    }
   ],
   "source": [
    "import scipy.linalg as splin\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rc(\"text\", usetex=True)\n",
    "%matplotlib inline \n",
    "\n",
    "n = 7\n",
    "A = np.random.randn(n, n)\n",
    "A = A + A.T\n",
    "print(np.linalg.norm(A))\n",
    "expA = splin.expm(A)\n",
    "print(expA.round(4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def qr_algorithm(A, num_iter, eps):\n",
    "    T = A.copy()\n",
    "    U = np.eye(A.shape[0])\n",
    "    conv = [(T, U)]\n",
    "    for i in range(num_iter):\n",
    "        Q, R = np.linalg.qr(T)\n",
    "        T = R @ Q\n",
    "        U = U @ Q\n",
    "        conv.append((T, U))\n",
    "        if np.sum(np.abs(np.tril(T, k=-1))) < eps:\n",
    "            break\n",
    "    return T, U, conv[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0675878641639698e-13\n"
     ]
    }
   ],
   "source": [
    "T, U, conv = qr_algorithm(A, 2000, 1e-6)\n",
    "print(np.linalg.norm(A - U @ T @ U.T))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  88.1374  -20.1911   -4.2125    3.0542    6.2766   -1.0849  -28.8463]\n",
      " [ -20.1911   18.0428   -5.4419   41.6549  -13.4495  -18.449    24.4322]\n",
      " [  -4.2125   -5.4419   21.6938  -53.9761   40.6681    2.6968  -33.6354]\n",
      " [   3.0542   41.6549  -53.9761  204.2881 -105.3993  -51.6985  108.6377]\n",
      " [   6.2766  -13.4495   40.6681 -105.3993   80.2581    6.1887  -70.7844]\n",
      " [  -1.0849  -18.449     2.6968  -51.6985    6.1887   30.0549  -17.7194]\n",
      " [ -28.8463   24.4322  -33.6354  108.6377  -70.7844  -17.7194   74.2981]]\n",
      "1.3182919857724327e-05\n"
     ]
    }
   ],
   "source": [
    "cur_term = np.eye(n)\n",
    "expA_approx = np.eye(n)\n",
    "K = 50\n",
    "denom = 1.\n",
    "conv = [expA_approx]\n",
    "for i in range(1, K):\n",
    "    denom = denom * i\n",
    "    cur_term = U @ T**i @ U.T / denom\n",
    "    expA_approx = expA_approx + cur_term\n",
    "    conv.append(expA_approx)\n",
    "print(expA_approx.round(4))\n",
    "print(np.linalg.norm(expA_approx - expA))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoIAAAG/CAYAAADM/Lb7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xl8VvWd9//3JwuELQlhU3ZCFHcgBNdaQQGdKk61AZdOa9sRYtuZTqdV0d7tPXP/Or0tdpkuc1cJTut0UQS0tnRTwKLWpRLAfUMCyKbIkgAqaz6/P64TG7NfW865cr2ej0ceIec6OefT+ZbOm+9q7i4AAABkn5ywCwAAAEA4CIIAAABZiiAIAACQpQiCAAAAWYogCAAAkKUIggAAAFmKIAgAAJClCIIAAABZiiAIAACQpfLCLiBTDBw40EePHp3Wd7z77rvq06dPWt+BxNE+0UXbRBvtE220T3Ql0zZr1qzZ5e6DOrqPINhJo0ePVk1NTVrfsWrVKk2ZMiWt70DiaJ/oom2ijfaJNtonupJpGzPb3Jn7GBoGAADIUgRBAACALEUQ7ICZzTSz6vr6+rBLAQAASCmCYAfcfZm7zy0qKgq7FAAAgJQiCAIAAGQpgiAAAECWIggCAABkqazbR9DMSiVVSlorqVxStbvXhVsVAABA18u6IChpibtPkiQzq5G0UNKscEsCAADoehkZBM2sUlKdu69o5bNySRWSaiWVSqptvC/4bE/jve5eZ2bTuqZqAACAaMm4IBgEt1Z78YJh3/nuPr3JtSVmVuvujcGw+TDwHjMrd/e16awbAAAgajJmsYiZlZrZAsXC3J42bquStKDZtQWS5gd/LklTeQAAABknY4Kgu9e6e5W7V7dzW+MikKZqgutSLEAWN/uccAgAALJSxgTBjphZsVrpLWxcERwMG9eqleDHsDAAAMhGGTdHsB0l0t+CXytK3X2FmX1wIQiHLRachGHfwSPasr9Bm3e/q4L8XBXk56pXfq7yc01NawYAAEiV7hQEmw/5tmWOmd2sv+0jOCd9JXXe0xt26xtPvC89sepD13NM6pWfq149cj8UEAf07aHjCgt0XFGBji8q0HFFvYLvBerXM4/wCAAAOmTuHnYNcTOzDZKqmm4fE2wNs8bdWyQgM3NJ01vbbqaD98yVNFeShgwZMmnRokXJFd6OvQcb9OJb7yknv6cOH5MON0iHjrkOH5OOHHMdalDsenCt/rBr78EG7Tvc8lkFuVJxgamkwDS4V47GFudobHGujutjyiEgJuzAgQPq27dv2GWgFbRNtNE+0Ub7RFcybTN16tQ17l7R0X3dqUdQUmyuYKpOCgkWplRLUkVFhU+ZMiUVj21T/1WrFO87Dh09pp37DumtfQe1o/6g3qp/X2/VH9Jb+97XjvqDWrPzgFZtjaXFwoI8TRjZX+UjizVxZH9NGFGsol75afhP0j2tSqB90DVom2ijfaKN9omurmib7hQEa4PvJWqyV2CwiKTp53Exs5mSZpaVlSVXXZr0zMvViJLeGlHSu9XPGxpctbsOaO3mOq3bslfr3qzTD1euV2NHcNngviofWaxJo/pr2slDNKBvzy6sHgAAhKnbBMHglJBatb49TF2woXQiz10maVlFRUUk5hLGKyfHVDa4n8oG99PsySMkSfsPHtHzW+u1dvNerdtSp+Uvv63FNVuVl/OipowbpCvLh+vCkwarID835OoBAEA6dZsgGFih2PFyTbeDKVcSK4Oj3iOYiH4F+TqvbKDOKxsoSXJ3vbJjv37z7DY9+Ow2rXhlp/oV5OmyM47XleXDVTGqP4tPAADohjI1CJao9VXC8yQtUTCvL1AVfCUk03sEO8PMdMrQQp0ytFA3X3KSntywS79eu00Prtuue5/ZohElvXTFxOG6cuIwjR7YJ+xyAQBAimRMEAzm+t2q2KbRxZLmm9l0Scvdfan0wfDwPDObL2l1cO/8RIeFs1Fujun8Ewbp/BMG6ZsfP6qHXnpLD6zdph8/sl4/WrleE0cW67pzRmvm+KHKzaGXEACATJYxQTBYCTyvE/etVctj5hLWHYeGO6tPzzxdWT5cV5YP11v1B/WbZ7dpyZqt+vJ9z+qOVRt008XjdNHJgxk2BgAgQ3WbI+bSxd2XufvcoqKisEsJ1XFFBaq6YKwe/vJH9aNrJurQ0WO6/uc1+sQdT+rp2t1hlwcAABJAEERccnJMl48fquVfuUD/94rTta3ufV1d/bQ+/dNn9OK2+rDLAwAAcSAIdsDMZppZdX09Iaep/NwcXXvWSD1601Td+ncn6fmtdbrsx3/RF+9Zq9p3DoRdHgAA6ASCYAcYGm5fQX6uqi4Yq8dunqp/vrBMf351p6b/52O65f7ntaP+/bDLAwAA7SAIIiUKC/L11Rnj9OhNU/Wps0fpgbXbdOF3H9U9f31TmXieNQAA2YAg2AGGhuMzqF9P/fvlp2rlVy9Q+ahife3XL2jOz2u068ChsEsDAADNEAQ7wNBwYkaU9NYvPneWvnHZKXps/S5d/J+PacXLb4ddFgAAaIIgiLTJyTH940fGaNk/fUSDCwt0/c9rdOsDL+jdQ0fDLg0AAIggiC4w7rh+evCL56rqglItWv2mLv3R41r75t6wywIAIOsRBNEleubl6ta/O1mL5pytI8dcs+58St9f/rqOHGsIuzQAALIWQbADLBZJrbNKB+iPXz5ffz9+qH60cr0q73iSfQcBAAgJQbADLBZJvcKCfH3/qgn6f9eWa/Oe93Tpj/6i3z+/I+yyAADIOgRBhObSM47XQ1/+qE4dWqgv3rNWCx7dwJ6DAAB0IYIgQjWksEC/vP4sXXr68brtj6/qG795UUeZNwgAQJfIC7sAoCA/Vz++ZqKG9++lBY/VanvdQf34monq05P/egIAkE70CCIScnJMt37sZH3z46dp1Ws7dVX1U9q572DYZQEA0K0RBDvAquGu9amzR+mu6ypU+867uuInT+r1t/eHXRIAAN0WQbADrBrueheeNET3zT1Hh4816BN3PKknN+wKuyQAALolgiAi6fThRfr1F87VcYUFuu6nz+iBtVvDLgkAgG6HIIjIGt6/t5Z+/lxVjCrRVxY/px+tXM/2MgAApBBBEJFW1Ctf//O5M3XlxGH6/vLX9fUHXyQMAgCQIuzPgcjrkZej780er0GFPbXg0Vrl5+bo32aeIjMLuzQAADIaQRAZwcx0yyUn6chR10+f2KiC/FzNu2QcYRAAgCQQBJExzEzfuOxkHTp6THc+ukG98nP1L9NOCLssAAAyFkGwA2Y2U9LMsrKysEuBYmHwm39/mg4dbdB/rnhdBfk5qrpgbNhlAQCQkVgs0gH2EYyenBzT/E+cocvOiJ1PfPcTG8MuCQCAjESPIDJSbo7pP6+aoMNHG/Tvy15Wz/xcXXPmyLDLAgAgo9AjiIyVn5ujH187URecOEhf+/UL+vU6Np0GACAeBEFktJ55uVrwqUk6p3SAvrr4Of3hhR1hlwQAQMYgCCLjFeTn6q7rKlQ+sr++dO86rXj57bBLAgAgI2RtEDSz5WHXgNTp3SNPP/vsZJ06tFBf+NVaPb7+nbBLAgAg8rIuCJrZNDObK2la2LUgtfoVxI6jGzu4r+b8vEYvbqsPuyQAACItMkHQzCrNrNVwZmblZja3McS1dV9nuPsKd69OvFJEWXHvHvr5587UgD49df3/1GjnvoNhlwQAQGRFIggGwW5hG5+VSprv7tVNQlxVcB1oYVC/nlr46QrtO3hEc35eo4NHjoVdEgAAkRRqEDSzUjNbIKlU0p42bquStKDZtQWS5qezNmS2U4YW6odXT9Tz2+p109Ln5e5hlwQAQOSEuqG0u9cqFvRkZvPauK1SLYNgjaQPFnsEc/7aO2dsubuvSKJUZKDppwzRzRefpPl/elUnDO6rL13EucQAADQV6ZNFzKxYrfQWunudmcnMSt29ljl/aMsNF5Rq/c79+v7y1zV2UF9desbxYZcEAEBkRGKOYDtKpFjwa+Nz5gmiXWam2648XRWj+uurS57VC1tZSQwAQCOLytwpM9sgqarpEK6ZlUta4+7Wyv0uaXq8Q77BM6cpNsfwdrUzbBwMOc+VpCFDhkxatGhRPK+K24EDB9S3b9+0viNb7Tvk+j9Pva9jLv3bOQXqXxD/v4Fon+iibaKN9ok22ie6kmmbqVOnrnH3io7ui/TQcDq4+1pJaxULgR3dWy2pWpIqKip8ypQpaa1t1apVSvc7stkJZ+xT5R1P6mdv9NB9c89Rrx65cf0+7RNdtE200T7RRvtEV1e0TdSHhiV9MFcwrHfPNLPq+nqGFDPdycfHVhK/sK1eNy59jpXEAICsF/UgWBt8L2l6sUkwrFWaufsyd59bVFSU7lehC0w7ZYhuueQk/f75HfrhyvVhlwMAQKgiHQSDRSK1kpr3CJZIqgu2n0kregS7n7kfLVXlpOH6wYr1+t3z28MuBwCA0EQ6CAZWSGo+2bE8uJ529Ah2P2amb11xmiaP7q+vLn6OM4kBAFkrSkGwRC17/iRpnqRZza5VBdeBhPTMy9Wd/zBJJX166Iv3rNX+g0fCLgkAgC4X9hFzxWY238yWKBYC55vZAjOrbLwnGB6eF9xXaWY3K3b2cNqHhYMaGRrupgb07akfXTNRW/e+r6/9+kUWjwAAsk7YR8zVqRM9e022fOly7r5M0rKKioo5Ybwf6TV5dIm+Mv1Efeeh13Te2AG6+syRYZcEAECXidLQMBCKz18wVuefMFD/9tuX9Npb+8MuBwCALpNwEDSzCWY2OnWlRBNDw91fTo7p+7MnqF9Bvv7pnrV67/DRsEsCAKBLJNMjuERSuZldmKpioohVw9lhUL+e+sFVE/TGOwf07799KexyAADoEskEwfvd/QF3fyRl1QAh+sgJA/XFKWVaXLNVD67bFnY5AACkXTJBcEPjH8xsjplNSEE9QKi+PO0ETR7dX//r1y+o9p0DYZcDAEBaJRMEP9hrw90XSlpqZtcnX1K0MEcwu+Tl5uhH10xUfl6O/umedTp45FjYJQEAkDbJBMFJzX6e7+53JVNMFDFHMPscX9RL360cr5d37NNtf3gl7HIAAEibZIJglZkdM7OHzOxGSR9KSgwVI5NNO2WI/vEjY/Q/T23Wn17cEXY5AACkRVJBULFj4aolnSnpa2a228zuC4aIq1JRIBCWeZecpDOGF+nmpc9ry573wi4HAICUSzgIuvtCd6939/vdfba7l0iqkLRC0mxJc1NVZJiYI5i9euTl6L+uKZe79M/3rtPRBo6gAwB0Lyk9WcTdNwYBcYakW1L57LAwRzC7jRzQW7d94nQ9u6VOD6w/EnY5AACkVDqPmKtO47OBLnPZGUN19eQR+uPGI1qzeW/Y5QAAkDJpC4Luzlgquo3/denJKikw3bTkObaUAQB0G+nsEQS6jX4F+frcaT1Vu+tdfe/h18IuBwCAlMhr7aKZ3anYwo9kZ8ebpA3uflWSzwFCd+rAXF171kjd9ZeNuuS04zRpVEnYJQEAkJRWg6C739DVhUSVmc2UNLOsrCzsUhABX/vYyXr0tXd005Ln9Yd/OV8F+blhlwQAQMIYGu4Aq4bRVN+eebq98gzV7npX332IIWIAQGYjCAJxOq9soD551kj99xMbVbNpT9jlAACQsLbmCH5b0pgUvaPW3W9N0bOASLj1Yydr1Wvv6Kalz+sPXzpfvXowRAwAyDytBkFJtyl2fFwq0GWCbqdvzzx9p/IMXXvXX/Xdh1/TNy47JeySAACIW1uLReolsQ8g0I5zywbqU2eP0k+fiK0injyaVcQAgMzCHEEgCbf83Uka3r+XblrynN4/zEbTAIDMQhAEktCnZ55u/8R4bdr9nm5/6NWwywEAIC4JBUEzKzSzCWY22swKU11UlJjZTDOrrq9npBytO2fsAH36nFG6+8lNemYjU2IBAJmjU0HQzOaY2cNmtsfMdktaKel2SdWS1prZbjNbbWY3mtno9JXb9dhHEJ0x75JgiHjpc3rv8NGwywEAoFPaDIJBr9+3zewhxY6am+XuJe4+wN0nu/uM4Kus8ZqkdZJuMbP7zOzCrvoPAYStcYh48+73dPuf2GgaAJAZ2tpHcKKkuZJud/eNnX2Yu69UrLewsRex3N2/m5JKgYg7Z+wAXRcMEV8+YajKR/YPuyQAANrVokfQzMZIKnX3z8cTAptz94WSFprZlckUCGSSmy45SccVFujrv35RR481hF0OAADtahEE3X2ju9+fioe7e727P5CKZwGZoG/PPP3bzFP08o59+vlTm8MuBwCAdmXd9jFmVm5mc83sZjNbYmalYdeE7uWS047TlHGD9P3lr+vtfQfDLgcAgDalNAia2R2pfF6qmVmxpAp3r3b32yUtkLQ85LLQzZiZ/s/lp+rIsQb9f797OexyAABoU1JBsMl+ghea2fWSZifxrEozm9bGZ429eNMavyf4mlJJ85r8XCOpNAiIQMqMGtBHX5xapt8/v0OPvf5O2OUAANCqVlcNS5KZTZBUolh4GiupOPhz47Xm4ckU22YmbkGwWyhpViuflUqa7+7Tm1xbYma17l4bz3vcfa2ZTW9yqUJSnbvXJVI30J6qC0r14Lpt+t+/eVF/+vJHVZCfG3ZJAAB8SHs9grdLWqHY8Ok8SVWKBcK9im0Rs07SDYr1As6QVBZc6zQzKzWzBYoFy7aOZKgKamhqgaT58byrUbPwWCVpTiLPATrSMy9X3/z4adq0+z3d+eiGsMsBAKCFNnsE3X2GmVVKmiZpnrt/6Iw1M7sz2CKm6bUV8bw8CGVVwe/Oa+O2SrUMgjVqMrfPzOYqFlLbstzdP1Rb8Dv3ufvSeGoG4nFe2UBdPn6ofrJqgz4+YZhGD+wTdkkAAHygzSAoSUFIWmpm15tZjbs/2/TjVu6/JZXFBXP3WvQWunudmcnMSt291t2r43zuNEm1zcMhkA5fv/Rk/fnVnfrGb17Uzz93psws7JIAAJDUycUi7n6XpPoQjo0rCd7f1hy+uLd+MbNySXsaQ2DQ6wmkzeDCAn11xol6fP0u/f6FHWGXAwDAB9rtEWwqOGVkY7CIpL6j+1Mkpat5g4UnKyUVN+mVqZXU6vBwMHw8V5KGDBmiVatWpbKcFg4cOJD2dyBxybTPSHeNKszR1+9/Vrk7X1OvPHoFU4m/O9FG+0Qb7RNdXdE2nQ6CjRqHh81srZkVuvu+1JeVHsGcxE4fABsMOVdLUkVFhU+ZMiVNlcWsWrVK6X4HEpds+wwoq9PHf/KEVr8/RP975impKwz83Yk42ifaaJ/o6oq2SXgfwWChyIBgD8FCSUrXucJh7vNnZjPNrLq+vqs6QdFdjR9RrE+eNVJ3P7lRL23nv08AgPAltaF0cC7xI5Kmm9mdiu0FmEqNW72UNL3YJBjGtY9gItx9mbvPLSoqSverkAVumnGSSvr00NcffFENDQltuwkAQMqk5Ig5d79fsSHUlPbcBYtEalt5boliG0GnPQgCqVTUO19f+9jJWvdmnRat3hJ2OQCALJeys4bdfa1iCzFSbYViJ4A0VR5cTzuGhpFqV0wcprNLSzT/T69q14FDYZcDAMhiKQuCgaokfrdErfcozlPLo+eq9OEzg9OGoWGkmpnpPz5+mt49dFTfe/i1sMsBAGSxlAbBYIuZTjOzYjObb2ZLFAuB881sQdO9/YLh4XnBfZVmdrNiZw93ybAwPYJIh7LB/fTpc0brvtVb9MqOjFl4DwDoZloEQTMbk6rVv2ZWaGbXt/W5u9e5+zx3n+Xu5u5j3b2q+bFv7r42uG+pu9/elSeC0COIdPmXi05QYa98/cfvX5Y7C0cAAF2vRRBssnH0HWY2OtEHm9kcSbcGp5IAaKaod77+ddqJeuKN3Vr5ys6wywEAZKFWh4bdfZ2kWyTdYGYPBWcNF3b0MDObYGZ3mtlDkja4+60prhfoVq49a6TGDuqjb/3hFR0+2hB2OQCALNPmHEF3r3f3W9z9Ykl7Jd1lZm+Y2XozWx0ExIeDa7vNbLViizgWuPvFwf6CGY85gkin/Nwcff3SU7Rx17v6xdObwy4HAJBlOnXEXLBP4P2NP5tZkaRSSXWS9rh7t01J7r5M0rKKioo5YdeC7mnKuEE6/4SB+uGK13XlxGHq36dH2CUBALJEoquGXbEzey+SVN6ZYWMArTMzfeOyU3Tg0FH9YMXrYZcDAMgicQfBYBFInaTlis0jXKnY4pIrUlxbJDA0jK5w4pB+uvaskfrlX9/UGzv3h10OACBLxBUEzewTim3uPMndc929zN1zJE2X9HkzuzAdRYaJ7WPQVf512onq3SNX3/r9K2GXAgDIEvH2CE529xnBquIPBPv8zVAsEAJIwIC+PfXPF5bpz6+9o0dffyfscgAAWSDeILi6g89rEi0EgHTduaM1akBv/cfvXtbRY2wnAwBIr1SfNfwhyWxIDWSjnnm5uvXvTtb6nQd07+otYZcDAOjm4g2Cy9s6fi64vqbZ5fkJVRUhLBZBV7v41CE6a0yJ/nP566p//0jY5QAAurF4g+AmSUvN7FjzL0lLJdU2udYgqTLVBXc1FougqzVuJ7P3vcP6r0fWh10OAKAbizcI1koaG6wYbv6VE3x98LOabEINoPNOG1akWZOG6+4nN2nTrnfDLgcA0E3FGwTnuPvGOO6/Lc7nAwjcOGOc8nNzdNsf2U4GAJAe8QbBve192Hz+YPNtZgB03uDCAn1hylg99NLbenLDrrDLAQB0Q/EGwQUdfF6VaCEAWrr+/FINK+6lb/3+FTU0eNjlAAC6mXiD4HQzm9raB2Z2p6RpyZcULawaRpgK8nP11Rkn6qXt+/SHF3eEXQ4AoJuJNwhWSyprOgRsZhPMrHFpY7dbHMKqYYTt7ycM07gh/fS9h1/XETaZBgCkUFxB0N1vcPeFktaZ2fVmdqOklZJuCD6bnZYqgSyWm2O68eJx2rjrXS2p2Rp2OQCAbiTRk0V2S5ot6QZJs9x9ZepKAtDctJMHq3xksX648nUdPHIs7HIAAN1EXEHQzK40swsV21j6YXcvi11u/bQRAKlhZpp3yUl6e98h/c+Tm8IuBwDQTcTbI7hUsZXDF7n7dyUp6A1cZ2Z3mtnrqS4QQMxZpQM0Zdwg/WTVBo6eAwCkRCIni1Q03x8w2GT6FkkDU1UYgJZuunic6t8/ourHNoRdCgCgG4i7R9DdW91Hxd3rJC1OviQAbTl1aJFmjh+qn/5lk3buPxh2OQCADBfvquFbzGy0md1mZtc3Xjezi8ys0N1vSH2J4WIfQUTNV6efqCPHGvRfj7wRdikAgAwX72KRiZJWSDJJxY3Xg3mCk81sdCqLiwL2EUTUjB7YR1dNHqF7/vqm3tz9XtjlAAAyWLxDw7Pdvczdb5G0sekHQRgsT1llANr0pYtOUF6u6fvLXwu7FABABos3CNY0+XNrB5+WJFELgE4aUligz543Rr95brte2bEv7HIAABkq3iDYNPxZK5+XJVELgDjc8NGx6tczT999iF5BAEBi4g2C1mSRiDe5OMHMHpL0cMoqA9Cuot75umHKWK18dadWb9oTdjkAgAwU76rh+yX1N7M9khaa2Xoz263YecNL3P2RdBSZSmZWbmbTzKzSzBaYWWnYNQGJ+uy5YzS4X0/N/+Orcm9ttgYAAG2L+6xhd/+OpDGKnTV8e/C91N3vSnFt6bJSUo27L5W0RtKSkOsBEtarR66+dNEJqtm8V39+bWfY5QAAMkzcQVCS3L3e3Ve6+8Lge9Kb7AU9dNPa+KzczOYGPXlz27qvk8YEm19LEuNpyHhXTR6hUQN66/Y/vaaGBnoFAQCdl1AQTLUg2C1s47NSSfPdvdrdV7h7taSqRId0m4RASaqSNC+R5wBRkZ+bo69MP1GvvrVfy57fHnY5AIAMEmoQNLNSM1sgqVRt985VSVrQ7NoCSfOTfO/Nkpa7+4pEnwNExcwzhurk4wv1vYdf1+GjDWGXAwDIEHlhvtzdaxULejKztnrmKtUyCNZIWt74g5nNlTS2nVd9KPAF7709GGZe7u7TE6kfiIqcHNPNF4/TZ+9erQfWbtXVZ44MuyQAQAYINQh2xMyK1UpvobvXmZnMrNTda4Ph4s48r1RSpbvfHlxaLGlB43NSWjzQxaaMG6TxI4r1X39+Q1eWD1ePvEjM/AAARFirQdDM7pRUodZPD4mHSdrg7lcl+PslUot5fU2VSoonwJVKGtDs5zpCILoDM9OXp52gz/6MXkEAQOdYVPYeM7MNkqqaDuGaWbmkNe7e4hQTM3NJ0+Od42dmlfrbUXjTJd3m7mvbuHeupLmSNGTIkEmLFi2K51VxO3DggPr27ZvWdyBxmdA+7q5vPn1Q+w67vn1+L+XltHYAUPeTCW2TzWifaKN9oiuZtpk6deoad6/o6L5IDw2nQ7B/YKN2h5SDIedqSaqoqPApU6aksTJp1apVSvc7kLiMaZ/jd+qzd6/W7n5jddXk7OgVzJi2yVK0T7TRPtHVFW2TEZOIgrmCYb17pplV19cnvVUi0CWmjBuk8cOL9ONH3tCRY6wgBgC0ra05gt9W7PSQVKh191sT/d3ge4mkD+YJNgmGaZ/b5+7LJC2rqKiYk+53AakQmyt44gcriLOlVxAAEL+2hoZv09/m0SUr4dM7gtXBtZKa9wiWqIsWeZjZTEkzy8rK0v0qIGWa9gpeWT5c+bkZ0fkPAOhirQbB4Mi4qIyFrlBsBXPTBR3lwfW0o0cQmYheQQBAZ0Spm6BELXv+pNgRcLOaXeNoOKADzBUEAHQk7CPmis1svpktUSwEzjezBcEWL5I+2ENwXnBfZXA03Pyu2vuPxSLIVI29glv3vq8H1m4NuxwAQAS1GBo2s4skTZRU7e770vnyxpDXifvW6sNDw12GoWFkMuYKAgDa0+L/K7j7SkkbJd1lZg+Z2fVmVtj1pQFIlpnpX6adQK8gAKBVrXYPuPv97j7b3S+WtFfS0mwNhQwNI9NNHTdYZwwv0n/9mbmCAIAP63CcKAiFMyTNVuzs4KVmdp+ZXZn26iLA3Ze5+9yioqKwSwES0ngG8ZY97+vXa7eFXQ4AIEI6PWHI3evdfWEQCufxovFUAAAgAElEQVRKGmtmDweh8ML0lQggWY29gj/+83p6BQEAH0ho5ngQCr8ThMJbJE0ysxozu4NQCEQPvYIAgNYkvYTQ3TcGobBC0u2SZjQJhROSLzFczBFEd0GvIACguZTuJRGEwluCUFgt6QYzuy+V7+hqzBFEd0GvIACgubRtKubu69z9Bne/Kl3vABAfegUBAE2xuyyQRegVBAA0lXAQNLMJZjY6daVEE3ME0d1MHTdYpw0r1J2PbtCxBg+7HABAiJLpEVwiqby7rxJmjiC6GzPT5y8oU+2ud/WnF98KuxwAQIiSCYL3u/sD7v5IyqoB0CUuOe04lQ7so5+sekPu9AoCQLZKJghuaPyDmc3pDlvFANkiN8d0wwVj9dL2fXps/a6wywEAhCSZIPhBN4K7L1Ts6Lnrky8JQFf4+MRhOr6oQD/58xthlwIACEkyQXBSs5/nu/tdyRQDoOv0yMvR9eeX6q8b92jN5j1hlwMACEEyQbDKzI6Z2UNmdqOkD62m6C5DxawaRnd2zZkj1L93vn7y5w0d3wwA6HaSCoKSShQ7QeRMSV8zs91mdl8wRFyVigLDxqphdGe9e+Tps+eN0cpXd+rVt/aFXQ4AoIslHATdfaG717v7/e4+291LJFVIWiFptqS5qSoSQPpcd85o9emRqztW0SsIANkmHWcNL3T3GZJuSeWzAaRHUe98ffLsUVr23Ha9ufu9sMsBAHShdB4xV53GZwNIoX/8yBjl5eTozsfoFQSAbJK2IOjurK4AMsSQwgJVVgzX0pqt2rnvYNjlAAC6SDp7BAFkkKqPlupoQ4P++y8bwy4FANBF8lq7aGZ3KrbwI9mzp0zSBne/KsnnAEizUQP66LIzhuqXT2/WF6aUqah3ftglAQDSrNUg6O43dHUhUWVmMyXNLCsrC7sUIO0+P2Wsfvvcdv3PU5v0pYtOCLscAECaMTTcAfYRRDY5+fhCXXjSYP3siY167/DRsMsBAKQZQRDAh3xx6ljtfe+IFj2zJexSAABp1tYcwW9LGpOid9S6+60pehaANJs0qkRnjinRwsdr9Q9nj1KPPP69CADdVatBUNJtih0flwqcZg9kmC9MGavP/Gy1Hly3TbMnjwi7HABAmrS1WKReEvsAAlnqghMH6dShhbrz0Q36xKThys2xsEsCAKRBVo/5mNmCsGsAosjM9IUpZard9a4eeumtsMsBAKRJ1gZBM5smaW7YdQBRdclpx6l0YB/dsWqD3JPdUhQAEEUpCYJmttjMrmzy80VmVhjnMyqDcNbaZ+VmNtfMpjV+T7LeYsXmLtYl8xygO8vNMc35aKle2Favp2uZ6gsA3VGqegQXuPsDjT+4+0rFTibplCDYLWzjs1JJ89292t1XuHu1pKrgeqKmufvaJH4fyApXTBymAX16aOHjtWGXAgBIg1QFwVozGy1JZvawma2WNL2jXzKz0mCeXqnaXl1cJan5XL4FkuYnUmgQOlck8rtAtinIz9WnzxmtR17dqTd27g+7HABAirW1fUyHzGy9YoFqibs/YmYTguHhBZJWBCuP2+XutYoFPZnZvDZuq1TLIFgjaXmTWuZKGtvOq5a7+4rGIWF3Z0gY6KR/OHukfrLqDd31+EZ9+xNnhF0OACCFEg6Cku6XtFrS7KBXb6Nim1DP6kwI7IwguLXoLXT3OjOTmZW6e20wXNwZ0ySVmFnjsHVxECJXBKEUQDMD+vZU5aThWrJmq746Y5wG9esZdkkAgBRJOAi6+y3BH++XJDMbo1jQ+pqZlSvWC/f5JOsrCd7VVg9eqaROBzh3X9r0ZzNbEEeIBLLWP35kjO555k394qlN+sqMcWGXAwBIEWtrWwgzu0nSBsV6y/bF/WCzMe6+MY77N0iqcvcVTa6VS1rj7i12szUzlzS96f1xvKtYsa1j5kuaJ2lpaz2CQW/hXEkaMmTIpEWLFsX7qrgcOHBAffv2Tes7kLhsb58frj2o9XuP6XtTeqtnbrQ2mM72tok62ifaaJ/oSqZtpk6dusbdO1y4216PYLWkNZLGmNlSSfc1XRnckXhCYFcLehhvD77au69asf87qKKiwqdMmZLWulatWqV0vwOJy/b26T1qj2YveEo7+5TqU2ePCrucD8n2tok62ifaaJ/o6oq2aW/V8CTFtoXJdfermobAYGHIfcHXjfHuGRivoAcvFGY208yq6+s5cQ/ZbfLo/ho/vEg//ctGNTSwwTQAdAftBcFZ7v6d5hfNbKJiPYVjFVu9WyZppZmNT0N9jcO1Jc1qKG72edq4+zJ3n1tUVJTuVwGRZhbbYHrjrne14pW3wy4HAJAC7Q0Nt9UFdqukuqbjzkEwW2xmlYnMJ2xLsDq4VlLzHsGSoAZW+gJd6JJTj9Ow4l5a+HitZpx6XNjlAACS1F6PYFtjP9MkLf7QjbE5d7OV4CbPHVihlqeUlKuLNoVmaBj4m7zcHP3jR8Zo9aa9Wvfm3rDLAQAkqb0gaI2nhTS5cJGkIjXZzLlREAaTWUpYopY9f1JsVe+sZteqgutpx9Aw8GGzJ49Qv4I83fV4ZNeDAQA6qc0gGOwTWG1mTZcHNoavtnrj4ppBbmbFZjbfzJYoFgLnm9kCM6tsUkedpHnBfZVmdrNiZw93ybAwPYLAh/XtmadPnjVKf3xxh7bseS/scgAASehoQ+kbJN0f7NlXrNgCkbntzAOMq0ewMeR14r61ktbG8+xUcfdlkpZVVFTMCeP9QBR95tzRuuvxWv33Xzbq3y8/NexyAAAJam9oWMHxbRWKbap8u6RJ7n5Xa/ea2YWKrSYG0M0dV1SgyycM1eKaLap/70jY5QAAEtRuEGzk7uvcfaG7r2t63cw+YWY3mdlixYZrF6alyhAxNAy07vqPlOq9w8f0q2c2h10KACBBnQqCbXH3+yUtVWzxyEYzuyMlVUUIi0WA1p0ytFDnnzBQdz+xSYePNoRdDgAgAUkFQSl2lFzQWzjb3T+fiqIAZIbrzy/Vzv2H9NvntoddCgAgAUkHQQDZ66MnDNS4If101+O1cufYOQDINATBDjBHEGibmen688fo1bf26/H1u8IuBwAQJ4JgB5gjCLTv8glDNbhfTy18nBMfASDTEAQBJKVnXq6uO3e0Hl+/S6++lbKjxgEAXYAgCCBpnzxrpAryc3T3E5vCLgUAEAeCYAeYIwh0rLh3D10xcbh+vW6b9r57OOxyAACdlHAQNLMiM7sxlcVEEXMEgc75zLmjdehogxat3hJ2KQCATkqmR7BE0uRUFQIgs407rp/OKxugXzy1SUePscE0AGQChoYBpMxnzh2j7fUH9fDLb4ddCgCgEwiCAFLmwpMGa2RJb/3siY1hlwIA6ASCIICUyc0xffqcUVq9aa9e3MYCKwCIOoJgB1g1DMRn9uQR6t0jVz9jKxkAiDyCYAdYNQzEp7AgX5WThmvZc9u168ChsMsBALSDIAgg5a47d7QOH2vQPX99M+xSAADtIAgCSLmxg/rqghMH6ZdPb9bho2wlAwBRRRAEkBafPW+0du4/pD++uCPsUgAAbSAIAkiLj54wSKUD+7BoBAAiLNkgaCmpAkC3k5Nj+sx5o/Xsljqte3Nv2OUAAFqRcBB0942S5qWwFgDdzJXlw9WvZ57ufnJT2KUAAFqRVI9gEAa7NfYRBBLXt2eeZlWM0O+f36G39x0MuxwAQDPMEewA+wgCybnu3FE65q5fPb057FIAAM0QBAGk1agBfXTRSYP1q7++qUNHj4VdDgCgCYIggLT77HljtPvdw1r2HFvJAECUEAQBpN25YwfoxCF99bMnNsrdwy4HABAgCAJIOzPTZ84do5e271PNZraSAYCoSCgImlmhmU0ws9FmVpjqotLJzOabmZvZXjNbbmalYdcEZIMrJg5TUa983c0G0wAQGXmducnM5kiaJalCkkuqldT4z/pSM+sfXLtP0lJ335T6UlNmg7uzETbQxXr1yNXVZ47QXY9v1Pa69zW0uFfYJQFA1muzRzDo9fu2mT2kWPib5e4l7j7A3Se7+4zgq6zxmqR1km4xs/vM7MJ4CjGzSjOb1sZn5WY218ymNX6P59kAouFTZ4+Su+sXbCUDAJHQao+gmU2UNFfS7fFsGu3uKyWtDJ4xx8zK3f27Hf1eEOwWKtbr2PyzUknz3X16k2tLzKzW3Ws7W1sTxWZWKalO0nRJt7l7XQLPARCn4f17a/opQ7R49RZ9edoJ6pmXG3ZJAJDVWvQImtkYSaXu/vlkTg5x94WSFprZlW3dY2alZrZAUqmkPW3cViVpQbNrCyTNT7C0andf6u4rFBvKXpngcwAk4FNnj9budw/rjy+8FXYpAJD1WvQIBuEvJUfHuXu9pAfa+bxWsaAnM2vr3OJKtQyCNZKWN/5gZnMljW2nlOVB8FPT3j93XxsMOxfTKwh0jXPHDlDpwD76xdOb9fGJw8IuBwCyWqcWi4TFzIrVSm+hu9eZmcys1N1r3b26k88rl7TQ3Sc1f17KigbQrpwc0yfPHqVv/u5lvbS9XqcO5fhGAAhLSvcRNLM7Uvk8SSVSu0Et3q1fatWkdzGYm7g0sdIAJKqyfLgK8nP0y6ffDLsUAMhqlswu/8EegqWKBbbGRR0DEnzWBklVjUO4wbVySWta2+7FzFzS9Kb3d/I90/S3ADlW7SwWCYac50rSkCFDJi1atCieV8XtwIED6tu3b1rfgcTRPqn10xcP6ekdR/WDKb3VOz+5HZ1om2ijfaKN9omuZNpm6tSpa9y9oqP72hwaNrMJ+lvAGyupcZi28Vpx819RbJuZSIsnOAZDztWSVFFR4VOmTElXWZKkVatWKd3vQOJon9QaeEK9LvvxX7Sz92h95rwxST2Ltok22ifaaJ/o6oq2aW+O4O2SpikW7hr/uV4bfK1ULAwuUGz+Xp1iC0wWp6NIFnMA3c9pw4o0YUSxfvH0Zl137miZsc87AHS1NucIuvsMSbMV29+vv7vnBJtHz3D32ZJq3H2hu9/v7iuDFcBxDdN2QuM+gSVNLwaLSJp+njZmNtPMquvr69P9KiDrfPqcUdrwzrt6asPusEsBgKzU7mKRYL+9GyTNCoaKP/RxK/ffksrigl7AWrUchi6RVJfghtLx1rDM3ecWFbGyEUi1j51+vPr3zuekEQAISadWDbv7XZLq4z02LkVWKHbGcVPlSn3vY6voEQTSpyA/V7Mnj9DDL7+tt+oPhl0OAGSdTm8f4+4b3f0RM5sQnD6SaiVq2fMnSfPU8ui5quB62tEjCKTXJ88cpQZ33fsMW8kAQFeLex9Bd382OH1kbbB9TMLMrNjM5pvZEsVC4HwzWxCcBdz4vjpJ84L7Ks3sZsW2qUn7sDCA9Bs5oLemnDhI9z7zpo4cawi7HADIKgmfLOLuC81sjJlVKLZwZJ+ZXenubR4p18oz6tSJnj13XytpbaK1JsPMZkqaWVZWFsbrgazwqXNG6XN312j5y2/rY6cfH3Y5AJA1kjpZpHG4WNJ0M7tTsRXG3QpDw0D6XXDiYA3v30s/f2pT2KUAQFZJyRFz7n6/YhsvtzbHDwDalZtj+uRZo/R07R6tf3t/2OUAQNZI2VnDwfDtylQ9LypYNQx0jdkVw9UjN0e/ZCsZAOgyKQuCgaoUPy90DA0DXWNA35667Izjdf/abXr30NGwywGArJDSIBisJgaAhPzDOaN04NBRPfjstrBLAYCs0CIIBiuBr0zFw82s0MyuT8WzwsLQMNB1Jo4o1qlDC/WLpzbLvcXhRQCAFGsRBINevY1mdoeZjU70wWY2R9KtwakkGYuhYaDrmJk+dfYovfrWftVs3ht2OQDQ7bU6NOzu6yTdIukGM3vIzK7vzObRwakjd5rZQ5I2uPutKa4XQDd3+YSh6leQp188xaIRAEi3NjeUdvd6xcKgzOwTku4ys3JJLqlO0p7g1rGS+kuqlVQjaUEQJAEgbr175Kly0nD98unNemf/KRrUr2fYJQFAt9Wpk0WCfQLvb/zZzIoklSoIhEFoBICU+IezR+lnT2zS4pot+uJUTvUBgHRJaNWwu9e7+7rgZJFuHQJZLAJ0vbGD+uq8sgH61dObdayBRSMAkC4JBcFgNfCFwdzBCzszfzBTsVgECMcnzxql7fUH9djr74RdCgB0W3EHwWA1cJ2kFYrNIVyh2CrjK1JcG4AsNu3kIRrYt4fueebNsEsBgG4rriAYLBqZJWmSu+e4e5m750iaLunzZnZhOooEkH165OVoVsUIPfLqTr1VfzDscgCgW4q3R3Cyu89ovirY3de6+wzFAiEApMTVk0foWINrcc2WsEsBgG4p3iC4uoPPaxItBACaGzWgjz5SNlD3rd7CohEASIOUnjXcXDInk0QFq4aBcF1z5khtq3tfj61n0QgApFq8QXB5W+cQB9fXNLs8P6GqIoRVw0C4pp8SLBr5K4tGACDVOrWhdBObJBWbWWtjNCZJTT4zxU4hAYCE9cjLUeWkEVr4eK3eqj+o44oKwi4JALqNeHsEayWNdffcVr5ygq8PflaT00gAIFGNi0aWsGgEAFIq3iA4x903xnH/bXE+HwBaGD0wtmhkEYtGACCl4g2Ce9v7sPn8webbzABAolg0AgCpF28QXNDB51WJFgIA7Zl+yhAN6NND97JoBABSJt4gON3Mprb2gZndKWla8iUBQEs98nJUWTFcK1/dqbf3cdIIAKRCvEGwWlJZ0yFgM5tgZuuDH7vd4hD2EQSi45rJI2Mnjaxm0QgApEJcQdDdb3D3hZLWmdn1ZnajpJWSbgg+m52WKkPEPoJAdIwe2EfnlQ1g0QgApEiiJ4vsljRb0g2SZrn7ytSVBABta1w08jiLRgAgaXEFQTO70swuVGxj6YfdvSx2ufXTRgAg1Wacclxs0cgzLBoBgGTF2yO4VLGVwxe5+3clKegNXGdmd5rZ66kuEACaalw0suIVFo0AQLISOVmkovn+gMEm07dIGpiqwgCgLVcHi0Y4aQQAkhN3j6C7t7p81t3rJC1OvqT0M7PKpl9h1wMgPmMG9tG5Ywfo3me2qIFFIwCQsHhXDd9iZqPN7DYzu77xupldZGaF7n5D6ktMLTO7WZLcfamkFZJuDbciAIngpBEASF68i0UmKhaeTFJx4/VgnuBkMxudaCFB71yrG1KbWbmZzTWzaY3fE3xHsaRbgxAod69z90mJ1gwgPBefyqIRAEhWXpz3zw5WCsvMPtH0A3dfGawe3hRvEUGwWyhpViuflUqa7+7Tm1xbYma17l4b56sqJNUGw8F1ksoVG+6O9zkAQtYjL0eVk4brrr9s1N8NKgi7HADISPHOEaxp8ufWJuaUxPMwMys1swWSSiXtaeO2KrU843iBpPnxvCtQqlj4W+HuKxQ7KWV5As8BEAFXnxlbNPL4tqNhlwIAGSneHsGm4c9a+bwsrofFeuKqJMnM5rVxW6VaBsEaNQlwZjZX0th2XrU8CH61kmqDhS1y97ogjJbSKwhknsZFI49u3aOGBldOTmv/swQAaEu8QdDM7Hp3v0tNQqGZTVCshy6RXrr2XlasVnoLgwCnxgDn7tWdfGRrYa8u2ToBhOeaM0fqn+/drcff2KULThwUdjkAkFHiXTV8v6T+ZrZH0kIzW29muxU7b3iJuz+S4vpKgve2FdZK43lY0OtXFwTMxqCZyFxDABEx49Qh6psvLV7NnoIAEC9zj38PLjMrUmzhRalivWw1be0vGMczN0iqCoZwG6+VS1rj7i3Ge8zMJU1ven8n31Oq2HD0BsWGkxe0FQSDIee5kjRkyJBJixYtiudVcTtw4ID69u2b1ncgcbRPdN39/AE9vsP0g6m91a8Hw8NRw9+daKN9oiuZtpk6deoad6/o6L54h4YlSUHoWxl8ZZQg9LU1H7H5vdWKLShRRUWFT5kyJY2VSatWrVK634HE0T7RtWX/I1q1/X293WuUZp4f10ABugB/d6KN9omurmibeFcNh6JxKDekd880s+r6+qQ6PAGk0Yh+ORo/oliLa7YokVEOAMhWUQ+CjUO2H9qWpkkwTPvcPndf5u5zi4qK0v0qAEm4qmKEXn/7gJ7dwvovAOisSAfBYJFIrZqcYhIokVTXFYs86BEEMsPM8cerV36u7mPRCAB0WqSDYGCFYgtTmioPrqcdPYJAZuhXkK9Lzzhey57brncPscE0AHRGq4tFzOxOxcJXspNtTNIGd7+qE/eWqGXPnxRb2LFEwaKNQFXwBQAfuGryCC1ds1W/f2GHZleMCLscAIi8VoOgu9/QFS8P5vrdqtg2NMWS5pvZdMVOAlka1FJnZvPMbL6k1cG987tq7z8zmylpZllZXIemAAhBxaj+Kh3UR/et3kIQBIBOSGj7mFQJ5gB2uJWLu6+VtDb9FbX67mWSllVUVMwJ4/0AOs/MdFXFCN32x1f1xs79KhvcL+ySACDSMmGOIAB02pXlw5WXY1pcszXsUgAg8tqaI/htSWNS9I5ad781Rc/qcgwNA5llUL+euujkwXpg7VbdOGOceuTx710AaEtbQ8O3qdnefUnYk6LnhIKhYSDzXDV5hB566W098urbuuS048MuBwAiq63FIvWS2DgPQEb66AmDNKSwp+5bvYUgCADtYMwEQLeTl5ujWZNG6NHX39GO+vfDLgcAIosg2AFOFgEy0+yKEWpwaSmLRgCgTQTBDnCyCJCZRg7orXPHDtDiNVvU0JDs3vgA0D0RBAF0W1dNHqEte97X07W7wy4FACKJIAig27r41ONUWJCnRau3hF0KAEQSQbADzBEEMldBfq6umDhMf3rpLdW/dyTscgAgcgiCHWCOIJDZZk8eocNHG/Tgs9vCLgUAIocgCKBbO3VokU4bVqhFq7fInUUjANAUQRBAt3fV5JF6Zcc+vbhtX9ilAECkEAQBdHuXjx+qnnk5uq/mzbBLAYBIIQgC6PaKeuXrY6cfr9+s2673Dx8LuxwAiAyCYAdYNQx0D1dNHqH9h47qjy/uCLsUAIgMgmAHWDUMdA9njSnRmIF9dO8zDA8DQCOCIICsYGa6evIIrd60V6+/vT/scgAgEgiCALJG5aTh6pGbo3v+Sq8gAEgEQQBZZEDfnrr4tOP0wNqtOniERSMAQBAEkFWuPXOk9h08qt89z6IRACAIAsgqZ5eWqHRQH93z181hlwIAoSMIAsgqZqZrzxyptW/W6ZUdnDQCILsRBDvAPoJA91M5abh65LFoBAAIgh1gH0Gg+ynu3UOXnn68Hly3Te8dPhp2OQAQGoIggKx07Vkjtf/QUS17bnvYpQBAaAiCALJSxaj+OmFwX4aHAWQ1giCArGRmuvaskXpua71e3MYcYADZiSAIIGtdOXG4eubl6Ff0CgLIUlkXBM1sr5l5s6+bw64LQNcr6p2vmeOH6rfPbtOBQywaAZB9sioImlmxpFnubo1fkua5++1h1wYgHNeeNVLvHj6m3zy7LexSAKDLRSYImlmlmU1r47NyM5trZtMavyf6Hndf0fSdkpYm+iwAmW/iiGKddFw/3fPXN+XuYZcDAF0qEkEwCHYL2/isVNJ8d6929xXuXi2pKrgeF3eva/LcYkkl7l6baN0AMp+Z6ZNnjdRL2/fp+a0sGgGQXUINgmZWamYLJJVK2tPGbVWSFjS7tkDS/CRff6ukxUk+A0A38PcTh6lXfi5byQDIOnlhvjzojauSJDOb18ZtlWoZBGskLW/8wczmShrbzquWNx0SDkxz97beCSCLFBbk6/LxQ/Xb57brf112sgoL8sMuCQC6RKhBsCPB8G2L3kJ3rzMzmVmpu9cGw8XxPHda82cCyG7XnjVS99Vs0YPrtunT54wOuxwA6BKRmCPYjhLpw3P7mol7nmCgXFJbzwSQhc4YXqRThxayaARAVol0j6Ck4jQ+u8NFIsGQ81xJGjJkiFatWpXGcqQDBw6k/R1IHO0TXalqm8n9j+julw7rvx98RGX9c5MvDJL4uxN1tE90dUXbRD0IpkVn9w0MhpyrJamiosKnTJmSzrK0atUqpfsdSBztE12papuKQ0e15Fsr9MqRgbp+yvjkC4Mk/u5EHe0TXV3RNlEfGpb0wVzBsN4908yq6+vZVgLo7vr2zNPfTxym3z2/XfXvHQm7HABIu6gHwcbh25KmF5sEw7TvAejuy9x9blFRUbpfBSACrj1zpA4dbdD9a7eGXQoApF2kg2CwSKRWLecKlkiqYzNoAKl22rAijR9epF8+vVmHjzaEXQ4ApFWkg2BghaSKZtfKg+tpx9AwkH3++cITVLvrXX1v+WthlwIAaRWlIFii1lcJz5M0q9m1quB62jE0DGSfaacM0bVnjdSCR2v1l/W7wi4HANIm7CPmis1svpktUSwEzjezBWZW2XhPMDw8L7iv0sxuVuzs4S4ZFqZHEMhO37j0FJUN7qt/Xfysdh84FHY5AJAWYR8xV6dO9Oy5+1pJa9NfUavvXiZpWUVFxZww3g8gHL165OpHV0/Ux//fE7p56fO667oKmVnYZQFASkVpaBgAIuWUoYW69WMnaeWrO/XzpzaHXQ4ApBxBsAMMDQPZ7TPnjtbUcYP0rT+8olff2hd2OQCQUgTBDrBYBMhuZqbvzBqvwoJ8fenedTp45FjYJQFAyhAEAaADA/v21Pdnj9frbx/Qt37/StjlAEDKEAQBoBM+euIgzTl/jH7x9GY9/NJbYZcDAClBEOwAcwQBNLrx4nE6dWihbr7/eb1VfzDscgAgaQTBDjBHEECjnnm5+tE1E3XoSIO+svhZHWvwsEsCgKQQBAEgDmMH9dW/X36KntywW9WPcdw5gMxGEASAOM2uGKFLTz9e33v4NT27pS7scgAgYQRBAIiTmen/XnG6hhQW6F8WrdMjr76td/ZzDB2AzBPqEXOZwMxmSppZVlYWdikAIqSod75+cPUEffq/n9Hn7q6RJB1XWKDThxfpjGFFOm14kU4fVqSBfXuGXCkAtAGVQBIAAArySURBVI0g2AHOGgbQlsmjS7T669P08vZ9en5rnV7YVq8XttVr+ctvf3DP0KJYODx9WJGGFBYoN8eUm2PKsQ9/z82RcnNylGumHJPUyrHG1srFTDj++NU9x1RQuzvsMtAG2qfrjSjprWHFvcIuQxJBEACS0rdnns4cU6Izx5R8cG3/wSN6afs+vbC1/oNw+NBLb7fzlCzwzNNhV4D20D5dat4lJ+nzU8aGXYYkgiAApFy/gnydXTpAZ5cO+ODavoNHVP/eER1rcB1zV0Pw/ViDq6FBwZ8bdKxBavCW29K0ckmuzNi+5rlnn9P4CePDLgNtoH263siS3mGX8AGCIAB0gcKCfBUW5IddRigOb8nVuWMHhl0G2kD7ZDdWDXeAk0UAAEB3RRDsACeLAACA7oogCAAAkKUIggAAAFmKIAgAAJClCIIAAABZiiAIAACQpQiCAAAAWYogCAAAkKUIgh1gQ2kAANBdEQQ7wIbSAACguyIIAgAAZCmCIAAAQJYiCAIAAGQpgiAAAECWMncPu4aMYGbvSNqc5tcMlLQrze9A4mif6KJtoo32iTbaJ7qSaZtR7j6oo5sIghFiZjXuXhF2HWgd7RNdtE200T7RRvtEV1e0DUPDAAAAWYogCAAAkKUIgtFSHXYBaBftE120TbTRPtFG+0RX2tuGOYIAAABZih5BAACALEUQBAAAyFJ5YRcAyczKJVVIqpVUKqnW3VeEW1X2MrNKSXWttQFtFZ6gXUoljQ2+L3D3pc3uoX1CYmbTJE2XtFuxNlrj7tXN7qF9QmZmxZLmu3tVs+u0TQiC9pgraamkPZJK/v/27v6ojaSJ4/ivq54AdDwZiAzAjuBEBsKOAJQBqovAJTLAjsAHGcBFYIsMrAwOlEHfH9ML42FlwEiaVe33U0WZfWF3UFvQ9LyspImk6/z132R8SAQrM7Oh0pvyKNt3aWYLd19UbFovxS+zL5KOW44Rq0oiCVw0iV/88Jyb2V6TbBCfeuJ9I3efZvvmZjZw9/PYJj7dMFNKNh4Qm6r2lGIyi+2lpJMiCdxofOgarm8i6aLYd6HH/xTYAjMbmtmF0l9adytOI1b1DN39ttlw96XS657Hg/jUM2nZd1PsJz6VRUKx13KI2NR1JOkPSfvu/kfZ06ENx4dEsL6xpNti3/fYjy1x94W7T8qurAKxqiCqfx/j39xNHB/GNvGp66hl3zL7nPjUN5J03bKf2FTm7stfVPc2Gh8SwYriF9uTClRUO/JfcKiMWNUTr/EwPloRn7rc/TjvFg5jSV8l4tMF0X3/d8t+YtNh24gPiWBde9JjQFvwBuwOYlVRdJeUfxGPlCb1LER8OsXMTiXdNuMDRXy6YLDi9Sc29Q3NbBwfpzEmurHx+DBZpK6yqwvdRay6ZyLpU3xOfDogfoEdSalKmB0iPhWZ2bhl3FmD2NR1J0l5fGIiSLNv4/GhIghg50TF6S6rOKED3P0qliWZxqzhg9pt6rvoWlxVTUJlMTawHJu+1Yk6JIId0DIIHh1FrOqLMTGTfCmF7Bjx6YDoxrqQ9E++n/hU8eEl680Rm05ZKHUXP8Rkk/EhEayrmSFUruk0KI6jPmLVHTNJfxb7iE/33EgaxCQF4lNBVGS/P3MasanIzM5adjcTQ4baQnwYI1iRuy/NbKGnYwD29DgIHh1ArLoh1nqclgOniU89UaGdS/qzZUKPFJMUiE8V7yTtm9nHbN+BUrVpJumbu18RmzqahaLN7Kp4nZukb7GN9w6JYH03Sm/W/AfoQexHtxCrimJc4Cz/wddUm2If8aljoFSVKH8hNbMZm3gQny1rWxc1KlDvi+V+iE0F7r4ws0lLMjdSmnXf/MG70fjQNVzfVE8fZzaJ/ahjT+0ztYhVJdlyCgMzO4iPkaTj7Ico8akgqoBfWw5NJZ0Tn875f8s+YlPPXb4WYHT5TiSdZOdsND7m7uu4Dt4gxnF8lPRN6a/oWx72vV3x5vtL6fUfK1U3bpQe/J1P6ydWWxaxuV9xeOHu+9m5xKeSqNjuS/o3/p2XFSniU08kG1NJH5T+0P0s6aLpzic29cQfukOlJH2goucjztlYfEgEAQAAeoquYQAAgJ4iEQQAAOgpEkEAAICeIhEEAADoKRJBAACAniIRBAAA6CkSQQAAgJ4iEQQAAOgpEkEAO8XMLortcTxu7qVff2Zm1+tv2frsQhtzZjY0s2sz8zI+ALqNRBDArhkW20dKjwR8qfdKD3Xvsl1o4wN3X7j7UWzuTAILgEQQwA6J5w6XSd+wfC7nr7j7sbtbcd1Z/uD3bVl137Y2dl08C1VKz+gGsCNIBAHskpGkebMRieHyLReMa5xJOnju3HWqdd8NGklauPub4gFgu0gEAeySI/1ccXon6VultuBnZWwA7AASQQC7pOwGfnXyEd2xHp9fSrqPQ5cx2eG+OH9sZvM4Ns8npsS17uPzCzP7YWaDmDxxaWb38XXXeRfwc/fN25jtO8gmZMzN7LTl+7pvOe+gOG8c7fS2tr3BSDE+MF6Ds/iYreHaADaERBDATljRDXzg7rdvuOyJUjIpSRNJ+5IOs3uOJV0qJZvHSuMTnyROZjZXqk5eRNfopLj+nrIu7efuW4p2zOP+R5K+SrpomaE7iPZexHWHsd1cp9m+ivs139Ng1b1fIh8fGHH64O7n8X2N33JtAJv1v9oNAIA2keTkCdeepEGxrMq7Ynvm7i+uELr70syaCuNdy6STL5LO3X0a21eR9E0kNfsG8bVNYqfs/OZ7OZE0N7Oxu1+94L6lL5I+u3uTYN6Y2a1SUnpZfM9Td7+K++5LOjOzQSSoTcL2KbZvlZLCtxrpcRLPB3f/3LRlDdcGsEEkggA6KUt6JD0khtNmMkJ00Q6zpGOtoso1UEqkzorDd8X2RL/g7rdmJj1d+uYl7RhFO36q/rn7jZktlap6eSKYf/6juFxz7B8z+yrp5o0V1UaTBOdJoJg4AnQfXcMAdsWwSCzKBGjt94t/D5W6OPOP4/zEtopejMW7jHF69+Xx32hHW9VwoSK5/FXyFceapG2mVKVcx7p/I0V3c9u4RADdRUUQQCcVXcOt3cKShlFpa7yqa/gZD4nXa9YplKRo51Ap2foeu+erv+JF7RgqdeXmhtn1XyRen8MYy3cqaWZmZzGm79WypO9TdHmfKY1D3P+d6wHYLhJBAJ2Udw1HcnHbJHkx6WFadh//prKbt7n/bXS9zvRYRWvaM1hVeYu2jSQdZ2P12iZjtN63RZPoTZR1QccEkmZyyKtF+8/N7KPSk0x+V7l+4EIpcW++71HzOgDoHrqGAeyC90Wl72GpkrfKEphJdOfm1z2RNIou3lF2fOWSKFE9XCpV2prnID9J1p65b3neRNJptGOcVd2uXlMBNbPTWDrmLLvOgR6XfWmWvXlN127bEj75xBGSQKDDSAQB7IKyonao9Y4PPFdKLmfKul8jiTmM+18rzd5d6PnZsCdKVbHLuOaqql3rfUsxAeNIj8vBTJQqoservmaFv5XG8k2y65xnEzwGSsu9vHvFNfeUTWSJ12wRSeZLq54AKjF3f/4sAKikrRvYzK7z5VqwPmb2Q6lbex2ziQF0HBVBAF33UzfwioWlsUYkgUB/kAgC6LpyDNraxgfiZ5Fkswg00CMkggC67q6YoTvUZtcP7C13XzK5A+gXxggCAAD0FBVBAACAniIRBAAA6CkSQQAAgJ4iEQQAAOgpEkEAAICeIhEEAADoqf8Axf+HHGs1VxMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10, 7))\n",
    "plt.plot([np.linalg.norm(expA - approx, \"fro\") / np.linalg.norm(expA) for approx in conv])\n",
    "plt.yscale(\"log\")\n",
    "plt.xticks(fontsize=20)\n",
    "plt.yticks(fontsize=20)\n",
    "plt.xlabel(\"\\# iterations, $k$\", fontsize=20)\n",
    "plt.ylabel(\"$\\|\\exp(A) - S_k \\|_F / \\|\\exp(A)\\|_F$\", fontsize=20)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Что нас ждёт дальше?\n",
    "\n",
    "- Те же задачи, но для матриц большого и очень большого размера\n",
    "- Активное использование структуры матриц\n",
    "    - разреженность\n",
    "    - малоранговость\n",
    "    - Тёплицевы матрицы и циркулянты\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",
    "- Почему эта задача важна?\n",
    "- Основные обозначения\n",
    "- Как задача формально записывается? \n",
    "- Критерий качества"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Методы решения\n",
    "\n",
    "- Какие методы вы используете для её решения?\n",
    "- Что из материалов курса вам было нужно для решения задачи?\n",
    "- Описание алгоритмов"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Результаты\n",
    "\n",
    "- Графики\n",
    "    - Читаемые оси и подписи к осям\n",
    "    - Понятная легенда\n",
    "    - Отличимые линии\n",
    "- Таблицы\n",
    "    - Выделенные ячейки, на которые надо смотреть\n",
    "    - Понятные наименования строк/столбцов\n",
    "- Картинки"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Анонс следующего занятия\n",
    "\n",
    "- 2 лекции\n",
    "- Способы вычисления SVD – последняя тема про плотные неструктурированные матрицы \n",
    "- Разреженные матрицы и прямые методы решения линейных систем с разреженными матрицами"
   ]
  }
 ],
 "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
}