{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d8b05cac",
   "metadata": {},
   "source": [
    "# Matrix inversion or solve\n",
    "> Short note on matrix inversion and why to prefer solve. Uses Hilbert matrices!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9261f5da",
   "metadata": {},
   "source": [
    "The short advice is always to solve for $x$ in linear equations $Ax=y$ with (for example) LU-decomposition, and don't do $A^{-1}y$. By why _exactly_?\n",
    "\n",
    "This post came to be after reading [Why Shouldn't I Invert That Matrix](http://gregorygundersen.com/blog/2020/12/09/matrix-inversion/) by Gregory Gundersen. He does an exellent job on describing matrix inversion, LU-decomposition and (ill) conditioned matrices. I'll repeat some terms here.\n",
    "\n",
    "The argument is that solving is _faster_ and _more precise_ than inversion:\n",
    "\n",
    "* LU decomposition requires less floating point operations than (plain) matrix inversion.\n",
    "* Floating point operations are less precise, so more lead to larger errors.\n",
    "\n",
    "\n",
    "In this post I use matrices related to polynomial interpolation, including Hilbert matrices. Hilbert matrices are _ill-conditioned_, so we expect that inversion will be bad.\n",
    "\n",
    "\n",
    "## Precision and floats\n",
    "I'll not go into the computational complexity, read the above post for a nice explanation. We are not solving symbolically but with 32 (or 64) bit floats. For example $0.1$ cannot be represented as a float: [see this article](https://www.exploringbinary.com/why-0-point-1-does-not-exist-in-floating-point/). The more operations we execute the more potential there is for errors. \n",
    "\n",
    "In mathematics both methods are precise if we solve symbolically. In mathematics there is the _condition number_ of matrices, a higher condition number means that small changes in $A$ lead to large changes in $x$. This seems exactly the intuition that we need to explore under what conditions we can expect more errors, and the difference between two methods is more pronounced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "d2bee436",
   "metadata": {},
   "outputs": [],
   "source": [
    "#hide\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy.random import default_rng\n",
    "from sklearn.datasets import load_diabetes\n",
    "import math\n",
    "from numpy.random import default_rng"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44c7e088",
   "metadata": {},
   "source": [
    "## Solve $Ax=y$\n",
    "We can solve $Ax=y$ with the LU-decomposition $A=LU$ and then solve $Lz=x$ and then $Ux = z$. Since $L$ and $U$ are triangular matrices this is straightforward.\n",
    "With triangular matrices you can start with an equation with one variable, solve it and continue.\n",
    "There'll always be an equation on a single variable (that you didn't solve yet).\n",
    "\n",
    "Finding the decomposition is the tricky part. The NumPy function `solve` does this with the method from LAPACK. We create an example from polynomial interpolation of the function $f(x)=x^2$ on the points $\\{-\\frac{1}{2},0,\\frac{1}{2}\\}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0d727056",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = lambda x: x**2\n",
    "x = np.array([-.5, 0., .5])\n",
    "A = np.vander(x, increasing=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b58c796f",
   "metadata": {},
   "source": [
    "For a small number of points this gives the same results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1f70fbd4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0., 0.])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.solve(A, f(x)) -  np.linalg.inv(A) @ f(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d2250b9",
   "metadata": {},
   "source": [
    "But for larger matrices it starts to differ:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c3a94bd5",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-1, 1, 40)\n",
    "A = np.vander(x, increasing=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "931aed2b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.40100652e-01,  1.11382317e-01,  1.00881039e-01,  4.94520867e-01,\n",
       "        3.78189174e-01, -6.48834908e+00, -3.71602241e+00,  2.78155158e+01,\n",
       "        9.44989685e+00, -7.18503682e+01, -2.39184816e-01,  1.33213768e+02,\n",
       "       -2.37054964e+01, -1.84869925e+02,  2.41140115e+01,  1.66567935e+02,\n",
       "        2.94400070e+00, -6.93285557e+01, -2.55111985e+01,  1.67702176e+00,\n",
       "        3.48501630e+01, -2.03380240e+01, -3.16823477e+01,  4.88047696e+01,\n",
       "        1.50645974e+01, -3.87264745e+01,  3.73249968e-01,  1.76218177e+01,\n",
       "       -2.31686235e+00, -5.51101135e+00, -1.84854387e+00,  4.96125557e-01,\n",
       "        2.44568978e+00,  1.08034577e+00, -1.51065340e+00, -8.77657212e-01,\n",
       "        4.59726355e-01,  3.29110196e-01, -5.81216797e-02, -4.90222286e-02])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.solve(A, f(x)) - np.linalg.inv(A) @ f(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ceda6ce",
   "metadata": {},
   "source": [
    "Also the speeds are different:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "2702d37a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "111 µs ± 13.5 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit np.linalg.inv(A) @ f(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e859f065",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "93 µs ± 11.8 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit np.linalg.solve(A, f(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac8916ca",
   "metadata": {},
   "source": [
    "## Hilbert\n",
    "The difference between the two method is due to the number of floating point operations and the amount of error introduced. We will show this with the Hilbert matrix as an example. The Hilbert matrix is _ill conditioned_, small changes in the matrix will lead to large changes in the computed output. \n",
    "\n",
    "Note that the Vandermonde matrix above is almost the Hilbert matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bc2e440",
   "metadata": {},
   "source": [
    "Define the Hilbert matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "7dd4e29c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.        , 0.5       , 0.33333333, 0.25      , 0.2       ],\n",
       "       [0.5       , 0.33333333, 0.25      , 0.2       , 0.16666667],\n",
       "       [0.33333333, 0.25      , 0.2       , 0.16666667, 0.14285714],\n",
       "       [0.25      , 0.2       , 0.16666667, 0.14285714, 0.125     ],\n",
       "       [0.2       , 0.16666667, 0.14285714, 0.125     , 0.11111111]])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def Hilbert(n):\n",
    "    return np.array([ [(1/(i+j-1)) for j in range (1, n+1)] for i in range (1, n+1)])\n",
    "Hilbert(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad1a5287",
   "metadata": {},
   "source": [
    "The inverse of the Hilbert matrix is known in closed form\n",
    "(see https://en.wikipedia.org/wiki/Hilbert_matrix):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "b718a047",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[     25,    -300,    1050,   -1400,     630],\n",
       "       [   -300,    4800,  -18900,   26880,  -12600],\n",
       "       [   1050,  -18900,   79380, -117600,   56700],\n",
       "       [  -1400,   26880, -117600,  179200,  -88200],\n",
       "       [    630,  -12600,   56700,  -88200,   44100]])"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def Hilbert_inv(n):\n",
    "    return np.array([ [\n",
    "        (-1)**(i+j) * (i+j-1) * math.comb(n+i-1, n-j) * math.comb(n+j-1, n-i) * math.comb(i+j-2, i-1)**2\n",
    "        for j in range (1, n+1)] for i in range (1, n+1)])\n",
    "Hilbert_inv(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81511035",
   "metadata": {},
   "source": [
    "The inverse of the inverse should be the same as the original:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "ab29785f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.all(np.isclose(Hilbert(10), np.linalg.inv(Hilbert_inv(10))))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55c06486",
   "metadata": {},
   "source": [
    "## Error inversion vs solve\n",
    "Define the error as the mean squared difference between solving through taking the inverse and solving throuh LU-decomposition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "id": "c6ef0217",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = default_rng()\n",
    "\n",
    "def error(n):\n",
    "    # Create the matrix A and solution y\n",
    "    A = Hilbert(n)\n",
    "    x = rng.uniform(size=(n,1))\n",
    "    y = A @ x\n",
    "    \n",
    "    # Solve Ax=y\n",
    "    x_inverse = np.linalg.inv(A)@y\n",
    "    x_solve = np.linalg.solve(A, y)\n",
    "    e_inverse = np.mean((x - x_inverse)**2)\n",
    "    e_solve = np.mean((x - x_solve)**2)\n",
    "    return np.stack((e_inverse, e_solve))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf091399",
   "metadata": {},
   "source": [
    "You can see that the difference is quite pronounced for large $n$ (chart is in log-scale)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "id": "e5852189",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA0yUlEQVR4nO3deXxU1d348c83OyQhkIUtARL2RVZD2ERRUcAN3MpmRW2lPlb7+PSpv+pj61KxauvTp7XighVxY1NRkNWWqiiiEpZAWAJhzYRAQiAJCdnn/P64A4SQYEJmy8z3/XrNa2bO3Ln3Ozfwveeec+65YoxBKaWU7wvwdABKKaXcQxO+Ukr5CU34SinlJzThK6WUn9CEr5RSfiLI0wFcTGxsrElMTPR0GEop1Wxs2rTpuDEmrq7PvDrhJyYmkpqa6ukwlFKq2RCRQ/V9pk06SinlJzThK6WUn9CEr5RSfsKr2/DrUllZic1mo6yszNOheJ2wsDASEhIIDg72dChKKS/U7BK+zWYjMjKSxMRERMTT4XgNYwz5+fnYbDaSkpI8HY5Sygs1uyadsrIyYmJiNNnXIiLExMTomY9Sql7NLuEDmuzroftFKXUxzTLhK6W8g91uWLTxMN/sPY7drlOte7tm14bvDSIiIiguLvZ0GEp5lDGGJ5el8/53hwFIaNOCnyR34s7kBDpEtfBwdKouWsP3QtXV1Z4OQamLMsbw+6VWsp95ZVf+NmUQnaNb8pd/7mHUC//mnrd/YHV6DhVVdk+H2uxk5hYz95sDLlm3Jvwm+PLLLxkzZgx33HEHvXv3Zvr06RhjWL16NXfeeed5y910000AfP7554wYMYIhQ4Zw5513nj1TSExM5Le//S1Dhgzhww8/5OWXX6Zv374MGDCAKVOmAFBSUsJ9991HSkoKgwcPZunSpe7/0crvGWN4cumOs8n+8Qm9mTgonvn3D2fdo1fz4Jju7Mop4oH3NzPi+bX8ceUuMnP1jPhijDF8tSePGXN/YOxfvuLF1bvJPeX8ARjNuknnmc92sPNIkVPX2bdjK566uV+Dl9+yZQs7duygY8eOjBo1ivXr1zN27FhmzpxJSUkJ4eHhLFq0iClTpnD8+HFmzZrFv/71L8LDw3nxxRf5y1/+wpNPPglATEwMmzdvBqBjx44cOHCA0NBQCgoKAHjuuee45pprmDt3LgUFBaSkpDB27FjCw8Odug+Uqo8xhqeX7eC97w5x/+gkHp/Q+7zBAp1jWvKbcb14ZGwP1u3NY9HGLOZ+c4A56/aT3KUNk4d24sYBHWgZ0qxTj9OUVlTzyZZs5q4/QGZuMXGRofz6up5MG9aZ2IhQp2/PZXtdRCYBNwKtgLeMMZ+LyBjgWWAHsNAY86Wrtu8uKSkpJCQkADBo0CAOHjzIFVdcwfjx4/nss8+44447WLFiBX/605/46quv2LlzJ6NGjQKgoqKCESNGnF3X5MmTz74eMGAA06dPZ9KkSUyaNAmwzg6WLVvGSy+9BFhDVA8fPkyfPn3c9GuVPzPG8MxnO3lnwyF+fkUS/3NDn3pHhgUFBnBN73Zc07sdeafKWbLZxqKNWTz60Tae+WwnNw/syOShnRiYEOWXo8tyCkt5b8Mh5v9wmILTlfTr2Iq//GQgNw7oQGhQoMu226iELyJzgZuAXGPMZTXKxwN/AwKBfxhjXjDGfAp8KiJtgJeAzwEDFANhgK2pwTemJu4qoaHnjsKBgYFUVVUBMGXKFF555RWio6NJTk4mMjISYwzXXXcdCxYsqHNdNWvqK1asYN26dXz22Wc899xzbN++HWMMH3/8Mb169XLtj1KqljPJft63B/nZFUk8cWP9yb62uMhQfnFVN2Ze2ZXUQydZ+EMWn2yxseCHw/RuH8lPkjtx6+B42oSHuPhXeN7WrALmfnOAldtzsBvD9X3bc98VSQxNbOOWA19j2/DnAeNrFohIIDAbmAD0BaaKSN8ai/zO8TnA18aYCcBvgWcuJeDm4qqrrmLz5s28+eabZ9vghw8fzvr168nMzASsNvk9e/Zc8F273U5WVhZXX301L774IoWFhRQXFzNu3Dj+/ve/Y4w1/G3Lli3u+0HKbxlj+MNyK9nfNyqJ3zUi2dckIgxNjOZ/fzKQH54Yy3O3XkZoUAB/WL6TYX9cy62vrufJpeksTs1i99Eiqqp9o8O3qtrO8m1HuO3V9UyavZ4vdudyz8hEvnr0al7/6eWkJEW77SynUTV8Y8w6EUmsVZwCZBpj9gOIyEJgoojsAl4AVhljNju+f+YveBKos4FKRGYCMwE6d+7cmPC8SmBgIDfddBPz5s3jnXfeASAuLo558+YxdepUysvLAZg1axY9e/Y877vV1dXcddddFBYWYozhV7/6Fa1bt+b3v/89jzzyCAMGDMBut5OUlMTy5cvd/tuU/zDG8OzyXby9/iD3jkrk9zddWrKvrVVYMNOHdWH6sC7syini0y3ZbMkq4ONNNt7dYE3nHhYcQJ8OrRgQH8Vl8VH0T4iie1wEQYHNY6xJwekKFm7M4t1vD3KksIwuMS15+ua+3JHciYhQz/RhyJnaYoO/YCX85WeadETkDmC8Mebnjvc/BYYBe4AZwEZgqzHmdRG5DRgHtAZe+7E2/OTkZFP7Bii7du3SNuuL0P2jnMUYw6wVu3jrmwPcMzKRp27u6/KaqN1u2H+8hPTsQrY7HjuyCympsIYqhwUH0LdDK/p78UEgM7eYed8e4ONN2ZRWVjOyWwz3jUri6t5tCQxwfU1eRDYZY5Lr+sxlhxljzMvAy7XKlgBLXLVNpZRzGGP440r3JnuAgAChe9sIureNYNLgeOD8g8A2WyHp2YV8tMnGOzXOBPp2aMVl8VG0axVGTHgIbcJDiAkPITo8hJjwUCLDgghwQrI9VVbJ0cIycgrLzj0XlZ73vrC0kpCgACYN6si9o5Lo06FVk7frLM5I+NlApxrvExxlSqlmyBjD86t28+bXB5gxoovbkn196joIVNsNB44XW2cBtiK2ZxewZHM2xeVVda4jMEBo0/LcQSA6IoTolo4DQoSjrGUIEWFBHC8uPz+hF5aRU1jK0cKys2caNcVGhNA+KoyENi0ZmhhNl5iWTBoc75JhlU3ljIS/EeghIklYiX4KMM0J61VKuZkxhhdW7WbOuv38dHgXnr6ln1cOmwwMELq3jaR720huHXyuvKyymvySCk4UV5BfUs7J0xXkF1dwosR65JdUcLKkgl1HisgvqaCwtLLebQQItI0Mo31UGD3bRXJlzzg6RIXRPqqF9dwqjLatQl06jNLZGjsscwEwBogVERvwlDHmLRF5CFiDNSxzrjFmh9MjVUq5lDGGF1bv5o11+7lreGf+MNE7k/3FhAUHEt+6BfGtGzaXT2W1nZOnKzhZUkl+STnFZVXERobSISqMuIhQr+obcIbGjtKZWk/5SmClUyJSSrmdMYYXV2fwxlf7mT6sM3+45bJml+wvRXBgAG0jw2gbGQZEejocl/Otw5dSqtGMMfx5TQavf7WPacM68+zEy5zSwam8j05ocQmee+455s+fT2BgIAEBAbzxxhsMGzaszmXHjBnDSy+9RHJynaOklLokmw6d4A/Ld3GkoJQYR8djTHgo0eEhxEaEEB0eSkzE+a8jQ4MuqLUbY3jp8wxe/XIfU1M6M0uTvU/ThN9IGzZsYPny5WzevJnQ0FCOHz9ORUWFp8NSfuJESQUvrNrF4lQbHaLCuKZXW06criC/uJy0kwWcKK7gVD0jVUICA84blRIbEUpZZTWr0o8yNaUTz03SZO/rNOE3Uk5ODrGxsWfn0ImNjQVg7dq1/OY3v6GqqoqhQ4fy2muvnTfPzuuvv86+ffv485//DMC8efNITU3llVde4f333+fll1+moqKCYcOG8eqrrxIY2Hx6/pXr2e2GxalZvLB6N8VlVfziyq786toehNdxxWZZZfXZUSnHi8vPjlI5XnLudX5xOfvzSigsrXQMveynyd4PNO+Ev+oxOLrduets3x8mvFDvx9dffz1/+MMf6NmzJ2PHjmXy5MkMGzaMe+65h7Vr19KzZ0/uvvtuXnvtNR555JGz37v99tsZMWLE2YS/aNEinnjiCXbt2sWiRYtYv349wcHBPPjgg3zwwQfcfffdzv1dqtnaeaSI3326nc2HC0hJjObZSZfRq339HYxhwYF0bN2Cjg0cqaL8h3baNlJERASbNm1izpw5xMXFMXnyZN544w2SkpLOzokzY8YM1q1bd9734uLi6Nq1K9999x35+fns3r2bUaNGsXbtWjZt2sTQoUMZNGgQa9euZf/+/Z74acrLnCqr5A+f7eTmV77hUP5pXrpzIIt+MfyiyV6pi2neNfyL1MRdKTAwkDFjxjBmzBj69+/P7Nmzf/xLWFMmL168mN69e3PrrbciIhhjmDFjBs8//7yLo1bNhTGGFdtzeHb5TnJPlTMtpTOPjutF65a+P32wci2t4TdSRkYGe/fuPft+69atdOvWjYMHD56d9vi9997jqquuuuC7t956K0uXLmXBggVnp0y+9tpr+eijj8jNzQXgxIkTHDp0yA2/RHmjA8dLuHvuDzw0fwuxEaEs+Y+RPHdrf032yimadw3fA4qLi3n44YcpKCggKCiI7t27M2fOHKZOncqdd955ttP2gQceuOC7bdq0oU+fPuzcuZOUlBQA+vbty6xZs7j++uux2+0EBwcze/ZsunTp4u6fpjyorLKaV7/cx+tf7iM0KIBnbunHXcO7uGV2ReU/Gj09sjvp9MiNp/un+fkiI5enlu7g8InTTBzUkSdu7OO48lOpxvPI9MhKeaO0rAKOFpUxplecxye9OlJQyrPLd7Iq/Shd48KZ//NhjOwe69GYlG/ThK98XkWVnZXbc3j724OkZRUAEBsRyl3DOzNtWGe316azC0pZvDGLN7/eT7Xd8Oi4Xvx8dJLHD0DK9zXLhG+M8YuJnRrLm5vnPCG3qIz3vz/M/O8Pc7y4nK6x4Tx9c1+6xITz7oaD/PVfe5n9RSY3D7BuVNE/IcplsZRXVfPPncdYnGrj6715GAPX9W3Hkzf1pVN0S5dtV6maml3CDwsLIz8/n5iYGE36NRhjyM/PJyzMv9t+jTFsPlzAO98eZOX2HKrshqt7xTFjZCJX9og7ezXp1b3bsj+vmHc3HOLD1CyWbMnm8i5tuHdUIuP6tSfYSdPi7jxSxOLULD7dmk3B6Uo6RoXx8NXduTO5kyZ65XbNrtO2srISm81GWVmZh6LyXmFhYSQkJBAcHOzpUNyuvKqa5Wk5vLPhINtshUSGBnFncifuHtGFxNjwi363qKySD1NtvPPtQQ6fOE2HqDDuGt6FqSmdiQ5v/HDIwtJKlm3NZnGqje3ZhYQEBnBdv3ZMTu7EqO6xOvJGudTFOm2bXcJXqqajhWV88P0h5n9/mPySCrrFhXPPyERuG5JQ5zwzF1NtN3yxO5e3vz3A+sx8QoMCmDQonntGJf7ofUntdsN3+/NZlJrF6vSjlFfZ6d0+kslDOzFpUDxtLuHAodSl0ISvfIoxhk2HTvL2twdZk36UamO4tndbZoxM5IrusU5p6ttz7BRvrz/IJ1tslFXaGd41mntHJTG2T7vzaujZBaV8lGrjw01Z2E6WEhkWxKRB8Uwe2ol+HVtps6NyO034yieUV1WzdOsR3vn2IDuOFNEqLIjJQzvx0+GJdI5xTXt4wekKFm7M4t1vD3KksIyENi2YMSKR9lFhLE7N4pvM4xgDo7rH8JPkTozr156wYB1tozxHE75q9r7ak8fTy3Zw4HgJPdtFMGNkIrcOjqdliHvGHVRV2/nnzmO8vf4gPxw8AUDHqDDuSO7EnZcnaAes8hpuv/BKRCYBNwKtgLeMMZ+LSDjwKlABfGmM+cAV21a+JbuglFmOi5OSYsN5+56hjOkV5/amkqDAACb078CE/h3YeaSIwtJKUpKitQNWNSsNTvgiMhe4Ccg1xlxWo3w88DcgEPiHMeYFY8ynwKci0gZ4CfgcuA34yBjzmYgsAjThq3pVVNn5xzf7+fvaTAzedXFS344X78BVyls1poY/D3gFePdMgYgEArOB6wAbsFFElhljdjoW+Z3jc4AE4MzdSqqbELPycV/vzeOpZTvYn1fCuH7t+P1NfUloo00mSjVVgxO+MWadiCTWKk4BMo0x+wFEZCEwUUR2AS8Aq4wxmx3L2rCS/lYuMi2ziMwEZgJ07ty5oeEpH3CkoJRZK3aycvtREmNa8va9Q7m6V1tPh6WUz2hqG348kFXjvQ0YBjwMjAWiRKS7MeZ1YAnwiojcCHxW3wqNMXOAOWB12jYxPuUEmbmnmP99Fr3aRzC6R5zTb51XUWXnrW8O8PLavdiN4b+v68n9V3bV0S5KOZlLOm2NMS8DL9cqKwHudcX2lOt8ssXGE5+kU1ZZjd1x+O0WF87oHnGM7hHLsK4xRDTyAqea1mce58ml6ezLK9G5ZZRysaYm/GygU433CY4y1cyVVVbzzGc7WPBDFimJ0bw8dTBFZZWs25PH13uPs3DjYeZ9e5CgAGFIlzaM7h7L6J5x9I+PatDIlaOFZTy7YicrtuXQJaYlb98zlKt7a/ONUq7UqHH4jjb85WdG6YhIELAHuBYr0W8EphljdjgjOB2H7xn784r55fwt7Mop4sEx3fj1dT0JqjWZWFllNZsPnWTd3uN8k5lHenYRAFEtghnVPYbRPeK4onvsBbX1iio7b68/wN/W7qXabvjl1d2Zqc03SjmNU8bhi8gCYAwQKyI24CljzFsi8hCwBmtY5lxnJXvlGcu3HeGxj7cTFCgXrXWHBQcysnus44YdvckvLmf9vny+3pPHN5nHWbn9KABJseFc0T2W0T1iCQkKYNaKXWTmFjO2Tzueulmbb5RyJ73SVgHWtAWzlu/ive8OMaRza16ZNuSSO2eNMezLK2bdnuN8k3mc7/bnc7rCGonbKboFT9/cj2v7tHNm+EopB73Fobqow/mneXD+JtKzi7h/dBL/b3zvJs0HLyJ0bxtJ97aR3HdFEhVVdjYfPsmRglJu6N9Bm2+U8hBN+H5udfpRHv0oDQHm/PRyru/X3unbCAkKYHjXGKevVynVOJrw/VRFlZ0XVu1m7voDDEyI4pVpQ7Q9XSkfpwnfD9lOnuah+VvYmlXAPSMTefyG3l4xR41SyrU04fuZtbuO8evFadjthlenD+GG/h08HZJSyk004fuJymo7L63J4I11++nboRWvTh/yo/d6VUr5Fk34fiCnsJSH528h9dBJpg3rzJM39dWRMkr5IU34Pu6L3bn894dplFVW87cpg5g4KN7TISmlPEQTvo8qr6rmxVUZzF1/gN7tI3ll2hC6t43wdFhKKQ/ShO+D9uUV8/D8LezMKWLGiC48fkMfbcJRSmnC9yXGGD7aZOOpZTsICQrgzbuTua6vTmGglLJowvcRp8oqeeKTdJalHWF412j+Onkw7aPCPB2WUsqLaML3AVsOn+RXC7dwpKCM/76uJw9e3b1Bc9IrpfyLJvxmzG43vLFuP//7eQbtWoWxaOZwkhOjPR2WUspLacJvpnKLyvj14jS+yTzODf3b8/ytA4hqGezpsJRSXkwTfjP0RUYuv1mcRklFFc/f1p8pQzshok04SqmL04TfjJRXVfOn1Rm89Y01tn7h1OH0aBfp6bCUUs2EJvxmYn9eMb9auIX07CLuHtGF/9Gx9UqpRtKE7+WMMSzZnM3vl6YTEhTgspuUKKV8n8sSvoh0BZ4AoowxdzjKxgDPAjuAhcaYL121fV9QWlHN/3yynU+2ZDMsKZq/ThlEh6hLu8+sUko16salIjJXRHJFJL1W+XgRyRCRTBF5DMAYs98Y87NaqzBAMRAG2JoSuK/LLSpj8pwNfLo1m/8a25P59w/XZK+UapLG3ql6HjC+ZoGIBAKzgQlAX2CqiPSt5/tfG2MmAL8Fnmnktv3GziNFTJy9nszcYub8NJn/HNtDL6RSSjVZoxK+MWYdcKJWcQqQ6ajRVwALgYn1fN/ueHkSCK1rGRGZKSKpIpKal5fXmPB8wtpdx7jj9W8xBj58YITOhaOUcprG1vDrEg9k1XhvA+JFJEZEXgcGi8jjACJym4i8AbwHvFLXyowxc4wxycaY5Li4OCeE1zwYY3jrmwPc/24q3eIiWPrQKPp1jPJ0WEopH+KyTltjTD7wQK2yJcASV22zuaqstvP0sh188P1hxvVrx/9NHkTLEB1ApZRyLmdklWygU433CY4y1QBFZZX88oPNfL33OL+4qiu/HdebAG2vV0q5gDMS/kagh4gkYSX6KcA0J6zX52WdOM198zZy4HgJL97en8lDO3s6JKWUD2tUwheRBcAYIFZEbMBTxpi3ROQhYA0QCMw1xuxweqQ+ZtOhE8x8dxNVdsO7P0thZLdYT4eklPJxjUr4xpip9ZSvBFY6JSI/sHRrNo9+tI0OUWHMvWco3eL0XrNKKdfTnkE3Msbwt7V7+eu/9pKSFM0bd11Om/AQT4ellPITmvDdpKyymsc+3sanW49w+5AE/njbZYQG6eRnSin30YTvBvnF5cx8bxObDp3k0XG9eHBMN52/XinldprwXWzvsVPc985GcovKmT1tCDcO6ODpkJRSfkoTvgt9s/c4//HBJkKDAlk4cziDO7fxdEhKKT+mCd9F9h47xX3zNtI1Lpx/zEgmoU1LT4eklPJzmvBdwG43PL5kOy1DA3n/58OIjahznjillHIrZ0yepmqZ/8NhUg+d5Ikb+miyV0p5DU34Tna0sIwXV+1mVPcY7rg8wdPhKKXUWZrwneypZelUVNt5blJ/HXqplPIqmvCdaHX6UdbsOMYjY3uSGBvu6XCUUuo8mvCdpKiskqeWpdOnQyt+PjrJ0+EopdQFNOE7yZ9W7ybvVDkv3Naf4EDdrUop76OZyQlSD57g/e8Oc8/IJAZ2au3pcJRSqk6a8JuovKqax5dsJ751C/77+p6eDkcppeqlF1410etf7mdvbjFv3zOU8FDdnUop76U1/CbIzD3F7C8yuXlgR67u3dbT4Sil1EVpwr9Edrvhf5ak0yIkkCdv6uvpcJRS6kdpwr9ECzdm8cPBEzxxYx/iInX6BKWU93NJwheRriLyloh8VKMsXETeEZE3RWS6K7brLrlFZTy/ahcjusZwp06foJRqJhqc8EVkrojkikh6rfLxIpIhIpki8hiAMWa/MeZntVZxG/CRMeZ+4JYmR+5BTy3bQXmVnT/eptMnKKWaj8bU8OcB42sWiEggMBuYAPQFpopIfQ3aCUCW43V148L0Hp/vOMqq9KP857U9SNLpE5RSzUiDE74xZh1wolZxCpDpqNFXAAuBifWswoaV9C+6XRGZKSKpIpKal5fX0PDc4lRZJU8u3UHv9pHMvLKrp8NRSqlGaWobfjznau1gJfV4EYkRkdeBwSLyuOOzJcDtIvIa8Fl9KzTGzDHGJBtjkuPi4poYnnP9eU0Gx06V8cLtA3T6BKVUs+OSK4WMMfnAA7XKSoB7XbE9d9h06CTvfXeIGSMSGaTTJyilmqGmVlOzgU413ic4ynxKRZWdx5dso0OrMH4zrpenw1FKqUvS1IS/EeghIkkiEgJMAZY1PSzv8sZX+9hzrJhZt15GhE6foJRqphozLHMBsAHoJSI2EfmZMaYKeAhYA+wCFhtjdrgmVM/Yl1fM3/+dyU0DOnBN73aeDkcppS5Zg6urxpip9ZSvBFY6LSIvYrcbHl+ynbDgAJ68WadPUEo1bzrU5CIWp2bxwwFr+oS2kWGeDkcppZpEE349ck+V8ceVuxjeNZqfJHf68S8opZSX04Rfj2eX76Ksys7ztw3Q6ROUUj5BE34dSiuqWbk9hxkjuuj0CUopn6EJvw47jhRSbTcMS4rxdChKKeU0mvDrsDWrAIABnaI8G4hSSjmRJvw6bLMV0jEqTEfmKKV8iib8OqTZChio8+UopXyMJvxaCk5XcCj/NAMSWns6FKWUcipN+LWk2QoBGKjt90opH6MJv5ZtWQWIQP94TfhKKd+iCb+WNFshXWPDiQwL9nQoSinlVJrwazDGaIetUspnacKv4WhRGXmnyhmoHbZKKR+kCb+GNMcFV1rDV0r5Ik34NaTZCgkOFPp0iPR0KEop5XSa8GtIyyqgd/tWhAYFejoUpZqHnDRY8gvY8CpUVXg6GvUjNOE72O2G7bZCHX+vVEOcOgqf/hLeuAp2LoU1j8OrwyFjFRjj6ehUPTThOxzIL+FUeZVeYavUxVSWwld/hpeHwLZFMPIh+O/dMO1DkABYMAXemwTHfOrW1j7DbQlfRMaIyNci8rqIjHHXdhvqTIftIO2wVepCxsC2D+HvyfDFLOh+DTz0A1w/C1q0hp7Xw4MbYPyLcGQrvH4FLP8vKDnu6chVDU1K+CIyV0RyRSS9Vvl4EckQkUwRecxRbIBiIAywNWW7rrDNVkjLkEC6xUV4OhSlvEvWD/CPsbDk5xAeA/esgMnvQ3TX85cLDIbhD8CvtkDKTNj0Drw8GL79u7bve4mm1vDnAeNrFohIIDAbmAD0BaaKSF/ga2PMBOC3wDNN3K7Tbc0qoH98FIEBejtDpQAoOAwf3gtvXQeFNpj4Ktz/JSRecfHvtYyGCS9aNf7Ow+Hz38Grw2D3Cm3f97AmJXxjzDrgRK3iFCDTGLPfGFMBLAQmGmPsjs9PAqH1rVNEZopIqoik5uXlNSW8BquosrMzp0jH3ysFUH4K/vWM1XyTsQqu/H/w8CYYPB0CGpEy4nrB9A9h+scQEAwLp8G7t8DR9B//rnKJIBesMx7IqvHeBgwTkduAcUBr4JX6vmyMmQPMAUhOTnZLdSDj6CkqquwMSNAROsqP2athy/vw71lQkgv9fwJjn4KohKatt8dY6HoVbJoHXzwHb4yGIXfD1b+DiDinhK4axhUJv07GmCXAEndtrzHSbAUAOqWC8l/7v4I1T8Cx7dBpGExdCAmXO2/9gcGQcj/0vwO++hP8MAe2fwxXPQrDHoCgek/6lRO5YpRONtCpxvsER5nXSssqIDo8hIQ2LTwdiroYY6A4z6qJKufIy4AFU62mlrJCuGMu3LfGucm+phZtYPzz8OB3kDgK/vkkzE6BHZ9CcS5UnNZ2fhdyRQ1/I9BDRJKwEv0UYJoLtuM022yFDEyIQkQ7bL1W+Sn4+OewZzVIILSKh9adIKpTrefOVhNEsN6PuE7VlZD1vbUf93wOxzMgJAKufRKG/9J9+y22B0xbBJlrrTOLD2fU+FCsmELCIdTxHFLPc2hEjdeREN0NYrpDUIh7focrlBVaB+JOKU5fdZMSvogsAMYAsSJiA54yxrwlIg8Ba4BAYK4xxmuvwigpr2Jv7inGX9be06Go+hRmw/zJkLsTRj1iXeBTaIPCLDj4DZw6AmfHBDiEt73wQHDmfXSSlSD8RUk+ZP7TSvKZ/4byQqsTtctIuHwG9L8TItp6Jrbu10LSVbBnlXX1bkUxVJQ4HsVQXuN9cS5U7D//89p/d7AqBDHdIK43tO1jPcf19s4DQVUF5O4AWypkb4bsVDi+x/o3/rjN6f9Om5TwjTFT6ylfCaxsyrrdJT27ELvRC6681pEtMH+K9R982mKrA7C26kooOmIdAAqyHM+Hreej262RJtXl55YPCoM+N8Og6ZB0JQT42NxJxli/e+8a2LPGSiYY6yDY52brIqmuV0NYK09HagkMsuJqLGOgqsxxUCiGsgI4ngl5uyB3NxxLh12fYV0CBAQEWWcAbXtDXB9rFFHbPlaZOw4ExsDJg5C9yXrYUq25iM7822wZCwnJVmd5/BDroOxkbuu09VbbHPew1RE6Xmj3CqsZp2UM/GwNtOtX93KBwdCmi/Woi90OJXmOs4LD1lnB9g+tR6sEGDgFBk2zaoXNVUWJ1fG6d43VVHPqiFXecTBc9VvoOQ46DGrcsEpvJwLBLawHjtE+HQefv0xlKRzfC3m7IXeX9ZyzDXYu47wDQUx36ywgtofVzxAaaT1CHM+hEeeXBTYgdZ4+ca7WfibJn863PgtqAR0HWR3Z8Zdbj9adrd/kQmK8uIMkOTnZpKamunQbv5y/ma2HC1j/2DXOXXF1JWz/CBKGQmx3567b1xkDG2ZbF+x0HGyNGIls59xtVJZBxkrYOh/2rbWaBjqPsGr9/SZZ/7G9md0OJ/bBvi+sJH/ga6umGBIB3a6GHuOgx/XO32++orLUajrJ3X3ujCBvt1UDpwE5MajFhQeBMwcGY7fOTE/sdyws1plE/BCIT7aSe9s+VkXFBURkkzEmuc7P/D3hj/7TvxkQ35rZ04c4b6WHNljziOTtgsAQuPJRq+3Z29oPvVF1Faz8DWx6G/rcAre+ASEtXbvNoiPWRGBbPoD8vRDcEvpOtGr9Xa7wfK24uspKTjlp5x5Ht1nNGGA1SfQcZz06j9R/Z01RXeXoOzhlPSqKobzIajY6+/5Urc9POT4vsl4bO3QY4Ki5J1s1eTdWIC6W8P26SSe/uJysE6XcNayepoDGOn3CGma25T2rc/C2f1i1yC+eg/SP4ea/WZeaq7qVFcKH98C+f1sHyGufck+ybdURrvgva5u2VNj6gfX3SltgnWYPnAaDpkKbRNfHUlVudU7XTO7Hdlht1WDVLNv3h4FTocNA66xEzyCdJzDImgyuRWtPR+ISfl3D/yIjl3vf3sjCmcMZ3jXm0ldkjNU08PnvrKQ18iGr3fRMD/uez2HFr61OxMvvhbFP++w/qEt28pA1Eid/L9z4F2v0iCdVllp9CFveh/1fAgYSR1u1/r4TnTN6oqLESuY5aZCz1XrO3QX2Kuvz0FbQfoCV2M88Ynv4Xiezciqt4ddjW1YhInBZfBM6bPMyrOabQ+utKxRv+r8LOxd7Xg9dvoMvn4fvXrVq/RP+ZCUOV4/9P7oddq+0OqLaJFqP1p29a5y6LdWaR72qAu76GLqO8XREVkdg/zusR6HNqu1vnQ+f/gesfNQ6XcdYbemm2roY7Lznusrt578vPXFuWGHLGCuhjxx7Lrm3TvR8c5LyKX6d8NNsBfRoG0FE6CXshorT8PVLsP5lq7Z388sw+Kf1/wcNjYBxz1kJZNmvrAtNek6AG19q+lwlF8RWAulLrHbw7E11LCBWM8aZA8B5jyQIj3X9geiMHZ/CJ7+AiHbWtLtxvdyz3caISrD6YUb/Bg5/ZzX55O6yatoSaI3yCAyp8f7Mc0Ct97XKI9qeS+6t4t23z5Xf8tuEb4xhm62AMb0u4YKTvf+ymmgKDlltqdc92/BJoDoOhvu/gO9fgy/+CLOHwTW/t4ZnNfVU/eh2a4KqbYutDqTYXjD+BRgw2Ro1dPJgjccB63nfv+FUzvnrCQ6/8EAQ3dXqfAqPbVqMZxgD3/wfrH3GOjOaMt9563YVEegywnoo1Qz5bcLPLijleHEFAxsz/r4oB1Y/Bjs/hZgeMOMz68KdxgoMgpEPW6NQVvwaVv/WGiVyy8tWh1xj1K7NB4ZCv1vh8nusDuKatcbIdtB52IXrqCy1LlQ6eRBOHDh3UDix3zogVJWeW7Z153NDy+Ivt2qnjR1FU1Vh/e4t78Flt1vzrHtTE5NSPspvE/6ZC64aNAe+vRo2/gPWPgvVFda0rqN+1fQZ/tp0gekfWSNCVj9m3RB65MOODt8fSaIXq823jG5cHMEtrKaUuppTjLEuaT++xxpbnJ0Kto2wwzHxqQRCu77nhqDFX26tp76zldKTsPhuOLDOmmd9zOPaTq2Um/htwk/LKiAkMIDe7X/k8vIjW+CzR6xRFN2ugRtecu4VmSJWu363a6whnev/Cjs+sTp/u197/rKNqc07M77IdtYjafS58lPHzl09mL0J0j+xDkBgNQl1HGzNuHjmTKBVvNWM9MFPrLOHSa9bQx2VUm7jt8Myp8zZQGmlnaW/HFX3AmWF8O/nYOObEB5nTena7zbXd6wd+BqWPwL5mVZtfdwfrTb2mrX5uN5Wkr+U2ryrnLny88wcIdmbrLMQe6X1eUR7q2lIAmDyB9bUuEopp9NhmbVU2w3p2UXcNiS+7gWMgbnjrZEYKffDNb+DMDfNtZM0Gh5YD1//r9WpuXOpddHNmdp88r1WJ6e3jegICLDGiMf2sOamAesioqPbz50FlBVaB7DmPGeNUs2YXyb8/XnFFJdX1X+HqyNbrKsdb/o/SL7PrbEBVgfmNU9YHZob/g7tLvOu2nxDBYVas/8l1FnZUEq5mV8m/LSzHbb11NozVllND30mujGqOrTtDRNnezYGpZTP8MvhEWlZBUSEBtE1NqLuBTJWWnOUhDdhugWllPIyfpnwt9kK6B8fRUBAHe3gJw9ZN07oNcH9gSmllAv5XcIvr6pmZ04RAy7WnAPQ6wb3BaWUUm7gdwl/d84pKqsNg+rrsM1YaV3EpCNJlFI+xm0JX0TCReQdEXlTRKa7a7u1pdkKABhQ1xW2pQXWrJfanKOU8kFNSvgiMldEckUkvVb5eBHJEJFMEXnMUXwb8JEx5n7glqZstynSsgqJjQilY1Qdc7dk/suai7z3je4PTCmlXKypNfx5wPiaBSISCMwGJgB9gaki0hdIALIci1U3cbuXbJutgIEJUUhdFy7tXmFdVRt/ufsDU0opF2tSwjfGrANO1CpOATKNMfuNMRXAQmAiYMNK+hfdrojMFJFUEUnNy8trSngXKC6vIjOvmAF1td9XVVg1/J7j9Y5CSimf5Io2/HjO1eTBSvTxwBLgdhF5Dfisvi8bY+YYY5KNMclxcQ2cY76BttsKMaaeC64OrbfmqdHROUopH+W2K22NMSXAve7aXl3OdtjWVcPPWGndINobbq+nlFIu4IoafjbQqcb7BEeZx22zFdA5uiXR4SHnf2CMNf6+29WNv5mHUko1E65I+BuBHiKSJCIhwBRgmQu202hpWYUMqOsOV0e3Q2GWNucopXxaU4dlLgA2AL1ExCYiPzPGVAEPAWuAXcBiY8yOpofaNMeLy8kuKGVQXePvM1YBAj3HuTsspZRymya14Rtj6rxlkTFmJbCyKet2tm0Xbb9fAZ1SIOISbmiulFLNhN9MrbA1q5AAgcvia93SsDAbctL06lqllM/zm4S/zVZAz3aRtAypdVKT4TgR0fZ7pZSP84uEb4whLaug7g7bjFUQ3Q1ie7o/MKWUciO/SPi2k6WcPF3JwNodtmVFcGCd1ZzjbfeIVUopJ/OLhL81qwDgwnvY7lsL9kqdLE0p5Rf8IuFvsxUQEhRAr/aR53+QsQpaRENCimcCU0opN/KLhJ+WVUi/jq0IDqzxc6srYc8aa7K0QL+8l7tSys/4fMKvthvSjxRe2Jxz+DsoK9DhmEopv+HzCT8zt5jTFdUXzpCZsRICQ6HbNZ4JTCml3MznE36ao8P2vCtsjbFudtL1KgiN8EhcSinlbr6f8G0FRIYFkRQTfq4wdxcUHNKLrZRSfsUvEv6AhCgCAmqMsz9zdW3P8XV/SSmlfJBPJ/yyymp255y6cMK0jFXQcQi06uCRuJRSyhN8OuHvzCmiym7OH6Fz6ihkp0Jvbc5RSvkXn074285cYVtzhM6e1daztt8rpfyMbyd8WyFtI0Np3yrsXOHuldC6C7Tt67nAlFLKA3w64W+1FTAgoTVyZmK0ihLY/6VVu9fJ0pRSfsZnE35RWSX780oYVLM5Z9+/obpc2++VUn7JbQlfRMaIyNci8rqIjHH19rbbCoFaF1xlrIKwKOg8wtWbV0opr9OghC8ic0UkV0TSa5WPF5EMEckUkcd+ZDUGKAbCANulhdtwaWfvYeuo4durrQ7bHuMgMNjVm1dKKa/T0Gki5wGvAO+eKRCRQGA2cB1WAt8oIsuAQOD5Wt+/D/jaGPOViLQD/gJMb1roF5eWVUBiTEtatwyxCrJ+gNP5OlmaUspvNSjhG2PWiUhireIUINMYsx9ARBYCE40xzwM3XWR1J4HQS4i1UbbZChmaGH2uIGMlBARD97Gu3rRSSnmlpkwEHw9k1XhvA4bVt7CI3AaMA1pjnS3Ut9xMYCZA586dLymw3KIycgrLzr+lYcZKSBoNYa0uaZ1KKdXcue3OH8aYJcCSBiw3B5gDkJycbC5lW2mODtuBZ9rvj++F/EwY9sClrE4ppXxCU0bpZAOdarxPcJR53DZbAYEBQr+OjoS/e4X1rO33Sik/1pSEvxHoISJJIhICTAGWOSesptmaVUDPdpG0CAm0CjJWQfsBEJXg2cCUUsqDGjoscwGwAeglIjYR+Zkxpgp4CFgD7AIWG2N2uC7UhiurrGZw59bWm+I8yPoeet/o0ZiUUsrTGjpKZ2o95SuBlU6NyAk+fGAkxjia//euAYw25yil/J7PTq1wdv6c3SuhVYLVpKOUUn7MZxM+AJWl1vw5vSboZGlKKb/n2wl//1dQVaqTpSmlFL6e8DNWQGgr6HKFpyNRSimP892Eb7dDxmprKoWgEE9Ho5RSHue7CT97E5Tk6q0MlVLKwXcTfsZKkEDooZOlKaUU+HrCTxwFLdp4OhKllPIKvpnw8/dB3m5tzlFKqRp8M+FnrLKe9epapZQ6y3cTftt+0CbR05EopZTXcNt8+G5jDHQYAK0v7eYpSinlq3wv4YvA+Nq31FVKKeWbTTpKKaUuoAlfKaX8hCZ8pZTyE5rwlVLKT2jCV0opP6EJXyml/IQmfKWU8hOa8JVSyk+IMcbTMdRLRPKAQ56O4yJigeOeDqKBmkusGqdzNZc4ofnE6u1xdjHGxNX1gVcnfG8nIqnGmGRPx9EQzSVWjdO5mkuc0HxibS5x1kWbdJRSyk9owldKKT+hCb9p5ng6gEZoLrFqnM7VXOKE5hNrc4nzAtqGr5RSfkJr+Eop5Sc04SullJ/QhN9AItJJRL4QkZ0iskNE/tNR/rSIZIvIVsfD43dOF5GDIrLdEU+qoyxaRP4pInsdz208HGOvGvtsq4gUicgj3rI/RWSuiOSKSHqNsjr3oVheFpFMEdkmIkM8HOefRWS3I5ZPRKS1ozxRREpr7NvXPRxnvX9rEXncsT8zRGSch+NcVCPGgyKy1VHusf15yYwx+mjAA+gADHG8jgT2AH2Bp4HfeDq+WrEeBGJrlf0JeMzx+jHgRU/HWSO2QOAo0MVb9idwJTAESP+xfQjcAKwCBBgOfO/hOK8HghyvX6wRZ2LN5bxgf9b5t3b8v0oDQoEkYB8Q6Kk4a33+v8CTnt6fl/rQGn4DGWNyjDGbHa9PAbuAeM9G1SgTgXccr98BJnkulAtcC+wzxnjNVdXGmHXAiVrF9e3DicC7xvId0FpEOngqTmPM58aYKsfb74AEd8RyMfXsz/pMBBYaY8qNMQeATCDFZcHVcLE4RUSAnwAL3BGLK2jCvwQikggMBr53FD3kOH2e6+mmEgcDfC4im0RkpqOsnTEmx/H6KNDOM6HVaQrn/yfytv15Rn37MB7IqrGcDe+pDNyHdfZxRpKIbBGRr0RktKeCqqGuv7W37s/RwDFjzN4aZd62Py9KE34jiUgE8DHwiDGmCHgN6AYMAnKwTvk87QpjzBBgAvBLEbmy5ofGOh/1ivG4IhIC3AJ86Cjyxv15AW/ah/URkSeAKuADR1EO0NkYMxj4NTBfRFp5Kj6ayd+6hqmcXzHxtv35ozThN4KIBGMl+w+MMUsAjDHHjDHVxhg78CZuOvW8GGNMtuM5F/gEK6ZjZ5oZHM+5novwPBOAzcaYY+Cd+7OG+vZhNtCpxnIJjjKPEZF7gJuA6Y6DE44mknzH601YbeM9PRXjRf7W3rg/g4DbgEVnyrxtfzaEJvwGcrTfvQXsMsb8pUZ5zbbaW4H02t91JxEJF5HIM6+xOvDSgWXADMdiM4ClnonwAufVmrxtf9ZS3z5cBtztGK0zHCis0fTjdiIyHvh/wC3GmNM1yuNEJNDxuivQA9jvmSgv+rdeBkwRkVARScKK8wd3x1fLWGC3McZ2psDb9meDeLrXuLk8gCuwTuG3AVsdjxuA94DtjvJlQAcPx9kVa4RDGrADeMJRHgOsBfYC/wKivWCfhgP5QFSNMq/Yn1gHoRygEqsN+Wf17UOs0TmzsWp424FkD8eZidUGfubf6euOZW93/JvYCmwGbvZwnPX+rYEnHPszA5jgyTgd5fOAB2ot67H9eakPnVpBKaX8hDbpKKWUn9CEr5RSfkITvlJK+QlN+Eop5Sc04SullJ/QhK+UUn5CE75SSvmJ/w918xKYSOtkHAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax =  plt.subplots()\n",
    "xs = np.arange(10, 200, 10)\n",
    "ys = [error(x) for x in xs]\n",
    "\n",
    "series = plt.plot(xs, ys)\n",
    "plt.legend(iter(series), ('Inverse', 'Solve'))\n",
    "ax.set_yscale(\"log\", nonpositive='clip')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7289414",
   "metadata": {},
   "source": [
    "## Error\n",
    "Additionally we can look at the quality of the inverse and you can see that the difference increases sharply for larger matrices (chart is in log-scale):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "cee36266",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD4CAYAAAAEhuazAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAkzElEQVR4nO3dd3xV9f3H8dfHsMMSwswgYRM2hOGqOFBExa0MNy11UEdbR7UtqHW0aq0KahEQqwgqWg0YxYWCWpU9khCIYSRh75n9/f1xr79fmh9gTG5ycu99Px8PHuR8703u+3jgzfXc8/0ec84hIiLh5SSvA4iISPVT+YuIhCGVv4hIGFL5i4iEIZW/iEgYquV1gPKKiopy8fHxXscQEQkqS5cu3eWca1F2PGjKPz4+niVLlngdQ0QkqJjZpmON67SPiEgYUvmLiIQhlb+ISBhS+YuIhCGVv4hIGFL5i4iEIZW/iEgYCnj5m1l7M5tmZnOON2Zm3czsJTObY2a3BjqDiEgo2Hkwn4nJqRzMKwz4zy5X+ZvZdDPbYWZryowPM7MMM8s0s/sBnHNZzrmxpZ9Xdsw5l+6cuwW4Gjit8rshIhI6CotLmLooi7Of+oKZ321iyca9AX+N8r7znwEMKz1gZhHAZOACIBEYZWaJ5X1hMxsBfACklPd7RERC3TeZuxj+7CL+8kE6fdudzEd3/YKzurYM+OuUa3kH59xCM4svMzwQyHTOZQGY2WzgEiCtnD8zGUg2sw+AN471HDMbB4wDiIuLK8+PFREJSrn7jvLYB+l8sHorsc3qM+W6/gxNbIWZVcnrVWZtn2ggu9R2DjDIzJoDjwJ9zewPzrnHy44B/wEuB+pygnf+zrkpwBSApKQk3W9SREJOXmExUxdlMWlBJs7B3ed25tdntqde7Ygqfd2AL+zmnNsN3PJTY8AXgX5tEZFg8ln6dh6am8bmPUcY1r01D17YjdhmDarltStT/rlAbKntGP+YiIicwIZdh3l4bioLMnbSoUUkr40dyBmd/t+qy1WqMuW/GOhkZgn4Sn8kMDogqUREQtCRgiImL8jk5YUbqB1hPDi8GzecGk+dWtU/5apc5W9ms4AhQJSZ5QATnHPTzGw8MB+IAKY751KrLKmISJByzjFv1VYeS0ln6/48Lu8bzf0XdKVl43qeZSrv1T6jjjOegi7VFBE5roxtB5mQvIZvs/aQ2KYxz4/qS1J8M69jBc+dvEREgsn+o4X849N1/Os/m2hYtxaPXNqD0QPjiDipai7d/LlU/iIiAVRS4pizLIe/fbSW3YcLGDUwjt+f14VmkXW8jvZfVP4iIgGyKmcff34/lRXZ++gX15QZNw2kR3QTr2Mdk8pfRKSSdh/K58n5Gby5JJvmkXV5+qreXNY3mpNqyCmeY1H5i4hUUFFxCTO/28zTH2dwuKCYsaclcMe5nWhcr7bX0X6Syl9EpAK+37CHP7+/hrXbDnJax+ZMvLg7nVo18jpWuan8RUR+hm3783j8w3TeX7GFtk3q8cKYflzQo3WVLcBWVVT+IiLlUFBUwvSvN/DcZ+spKnH85uyO3DakI/XrVO0CbFVF5S8i8hO+XLeTh5JTydp1mHO7teRPFyXSrnmk17EqReUvInIc2XuO8PC8ND5J205CVCSv3DSAs7oE/sYqXlD5i4iUcbSgmBe//IGXvvyBCDPuHdaFsacnULdWcJ7iORaVv4iIn3OO+anbeGReOrn7jnJx77Y8MLwrbZrU9zpawKn8RUSAzB2HeGhuKovW76JLq0bM+tVgTunQ3OtYVUblLyJh7WBeIc9/nsn0rzZQv04EEy5O5LrB7agVUf1r7Fcnlb+IhCXnHO+tyOWxlLXsPJjPNUmx3DOsC1EN63odrVqo/EUk7KRu2c+E91NZsmkvvWOa8PL1SfSJbep1rGoV8PI3s/bAg0AT59yV/rFI4AWgAN+N298rve2cmxnoHCIiZe07UsBTH2fwxnebadqgDn+9oidX9Y+t0QuwVZXy3sZxOnARsMM516PU+DDgWXy3cZzqnHvCOZcFjDWzOaV+xOXAHOfcXDN7EzipzLbKX0SqTHGJY/bizTw1P4P9Rwu5/pR47j63M00a1PwF2KpKed/5zwAmAf/6ccDMIoDJwFAgB1hsZsnOubRjfH8MsNr/dfExtkVEqsTSTXuZkLyGNbkHGJjQjIdGdKdbm8Zex/Jcee/hu9DM4ssMDwQy/e/0MbPZwCXAsco/B1/hr8D3rr/s9jGZ2ThgHEBcXFx5ooqIALDjYB5//TCDd5bl0KpxXZ4d2YcRvdsG3QJsVaUy5/yjgexS2znAIDNrDjwK9DWzPzjnHgfeBSaZ2YXAXHzn/EtvH5NzbgowBSApKclVIquIhInC4hJe/WYj//h0PflFxdw6pAPjz+pIZF1d31JawP9rOOd2A7eUGTsM3FTmqWW3RUQq5evMXUxMTmX9jkOc2bkFEy5OpH2Lhl7HqpEqU/65QGyp7Rj/mIhItcrdd5RHP0gjZfU2YpvV5+Xrkzi3W0ud4jmBypT/YqCTmSXgK/2RwOiApBIRKYe8wmJeXpjF5C8yAfjt0M6M+0V76tUOnQXYqkp5L/WcBQwBoswsB5jgnJtmZuOB+fgu9ZzunEutsqQiIn7OOT5L38HD89LYvOcIF/RozYMXdiPm5AZeRwsa5b3aZ9RxxlOAlIAmEhE5gQ27DvPw3FQWZOykQ4tIXh87iNM7RXkdK+jo428RCQqH84uYvCCTqYs2UKfWSTw4vBs3nBpPnVqhvQBbVVH5i0iN5pxj3qqtPPpBOtsO5HF5v2juH9aVlo3reR0tqKn8RaTGWrvtABOTU/k2aw/d2zZm8pi+9G/XzOtYIUHlLyI1zv6jhTzzyTpe+3YTjerV4i+X9mDUwDgiwnABtqqi8heRGqOkxDFnaQ5//Wgte44UMHpgHL8/rwsnR9bxOlrIUfmLSI2wMnsff05OZWX2Pvq3O5lXRwykR3QTr2OFLJW/iHhq96F8npyfwZtLsmkeWZe/X92by/pGa3ZuFVP5i4gniopLmPndZp7+OIMjBcX88vQE7jinE43qhe8a+9VJ5S8i1e67rN1MSE5l7baDnN4xiokjEunYspHXscKKyl9Eqs22/Xk8lpJO8sotRDetz4tj+jGsR2ud4vGAyl9Eqlx+UTHTv9rI85+vp6jEccc5nbj1zA7Ur6MF2Lyi8heRKvVFxg4empvGhl2HObdbK/58USJxzbUAm9dU/iJSJTbvPsLD89L4NH07CVGRzLhpAEO6tPQ6lvip/EUkoI4WFPPiF5m8tDCLWicZ9w3rys2nx1O3lk7x1CQqfxEJCOcc81O38ci8dHL3HWVE77Y8MLwbrZtoAbaaqFrK38zigOeAPcA659wTZnYGMMafIdE5d2p1ZBGRwMvccZCJyWl8lbmLrq0bMXvcYAa3b+51LDmBCpe/mU0HLgJ2OOd6lBofBjyL7+5eU51zTwA9gTnOudfN7E0A59wiYJGZXYrvlpAiEmQO5hXy3GfreeXrjTSoE8HEixO5dnA7akVojf2arjLv/GcAk4B//ThgZhHAZGAokAMsNrNk4FtgjpndDLxW5ueMBsZWIoeIVDPnHP9ensvjH65l16F8rkmK5Z7zu9C8YV2vo0k5Vbj8nXMLzSy+zPBAINM5lwVgZrOBS4BCfPf9XWhmc4BX/I/HAfudcwcrmkNEqtea3P1MSE5l6aa99I5tytTrk+gd29TrWPIzBfqcfzSQXWo7BxgEvARMNLPRwMZSj4/F/w/BsZjZOGAcQFxcXICjisjPsfdwAU99nMEb32+mWYM6/O2KXlzZP4aTtMZ+UKqWD3ydc2uAK48xPuEnvm8KMAUgKSnJVU06ETmR4hLHrO8389THGRzMK+KGU+K5e2hnmtTXAmzBLNDlnwvEltqO8Y+JSBBaumkPf34/ldQtBxiU0IyHLulO19aNvY4lARDo8l8MdDKzBHylPxLfB7oiEkR2HMzjiQ/X8u6yXFo3rsfzo/pyUa82WoAthFTmUs9ZwBAgysxy8H2gO83MxgPz8V3qOd05lxqQpCJS5QqLS3j1m43849P15BcVc9uQDtx+Vkci62o+aKipzNU+o44zngKkVDiRiHji68xdTEhOJXPHIYZ0acGEi7uTEBXpdSypIvrnXCTMZe85wqMfpPNR6jbimjVg6vVJnNOtpU7xhDiVv0iYyiss5qUvf+DFL37gJDN+f15nfnlGe+rV1gJs4UDlLxJmyi7AdlGvNjwwvBttm9b3OppUI5W/SBhZv/0gD831LcDWpVUjZv1qMKd00AJs4UjlLxIGDuQV8uyn63n1Gy3AJj4qf5EQVlLimLMsh799tJbdhwsYOSCW35+nBdhE5S8SslZm72NCciorsvfRN64p028cQK+Ypl7HkhpC5S8SYnYdyufJjzJ4a2k2zSPr8vRVvbmsb7QWYJP/ovIXCRGFxSW89p9NPPPpOo4WFPOrM9rzm7M70qieFmCT/0/lLxICvsncxcS5qazbfogzOkUx4eLudGzZ0OtYUoOp/EWCWO6+ozz6QRopq7cR26w+U67rz9DEVpqdKz9J5S8ShPIKi/nnl1m8+GUmAL8d2plxv9DsXCk/lb9IEHHO8XHadh6Zl0bO3qNc2LMNfxjelZiTG3gdTYKMyl8kSGTuOMRDc1NZtH4XnVs15I1fDuLUjlFex5IgpfIXqeEO5hXy3GfreeXrjdSvE8EE/+zc2pqdK5Wg8hepoUpKHP9enssTH61l16F8ru4fyz3DuhCl2bkSANVS/mZ2BjDG/3qJzrlT/eORwJfAROfcvOrIIhIMVufsZ0LyGpZt3kef2KZMvT6J3rFNvY4lIaQyt3GcDlwE7HDO9Sg1Pgx4Ft9tHKc6555wzi0CFpnZpfju8/uj+4C3KppBJNTsPpTPUx9nMHtxNs0j6/Dklb24ol+MZudKwFXmnf8MYBLwrx8HzCwCmAwMBXKAxWaW7JxL8z9lNDDW/9yhQBpQrxIZREJCUXEJr3+7ib9/so4jBcWMPS2BO87tRGPNzpUqUpl7+C40s/gywwOBTOdcFoCZzQYuAdLMLA7Y75w76H/uECASSASOmlmKc66konlEgtV/ftjNxORUMrYf5PSOUUwckUjHlo28jiUhLtDn/KOB7FLbOcAg/9djgVd+fMA59yCAmd0I7DpW8ZvZOGAcQFxcXICjinhry76jPJqSzgerthLdtD4vXduP87u31uxcqRbVdrWPc27CccZnnOB7pgBTAJKSklzVJBOpXnmFxby8MIvJX2TiHNx1biduObODZudKtQp0+ecCsaW2Y/xjImHPOcen6Tt4ZF4am/cc4YIerXlgeDdim2l2rlS/QJf/YqCTmSXgK/2R+D7kFQlrP+w8xMNz0/hy3U46tmzI62MHcXonzc4V71TmUs9Z+D60jTKzHGCCc26amY0H5uO71HO6cy41IElFgtCh/CKe/2w907/eQL1aEfzpokSuP0Wzc8V7lbnaZ9RxxlOAlAonEgkBzvlm5z7+4Vp2Hsznqv4x3DusKy0aaXau1Axa3kEkwNbk7mdCcipLN+2ld0wTXr4+iT6anSs1jMpfJED2HC7gyfkZzF68mWYN6vC3K3pxZX/NzpWaSeUvUklFxSW88f1mnv54HYfyi7jp1ATuPLcTTeprdq7UXCp/kUr4Lms3E5JTWbvtIKd2aM7EEd3p3Eqzc6XmU/mLVMDW/Ud5LGUtc1duIbppfV4c049hPTQ7V4KHyl/kZ8gvKmbqog1M+jyTYue445xO3HpmB+rX0excCS4qf5Fy+ix9Ow/PS2PT7iOc370Vf7wwUbNzJWip/EV+QtbOQzwyL40FGTvp0CKS18YO5IxOLbyOJVIpKn+R4ziUX8SkzzOZ9lUWdWtF8McLu3HDqfGanSshQeUvUoZzjvdXbOHxD9PZfiCfK/vHcO+wLrRspPsOSehQ+YuUsiZ3PxOTU1myaS+9Yprw4rX96Rd3stexRAJO5S8C7D1cwFMfZzDr+800bVCHJy7vydVJsZqdKyFL5S9hrbjE+WfnZnAwr4jrT4nn7nM706SBZudKaFP5S9j6fsMeJiSnkr71AIPbN2PiiO50bd3Y61gi1ULlL2Fn2/48Hv8wnfdXbKFtk3pMHt2P4T01O1fCi8pfwkZ+UTHTvvLNzi0qcdxxdkduGdKBBnX010DCT7X8qTezk4BHgMbAEufcq2bWHngQaOKcu7I6ckj4WrB2Bw/PS2PDrsMMTWzFny5MJK65ZudK+KrwbBUzm25mO8xsTZnxYWaWYWaZZna/f/gSfDdzLwRyAJxzWc65sRV9fZHy2LjrMDfPWMxNMxZjBq/ePJCXr09S8UvYq8w7/xnAJOBfPw6YWQQwGRiKr+QXm1ky0AX4xjn3TzObA3xWidcV+UmH84uYvCCTqYs2UDvCeGB4V248NYE6tTQ7VwQqdw/fhWYWX2Z4IJDpnMsCMLPZ+N71ZwMF/ucUl/c1zGwcMA4gLi6uolEljDjnSF65hcdT1rLtQB6X943m/gu60rKxZueKlBboc/7R+Ir+RznAIOBZ4HkzOwNYCGBmzYFHgb5m9gfn3ONlf5hzbgowBSApKckFOKuEmLQtB5g4N5XvN+yhR3RjJo/pS/92zbyOJVIjVcsHvs65I8DYMmO7gVuq4/UltO07UsDTH69j5nebaFK/No9d1pNrBsQSodm5IscV6PLPBWJLbcf4x0QCrrjEMXvxZp6an8H+o4VcN7gdvx3aRbNzRcoh0OW/GOhkZgn4Sn8kMDrAryHCko2+2bmpWw4wKME3O7dbG83OFSmvCpe/mc0ChgBRZpYDTHDOTTOz8cB8IAKY7pxLDUhSEWD7gTye+HAt/16eS5sm9Xh+VF8u6tVGs3NFfqbKXO0z6jjjKUBKhROJHENBUQnTv97A85+tp7DYMf6sjtx2lmbnilSU/uZIjbcgYwePzE0ja9dhzu3Wkj9dlEi75pFexxIJaip/qbE27T7MI/PS+DR9BwlRkbxy0wDO6tLS61giIUHlLzXOkQLf7NyXF/pm595/QVduPk2zc0UCSeUvNYZzjnmrtvJYSjpb9+dxmX92bivNzhUJOJW/1Ahrtx1gYnIq32btIbFNY54f1ZekeM3OFakqKn/x1P4jhfz9kwxe+3YTjevX5tHLejByQJxm54pUMZW/eKK4xPHWkmyenJ/BviMFjBnUjt+d15mmDep4HU0kLKj8pdpl7TzEXW+uYFXOfgbG+2bnJrbV7FyR6qTyl2qVsnor985ZRe0I49mRfRjRu61m54p4QOUv1aKwuITHU9Yy/esN9I1ryuTR/WjbtL7XsUTClspfqtzW/UcZ/8Zylm7ay42nxvPA8G66Zl/EYyp/qVJfrd/FnbOXk1dYzKTRfbmoV1uvI4kIKn+pIiUljkkLMnnm03V0atmQF8b0p2PLhl7HEhE/lb8E3N7DBdz91gq+yNjJpX3a8tjlPbX6pkgNo7+RElArsvdx+8xl7DyYz18u7cGYQXG6mkekBlL5S0A453j92008PC+Nlo3qMefWU+gV09TrWCJyHNVW/mYWCXwJTHTOzTOzOOA5YA+wzjn3RHVlkcA6nF/EA/9ezfsrtnBWlxY8c00fzdQVqeEqfL2dmU03sx1mtqbM+DAzyzCzTDO7v9RD9wFvldruCcxxzt0M9K1oDvFW5o6DXDr5a+au3MI953dh2g0DVPwiQaAy7/xnAJOAf/04YGYRwGRgKJADLDazZCAaSANKr837LTDHzG4GXqtEDvFI8sot3P/OKurXjuC1sYM4rWOU15FEpJwqcw/fhWYWX2Z4IJDpnMsCMLPZwCVAQyASSASOmlkKcBO+m74vNLM5wCtlX8PMxgHjAOLi4ioaVQIsv6iYxz5I59X/bCKp3clMGt2P1k205r5IMAn0Of9oILvUdg4wyDk3HsDMbgR2OedKzOwjYKKZjQY2HuuHOeemAFMAkpKSXICzSgXk7jvKbTOXsTJ7H788PYH7LuhK7QjN1hUJNtV6tY9zbkapr9cAV1bn60vlfLluJ3fNXk5hseOla/sxrEcbryOJSAUFuvxzgdhS2zH+MQlixSWOZz9bz/Ofr6dLq0a8eG1/EqIivY4lIpUQ6PJfDHQyswR8pT8SGB3g15BqtPtQPne9uYJF63dxRb8Y/nJpD+rXifA6lohUUoXL38xmAUOAKDPLwffh7TQzGw/MByKA6c651IAklWq3dNNexr+xjN2HC3ji8p5cMyBWs3VFQkRlrvYZdZzxFCClwonEc845ZnyzkUc/SKdN03q8e+up9Ihu4nUsEQkgLe8g/+VQfhH3vbOKD1Zt5dxurXj6qt40aVDb61giEmAqf/lf67Yf5JbXl7Jx12Huv6Arv/5Fe53mEQlRKn8B4N/Lc3jg3TVE1q3FG78azOD2zb2OJCJVSOUf5vIKi3l4XhpvfLeZgQnNmDSqLy0ba7auSKhT+Yex7D1HuG3mMlbn7ufXZ7bnnvO6UEuzdUXCgso/TH2+djt3v7mSEueYcl1/zuve2utIIlKNVP5hprjE8fdPMpi84AcS2zTmxWv70a65ZuuKhBuVfxjZeTCfO2cv55sfdjNyQCwTR3SnXm3N1hUJRyr/MLF44x5un7mM/UcLefLKXlyVFPvT3yQiIUvlH+Kcc0xdtIEnPlpL7Mn1efXmgXRr09jrWCLiMZV/CDuQV8g9b69kfup2hnVvzd+u6kXjepqtKyIq/5CVtuUAt81cSvbeo/zxwm6MPT1Bs3VF5H+p/EPQ20uy+eN7a2jaoDazxw1mQHwzryOJSA2j8g8heYXFTExOZfbibE7t0JxnR/alRaO6XscSkRpI5R8iNu0+zK2vLyNt6wHGn9WRu4d2JuIkneYRkWNT+YeAj1O38bu3V3KSGdNvTOLsrq28jiQiNVy1lL+ZdQPuBKKAz5xzL5rZEOARIBWY7Zz7ojqyhJKi4hKenJ/BPxdm0TO6CS+M6UdsswZexxKRIFDhVbzMbLqZ7TCzNWXGh5lZhpllmtn9AM65dOfcLcDVwGn+pzrgEFAPyKlojnC140Aeo6d+xz8XZjFmUBxv33KKil9Eyq0ySzjOAIaVHjCzCGAycAGQCIwys0T/YyOAD/i/Wzwucs5dANwHPFSJHGHnPz/sZvhzX7E6Zz/PXNObRy/rqWUaRORnqXD5O+cWAnvKDA8EMp1zWc65AmA2cIn/+cn+sh/j3y7xf89e4JiXpJjZODNbYmZLdu7cWdGoIaOkxPHiFz8wZuq3NK5fi/duP43L+sZ4HUtEglCgz/lHA9mltnOAQf7z+5fjK/kUADO7HDgfaApMOtYPc85NAaYAJCUluQBnDSr7jxTyu7dX8mn6di7s1Ya/XtGLhnX1eb2IVEy1tIf/w9wvyoy9C7xbHa8f7Nbk7ufWmUvZui+PCRcncuOp8ZqtKyKVEujyzwVKLxcZ4x+TCnDOMXtxNhOSU2keWYc3f30K/dud7HUsEQkBgS7/xUAnM0vAV/ojgdEBfo2wcLSgmD++t4Z3luVwRqco/nFNH5o31GxdEQmMCpe/mc0ChgBRZpYDTHDOTTOz8cB8IAKY7pxLDUjSMJK18xC3zVxGxvaD3HlOJ+44p5Nm64pIQFW4/J1zo44znsL/Xc4pP9OHq7dyz5xV1I4wXrlxAEO6tPQ6koiEIF0uUkMUFpfwxIdrmfbVBvrENmXymH5EN63vdSwRCVEq/xpg2/48bn9jGUs37eWGU9rx4IWJ1KlVmfl3IiInpvL32NeZu7hj1nKOFhbz3Ki+jOjd1utIIhIGVP4eKSlxvPBFJn//ZB3tWzTkzWv70bFlI69jiUiYUPl7YN+RAu5+cwULMnZySZ+2PHZZTyI1W1dEqpEap5qtytnHra8vY+fBfB65tAfXDorTbF0RqXYq/2rinOP17zbzyNw0WjSqy9u3nELv2KZexxKRMKXyrwZHCop44N3VvLdiC0O6tOCZq/twcmQdr2OJSBhT+VexzB2HuPX1pWTuPMTvhnbm9rM6cpJm64qIx1T+VWjuyi3c984q6teO4LWbB3F6pyivI4mIACr/KlFQVMJjKenM+GYj/dudzKTRfWnTRLN1RaTmUPkHWO6+o9w+cxkrsvcx9vQE7r+gK7UjNFtXRGoWlX8AfbluJ3fNXk5hseOFMf0Y3rON15FERI5J5R8AxSWO5z5bz3Ofr6dLq0a8MKYf7Vs09DqWiMhxqfwrac/hAu6cvZxF63dxeb9oHr20J/XrRHgdS0TkhFT+lbBs815un7mM3YcLePzynowcEKvZuiISFKql/M3sUuBCoDEwzTn38bHGqiNLIDjnePWbjTyakk7rJvV455ZT6RnTxOtYIiLlVuHLUMxsupntMLM1ZcaHmVmGmWWa2f0Azrn3nHO/Am4BrjneWDA4lF/Eb2YtZ+LcNM7s3IJ5489Q8YtI0KnMNYgzgGGlB8wsApgMXAAkAqPMLLHUU/7of5yfGKuR1m0/yIhJX5Gyeiv3DuvClOuSaNKgttexRER+tsrcw3ehmcWXGR4IZDrnsgDMbDZwiZmlA08AHzrnlvkfs7JjZZnZOGAcQFxcXEWjBsR7y3P5w7uriaxbi5m/HMwpHZp7mkdEpDICfc4/GsgutZ0DDAJ+A5wLNDGzjs65l44z9l+cc1OAKQBJSUkuwFnLJb+omIfnpjHzu80MjG/GpNF9adm4nhdRREQCplo+8HXOPQc891NjNU32niPc/sYyVuXs59e/aM8953ehlmbrikgICHT55wKxpbZj/GNB5/O127n7zZWUlDj+eV1/zu/e2utIIiIBE+jyXwx0MrMEfKU/Ehgd4NeoUsUljmc+WcekBZl0a9OYl67tR7vmkV7HEhEJqAqXv5nNAoYAUWaWA0xwzk0zs/HAfCACmO6cSw1I0mqw61A+d8xazjc/7ObqpBgevqQH9Wprtq6IhJ7KXO0z6jjjKUBKhRN5ZMnGPdz+xjL2HSnkb1f04uoBsT/9TSIiQSrsl3dwzjHtqw08/uFaYk6uz7u3DaB7W03aEpHQFtblfyCvkHvfXsVHqds4L7EVT13dm8b1NGlLREJf2JZ/+tYD3Pr6UrL3HuXB4d345RkJWpRNRMJGWJb/nKU5/PG91TSuV5tZvxrMwIRmXkcSEalWYVX+eYXFTExOZfbibE5p35xnR/WhZSPN1hWR8BM25b959xFunbmU1C0HuG1IB347tLNm64pI2AqL8v8kbTu/fWsFBky7IYlzurXyOpKIiKdCvvyf/XQ9z3y6jh7RjXlxTH9imzXwOpKIiOdCvvzjoxowelAcf74oUbN1RUT8Qr78L+kTzSV9or2OISJSo+gTTxGRMKTyFxEJQyp/EZEwpPIXEQlDKn8RkTCk8hcRCUMqfxGRMKTyFxEJQ+ac8zpDuZjZTmBTBb89CtgVwDhe0r7UPKGyH6B9qakqsy/tnHMtyg4GTflXhpktcc4leZ0jELQvNU+o7AdoX2qqqtgXnfYREQlDKn8RkTAULuU/xesAAaR9qXlCZT9A+1JTBXxfwuKcv4iI/LdweecvIiKlqPxFRMJQyJW/mU03sx1mtqbUWDMz+8TM1vt/P9nLjOV1nH2ZaGa5ZrbC/2u4lxnLw8xizWyBmaWZWaqZ3ekfD7rjcoJ9CcbjUs/Mvjezlf59ecg/nmBm35lZppm9aWZ1vM56IifYjxlmtqHUMenjcdRyM7MIM1tuZvP82wE/JiFX/sAMYFiZsfuBz5xznYDP/NvBYAb/f18AnnHO9fH/SqnmTBVRBPzOOZcIDAZuN7NEgvO4HG9fIPiOSz5wtnOuN9AHGGZmg4G/4tuXjsBeYKx3EcvlePsBcE+pY7LCq4AVcCeQXmo74Mck5MrfObcQ2FNm+BLgVf/XrwKXVmemijrOvgQd59xW59wy/9cH8f2hjiYIj8sJ9iXoOJ9D/s3a/l8OOBuY4x+v8cflBPsRlMwsBrgQmOrfNqrgmIRc+R9HK+fcVv/X24BWXoYJgPFmtsp/WqjGnyopzczigb7AdwT5cSmzLxCEx8V/emEFsAP4BPgB2OecK/I/JYcg+Met7H445348Jo/6j8kzZlbXu4Q/yz+Ae4ES/3ZzquCYhEv5/y/nu7Y1aN8VAC8CHfD97+1W4GlP0/wMZtYQeAe4yzl3oPRjwXZcjrEvQXlcnHPFzrk+QAwwEOjqbaKKKbsfZtYD+AO+/RkANAPu8y5h+ZjZRcAO59zSqn6tcCn/7WbWBsD/+w6P81SYc267/w96CfAyvr+wNZ6Z1cZXljOdc+/6h4PyuBxrX4L1uPzIObcPWACcAjQ1s1r+h2KAXK9y/Vyl9mOY/xSdc87lA68QHMfkNGCEmW0EZuM73fMsVXBMwqX8k4Eb/F/fALzvYZZK+bEs/S4D1hzvuTWF/5zlNCDdOff3Ug8F3XE53r4E6XFpYWZN/V/XB4bi+wxjAXCl/2k1/rgcZz/WlnpjYfjOkdf4Y+Kc+4NzLsY5Fw+MBD53zo2hCo5JyM3wNbNZwBB8S6BuByYA7wFvAXH4loW+2jlX4z9IPc6+DMF3asEBG4FflzpvXiOZ2enAImA1/3ce8wF858qD6ricYF9GEXzHpRe+Dw8j8L0RfMs597CZtcf3rrMZsBy41v/uuUY6wX58DrQADFgB3FLqg+Eaz8yGAL93zl1UFcck5MpfRER+Wric9hERkVJU/iIiYUjlLyIShlT+IiJhSOUvIhKGVP4iImFI5S8iEob+BwV7soPIh2Z9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def error_inv(n):\n",
    "    return np.mean((Hilbert_inv(n) - np.linalg.inv(Hilbert(n)))**2)\n",
    "\n",
    "fig, ax =  plt.subplots()\n",
    "xs = np.arange(10, 50, 10)\n",
    "ys = [error_inv(x) for x in xs]\n",
    "\n",
    "plt.plot(xs, ys)\n",
    "ax.set_yscale(\"log\", nonpositive='clip')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d0d459d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}