seminars/seminar7/seminar7.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Семинар 7\n",
    "\n",
    "# Введение в итерационные методы"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Степенной метод: напоминание\n",
    "\n",
    "- Задача поиска максимального по модулю собственного значения\n",
    "- Сходимость линейная, коэффициент $\\frac{|\\lambda_2|}{|\\lambda_1|}$\n",
    "- Может не сходится. В каком случае?\n",
    "- Сложность пропорциональна сложности умножения матрицы на вектор"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Сдвиги для ускорения сходимости степенного метода\n",
    "\n",
    "- Метод обратной итерации со сдвигами\n",
    "\n",
    "\\begin{align*}\n",
    "& x_{k+1} = (A - \\sigma I)^{-1} x_k\\\\\n",
    "& x_{k+1} := \\frac{x_{k+1}}{\\|x_{k+1}\\|}\n",
    "\\end{align*}\n",
    "\n",
    "- Куда сойдётся?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "sns.set_context(\"talk\")\n",
    "import scipy.linalg as splin\n",
    "\n",
    "def shift_inverse_iteration(A, x0, num_iter, shift=0.):\n",
    "    x = x0.copy()\n",
    "    conv = [x]\n",
    "    for i in range(num_iter):\n",
    "        x = np.linalg.solve(A - shift * np.eye(n), x)\n",
    "        x = x / np.linalg.norm(x)\n",
    "        conv.append(x)\n",
    "    return x, conv\n",
    "\n",
    "def shift_inverse_iteration_fast(A, x0, num_iter, shift=0.):\n",
    "    x = x0.copy()\n",
    "    conv = [x]\n",
    "    lu_factor = splin.lu_factor(A - shift * np.eye(n))\n",
    "    for i in range(num_iter):\n",
    "        x = splin.lu_solve(lu_factor, x)\n",
    "        x = x / np.linalg.norm(x)\n",
    "        conv.append(x)\n",
    "    return x, conv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.00129089  1.14452777  1.72651117  4.27990214  5.59525073  8.50188201\n",
      " 12.20387858 13.17684228 15.93954198 32.29176751]\n"
     ]
    }
   ],
   "source": [
    "n = 10\n",
    "A = np.random.randn(n, n)\n",
    "A = A.T @ A + np.eye(n)\n",
    "# A = A + A.T\n",
    "true_spectrum = np.linalg.eigvals(A)\n",
    "print(np.sort(true_spectrum))\n",
    "x0 = np.random.randn(n)\n",
    "num_iter = 40"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0012950263476652\n"
     ]
    }
   ],
   "source": [
    "# smallest_eigv, conv_sm = shift_inverse_iteration(A, x0, num_iter)\n",
    "# print(smallest_eigv @ A @ smallest_eigv)\n",
    "smallest_eigv, conv_sm = shift_inverse_iteration_fast(A, x0, num_iter)\n",
    "print(smallest_eigv @ A @ smallest_eigv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1a26fbf208>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAEMCAYAAAAWDss+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl4XPV97/H3V4vlRbJkbRgwsfEmmZ2AYxYbyyR5GiCscZNQUkjapje5JTe3TUO6JLeE55Z7S8u9NyQ3t4WSQto8acliCMQ0oWBhbLMECBisWPIONgYttmRbtmwt3/vHOSNG45E8kmfmjEaf1/PMM+NzfjPz08/Lx7/zW465OyIiIplQEHUFREQkfylkREQkYxQyIiKSMQoZERHJGIWMiIhkjEJGREQyRiEjIiIZo5AREZGMUciIiEjGKGRERCRjFDIiIpIxChkREckYhYyIiGSMQkZERDKmKOoKZIuZ9RGE6oGo6yIiMo5MBwbcfUx5kZchY2aNSQ4XApSXl5dntzYiIuNXV1cXnMRVr7wMmWH0l5eXF3Z2dkZdDxGRcaOiooKurq4xXwHKy5Bx94bEY2bWCagXIyKSRRr4FxGRjFHIiIhIxihkREQkY3I6ZMysxMz+ycx2m1mXma0xs7OjrpeIiKQmp0OGYGLCduASoBJ4HHg00hqJiEjKzN2jrkPKzGwS0APUuHvHKN/bWV5eXj7aKczuzi82vUvbwaOsqK9l1oypo3q/iMh4Fk5h7nL3irG8P609GTObZWbfMrN1ZnbIzNzMGoYpW2pm95nZXjM7YmYvm9l1J/iKy4DW0QbMyTAzvvrjjXzjsU28sbsrW18rIpIX0n25bD5wM3AIePoEZVcBtwBfB64BmoBVZnZ1ssJmVgH8A/AXaattimrKSgBoP3Q0218tIjKupXsx5lp3rwUwsxuApD2TMEg+Atzk7qvCY2uAucC9wOqE8pMJxmKecPfvpbnOJ1RdWsL2tm7aDipkRERGI60h4+4DKRa9EegCHot7r5vZw8D9ZnaWuzcBmFkh8K/AHuBPh/vAcEX/SMa82j/Wk2k7dGysHyEiMiFFNbvsHKApSShtjDsf8wAwBfisRzRLoaZUl8tERMYiqr3LqoCWJMf3xZ3HzGYDnyOYUbbfzGLlrnL35+LfGD/zYZhdmJcS7sQ8WtWlkwB0uUxEZJSi3CBzpF6JA7j7LsBGKJcVGvgXERmbqC6XdRD2VhJUhs/7kpyLTHV4uazt4FHG07oiEZGoRRUym4BFZpb4/eeGz29muT4jivVkjvYNcOhoX8S1EREZP6IKmVVABXBtwvFbgebYzLJcEevJALRrhpmISMrSPiZjZivDl4vD5+VmVg10u/uT4bHVwBrgQTOrAnYAtxEMzl+f7jqdrKpw4B+CS2ZnVk+LsDYiIuNHJgb+f5Tw6zvD513AHBhcE3MDcHf4qCBY8X+Tuz+egTqdlJKiQsqnFNN1pFeD/yIio5D2kHH3lGaDufsB4PbwkfOqSyfRdaRX05hFREYhyinMGePuDYnHwh0BTmrV/7a2bvVkRERGIS9DZpjFmKUn85nVWvUvIjJquX7TspwRv1ZGRERSk5c9mUxdLgNtkikiMhrqyaRocJNM9WRERFKWlz2ZTIzJvN+TCbaWidusU0REhqGeTIpiYzLH+gY4qK1lRERSkpc9mUyMyVSXDV31P31y8Vg/SkRkwlBPJkVV0+L2L9O4jIhIShQyKZpUVEDF1KD30qa1MiIiKcnLy2WZGPiHYIZZ5+Fe9WRERFKknswovL/qX2tlRERSkZc9mUwM/ANUl2nVv4jIaKgnMwo12r9MRGRUFDKjEJvGrIF/EZHUKGRGQVvLiIiMTl6OyWRqdllsTKb90DFtLSMikgL1ZEYh1pM51j/AgSPaWkZE5ETysieTqdllsU0yIRiXKZ+qrWVEREainswoVE4bun+ZiIiMTCEzCsWFBYNBo2nMIiInppAZpepShYyISKpyPmTM7Jtm1mRmA2b26ajrE9taRpfLREROLOdDBtgCfBl4KeqKwPuD/+rJiIicWM6HjLv/i7s/BfREXRdQT0ZEZDTSGjJmNsvMvmVm68zskJm5mTUMU7bUzO4zs71mdsTMXjaz69JZn0yoKdNOzCIiqUp3T2Y+cDNwCHj6BGVXAbcAXweuAZqAVWZ2dZrrlFbV2iRTRCRl6V6MudbdawHM7AYgac8kDJKPADe5+6rw2BpgLnAvsDrN9Uqb+Nll2lpGRGRkae3JuPtAikVvBLqAx+Le68DDQL2ZnTXa7zazzpEenORq/5jY5bLefqfrSG86PlJEJG9FNfB/DtCUJJQ2xp0HwMyKzWwyQV2LzWyymUU2YSG2fxlo8F9E5ESi+se6CtiX5Pi+uPMxDwBHgGXA98PXVyS+0d0rYg/gtSSP/nRUvHLaJGJXyHRfGRGRkUU5hdlTOefun3V3S3g0Zr56yRUVFlA5NTYuoxlmIiIjiSpkOhjaW4mpDJ+T9XJyhtbKiIikJqqQ2QQsSjK2cm74/GaW6zMqWvUvIpKaqEJmFVABXJtw/Fag2d2bsl+l1MWmMasnIyIysrTftMzMVoYvF4fPy82sGuh29yfDY6uBNcCDZlYF7ABuA5YC16e7TummnoyISGoycWfMHyX8+s7weRcwB4I1MeFizbvDRwXBiv+b3P3xDNQprbTqX0QkNWkPGXdPaQm8ux8Abg8f44oG/kVEUpOJnkzk3L0h8VgmVv13HDrGwIBTUKCtZUREksnLkDGzxiSHS9P1+bGeTN+A03mkd/CWzCIiMlTO308mF8V6MqBxGRGRkeRlTybTl8sqp02iwGDAof3gURaeUpaOjxURyTvqyYxBYYENXiLT/mUiIsPLy55MpsdkIBiXaT90TDPMRERGoJ7MGMXGZdSTEREZXl72ZDI9JgPv31em/aB2YhYRGY56MmNUra1lREROSCEzRtokU0TkxPLyclk2Bv61SaaIyImpJzNGsVX/Hd3B1jIiInK8vOzJZGXgP+zJ9A84+w8fo6q05ATvEBGZeNSTGaPq0vitZTTDTEQkGYXMGM2YGmwtAxr8FxEZjkJmjAoLbPASmQb/RUSSy8sxmWzMLoPgklnbwaPqyYiIDEM9mZOgacwiIiPLy55MNmaXQdyCTIWMiEhS6smchNj+ZbpcJiKSnELmJLx/uUxTmEVEklHInIRq9WREREaU8yFjZqeY2S/N7LCZbTSzi6KuU0ysJ7Ov+yj92lpGROQ4OR8ywP8DtgJVwHeAn5hZcbRVCsR6MgMO+w/rkpmISKKcDhkzKwM+DnzT3Y+4+/3hqcsirNag2Owy0CUzEZFk0hoyZjbLzL5lZuvM7JCZuZk1DFO21MzuM7O9ZnbEzF42s+sSii0AOt39vbhjbwBnpbPeYzVj6iQKw71ltFZGROR46e7JzAduBg4BT5+g7CrgFuDrwDVAE7DKzK6OKzMNOJDwvgNkYPX+WBQUGFXTdPMyEZHhpHsx5lp3rwUwsxuAxJ4J4bmrgY8AN7n7qvDYGmAucC+wOizaDZQlvH06QYjlhJqyEloPHlVPRkQkibT2ZNx9IMWiNwJdwGNx73XgYaDezGKXw7YAM8zslLj3nkPQ6xnCzDpHepDm1f4x1aVaKyMiMpyoBv7PAZqShNLGuPO4+0HgCeAbZjbZzP4AMGBD1mp6AlorIyIyvKj2LqsCWpIc3xd3PuaLwL8A+wmmMn/C3XsT3+juFbHXw+zCvBQoHGN9h6VNMkVEhhflBpkjrV4cPBfOLPto5qszNoObZKonIyJynKgul3UwtLcSUxk+70tyLifFejLvHeghGFYSEZGYqEJmE7DIzBK//9zw+c0s12fMzj5tOgD7D/eybmt7xLUREcktUYXMKqACuDbh+K1As7sfN3ssV82vLWPJmUEH7KH1O6OtjIhIjkn7mIyZrQxfLg6fl5tZNdDt7k+Gx1YDa4AHzawK2AHcRjA4f32665Rpn7t8Di/u2Mczza3sbO9mTvW0qKskIpITMjHw/6OEX98ZPu8C5kCwJiZcrHl3+KggWPtyk7s/noE6ZdRHFp3C6RVT2NN5hO8/v4v/dm1O7HojIhK5tF8uc3cb5jEnodwBd7/d3We6+2R3/6C7P5ru+mRDUWEBv3vpbAB+9PLbHDraF3GNRERyQ07vwjxW7t6Q+CDDW9F8evEZTC4u4ODRPn766u5MfpWIyLiRlyFjZo2JDzK8qWbF1EnceOHpADy0YScDuomZiEh+hkxUbrtsDgDb27p5TtOZRUQiXfGfMeHlsSEyuUlmTP3M6Vw6t4rnt3fw8IadLF9Yk8mvExHJeerJpFmsN/PM5lZ2tHdHWxkRkYjlZchEMSYT85FFtZxeMQWA7z+/MxtfKSKSs/IyZKJUVFjArYPTmXdrOrOITGh5GTJRTGGO96lwOvOho3385BVNZxaRiSsvQyZqwXTmWQA8rOnMIjKBKWQy5LOx6czt3azd0hZtZUREIpKXU5iHuTNmVgb+Y+pmlnHZvCo2bOvgoQ07aairzebXi4jkBPVkMig2nbmxuY3tbVkbEhIRyRl5GTJRD/zHxHZnBrjz8Sb6+geyXQURkUjlZcjkisIC446P1QGwtqWNu1dvjrhGIiLZpZDJsOsvOJ0/vGIuAN9bv4MfvvRWxDUSEckehUwWfO1j9Xy4Phj4/8ajb/L8to6IayQikh0KmSwoLDC+dfOF1J1SRt+A88UfvMJO7WsmIhNAXoZMlHuXDae0pIh/vO1iqqZNovNwL7//8K/oOtIbZZVERDIuL0MmV51ROZV/+N2LmFRYwLa2br70w19rxpmI5LW8DJlcmcKczMVzKrn7pnOBYMbZf//5byKukYhI5uRlyOS6lRfN4gvL5wHBrZp/8OKuiGskIpIZCpmI3PFbdXz0rFOAYMbZA2u3466NNEUkv+R0yJjZN82sycwGzOzTUdcnnQoKjP/zqQs4/4wKBhz+evVv+ON/e42e3v6oqyYikjY5HTLAFuDLwEtRVyQTppUU8W9/eAk3Xng6AI++9g4r/34DezqPRFwzEZH0yOmQcfd/cfengJ6o65Ipk4sL+V+fPJ+vX7OIAoM39xzgum+v44XtWrApIuNfSiFjZrPM7Ftmts7MDpmZm1nDMGVLzew+M9trZkfM7GUzuy6ttc4zZsYfLJvL939vCRVTi+noPsZn/vFFHt6wU+M0IjKupdqTmQ/cTDAN+OkTlF0F3AJ8HbgGaAJWmdnVY63kRLF0QTU/+6Ol1M8Mdgb4q59t4ms/2ahxGhEZt1INmbXuXuvuHwP+abhCYZB8BPgDd3/Q3Z8BbgOeB+5NKPtk2CtK9rhjjD/PuPeBqqn89D9fxjXnngrAIy/v5tpvr+NFXT4TkXEopZBx91SXpd8IdAGPxb3XgYeBejM7K+74Ve5eOszjnlH8DACYWedID6B8tJ8ZlamTivjO71zIV3+rjgKDLa2H+NT9L/Anj7xG28GjUVdPRCRl6R74PwdoShJKG+POp8zMis1sMkE9i81sspnl9GSFdDEz/mjFfH52+1LOP6MCgJ++uocP39vIP7+wi/4BjdWISO5L9z/YVcC+JMf3xZ0fjQeAI8Ay4Pvh6yuSFXT3itgDeC3JY1wObJxzejmrvngZd994LuVTijnQ08c3Hn2TG7+7ntff7oy6eiIiI8pEr2Ck/2KP6r/f7v5Zd7eER+PJVW/8KSgwfmfJB3jmK8v55MWzANi4u4sbvruev1z1Bu9oXY2I5Kh0h0wHyXsrleFzsl6OpKiqtIR7Vp7Pj79wKfUzy3CHH7z4FsvuWcN/+eGv2bhbPRsRyS3pDplNwKIk4ybnhs9vpvn7JqSL51TyxJeW8o2Pn0VNWQn9A87PXn+H676znk/+/fP8YtO7GrMRkZyQ7pBZBVQA1yYcvxVodvemNH/fhFVUWMDvLz2TdV9bwd/99vnUzywD4KWd+/hP//wKV97byMMbdtJ9tC/imorIRGaprig3s5Xhy8XAHcCdBD2Xbnd/MixjBIs1zwvL7CBYJ3MrcL27P57Oyo9Q18Ykh5eWl5cXdnbm5yUld2fDtg7+8bntrGluGzw+ubiAFXW1XH3uqVxZX8u0kqIIayki401FRQVdXV1d4aSqURtNyAxXcJe7z4krNx24G1hJ0KtpAu5y90fHUsGxmIghE29r60EeXLeTn766m6N9788mn1xcQMPCWq4+71Q+rMARkRRkLWTGk4keMjFdh3t56jfvsfqNvTy3pY3e/vd/r0uKCmioq+HK+lqWLajhtIopEdZURHKVQiZFZtZZXl5ePpFCJl7XkV7+oykWOO0c6x+6XnZ+bSnLFlRzxcIalpxZydRJ6uWIiEImKfVkRnagJwicX256j/Xb2jnYM3RywKTCAi6eM4NlC2q4dF4V55w2naLCCbHRgogkUMgkoZBJXV//AK/v7mRtSzvPbWnjtbc7SZz9XFZSxOIzK7l0bhWXzqti0anTKSywaCosIlmlkEnRRL9clqquw71s2NbO2i1tbNjWwa6Ow8eVmT65iCVzq7h0bhWXzK2ifmYZBQodkbykkEmRQmZs9nQe4fltHTy/rYMXtnckvTV0xdRilpxZySVh6NSdotARyRcKmSR0uSwz3J239x3h+e3tYejs490Dx98Ze8bUYpacGVxau3ReFQtqSwmWUInIeKOQSUIhkx3uzq6Ow7ywvYPntwe9ndYk97upLi3hkrmVQejMreLM6mkKHZFxQiGTIl0uyzx3Z2cYOhvCS2zth44PnZnTJ3PpvCoum1fF5fOrtUZHJIcpZFKkkMk+d2db26FgTCfs6ew/3HtcuTlVU7lsfjWXz6vm0nlVVE6bFEFtRSQZhUyKFDLRGxhwmt87yPPbgp7Oi9s7OJhkA89Fp05n6fwqli6o4UNzKpkyqTCC2ooIKGSS0pjM+NDXP8Abe7rYsK2DDdva+dXO/RzrG7oTwaTCAj44u4Kl86u5fH41555eroWhIlmkkElCITM+9fT28+qu/azf1s66rR28sTvJwtDJRVw6t4plC2tYNr+a2VVTNYlAJIMUMinS5bLxp+twL89v72D91nbWb21ne3v3cWXOqJzC0vk1XLGgmsvmVVM+tTiCmorkL4VMihQy49+eziOs39LOuq3BY1/3sSHnCwzOm1UxuNHnhWdU6NKayElSyKRIIZNfBgacpr0HeG5LsOfayzv3H7ezdFlJEZfNr2LZghqWL6zhjMqpEdVWZPxSyKRIIZPfjhzr58UdHYOh0/LeoePKnFk9jWULqlm+sIZL5lbppm0iKVDIpEghM7Hs7TrCc1vaWdvSxrqt7XQmrM8pLjQunl3J8roarlhQw6JTyzSBQCQJhUwSml0m8foHnDf3dLG2pY21W9p49a1O+hOmrdWUlXDFghquWFjNFQtqmKEFoSKAQiYphYyM5EBPLxu2tvNsS9DTSdxZ2gzOn1XB8oU1LK+r4fxZFbp/jkxYCpkU6XKZJBNsfdPN2pY2nm1p44XtHRxNWBBaPqV4cCxn+cIaaqdPjqi2ItmnkEmRQkZS0dPbz0s79vFsGDpbW4+fQHDWqdNpqAsC54OzZ1CsadKSxxQyKVLIyFjs3n+YtS3tPNvSyvqtHRxK2GutrKSIy+dXB6FTV8Op5dpRWvJL3oaMmZUAfw98FCgDXgVud/dNY/w8hYyclGN9A7z61n4am4Nezm/2HjiuTN0pZYOBc/HsSiYVqZcj41s+h8w04E+AfwL2Al8GvujuC8b4eQoZSav3DvTwbBg4a7e0cbBnaC+ntKSIy+ZVsaK+lgb1cmScytuQSWRmk4AeoMbdO8bwfoWMZExf/wC/fruTxuZWGpvb2PSOejmSH7ISMmY2C/gqcBFwATANWOHujUnKlgJ3A78NVACbgLvc/WdjqWDc5zYA/+ruM8f4foWMZE3rgR6ebWmjsaWN51raOJCkl3P5/Coa6tTLkdyWrZBpAB4hGBc5ClzH8CHzFPBB4A5gB/BZ4BbgWndfPaZKmlUALwJ/4+7fG+NnKGQkEqn0cupnlrG8roYVdbVcpBlrkkOyFTIF7j4Qvr4BWEWSkDGzq4GfAze5+6rwmAHPAVXuviiu7JPAsmG+8i53vycsNxn4d+AVd//K6H68IXVTyEhOaD3QQ2NLG882Jx/LKSspYumCalbU1bK8roZTtC5HIpT1MZkThMwDBJfJKmOhFB7/PHA/cLa7N43iuwqBnwDdwGd8hMqa2YnSo7y8vByFjOSSvv4BXn3r/V5OU5IZa2edOp0V9UEv5wLdvkCyLNdC5nnA3f2yhONLgBeAT7n7I6P4ru8BpwMfd/feE5RVyMi4F5uxtqa5lXVb2jmYsC5n+uQirlhYM9jLqS4tiaimMlGcbMike6/zKqAlyfF9cedTYmazgc8RzCjbH7dD7lXu/lxi+fgGGG7vMqAw1e8XicIp0yfzycVn8MnFZ9DbP8Aru/azprmVZ5vb2PzuQQ709PHExr08sXEvAOfNKqehrpYVdTWcpz3WJAdl4oYaI3WNUu42ufsuQH9jZMIqLizgkrlVXDK3ij+/ahHvdB7h2ZY21mxuZd3Wdg4f62fj7i427u7ivqe3MGNqMcsX1rCivlY7SUvOSHfIdJC8t1IZPu9Lck5EUnBaxRRu/tAHuPlDH+BoXz8v79xPY3Mrz2xuZVtbN/sP9/Loa+/w6GvvUGBwwRkVrKirZUV9LWedOp0C9XIkAukOmU3AJ+Jno4XODZ/fTPP3iUxIJUWFXD6/msvnV/OX15zF2/sO09jcyprmNjZsa6enN5hQ8Opbndz7VAs1ZSU0hL2cpQuqmT65OOofQSaIdA/8XwM8Adzg7o/FHV8L1Lp7/UnXOLU6NiY5rPvJyITQ09vPizv2sWZzK2uaW9nVcXjI+cIC4+LZM1hRX8uKuloWnlKqu4LKsLI2u8zMVoYvFxMstLyToOfS7e5PhmUMeBo4j/cXY94G3Apc7+6Pj6WSo6WQEXnfjvbuwcB5cfs+jvUPvV/OaeWTaQgD57J5VUwrycRQrYxX2QyZ4Qrucvc5ceWmE2wrs5JgW5kmgsWVj46lgmOhkBFJ7vCxPjZs7WBNuC4n8a6gkwoLWDK3cnDG2tya0ohqKrliwmyQORoKGZETc3da3jsUjuW08vLO/fQNDP33YHbV1MHJA0vOrGRysVYBTDQKmRRpWxmRkR3o6WXdlnbWbG6lsaWNtoNHh5yfXFzA5fOqw0trNcyaMTWimko2KWSSUE9G5OQMDDhNew/wzOZWGptb+fXbnST+U7GgtnTwXjmL51RqU888pZBJQiEjkl77uo/x3JZgIeizLW3sPzx0l6fSkiKWhZt6NtTVUKtNPfOGQiZFulwmkh79A85rb3eGl9VaeXPP8Zt6nn3a9HAsp4YLzpih7W7GMYVMihQyIpkRu3VBY3Mrz7Ucv6lnRWy7m7parlhYQ6W2uxlXFDJJ6HKZSDTiN/Vcs7mVlvcODTlvBhdqu5txRSGThEJGJDfs3n+YxuZgLGd9uN1NvJqyElaEdwS9XNvd5CSFTIp0uUwkWifa7qaowLh4zgyuDHcfmF+r7W5ygUImRQoZkdzh7sF2N2Ev58UdHfT2D/236PSKKayor+HK+lounVvNlElaCBoFhUyKFDIiuav7aB/rt7azpjmYQLC3q2fI+UlFBVw6t4oVdTVcWX8KH6jSQtBsUcgkoTEZkfHL3dn87sFgf7XNbbzy1n76E7a7mVszLZg8UFfL4jNnUFKkXk6mKGSSUMiI5I+uw72s3dI2eBvqju5jQ85PmxTcWye2+8Cp5VMiqml+UsikSJfLRMa/gQHnjT1dg9vdvL6767gyi06dHsxYq6/lwjMqKNJ2NydFIZMihYxI/mk7eJRnW4JeztqWNg72DF0IWj6lmCsW1rCiroblC2uoKi2JqKbjl0ImRQoZkfzW1x/ccjrWy9n87sEh583g/FnBQtAr62s5+zQtBE2FQiZFChmRieWdziPhzgNtrN/azpHe/iHnq0tLaKgLpkgv1ULQYSlkUqSQEZm4enr7eWnHvsHtbnYmWQh60ewZrKgPejkLtBB0kEImCc0uE5GR7GjvHtx54MXt+zjWP3S7m9Mrpgz2ci6dV8XUSUUR1TR6CpkkFDIikqruo31s2NYx2MsZaSHoivpaZldNi6im0VDIpEiXy0TkRNyd5vcOsmZzMGPtlV1aCKqQSZFCRkRGq+twL89tbeOZzSdeCLqirpaZ5fl3R1CFTIoUMiJyMgYGnI17uoI7gk6ghaB5HTJmdj9wLTAN2AX8ubs/McbPUsiISNq0HTzK2pY2nsnzhaD5HjL1wA53P2pmi4GngDnuPuqkUMiISKakuhA0dq+c8bQQNK9DJp6ZXQhsAC52901jeL9CRkSyIp8WgmYlZMxsFvBV4CLgAoLLVyvcvTFJ2VLgbuC3gQpgE3CXu/9sTBU0+y7wOWAy8HPgWh9DMipkRCQKsYWgsV7OeFsImq2QaQAeAV4FjgLXMXzIPAV8ELgD2AF8FriFIBxWj6mSZoXACuAsd79vjJ+hkBGRyI23haDZCpkCdx8IX98ArCJJyJjZ1QS9jZvcfVV4zIDngCp3XxRX9klg2TBfeZe735OkHk8A33b3X6TwsyW+VyEjIjklfiFo4+ZW3snBhaBZH5M5Qcg8QHCZrDIWSuHxzwP3A2e7e9NYKhp+zpPA4+7+3STnTpQe5eXl5ShkRCQXuTst7x3imbCXkysLQU82ZNLdDzsHaIoPmNDG+POpfFA4tnMD8CjQA9xIcMnsa+mpqohI7jAz6maWUTezjC82zBtcCLpmcxvPtrTSfugY29u62d62gwfX7Rg3dwRNd8hUAS1Jju+LO58qB34P+A5gwFbgZnffmLRwXMoOt3cZkN/7P4hI3iifWszHzzuNj593WtKFoN3H+vll03v8suk9AOpnlgVTpHNsIWgmRpRGuv6W8rU5d+8Grjz56oiIjG8FBcYFZ1RwwRkV/PFHFya9I+jmdw+y+d2DfLdxW04tBE13yHSQvLdSGT7vS3JORERGoaashJUXzWLlRbPo7R/g1V37WdPcNrgQtOtIL4+//g6Pv/7O4ELQ3zp7Jl9smJf1uqY7ZDYBn4ifjRZHg029AAAGaUlEQVQ6N3x+M83fJyIyoRUXFrBkbhVL5lbxZ1fVs6fzCI0JC0Ffe7uTsslFeREyq4DfJ9hv7LG447cCzSczs0xERE7s9Iop3LJkNrcsmU1Pbz+/2hksBD1vVnkk9Uk5ZMxsZfhycfi83MyqgW53fzI8thpYAzxoZlUEizFvIxh0vz49VRYRkVRMLi5k2YIali2oiawOKa+TMbPhCu5y9zlx5aYTbCuzkmBbmSaCxZWPnlxVT44WY4qIjF7W1sm4e0qb6bj7AeD28BGJYaYwl2a7HiIiE11uTKQWEZG8FO3Oaxni7g2Jx8JtZ6IZ+RIRmaDyMmR0uUxEJDfocpmIiGTMuLkz5skyswHAyst1xUxEJFVdXV0A7u5j6pRMpJDpI+i5HRjD22PJ1JW+GuU1tdfoqL1GR+01eifTZtOBAXcf0/DKhAmZkxG7V81Y54lPNGqv0VF7jY7aa/SibDONyYiISMYoZEREJGMUMiIikjEKGRERyRiFjIiIZIxCRkREMkYhIyIiGaN1MiIikjHqyYiISMYoZEREJGMUMiIikjEKmRGYWamZ3Wdme83siJm9bGbXRV2vXGBms8zsW2a2zswOmZmbWcMwZX/HzF43sx4z221m/9PMJme5ypExsw+b2UNm1mxmh8M2+KmZnZuk7EfN7IXwz1urmf2DmU2oPbrM7DIz+4WZ7Qn/zLSZ2TNmdlWSshO+vZIxszvDv5OvJTmX1TZTyIxsFXAL8HXgGqAJWGVmV0daq9wwH7gZOAQ8PVwhM/sM8ANgPXAVcDfwR8BDma9izvgC8AHgfxO0wZ+Ev/6VmV0SKxSG9GrgbeBa4E+B64Cfm9lE+rs6A2gGvgJ8DPhD4Ciw2sw+HSuk9krOzM4Gvga8l+RcA9luM3fXI8kDuBpw4Ma4YwasA34Tdf2ifgAFca9vCNuqIaFMIbAXeCzh+OfD8kui/jmy1Fa1SY5VAPuBn8Qdewn4dULbfjRsq09F/XNE3IZF4T+Mz6i9RmynAuAF4NtAI/Bawvmst9mETfsU3Ehw74XHYgc8+B15GKg3s7OiqlgucPeBFIpdAswkaLN4PwB6gU+ku165yN1bkxzrBLYAswDM7HRgMfDP8W3r7k8Be5ggbTUcd+8j+PvYC2qvEfwxwZ+pv0w8EVWbKWSGdw7QlOQf041x52VksTZ6M/6gux8GtjGB29DMagh+/ljbJG2r0BtMwLYyswIzKzKz08zsm8BCgkuOoPY6jpnNBe4Cbnf3ZDdnjKTNFDLDqwL2JTm+L+68jCzWRsO144RsQzMz4H6Cv39/Fx5WWx3vEYKeyx7gvwKfdPd/D8+pveKEf6YeAH7h7o8OUyySNlPIjGyk7RC0VULqhmuridqGf0swjvUFd/9Nwjm11fvuAD5EMDC9GnjEzG5OKKP2CnweuBj4Ugpls9pmY7pn8wTRQfJkrwyfk/1vQIbqCJ+r4l7HVAI7slud6JnZXxPMmvqyuz8Udyq+rRJVMgH/vLn7dmB7+MvHzexx4P+a2b+h9hpkZtXAPcD/ALrjpiMXAYXhr3uIqM3UkxneJmBRkml9sbUNya5rylCbwuch13rNbCowjwnWhmZ2F/AXwB3ufl/C6aRtFTqXCdZWw3iJYHpzDWqveLOAcoKQ2R/3uJygffYDdxJRmylkhreKYJrptQnHbwWa3b0p+1Uad14A3gV+N+H4zUAx8NOs1ygiZvZXwDeAb7j73yaed/fdwMvALfH/sTGzDwOnM4HaKplwzKEB6AQ61F5DbAVWJHm8TjDBZgVwf1Rtpstlw1sNrAEeNLMqgks7twFLgeujrFiuMLOV4cvF4fPysOve7e5Punufmf0Z8JCZfQf4MbAI+Bvgx+7+QvZrnX1m9hWC/0k+AfxH/AJM4Ki7/zp8/TXgl8APzex+4DSCtnoR+FH2ahwtM/sBsAt4BWgHTiX4u3cl8KVwOjOovQBw90MEa2KGMLPO8Hz8uey3WdSLh3L5AUwHvkPwv/Ee4FXghqjrlSsPgoHCZI+dCeU+QzBF8ijBTKF7gClR1z+L7dQ4irb6WPgXvgdoI5gxNCPqnyHL7XU78DzBGEJf+PwL4NokZSd8e43Qjo0kLMaMos10PxkREckYjcmIiEjGKGRERCRjFDIiIpIxChkREckYhYyIiGSMQkZERDJGISMiIhmjkBERkYxRyIiISMb8fzJ03R3R7XRSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy([np.linalg.norm(A @ x - (x @ A @ x) * x) for x in conv_sm])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.7309762959740826\n"
     ]
    }
   ],
   "source": [
    "shift = 3\n",
    "middle_eigv, conv_m = shift_inverse_iteration(A, x0, num_iter, shift)\n",
    "print(middle_eigv @ A @ middle_eigv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1a273a4860>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEMCAYAAADqG+D0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAF6NJREFUeJzt3X+QbGdd5/H3t+fHnZl7yUzuDaghBVlkKQmkLEtZKWSXuEIVJAIhUrtm+WlRWKhYriuFlgu1wNayri4qiGUZihXUlD8Qr1TYy2KERIyIGBECuWzCjxAJxkDuzZ279+fcmX72j3POzJkzp7tP93RP306/X1Vd3XP69DlPP/fHZ57nOed5IqWEJEn9ao27AJKkyWSASJIGYoBIkgZigEiSBmKASJIGYoBIkgZigEiSBmKASJIGYoBIkgZigEiSBmKASJIGYoBIkgZigEiSBmKASJIGMjvuAuxWRKyTBeHJcZdFkibIJUA7pTRwDsSkrwcSEW0glpeXx10USZoYq6urACmlNHBP1MS3QICTy8vLyydOnBh3OSRpYqysrLC6urqrnhvHQCRJAzFAJEkDMUAkSQMxQCRJAzFAJEkDMUAkSQN5NFzGO7DP/OMjfPmhUzz+0kV+4MmXjbs4kjRRproF8oE7H+CNH7yLP/j0P467KJI0caY6QBbnZgA4d2FjzCWRpMkzlgCJiH0R8TsR8UBErEbEbRHxtL0ux+J89vXPrBkgktSvcbVAZoGvAs8EDgK3AH+214VYms+GgM7aApGkvo0lQFJKp1NK/zWl9EBKaQN4N/CdEXFoL8uxkHdhnbUFIkl9axQgEXFFRLwzIu6IiFMRkSLimg77HoiId0XEgxFxNiLujIgX9TjFs4BvppSO9Vn+XXEMRJIG17QF8mTgRuAU8LEe+x4GXga8CbgOOAocjohr63aOiBXgt4FfbFiWoSnGQOzCkqT+Nb0P5BMppccBRMT1QG2LIg+J5wI3pJQO59tuA54EvAM4Utl/gWzs48Mppf810DfYhcW57Os7iC5J/WsUICmldsPjvQRYBT5U+myKiPcDN0XEVSmlowARMQP8IfAN4A2dDhgRvRb6GHglqcV5u7AkaVDDHkR/OnC0JnDuKr1feA+wCLw6jWlZxGIM5MJG4sJG04yUJMHwpzI5BNxbs/146X0i4onAjwHngEciotjvBSmlvyp/MKW00u2EeQtloFZIESCQtULmZqb6vkpJ6sso5sLq1ppIACml+4Host+eKAbRIbuU9zELc2MsjSRNlmEHyDHyVkbFwfz5eM17jUXE7TWbDwx6vMX5ra/vlViS1J9h99ncDTw1IqrHvTp//sKQz7cr5S4sA0SS+jPsFshh4DXACyldiQW8ErinuAJrUCmla6rbhjUG4t3oktSfxgESES/NXz4jf35ORFwGnE4pfSTfdgS4DXhvPi3JfcCrgGcDLx5OkYdn32xpDMQWiCT1pZ8WyAcqP78lf74fuBI27/m4Hnh7/lghuxP9hpTSLbsq6Qi0WsHi3AxnL2zYApGkPjUOkJRSo6umUkongdfnj6Ea9iA6ZDcTnr2wYQtEkvo09Tc+LDojryQNZKLWRB/2IDrAwlyWoU5nIkn9sQWSz4flhIqS1J+pD5ClOVcllKRBTH2ALOQtEANEkvozUWMgI7kKqxgDsQtLkvoy9S2QzauwbIFIUl8mqgUyiquwigkVHUSXpP7YAplzVUJJGoQBkq8JYheWJPXHAPFOdEkayNQHyMKcNxJK0iAmahB9FJfxLuWD6I6BSFJ/pr4F4hiIJA1mologI7mM1zEQSRrI1LdAFjYv422PuSSSNFmmPkCKMZC1jTbrG4aIJDU19QFSdGGB4yCS1A8DZH6rCgwQSWpu6gNkodQCObdmF5YkNTVRV2GNZjr3rQA5c2F9N4eSpKky9S2QYhAdvJRXkvoxUS2QUdwHsm/WMRBJGsTUt0BarWChWJXQAJGkxqY+QKB8N7qD6JLUlAHC1jjImTUH0SWpKQME7MKSpAEYIMDifN6FZYBIUmMGCI6BSNIgJuoy3lHcSAilVQm9kVCSGrMFAizlXVjnvJFQkhqbqBbIKG4khFIXlmMgktSYLRDKg+iOgUhSUwYIW2MgzoUlSc0ZIGyNgZx1EF2SGjNAKF/GawtEkpoyQCh1YTkGIkmNGSBsDaI7lYkkNWeAsNWF5WSKktScAUJpEN0xEElqzABhawzknGMgktSYAcJWF9baRpv1DUNEkpqYqKlMRjWZYjGIDnBuvc2BGXNVknrxf0q2xkDAgXRJamqiWiCjmkyxGAMBOOeaIJLUiC0QtsZAwBl5JakpA4TtYyAGiCQ1Y4AAC7OOgUhSvwwQoNUKFuayqnA6E0lqxgDJbc3I6yC6JDVhgORc1laS+mOA5BbmDRBJ6ocBktuaUNFBdElqwgDJOQYiSf0xQHILjoFIUl8MkNzinKsSSlI/DJBcMQbijYSS1IwBklvcvArLMRBJasIAyW2OgbisrSQ1MlHTuY9qQSlwDESS+mULJOed6JLUn4lqgYxqQSnYGgM5YxeWJDViCyRXBIhdWJLUjAGSW3QQXZL6YoDkHAORpP4YILnN+0BsgUhSIwZIrmiBrG20Wd/wZkJJ6sUAyRUtEIBz6waIJPVigOSKFgjYjSVJTRgguYVSgHgpryT1ZoDklkpdWN5MKEm9GSC58hiIl/JKUm8GSG5h1jEQSeqHAZJrtYJ9s1l1OAYiSb0ZICVLTqgoSY0ZICVOZyJJzRkgJQvzBogkNWWAlGyuSmgXliT1ZICU2IUlSc0ZICWuSihJzRkgJZtdWLZAJKknA6TENUEkqTkDpMQxEElqzgApcQxEkpobW4BExFsj4mhEtCPiR8dVjjLHQCSpuXG2QL4E/Azw6TGWYRu7sCSpubEFSErp91NKtwLnxlWGKgfRJam5RgESEVdExDsj4o6IOBURKSKu6bDvgYh4V0Q8GBFnI+LOiHjRUEs9Igt2YUlSY01bIE8GbgROAR/rse9h4GXAm4DrgKPA4Yi4dtBC7hVn45Wk5mYb7veJlNLjACLieqC2RZGHxHOBG1JKh/NttwFPAt4BHOm3gBFxoscuy/0esxPHQCSpuUYtkJRSu+HxXgKsAh8qfTYB7we+KyKu6ruEe8jZeCWpuaYtkKaeDhytCZy7yu8DRMQcMEMWYnMRsQCsVT+bUlrpdsK8hTKUVkjRAllbb7PRTsy0YhiHlaRHpWFfhXUIOF6z/Xjp/cJ7gLPAvwZ+N3/9b4Zcnr4UYyBgK0SSehnFZbypyXsppVenlKLyuH0E5WmsaIGAl/JKUi/D7sI6xvZWRuFg/lzXOmksIm6v2XxgN8csWygFiJfySlJ3w26B3A08NSKqx706f/7CkM83VIt2YUlSY8NugRwGXgO8kNKVWMArgXtSSkd3c/CU0jXVbcMcRN82BmIXliR11ThAIuKl+ctn5M/PiYjLgNMppY/k244AtwHvjYhDwH3Aq4BnAy8eTpFHZ2F2K0C8mVCSuuunBfKBys9vyZ/vB66E7J6P/EbDt+ePFbLLdm9IKd2yq5LugVYr2Dfb4vx62zEQSeqhcYCklBrdFJFSOgm8Pn8M1agH0SEbBzm/3nYMRJJ6cEGpis3pTOzCkqSuhj2IPlKjHkSH0qqEtkAkqStbIBWbqxLaApGkrgyQCmfklaRmDJCKRWfklaRGJmoMZE+uwnIQXZIasQVS4broktTMRLVA9uQqLMdAJKkRWyAVCwaIJDVigFQUXVhOZSJJ3RkgFUt5C8TJFCWpu4kaA9mrubDAQXRJ6sUWSEUxBmIXliR1N1EtEK/CkqSLhy2QiiXvRJekRgyQioV5B9ElqQkDpKLowlpbb7PRTmMujSRdvAyQiiJAwIF0SerGAKkoLuMFx0EkqZuJugprL2fjBe8FkaRubIFU2AKRpGYmqgWyl/eBgC0QSerGFkjFwpwtEElqwgCpmGkF+2azajFAJKkzA6SGEypKUm8GSA3XRZek3gyQGk6oKEm9GSA1nNJdknozQGosOaGiJPU0UfeB7MWd6FAaRLcFIkkd2QKpseAguiT1NFEtkL24Ex22BtEdA5GkzmyB1HBVQknqzQCpUXRhOYguSZ0ZIDWKQXS7sCSpMwOkhneiS1JvBkgN70SXpN4MkBpOpihJvRkgNWyBSFJvBkgN70SXpN4MkBoOoktSbxN1J/pez4V1fr3NRjsx04phn0KSJp4tkBqLpXXRvRdEkupNVAtkz+bCmt8KkLMXNti/b6KqSZL2hC2QGuUWiOMgklTPAKlRboHYhSVJ9QyQGuUWiBMqSlI9A6TGwtz2MRBJ0k4GSI2ZVjA/m1WNASJJ9QyQDjZXJbQLS5JqGSAdFKsSOgYiSfUMkA6cUFGSujNAOigG0r2MV5LqGSAduCaIJHVngHSw5JTuktSVAdJB0YXlILok1TNAOlh0DESSujJAOvAqLEnqzgDpwEF0Sepuoha62KsVCcF10SWpF1sgHbguuiR1N1EtkL1akRAcA5GkXmyBdLBgF5YkdWWAdLDkbLyS1JUB0kExiH7GFogk1TJAOnAQXZK6M0A6KKYyOb/ept1OYy6NJF18DJAOii4sgHPrtkIkqcoA6WCpFCBOqChJOxkgHRRjIOA4iCTVMUA6WCgFiDPyStJOBkgH5TEQbyaUpJ0MkA7KXVif/MoxNrwSS5K2MUA6mGkFly8vAPBLH/m/PP/XP8GRzz/oJb2SlDNAurj5tc/kuqu/A4AvffMUP3nzZ7juN+7g1qMPkZJBImm6xaT/RxgRJ5aXl5dPnDgxsnPc/U+r/NqtX+IvvvjQ5rbvvmKZn33eU3jOUx5LRIzs3JI0CisrK6yurq6mlFYGPYYB0ofPff0Ev3rrvfzlvd/a3HZw/zxXHlriysv2c+Wh/Tzx0BL/4rL9PPHQfpYX50ZeJkkahAHC3gZI4c6vHedXb72XT37lWNf9VpbmuOzAPg7un+eyA/Mc3D/Pof37OJS/Xl6cY2l+hsW5WZbmZ7LHvlkW52aYadmqkTQ6BgjjCZDC0X86yT0PneS+h89w/7HTfO3h09z38GlOnlvf9bH3zbaYn2kxOxPMzbTyRzA702K2lW1rtYJWwEwErVbkz9CKyB8Q+TMEEWTb8teRvyZ7m4DN7rjs9c5tlD4TWx+t/Bxbr2Nrj+J41X2rx6p2CUblnOVjVMu24/ydylU6ad0+27ZFtYw1+2z7zvV1WC47EaX3tx8v8j+M2j+P2HmuuuOU6z77c48dxyifo/YYO8q29XeI6nbyc5T+3Kt/z7Z9x81zbp27Varn4u9tp2MU59r8HjXby59rxfZz2e08nACZqBUJLzZXXX4JV11+yY7tJ86scd/Dp3ngkbMcO3We46fXePj0GsdPreWvs20nz16g00Vd59fbnF9vj/gbSNNrZ7BuD7MikIhyMJWDsgimSkh12lb6bBGW5Z/LZSj/8lcN6LrP3/SK72N+du+viTJARmBlaZ7vecI83/OES7vul1Li/Hqbs2sbnLmwwdm1dc6sbXBmbYOzaxusbbRZ30hc2GhzYaPNert4nVjfaNNO0E6JjXainRLtdqKdYCNlP5O/306Q8teQPWdvp7wckPLnvGTZti77FNsrT6SUSq/rt1M5xuYxU822mnJtbU879ilv336c1PE85XIV+1XLXT4+O7aV62vrc+Vz1dZF2ipvqu5b+sWi9v3N85d/rtR95TzlP/fN8uw45vbv9miW8n8rpS1jK8tupTGV3QAZo4hgYW6GhbkZukeNNB4p7QyXIoiy97cH0OZ7XYIpUfxyU7NP2n7MdvlzlWO02zt/MdjaP20/JtBu1x+rXfpebDt3tl+77txp65eEdqnc7cp526X9i188inJsHnfbMUu/6OX7bbSLXwTStuOV62lmTF1yBoikjmKzuwZKIzgSMMYbCSPi2yLizyPiTETcFRHfO66ySJL6N8470X8L+DJwCHg38MGI8MYJSZoQYwmQiHgM8MPAW1NKZ1NKN+VvPWsc5ZEk9a9RgETEFRHxzoi4IyJORUSKiGs67HsgIt4VEQ9GxNmIuDMiXlTZ7V8CJ1JKD5W2fR64aqBvIUnac01bIE8GbgROAR/rse9h4GXAm4DrgKPA4Yi4trTPfuBk5XMngQMNyyNJGrOmV2F9IqX0OICIuB6otijI37sWeC5wQ0rpcL7tNuBJwDuAI/mup4HHVD5+CVlASZImQKMASSk1vSX6JcAq8KHSZ1NEvB+4KSKuSikdBb4EXBoR31bqxno68D+rB4yIXnOULK+urrKyMvDd+JI0dVZXVyH7xX1gw74P5OnA0ZrAuavy/v+LiA8Db46INwAvJ7vI/JMDnjetrq5Wu8SaWs6fVwf8/LSxvvpjffXH+urPburrEmBX8yUNO0AOAffWbD9eer/wE8DvA4+QXc77IymlC9UP7mairyaKFs6oz/NoYX31x/rqj/XVn3HX1yjuRO82KcvWtENZ19XzRnB+SdIeGPZ9IMfY3sooHMyfj9e8J0maQMMOkLuBp0ZE9bhX589fGPL5JEljMuwAOQysAC+sbH8lcE9+BZYk6VGg8RhIRLw0f/mM/Pk5EXEZcDql9JF82xHgNuC9EXEIuA94FfBs4MXDKbIk6WLQeEnbiOi04/0ppStL+10CvB14KVlr5CjwtpTSn+2uqKMx7qsYJo311R/rqz/WV3/GXV8Tvyb6bo37D2DSWF/9sb76Y331Z9z1NfUBIkkazDjXA5EkTTADRJI0EANEkjSQqQyQhoteTaU+Fw/7DxHxuYg4FxEPRMQvRcTCHhd5bCLihyLifRFxT0ScyevgTyPi6pp9nxcRn8r/vn0zIn47IqZuoDginhURH42Ib+R/b74VER+PiBfU7GudVUTEW/J/k5+teW/P62sqA4Rmi15Nq0aLh0XEy4Gbgb8GXkB26fZPAe8bfREvGq8DngD8Glkd/Kf857+LiGcWO+UBfAT4OtlNtm8gW1Pnf9fM2vBodylwD/BzwPOBHwfOA0ci4keLnayznSLiacDPAw/VvHcN46ivlNJUPYBrySZ1fElpWwB3AF8cd/nG/QBapdfX53V1TWWfGeBB4EOV7a/N9//+cX+PPaqrx9VsWyGbYfqDpW2fBv6hUrfPy+vq34/7e4z7QXZD89eBj1tnHeuoBXwK+A3gduCzlffHUl/TmOS1i14B7we+KyKmel321GzxsGcC305WZ2U3AxeAHxl2uS5GKaVv1mw7QbZg2hUAEfF4stkbfq9ctymlW4FvMCV11U1KaZ3s3+QFsM46+Fmyv1P/ufrGOOtrGgOkyaJX6q6oo22TY6aUzgBfYYrrMCIeS/b9i7qpravc55nSuoqIVkTMRsTlEfFW4ClkXYFgnW0TEU8C3ga8PqVUt3De2OprGgPkEPXTytcteqV6RR11qseprMOICOAmsn9XxfLM1lW9PyZrcXwD+I/Av0sp/Z/8Pessl/+deg/w0dR5Oqix1dc0Bgg0XPRKPXWqq2mtw18hGzd6XUrpi5X3rKvt3gj8K7KB3iPAH0fEjZV9rLNsXPH7gJ9usO+e19coViS82Lno1e4dy58PlV4XDpLNwjxVIuK/kV1Z9DMppfeV3irXVdVBpvTvW0rpq8BX8x9viYhbgN+MiD/COgMgn+38l4H/DpwuXZI7C8zkP59jjPU1jS0QF73avbvz5219qxGxBHwnU1aHEfE24BeBN6aU3lV5u7auclczZXXVxafJLvF9LNZZ4QpgmSxAHik9foCsbh4B3sIY62saA8RFr3bvU8A/A6+obL8RmAP+dM9LNCYR8V+ANwNvTin9SvX9lNIDwJ3Ay8q/tETEDwGPZ4rqqpO8n/8a4ARwzDrb9GXgB2senyO7WOUHgZvGWV/T2IXlolc99Fo8LKW0HhG/ALwvIt4N/AnwVOB/AH+SUvrU3pd670XEz5H9Bvhh4C/KNw8C51NK/5C//nngz4E/iIibgMvJ6upvgQ/sXYnHLyJuBu4H/h54GPgOsn9//xb46fySXrDOSCmdIrvnY5vSFO7l98ZTX+O+QWZMN+VcAryb7Lfoc8BngOvHXa6L5UE26Fb3+Fplv5eTXSZ4nuxqml8GFsdd/j2sp9v7qKvn5/+YzwHfIruy5tJxf4cx1Nnrgb8h67dfz58/CrywZl/rrL4Ob6dyI+G46sv1QCRJA5nGMRBJ0hAYIJKkgRggkqSBGCCSpIEYIJKkgRggkqSBGCCSpIEYIJKkgRggkqSB/H9sMhPpxo6LuQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy([np.linalg.norm(A @ x - (x @ A @ x) * x) for x in conv_m])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "32.29176750885745\n"
     ]
    }
   ],
   "source": [
    "shift = 30\n",
    "large_eigv, conv_l = shift_inverse_iteration(A, x0, num_iter, shift)\n",
    "print(large_eigv @ A @ large_eigv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1a273a4940>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAEMCAYAAABtKgnyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmYk+W5x/HvncwGDGQERGVRqqiAIApqXSu2WgX3ilbErR6P1qq19RSr1p5ztB5b1x6pPVYtFfddtFrcqliLawFlG0VUQEEUBWYQGIaZyX3+SAZCyAyzZPJm+X2u670yvu+T5M4jM788yZ035u6IiIgEIRR0ASIiUrgUQiIiEhiFkIiIBEYhJCIigVEIiYhIYBRCIiISGIWQiIgERiEkIiKBUQiJiEhgFEIiIhIYhZCIiARGISQiIoFRCImISGAUQiIiEpiioAvIBmZWTyyQVwddi4hIDukGRN29zVli+j4hMLMoYJFIJOhSRERyRnV1NYC7e5tfVdNKKGZ1JBKJVFVVBV2HiEjOqKiooLq6ul2vIOk9IRERCYxCSEREApOzIWRmpWZ2t5ktMbNqM5tqZnsEXZeIiLRczoYQsfezPgH2B7oDzwBPBVqRiIi0St50x5lZCbAe2NbdV7TyulVqTBARaZ14Y0K1u1e09TYythIys75mdquZTTOzNWbmZjayibHlZjbBzJaZWY2ZTTez47ZyFwcCy1sbQCIiEpxMvhw3ABgLrAFe3srYycA44CrgaKASmGxmo1MNNrMK4A7gyrRV20LrNtRz7bOVfLpiXabvWkQk52Xyc0KvuXsvADM7AUi5sokHzeHAD9x9cnzfVGBn4GZgStL4MmLvBT3r7n/puPJT+9Hd/+LthSuZs7Sah/59f0Ihy3QJIiI5K2MrIXePtnDoiUA18HTCdR24BxhoZoMb95tZGHgYWAr8oqkbNLOq5jagzadK+PHIXQB4e+FK7n1zUVtvRkSkIGVjd9wQoDJFaM1OON7oLqATcLYH1GFx2O69OGWfvgBc//x8Fn29NogyRERyUjaGUA9gZYr9KxOOY2Y7AT8CvgOsijc7rDGzQ5Kv6O4VzW3EVl5tdtUxg9khUkZNXQOXPT6baDQ/Og5FRDpaNoYQQHN/xR3A3Re7u7l7J3cvT9j+maEaN+pWVszvTtoTgHcWrWTSG4syXYKISE7KxhBaQXy1k6R7/DLVKilwh+62Lafu2w+AG174gIV6WU5EZKuyMYTmAYPMLLm2ofHLuRmup8V+dfQgekfKWF8XZfxjs2jQy3IiIs3KxhCaDFQAxybtPxOY7+6VmS+pZbqWFXP9mNjLctMXr+Lu1xcGXJGISHbLaAiZ2RgzGwMcEN91aHzfqIRhU4CpwEQzO8fMDjOzScDBwPhM1tsWh+y6LWP32xGAG1+Yz8dfrQm4IhGR7JXRc8eZWVN3ttjd+yeM6wZcB4whtiqqBK5x9w45QWm6zx23praeI3//Gkurahi+YwWP/fhAwvoQq4jkmZw6dxxAvJst1dY/adxqd7/I3bd39zJ3H95RAdQRykuLuCH+stzMT6uYOO2TgCsSEclO2fieUF44aEBPTt8/9rLcTS9+yNyl7fookohIXlIIdaArRg2iX/dObKiPcs6kf7G0qibokkREsopCqAN1KS1i4ln70rWsiOXf1HL2X96huqYu6LJERLKGQqiD7bZdV+44YwTFYWPB8jWcf990ausbgi5LRCQrKIQy4MBdenLjmGEAvPXJSn75+Gzy5RttRUTaQyGUISfs3YfxR+4OwFPvfc5NL84PuCIRkeAphDLoJyN32fhB1j9O/ZgH3/404IpERIKlEMogM+M3x+/BYbtvC8Cvn57L1A+WB1yViEhwFEIZVhQOcdtpwxnaJ0JD1LnwwZnMWaLPEIlIYVIIBaBLaRETz96HPhWdWLehgfPum87q9WrdFpHCoxAKSK+uZdxzzr50Kg6zrHo91z6btScHFxHpMAqhAA3o1ZVfHhXrmHt0+hKmztf7QyJSWBRCATvzgP58+1uxL4294ok5OqOCiBQUhVDAQiHjxjHD6FwS5ovV6/mNXpYTkQKiEMoCO/bozOWjBgLw+IwlvPLBlwFXJCKSGQqhLHH6t3di/53jL8s9OYfqdXpZTkTyn0IoSyS+LPfl6lqufnZe0CWJiHQ4hVAW6de9M1eMHgTAkzOX8vdKvSwnIvlNIZRlxu23Iwfu0gOAKyfPoWrdhoArEhHpOAqhLBMKGdeftCddSsIs/6aWq59Rt5yI5K+cDiEzu9rMKs0samanBl1PuvTr3pkrj469LDf53aU8P/eLgCsSEekYOR1CwALgEuCdoAtJt9P225GDB/QE4PInZ7OsuibgikRE0i+nQ8jd73f3l4D1QdeSbmbGTScPY5vOxVStq+OSh9+jIapvYxWR/JKWEDKzvmZ2q5lNM7M1ZuZmNrKJseVmNsHMlplZjZlNN7Pj0lFHvtk+UsZNJ8e+FvydhSv5wysLAq5IRCS90rUSGgCMBdYAL29l7GRgHHAVcDRQCUw2s9FpqiWvfG/QdvzooP4ATHh5AW9/siLYgkRE0ihdIfSau/dy96OAu5saFA+aw4Fz3X2iu78CnAW8CdycNPa5+Koq1XZZmurOCZePGsgevbsRdfjZI++xaq3atkUkP6QlhNw92sKhJwLVwNMJ13XgHmCgmQ1O2D/K3cub2G5oTX1mVtXcBkRac3uZVloU5g9j96ZzSey7hy57YjaxaRMRyW2ZbkwYAlSmCK3ZCcdbzMyKzayM2OMoNrMyM8vpZoum7LxtOdccH5uelyq/5L63FgdckYhI+2X6D3YPYGWK/SsTjrfGXUANcAhwb/zn7yQPcveK5jZiq7Osd9LwPhy/V28Arv3b+1R+vjrgikRE2ieIVUNzryO16jUmdz/b3S1pe7V95WUvM+PaE4awU4/ObKiPcvFDM1m3oT7oskRE2izTIbSC1Kud7vHLVKskSdC1rJgJp+5NUcj4+Ku1XP1XndZHRHJXpkNoHjAoxfs2Q+OXczNcT04a1q+Cy47aHYBHpn/GEzOWBFyRiEjbZDqEJgMVwLFJ+88E5ru7nta30LkH78x3B/YC4FdPzeH9ZXp/SERyT9pCyMzGmNkY4ID4rkPj+0YlDJsCTAUmmtk5ZnaYmU0CDgbGp6uWQhAKGbecMoy+23RifV2UC+6fwer1+jZWEcktlq7Pm5hZUze02N37J4zrBlwHjCG2KqoErnH3p9JSSBuYWVUkEolUVVUFVUKbzVlSzUl/eoMN9VGO3GM7/nT6CMws6LJEpABUVFRQXV1dHe8ybpO0rYRSdKk1bv2Txq1294vcfXt3L3P34UEGUK4b2jfCNcftAcAL877krn9+EnBFIiItl5cf7Cw0P9y3HyeP6AvA9c/P5y2dX05EcoRCKA+YGb85YQiDduhGQ9S56MF3Wb46777dQkTykEIoT5QVh/nT6cPpWlbE12tqufDBmdQ1tPSUfiIiwVAI5ZGdenThllP2AuBfi1Zxw/MfBFyRiEjzFEJ55ojB2/GTkbsAcNc/FzJlzrKAKxIRaZpCKA9desRuHLhL7OxI4x+bxUfLvwm4IhGR1BRCeagoHGLC2L3ZIVLG2g0NnHffDL7RB1lFJAsphPJUz/JSbj99BCXhEJ98tZZLH51FNKovwhOR7KIQymN79avgmuNjH2R9qfJLbv/HxwFXJCKyOYVQnjt1vx0Zu9+OANz04nxenb884IpERDZRCBWA/z5uMMP6VeAOlzz8Hp+uWBd0SSIigEKoIJQWxT7I2rO8hOqaOs6/fwY1GxqCLktERCFUKHaIdOK204YTDhnvL1vNFU/OJl1nUBcRaSuFUAHZf+ceXDl6EABPvfc5k95YFGxBIlLwFEIF5pyD+nP8Xr0BuPZv7/O2zrgtIgFSCBUYM+N3P9iTgdt3pSHqXPjgTJZV1wRdlogUKIVQAepUEuaOM0YQ6VTM12s28OP7Z7K+To0KIpJ5CqECtVOPLkwYuzdmMOuzKv7z6blqVBCRjFMIFbBDd9uW8UfuDsCj05fwwNufBlyRiBQahVCBu+DQXRg1ZHsArn5mHtMXrQy4IhEpJAqhAmdm3HjyMHbtVU5dg3PBAzP5Ul8NLiIZkvMhZGYHmFnUzC4PupZcVV5axJ1n7kPXsiK++qaWC+6fQW29GhVEpOPldAiZWQj4X+CdoGvJdd/q2YVbT90LM5j5aRVXP1MZdEkiUgByOoSA84HXgQ+CLiQffHfgdvz88N0AePDtT3noHTUqiEjHancImVlfM7vVzKaZ2RozczMb2cTYcjObYGbLzKzGzKab2XFtvN8ewCXA1e0oX5JcdNgADh+0HQD/9fQ8Zn66KuCKRCSfpWMlNAAYC6wBXt7K2MnAOOAq4GigEphsZqPbcL+/BW5x9+o2XFeaEAoZv//hMHbetgsbGqJccP8Mln+jRgUR6RjpCKHX3L2Xux8F3N3UoHjQHA6c6+4T3f0V4CzgTeDmpLHPxVdVqbbLzGwEMBz4c0sKNLOq5jYg0sbHnpe6lhVz5xn7UF5axJera7nwgZlsqI8GXZaI5KF2h5C7t/Sv04lANfB0wnUduAcYaGaDE/aPcvfyJrYbgEOAwcByM/saOBX4tZn9pb2PR2IG9CrnllOGAfCvRau49m9qVBCR9MtkY8IQoDJFaM1OON5SfwZ2A/aKb38FbgV+kWqwu1c0txELR0ny/T2256ff2xWAe99czKPTPwu4IhHJN5kMoR5Aqo/jr0w43iLuvsbdlzRuwDpgtbvr4/5p9rPv7cr3BvYC4Kqn5jLrs6qAKxKRfJLpFu3mzpDZ5rNnuvvZ7v67tl5fmhYKGb8/dS927tmFDfVRzr9vBl99Uxt0WSKSJzIZQitIvdrpHr/UKiZLdSsr5o4zRtClJMwXq9dz4YMzqWtQo4KItF8mQ2geMCh+loNEQ+OXczNYi7TSrtt15eZT9gLgnYUr+Z+/vR9wRSKSDzIZQpOBCuDYpP1nAvPdXe1XWe6oIdtz8XcHADDpjUU8MWNJwBWJSK4rSseNmNmY+I/7xi8PNbOewFp3fy6+bwowFZgYP9vBQmKfEzoYOD4ddUjH+9nhuzF3aTVT53/FFZPnsNt2XRnaVx+zEpG2sXR8m6aZNXUji929f8K4bsB1wBhiq6JK4Bp3f6rdRbSDmVVFIpFIVZU6v1qiuqaO42+bxqIV6+gdKeOvFx9Mz/LSoMsSkQyrqKigurq6Ov5RlzZJSwjlOoVQ63345Tec8MfXWbehgW9/qzv3n/ttisO5fj5cEWmNdISQ/mpIm+y2XVduPjl2RoW3F67kuilqVBCR1lMISZuNGroDFx62CwB3v76IJ2eqUUFEWkchJO1y6RG7M3L3bQG44sk5zF2qMyCJSMsphKRdwiHj1h/uzU49OlMbP6PCijU6o4KItIxCSNot0jn21Q+dS8IsrarhogffpV5nVBCRFlAISVrsvn1XbhwTa1R485MVXDdF37guIlunEJK0OXrPHbhgZKxR4S+vL2Tyu2pUEJHmKYQkrX7x/d35zm6xRoXLn1Cjgog0TyEkaRUOGRNO3Ysdu29qVFi5dkPQZYlIllIISdpVdC7hzjNH0Km4sVFhphoVRCQlhZB0iIHbd+PGk/cE4I2PV/C759SoICJbUghJhzlmz978+NBYo8Kfpy3k6feWBlyRiGQbhZB0qPFH7s4hu/YE4JdPzFajgohsRiEkHSocMv4wdm927N6Z9XVqVBCRzSmEpMNVdC7hjjM2NSpc/JAaFUQkRiEkGTFoh02NCq9/tILrn1ejgogohCSDjtmzN+cfujMAd/1TjQoiohCSDLvsyIGbNSrM+1yNCiKFTCEkGdXYqNCve6eNjQqr1KggUrAUQpJxFZ1LuPOMfehUHGbJqhouUqOCSMHK6RAyszVJm5vZSUHXJVs3aIduXD9mU6PCDS/MD7giEQlCToeQu5c3bsCBQA3wYsBlSQsdN6w3538n1qhw52ufqFFBpADldAglGQc87e7fBF2ItNz4I3fn4AGbGhUqP18dcEUikklpCSEz62tmt5rZtISXxUY2MbbczCaY2TIzqzGz6WZ2XDvv34CxwAPtuR3JvKJwaLNGhfPum65GBZECkq6V0ABiIbAGeHkrYycTW7VcBRwNVAKTzWx0O+7/O0Bn4IV23IYEZJsuJdxx+j6UFYdYsqqGix96V40KIgUiXSH0mrv3cvejgLubGhQPmsOBc919oru/ApwFvAncnDT2uRSNB43bZUk3PQ541N3rmrjfquY2INKOxy5pMLh3N24YMwyAaR99rUYFkQKRlhBy95Y+bT0RqAaeTriuA/cAA81scML+UYmNB0nbDY3jzKwEGINeist5xw3rzXkJjQp/nfV5wBWJSEfLdGPCEKAyRWjNTjjeWqOJBdsbTQ1w94rmtvj1JQtcduTuHDSgR+znx2epUUEkz2U6hHoAK1PsX5lwvLXGAQ/GV1SS44rCIW4bO5y+28TPqHC/GhVE8lkQLdrNhUWrg8TdT3b3X7WjHsky23SJffVDWXGIz1bW8NOH1aggkq8yHUIrSL3a6R6/TLVKkgK0R+8I158UO6PCPxd8zY0vqlFBJB9lOoTmAYPMLPl+h8Yv52a4Hslix+/Vh38/5FsA3PGPT3hGjQoieSfTITQZqACOTdp/JjDf3SszXI9kuV8eNTChUUFnVBDJN2kLITMbY2ZjgAPiuw6N7xuVMGwKMBWYaGbnmNlhZjYJOBgYn65aJH8kNirU1DVw/v3TqVqnRgWRfGHpaiozs6ZuaLG7908Y1w24jthneyqInTHhGnd/Ki2FtIGZVUUikUhVVVVQJchWzPu8mpNuf4P1dVEO2bUnk360H+GQBV2WSEGrqKigurq6Ov5RlzZJ20rI3a2JrX/SuNXufpG7b+/uZe4+PMgAktyQ3KhwwwsfBFyRiKRDPp1FW/JccqPCs7PVqCCS6xRCklMSGxXGPzab95epUUEklymEJKfEvvphOH0qYo0K592nRgWRXKYQkpzTvUsJd5656YwKFz/0Lg1RnbVJJBcphCQnbXFGBX31g0hOUghJzjp+rz6ce3CsUeFP//hYjQoiOUghJDnt8lEDOXAXNSqI5CqFkOS0onCI205To4JIrlIISc7rHv/qh9IiNSqI5BqFkOSFIX3UqCCSixRCkjdO2LsP/6ZGBZGcohCSvHJFUqPCB1+oUUEkmymEJK/Ezqiw96ZGhXtnqFFBJIsphCTv9Cgv3dio8OnKdWpUEMliCiHJS2pUEMkNCiHJW2pUEMl+CiHJa1eMGsgBO6tRQSRbKYQkr8XOqKBGBZFspRCSvJfcqPDTh99To4JIllAISUEY0ifC704aCsBrH37FTS+qUUEkG2R1CJnZ1WZWaWZRMzs1xfHtzOxFM1tnZrPNbEQQdUpuOHHvvpxzUKxR4fZXP+Zvs5cFXJGIZHUIAQuAS4B3mjh+O/AR0AO4DXjCzIozVJvkoCtHb2pU+MVjs9SoIBKwrA4hd7/f3V8C1icfM7OuwDHA1e5e4+53xg8dmMkaJbeoUUEku7QohMysr5ndambTzGyNmbmZjWxibLmZTTCzZWZWY2bTzey4tFYdsytQ5e5fJuybAwzugPuSPKJGBZHs0dKV0ABgLLAGeHkrYycD44CrgKOBSmCymY1ua5FN6AIkv5ayGihP8/1IHlKjgkh2aGkIvebuvdz9KODupgbFg+Zw4Fx3n+jurwBnAW8CNyeNfS6+qkq1XdaCmtYCXZP2dSMWlCJbpUYFkeC1KITcPdrC2zsRqAaeTriuA/cAA81scML+Ue5e3sR2QwvuawGwjZltl7BvCLGV12bMrKq5DYi08PFJnrly9ED237k7AOMfn8X8L74JuCKRwpLuxoQhQGWK0JqdcLzFzKzYzMqI1VlsZmVmFgJw92+AZ4Ffx/efCxjwRrsegRSUonCIP542nD4VnVi3oYHz7ptO9bq6oMsSKRjpDqEewMoU+1cmHG+Nu4Aa4BDg3vjP30k4fgGwO7CKWCv3Se6+xV8Qd69obiO2epMCldiosHjFOn76sL76QSRTOqJFu7nf3lb9Zrv72e5uSdurCce/dPcj3L2Tuw919xltLVoK25A+EX77g1ijwj8+/Iqb1aggkhHpDqEVpF7tdI9fploliWSFHwzvy48O6g/A/736MVPmqFFBpKOlO4TmAYMa37dJMDR+OTfN9yeSVleOHrSxUeEXj6lRQaSjpTuEJgMVwLFJ+88E5rv7Fp1rItmkWI0KIhnV4hAyszFmNgY4IL7r0Pi+UQnDpgBTgYlmdo6ZHWZmk4CDgfHpKlqkI6lRQSRzLPYxnhYMNGtq4GJ3758wrhtwHTCG2KqoErjG3Z9qX6kdx8yqIpFIpKqqKuhSJIs8OXMJlz46C4ALD9uF8UcODLgikexSUVFBdXV1dbzLuE1avBJK0aXWuPVPGrfa3S9y9+3dvczdh2dzAIk0JbFR4Y9TP+Y5NSqIpF1Wn0VbJGiJjQr/oUYFkbRTCIk0Q40KIh1LISSyFT3KS/nT6ZsaFS55RI0KIumiEBJpgaF9N51R4dX5X3HLSzqjgkg6KIREWkiNCiLppxASaQU1Koikl0JIpBWKwyFuO204vSNlalQQSQOFkEgr9Swv5Y4z9lGjgkgaKIRE2kCNCiLpoRASaSM1Koi0n0JIpB2SGxU+/FKNCiKtoRASaYfGMypsbFS4V40KIq2hEBJppx7xRoWSohCL1Kgg0ioKIZE0GNo3wm9P3NSo8PuXPgy4IpHcoBASSZOTRvTl7AP7A3Db1I94fq4aFUS2RiEkkka/OnoQ+30r3qjw6CwWqFFBpFkKIZE0Kg6H+L9xw9khUsbaDQ2cd98MqmvUqCDSFIWQSJr1jH/1Q0lRiIVfr+Xnj7xHVI0KIikphEQ6wLB+FfzPCUMAeOWD5fzv39WoIJKKQkikg5y8Tz/OOmAnACa88hEvzPsi4IpEsk9Wh5CZXW1mlWYWNbNTUxy/08yWmdlqM5tjZscEUadIU646ZjD79Y81Klz6yHt8tFyNCiKJsjqEgAXAJcA7TRy/Bejv7t2Ac4D7zawiU8WJbE1xOMQfxw1n+27xRoV7Z7B6vRoVRBpldQi5+/3u/hKwvonjH7h7bfw/64FSoE+m6hNpiW27lvKnM0ZQEg7xyddr+fnDalQQadSiEDKzvmZ2q5lNM7M1ZuZmNrKJseVmNiH+MlmNmU03s+PSWvXm9/d/ZlYDzAReBio76r5E2mqvfhVcG29UePmD5dz68oKAKxLJDi1dCQ0AxgJriP2hb85kYBxwFXA0sVCYbGaj21pkc9z9J0A5cATworvrKaZkpVP27ccZ+8caFW59eQEvqlFBpMUh9Jq793L3o4C7mxoUD5rDgXPdfaK7vwKcBbwJ3Jw09rn4qirVdllrHoS7N7j734Hvm9mRrbmuSCb9+pjB7Nt/GwAufXQWHy1fE3BFIsFqUQi5e7SFt3ciUA08nXBdB+4BBprZ4IT9o9y9vInthlY8hkRhYJfknWZW1dwGRNp4fyKtUlK0qVFhTW095903XY0KUtDS3ZgwBKhMEVqzE463mJkVm1kZsTqLzazMzELxY+Vmdnr8ssjMTgYOA6a18zGIdKheXcu4/fThsUaFr9Zy6SOz1KggBSvdIdQDWJli/8qE461xF1ADHALcG//5O/FjTqwtewmwArgcGOvus5NvxN0rmtuIrd5EMmbvHbfhNyfsAcDf3/+SCa+oUUEKU0e0aDf3lK5VT/fc/Wx3t6Tt1fixte7+3XiQRNx9hLtPbk/hIpn0w3135PT9dwTgf/++QGdUkIKU7hBaQerVTvf4ZapVkkjB+s9j9tjYqPDzR97jgy9WB1yRSGalO4TmAYMa37dJMDR+OTfN9yeS00qKQtx++gh6R8pYt6GBc++Zzoo1tVu/okieSHcITQYqgGOT9p8JzHd3fZBUJEnP8lLuOmsfOhWHWbKqhgsemMmG+pY2pIrkthaHkJmNMbMxwAHxXYfG941KGDYFmApMNLNzzOwwM5sEHAyMT1fRIvlmj94RbjllGADvLFzJf/11HvrctRQCa+k/dDNrauBid++fMK4bcB0whtiqqBK4xt2fal+pHcfMqiKRSKSqqiroUqTA/f6lDzee0uea4/fgzAP6B1uQSDMqKiqorq6ujncZt0mLQyifKYQkW0SjzoUPzuS5uV8QDhn3nbMfBw7oGXRZIimlI4Sy+izaIoUmFDJuPmUYg3boRkPU+cmDM1m8Ym3QZYl0GIWQSJbpXFLEXWeOoEeXEqrW1fFv90znG53aR/KUQkgkC/XdpjN3nDGC4rDx0fI1jLn9TZ6cuYTa+oagSxNJK70nhN4Tkuz16PTPuOzxTWei6lleyhn778S4/XekZ3lpgJWJqDEhbRRCks3e/XQVf562kOfnfkFD/ESnJeEQx+3Vmx8d1J89em/9JPD1DVGqaupYuXbDxm3F2g2srqkjHDJKi0KUFoVjl8UhyorClBaHCJmxoT5KbX0D6+til7X1UWrjP4dCtvF6ZcXx6xeFKC0OUxwyahs2ja2tj8avG/sZ2Dh24/Xi91saDlEX9Y1ja+M11NZFWV/fQDTqm8Ym1F5WHKakKETUPX6/W9ZeVx+lJKHOzWsPY8Zm95V4Oxvqo4RDoc2vU7zp/otCtqnW+i0fO7DFPDX+XBIOsaEhusU8rU+er8R5aryNcIj6qG8xT421N0QbH/OWNTfWM6xfBcXh1r04phBKE4WQ5ILPq2q4983FPPTOp1TXbHqPaM++EcpLi6hvcDY0RKlriFLf4NQ1RKmLRvlmfT3VNXXoV12aM+s/v0+kc3GrrqMQShOFkOSSmg0NPPnuEu5+fVGbvxSvW1kRPcpL6dapGN+4athytdLgvsVKo/GZdEk4RNTZbHxtwrPvDQ2xZ99lTax2gCbvt7Y+Skk4tMWz9pL4bYWN2KqhidVO2GyLlUZj7UXhEHUJK7TE662va8Bh40oweZVUHDYanE11Jq3U6hq8yVVSSVEIb5yvpBVWbd2m+UpelZYWhSmLz9eG5FVhwv0XhzetSjd77PH5qmvwlCva9fUNuMMHvzmKsuJwq/4dKYTSRCH7jH5UAAAF30lEQVQkucjd+eeCr3nrkxWEzCgOhygKGyXhEMVho7goRHEoROfSMN27lNCjSynbdClmm84lLX7Zxd0xszbXl2vXDfK+g7xufdQpClmrb0MhlCYKIRGR1tOHVUVEJKcphEREJDAKIRERCYxCSEREAqMQEhGRwCiEREQkMGrRBswsClgksvXTn4iISEx1dTWAu3ubFzQKIcDM6omtCle38SYa06s6PRXlPc1X62i+Wkfz1Trtma9uQNTdi9p65wqhNDCzKoD2fGCrkGi+Wkfz1Tqar9YJer70npCIiARGISQiIoFRCImISGAUQiIiEhiFkIiIBEYhJCIigVEIiYhIYPQ5IRERCYxWQiIiEhiFkIiIBEYhJCIigVEItZGZlZvZBDNbZmY1ZjbdzI4Luq5sYGZ9zexWM5tmZmvMzM1sZBNjTzOzWWa23syWmNnvzKwswyUHxsy+Z2aTzGy+ma2Lz8GTZjY0xdgjzOyt+L+35WZ2h5kV3PnRzOxAM3vBzJbG/918ZWavmNmoFGM1Z0nM7L/jv5PvpTiW8flSCLXdZGAccBVwNFAJTDaz0YFWlR0GAGOBNcDLTQ0ys9OBB4DXgVHAdcCFwKSOLzFr/BjYEfg9sTm4NP7f/zKz/RsHxUN8CvAZcCzwC+A44G9mVmi/x9sA84H/AI4CzgNqgSlmdmrjIM3ZlsxsD+CXwJcpjo0kiPlyd22t3IDRgAMnJuwzYBrwftD1Bb0BoYSfT4jP1cikMWFgGfB00v5/j4//dtCPI0Nz1SvFvgpgFfBEwr53gHeT5vaI+Fz9MOjHEfQGFBH74/mK5qzJOQoBbwF/AF4F3ks6Hsh8FeSzgTQ4kdh3bzzduMNj/8fuAQaa2eCgCssG7h5twbD9ge2JzVmiB4A64KR015WN3H15in1VwAKgL4CZ9QH2Be5LnFt3fwlYSoHMVXPcvZ7Y72QdaM6a8HNi/6Z+lXwgyPlSCLXNEKAyxR/b2QnHpXmNczQ3cae7rwM+poDn0My2Jfb4G+cm5VzFzaFA58rMQmZWZGa9zexqYDdiL2uC5mwzZrYzcA1wkbun+vLOwOZLIdQ2PYCVKfavTDguzWuco6bmsSDn0MwMuJPY7+ZN8d2aq9QeJbbyWQr8DDjF3Z+PH9OcxcX/Td0FvODuTzUxLLD5Ugi1XXOnmtBpKFquqbkq1Dm8kdj7aD929/eTjmmuNncZsB+xN8+nAI+a2dikMZqz2Pus+wAXt2Bsxuerzd8LXuBWkPqZQff4ZapnE7K5FfHLHgk/N+oOLMxsOcEzs/8h1vF1ibtPSjiUOFfJulOg/97c/RPgk/h/PmNmzwB/NLNH0JwBYGY9gRuA3wJrE9qti4Bw/L/XE+B8aSXUNvOAQSnaFhs/25HqdVXZ3Lz45WavNZtZZ2AXCmwOzewa4ErgMnefkHQ45VzFDaXA5qoZ7xBr394WzVmjvkCEWAitStgOIjY3q4D/JsD5Ugi1zWRibbTHJu0/E5jv7pWZLynnvAV8AZyRtH8sUAw8mfGKAmJm/wX8Gvi1u9+YfNzdlwDTgXGJT3zM7HtAHwporpoSf99jJFAFrNCcbfQRcFiKbRaxBqDDgDuDnC+9HNc2U4CpwEQz60HspaOzgIOB44MsLFuY2Zj4j/vGLw+NvzSw1t2fc/d6M7scmGRmtwGPA4OA64HH3f2tzFedeWb2H8SeiT4L/D3xA6pArbu/G//5l8CLwENmdifQm9hcvQ08lrmKg2dmDwCLgRnA18AOxH7/vgtcHG/XBs0Z7r6G2GeCNmNmVfHjiceCma+gP0CVqxvQDbiN2LP59cBM4ISg68qWjdgbmam2RUnjTifWAlpLrMvpBqBT0PVncJ5ebcVcHRX/g7Ae+IpYx9M2QT+GAObsIuBNYu9j1McvXwCOTTFWc5Z6Dl8l6cOqQc2Xvk9IREQCo/eEREQkMAohEREJjEJIREQCoxASEZHAKIRERCQwCiEREQmMQkhERAKjEBIRkcAohEREJDD/D7oiKsDuZC9KAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "res = np.array([np.linalg.norm(A @ x - (x @ A @ x) * x) for x in conv_l])\n",
    "plt.semilogy(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Метод Релея\n",
    "\n",
    "- Адаптивный сдвиг на каждой итерации\n",
    "\n",
    "$$ \\sigma_k = \\frac{(Ax_k, x_k)}{(x_k, x_k)} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def rayleigh_quotient_iteration(A, x0, num_iter):\n",
    "    x = x0.copy()\n",
    "    x = x / np.linalg.norm(x)\n",
    "    conv = [x]\n",
    "    \n",
    "    for i in range(num_iter):\n",
    "        lam = x @ A @ x\n",
    "        x = np.linalg.solve(A - lam * np.eye(n), x)\n",
    "        x = x / np.linalg.norm(x)\n",
    "        conv.append(x)\n",
    "    return x, conv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "eigv_rq, conv_rq = rayleigh_quotient_iteration(A, x0**2, num_iter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.279902142247673\n",
      "[ 1.00129089  1.14452777  1.72651117  4.27990214  5.59525073  8.50188201\n",
      " 12.20387858 13.17684228 15.93954198 32.29176751]\n",
      "[-0.66949928  0.47543029  0.88485231 -0.14937946 -1.50683479  0.50107061\n",
      " -1.98957823 -0.20772708 -1.4726692   0.74007251]\n"
     ]
    }
   ],
   "source": [
    "print(eigv_rq @ A @ eigv_rq)\n",
    "print(np.sort(np.linalg.eigvals(A)))\n",
    "print(x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1a27370470>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAEMCAYAAABtKgnyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8FHX6wPHPs+mQRkioQQEpigiK4oEKAoJSpIrSVBS5886znQXPwgnqqT/r6XkWPDxEURE1IB5NBRVPUBE7ih0FkRJIQiB9v78/Zjdslk1I2d2Z3Tzv12teOzvz3ZknQ8iz35lvEWMMSimllB1cdgeglFKq8dIkpJRSyjaahJRSStlGk5BSSinbaBJSSillG01CSimlbKNJSCmllG00CSmllLKNJiGllFK20SSklFLKNpqElFJK2UaTkFJKKdtoElJKKWUbTUJKKaVsE2t3AE4gIuVYCbnA7liUUiqCpAJuY0y9c4nofEIgIm5A0tLS7A5FKaUiRn5+PoAxxtT7rprWhCwFaWlpaXl5eXbHoZRSESM9PZ38/PwG3UHSZ0JKKaVso0lIKaWUbSI2CYlIgoj8R0S2iki+iKwRkWPtjksppVTtRWwSwnqe9QPQB8gAlgKLbY1IKaVUnURN6zgRiQeKgSxjTG4dP5unDROUUqpuPA0T8o0x6fU9Rthax4lINnA9cCJwPNAUGGiMeStA2WTgTuBcIB34ErjNGPNqDac4BdhZ1wSklAps//79FBQUUF5ejtvttjscFUYul4vExEQyMzMRkdCeK6RHr6oTMAkoBN48TNkcYApwCzAC2ATkiMjwQIVFJB14ArgpaNHWVul+WHEj7P0p7KdWKhTcbjdbt27l559/pqCggLKyMrtDUmFWVlbG7t272bZtG6G+WxbOfkLvGGNaAIjIGGBUoEKeRDMYGGeMyfFsWwN0BO4HlvmVT8R6FvSaMeap0IVfjQXnwZZ34ddP4KLXwBUT9hCUCqb8/Hz27dtHZmYmzZs3x+WK5EfHqr727NnDjh072L17N1lZWSE7T9h+u4wxta3PjwXygSU+nzXA08DRItLNu11EYoAXgG3AddUdUETyalqA+g+VcNrV1uvP78G6R+p9GKWcorCwkPj4eDIzMzUBNWIZGRkkJCRQXFwc0vM48TesO7ApQNL6zGe/15NAEnCRsauFRechcNI0a331HfDbF7aEoVSwuN1uYmNjQ/4sQDlfTExMyJ8HOjEJNQf2BNi+x2c/InIkcDHQH9grIoWepZ//B40x6TUtWDWv+jvzDsjoCBWlkHMplJc06HBKKdVYODEJAdRUqzEAxpgtxhgxxiQZY5J9lrVhivGg+KYwdg6IC3Z8AWvuDHsISikViZyYhHLx1Hb8ZHheA9WS7NeuN5x2jbX+v4dgyzp741FKqQjgxCT0JXCMiPjHdpzn1bkPXU6/AVr1AIx1W65kn90RKaX8zJs3DxGpXGJjY8nOzmbatGls3779kPKFhYXccccd9OzZk6ZNm5KSksLJJ5/MI488os3Xg8CJSSgHq4PqSL/tFwKbjTGbwh9SLcXGw7g5EJMAeVtgZfi7LSmlamf+/PmsW7eOlStXcsEFFzB//nwGDRpUJbHs2LGDPn36cPfddzNy5EheffVVXnrpJfr27cs111zDWWedRVFRkY0/ReQL63xCIjLes9rb83q6iGQC+40xyz3blgFrgLki0hz4EZgKnAaMDme89dLiGDjjb7DqZtg4H7oOh67D7I5KKeXnuOOO4/jjjwfgjDPOYOfOnTz11FOsXbuWQYMGAXDhhRfyzTff8M4779CnT5/Kz5511ln079+f8ePHM2PGDP75z3/a8jNEg3DXhBZ5lhme97M87x/zFvA0tR6D1f/nTmA50AOr8+rScAZbb30ug/aeRnqvXgH7d9sbj1LqsE488UQAdu7cCcCGDRtYtWoV06dPr5KAvM455xyGDh3KnDlz2LVrV1hjjSZhrQkZY2rV8cAYUwBc7lkij8sFYx6Fx06F/btg6VUw4VnQfhcqgpVXuNmeH9qOi3XROi2R2JjgfY/+6aefAOjSpQsAr7/+OgCjRgUc3AWA0aNHs2LFCt5++23Gjx9fbTlVPZ3eO1TSj4Bh/weL/wRfvwbvPgj9rrE7KqXqbXt+Mf3uWWN3GJXWzhhIu4wm9f58RUUF5eXlHDhwgNWrV/P4448zefJkevXqBcDPP/8MQIcOHao9Rvv27YGDCUzVnSahUOo5Cb5fDZ8vgjdnQ3JLOGGK3VEppYCTTjqpyvt+/foxb968Oh3DO1CLji5Rf5qEQkkERv8LCnfCj29bz4eaZkKXs+yOTKk6a52WyNoZA+0Oo1LrtMQGfX7BggV06dKF/Px8nn76aZ555hmuvvpq/vWvfwFwxBFHAPDjjz/StWvXgMfYsmULAO3atWtQLI2ZJqFQi02wngfNGwG/fQYvToWpr0K7k+2OTKk6iY1xNej2l9N069atSuu4/Px8HnvsMS666CJ69+7NkCFDuOmmm1iyZAlDhw4NeIzFixcTFxfH6aefHs7Qo4oT+wlFn8RUOP9laNYByotgwbmw82u7o1JK+XjwwQeJjY1l5syZgHW7bvDgwcydO5f169cfUv7ll19m5cqVTJ48mZYtW4Y73KihSShcklvABa9A0ywozoNnx0H+VrujUkp5dOzYkUsvvZSVK1eybp017Nb8+fPp1KkTgwcPZubMmbz55pusWrWKq6++mkmTJtGzZ08eeUSncGkITUJBUOGu5SwSGR2tGlF8ChRsg2fPgQPOHApPqcZo5syZpKSk8Le//Q2A1q1b8/777zNjxgwWL17MyJEjOeuss3jooYeYOnUq69evJzk52eaoI5vYNQ2Pk4hIXlpaWlpeXl6dP/vtjn1c8fzH3DryWPoeFWjc1QB+eNtKQO4yaNcHLsiB+Oi5164im/dh+5FHHmlzJM60c+dO+vbtS0xMDGvXro3qW3GH+11IT08nPz8/3zMlTr1oTagBjDHclPM5X/+2j6te+JjdhbWcR6jj6dYYcwj8sh5y/gD6ZUCpiNCiRQtWrFjB3r17Oeuss6jPl1d1kCahBhAR7h3fk+SEWHbuK+EvCz/BXdtbc93HWZ1ZAb5aao0zp5SKCJ07d2bXrl188sknpKfXuxKg0CTUYO0zm3L3OdYsE2u/3c2jb31X+w//7lI4dpy1vvIm2PNjCCJUSinn0iQUBGf3aMMFfax7pg+8/g3rf8it/YdH3A8praG0EHL+CO6KEEWplFLOo0koSG4ecQzHtknFbeDK5+vwfKhJBoz2NPH8ZT28p0PCK6UaD01CQZIYF8O/Jveq3/OhToPhpEus9TV/h9+cO3msUkoFkyahIGrQ86Ezb4eMo6CiFF75A5TXsiallFIRTJNQkPk/H3q/ts+H4pvC2CdAXLDzS6tGpJRSUU6TUAhUeT5Ul/5D7XpDv2ut9f89DFvWhS5IpZRyAE1CIeD7fGhHQR2fD/WfAa17AgZyLoWSfSGNVSml7KRJKET8nw+t2vRb7T4YGw9j50BMAuRtsfoPKaWC5r333mPkyJFkZ2eTkJBA69at6d+/P3fffXeVciLC1VdffdjjzZo165BJ7XJzcznvvPPIyspCRLjooosoLi5m1qxZvPXWW7WKc968eYhIlVlbH3/88TpPvOd0EZ2ERGS2iGwSEbeITLQ7Hn9n92jD6V2yAHh547baf7DF0TD4Vmt943z4+r8hiE6pxue1116jX79+FBcXc//997Ny5Uruu+8+evTowcsvv1yvY06fPr1y1G2v22+/naVLl/LII4+wbt06Zs6cSXFxMbNnz651EhoxYgTr1q2jdevWlduiMQlF+qR23wJXAbfbHUh1xp+Yzdvf7GLN1zvZs7+UjKbxtfvg7/4Em5fDT2th8WXwx7WQfkRog1Uqyt1777106tSJFStWEBMTU7l9ypQpuN3ueh0zOzub7OzsKtu+/PJLOnfuzIQJEyq31XWMuaysLLKysuoVU12UlJSQkJAQ8vNUJ6JrQsaYZ40xrwPFdsdSnSHdWpKSEEu527D0019r/0GXC8Y9CU0yrfmHFl0E5aUhi1OpxiA3N5esrKwqCcjL5Qr853DevHkcffTRNGnShJ49e/Laa69V2e97O+6nn35CRHjjjTf4/PPPERFEhHnz5tGsWTMAZs+eXbl91qxZ1cbqfzuuffv2fPrpp7z99tuVnx8wYEBl+a1btzJt2jRat25NfHw8Xbp04R//+EeVY7711luICAsWLODyyy+nZcuWJCY2bJr0hgpKTUhEsoHrgROB44GmwEBjzFsByiYDdwLnAunAl8BtxphXgxGL0yTGxTD8uNYs3PALr2zcytRT2tf+w6mt4Zwn4ZlxsO0jeONWGHpXyGJVqkYV5dY8WE6R2hZi6vYnrE+fPsydO5drr72WyZMn07NnT2Jjqz/GkiVLWL9+PbfffjvJycncc889jB07ls2bN9OxY8dDyrdu3Zp169Zx2WWXUVhYyPz51sDExxxzDNnZ2QwZMoRLLrmE6dOnAxxSg6pJTk4OEyZMIDk5mUcffRSA1NRUAH799VdOPvlkkpOTufPOO2nXrh2vv/461113Hbm5udx+e9WbRTNmzGDQoEHMmzePffvsbfwUrNtxnYBJwEbgTWBUDWVzgF7ADOBH4CIgR0RGGmOWBSkeRxnXqy0LN/zCp1vz+W5nIZ1a1GESrKMGQf/r4Z17YP2jcOQpcMzI0AWrVHUKtsFDPeyO4qCrPoNmdZvz6K677uL777/ngQce4IEHHiApKYlTTjmFMWPGcOmllxIXF1elfGlpKW+++SZNmzYFoFevXrRp04aFCxdy4403HnL8hIQE+vTpQ2pqKm63mz59+lTuO+mkkwAr8fhur60TTjiBJk2akJqaesjnZ82axYEDB/joo48qnyENHjyY0tJS7r33Xv7yl7+QkZFRWb5bt24888wzdY4hFIJ1O+4dY0wLY8xQ4D/VFRKR4cBgYLoxZq4xZjUwFVgH3O9XdrmIFFazzAhS3GHRu30G2c2SAMj5uB5Teg/4K7TvZ60v/rOOtq1UPWVlZbFmzRo2btzIPffcw/Dhw9m4cSNXXHEFffv2pbi46p39QYMGVSYggJYtW9KiRYvKyd6cYtmyZZxxxhlkZWVRXl5euQwfPpySkhLWr19fpfzYsWNtivRQQakJGWNq+0RvLJAPLPH5rBGRp4E5ItLNGLPJs31YMGIDa+bUwxRJC9a5AnG5hHEntOXh1d+Rs3Eb1w7pisslh/9g5QFi4Jy58PhpsH+n9XzoklUQa9/DRNUIpba1ah9Okdq23h894YQTOOGEEwAoKipi+vTpPPfcczz11FNcdtllleWaNz90tuSEhIRDkpXdduzYwSuvvHJITc5r9+7dVd77trizW7hbx3UHNgVIWp/57q/twUQkDojBqtHFiUgiUFqHpBg2Y3tl8/Dq7/g1v5j1P+ZyylGZdTtASksYPxfmj4btn8DKm2HEfaEJVqlAYmLrfPsrEiQlJXHDDTfw3HPP8cUXkTl4cGZmJr169WL27NkB93fo0KHKe/9+TXYKdxJqDnwTYPsen/118STW7TyAfsB8YCDwlm+hw81/7qkphbQ21CGzKScckc7HP+fxysZtdU9CAB36w4AbrXHlPnzSej7UfVzwg1UqSm3fvj1gLeCrr74CQltD8DaDLioqatAxAn1+2LBhvPHGG3Tp0qWysUKksKOfUE3j19RybBtPYWMuwmrYEBHG9crm45/zWP75dm4bfSxN4utx+ftdC1vegx/WwKtXWkP8ND8q+MEqFYWGDh1Ku3btGDVqFF26dKGsrIwPPviA++67jxYtWlS2WguFpKQkOnbsyNKlSxkyZAjp6em0adOGNm3a1PoY3bt3Z8GCBSxatIgOHTqQkpJC165duf3221m1ahWnnnoqV155JZ07d6awsJDvvvuOpUuXsmrVqoDN0p0g3P2Ecglc2/E229gTYF/UGNmjNXExwv7SClZ9uaN+B3HFWP2HUlpD6T54cSqU1f+blVKNyc0330xKSgr33nsvZ599NmeffTZz587lvPPO48MPPwz5s5I5c+YQFxfHiBEj6N27N3PmzKnT52+99Vb69evHxRdfTO/evbn00ksBaNu2LRs2bKBfv37ccccdnHnmmVx88cW8+OKLDBgwoNo+UE4gxtSp8nH4A4qMwWqGfUg/IRH5N3AO0Nz3uY2ITMe6tXast2FCOIlIXlpaWlpdezTXxx+f+YgVX/5G/y5ZzJ92cv0P9NP/4OmRYCqgxwTPNBDOuc+rIpe35deRR0bf8x9VN4f7XUhPTyc/Pz//cI88ahLu9JiD1UHVv6PLhcBmOxJQuI3rZbXoeffbXewoaEALm/anwpDbrPXPFsL7jwchOqWUCq+gJSERGS8i44G+nk2ne7b5NrVeBqwB5orINBEZKCLzgNOwRlyIegO6tqBZkzjcBpZ80sDe533/DN3HW+srb4Yf32l4gEopFUbBrAkt8izejqSzPO8f8xYw1r2/McALWEP3LAd6AOOMMUuDGItjxce6GNnTehD58kfbaNDtUBEY9U9oeZx1W27RRZD3S3ACVUqpMAhaEjLGSDVLe79yBcaYy40xrYwxicaYXsaYxcGKIxKM62WNF7V5xz42bS9o2MHim8DEZyGpGRzIhYVTtKGCUipiOLfJRBTrmZ1GxyxrKJCcuswzVJ1m7WH8UyAu2P4pvPYXCHKDE6WUCgVNQjYQEc7x1IYWf/Ir5RVBGODhqEEweJa1/unz8P4TDT+mapREpN5z66joUlFREfLm3ZqEbDL6eOu50O7CEtZ+t/swpWvplCvhWM8ICitvgp/eDc5xVaOSkJBASUkJFRUVdoeibOR2uykrK6txqotg0CRkk+xmTejT0eqj+0owbsmB1VBh9CPQsrvVUOHFqdpQQdVZamoqxhi2b9+uiaiRMsawe/du3G43aWkhHdEs4qf3jmijj2/L+h/28M43u4J30PimMOFZmDMADuyGhefDtBUQlxS8c6iolpSURMuWLdmxYweFhYUkJCQ4use9Cr7y8nJKS0tJTU2lSZMmIT2X/mbZ6OhWKQDkF5VRUFwWvANndPBpqPAJLLlcGyqoOsnIyKB9+/ZkZGQ4dswxFTrx8fG0atUqLFM+aE3IRtnNDn7D2La3iNTWgecCqZdOZ8CQ22HVzfDFS9CqO5z2l+AdX0W9pKQkkpK0Bq1CS2tCNspMjich1von2Lo3BH17+v4Zek6y1t+YDd+sDP45lFKqATQJ2UhEaOuZ9nvb3gOhOAGc/Q9oexJg4OXpsGtz8M+jlFL1pEnIZt5bciGpCQHEJVoNFZJbQUkBPD8JikI/WrhSStWGJiGbZXtqQiFLQgCprWHicxCTAHu+h5cvAbc2vVVK2U+TkM0qk1BeCG7HVTnRiTDqYWv9uzfgjVtDez6llKoFTUI2a5sehpqQV8+J0Pdya/29f8KnL4T+nEopVQNNQjbzPhPKO1BGYUl56E845DY46gxr/dUrYetHoT+nUkpVQ5OQzdo1O9gPY1s4akOuGBg/FzKOgooSeGEyFGwP/XmVUioATUI2y0xOIL6yr1CInwt5JTWDSS9AQhoU/mYlIp2DSCllA01CNnO5hOxwPhfyyupi1YjEBb9utG7N6dA+Sqkw0yTkAJUdVvPCXBvpPMR6RgTw+Yvwv3+E9/xKqUZPk5ADHOwrFKbbcb76Xg49J1vrb8yGzSvCH4NSqtGK+CQkIn1FxC0if7U7lvoK+agJNRGBsx+E7N5UDu2z86vwx6GUapQiOgmJiAv4B/CB3bE0RFhGTaiJd2iflDZQug+enwgH9tgTi1KqUYnoJARcCvwP+NruQBrCm4T27C/lQGkY+goFktIKJj0HsYmw9ydYNBUqgjjHkVJKBdDgJCQi2SLykIi8KyKFImJEZEA1ZZNF5GER2S4iRSKyQURG1fO8zYGrgNkNCN8R2qZXnVfINm1OgDGPWus/vgMrbrQvFqVUoxCMmlAnYBJQCLx5mLI5wBTgFmAEsAnIEZHh9TjvXcADxpj8enzWUVqkJBAXI4CNt+S8up8D/a6z1j98EjY8ZW88SqmoFowk9I4xpoUxZijwn+oKeRLNYGC6MWauMWY1MBVYB9zvV3a5p1YVaJkhIicCvYB/ByF+27lc4jOGnA0t5PwNvBm6jrDWl10PP661Nx6lVNRqcBIyxrhrWXQskA8s8fmsAZ4GjhaRbj7bhxljkqtZ7gH6Ad2AnSKyG5gIzBSRgF/bRSSvpgVIq9cPH0S2tpDz53LBuCegRTdwl8OLF8KeH+2OSikVhcLZMKE7sClA0vrMZ39t/RvoAhzvWV4FHgKua2iQdqmsCYW7w2p1ElJg0vPQpDkU7bGG9inZZ3dUSqkoE84k1BwI1O53j8/+WjHGFBpjtnoX4ABQYIwJ2K7YGJNe04JVQ7OV7c20A2nWHs6bD65Y2LkJXrkU3LWt+Cql1OGFu4l2TYOT1XvgMmPMRcaYu+v7eSfIzvAM3eOEZ0K+2p8Gw++z1jf/F9bcYW88SqmoEs4klEvg2k6G57VR9470PhPaXVhKUanDpt4+6WI4+Q/W+tr74fOX7I1HKRU1wpmEvgSO8Yxy4Os4z+sXYYzFcbJ95xUK9VTf9XHWXdDhdGt9yZ9hm06Gp5RquHAmoRwgHRjpt/1CYLMxZlMYY3GcFimJxLoc0lcokJhYOHceZHSE8mJ4YYpOhqeUarCgJCERGS8i44G+nk2ne7YN8ym2DFgDzBWRaSIyUETmAacB1wcjjkgW4xLa2DGvUF00yfBMhpcK+7bDwik6GZ5SqkGCVRNa5FlmeN7P8rx/zFvA0ydoDPACcCewHOgBjDPGLA1SHBHNkS3k/GV1hXPmAmLdktPJ8JRSDRCUJGSMkWqW9n7lCowxlxtjWhljEo0xvYwxi4MRQzSwdV6huuhypk6Gp5QKikgfRTuqeFvIhX2G1fo45QroOclaf2M2bF5ubzxKqYikSchB2jr9mZAvETj7H1Unw9vRqNuWKKXqQZOQg3hvx+3aV0JxmcP6CgUSlwgTFkBqWygttCbD259rd1RKqQiiSchBsjN85hWKhFtyACktYeJzEJsEeVt0MjylVJ1oEnKQlikJzu4rVJ02xx+cDO+ntbB8Rs3llVLKQ5OQg8TGuGiVlgjYPMNqfXQfB6ffYK1veAo+eNLeeJRSEUGTkMNETDPtQE7/Kxzjma19+Q3ww9v2xqOUcjxNQg7jqMnt6srlgrGPQ8vjwFRYk+Hlfm93VEopB9Mk5DARXRMCiG8Kk56DJplQnAfPT4Ji26drUko5lCYhh4moDqvVST8CJi4AVxzs3mz1IXJHQJNzpVTYaRJyGG+H1R0FJZSUR/Af7iP6wEjPcD7froI3ZtkajlLKmTQJOYzvvEK/5hXbGEkQnHA+9Pmztf7ew/DJ8/bGo5RyHE1CDtM6LZGYyr5CEfpcyNeQ2+CoM6z1pVfCLx/YG49SylE0CTlMbIyLVqlWX6GIbCHnLyYWxj8FzTtDRak1GV7+VrujUko5hCYhB/Lekou4DqvVSUqHyQshMQ3274QXJkNpFNTylFINpknIgdpGejPtQJofZU0PLjGw/VNY/CedDE8ppUnIiSK6w2pNjhoEQ++y1jcthnfutTcepZTtNAk5UERM811fJ/8Bek211tf8HTYtsTcepZStNAk5kDcJ7dhXTGm52+ZogkwEht8HR55qvc/5I2z/zN6YlFK2iegkJCKFfosRkXPsjquhstOt23HGwPb8KKwNxcbDefOtkRXKDlhD+xTutDsqpZQNIjoJGWOSvQtwClAErLI5rAZrlZaIp6tQdN6SA2iaCZNegPhkKNgKC8+H8hK7o1JKhVlEJyE/U4Alxph9dgfSUPGxvn2FoqiFnL+Wx8K4OYDAL+/Da9doizmlGpmgJCERyRaRh0TkXZ/bYgOqKZssIg+LyHYRKRKRDSIyqoHnF2ASsKAhx3GSqG0h5+/oETDoFmv9k2dh/aP2xqOUCqtg1YQ6YSWBQuDNw5TNwaq13AKMADYBOSIyvAHn7w80AVY24BiOEnUdVmvS71roPt5aX3ULfPu6vfEopcImWEnoHWNMC2PMUOA/1RXyJJrBwHRjzFxjzGpgKrAOuN+v7PIADQ+8ywy/Q08BXjTGlAXp57Fd22hupu1PBEY/Am1OAOOGl6bBrm/sjkopFQZBSULGmNq2Ix4L5AOVnUOMMQZ4GjhaRLr5bB/m2/DAb7nHW05E4oHx1HArTkTyalqAtLr9xKEX8ZPb1VVcEkx8DpJbQUkBPD8BDuyxOyqlVIiFu2FCd2BTgKT1mc/+uhqOldjea0hgTuN9JvRbQRT2FapOahsrEcUkwJ4f4KWLoaLc7qiUUiEU7iTUHAj09XaPz/66mgI856lRBWSMSa9pwUpijuKtCbkN/JYf4fMK1UX2iTD6X9b6D2/ByptsDUcpFVp2NNGuqQ1undvnGmPONcbc3IB4HKl1WhLi7SuU10huyXn1OBdOu8Za/+AJ+GiereEopUIn3Ekol8C1nQzPqz4E8IiPddEyJYrmFaqrQTOhq6fB5H+vhZ/etTcepVRIhDsJfQkcIyL+5z3O8/pFmONxtKgeyPRwXC6rI2uLbuAuh4UXwN6f7I5KKRVk4U5COUA6MNJv+4XAZmPMpjDH42gtUhMA2F3YSIezSUiBSc9Dk+ZQtMcaY64k4gfEUEr5CFoSEpHxIjIe6OvZdLpn2zCfYsuANcBcEZkmIgNFZB5wGnB9sGKJFpnJVhLKbaxJCKBZe2uwU1cs7NwEr/wB3I2ktaBSjUAwa0KLPIu3I+ksz/vHvAU8LdjGAC8AdwLLgR7AOGPM0iDGEhWaN/XWhEptjsRm7U+DEZ6+zJuXwerb7Y1HKRU0QUtCxhipZmnvV67AGHO5MaaVMSbRGNPLGLM4WHFEk8yUeKCR14S8TrwITr7UWn/3Afhska3hKKWCI5pG0Y46WhPyc9ad0HGgtb7kz7D1I3vjUUo1mCYhB8tMtmpChSXlFJdV2ByNA8TEwrn/gYyjoKIEXpgMBb/aHZVSqgE0CTmYt2ECNOIWcv6SmlmT4SWkQeH6EMiiAAAVDUlEQVRvViIqa4RN2JWKEpqEHKy5pyYEkKu35A7K6gLjnwJxwa8fW7fmdDI8pSKSJiEHS06IJSHW+ifSmpCfzoPhzDus9S9ehrX311xeKeVImoQcTER8+gppTegQfS6D48+31lffDl//1954lFJ1pknI4byNE3ZpTehQInD2A9Dud9b7l38Pv+nIT0pFEk1CDtdca0I1i02ACc9CWjso228N7bN/t91RKaVqSZOQw3lrQrn7tSZUreQW1mR4cU0g/2drsNNyTdpKRQJNQg7nrQlpw4TDaN0Dxj5hrf/8Hiy7VlvMKRUBNAk5nDZMqINuo2CgZ37DjfPh/SfsjUcpdViahBzOeztOa0K11P96OHastb7yRvjuTXvjUUrVSJOQw3lrQnv2l1Lh1ttLhyUCox+F1j3BuGHRxbD7W7ujUkpVQ5OQw3lHTXAb2HtAb8nVSnwTmPg8JLeEknx4fiIU7bU7KqVUAJqEHM47kjboc6E6SWsLExZATALkfgcvTYOKcrujUkr50STkcBlN4xGx1vW5UB216w2jHrbWv18Nr8+0Nx6l1CE0CTlcjEvIaKKNE+qt50Q49Sprff2jVqs5pZRjaBKKAJnJOrldg5xxK3QZaq2/dg1sec/eeJRSlTQJRQBv4wSd5rueXDEw7knIOgbcZbDwfNi7xe6olFI4PAmJyGwR2SQibhGZGGB/SxFZJSIHROQzETnRjjhDLVNHTWi4xFSY9Lw1Kd6BXGsyvJJCu6NSqtFzdBICvgWuAj6oZv9jwHdAc+AR4GURiQtTbGFzsCakt+MaJKMDnDcfXLGw4wvIuRTcbrujUqpRc3QSMsY8a4x5HSj23yciKcDZwGxjTJExZo5n1ynhjDEcKmtC+zUJNViH/jDsHmv969dgzd/tjUepRq5WSUhEskXkIRF5V0QKRcSIyIBqyiaLyMMisl1EikRkg4iMCmrUls5AnjFmh8+2z4FuITiXrSqH7tmnt+OCovcl0Pv31vra++Dzl+yNR6lGrLY1oU7AJKAQONxgXDnAFOAWYASwCcgRkeH1DbIaTYECv20FQHKQz2O7ykFM95dgdGTo4Bh6l1UrAljyZ9j2kb3xKNVI1TYJvWOMaWGMGQr8p7pCnkQzGJhujJlrjFkNTAXWAff7lV3uqVUFWmbUIqb9QIrftlSsRBlVvNM5FJe52V9aYXM0USImDs59Gpp1gPJieGEKFGy3OyqlGp1aJSFjTG2f3o4F8oElPp81wNPA0SLSzWf7MGNMcjXLPbU417dAMxFp6bOtO1bNqwoRyatpAdJq+fPZwns7DrSZdlA1yYDJCyEhFfZtt1rMlRXZHZVSjUqwGyZ0BzYFSFqf+eyvNRGJE5FErDjjRCRRRFwAxph9wGvATM/26YAAUdcT0Xf8OG2mHWRZXeGcuYDArxvh1St0MjylwijYSag5sCfA9j0+++viSaAI6AfM96z399n/J6ArsBerKfc5xpgy/4MYY9JrWrBqb46VFB9D0/gYQEdNCIkuZ8KZt1vrny+Cdx+0Nx6lGpFQNNGu6Wtknb5iGmMuMsaI3/KWz/4dxpghxpgkY8xxxpiofbqcmaIdVkOq7+Vw/BRr/c3b4Otl9sajVCMR7CSUS+DaTobnNVAtSdVC86baYTWkRODsByH7ZMDAK7+HHV/aHZVSUS/YSehL4Bjvcxsfx3levwjy+RoNHbonDGITYMKzkJoNpYXWZHj7c+2OSqmoFuwklAOkAyP9tl8IbDbGHNJyTdWOt5m21oRCLKUlTHoO4ppA3s/w4oVQrtdcqVCpdRISkfEiMh7o69l0umfbMJ9iy4A1wFwRmSYiA0VkHnAacH2wgm6MsjzNtHdpTSj0WveEMY9Z61veheXXa4s5pUIktg5lF/m9n+V53QK0B6tPkIiMAe70LOlY/XbGGWOWNijSRu5gTUiTUFgcOwZ23Qhv3QUfzYMWx8Lv/mB3VEpFnVonIWOM1LJcAXC5Z1FBcnDoHr01FDb9Z8DOTbBpCaz4K2R2hqMG2h2VUlHF0aNoq4O80znkHSijrEKnHwgLl8u6LdeqB5gKWDQVcr+3OyqlooomoQjhrQkB7NHaUPjEN7Umw2vaAorzrRZzxY7u26xURNEkFCF8x4/bpVM6hFdaNkxcADHxsPsbeGkauHUgWaWCQZNQhEhLiiPWZT2W0+dCNmh3Mox82Fr/7g14/W/2xqNUlNAkFCFEpPK5kE5uZ5PjJ8EpV1rr6x6Bj5+1Nx6looAmoQjiHU07d78mIdsMngWdz7LWl14NP6+3MxqlIp4moQhycBBTvR1nG1cMnPNvyDoa3GWw8HxrZAWlVL1oEoogmZ5BTHX8OJslplot5pKawf5d8PxkKIm6CX2VCgtNQhFEa0IOktERzpsPrljY8Tks/iO4tf+WUnWlSSiCHJzOQWtCjtChPwz7P2v9q6XWED9KqTrRJBRBMnUkbefpPR1OusRaf+ce+OIVe+NRKsJoEoog3ibauftLMDqqs3MM+z9o389aX3wZ/PqxvfEoFUE0CUUQb02orMJQUFRuczSqUkyc9XyoWQcoL7IaKuz7ze6olIoImoQiiO/4cTqvkMM0ybBazMWnwL5f4YUpUFZsd1RKOZ4moQiS0fTg+HHaOMGBWhxj9SFCYNsGWHqVToan1GFoEoog8bEu0pLiAG2m7Vhdh8LgW631z16A9x62Nx6lHE6TUITJ9GmcoBzq1KuhxwRr/fVb4ZuV9sajlINpEoow3mm+dRBTBxOxRtxueyJg4KVLYOfXdkellCNpEoowWd4kpNM5OFtcIkx8DlJaQ+k+azK8A3vsjkopx9EkFGF0OocIktLKSkSxibD3R3jxQqgoszsqpRzF0UlIRGaLyCYRcYvIxAD754jIdhEpEJHPReRsO+IMp4PTOWhNKCK07QWj/2Wt/7QWVvzV3niUchhHJyHgW+Aq4INq9j8AtDfGpALTgGdFJD1cwdkhM0VH0o44x42HftdZ6x/+Gz6ca288SjmIo5OQMeZZY8zrQMBef8aYr40x3r/G5UAC0DZc8dmhsiakTbQjy8Cb4WhPRX35DPhxrb3xKOUQtUpCIpItIg+JyLsiUigiRkQGVFM2WUQe9twmKxKRDSIyKqhRVz3foyJSBGwE3gQ2hepcTpDlqQkVlpRTXFZhczSq1lwuGPsEtDgW3OXw4gWw50e7o1LKdrWtCXUCJgGFWH/oa5IDTAFuAUZgJYUcERle3yBrYoy5DEgGhgCrTJSP7OmtCYHekos4CcnW0D5NmkPRXnh+EhQX2B2VUraqbRJ6xxjTwhgzFPhPdYU8iWYwMN0YM9cYsxqYCqwD7vcru9xTqwq0zKjLD2GMqTDGvAGcKSJnBYgrr6YFSKvL+ezkndgO9JZcRGp2JEx4FlxxsOsreOX34NYarWq8apWEjDG1nTJyLJAPLPH5rAGeBo4WkW4+24cZY5KrWe6pw8/gKwY4qp6fjQhN42NIiLX+2bQmFKGOPAVGeL6TfbMC3rzN3niUslGwGyZ0BzYFSFqf+eyvNRGJE5FErDjjRCRRRFyefckicr7nNVZEzgUGAu/6H8cYk17TgpU4I4KI6OR20eDEqfC7P1nr//sHfLrQ3niUskmwk1BzIFC38D0+++viSaAI6AfM96z39+wzWM2ytwK5wF+BScaYzwIcJ6p4x4/T6Rwi3Jl3QMeB1vqrV8Av1fVEUCp6haKJdk0NA+rUaMAYc5ExRvyWtzz79htjBnlqM2nGmBONMTkNCTxSaE0oSsTEwrn/geadoKLEaqiwd4vdUSkVVsFOQrkEru1keF518KwgqBy6R2tCkS+pGUx+0Xo9sBuemwDFEXN3WKkGC3YS+hI4xvvcxsdxntcvgny+RqmyJqTTOUSH5kdVbTG36CKo0OnbVeMQ7CSUA6QDI/22XwhsNsZEdUfScDk4nYPejosa7U+DUZ4J8L5fbY2qEN1d3pQCILa2BUVkvGe1t+f1dBHJBPYbY5Z7ti0D1gBzRaQ58CNWP6HTgNHBCVnpxHZR6vjJkPsdrL0fNsyFzM7Q5092R6VUSNU6CQGL/N7P8rxuAdqD1SdIRMYAd3qWdKwRE8YZY5Y2KFJVyXs7bs/+UirchhiX2ByRCpqBt1iJaNMSWHEjNOtgTRmuVJSq9e24AK3UvEt7v3IFxpjLjTGtjDGJxphexpjFQY+8EfM2THAb2HtAb8lFFZcLxjwObXphzco6DX773O6olAoZR4+irQLz1oRAm2lHpfgmMOkFSGsHZfutFnMF2+2OSqmQ0CQUgZo1icd7B06baUeplJZWIopPhoJt8Pip8MZsyPvF7siUCipNQhEoxiVkNNW+QlGvVXc4d541PfiBXHj3AXioB7wwBX54W1vPqahQl4YJykGaN01gd2Epu/V2XHTrPASu+hQ+mgcbnoLCHfD1a9aS2RVO/j30nAgJKeB2g7sMKsqgotSat6iiFEoKrSRWZdljvRbngSsWYhOsZBebADEJB9+LyxrNobwEyos9S6nntQRcMVU/G5sIMfGe1zjr/N6y/q9w6Gd9j+Eur/6z7gqIjff7rE/8xn2wbJX4S6yYYnw/63du5NDzVZQc/Nn9r1ds4sFYXLGB4/W+iquGn9l7vUqquV7GL26/1yrXy+8Y7vKDMcYE+JljE+DkP1i3gsNMk1CEykyJZ/MOyNWaUPRLaQUD/gqnXQNfL4UPnoSf18HuzbDsOlh+g1XO6JQQqgFOutiW02oSilDeye30dlwjEhsP3c+xlu2fwYdPwmeLoLyodp9PSIMmGdaket4lKd2qVfjXFrzf+k3FobUjb+3hsDWOMp9v395v8D6fh6rf2Kt8vtSqGfjXrqrU0PxqWb7vA9U4AtbQAhzDuAPXkmISrGOYisPUOPxqR5XnTwBMgFqS99xlnp+5hhpaRYB/o8paaQ01NIk5+DNXVzuNTQzVb26NNAlFKB3EtJFr3QNG/RMGz4ZfN1p/dGPiraF/YmJ91uMgvikkZVh/kJRyGE1CEUoHMVWAVbPpNNjuKJSqN20dF6GyvOPHaU1IKRXBtCYUoXxrQsYYRKyOQyXlFeQfKCO/yFrc2opXKVULvY5IJzYm/PUSTUIRyjuSdkm5m2EPrSW/qIy8A2UUlWkLKaVU3X0+60xSNAmp2mqbnkSMS6hwG77+bZ/d4SilVL1oEopQWSkJPDalF19syyetSTzpSXGkJcWR3sRa0pLiSU2KJdalj/2UUodn12D8YnToD0QkLy0tLS0vL8/uUJRSKmKkp6eTn5+fb4xJr+8x9GuyUkop22gSUkopZRtNQkoppWyjSUgppZRtNAkppZSyjSYhpZRSttEm2oCIuAFJS0uzOxSllIoY+fn5AMYYU+8KjSYhQETKsWqFBfU8hDd75Qcnoqin16tu9HrVjV6vumnI9UoF3MaYeg98oEkoCEQkD6AhHbYaE71edaPXq270etWN3ddLnwkppZSyjSYhpZRSttEkpJRSyjaahJRSStlGk5BSSinbaBJSSillG01CSimlbKP9hJRSStlGa0JKKaVso0lIKaWUbTQJKaWUso0moXoSkWQReVhEtotIkYhsEJFRdsflBCKSLSIPici7IlIoIkZEBlRTdrKIfCoixSKyVUTuFpHEMIdsGxE5Q0TmichmETnguQaviMhxAcoOEZH1nt+3nSLyhIg0uvHRROQUEVkpIts8vze7RGS1iAwLUFavmR8RmeX5P/lJgH1hv16ahOovB5gC3AKMADYBOSIy3NaonKETMAkoBN6srpCInA8sAP4HDAPuBP4MzAt9iI7xR+AI4EGsa3CN5/2HItLHW8iTxJcBvwAjgeuAUcB/RaSx/T9uBmwGrgWGAn8ASoBlIjLRW0iv2aFE5FjgBmBHgH0DsON6GWN0qeMCDAcMMNZnmwDvAl/ZHZ/dC+DyWR/juVYD/MrEANuBJX7bf+8p/zu7f44wXasWAbalA3uBl322fQB87Hdth3iu1QS7fw67FyAW64/nar1m1V4jF7Ae+CfwFvCJ335brlej/DYQBGOx5t5Y4t1grH+xp4GjRaSbXYE5gTHGXYtifYBWWNfM1wKgDDgn2HE5kTFmZ4BtecC3QDaAiLQFegPP+F5bY8zrwDYaybWqiTGmHOv/ZBnoNavGX7B+p27232Hn9dIkVD/dgU0B/th+5rNf1cx7jb7w3WiMOQB8TyO+hiKShfXze69NwGvl8TmN9FqJiEtEYkWkjYjMBrpg3dYEvWZViEhH4DbgcmNMoMk7bbtemoTqpzmwJ8D2PT77Vc2816i669gor6GICDAH6//mfZ7Neq0CexGr5rMNuBo4zxizwrNPr5mH53fqSWClMWZxNcVsu16ahOqvpqEmdBiK2qvuWjXWa3gv1nO0PxpjvvLbp9eqqhnAyVgPz5cBL4rIJL8yes2s56wnAVfUomzYr1e95wVv5HIJ/M0gw/Ma6NuEqirX89rcZ90rA/gxvOHYT0T+jtXi6ypjzDyfXb7Xyl8GjfT3zRjzA/CD5+1SEVkK/EtEFqLXDAARyQTuAe4C9vs0t44FYjzvi7HxemlNqH6+BI4J0GzR27cj0H1VVdWXntcq95pFpAlwFI3sGorIbcBNwAxjzMN+uwNeK4/jaGTXqgYfYDXfzkKvmVc2kIaVhPb6LKdiXZu9wCxsvF6ahOonB6sZ7Ui/7RcCm40xm8IfUsRZD/wGXOC3fRIQB7wS9ohsIiK3AjOBmcaYe/33G2O2AhuAKb5ffETkDKAtjehaVcfz3GMAkAfk6jWr9B0wMMDyKVYDoIHAHDuvl96Oq59lwBpgrog0x7p1NBU4DRhtZ2BOISLjPau9Pa+ne24N7DfGLDfGlIvIX4F5IvII8BJwDPB/wEvGmPXhjzr8RORarG+irwFv+HZQBUqMMR971m8AVgHPi8gcoA3WtXofWBS+iO0nIguALcBHwG6gNdb/v0HAFZ7m2qDXDGNMIVafoCpEJM+z33efPdfL7g5UkboAqcAjWN/mi4GNwBi743LKgvUgM9Dyk1+587GagJZgtXK6B0iyO/4wXqe36nCthnr+IBQDu7BaPDWz+2ew4ZpdDqzDeo5R7nldCYwMUFavWeBr+BZ+nVXtul46n5BSSinb6DMhpZRSttEkpJRSyjaahJRSStlGk5BSSinbaBJSSillG01CSimlbKNJSCmllG00CSmllLKNJiGllFK2+X9elhBFjty8pwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "res = np.array([np.linalg.norm(A @ x - (x @ A @ x) * x) for x in conv_rq])\n",
    "plt.semilogy(res, label=\"RQ\")\n",
    "plt.semilogy([np.linalg.norm(A @ x - (x @ A @ x) * x) for x in conv_l], label=\"Shift iter\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Как найти другие собственные значения?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1.71428299e-15  1.00129089e+00  1.14452777e+00  1.72651117e+00\n",
      "  5.59525073e+00  8.50188201e+00  1.22038786e+01  1.31768423e+01\n",
      "  1.59395420e+01  3.22917675e+01]\n"
     ]
    }
   ],
   "source": [
    "B = A - eigv_rq @ (A @ eigv_rq) * np.outer(eigv_rq, eigv_rq)\n",
    "print(np.sort(np.linalg.eigvals(B)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.726511169331873\n"
     ]
    }
   ],
   "source": [
    "eigv_rq_B, conv_rq_B = rayleigh_quotient_iteration(B, x0, num_iter)\n",
    "print(eigv_rq_B @ B @ eigv_rq_B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1.56067399e-15  3.01144918e-15  1.00129089e+00  1.14452777e+00\n",
      "  5.59525073e+00  8.50188201e+00  1.22038786e+01  1.31768423e+01\n",
      "  1.59395420e+01  3.22917675e+01]\n"
     ]
    }
   ],
   "source": [
    "C = B - eigv_rq_B @ (A @ eigv_rq_B) * np.outer(eigv_rq_B, eigv_rq_B)\n",
    "print(np.sort(np.linalg.eigvals(C)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0012908945395032\n"
     ]
    }
   ],
   "source": [
    "eigv_rq_C, conv_rq_C = rayleigh_quotient_iteration(C, x0, num_iter)\n",
    "print(eigv_rq_C @ C @ eigv_rq_C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### И так далее для других собственных векторов....\n",
    "\n",
    "- Какая правильная стратегия использования метода Релея?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Оптимальный шаг в методе Ричардсона\n",
    "\n",
    "- Решение линейной системы \n",
    "\n",
    "$$ Ax = b $$\n",
    "\n",
    "- Метод Ричардсона \n",
    "\n",
    "$$ x_{k+1} = x_k - \\tau (Ax_k - b) = (I - \\tau A)x_k + \\tau b$$\n",
    "\n",
    "- Оптимальный $\\tau$ для $A = A^* \\succ 0$\n",
    "\n",
    "$$ \\tau^* = \\frac{2}{\\lambda_{\\max} + \\lambda_{\\min}} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Как его получить?\n",
    "\n",
    "- Найдём минимум $\\| I - \\tau A \\|_2 = \\lambda_{\\max}(I - \\tau A) = \\max (|1 - \\tau \\lambda_{\\max}(A)|, |1 - \\tau \\lambda_{\\min} (A)|)$\n",
    "- Две \"v\" образные функции\n",
    "- Нарисуем функцию максимум от них\n",
    "- Минимум максимума – это точка в которое они равны!\n",
    "- $1 - \\tau^* \\lambda_{\\max}(A) = -1 + \\tau^* \\lambda_{\\min}(A)$\n",
    "- Получаем ответ!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Метод сопряжённых градиентов: напоминание\n",
    "\n",
    "- Главный метод решения большой разреженной линейной системы с симметричной положительно-определённой матрицей\n",
    "- Требует знания процедуры умножения матрицы на вектор\n",
    "- Чувствителен к кластеризации спектра – чем менее равномерен спектр, тем лучше сходимость"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A is normal matrix: ||AA* - A*A|| = 0.0\n",
      "Condition number = 396.57461187799817\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "n = 1000\n",
    "\n",
    "\n",
    "# Random\n",
    "A = np.random.randn(n, n)\n",
    "mu = 10\n",
    "A = A.T.dot(A) + mu * np.eye(n)\n",
    "\n",
    "# Clustered eigenvalues\n",
    "# A = np.diagflat([np.ones(n//4), 10 * np.ones(n//4), 100*np.ones(n//4), 1000* np.ones(n//4)])\n",
    "# U = np.random.rand(n, n)\n",
    "# Q, _ = np.linalg.qr(U)\n",
    "# A = Q.dot(A).dot(Q.T)\n",
    "\n",
    "# A = (A + A.T) * 0.5\n",
    "print(\"A is normal matrix: ||AA* - A*A|| =\", np.linalg.norm(A.dot(A.T) - A.T.dot(A)))\n",
    "print(\"Condition number = {}\".format(np.linalg.cond(A)))\n",
    "\n",
    "b = np.random.randn(n)\n",
    "\n",
    "# Hilbert matrix\n",
    "# A = np.array([[1.0 / (i+j - 1) for i in range(1, n+1)] for j in range(1, n+1)])\n",
    "# b = np.ones(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAakAAAEPCAYAAAD4aTuoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8VFXex/HPLwVCTShBCIReBSlS7Ipl7WVtiK6udd1mXZ917evqro9bUB91q6trdwXFrmtXdC2A0nsnEAg1ARJC2nn+uDdhiBOSTCa5M5Pv+/XKa5hzT+78ziTML+feU8w5h4iISCxKCjoAERGRmihJiYhIzFKSEhGRmKUkJSIiMUtJSkREYpaSlIiIxCwlKRERiVlKUiIiErOUpEREJGYpSYmISMxSkhIRkZilJCUiIjFLSUpERGKWkpSIiMSslKADiBdmVoaX1HcEHYuISBxpD1Q45yLKN6b9pOrGzCoAS09PDzoUEZG4UVBQAOCccxFduVNPKgwz+yRMcUV6enpyfn5+U4cjIhK3MjIyKCgoiPgKlO5JiYhIzFJPKgzn3PjqZWaWD+han4hIE1JPSkREYpaSlIiIxCwlKRERiVlKUiIiErM0cCKMGoagt23qOEREgpazrYh/TFvJnacfSIuUpu/XKEmJiEhYc3LyueixrygsKScroxU/Hd+vyWNQkgpDQ9BFpLn7bNlmrnl+FoUl5bRrmULfzDaBxKEkJSIiVTbv3MM9by7kjTm5ALRtmcKzVx3CiOyMQOJRkhIREQp2l/L4Zyv5+7SV7CmrAGBs7w7cc9YwhnRrH1hcSlIiIs1Y3o5iXpi+lme+XMPWwpKq8rtOP5DLDu9NUpIFGJ2SlIhIs7SxwEtOf/10BSV+z6lNi2QuPqwXE8Zk0y8zNgY0K0mJiDQjizfu4LFpq3h19nrKK7ytmtJbpXLhuJ5cfkRvDmifFnCE+1KSEhFJcDuKS3lr7gb+PX0tc9YVVJV3bNOC74/szjXH9adjmxYBRlgzJSkRkQRUWl7Bx4s38fqcXN5bmFd1SQ+gV6fWXHVUXy4Ykx3IBN36UJIKQytOiEi8Ki2v4LXZuTz84TLWbiuqKm+VmszJw7py4biejO3dAbNgB0TUlZKUiEicK69wfL1qKx8t2sSbczewcUcxAGYwrndHzh+TzcnDutK2Zfx95MdfxE1AK06ISDwoKCrlsc9WMuWbHPJ27KkqTzI4fXgW158wIGZG6UVKSUpEJI4U7injg0V5vDEnl0+Xbqa03FUdG9EjnaMHZjJhTDbZHVsHGGX0KEmJiMS44tJyPl68iTfm5vLR4k0Ul+4dBJGWmsQlh/bi4kN70atTMOvrNSYlKRGRGFRSVsFnyzbzxpxc3l+YR2FJedWxFslJjB+Uyekjsjh+cBfaxOG9prpK3JaJiMSZkrIKvlq5lfcX5vH6nFwKdpdWHUtJMo4c0JkzhmfxvaEH0D4tNcBIm46SlIhIgCpH5r0xJ5d35m8kv2hvYkoyOLRvJ84YkcXJQ7vSIUYn3DYmJSkRkSbmnGPOugJen53Lm3Nz2bRz78g8MzioezpnDM/irJFZdImxZYqampKUiEgTWZa3k1dnr+eNORv2mWgLMDI7gzNHZHHa8G4xt35ekJSkREQaUc62Il6fk8vrs3NZkrdzn2ODDmjHmSOzOGN4Fj07JcaQ8WhTkhIRibJVWwqZMjOHDxblsTRv1z7HenRoxZkjsjhzZBaDuwa3mWC8UJISEYmClZt38cacDXy8ZBOzc/L3OdalXUvO8C/ljcrOiJt182KBklQYWmBWRGpTWl7BN2u288mSzXyyZBOLN+57KS8rPY1zDu7BUQM6M6Z3R5ID3uE2XilJiYjUw7drt/PW3A28FbKQa6Ws9DROGtaV4wZ34fB+nZWYokBJKgwtMCsiodbn7+adeRv4z/yNzFyzfZ9jB3Zrz/hBmRw7uAuje3YgSYkpqpSkRETCmL++gP8u38I78zd+5x7TkG7tOX5wFyaMydaovEamJCUi4svN382LM3L4cHEe89fv2OdYxzYtOH5wF045qCvHDuqiwQ9NRElKRJq1gqJSPlm6iWe+XPOdS3nd0tM4fkgXTh3WjXF9OpKSHNtbrSciJSkRaXZ2Fpfy/sI83pm/kfcX5u1zrHPbFpxzcA/GD8zksH6d1GMKmJKUiDQL5RWO9xdu5N0FeXywKI+dxWVVx1KTjWMGZnLRIT05akAmqeoxxQwlKRFJWHvKyvl0yWYmz8zhixVbKQrZkyk12TisX2dOP6gbJw49gIzWzW+F8XigJCUiCaWiwjFzzXZembWet+bmsiOkxwRwSJ+OnDa8GycP60qXdlrINdYpSYlIQli+aSevzFrPq7NyWZ+/e59jR/TvxDmjejCuT0eyO2rIeDxRkhKRuDV/fQEvfbOOr1Zu/c6yRAd1T+f7o7pzura+iGtKUiISV4pKyvho8Sae+mI1M1bvO2S8R4dWfH9kd74/Kov+XdoFFKFEk5KUiMSFyl7TC9PXsqesoqp84AFtOXlYN28h114dNGQ8wShJiUhMqqhwzMrZzufLtvLR4jzmrCuoOpacZBw7qAtXHNFbc5kSnJKUiMSUbYUlTJ6Zw+QZOazcUrjPse4Zrbjs8N5MGJtNeqvUgCKUppTQScrM/gKcgbd6+U5gCnCzc64k0MBEZB/5RSU8P30tny3dwjdrt1MScjmvf5e2HNGvE+MHd2H8wEz1mpqZhE5SwKPAL51zhWaWCUwGfgXcG2xYIuKcY866AqbMzOGVWev3mWjbLi2Fs0d15+xR3RnVs0OAUUrQEjpJOecWVi8C+gcRi4h4Nu0o5rXZuUz5JoelebuqylukJHHmiCyOH9yFIwd0pl2aLudJjCQpM7sVOBgYDfQB1jjnetdQNwm4Hvgx0BvYjNdDuss5Vxim/i3A7Xjbv28Fbo5+C0RkfyoqHP9dsYUXpq/lP/M3UuH2HhvctR3nje7B2aO606lty+CClJgUE0kKuA/YBnwLZNRS90HgOuAVYBIwxH8+ysxOcM5VhFZ2zt0P3G9mQ4BLgI1Rjl1EarBmayHPfb2Wd+ZvIGfb3lUg0lulcsaIblwwpifDurfXfSapUdSSlJmlAsOAIufcknp+ez/n3Er/PPPxej3hXmMocC0w1Tl3bkj5KuBhYCLwfLjvdc4tMrPZwJPACfWMT0TqqKSsgo8W5/GXT1YwN2TYOMDoXh24aFxPTh/RjZYpyQFFKPGk3knKzCYA5wE/cc5t88v6Ae8A/fznrwETnHNlNZ4oRGWCqoMLAQMeqlb+GHA/cDE1JClfMjCgjq8lIvVQOdn2tdnr2V5UWlV+QPuWTBzbk9OGd2PgAVoFQuonkp7UFUBWZYLyTcIbkPAR0Ak4C7gcL3lE01igApgeWuicK/Z7SWMry8ysLXA+3mXBArxe3p3Au+FObGb5tbx2euRhiySmdduL+GBhHi/OXMeiDXu3W08yOKJ/Zy49rDfHDu5CcpIu50lkIklSBwLvVz4xs/bAqcBk59xE/7LfbBonSWUBW5xze8IcWw8cbmYt/HlQDrgIL4G2ADYBLwO/jnJMIs1KWXkFb8zN5c05G/h4yaZ9BkH0y2zDeaOzOXtUd7qma1FXabhIklQmsCHk+WH+ef4N4JwrNbP38S7NRVtrIFyCAigOqVPij/T7Xl1P7Jzb74ANv6el3pQ0W5t2FvPPz1bx2uz15O3Y+9+wfVoKJw/ryoXjejIyO0ODICSqIklSO9n3w/oYvF7L5yFlxUBjXHwuArrUcCwtpI6IRMHGgmI+XrKJDxfl8enSzZSWe90mMziyf2fOG92DUw/qpu3WpdFEkqSWAaeYWUu85HQ+MNc5tyWkTi+8y2vRlgscaGYtw1zy6453KbDBSx6Z2SdhisOOOBRJRGu3FvF/Hy7jlVnr9rmc16ZFMlce2YfzRmfTs5M2D5TGF0mS+gfwL7xkVYo3ofaGanUOARY0KLLwZgAnAuOAzyoLzSwNGAlMa4TXFGk2Fm/cwb8+X83UWeuqek3t0lI4emAmJx54AOMHdiG9tVaCkKZT7yTlnHvKzAYBV/tFj/pfAJjZcXiJ6y/RCLCaF4Hb8JLiZyHlP8K7F/VcNF7EOTe+epnuSUmi2rWnjHfmbeCF6Wv5du3eQa7dM1px7XH9OXd0D13Ok8CYc672WvU5oVkLoBVQWNd5UmZ2Cd4lQvAm67bAG5UH3hJJz4TUfQS4Bm9o+dvsXXHiv8Bx1VeciBYzy09PT0/Pz69tpLpIfFixeRf/nr6WyTPXUbB777ymvpltuPzw3pw/Jpu0VE24lYbJyMigoKCgoLbBaTWJepKKKAjvHtAxNRz+NLRnY2bJeD2pq/F6bFvwelh3Oed2hTtBhPFUd2R6enqykpTEs+LScj5evMnbFmPZ3tvILVOSOGZgJhPHZTN+YBeSNK9JoqShSSriZZHMbDjePKQhQBvn3Al+eW+8e0bvO+e21+Vc4S6v7aduOV4va1JtdUXEs2jDDibPzOHFGTn7bInRPaMVF4zN5qJDetJZi7tKDIooSZnZPXj3hiovVId2x5KAF/B6O480KLqA6J6UJIKKCscHi/L466crmBVyrynJ4JiBmVw4ricnDDlAvSaJaZGs3TcRuANveaFfARcAt1Qed86tNLOZwJnEaZISiWfOOd5dkMef3lvC8k17r4D37dyG88dkc8HYbDq2aRFghCJ1F0lP6jpgOXCWc67EzM4OU2cRML4hgQVJ86QkHu0sLuX5r9fy7Ndr9tkWY/ygTH58dD8O7dtRq0FI3IkkSR0EPFnLpNlc4IDIQhKR+sjN383b8zbw109WsLVw73/LYwdlcvPJgxnSrX2A0Yk0TCRJyvBWIt+fA9i7ll7c0T0piQdzcvL594wcXpyxtmpViJYpSZw0tCs/O7Yfg7sqOUn8i3RZpMNrOugPET+SxllxQqTZm7F6Gw9/uGyfIeTtWqZwxsgsbjh+AF3aa/VxSRyRJKnJwG/N7CbnXLhh4Lfi7S31fw2KTESqFO4p48UZObwyaz3z1u/d7bZvZhsmjs3mh4f11sRbSUiRJKmH8BaV/YO/S68DMLM/AUcBY4Cv8Nb4i0saOCGxYkPBbh58fylvz9vIrj17F3AZmZ3B9ccPYPygTA2GkIQWydp9u83sWLye0g/wtmQH+AXevapngWvquiSSiHzX7pJy/jFtJX/+ZDklZd4t4CSDEw/sykWH9OSoAZ2VnKRZaNCySGbWEW/L9k54W7RPd85tjlJsMUVr90lTyNlWxEvfrOOJz1exM6Tn9PNj+3HRIb3ontEqwOhE6i+wZZEAnHPb8Cb1ikgDbNpRzN1vLODteRurylKTjSuO6MOVR/bRYAhpthqUpESkYXbtKeOxaSv5x7SV7C711tRrn5bCJYf14ooj+tBJ6+lJMxfJskhP1LGqc85dWd/zizQHGwp28/dPV/LSN+uqBkS0aZHM9ScM4Ioj+pCi/ZtEgMh6UpfVctzhTfh1QFwmKY3uk8ayakshv3trEdOWbqak3BsQkZaaxMWH9OK6EwbQPk273oqEiiRJ9amhPANvEMWdwBeELDor0tyVlld4o/U+Xl61VUbHNi34yTF9mTAmm4zWWvBVJJxIhqCvqeHQGmCOmb0LzAU+AB5vQGyB0bJIEi17ysr5y8crePKL1VW733Zs04KrjurDReN6KjmJ1CLqAyecczlm9gZwPXGapESi4bNlm7nrtQWs2lIIgBlcfngfrj9+AOmtdVlPpC4aa3RfHjCgkc4tEtNWbt7Fr19fULW2XpLBlUf24cJxPembqVubIvUR9STlLzB7HN7kXpFmY8XmXTz++Somz8ihzF+WfFTPDO49axjDuutKsUgkIhmCfvR+zpUNXA6MBP7ZgLhE4saarYXc9/Yi3luYR+UCLtkdW3HjCQM5e1R3LV8k0gCR9KQ+wV9UtgYGTAN+GUlAIvGivMLx9JermfTe0qq5TtkdW3H54X34waE9aZmiVclFGiqSJHUP4ZNUBbAdb/2+6Q2KKmCaJyX745zjy5Vbue/tRcxfvwOAzHYtufP0Azl1WFdNxBWJokiGoN/dCHGIxDznHNOWbeFvn6zgy5Vbq8onjs3mllMGazi5SCPQ2n1haJ6UVLcwdwc3vzynqucE3p5Od54+hNG9OgYYmUhiU5IS2Y/S8gp+99YinvlqDeX+iL3D+nbi6mP6Mn6gNhwUaWy1JikzWxnhuZ1zrl+E3ysSuJmrt3HTlDms2VoEQO9Orfn9ucM5pG+ngCMTaT7q0pNKYv+j+WqiPzElLq3ZWsivX1/AJ0v27t955ZF9uPnkQRqxJ9LEak1SzrneTRCHSOCcczz1xWp+/58lVXs7dW2fxh/PH85RAzIDjk6kedI9KRFg085ifvHiHD5f7i1l1LFNC245ZTDnjOquIeUiAVKSkmatqKSMSe8t5ekvV1Na7l3VPnZQJn84bwSZ7bQrrkjQIk5SZtYSb/+o7kDY/83OuacjPX+QNJk38Tnn+HjJJn731iJWbPZWKW+XlsI1x/bn6qP7atSeSIyIKEmZ2RXAH4AONVXBG2wRl0lKEtv2whLuen0Bb8zJBbwtNC49rDc3njBQW2iIxJhIFpg9GW/x2AXA74BJwKvAdGA8cCIwBXg7alE2MU3mTUzOOd6Zv5E7X53P1sISAA7v14nbTh2iVcpFYlQkPambgK3A4c65nWY2CZjtnLsfuN/MrgT+BjwSxThFGmTrrj3cMnUe7y/Mqyq75ZTB/FiX9kRiWiRJ6mDgNefczpCyquFPzrnHzewS4HbglAbGJ9IgFRWOf36+koc/XF61Uvnh/Tpx88mDGZmdEXB0IlKbSJJUG2BDyPNioH21OjOBKyINSiQa1m0v4s5X5/OxPyk3vVUqt54ymInjegYcmYjUVSRJaiMQOrNxAzCoWp10QFPzJTAvf7OOu16bT2GJNyn31IO68r9nD9fACJE4E0mSWsC+SekzYKKZHeWc+8zMhgET/HoiTaq8wvGPaSv5/X8WA9ChdSp3nXEgZ4/qEXBkIhKJSJLUO8BDZpblnMvFG4p+PvCJmW0DOuINQf9t9MIUqd3yTbu4+aU5fLs2H4BxfTryyIWjOKB9WsCRiUikIlnv5e94E3i3ADjnFgLH4yWvLcB7wCnOubgdgi7xxTnHtKWbOecv/61KUBPHZvPU5eOUoETiXCQ785YCedXKvgJOj1ZQInW1cvMurn1hFgtyvc0IO7RO5c8XHczh/TsHHJmIREMkk3kznHP5jRGMSH18s2YbP3x8etXgiDG9OvD784bTL1MrWIkkikjuSW0ws9eBp4D/OOcqohyTyH6VlFXwyEfL+PPHy6lw0DIliQcmjOS04d2CDk1EoiySJLUab6DEecAmM3sWeNo5Ny+agQVJC8zGrrwdxfzsuW/5Zs12AHp1as0/LhnDoK7tAo5MRBpDvQdOOOeGAIfgLX2UirdM0mwz+8bMrjMz3QyQRjF/fQET/v5lVYK6+NCevHXdUUpQIgnMnItkZ3j/m81SgTOBS4GT8XpmpXgj/Z50zr0ajSBjgZnlp6enp+fn63ZcECbPyOGO1+ZTUlZBi5QkHp44ipOHdQ06LBGpRUZGBgUFBQXOuYjWIWtQktrnRF4P6mLgEmAUUOGcS5hNFZWkgrGnrJy7X1/AC9NzAOjdqTWTJoxkdK+adokRkVjS0CQVzSSyFW+ViUXAsCifW5qhLbv2cP2/Z/Hf5VsBb2HYv10ymvZpWtpIpLlocCIxs8F4l/suBrLwVptYjjf6TyQib8zJ5devL2Cbv+/TRYf05PZTh9Cmpf72EWlOIt2ZtwNwIV5yGoOXmHYAj+Pdi/oiahFKs7KnrJx731zIs1+tBaB9Wgr3fn8YZ43sHnBkIhKESCbzvoS3ukQq3hbxHwBPAq8454qjGp00K+UVjv+ZMrdqW/cRPdL52yWj6ZbeKuDIRCQokfSkzgGW4F3Oe8Y5tz66IUlztLuknOv/PYv3/J1zf3x0X3550iBSkiNZXlJEEkUkSepwf60+kaiYt66Anz73Deu27wbg6qP7csspg7Wtu4hEtMCsEpREzeycfH74+NfsKC7DDO487UCuOLJP0GGJSIyIeKiUmSXjbX7YgRp24XXOTYv0/JL43pyby02T57CnrIIOrVN57IdjGNO7Y9BhiUgMiXR0353AjXjbxO+PtpCX76iocNw6dR4vzvQm6HZtn8ZTV4zT8kYi8h2RjO67GfgNUAA8A+QAZVGOSxLY7a/Or0pQI7Iz+PNFo+jRoXXAUYlILIqkJ/UjYD1wsHNuc5TjkQS2vbCEe99cyNRZ3oDQ4wd34c8/OJi0VHW4RSS8SJJUNvCYEpTUx87iUq54agaz/O3djx2UyWM/HENSkkbwiUjNIpmEkkccrMtnZi3N7DEzW2lmu8xsmZndEHRczdHijTs46cFpVQnqf04cyD8vHasEJSK1iiTZTAbONrOWzrk90Q4oilKAjcCJwEpgOPCumW1wzr0YaGTNSGl5BbdOnUduQTEtUpL47VnDmDA2O+iwRCRORNKTugvYALxkZjE7ocU5V+icu9M5t9w5V+Gcmw28BRwRdGzNxbrtRXz/z/+t6kE9PHGkEpSI1EskSWoB0Bs4FVhuZtv8S2rVv1bU9YRmdquZTfG/z5nZ6v3UTTKzG81ssZkVm1mOmU0ysza1vEYKcCQwt65xSeSKS8u55vlZLMjdAcBN3xvIycO6BRyViMSbSC73JeENOV8bUhbu5kJ9bjjcB2wDvgVq2xjrQeA64BVgEjDEfz7KzE5wzlXU8H0P4w2bf7oecUkENhTs5qqnZlYlqEcvGsXpw7MCjkpE4lEkyyL1boQ4+jnnVgKY2XygbbhKZjYUuBaY6pw7N6R8FV4Smgg8H+b7JuH1oo5zzpVEP3ypNG3pZv5nyhw27dxDcpLxq5MHKUGJSMRiYonpygRVBxfi9dAeqlb+GFCEt/HiPszsIbzBE8c757Y0JE7Zv+WbdnHlUzPYtHMP7dJSePLysVx9dL+gwxKROBaNnXk7AG2dczlRiKc2Y4EKYHpooXOu2Mxm+8dDY3sYOA44trZ5XWaWX8tr17YEVLO2akshlz4xndJyR7u0FF79+RH0ywzbIRYRqbOIelJm1tYfrLAR2AKsCjl2iJm9bWYHRyvIEFnAlhqGvq8HOptZCz+OXniXBvsDq/y5UrvM7J1GiKtZW7xxB+f+9QvW5++mdYtkHr90rBKUiERFJGv3pQOfA0OB2XhJakhIlXnAUXiX5r6NQoyhWgM1zc0qDqlT4pxbQz0Gbzjn9jtgw+9pqTdVzfJNO/nBY1+zrbCEDq1TeeKysYzq2SHosEQkQUTSk7odL0Fd5pw7GJgSetA5VwR8Chzf8PC+owhoWcOxtJA60gQKikq56qmZbPUT1AtXH6oEJSJRFen28e865/Y3lHsN1e4PRUkucGANq110x7sU2ODRe2b2SZhiXb8KUVpewU1TZrN6axEpScYTl41lcNf2QYclIgkmkp5UD2qfELuLxrk0NgMv5nGhhWaWBowEZjbCa0o1FRWOX700lw8WbQLg12cOVQ9KRBpFJD2pnUCXWur0wbtXFW0vArcBNwCfhZT/CO9e1HPReBHn3PjqZbon5dlTVs5Pn/2WjxZ7Ceq64/pzyaG9Ao5KRBJVJElqBnC6mbVzzu2sftDMuuEtmfRmXU9oZpcAlZ90mUALM7vDf77GOfcMgHNunpn9GbjGzKYCb7N3xYlPCTORV6LHOcdtU+dXJaiLDunJjd8bGHBUIpLIIklS/we8A7xtZleHHjCzIXgTa9PwVoCoqyuBY6qV3es/foq3A3ClG4DVwNXAaXg9tkeAu/azJFK96J5UeI9+tJyXv10HwNVH9+XWUwZjpu02RKTxRLIs0rtmdjdwNzAfKAUwsy1AB7xh379yzn1Rj3OOr0fdcrw1+ybVOWhpsGe/WsOk95cCcOaILCUoEWkSEa044Zy7x8w+w7vMdijQCXB4l98edM59FL0Qm57uSe1ry649/PathQD0y2zDfeccpAQlIk0i4mWRnHMfAx9HMRaJQbv2lPGjp2dSXFqBGTx31aG0bRnzGzOLSILQp00YuiflKSuv4LoXZjFrbT5m8NAFI+manlb7N4qIRElMrIIusenpL9dUjeS758yhnDWye8ARiUhzE8nafRV495/2xwE7gEXAVODRGhaFjUm6JwWvzV7PfW8vAuCUYV255LDewQYkIs1SJD2paXgrThjethlr8LbOWOM/N7xFZtcDo4A/AP+tbXt3iR0zV2/jl1PmUlbhGNKtPfefOzzokESkmYokSV2I16P4N96Oun2dc4c55/oC/fzy9sD3gAOAJ4CDgZujE7I0pg0Fu7nm+VmUlFfQPaMVz145jvRWqUGHJSLNlDlX25W7at9g9jRwoHNuzH7qzAQWOOcuNbNkYAFQ5pwb1qBom0gNAyeOTE9PT87Pr21vxPjlnOO0hz9n4YYdtG6RzJSfHMbQrGZzhVNEGkFGRgYFBQUFtW2HVJNIelInAe/XUud94GSomnw7DW89P4lhz09fy8INOwC4/9zhSlAiErhIhqC3w7uctz/pfr1K2yJ4ncA0x4ETM1Zv4/ZX5gMwplcHzhyRFXBEIiKR9aQWAxeYWdhPMTPrAVyAN7KvUjawNYLXkiYwJyefK/41A4DenVrzlx8cHHBEIiKeSHpSk/AWfP3WzB4B/gvk4Q2SOBK4FsgAHgAwsxTgBPbdWkNiyN1vLGDnnjLap6Xwz0vH0qW9JuyKSGyIZIHZ5/xe1O+Ae6odNqAMuN05V7m3UwZwF/B1QwKVxvH6nFxmrfUGg9x95lD6d2l2C2uISAyLdIHZP5rZS8AP8HbETcebvDsLeN45tzKk7hbg71GItck0l2WRvli+hVte9jZZPrRvR87QfSgRiTENWWB2FfDbKMYiTShvRzE/fvYbikrKaZeWwkMXjCI1WatkiUhs0QKzYST66L6KCsftr8xnZ3EZ6a1Smfzjw7RwrIjEpFqTlJkd7f9zunOuOOR5rZxz0yKOTBrN7a/O44NFeQD85syhDOrarpbvEBEJRl16Up/gLRg7BFga8rwukiOKShrNK7OHD75rAAAQCklEQVTW8cL0HAAuHJfNWSN1H0pEYlddktQ9eElpS7XnEmdWbynkVy/PA6Bbehp3nHagdtgVkZhWa5Jyzt29v+cSP/73nUWUlFWQmmw8d9UhtNEOuyIS4zScq5n4YsUW3l3g3YeaNGEkfTMTbkS9iCSgOv0p7Q+WWO2cW1vH+iOAEc65pxsSXFASbZ5UeYXjt296q1Qd3DODM4Z3CzgiEZG6qWtP6mPgstACM/uVmdW0Ht/3gX81IC6JotDVze86Y6juQ4lI3KjrTYlwn2ppeEseJZxEmif1wcI87nrNW9387FHdGZmdkD8yEUlQunOewGas3sZVT88EoGv7NG49ZXDAEYmI1I8GTiSoopIyfvrstwC0bZnCE5dpdXMRiT9KUgnIOcclj09ny649ADwwYQQHZtW2T6WISOxRkkpAz329lm/WbAfgnFHdOX7IAQFHJCISmfrck9IqE3FgR3Epk95bAsAhfToyacIIjeYTkbhVnyR1t5ndXb3QzMqjF4401D8+Xcn2olJat0jmkYtGKUGJSFyrT5Kq76dd3Pa84nUy7+6Scl6Y7s23vvLIPnRpp4ESIhLf6pSknHO6dxUH7n59AVsLS0hJMi4Ymx10OCIiDaZ5UmHE42TenG1FTP7G24LjlycNokeH1gFHJCLScOohJYj7/7MY5yCjdSqXH9En6HBERKJCSSoBvDV3A2/N3QDAzScNpkWKfqwikhj0aRbntheWcMvUuQD07dyG88f0CDgiEZHoUZKKc6/MWs/O4jIAnrhsLKnJ+pGKSOLQJ1ocK69w/HuGN+T8B4f0pHfnNgFHJCISXUpScey+txexNG8XABPH9gw4GhGR6FOSilOfLdvM45+vAmDi2GwO6hGzo+NFRCKmJBWHikrKuO2VeQAM6NKW208bEnBEIiKNQ0kqDj371Rpytu2mRXISf714NO3SUoMOSUSkUShJxZkNBbv5+6crAbhgbDb9u8T8koIiIhFTkooz97+zmK2FJbRPS+Hqo/sGHY6ISKPS2n1hxOoq6HNy8nltdi7grc+X3VHr84lIYlNPKk6UlVfwq5e9lSXatEjm5GHdAo5IRKTxqScVRiyugv7kF6tZvHEnAA9eMJLMdi2DCkVEpMmoJxUHnHNVmxmeNTKLE4d2DTgiEZGmoSQVBx7+cDkrNhcCaBsOEWlWlKRi3LK8nTz80TIAzjm4OyO0soSINCNKUjGsqKSMO16dT3mFo3PbFvzmzKGYWdBhiYg0GSWpGPbHd5fw9aptAFxzbH+tLCEizY6SVIzaUVzK5Bk5AFx5ZB8uPbx3sAGJiARASSpG3fvGQgpLymmVmsx1xw3QZT4RaZaUpGLQrLXbmfLNOgB+dHRf0lvrMp+INE9KUjHGOcetU71tODLbteT64wcEHJGISHCUpGLI7pJyzvnrF1UrS9z0vYEkJ+kyn4g0X0pSMeSd+RuYtTYfgAvGZDNxnLaEF5HmTUkqRoQufXRk/878/rzhAUckIhK8hE5SZvZzM5tuZsU1bL8RM26aMocZq7cDMHFcdsDRiIjEhkRfBX0DcD8wFjgs4FhqNHlGDlO/XQ/A0QMzOe0gbcMhIgIJnqScc1MBzCxmb+58uCiP21/1RvON6pnBYz8crTlRIiK+mLjcZ2a3mtkUM1tpZs7MVu+nbpKZ3Whmi/3LeDlmNsnM2jRhyFGRs62Inz//LaXljgFd2vL0FeNomZIcdFgiIjEjJpIUcB9wHLAC2F5L3QeBB4CFwLXAFOA64A0zi5X21GpHcSm/mDyb4tIK0lKT+OvFo7U2n4hINbFyua+fc24lgJnNB9qGq2RmQ/ES01Tn3Lkh5auAh4GJwPONH27DfLF8C798aS7r83cD8OCEkfTvErbJIiLNWkz0PCoTVB1cCBjwULXyx4Ai4OJoxhVtu0vK+dO7S7jon1+zPn83KUnGAxNGcIoGSoiIhBUrPam6GgtUANNDC51zxWY22z9excxS8NqYAiSZWRpQ4ZwrqX5iM8uv5bUj2m1wyUZv08LVWwpZtmkXJWUVAHRLT+PRi0YxulfHSE4rItIsxFuSygK2OOf2hDm2HjjczFqEJKE7gF+H1NkNfAqMb9QoQ5SWV/DW3A1Vz83gvIN7cPeZQ2nTMt7efhGRphVvn5KtgXAJCqA4pE4JgHPubuDuupzYOZexv+N+T6vevanendtw/uge9O7chuyOrTm0b0e6tEur72lERJqleEtSRUCXGo6lhdSJGW1bpvDH80cEHYaISFyKtySVCxxoZi3DXPLrjncp8Dv3m+qrhiWUNPxORKSJxcTovnqYgRfzuNBCf0DESGBmEEGJiEjjiLee1IvAbcANwGch5T/Cuxf1XDRexDk3vnpZpPekREQkcjGRpMzsEqCX/zQTaGFmd/jP1zjnngFwzs0zsz8D15jZVOBtYAjeihOfEgcTeUVEpO5iIkkBVwLHVCu713/8FHgmpPwGYDVwNXAasAV4BLjLOVcRjWB0T0pEJDbERJIKd3ltP3XLgUn+l4iIJDBzzgUdQ1wwswrA0tN1W0pEpK4KCgoAnHMuooF6SlJ1ZGZleCMLd0Tw7ZWZrSB6EcU8tbl5UJubh4a0uT3ecnQRXblTkmoClesC1raqRSJRm5sHtbl5CLLN8TZPSkREmhElKRERiVlKUiIiErOUpEREJGYpSYmISMxSkhIRkZilJCUiIjFL86RERCRmqSclIiIxS0lKRERilpKUiIjELCUpERGJWUpSjcTMkszsRjNbbGbFZpZjZpPMrE3QsdWVmQ00s3vM7Csz22xmO81stpndHq4dZjbIzF41s+1mVmhmn5nZcTWcO93MHjGz9f77s8DMfmpm1vgtqx8za21mq8zMmdmjYY4nRLvNrKOZ/cnMlvuxbTazj83sqGr1DjGzD/zfhx1m9h8zG1nDObPM7Gn/XLvNbKaZnd80Ldo/M2trZreZ2Ty/LVvM7Aszu6z6zyPe2mxmt5rZFDNb6f/erq6lfqO0z8xa+p8hq8xsj5mtMLM7zCy1zo1xzumrEb6A/wMcMBX4EfAAUAp8BCQFHV8d23A/sBN4DrgW+Anwot+uOUCrkLr9gK1AHnAr8DNglt/mE6qdtwUw3T/2gP/+TPXPe3fQ7Q7zPvzJfx8c8Gi1YwnRbqAXsArY7P/crwBuBP4FTAypdyhQDKzwj9/o/3sncFC1c3YEVgK7gHvwdtP+xG/v5QG3Nwn4DCgHnvBjuwH42o/v9/HcZv/1tgLvA9uA1fup22jtA171jz0OXOU/OuDJOrclyF+URP0ChgIVwMvVyq/1f0AXBR1jHdsxBkgPU/5bvx3XhJRN9v/DjwwpawusAZbgT3fwy3/mf/+11c77MlAC9Aq67SExHQyUAb8gfJJKiHb7H9g5QLda6k3H21Ote0hZd7/svWp1/+C394yQsmT/HFuBtgG29zA/tgerlbfwP4Tz47nNQN+Qf89n/0mqUdoHnOrXnVTtHJP88sPr1JagfkkS+SvkQ/yoauVpQCHwdtAxNrB9B/nt+5v/vA3eX2Ifhql7p193XEjZ5/77kFat7lF+3ZuDbqMfTzLwDfAm0JtqSSpR2g0cTUjyBFKB1mHq9ffrPR7m2ON4f5h1DSlbBywPU/cS/zwTAvzZnuTH8Mswx6YD6xOlzftLUo3ZPuBZvyy7Wt1sv/wvdYlf96Qax1i8H+700ELnXDEw2z8ez3r4j3n+43CgJfBlmLpf+Y9jwbtXh9c7meW/H6Gm471vsfL+3AgMBq6p4XiitPtU/3Gtmb0B7AYKzWypmV0cUq8yvpraa8BoADPrhvfX+Fc11A09XxCmA/nAzWZ2vpn19O8t/i9eG+726yVSm8NpzPaNxUv2OaEV/ee51PG9UJJqHFnAFufcnjDH1gOdzaxFE8cUFWaWDNyFdwnseb84y39cH+ZbKsu6+48dgFbh6vrv19aQuoExsz7Ab4B7nHOra6iWKO0e5D8+hnfP4VLgSrxLkM+Y2eX+8fq0tz51m5xzbjtwJt79msl4l2cXAz8HznXOPeZXTZg216Ax25dVQ93K+nV6LyLac15q1RoIl6DAuzxUWaekacKJqofwbrTe5pxb4pe19h/Dtbm4Wp391a2s37qGY03pr3gDCR7YT51EaXc7/3EncKxzrgTAzF7Buz9zn5k9RfTaW71uUHbhXQp7HfgCL0H/HHjezM5yzr1P4rW5usZsX22fg3V6L5SkGkcR0KWGY2khdeKKmd2Ld+nrH865/w05VNmWlmG+rXp791e3sn6g741/ietE4GjnXOl+qiZKu3f7jy9UJijwehtm9jrwQ7zeVrTaG/j/ATM7CC8x3eic+1tI+Qt4iesxM+tHArW5Bo3ZvqIa6lbWr9N7oct9jSMX75JeuB9Qd7xLgXHVizKzu4E78IYk/6Ta4Vz/MVz3vbKsstu/He9D8Tt1/ferEzVfImh0fgwPAG8DG82sv5n1xxuiDZDul2WQOO1e5z9uDHNsg//Ygfq1tz51g3Aj3gfllNBC51wR8Bbez7s3idXmcBqzfbk11K2sX6f3QkmqcczAe2/HhRaaWRowEpgZRFCRMrNfA78Gngaucv4QnRDz8Lr1h4X59kP9x5kAzrkK4FtgVJgkPg7vfQvy/WkFZAKnActCvj7xj1/sP7+KxGl35QCfHmGOVZZtwvu9hprb6/BGQ+Kc24D3IXRoDXUh2J9z5YdncphjKSGPidTmcBqzfTOA7maWHVrRf55FXd+Lph4O2Ry+8IZo72+e1MVBx1iPttzlx/w0+5mEjPcXaTkwIqSscr7QUvadL/Rzap4vVAr0CbC9qcB5Yb5+6sf8jv98YKK0G6+XtAOvRxU6z6Ub3n2bpSFlM/y6WSFlWX7ZB9XO+0dqnlOzHWgX4M/5QcIM+wcqe8jbgJREaDO1z5NqlPbh/aG3v3lSR9Yp/qB+SRL9C3iEvStOXOX/YErx/iKPlxUnKj9U1+Ddl7i42tf3Qur29/9j5wG3sHflhTLgpGrnbYH3V1Sp/75cxd6VF+4Nut01vBe9CT+ZNyHajbdygPM/0H7ht2UN3uCeE0PqHY7Xe1yBt0LDDf6/dxGSqP26nYDVeAMyfuO/xsf+61wZcHt74Y2orACewbuEfRveYBkH/Cye24w3b+kO/ysPL4FUPr+kWt1Gax/whn/sn3gjRv/pP3+mzm0J8hclkb/w/rq4CW/VgT143eQHCHCWfQRteNL/harp65Nq9YcAr+HNPynCm7x6Qg3nzgAexfurdQ+wEG9QhjV2uyJ8L3oTJkklUruBc/DmuxT6H0LvAUeEqXcY8KH/IbYTeBc4uIZzdsdLAlvwRnR9C1wQdFv92PoBT+H1IEvxeg7TgHPivc3sXaqo1v+3jdk+vPt+v8VLbHvwRoveCaTWtS3amVdERGKWBk6IiEjMUpISEZGYpSQlIiIxS0lKRERilpKUiIjELCUpERGJWUpSIiISs5SkREQkZilJiYhIzFKSEhGRmPX/1Mm54dQIEjQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "eigs = np.linalg.eigvals(A)\n",
    "# print(eigs)\n",
    "plt.semilogy(np.unique(eigs))\n",
    "plt.ylabel(\"Eigenvalues\", fontsize=20)\n",
    "plt.xticks(fontsize=18)\n",
    "_ = plt.yticks(fontsize=18)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.059164723442146e-09\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "import scipy.sparse.linalg as spsplin\n",
    "\n",
    "res_hist = []\n",
    "def callback(x):\n",
    "    res_hist.append(np.linalg.norm(A @ x - b))\n",
    "    \n",
    "x, info = spsplin.cg(A, b, tol=1e-10, callback=callback)\n",
    "print(np.linalg.norm(A @ x - b))\n",
    "print(info)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "226\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAEMCAYAAAD9OXA9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8VFX6x/HPk0ICBBIIXRCkhyYgTVAE2yqICuIiFrD3squu7rq6q+5a1p/rKuKqWLFhD4iUVREUpAvSIkV6LwmEQAKE5Pz+mIGN2QQyYSZ3Jvm+X695TThz5s5zLzPzzDn3nHPNOYeIiEioRHkdgIiIlG9KNCIiElJKNCIiElJKNCIiElJKNCIiElJKNCIiElJKNCIiElJKNCIiElJhnWjMLM7M3jKzTWaWaWZTzayt13GJiEjJhXWiAWKANUAPoCYwHhjraUQiIhIQi6QlaMysEnAAqO2cS/c6HhEROb6gtmjMrKGZvWBmM8xsn5k5M+tTTN0EMxthZlvNLMfM5pvZxcd5iZ7ADiUZEZHIEeyus+bAUGAfMOU4dVOBq4CHgf5AGpBqZv2KqmxmScCrwENBi1ZEREIuqF1nZhblnMv3/30pvmTS1zk3rVC9fsAEYJBzLtVfZsB0INk5l1KofjwwGfjROXdf0AIWEZGQiwnmxo4kmRIYCGQC4wo815nZaGCUmbVxzqUBmFk08CGwGbi/uA2a2Z7jvGYi4IC9JYxRRESgOpDvnCt1vghqoglAOyCtiMS0uODj/r9fAyoDl7sTb35ZYmJi4gluQ0SkwsjMzIQTPM3iVaJJBlYWUZ5R4HHMrDFwHb6RZrt9vWsAXOicm17wic65pGO9oJntSUxMTNyz53gNHxEROSIpKYnMzMwT6gnyKtGArxvrmI8559YDdox6IiIS5ryasJmOv9VSSE3/fUYRj4mISATyKtEsA1LMrPDrt/ffLy3jeEREJES8SjSpQBIwoFD5MGDFkRFnIiIS+YJ+jsbMBvv/7Oq/P8vMagH7nXOT/GUTganAG2aWDKwFhgNnAJcEOyYREfFO0Nc6M7PiNrjeOdekQL3qwJPAYHytmzTgcedcSBbN1KgzEZHA+UedZR5vZO+xBL1F45wr0Sgx59xe4E7/LazNWLWLtbv2cX7betStHu91OCIiEcXL4c0R451Z6/gqbTuPjFtG3epxNE6uyoBTG3D5aQ2Jj432OjwRkbAW7tejCQu1qsVRPd6Xk7fvPcjctRk8MnYpl4z8gZXbszyOTkQkvEXU9WhOxImeo8nNy+fH9bvZvDuHmavTGfvTZvLyHWZwXkpdHrigFc3rVAty1CIi3grGORolmlL6cX0Gv/9oERsysgGIjjL6tqrN0G4nc05K3aC8hoiI15RoAhCKUWeH8/KZtHQb//xqBevSs4+WX9C2Ho9d0lYDB0Qk4inRBCCUw5sPHc5n8rJtfDh3AzNX+y7+WS0+hnNT6tKnVW0uPrUBBRYEFRGJGEo0ASiLeTTOOcb9tIXHxi9jd3bu0fLz29Tlz/1TaJxcNWSvLSISCko0ASjLCZvp+w7y0fyNzFmTwXcrdwIQZXBmi9r0a1+P89rUo2bVSiGPQ0TkRCnRBMCLlQGcc6Qu3Mw/v1rJ5j05R8tjoozb+zTjzrNbUClGI8xFJHwp0QTAyyVocvPy+SZtOxOXbuPbn7ez/1AeAKc2TOSlqzrTsEaVMo9JRKQklGgCEC5rneUcyuP5b1YyavoanIPEyrE8P6QjfVvX8TQuEZGiBCPRqN+mjFWuFM2f+qXwzvXdqFm1Epk5uVz39jwe/WIZ2YcOex2eiEjQKdF45MwWtZlw9xmc1rgGAG/PXMeFL0xnzpp0jyMTEQkuJRoP1U+szIc39+APv2lFbLSxPj2bIaNm88Cni8jYf8jr8EREgkLnaMLEim1ZPPDpIhZtygQgqUosfx3QhoGdGnocmYhUZDpHU460qleNz2/vxd8uaUu1+Bj2ZOfy+48W8dxXK6goPwZEpHxSogkj0VHGNac34dv7+tC7ZW0ARnz7Cze9M5892epKE5HIpEQThmpXi+ON4V24omsjAL75eQf9R8xgwYbdHkcmIhI4JZowFRsdxdOXdeDZy0+lcmw0m/fk8NtXZjHq+9Xk5asrTUQihwYDRIBV27O444MFrNy+D4BWdatxW59mXNCuni4lLSIhpZUBAhDJiQZ8Kwo8/uUyxszdeLTspKTKjBjakdMa1/QwMhEpz5RoAhDpieaIpZszeWnqL3ydtp3D/ktJdzgpkeE9mzCw00m67o2IBJUSTQDKS6I54pcdWdz5wUKWb8s6WtajaU3+fml7mtdJ8DAyESlPlGgCUN4SDfhWhZ69Jp0P5mxg0tJtAMRGG7f1ac7tfZrp/I2InLAKMWHTzOqa2Vdmlm1mi83sNK9jChex0VGc2aI2L199Gm9d15VGNSuTm+cYMWUV/V6YzszVu7wOUUQk/BMN8DLwC5AMjAQ+M7NYb0MKP31b1eGr353F7X2aERNlrNm1nytfm8N9H2vdNBHxVlh3nZlZNSAdaOSc2+4vWwcMd859F+C2yl3XWXGWb9vLQ58vYcEG377WqBLLPee04PIujagaF+NxdCISScKu68zMGprZC2Y2w8z2mZkzsz7F1E0wsxFmttXMcsxsvpldXKhaC2DPkSTjtwRoE8y4y5vW9arz6a09+ful7agWH8Pu7FweHZ9Gr398y0fzNpCvCZ8iUoaC3XXWHBgK7AOmHKduKnAV8DDQH0gDUs2sX4E6VYG9hZ63F9CwquOIijKu7tGYKfeexdBujYiLiWJPdi4PfraEIaNmsXJ71vE3IiISBEHtOjOzKOdcvv/vS/Elk77OuWmF6vUDJgCDnHOp/jIDpgPJzrkUf1lnYJJzrm6B544HJjrnXi60zeP1iSUmJiZSEbrOirIz6yBPTEhj7E9bAIiJMm7q3ZS7z25B5UoanSYiRQu7rrMjSaYEBgKZwLgCz3XAaKC1mR3pGlsF1DCzugWe2w5f60cCULtaHM9f0Yn3buhOk+QqHM53vDxtNec//x3frdzpdXgiUo55NeqsHZBWRGJaXOBxnHNZwJfAI2YWb2Y3AgbMLLxB51zSsW74EluFd0aLWkz+XW/uPqcFsdHGxowchr85l999uJD0fQe9Dk9EyiGvEk0ykFFEeUaBx4+4DWgF7AbuAS5zzuWGNrzyLT42mnvPa8mke3rT7RTfOmljf9pCn2encd/Hi/hlh87fiEjweDmP5lgnh44+5pzb7pw7zzlX2TnX3jn3YxnEViE0r5PAhzf14OlB7akeH0PWgcN8tmATF4/8gf8s2+Z1eCJSTniVaNL5davliCPLEBfV2pEQiIoyruh2MlPv78NjF7elQWI82YfyuOXdH3lxyipdRlpETphXiWYZkGJmhV+/vf9+aRnHU+ElJ8QxvGcTxt15Bl0a1wDgn1+v5M4PFpJ96LDH0YlIJPMq0aQCScCAQuXDgBXOOY0q80jtanG8f1N3hnTxXUZ6wpKt9B8xg9lr0j2OTEQiVdDXIzGzwf4/u/rvzzKzWsB+59wkf9lEYCrwhpklA2uB4cAZwCXBjkkCExcTzdOXtSelfjX+PuFn1u7azxWjZnNl95P544WtqR6vpeZEpOSCvtaZmRW3wfXOuSYF6lUHngQG42vdpAGPO+fGBjWg/75ehVnrLJiWb9vLg58tYdFG33GrWz2Opy/rQN9WdTyOTETKgq5HEwAlmtLLy3e8PXMdz/5nBTm5eZjBXWe34I6+zYiL0aoCIuWZEk0AlGhO3MaMbO75cOHRVaFPSqrM45e05ZyUusd5pohEqrBbgkbKt0Y1qzDm5h7ceMYpxEQZm/fkcMPo+dz38SIyczSHVkSKphaNlMq6Xfv50+dLmOUfjVY/MZ6nL+vAWS1rexyZiASTus4CoEQTfPn5jvfmrOepicvJyc0DoE+r2txzTgs6nVzD4+hEJBjUdSaeiooyhp3ehMm/O5NuTXyLOkxbsZPfvjqLLxZt8Tg6EQkXSjRywhonV+XDm3sw6prTaF4ngdw8x91jFjLup81ehyYiYUCJRoIiKso4v209Pru1J6f5l7C5/5NFTFuxw+PIRMRrSjQSVIlVYnlzeFda+Fs21789j5HfrmL/Qa2XJlJRaTCAhMS2zAPcMHoey7bsBSCpSizX9mzC8NObUKNqJY+jE5GS0qizACjRlL2cQ3n8Y/JyxszdwMHDvoupVqkUzbU9m3D3OS2Ij9WqAiLhTokmAEo03tmZdZA3f1jLe7PWk+XvQmteJ4FHLmpD7xa1MDOPIxSR4ijRBECJxnt7D+Ty2vdreHnaag7n+95357epyzODO5BURd1pIuFI82gkolSPj+W+81sx9o5enN7Ud4HVr9K20++F6cxfp4uqipRXatGIJ5xzfLZgM4+MXUpObh7RUcY1PRpzR9/m1K4W53V4IuKnrrMAKNGEp1927OPODxawfFsWANFRxgXt6vHkpe1JrKILrIl4TV1nEvGa10lg3J29+HO/FBIrx5KX75iweCuDXv6BDenZXocnIkGgFo2EjQO5eXy2YBOPfZHGobx8kqtWYtSwLkdXGhCRsqcWjZQr8bHRXNW9Me/d2J2kKrGk7z/E0FGz+XDuBirKDyKR8kiJRsJOt1Nqknp7L5rWrsqhvHz++PkSbnpnPmt27vM6NBEpBXWdSdjaeyCXBz5ZzORl2wAwg1MbJtG7RS1u79tcKwuIlAGNOguAEk1kcs4xfvFWnpm8nE27c46WX9iuHiOv7Ex0lFYVEAklJZoAKNFEtsN5+Xy3ciff/LydMXM3AjCo80n847IOxEarB1gkVIKRaGKCGZBIqMRER3FOSl3OSalL9fhYXv1+DZ8v2MzqHfu47LSG/LZLI3WliYSpsP4paGZxZvaWmW0ys0wzm2pmbb2OS7z1xwtbc/fZzQFYtCmTv4xbxgXPf8/ctVrGRiQchXWiwdfiWgP0AGoC44GxnkYknjMz7j2/Fe/d0J3LOjckJspYl57N0Ndm89r3a8jPrxjdwSKRIqLO0ZhZJeAAUNs5lx7gc3WOppxatT2Lu8YsPLqMTZfGNRh5ZWfqJcZ7HJlI5CvTCZtm1tDMXjCzGWa2z8ycmfUppm6CmY0ws61mlmNm883s4tIGWUBPYEegSUbKtxZ1q5F6ey+GdGkEwPz1uxkyahabdmsJG5FwEEjXWXNgKLAPmHKcuqnAVcDDQH8gDUg1s36lCRLAzJKAV4GHSrsNKb8qV4rmH4M78MbwLsTHRrE+PZtB/57J4k1qwYp4rcRdZ2YW5ZzL9/99Kb5k0tc5N61QvX7ABGCQcy7VX2bAdCDZOZdSoO4k4MxiXvJx59wz/nrxwGTgR+fcfSXfvV/Fpa6zCmL2mnRuGj2frIOHiY4yBnU6iT/8phV1qqsrTSRQZdp1diTJlMBAIBMYV+C5DhgNtDazNgXKL3TOJRRzO5JkooEPgc3A/cW9qJntOdYNSCzpvkpk69E0mc9v70mT5Crk5Ts++XET5zz3HR/M2aCBAiIeCMWos3ZAWhGJaXGBxwPxGlAZuNZF0sgF8VSLutX4z+978/dL21GjSixZBw7zUOoSrhg1m192aM00kbIUikSTDBQ1oSGjwOMlYmaNgeuA3sBu/yCEfWb2P91tzrmkY93wtbKkAomLiebqHo355t6zGNTpJADmrsug3wvTmbx0m8fRiVQcoZpHc6yWR4lbJc659c45c85VLtStNj0IMUoFkZwQx3NDOvLO9d1oVLMyh/Lyuffjn1i5Pcvr0EQqhFAkmnSKbrXU9N9r+rZ4onfL2oy9vRcNEuPJPpTHJSN/4G9fprEt84DXoYmUa6FINMuAFDMrvO32/vulIXhNkRJJTojjlWtOo0aVWHJy83hjxlrOfOZbfvvKLD6apwusiYRCKBJNKpAEDChUPgxY4ZxLC8FripRYh4ZJfP9AXx68oDW1EiqRm+eYuy6DBz9bwi3v/siWPTnH34iIlFhAS9CY2WD/n12BB4BH8bVg9jvnJvnrGL4JnR38ddYCw/Elmkucc+ODFXwgNI9GipJzKI9pK3bwxaItTPIPEIiLieLyLg258YymNKlV1eMIRbxV5tejMbPiKq93zjUpUK868CQwGF/rJg3fBEzPFsRUopFjcc4x9qfNPDlxOTuzDgKQEBfDezd2p2OjUn++RCKeLnwWACUaKYn9Bw/z4byNvPLdanZmHSSxcixvDO9ClyY1j/9kkXKoTFcGEKkIqsbFcMMZp/DxLadTu1ocmTm5XPnaHD6Yo4ECIqWlRCNShFNqVeXTW0+nWe2qHMrL56HUJdz5wUIyc3K9Dk0k4ijRiBSjcXJVPr+9F/3a1wNgwpKt9HthOqs00VMkIEo0IseQWDmWl67szFOD2hMfG8XmPTkMfW0OizbqXJ9ISWkwgEgJLd2cyVWvzznafdb55CSu6t6YSzo2ICZav9mkfNKoswAo0UgwLNmUye8//ulXK0C3rleNP/VLoXeLWvimkYmUH0o0AVCikWBxzjF7TQbvzVnPxCVbOfIRSqlfnVvPasqADg2IilLCkfJBiSYASjQSCos37eGpicuZtSb9aNmQLo14+rL2at1IuaBEEwAlGgmlRRv38OK3v/DNz9sBOL9NXa7u0ZjeLWt7HJnIiVGiCYASjYSac44/fb6ED+dtPFp299nN+f15LdW6kYilRBMAJRopC3n5jk9/3MiYuRv5yT8E+ozmtfhz/xRS6lf3ODqRwCnRBECJRspSbl4+D32+hE9+3ASAGVzRtRF/7t+GhLgYj6MTKTmtdSYSpmKjo3hmcAfeurYrzWpXxTkYM3cjA1/6gbQte70OT6RMqUUjEmK5efmMnrmOf0xeTm6ewwzOTalL31Z1GNT5JOJjo70OUaRY6joLgBKNeO3H9bt58LPFv5rsWbd6HE8P6kDf1nU8jEykeEo0AVCikXBwOC+f8Yu38O3ynfxn6TYO5eUTE2W8OLQTF7av73V4Iv9DiSYASjQSbjZmZHPTO/NZvi0LM3jwgtbc0ruphkJLWNFgAJEI1qhmFcbc1IOOjZJwDp6etJzfffQTOYfyvA5NJKiUaEQ8VKNqJT68uQeDOp8EwLifttDjqSk8M3k5h/PyPY5OJDiUaEQ8Fh8bzT8vP5WH+6cQG21k5uTy72mruf39BWQfOux1eCInTOdoRMLIzqyDvPLdat6YsRaAk5Iq83+DO9CzeS2PI5OKSudoRMqZ2tXieOSiNjw6oA2Von1X9Lzmzbm89cNaKsqPQil/lGhEwtC1vU5h8u/OpN1J1cnLdzw2Po2HUpew/6C60iTyKNGIhKmmtRP45Jae9O/gm18zZu5Gejw5hZHfriI/X60biRwRkWjM7HQzyzezP3odi0hZqlwpmpFDO3H/+S2Ji4ki6+Bhnv1qJde8OYcf1+/2OjyREgn7wQBmFgXMAhww1jn3dCm3o8EAEtF27z/Ec1+v5N3Z64+W/aZtXR6/pB11q8d7GJmUZxVlMMAtwA/Acq8DEfFSjaqV+Nul7Rh1zWm0PykRgP8s2865z33HmLkbNFhAwlaJEo2ZNTSzF8xshpntMzNnZn2KqZtgZiPMbKuZ5ZjZfDO7uDTBmVkycA/wWGmeL1Iend+2Hl/c2Yvnh3SkRpVYsg4c5k+fL+GuMQs1WEDCUklbNM2BocA+YMpx6qYCVwEPA/2BNCDVzPqVIr6ngOecc5mleK5IuWVmXNrpJL659ywu8g8W+HLxVq5+Yw57D+R6HJ3Ir5U00XzvnKvjnLsAeKu4Sv5kci5wo3PuDefct8BwfOdY/lmo7iR/66io2wNmdhrQGXi9JAGa2Z5j3YDEEu6rSMRITojjxaGdeOzitpjBwg17GDpqNuvT93sdmshRJUo0zrmSLro0EMgExhV4rgNGA63NrE2B8gudcwnF3J4BzgTaADvMbBdwBfCImb1ZwlhEKgQzY3jPJvzz8lOJMli2ZS/9R8xg/KItXocmAkCwL17eDkgrIjEtLvh4Cbf1OvBpgX8/B/wCPFtU5eONiFCrRsq7QZ0bUrd6PL/76Cd2Zh3krjELSV24mWGnN6b7KclUrqQreYo3gj3qLBnIKKI8o8DjJeKc2+ec23TkBmQDe51zRW1fRIBezWsx6Z4z6d2yNgDfLt/BtW/N47S/f807s9ZpZJp4IhTDm4/1Ti71u9w5d21p59CIVCS1EuJ4+9quvHL1aXRpXAMzyD6Ux1/GLeP6t+exM+ug1yFKBRPsRJNO0a2Wmv57tUZEykBUlHFBu3p8eltP5v/5XAZ28l3vZuqKnVzw/PdM+Xm7xxFKRRLsRLMMSPHP5i+ovf9+aZBfT0SOIzkhjn8N6ciIoZ2oFh9D+v5D3DB6Pk9N+pk8rZkmZSDYiSYVSAIGFCofBqxwzpV0IICIBNnFpzZg8u960/0UXwfDq9+tYdDLM5m2YofHkUl5V+K1zsxssP/PrsADwKP4WjD7nXOT/HUM34TODv46a/HNoxkGXOKcGx/M4AOhtc5EfA7n5fPExJ9564d1R8vOTanDkwPbU0drpkkhwVjrLJBEU1zF9c65JgXqVQeeBAbja92kAY8758aWNshgUKIR+bW5azP419crmbUmHYD6ifG8OLQTXZrUPM4zpSIp00QT6ZRoRP6Xc47xi7fy0OdL2OdfJ63bKTX5c78UTm1U6u8VKUeUaAKgRCNSvBXbsrj3459YtmXv0bLeLWtzR59mdG9a4ulvUg4p0QRAiUbk2JxzzPhlF09M+Jnl27IAiDL4w29ac0vvpkRFmccRiheUaAKgRCNSMnn5jik/b+eFKauOtnA6Nkri+SEdaVKrqsfRSVmrKBc+E5EyFB1lnN+2Hp/d1pOh3RoB8NPGPQx7cy679x/yODqJREo0IlKk+NhonhrUgfdv7E58bBQbMrIZ+tpsJi7ZqomeEhAlGhE5pl7Na/Hs5acCsHxbFre/v4A+z05lqiZ6Sgkp0YjIcV3UoQGf3no656bUxQw2ZuRw0+j5TF66zevQJAJoMICIBGT1zn3c8f4Clm/LIjrKGHFFJ/r7Lyct5Y8GA4hImWtWO4EPbupBm/rVyct33DVmAU9MSCPrQK7XoUmYUqIRkYDVrFqJD27qToeGieQ7eG36Ws7553eMXbhZF1eT/6GuMxEptQO5ebw+fQ0jp/7CgVzfFdy7NK7BbX2a0adVHaI1yTPiacJmAJRoREJn854cnpiQxsQl/x0c0KNpTUYN60L1+FgPI5MTpUQTACUakdD74ZddvPLdaqav2gVAy7oJ/OWitvRqnozvKiISaZRoAqBEI1J2PvtxEw98tvjoxM52J1Xnjj7NuaBdPSWcCKNEEwAlGpGytXDDbp6auJy56zKOlp3ZohZPDmxPo5pVPIxMAqFEEwAlGpGy55xjwYY9/HvqL0xZ7ltJoGqlaP7cvw1DuzVS6yYCKNEEQIlGxFtfLdvGn8cuZWfWQcB3+egXruhE1bgYjyOTY9GETRGJGOe3rcfXv+/NpR0bAPDNzzsY+tpsrQhdASjRiEiZSapSieev6MQTA9sRZbB4UyZXvzGHzGytKlCeKdGISJm7qntjRl7ZmegoY9mWvVzz5hwyc5RsyislGhHxRL/29Xl+SMejLZvBL89k6oodWsKmHFKiERHPDDi1Af/yJ5tVO/Zx3VvzuGLUbFZtz/I6NAkiJRoR8dQlHU9i3B1n0Kt5MgBz1mYw6N8zmbl6l8eRSbCEdaIxn0fMbIuZ7TWzaV7HJCLB175hIu/f2IN3ru9GnWpxZB08zLVvzWOaruJZLoT1PBozuxsYAAwHtgEdnXMLSrktzaMRiQCb9+RwzRtzWLNzP1EGTWsn0L99fa7vdQqJVbRAZ1kr1xM2zSwa2AT0dM6tDcL2lGhEIsSOvQe48vU5/LJj39GymlUr8fyQjvRuWdvDyCqeMpuwaWYNzewFM5thZvvMzJlZn2LqJpjZCDPbamY5ZjbfzC4uRWyNgHjgGjPbYWZpZnZ5KbYjIhGmTvV4vrzrDN66riu3ntWMhLgYMvYfYvhbcxn1/WqNTIswJT1H0xwYCuwDphynbipwFfAw0B9IA1LNrF+AsTUAkoB6wMnAzcCbZtY8wO2ISASKj42mb6s6/PHC1vzn973p2CgJ5+DJicu59+NF7MnWigKRokRdZ2YW5ZzL9/99Kb5k0tc5N61QvX7ABGCQcy7VX2bAdCDZOZdSoO4k4MxiXvJx4GtgAXCyc26j/zlfAF84514PZCf9z1XXmUgEO5Cbx/2fLOLLxVsBqJVQib8MaMuADvW1OGcIBaPrrESr2R1JMiUwEMgExhV4rjOz0cAoM2vjnEvzl194rA2ZWQJQ4qnCZna8DJJY0m2JSPiJj41mxBWd6HZKTZ6ZvIJd+w5x95iFpC7YxGMXt+PkZF16IFwFe3hzOyCtiMS0uMDjJeKc24cvYT1gZrFm1h04C5galEhFJOJERRnDTm/C1/f25tyUugBMXbGTs/85jfs/WcTkpVvJzSvp72IpK8FONMlARhHlGQUeD8RtwCnAbuBd4Hrn3OqiKjrnko51w9fSEpFyoH5iZV4bdhqvXN2ZBonxHM53fPrjJm59bwEDXpzBoo3qIg8noZiweayTPgENFXHO7XLOXeScS3DOtXTOfXaCsYlIOWFmXNCuPt/e34e/XdKWns2SiTJYvi2Lwa/M5J1Z67wOUfyCnWjSKbrVUtN/X1RrR0Sk1OJjo7nm9CZ8cFMPUm/vRfM6CeTmOf4ybhkjv13ldXhC8BPNMiDFzApvt73/fmmQX09E5KhTGyUx7o5enNfGd/7m2a9W8tdxS8k6oEsQeCnYiSYV39yXAYXKhwErjow4ExEJlapxMbx0ZWfOaV0HgNGz1nPGP6by+Pg0vlq2jYOH8zyOsOIp8RI0ZjbY/2dX4AHgUXwtmP3OuUn+OoZvQmcHf521+NYpGwZc4pwbH8zgA6F5NCIVS25ePq9PX8sLU1ZyIPe/I9Ga10ngX7/tSPuGmvFQEmW61pmZFVdxvXOuSYF61YEngcH4WjdpwOPOubGlDTIYlGhEKqZd+w7y4dwNTFm+g8WbMsnLd8THRvHBTT3ofHINr8MLe+V6Uc1gU6IRkSXdW0NNAAAN4ElEQVSbMrnt/R/ZtDuHGlViufucFvy2SyOqxpVo7nqFVGaLaoqIlAftGyby7g3dqVm1Eruzc3lsfBoXvPA9c9dqQGwoKdGISIVySq2qpN7ek6HdTqZSTBQbM3L47auzGDpqNlszc7wOr1xSohGRCqdxclWeGtSeL+86g1P9gwJmrUnn5nd+5ECuRqUFmxKNiFRYLetWY+wdvXjpys5EGSzZnMmt7/3IzqyDXodWrijRiEiFZmb071CfBy9oDcC0FTv5zfPfM3HJVo8jKz+UaEREgJt7N+VfQ06lerzvap63v7+Au8YsJGO/LrB2ojS8WUSkgG2ZB3jws8V8t3InANXiYri5d1Nu69OMmOiK99tc82gCoEQjIiXlnOOjeRt5YsLPZB08DEDPZsm8dGVnalSt5HF0ZUvzaEREQsDMuKLbyUx/sC/X9WoCwMzV6QwYOYNvl2/XIp0BUotGROQ4vly8hT98spgc/9DnanExjBrWhdObBXotx8ijFo2ISBm4qEMDPrutJ6c2SiLKIOvgYW4YPY+JS7ZSUX6snwi1aEREArBpdzZDXp3N5j2+VQQ6n5zEn/ql0LVJzeM8MzKpRSMiUsYa1qjCR7f0OHq9mwUb9nD5K7O4cfR8NqRnexxdeFKLRkSklGavSeepSctZtNH3vVItPoZHB7Tl/LZ1qRYf63F0waHhzQFQohGRUHDOMWnpNv4ybhm79vmWromNNrqfksx957ekU4Rf80ZdZyIiHjMz+rWvz4S7z+C8NnWJjTZy8xwzftnF8Dfnsj59v9chek4tGhGRINp38DAzVu3k4bFL2bXvEC3rJvD2dd1okFTZ69BKRS0aEZEwkxAXwwXt6vPSlZ2JiTJWbt/HgBdnMHtNuteheUaJRkQkBLo3Teat67qSWDmW9P2HuPr1Ofzp88VMXbGD/PyK0ZN0hLrORERCaEN6Nje/O5/l27KOlrWsm8CtZzVjwKkNiA3zhTo16iwASjQi4pXsQ4d5b/Z6vvl5B3PXZhwtT6lfnReHdqR5nWoeRndsSjQBUKIRkXCwbEsmL09bzcQlW8l3UDk2mjev7Rq266ZpMICISIRp2yCRkVd25uNbTqdBYjw5uXlc//Y8Zqza5XVoIaNEIyLigS5NavLJbT1pVLMyObl5DH9rLs99tYJNu8vfMjZhnWjMrLOZzTSzvWa22syGeR2TiEiwnJRUmY9vOZ2U+tXJy3eM+PYXzn72OyYv3eZ1aEEV1okGeAeYACQBlwMjzayFtyGJiARP/cTKfHrr6dzcuym1EipxKC+fOz9YwEtTfyH70GGvwwuKsB4MYGZZQAfn3Fr/v+cCf3POjS/FtjQYQETCWsb+Q1z1+hx+3roXgCbJVXjvxu40rFHFs5jKbDCAmTU0sxfMbIaZ7TMzZ2Z9iqmbYGYjzGyrmeWY2Xwzu7iU8b0EXGNmMWbWDWgIzCnltkREwlrNqpX48OYe3HjGKVSKiWJdeja/fWUWy7Zkeh3aCSlp11lzYCiwD5hynLqpwFXAw0B/IA1INbN+pYjvS+B64AAwE3jIObejFNsREYkIiZVjefiiNrx7fTeqVopmS+YBBv57JuN+2ux1aKVWoq4zM4tyzuX7/74UXzLp65ybVqheP3znVAY551L9ZQZMB5KdcykF6k4CzizmJR8HXgfWAtcCXwApwGTgMufc/7RqzOx4fWKJiYmJqOtMRCLFsi2Z3P7+AtanZxNl8MhFbRjStRFVKsWUWQxl1nV2JMmUwEAgExhX4LkOGA20NrM2BcovdM4lFHN7BmgG7HXOpTrn8pxzS4FpQO8SxiIiEtHaNkjkizvPoGOjJPIdPDY+je5PTGH8oi1ehxaQYI86awekFZGYFhd4vKRWAtXMbID5tAbOBpYUVdk5l3SsG74EKCISURIrx/LODd3866IZWQcPc9eYhdz5wQKWbo6Mr7VgJ5pkIKOI8owCj5eIcy4TGAL8DdgLfAWMdM5NPtEgRUQiSfX4WF4c2olZfzqHXs19X6NfLt7KgJEz+Ou4pew7GN7DoEMxj+ZYJ30CGkvtnPuPc66jc66ac+5k59yTJxibiEjEqpUQx+jruvHs5afSsm4CzsHoWes577nvmLoifMdJBTvRpFN0q6Wm/76o1o6IiJRQTHQUg09ryIS7z+SPF7YmLiaKrZkHuO6teTz6xbKwbN0EO9EsA1LMrPB22/vvlwb59UREKqTY6ChuPasZX/2+N12b1ADg7ZnrOPvZaYz6fjV7D+R6HOF/BTvRpOJbLmZAofJhwArnXFqQX09EpEJrnFyVMTf14A+/aUV8bBQ7sg7y5MTlXPj8dBZvCo/pHCVegsbMBvv/7Ao8ADyKrwWz3zk3yV/H8E3o7OCvsxYYji/RXFKapWOCRUvQiEh5t3lPDq99v4YP523gQG4+cTFRvH9jd7o0qXn8JxejTC98ZmbFVVzvnGtSoF514ElgML7WTRrwuHNubGmDDAYlGhGpKFZsy+Kmd+azISObGlViGXNzD1rXq16qbekKmwFQohGRimR9+n4G/nsmGfsPUSk6invPb8lNZzYlOsoC2o6usCkiIkVqnFyV0dd14+SaVTiUl8+IKavYmpnjSSxlt2COiIiUqfYNE5l0z5k8PWk5repV8+xyA+o6ExGRYqnrTEREwp4SjYiIhJQSjYiIhJQSjYiIhJQSjYiIhJQSjYiIhFRFGt6cD1hiYqLXoYiIRIzMzEwA55wrdcOkIiWaw/hacHtL8fQj2Skyrpsaejoev6bj8V86Fr9WHo5HdSDfOVfqCf4VJtGcCDPbA3AiE5bKEx2PX9Px+C8di1/T8fDRORoREQkpJRoREQkpJRoREQkpJRoREQkpJRoREQkpJRoREQkpJRoREQkpzaMREZGQUotGRERCSolGRERCSolGRERCSonmGMwswcxGmNlWM8sxs/lmdrHXcYWamfUxM1fMrXWhuueZ2Wz/8dlhZq+aWcSu62RmDc3sBTObYWb7/Pvcp5i6V5rZIjM7YGabzOxpM4svol5dMxttZrvMbL+ZTTezniHfmRNU0mNhZuuKea88XUTdSD0W55jZ22a2wsyy/f/fn5tZ+yLqlugzUZG+X5Roji0VuAp4GOgPpAGpZtbP06jKzoPA6YVu64486P/SmQhsBAYA9wMXAxPMLFLfW82BocA+YEpxlczsauB94AfgQuBJ4A7g7UL14v3bOQu4CxgIZAFTzKxT8MMPqhIdC7/v+d/3yksFK0T4sbgVOBn4F77/73v9/55nZj2OVArwM1Fxvl+cc7oVcQP6AQ4YWKDMgBnAz17HF+J97+Pf90uPU28usBCIKlB2nv+5Q7zej1Lue8F9udS/L30K1YkGtgLjCpXf5K/fvUDZ7f6yzgXK4oA1wCSv9/dEj4X/sXXA2BJsL5KPRZ0iypKA3cBnBcpK9JmoaN8vkfqrsywMxHcNiXFHCpzv3TAaaG1mbbwKLByY2UlAV+Bd51z+kXLn3NfAZuAyr2I7EQX35Rh6APXwvRcKeh/I5df7PhBY4pxbUOA1DgJjgPPMrNqJRRw6JTwWgYjkY7GjiLI9wCqgIQT8mahQ3y9KNMVrB6QV8WFbXODx8u5VMztsZplm9qWZnVbgsSP7v7SI5y2hfB+fIvfdOZcNrObX+96ucD2/xfhaRimhCNADZ/vP4xwysyVmdpuZWaE65epYmFltfr1PgXwmKtT3ixJN8ZKBjCLKMwo8Xl5lAs8DNwN9gT8AbYAfzKy7v86R/S/uGJXn4xPIvleE99GXwN34uoOGACuBfwPPFapXbo6FP4mOwvcd+qy/WO+LYpT60pwVxLGWTSi3Syo45xbi62c+YrqZfYHvl9oTwLkFqxe3mRCFF05Kuu/l+n3knLuzUFGqmb0P3G1mzzvn1hesfqxNBT+6kPk/fOetrnPO/VzoMb0vClGLpnjpFP2roqb/vqhfI+WWc24b8BW+8xPgOz5Q/DEqz8cnkH2vqO+j0fi+X7oVKCsXx8LMngDuA+5xzr1d4CG9L4qhRFO8ZUBKEUMSj4ybL6oftryL4r+/tJb574vqS25P+T4+Re67mVUBmvHrfV9WuJ5feyAPWB6KAMPAkc9NwXMQEX8szOxx4CHgAefciEIPB/KZqFDfL0o0xUvFN3xxQKHyYcAK51xa2YfkHTOrh2+Y5mwA59wmYD5wVcEPi5mdA5wEfO5FnGVkNrANuKZQ+VAgll/veyrQ3sw6Hikws0r+ut845/aGOFavDMOXZOYVKIvoY2FmfwUeAR5xzv1f4ccD/ExUrO8Xr8dXh+sN35j2b4FdwPX4Toq/je/DM8Dr+EK87+8Df8M3BLMPvslqa4FsoEuBemcDh4GPgHPwffFuwfdFHO31fpzA/g/23/6BrwX3V/+/LyxQZ7j/sZH+Y3QbsBf4pNC24vFNxFuD70T5efhOnucAp3m9ryd6LPAliQ/9//d98Q3hTfXXfaa8HAt8XWUOGI+v+7jgrVOBeiX6TFS07xfPAwjnG1Dd/0WyDTgALOA4kxjLww34I/ATsAffvJBt/i+TdkXUvQCY4z8+O4HXgBpe78MJ7r8r5rauUL2r8Q1bPYhvnsQzQOUitlcPeBdfv3s2vkl5Z3i9n8E4Fv4v2m/wTWA9hG+m/0xgeDHbi8hjAUwL4H1Ros9ERfp+0fVoREQkpHSORkREQkqJRkREQkqJRkREQkqJRkREQkqJRkREQkqJRkREQkqJRkREQkqJRkREQkqJRkREQur/AUg3eQ6bCRH/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy(res_hist)\n",
    "print(len(res_hist))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Регуляризация задачи наименьших квадратов\n",
    "\n",
    "$$ \\| Ax - b \\|_2^2 + \\lambda \\|x\\|_2^2 \\to \\min_x $$ \n",
    "\n",
    "- Решение записывается аналогично\n",
    "\n",
    "$$ (A^{\\top}A + \\lambda I)x = A^{\\top}b $$\n",
    "\n",
    "- Каким свойством теперь обладает матрица по сравнению с исходной постановкой без регуляризации?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "m, n = 1000, 10\n",
    "A = np.random.randn(m, n)\n",
    "x_true = np.random.randn(n)\n",
    "b = A @ x_true + 0.01 * np.random.randn(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.69988931780442\n"
     ]
    }
   ],
   "source": [
    "lam = 1\n",
    "x = spsplin.lsmr(A, b, damp=lam)[0]\n",
    "print(np.linalg.norm(A @ x - b)**2 + lam * np.linalg.norm(x)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.699889338125373\n"
     ]
    }
   ],
   "source": [
    "lin_op = spsplin.LinearOperator((n, n), \n",
    "        matvec=lambda x: lam * x + A.T @ (A @ x))\n",
    "x_cg, _ = spsplin.cg(lin_op, A.T @ b)\n",
    "print(np.linalg.norm(A @ x_cg - b)**2 + lam * np.linalg.norm(x_cg)**2)"
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}