lectures/lecture8/lecture-8.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 8: Symmetric eigenvalue problem and SVD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Recap of the first part\n",
    "\n",
    "- QR decomposition and Gram-Schmidt algorithm\n",
    "- Schur decomposition and QR-algorithm (basic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Plan for today\n",
    "\n",
    "Today we will talk about:\n",
    "\n",
    "- Algorithms for the symmetric eigenvalue problems\n",
    "    - QR algorithm (in more details)\n",
    "    - Divide-and-Conquer\n",
    "    - bisection\n",
    "- Algorithms for SVD computation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Schur form computation\n",
    "\n",
    "- Recall that we are trying to avoid $\\mathcal{O}(n^4)$ complexity. \n",
    "- The idea is to make a matrix have a simpler structure so that each step of QR algorithm becomes cheaper.\n",
    "\n",
    "- In case of a general matrix we can use the **Hessenberg form**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Hessenberg form\n",
    "\n",
    "The matrix $A$ is said to be in the Hessenberg form, if\n",
    "\n",
    "$$a_{ij} = 0, \\quad \\mbox{if } i \\geq j+2.$$\n",
    "\n",
    "$$H = \\begin{bmatrix} * & * & * & * & * \\\\ * & * & * & * & * \\\\ 0 & * & * & * & *\\\\ 0 & 0 & * & * & *\\\\ 0 & 0 & 0 & * & * \\\\ \\end{bmatrix}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Reduction any matrix to Hessenberg form\n",
    "\n",
    "- By applying Householder reflections we can reduce any matrix to the Hessenberg form\n",
    "\n",
    "$$U^* A U = H$$\n",
    "\n",
    "- The only difference with Schur decomposition is that we have to map the first column to the vector with two non-zeros, and the first element is not changed.\n",
    "\n",
    "- The computational cost of such reduction is $\\mathcal{O}(n^3)$ operations.\n",
    "\n",
    "- In a Hessenberg form, computation of one iteration of the QR algorithm costs $\\mathcal{O}(n^2)$ operations (e.g. using Givens rotations, how?), and the Hessenberg form is preserved by the QR iteration (check why)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Symmetric (Hermitian) case\n",
    "\n",
    "- In the symmetric case, we have $A = A^*$, then $H = H^*$ and the upper Hessenberg form becomes tridiagonal matrix.\n",
    "\n",
    "- From now on we will talk about the case of symmetric tridiagonal form.\n",
    "\n",
    "- Any symmetric (Hermitian) matrix can be reduced to the tridiagonal form by Householder reflections.\n",
    "\n",
    "- **Key point** is that tridiagonal form is preserved by the QR algorithm, and the cost of one step can be reduced to $\\mathcal{O}(n)$!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR algorithm: iterations\n",
    "\n",
    "- The iterations of the QR algorithm have the following form:\n",
    "\n",
    "$$A_k = Q_k R_k, \\quad A_{k+1} = R_k Q_k.$$\n",
    "\n",
    "- If $A_0 = A$ is <font color='red'> tridiagonal symmetric matrix </font>, this form is preserved by the QR algorithm.\n",
    "\n",
    "Let us see.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "DeviceArray([-1.47168069, -0.75157939,  0.        ,  0.        ,\n",
       "              0.        ,  0.        ,  0.        ,  0.        ,\n",
       "              0.        ,  0.        ,  0.        ,  0.        ,\n",
       "              0.        ,  0.        ,  0.        ,  0.        ,\n",
       "              0.        ,  0.        ,  0.        ,  0.        ],            dtype=float64)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD8CAYAAACxd9IeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAIxUlEQVR4nO3dwYuU9x3H8c+nCb2EHAyOIm0tJUghly7dxUugGEqCzcXkUKgnD4HNof4B3pJjLiGnEthQ0UssvUg8hDTixUsP3YCkW2gxhG1rFHfFS28l5tuDj7DYXXec55l5Zp7P+wXLzDyOM7/H5e0zs359xlUlAMP3vb4XAGA2iB0IQexACGIHQhA7EILYgRC9xm77pO1/2P7K9rk+19Il25u2/2r7hu31vtfThu3ztrdsb+zY9oLtq7ZvNpcH+lzjJPbYr3dtf9N8327Yfr3PNXatt9htPyPpd5J+JeklSadtv9TXeqbglapaqqqVvhfS0gVJJx/bdk7Stao6Julac3vRXND/75ckfdB835aq6tMZr2mq+jyyH5f0VVV9XVX/lfQHSad6XA92UVXXJd1/bPMpSReb6xclvTHLNXVhj/0atD5j/4Gkf++4favZNgQl6XPbX9he7XsxU3C4qu5IUnN5qOf1dOms7S+bl/kL9/bkSfqM3btsG8rs7stV9XM9fIvyW9u/6HtBGMuHkl6UtCTpjqT3e11Nx/qM/ZakH+24/UNJt3taS6eq6nZzuSXpsh6+ZRmSu7aPSFJzudXzejpRVXer6kFVfSfpIw3s+9Zn7H+RdMz2T2x/X9JvJF3pcT2dsP2c7ecfXZf0mqSNJ/+uhXNF0pnm+hlJn/S4ls48+gus8aYG9n17tq8nrqpvbZ+V9CdJz0g6X1V/62s9HTos6bJt6eGf78dV9Vm/S5qc7UuSTkg6aPuWpHckvSfpj7bfkvQvSb/ub4WT2WO/Tthe0sO3k5uS3u5rfdNg/osrkIEJOiAEsQMhiB0IQexACGIHQvQe+0DHSSUNd9/Yr8XUe+yShvwHPNR9Y78W0DzEDmAGZjpUY3sqT7a8vDyNh21te3tbo9Go72V0jv2aX5ubm7p3795u/8msv3HZLq2vL/TJYIDOrKzsfa6UVi/jh3paKWCIJo494LRSwKC0ObJzWilggbSJfazTStletb2+6GdZBRZdmx/QjXVaqapak7QmTe+n8QD21+bIPtjTSgFD1Cb2QZ5WChiqiV/GD/i0UsAgtRqqaT4xo/dPzWjO9zYWTsOFVMzGAyGIHQhB7EAIYgdCEDsQgtiBEMQOhCB2IASxAyGIHQgx09iXl5dVVWN9TYvtsb+AIeHIDoQgdiAEsQMhiB0IQexACGIHQhA7EILYgRDEDoQgdiDE3H5k89OMzE5rtHXcx+WMtVgEHNmBEMQOhCB2IASxAyGIHQhB7EAIYgdCEDsQgtiBEMQOhJjbcdmn0fdo7dM8JqO16AtHdiAEsQMhiB0IQexACGIHQhA7EILYgRDEDoQgdiDEICbongbTdkjFkR0IQexAiFYv421vSvqPpAeSvq2qlS4WBaB7Xbxnf6Wq7nXwOACmiJfxQIi2sZekz21/YXu1iwUBmI62L+Nfrqrbtg9Jumr771V1fecdmr8EViXp6NGjLZ8OwKRaHdmr6nZzuSXpsqTju9xnrapWqmplNBq1eToALUwcu+3nbD//6Lqk1yRtdLUwAN1q8zL+sKTLzUTYs5I+rqrPOlkVgM5NHHtVfS3pZx2uZe6MO67a9+fDS4zWYn/80xsQgtiBEMQOhCB2IASxAyGIHQhB7EAIYgdCEDsQgtiBEHFnl52Gvs9Y+7SPy2htJo7sQAhiB0IQOxCC2IEQxA6EIHYgBLEDIYgdCEHsQAhiB0IwLjtjjNaiLxzZgRDEDoQgdiAEsQMhiB0IQexACGIHQhA7EILYgRDEDoRgXHaOMVqLLnFkB0IQOxCC2IEQxA6EIHYgBLEDIYgdCEHsQAhiB0IQOxCCcdmBYLQW++HIDoTYN3bb521v2d7Yse0F21dt32wuD0x3mQDaGufIfkHSyce2nZN0raqOSbrW3AYwx/aNvaquS7r/2OZTki421y9KeqPbZQHo2qTv2Q9X1R1Jai4PdbckANMw9R/Q2V61vW57fXt7e9pPB2APk8Z+1/YRSWout/a6Y1WtVdVKVa2MRqMJnw5AW5PGfkXSmeb6GUmfdLMcANMyzj+9XZL0Z0k/tX3L9luS3pP0qu2bkl5tbgOYY/tO0FXV6T1+6ZcdrwXAFDEuG4jR2kyMywIhiB0IQexACGIHQhA7EILYgRDEDoQgdiAEsQMhiB0IwbgsnojR2uHgyA6EIHYgBLEDIYgdCEHsQAhiB0IQOxCC2IEQxA6EIHYgBOOy6MwijdYmjtVyZAdCEDsQgtiBEMQOhCB2IASxAyGIHQhB7EAIYgdCEDsQgnFZ9KLv0drEM9ZyZAdCEDsQgtiBEMQOhCB2IASxAyGIHQhB7EAIYgdCMEGHuce0XTc4sgMhiB0IsW/sts/b3rK9sWPbu7a/sX2j+Xp9ussE0NY4R/YLkk7usv2Dqlpqvj7tdlkAurZv7FV1XdL9GawFwBS1ec9+1vaXzcv8A3vdyfaq7XXb69vb2y2eDkAbk8b+oaQXJS1JuiPp/b3uWFVrVbVSVSuj0WjCpwPQ1kSxV9XdqnpQVd9J+kjS8W6XBaBrE8Vu+8iOm29K2tjrvgDmw74TdLYvSToh6aDtW5LekXTC9pKkkrQp6e3pLRFAF/aNvapO77L591NYC9DauOOqfX8+vDT70Vom6IAQxA6EIHYgBLEDIYgdCEHsQAhiB0IQOxCC2IEQxA6E4OyyiNT3GWuf9nG7GK3lyA6EIHYgBLEDIYgdCEHsQAhiB0IQOxCC2IEQxA6EIHYgBOOywD4WbbR2LxzZgRDEDoQgdiAEsQMhiB0IQexACGIHQhA7EILYgRDEDoRgXBbo0DyM1u6FIzsQgtiBEMQOhCB2IASxAyGIHQhB7EAIYgdCEDsQgtiBEO7iQ97HfjJ7W9I/H9t8UNK9mS1itoa6b+zX/PpxVY12+4WZxr7rAuz1qlrpdRFTMtR9Y78WEy/jgRDEDoSYh9jX+l7AFA1139ivBdT7e3YAszEPR3YAM0DsQAhiB0IQOxCC2IEQ/wNApy3htX+LEgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD8CAYAAACxd9IeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAIxUlEQVR4nO3dwYuU9x3H8c+nCb2EHAyOIm0tJUghly7dxUugGEqCzcXkUKgnD4HNof4B3pJjLiGnEthQ0UssvUg8hDTixUsP3YCkW2gxhG1rFHfFS28l5tuDj7DYXXec55l5Zp7P+wXLzDyOM7/H5e0zs359xlUlAMP3vb4XAGA2iB0IQexACGIHQhA7EILYgRC9xm77pO1/2P7K9rk+19Il25u2/2r7hu31vtfThu3ztrdsb+zY9oLtq7ZvNpcH+lzjJPbYr3dtf9N8327Yfr3PNXatt9htPyPpd5J+JeklSadtv9TXeqbglapaqqqVvhfS0gVJJx/bdk7Stao6Julac3vRXND/75ckfdB835aq6tMZr2mq+jyyH5f0VVV9XVX/lfQHSad6XA92UVXXJd1/bPMpSReb6xclvTHLNXVhj/0atD5j/4Gkf++4favZNgQl6XPbX9he7XsxU3C4qu5IUnN5qOf1dOms7S+bl/kL9/bkSfqM3btsG8rs7stV9XM9fIvyW9u/6HtBGMuHkl6UtCTpjqT3e11Nx/qM/ZakH+24/UNJt3taS6eq6nZzuSXpsh6+ZRmSu7aPSFJzudXzejpRVXer6kFVfSfpIw3s+9Zn7H+RdMz2T2x/X9JvJF3pcT2dsP2c7ecfXZf0mqSNJ/+uhXNF0pnm+hlJn/S4ls48+gus8aYG9n17tq8nrqpvbZ+V9CdJz0g6X1V/62s9HTos6bJt6eGf78dV9Vm/S5qc7UuSTkg6aPuWpHckvSfpj7bfkvQvSb/ub4WT2WO/Tthe0sO3k5uS3u5rfdNg/osrkIEJOiAEsQMhiB0IQexACGIHQvQe+0DHSSUNd9/Yr8XUe+yShvwHPNR9Y78W0DzEDmAGZjpUY3sqT7a8vDyNh21te3tbo9Go72V0jv2aX5ubm7p3795u/8msv3HZLq2vL/TJYIDOrKzsfa6UVi/jh3paKWCIJo494LRSwKC0ObJzWilggbSJfazTStletb2+6GdZBRZdmx/QjXVaqapak7QmTe+n8QD21+bIPtjTSgFD1Cb2QZ5WChiqiV/GD/i0UsAgtRqqaT4xo/dPzWjO9zYWTsOFVMzGAyGIHQhB7EAIYgdCEDsQgtiBEMQOhCB2IASxAyGIHQgx09iXl5dVVWN9TYvtsb+AIeHIDoQgdiAEsQMhiB0IQexACGIHQhA7EILYgRDEDoQgdiDE3H5k89OMzE5rtHXcx+WMtVgEHNmBEMQOhCB2IASxAyGIHQhB7EAIYgdCEDsQgtiBEMQOhJjbcdmn0fdo7dM8JqO16AtHdiAEsQMhiB0IQexACGIHQhA7EILYgRDEDoQgdiDEICbongbTdkjFkR0IQexAiFYv421vSvqPpAeSvq2qlS4WBaB7Xbxnf6Wq7nXwOACmiJfxQIi2sZekz21/YXu1iwUBmI62L+Nfrqrbtg9Jumr771V1fecdmr8EViXp6NGjLZ8OwKRaHdmr6nZzuSXpsqTju9xnrapWqmplNBq1eToALUwcu+3nbD//6Lqk1yRtdLUwAN1q8zL+sKTLzUTYs5I+rqrPOlkVgM5NHHtVfS3pZx2uZe6MO67a9+fDS4zWYn/80xsQgtiBEMQOhCB2IASxAyGIHQhB7EAIYgdCEDsQgtiBEHFnl52Gvs9Y+7SPy2htJo7sQAhiB0IQOxCC2IEQxA6EIHYgBLEDIYgdCEHsQAhiB0IwLjtjjNaiLxzZgRDEDoQgdiAEsQMhiB0IQexACGIHQhA7EILYgRDEDoRgXHaOMVqLLnFkB0IQOxCC2IEQxA6EIHYgBLEDIYgdCEHsQAhiB0IQOxCCcdmBYLQW++HIDoTYN3bb521v2d7Yse0F21dt32wuD0x3mQDaGufIfkHSyce2nZN0raqOSbrW3AYwx/aNvaquS7r/2OZTki421y9KeqPbZQHo2qTv2Q9X1R1Jai4PdbckANMw9R/Q2V61vW57fXt7e9pPB2APk8Z+1/YRSWout/a6Y1WtVdVKVa2MRqMJnw5AW5PGfkXSmeb6GUmfdLMcANMyzj+9XZL0Z0k/tX3L9luS3pP0qu2bkl5tbgOYY/tO0FXV6T1+6ZcdrwXAFDEuG4jR2kyMywIhiB0IQexACGIHQhA7EILYgRDEDoQgdiAEsQMhiB0IwbgsnojR2uHgyA6EIHYgBLEDIYgdCEHsQAhiB0IQOxCC2IEQxA6EIHYgBOOy6MwijdYmjtVyZAdCEDsQgtiBEMQOhCB2IASxAyGIHQhB7EAIYgdCEDsQgnFZ9KLv0drEM9ZyZAdCEDsQgtiBEMQOhCB2IASxAyGIHQhB7EAIYgdCMEGHuce0XTc4sgMhiB0IsW/sts/b3rK9sWPbu7a/sX2j+Xp9ussE0NY4R/YLkk7usv2Dqlpqvj7tdlkAurZv7FV1XdL9GawFwBS1ec9+1vaXzcv8A3vdyfaq7XXb69vb2y2eDkAbk8b+oaQXJS1JuiPp/b3uWFVrVbVSVSuj0WjCpwPQ1kSxV9XdqnpQVd9J+kjS8W6XBaBrE8Vu+8iOm29K2tjrvgDmw74TdLYvSToh6aDtW5LekXTC9pKkkrQp6e3pLRFAF/aNvapO77L591NYC9DauOOqfX8+vDT70Vom6IAQxA6EIHYgBLEDIYgdCEHsQAhiB0IQOxCC2IEQxA6E4OyyiNT3GWuf9nG7GK3lyA6EIHYgBLEDIYgdCEHsQAhiB0IQOxCC2IEQxA6EIHYgBOOywD4WbbR2LxzZgRDEDoQgdiAEsQMhiB0IQexACGIHQhA7EILYgRDEDoRgXBbo0DyM1u6FIzsQgtiBEMQOhCB2IASxAyGIHQhB7EAIYgdCEDsQgtiBEO7iQ97HfjJ7W9I/H9t8UNK9mS1itoa6b+zX/PpxVY12+4WZxr7rAuz1qlrpdRFTMtR9Y78WEy/jgRDEDoSYh9jX+l7AFA1139ivBdT7e3YAszEPR3YAM0DsQAhiB0IQOxCC2IEQ/wNApy3htX+LEgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import jax.numpy as jnp\n",
    "import jax\n",
    "import matplotlib.pyplot as plt\n",
    "from jax.config import config\n",
    "config.update(\"jax_enable_x64\", True)\n",
    "#Generate a random tridiagonal matrix\n",
    "\n",
    "n = 20\n",
    "d = jax.random.normal(jax.random.PRNGKey(0), (n, ))\n",
    "sub_diag = jax.random.normal(jax.random.PRNGKey(1), (n-1,))\n",
    "\n",
    "mat = jnp.diag(d) + jnp.diag(sub_diag, -1) + jnp.diag(sub_diag, 1)\n",
    "mat1 = jnp.abs(mat)\n",
    "mat1 = mat1/jnp.max(mat1.flatten())\n",
    "plt.spy(mat)\n",
    "q, r = jnp.linalg.qr(mat)\n",
    "plt.figure()\n",
    "b = r.dot(q)\n",
    "# b[abs(b) <= 1e-12] = 0\n",
    "b = jax.ops.index_update(b, jax.ops.index[abs(b) <= 1e-12], 0)\n",
    "plt.spy(b)\n",
    "#plt.figure()\n",
    "#plt.imshow(np.abs(r.dot(q)))\n",
    "b[0, :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Tridiagonal form\n",
    "\n",
    "- In the tridiagonal form, you do not have to compute the $Q$ matrix: you only have to compute the **triadiagonal part** that appears after the iterations \n",
    "\n",
    "$$A_k = Q_k R_k, \\quad A_{k+1}  = R_k Q_k,$$\n",
    "\n",
    "in the case when $A_k = A^*_k$ and is also tridiagonal.\n",
    "\n",
    "- Such matrix is defined by $\\mathcal{O}(n)$ parameters; computation of the QR is more complicated, but it is possible to compute $A_{k+1}$ directly without computing $Q_k$.\n",
    "\n",
    "- This is called **implicit QR-step**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Theorem on implicit QR iteration\n",
    "\n",
    "All the implicit QR algorithms are based on the following **theorem**: \n",
    "\n",
    "Let \n",
    "\n",
    "$$Q^* A Q = H$$ \n",
    "\n",
    "be an irreducible upper Hessenberg matrix. Then, the first column of the matrix $Q$ defines all of its other columns. \n",
    "It can be found from the equation\n",
    "\n",
    "$$A Q = Q H. $$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Convergence of the QR-algorithm\n",
    "\n",
    "- The convergence of the QR-algorithm is a very **delicate issue** (see E. E. Tyrtyshnikov, \"Brief introduction to numerical analysis\" for details). \n",
    "\n",
    "\n",
    "**Summary.** If we have a decomposition of the form\n",
    "\n",
    "$$A = X \\Lambda X^{-1}, \\quad A = \\begin{bmatrix}A_{11} & A_{12} \\\\ A_{21} & A_{22}\\end{bmatrix}$$\n",
    "\n",
    "and \n",
    "\n",
    "$$\n",
    "\\Lambda = \\begin{bmatrix} \\Lambda_1 & 0 \\\\ \n",
    "0 & \\Lambda_2 \\end{bmatrix}, \\quad \\lambda(\\Lambda_1)=\\{\\lambda_1,\\dots,\\lambda_m\\}, \\ \\lambda(\\Lambda_2)=\\{\\lambda_{m+1},\\dots,\\lambda_r\\},\n",
    "$$\n",
    "\n",
    "and there is a **gap** between the eigenvalues of $\\Lambda_1$ and $\\Lambda_2$ ($|\\lambda_1|\\geq \\dots \\geq |\\lambda_m| > |\\lambda_{m+1}| \\geq\\dots \\geq |\\lambda_r| >0$), then the $A^{(k)}_{21}$ block of $A_k$\n",
    "in the QR-iteration goes to zero with \n",
    "\n",
    "$$\\Vert A^{(k)}_{21} \\Vert \\leq  C q^k, \\quad q = \\left| \\frac{\\lambda_{m+1}}{\\lambda_{m}}  \\right |,$$\n",
    "\n",
    "where $m$ is the size of $\\Lambda_1$.\n",
    "\n",
    "So we need to increase the gap!  It can be done by the **QR algorithm with shifts**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR-algorithm with shifts\n",
    "\n",
    "$$A_{k} - s_k I = Q_k R_k, \\quad A_{k+1} = R_k Q_k + s_k I$$\n",
    "\n",
    "\n",
    "The convergence rate for a shifted version is then\n",
    "\n",
    "$$\\left| \\frac{\\lambda_{m+1} - s_k}{\\lambda_{m} - s_k}  \\right |,$$\n",
    "\n",
    "where $\\lambda_m$ is the $m$-th largest eigenvalue of the matrix in modulus. \n",
    "- If the shift is close to the eigenvalue, then the convergence speed is better.\n",
    "- There are different strategies to choose shifts. \n",
    "- Introducing shifts is a general strategy to improve convergence of iterative methods of finding eigenvalues. \n",
    "- In next slides we will illustrate how to choose shift on a simpler algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Shifts and power method\n",
    "\n",
    "- Remember the power method for the computation of the eigenvalues.\n",
    "\n",
    "$$x_{k+1} := A x_k, \\quad x_{k+1} := \\frac{x_{k+1}}{\\Vert x_{k+1} \\Vert}.$$\n",
    "\n",
    "- It converges to the eigenvector corresponding to the largest eigenvalue in modulus. \n",
    "\n",
    "- The convergence can be very slow.\n",
    "\n",
    "- Let us try to use shifting strategy. If we shift the matrix as \n",
    "\n",
    "$$  A := A - \\lambda_k I$$\n",
    "\n",
    "and the corresponding eigenvalue becomes small (but we need large). \n",
    "- That is not what we wanted!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Inverse iteration and Rayleigh quotient iteration\n",
    "\n",
    "- To make a small eigenvalue large, we need to **invert the matrix**, and that gives us **inverse iteration**\n",
    "\n",
    "$$x_{k+1} = (A - \\lambda I)^{-1} x_k,$$\n",
    "\n",
    "where $\\lambda$ is the shift which is approximation to the eigenvalue that we want. \n",
    "\n",
    "- As it was for the power method, the convergence is linear.\n",
    "\n",
    "- To accelerate convergence one can use the **Rayleigh quotient iteration** (inverse iteration with adaptive shifts) which is given by the selection of the **adaptive shift**:\n",
    "\n",
    "$$x_{k+1} = (A - \\lambda_k I)^{-1} x_k,$$\n",
    "\n",
    "$$\\lambda_k = \\frac{(Ax_k, x_k)}{(x_k, x_k)}$$\n",
    "\n",
    "- In the symmetric case $A = A^*$ the convergence is **locally cubic** and **locally quadratic** otherwise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Singular values and eigenvalues (1)\n",
    "\n",
    "- Now let us talk about singular values and eigenvalues. \n",
    "\n",
    "- SVD \n",
    "\n",
    "$$A = U \\Sigma V^*$$\n",
    "\n",
    "exists for any matrix.\n",
    "\n",
    "- It can be also viewed as a reduction of a given matrix to the diagonal form by means of two-sided unitary transformations:\n",
    "\n",
    "$$\\Sigma = U^* A V.$$\n",
    "\n",
    "- By two-sided Householder transformation we can reduce any matrix to the **bidiagonal form** $B$ (how?)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Singular values and eigenvalues (2)\n",
    "\n",
    "- **Implicit QR algorithm** (with shifts) gives the way of computing the eigenvalues (and Schur form).\n",
    "\n",
    "- But we cannot apply QR algorithm directly to the bidiagonal matrix, as it is not diagonalizable in general case.\n",
    "\n",
    "- However, the problem of the computation of the SVD can be reduced to the **symmetric eigenvalue problem** in two ways:\n",
    "\n",
    "\n",
    "1. Work with the tridiagonal matrix \n",
    "\n",
    "$$T = B^* B$$ \n",
    "\n",
    "2. Work with the extended matrix \n",
    "\n",
    "$$T = \\begin{bmatrix} 0 & B \\\\ B^* & 0 \\end{bmatrix}$$\n",
    "\n",
    "- The case 1 is OK if you **do not form T directly**!\n",
    "\n",
    "- Thus, the problem of computing singular values can be reduced to the problem of the computation of the eigenvalues of symmetric tridiagonal matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Algorithms for the SEV  (symmetric eigenvalue problem)\n",
    "\n",
    "Done:\n",
    "- QR algorithm: the \"gold standard\" of the eigenvalue computations\n",
    "- RQI-iteration: Rayleigh quotient iteration is implicitly performed at each step of the QR algorithm\n",
    "\n",
    "Next slides:\n",
    "- Divide-and-conquer\n",
    "- Bisection method\n",
    "- Jacobi method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Divide-and-conquer\n",
    "\n",
    "- Suppose we have a tridiagonal matrix, and we split it into two blocks:\n",
    "\n",
    "$$T = \\begin{bmatrix} T'_1 & B \\\\ B^{\\top} & T'_2 \\end{bmatrix}$$\n",
    "\n",
    "- We can write the matrix $T$ as\n",
    "\n",
    "$$T = \\begin{bmatrix} T_1 & 0 \\\\ 0 & T_2 \\end{bmatrix} + b_m v v^*$$\n",
    "\n",
    "where $vv^*$ is rank 1 matrix, $v = (0,\\dots,0,1,1,0,\\dots,0)^T$.\n",
    "\n",
    "- Suppose we have decomposed $T_1$ and $T_2$ already:\n",
    "\n",
    "$$T_1 = Q_1 \\Lambda_1 Q^*_1, \\quad T_2 = Q_2 \\Lambda_2 Q^*_2$$\n",
    "\n",
    "- Then (check),\n",
    "\n",
    "$$\\begin{bmatrix} Q^*_1 & 0 \\\\ 0 & Q^*_2 \\end{bmatrix} T\\begin{bmatrix} Q_1 & 0 \\\\ 0 & Q_2 \\end{bmatrix} = D + \\rho u u^{*}, \\quad D = \\begin{bmatrix} \\Lambda_1 & 0 \\\\ 0 & \\Lambda_2\\end{bmatrix}$$\n",
    "\n",
    "- I.e. we have reduced the problem to the problem of the computation of the eigenvalues of <font color='red'> diagonal plus low-rank matrix </font>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Diagonal-plus-low-rank matrix\n",
    "\n",
    "It is tricky to compute the eigenvalues of the matrix\n",
    "\n",
    "$$D + \\rho u u^* $$\n",
    "\n",
    "The characteristic polynomial has the form\n",
    "\n",
    "$$\\det(D + \\rho uu^* - \\lambda I) = \\det(D - \\lambda I)\\det(I + \\rho (D - \\lambda I)^{-1} uu^*) = 0.$$\n",
    "\n",
    "Then (prove!!) \n",
    "\n",
    "$$\\det(I + \\rho (D - \\lambda I)^{-1} uu^*) = 1 + \\rho \\sum_{i=1}^n \\frac{|u_i|^2}{d_i - \\lambda} = 0$$\n",
    "\n",
    "Hint: find $\\det(I + w u^*)$ using the fact that $\\text{det}(C) = \\prod_{i=1}^n\\lambda_i(C)$ and $\\text{trace}(C) = \\sum_{i=1}^n \\lambda_i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Characteristic equation\n",
    "\n",
    "$$1 + \\rho \\sum_{i=1}^n \\frac{|u_i|^2}{d_i - \\lambda} = 0$$\n",
    "\n",
    "How to find the roots?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt4AAAJMCAYAAADNOqEXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3hU1cI18HVCEhKK6YVI7733IoKgIMgVUURRRPQCKoq+XhtSREGuDbChEUUuHa5KMwS8dJUuCSVMQkmAFEhvkzKZzOzvDxI+VHqyz07OWb/nmSdtMlmcCZM1e/bZWxNCgIiIiIiI5HJRHYCIiIiIyAxYvImIiIiIdMDiTURERESkAxZvIiIiIiIdsHgTEREREemAxZuIiIiISAcs3kREREREOqgQxVvTtIaaps3TNM2iaZpV07TskvcXaZrWV3U+IiIiIqKyclUdQNO0cQC+AOBZ8qk8AG4AmpdcnAB2qUlHRERERFQ+lI54a5o2CsC3uFS6vwDQSAhRQwhRDUAwgCcB7FEYkYiIiIioXGiqtozXNC0QQDQAHwBThBBzlAQhIiIiItKByhHv53CpdMcA+EBhDiIiIiIi6VQW79Elb5cIIZwKcxARERERSaekeGua5gegScmHv2ma1l/TtC2apmVqmpavadoJTdP+rWmav4p8RERERETlTdWId5Mr3r8XwNaSt1VKPtcCwBsAIjVNa6FzNiIiIiKicqfk5EpN0wYBCC/5UACIAvCMEOKApmkuAO4DsBhAIAALgLZCiOKr3M54AOMBwNPTs1OdOnV0SP9nTqcTLi4VYjl0kqwy39dpaWnIyMhA06ZNL38uOTkZVqsVjRo1Upis4rqV+zs3NxcXLlxAvXr1ULVqVQBATk4OLl68iAYNGsDNzU1mVMMrKChAfHw8ateujWrVqpX7sa3M/7fLQ3FxMWJjYxEUFAQvLy/Y7XbExcUhODgYd9xxh+p45crs97WZqLqvT548mSaECLjmFYQQul8A3I9LhVsAsANofJXrDLniOo/c6DY7deokVNixY4eSn0v6q8z39auvviqqVav2p8+NHz9eBAcHK0pU8d3K/b1ixQoBQFgslsufW7p0qQAgTp06JSGduWzfvl0AuHyfLFmyRAAQp0+fLpfbr8z/t8tDbGysACAWL14shBDi5MmTAoBYtmyZ4mTlz+z3tZmouq8BHBLX6auqnvZZr3g/TAhx+q9XEEKEAThZ8uEAXVIRGVRRURHc3d1VxzAsu90OABzZloTHV67S4+vqqnxPPSLDU1W8k654P+Y61yv9mv5zSIgMhMVbrtLiwmMsB4u3XMXFl2Zy8vgSyaeqeMcCKCh5/2YmmavZ5YfIIFi85SoqKgLA4iILi7dcHPEm0o+S4i0urdu9s+TD5te5arOSt+ekBiIyOBZvuVgM5eLxlYsj3kT6UXlq79KSt0M0TWv81y9qmjYEQOkSDJt0S0VkQCzecrEYysXjK1dp8eaIN5F8Kov3agB/AHAFsFbTtC4AoGmaS8lyg9+VXO8AgDA1EYmMgcVbLs7xlosjsnLxiQ2RfpQ9vRVCODVNexDALgCtARzQNC0XlzbRqVZytRgAD5csz0JEt4nFWy7O8ZaLc5Dl4og3kX6UriIvhEgA0A7ATADHcal0CwARAN4G0FkIEa8uIZEx2O12Fm+J7HY7NE1DlSpVbnxlumUckZWLx5dIP8qf3gohrADeKbkQkQQc8ZbLbreztEjEYigXR7yJ9MN9U4lMgMVbLr6iIBeLt1ycykOkHxZvIhNg8ZarqKiIpVAiFm+5ePIqkX5YvIlMgMVbLk41kYvFWy5ONSHSD4s3kQlwRFYuFm+5OBVCLj6xIdIPizeRCdhsNlStWlV1DMPiHG+5iouLUaVKFWiapjqKIXHEm0g/LN5EJlBYWAgPDw/VMQyLryjIxVcU5OKIN5F+WLyJTIDFWy4WQ7l4fOXiiDeRfli8iUzAZrOxeEvEYigXj69cXNWESD8s3kQGJ4TgiLdkXDVGLhZvuXjyKpF+WLyJDK6oqAgAWLwlYjGUi8dXLk41IdIPizeRwRUWFgJg8ZaJxVAuHl+5eHIlkX5YvIkMjsVbPhZDuYqLizkaKxFHvIn0w+JNZHClxZvreMvD5QTl4hMbuex2O1xcXODiwkpAJBv/lxEZnM1mA8ARb5mKior4xEYiFm+5+IoCkX5YvIkMjlNN5OOqMXKxeMtVXFzM40ukExZvIoNj8ZbPZrNxxFsiFm+57HY7R7yJdMLiTWRwLN7ycYMiuVi85eKIN5F+WLyJDI7FW77CwkKOeEvEDYrk4og3kX5YvIkMjsVbPk41kYvHVy6eXEmkHxZvIoPjqiZyCSE41UQyFm+5OJWHSD8s3kQGx3W85SoqKgLA4ysTi7dcHPEm0g+LN5HBcaqJXKWvKLAYysPiLRdPriTSD4s3kcGxeMvF4ysfNyiSiydXEumHxZvI4FgM5eKIt3w2m42rmkjE4k2kHxZvIoMrKCgAwOItC4u3fJxqIhdfUSDSD4s3kcEVFBTA3d0dVapUUR3FkPiKgnws3nKxeBPph8WbyODy8/NRrVo11TEMiyPechUXF8PpdPL4SsSpPET6YfEmMjgWb7lYvOXi8ZWPO4MS6YfFm8jg8vPz4enpqTqGYXGqiVxcJ10+TuUh0g+LN5HBFRQUcMRbIo7IysXjKx9HvIn0w+JNZHDXm2oihNA5jfFwxFsuFm/5rnVyJR8fiMofizeRwV1rqommaQrSGA+LoVylx5cjsvL89eRKPjYQycPiTWRwnGoiF4u3XDy+8nGqCZF+WLyJDI6rmsjFqSZysXjLx5MrifTD4k1kcFzVRK7S4s3iIgdXNZGPI95E+mHxJjI4jnjLVVBQAAB8ciMJR7zlcjqdsNvtPL5EOmHxJjI4zvGWq7R4c6qJHCzectntdgA8eZVILyzeRAYmhOCIt2QFBQXw8PCAiwsfTmVg8ZaLq8YQ6Yt/KYgMzG63w+FwcBqERAUFBTy+ErF4y8U59ET6YvEmMrD8/HwA4Ii3RCzecnFEVq7S4s3jS6QPFm8iAyst3iyG8rB4y1W6agyPsRx8YkOkLxZvIgOzWq0AgJo1aypOYlxcrlEurpMuF6eaEOmLxZvIwPLy8gAANWrUUJzEuDjiLReXa5SLI95E+mLxJjKw0hHv6tWrK05iXFyuUS4u1ygXR7yJ9MXiTWRgpcWbI97ycMRbrsLCQri7u3O5Rkk44k2kLz6SERkYi7d8LN5y8fjKxVVNiPTF4k1kYJzjLR+LoVyFhYWcZiIRp5oQ6YvFm8jAOOItH4u3XDy+cnGqCZG+WLyJDIwnV8rH5QTlKigo4Ii3RKXFm8eYSB8s3kQGZrVa4eLiwj+qEnFEVq7CwkIeX4m4TjqRvli8iQwsLy8PNWrUgKZpqqMYktPphM1m43KCEnHEWy4WbyJ9sXgTGZjVauX8bom4nbl8fEVBLq6TTqQvFm8iA7NarZzfLVHpqjE8xvJwVRO5OOJNpC8WbyID44i3XCze8nHEW67S4s3lBIn0weJNZGC5ubm44447VMcwLBZv+XhypVyFhYVwdXWFq6ur6ihEpsDiTWRgOTk5LN4SsXjLx5Mr5eITGyJ9sXgTGRiLt1ws3vJxqolcnENPpC8WbyIDY/GWi8VbPhZDuXh8ifTF4k1kYCzecnFnULmcTienQkjG4k2kLxZvIoOy2Wyw2Wws3hJxxFuu0jWmeXzl4Rx6In2xeBMZVG5uLgCweEvE4i0Xj698HPEm0heLN5FB5eTkAGDxlonFUK78/HwAQLVq1RQnMS4WbyJ9sXgTGRSLt3x5eXnQNI3FRRI+sZGPxZtIXyzeRAbF4i1fXl4eqlevDk3TVEcxJI54y8fiTaQvFm8ig8rOzgYA1KxZU3ES4yot3iQHR7zl46oxRPpi8SYyqKysLACAj4+P4iTGdTPFWwihUxrjudGIN49t2XHEm0hfLN5EBpWRkQEA8PX1VZzEuHJzc6/5igKnn5TdtUa8eWzLD4s3kb5YvIkMKjMzEwDg5eWlOIlxXa94U9mVFm/O8ZaH63gT6YvFm8igMjIy4O3tjSpVqqiOYli5ubk8eVWi0qkmnOMtT35+Pp/YEOmIxZvIoDIyMji/W7KcnByOeEvEEW+57HY77HY7jy+RjipM8dY0rYamafGapomSy1jVmYgqs8zMTM7vloxTTeTicoJyFRQUAODxJdJThSneAGYBqK06BJFRZGRk3LB4c1WIsmHxlisvLw9Vq1bldClJuGoMkf4qRPHWNK0jgEkA9qvOQmQUN5pqwpUhysbhcCAvL49zvCXKz8/n/G6JrjXizccGInmUF29N01wAhJZ8+JzKLERGwqkmclmtVgDcoEgmq9XK4i0Rp/IQ6U958QbwIoDOAL4SQkSoDkNkBEKIm5pqQrcvNzcXAIu3TFarlcdXIhZvIv0pLd6apt0J4D0AyQCmqsxCZCRWqxUOh4OrmkjE4i0f59DLxeJNpD9XxT//cwA1ATwvhMhWnIXIMLhrpXw5OTkAwDneElmtVtSoUUN1DMNi8SbSn7IRb03THgAwHMBOIcQyVTmIjKi0eHPEWx6OeMvHEW+5WLyJ9KdkxFvTtOoAvgBgB/BCGW5nPIDxABAUFISdO3eWS75bYbValfxc0l9luq8PHz4MADh37tw1MyclJcFut1eaf5PebnR///777wCAU6dOwel0/u3rJ06cAADs378fiYmJUjIaXWpqKoKDg/92P1gsFgCXjm1CQkKZf05l+r9dnv744w8AwNGjR5Genn7586W/rxaLxXDHxaz3tRlV1Pta1VSTdwHUBfChEOLE7d6IEOIbAN8AQOfOncXdd99dPuluwc6dO6Hi55L+KtN9nZaWBgC455570KZNm6teZ/Xq1XBzc6s0/ya93ej+PnXqFADg3nvvRZ06df729dLy0q1bNzRp0kRKRqMrLi5G48aN/3Y/lJbtbt26oXHjxmX+OZXp/3Z5iomJAQD0798ftWrVuvz5M2fOAABatGhhuONi1vvajCrqfa37VBNN09oDmAwgHpcKOBGVM041kS8zMxMAj7FMnGoiF6eaEOlPxYj3pwCqAHgbgKZp2rXOnKla8jWnECJft3REBlBaCnlypTyZmZmoUqUK15mWxG63w2az8eRKiVi8ifSn4uTKeiVvlwDIvcql1NclH9/2VBQis8rIyIC7uzs8PT1VRzGsrKws+Pj4cJc/SbhBkXz5+flwdXWFm5ub6ihEplERNtAhonKWkpKCwMBAlkKJMjMzOc1EIhZv+fLy8jjaTaQz3Yu3EKK+EEK71uWKqz5d8rn6emckquySk5MRFBSkOoahsXjLVbpcI6eayMOdQYn0xxFvIgO6ePEii7dkLN5ycZ10+bhBEZH+WLyJDCg5ORnBwcGqYxhaVlYWvL29VccwrOzsS5sZ8xjLw+JNpD8WbyKDcTqdSElJ4Yi3ZBzxlisrKwsAi7dMLN5E+lO1gc41/WWeNxHdooyMDBQXF7N4SySEYPGWjMVbPqvVylfGiHTGEW8ig0lOTgYAFm+JrFYrHA4Hi7dELN7yccSbSH8s3kQGU1q8OZIlD0uhfFlZWXBzc+Na9BKxeBPpj8WbyGA44i0ft4uXr/TkVa5FLw+LN5H+WLyJDObixYsAWLxlYvGWLzMzk68oSCSEYPEmUoDFm8hgkpOT4ebmxlIoEYu3fFyuUS6bzQaHw8HiTaQzFm8igyndtZIv0ctTOsebxVseFm+5rFYrAO4MSqQ3Fm8ig+F28fKVjnizGMrD4i1X6c6gLN5E+mLxJjKYCxcucEUTyTIyMuDi4gIvLy/VUQyLxVsuFm8iNVi8iQzm/PnzqFu3ruoYhpaamgp/f3+4uPAhVBYWb7mys7MBgE8eiXTGvxpEBpKXl4eMjAzUqVNHdRRDS0lJQWBgoOoYhlVYWIjCwkIWb4m4Fj2RGizeRAYSHx8PABzxlozFW67S0ViWQnk44k2kBos3kYGcP38eAIu3bMnJySzeEnE0Vj4+uSFSg8WbyEBKR7w51UQujnjLxeItX+kx5og3kb5YvIkM5Pz589A0DXfeeafqKIZVWFiInJwcFm+JWLzly87OhoeHB6pWrao6CpGpsHgTGUh8fDxq1aoFNzc31VEMKzU1FQBYvCVi8ZYvKyuLo91ECrB4ExkIlxKULyUlBQCLt0ws3vJlZ2ezeBMpwOJNZCAs3vKxeMuXkZEBAPDx8VGcxLiys7P5xIZIARZvIoMQQiA+Pp4nVkrG4i1famoqqlevDk9PT9VRDItTTYjUYPEmMoiLFy+isLAQ9evXVx3F0Fi85eOqMfJxxJtIDRZvIoM4deoUAKBJkyaKkxhbcnIyPDw8UKNGDdVRDCs1NRUBAQGqYxgaR7yJ1GDxJjIIFm99lI7GapqmOophsXjLxxFvIjVYvIkM4vTp03Bzc+PJlZJx10r5UlNTeYwlKioqQkFBAUe8iRRg8SYyiFOnTqFBgwZwdXVVHcXQkpKSuEGRREIIpKSkcMRbIm4XT6QOizeRQZw6deqWp5kIISSlMa6EhATUrl1bdQzDys3NRVFREYu3RKXF+0Yj3nx8ICp/LN5EBiCEwOnTp2+peHOO8q2zWq3Iyspi8ZaIO4PKV7pB0bWKNx8biORh8SYygAsXLiA/Px+NGzdWHcXQEhISAIBrpUtUulwjR7zl4VQTInVYvIkM4OTJkwC4oolspcWbI97ylI54s3jLc6MRbyKSh8WbyACioqIAAK1atVKcxNhYvOXjVBP5bnaONxGVPxZvIgM4fvw4vL29ERISojqKoZUWb65qIg+nmsiXlpYGAPD391echMh8WLyJDODYsWNo06YNT4qSLCEhAf7+/vDw8FAdxbBSU1NRvXp1eHp6qo5iWGlpafDw8ED16tVVRyEyHRZvokpOCIHjx4+jdevWqqMYXkJCAk+slIyb58iXmpoKf39/PlEnUoDFm6iSS0hIQHZ2Ntq0aaM6iuFxDW/5Lly4gKCgINUxDC01NZVTeYgUYfEmquSOHz8OABzx1kF8fDznd0uWmJjIJzeSpaWlcX43kSIs3kSV3NGjRwGweMuWk5ODjIwM1K9fX3UUQ0tMTOSTG8k44k2kDos3USV36NAhNGjQAD4+PqqjGNqZM2cAAI0aNVKcxLhycnKQm5vL4i0ZizeROizeRJXcwYMH0aVLF9UxDC82NhYAi7dMiYmJALhco0w2mw25ubmcakKkCIs3USWWmpqKc+fOsXjrgCPe8pUWb87xlqd0DW+OeBOpweJNVIkdPHgQAFi8dXDmzBn4+/vjjjvuUB3FsDjiLV/pzqAc8SZSg8WbqBI7dOgQNE1Dx44dVUcxvDNnznC0W7LSnUG5A6s8HPEmUovFm6gSO3DgAFq0aIGaNWuqjmJ4LN7yJSYmwtfXl7tWSlQ64s3iTaQGizdRJeVwOPDbb7+hd+/eqqMYXlFREc6fP4+GDRuqjmJoXMNbPk41IVKLxZuokjp27Biys7Nx1113qY5ieOfOnYPT6eSIt2QJCQmc3y1ZWloaNE2Dr6+v6ihEpsTiTVRJ7d69GwBYvHUQExMDAGjatKniJMbGzXPkS01NhZ+fH6pUqaI6CpEpsXgTVVK7d+9G/fr1UadOHdVRDO/EiRMAgJYtWypOYlz5+flITk5GgwYNVEcxtIsXLyIwMFB1DCLTYvEmqoScTid2796NPn36qI5iClFRUQgJCYG3t7fqKIYVFxcHAJxHLxlfVSBSi8WbqBKKjIxEamoqBg4cqDqKKZw4cYKj3ZKV7gzK4i1XUlISizeRQizeRJXQ5s2bAQD33nuv4iTG53Q6YbFYWLwlY/GWz+Fw4OLFiyzeRAqxeBNVQps3b0bHjh0RFBSkOorhxcfHIy8vD61atVIdxdBiY2NRs2ZN+Pn5qY5iWCkpKXA4HNygiEghFm+iSiY7Oxt79uzBoEGDVEcxBZ5YqY/Y2Fg0aNAAmqapjmJYiYmJAMARbyKFWLyJKplNmzbB4XBgyJAhqqOYwvHjxwGweMsWGxvLaSaSlRZvjngTqcPiTVTJ/PjjjwgODkb37t1VRzGFw4cPo169etxwRCIhBOLi4li8JUtKSgLAEW8ilVi8iSqR/Px8hIeHY/jw4XBx4X9fPRw+fBgdO3ZUHcPQLl68iIKCAhZvyRITE+Hi4sJzQ4gU4l9uokpk8+bNyM/Px4gRI1RHMYWcnBycPHmyTMVbCFGOiYypdGfQJk2a3NL38djemsTERAQHB3PXSiKFWLyJKpGlS5ciKCiI28Tr5MiRIwBwW8WbJwnePIvFAgBo0aLFTV2fx/b2cA1vIvVYvIkqiZSUFPz888944okn4ObmpjqOKRw+fBjA7RVvunkWiwU1atRA7dq1VUcxtMTERJ5YSaQYizdRJbFixQoUFxfj6aefVh3FNP744w8EBwcjODhYdRRDs1gsaN68OUeyJRJC4Ny5c6hXr57qKESmxuJNVAkIIfD999+jS5cu3MhFR3v37uXqMTqwWCw3Pc2Ebk96ejqsVisaNGigOgqRqbF4E1UChw8fxtGjRzF27FjVUUwjMzMTp0+fRq9evVRHMbTc3FwkJiayeEsWGxsLACzeRIqxeBNVAvPnz0eNGjXw+OOPq45iGlFRUQCAnj17Kk5ibNHR0QBu/sRKuj1xcXEAWLyJVGPxJqrgEhISsGrVKjz77LPw9vZWHcc0jh8/Dnd3d55YKdmJEycAsHjLxuJNVDGweBNVcF988QWcTideeukl1VFM5fjx4+jUqRM8PDxURzG0yMhIVKtWDY0bN1YdxdDi4uLg7++PmjVrqo5CZGos3kQVWE5ODkJDQzFixAiOVOkoLy8PMTEx6N27t+oohhcREYG2bdtyUxfJ4uLi+BhCVAGweBNVYPPmzUNWVhbeeOMN1VFM5bfffkNxcTEGDBigOoqhCSEQGRmJ9u3bq45ieCzeRBUDizdRBZWWloZPPvkEI0aMQKdOnVTHMZWtW7fCzc2NI96SxcXFITs7Gx06dFAdxdAcDgfOnTuHhg0bqo5CZHos3kQV1L///W/k5eXh3XfflfYzhBDSbrsy27ZtG1q2bIlq1aqpjmJokZGRAMDiLVlCQgLsdvstF28+PhCVPxZvogrozJkz+OKLL/Dkk0+iZcuWUn4Gdwm8urS0NERGRnI1Ex1ERESgSpUqaN26teoohmaxWAAAzZs3v6nr87GBSB4Wb6IKRgiBSZMmwd3dHe+//77qOKYTHh4OIQS6du2qOorhHThwAK1atYKnp6fqKIZ2q8WbiORh8SaqYNauXYvNmzdj5syZCAkJUR3HdH7++WcEBwejadOmqqMYmtPpxP79+9GjRw/VUQwvOjoafn5+CAgIUB2FyPRYvIkqkIyMDLz44oto06YNXnzxRdVxTMdut2Pz5s0YMmQIXFz48CiTxWJBdnY2i7cOLBYLR7uJKgj+ZSGqQF544QWkpKTg+++/h6urq+o4pvPrr78iJycHQ4cOVR3F8Pbu3QsALN46iI6O5s6gRBWEsuKtaVpdTdNe1jRto6Zp5zVNs2malqtp2hFN0/6taVotVdmIVFi1ahVWrVqFGTNmcPlARdasWYNq1aph4MCBqqMY3t69e+Hn54cmTZqojmJo6enpSE1N5Yg3UQWhZEhN07Q6AM4CuPLU6RwA1QG0LbmM1zRthBBih/4JifQVHR2Nf/7zn+jRowfefPNN1XFMyW6344cffsCwYcNQvXp11XEMb8+ePejevTtX0JAsOjoaAE+sJKooVI14l+4NHAbgEQC+QggvANUA3A8gDoAPgHWapgWriUikj9zcXDz00EPw9PTEmjVrOMVEkW3btiE9PR2PPfaY6iiGd+HCBURHR6Nv376qoxjesWPHAEDasqREdGtU/YXPBNBBCHHkyk8KIYoAhGuadj+ACAB3AJgAYKb+EYnkKy4uxmOPPYaYmBhs3boVtWvXVh3JtFauXAkvLy/cd999qqMY3o4dl17I7N+/v+IkxhcZGQkvLy/Ur19fdRQigqIRbyFE9l9L91++Hg1gX8mHnOxKhiSEwIQJExAWFoYFCxagX79+qiOZVmFhIdauXYuHHnoIVatWVR3H8LZv3w5vb2+0b99edRTDi4yMRPv27Tmlh6iCqMirmqSXvK1y3WsRVVLTp0/HokWLMH36dEyYMEF1HFNbu3YtcnNz8fjjj6uOYgo7duxA3759UaUKH95lcjgcOHr0KDp06KA6ChGVqJDFW9M0VwC9Sj48rjILUXkTQmDGjBmYNWsWnnnmGbzzzjuqI5leaGgoGjZsyKkPOoiNjUVsbCyPtQ5OnjyJgoICvrJAVIFUyOIN4AUAwQCcAJYozkJUboQQePPNN/Huu+9i3LhxCA0N5UvAilksFuzatQvjx4/npjk6CAsLAwDcf//9ipMYX2RkJACweBNVIBVu+QRN09oCeL/kwy+EEFHXue54AOMBICgoCDt37pQf8C+sVquSn0v6K+t9XVxcjPnz5yMsLAzDhg3D6NGj8euvv5ZfwFuUmJgIu91u+t/fL7/8Eq6urmjatOmfjkVZ7+8TJ04AAPbv34+kpKQypjSOpUuXonbt2khISEBCQsJt3caVxzYxMbHMmYz6OL5+/Xq4uroiJSXllv59Fy5cAHDpSanRjotR72v6uwp7XwshKswFQC0AsQAEgEMAPG72ezt16iRU2LFjh5KfS/ory32dmZkp7rnnHgFATJkyRTidzvILdpuef/554e/vrzqGUnl5ecLHx0c8+uijf/taWf9vL1++XAAQMTExZbodI7FaraJq1ari5ZdfLtPtLFu2TAAQJ0+eLJdcRn0c79Onj+jatestf19sbKwAIBYvXiwhlVpGva/p71Td1wAOiev01Qrzuqqmab4AfgHQAMApAEOEEIVqUxGV3ZEjR9CtWzfs3r0bixcvxuzZszm9pIJYtGgRMjMzMWnSJNVRTGHr1q2w2WwYMmSI6iiGZ7fbcejQIfTo0UN1FCK6QoUo3pqmeQHYAqA1gPMABgghktWmIiobIQRCQ0PRrVs35ObmYuvWrXjqqadUx6ISxcXF+OSTT9CzZ0/07t1bdRxTWLNmDfz8/Lhxjg6OHDmCgsmFQQwAACAASURBVIICFm+iCkZ58dY0rTqATQA6A7iIS6X7vNpURGWTmpqKRx99FBMnTkTfvn0RGRmJu+66S3UsusJ///tfnD17Fm+88YbqKKZQUFCADRs24KGHHoKbm5vqOIa3d+9eAGDxJqpglBZvTdM8AWwE0BOX1u0eIIQ4pTITUVkIIbBixQq0aNEC69atw5w5cxAeHo7AwEDV0egKTqcTH3zwAVq0aIGhQ4eqjmMK4eHhsFqtGDlypOooprBv3z6EhISgTp06qqMQ0RWUrWqiaZo7gJ8A9AOQBeBecZ0VTIgqupMnT+KVV17Bpk2b0K1bN3z33Xdo1aqV6lh0FT/88AOOHDmCpUuXcglBnaxevRoBAQG4++67VUcxPCEEfv31V/Ts2ZPnkxBVMEr+4miaVgXACgCDAOQCGCyEOKwiC1FZZWRk4OWXX0arVq2we/duzJ07F7///jtLdwVVXFyMadOmoXXr1njsscdUxzGFvLw8/PzzzxgxYgRcXSvcKraGc+bMGcTHx3OTIqIKSNUjYC8AI0redwOw7jrPyuOFEF10SUV0C/Ly8vD1119j9uzZyM7OxjPPPIP33nsPQUFBqqPRdSxZsgQnT57EunXruGW5Tv773/8iPz+fT3R0smPHDgBAv379FCchor9SVbyvHGn3KLlcC5cUpAolNzcXCxYswCeffILU1FQMHDgQH3/8Mdq2bas6Gt1Abm4upk2bhq5du2LYsGGq45jGwoUL0axZM/Tp00d1FFPYvn07atWqhWbNmqmOQkR/oaR4CyF2AuDEM6pUUlNTMXXqVHz11VfIyMjAoEGDMG3aNPTs2VN1NLpJs2bNQlJSEn788UfOfdVJVFQU9uzZg48//pjHXAdCCOzYsQMDBgzg8SaqgDjZjug6hBD4/fff8fnnn+OHH36AEALDhg3DlClT0LVrV9Xx6BbExMRg3rx5ePrpp9G9e3fVcUxj4cKFcHNzw5gxY1RHMYXIyEgkJydjwIABqqMQ0VWweBNdxcWLF7F8+XL85z//wbFjx+Dt7Y1HHnkEc+bMQYMGDVTHo1skhMBLL72EatWqYc6cOarjmEZ+fj6WLl2K4cOHIyAgQHUcU9i0aRMAYPDgwYqTENHVsHgTlcjLy0NYWBiWLFmCzZs3w+FwoGvXrggNDcXo0aNx8OBBlu5KavXq1fjll18wf/58nvyqo8WLFyMjIwOTJk1SHcU0wsLC0KVLF/6eE1VQLN5kajk5Ofj555/xww8/YPPmzSgoKEBISAhee+01jBkzBi1atFAdkcooOTkZkyZNQteuXfHCCy+ojmMaDocDc+fORbdu3dC7d2/VcUwhLS0N+/btw/Tp01VHIaJrYPEm0zlz5gzCw8MRHh6OrVu3oqioCLVq1cK4ceMwYsQI3HXXXVxmziCEEHjuuedgtVqxePFiriGto/Xr1+PMmTP44IMPeJKfTtavXw8hBB544AHVUYjoGvhXiAwvLy8Pu3fvxubNmxEeHo5Tp04BABo2bIgXXngBI0aMQI8ePbiDoQGtWrUKa9euvbw9POlDCIEPP/wQDRs2xIMPPqg6jmmsXr0aDRs2RMeOHVVHIaJrYPEmw8nOzsbvv/+O3bt3Y9euXTh06BCKi4vh4eGBfv364cUXX8SgQYPQpEkT1VFJojNnzmDixIno0aMHXn31VdVxTCU8PBz79+/HN998w1ePdJKamort27fjjTfe4CsMRBUYizdVakIInD17FgcOHMD+/fuxa9cuREZGwul0ws3NDV26dMFrr72Gu+++G3369IGnp6fqyKQDm82GRx99FC4uLlixYgXLn46EEJg2bRoaNmyIsWPHqo5jGj/++CMcDgceffRR1VGI6DpYvKlSSUlJwcGDB3Hw4EEcOHAABw8eRFpaGgDAw8MD3bt3x9SpU9G3b190794d1apVU5yYVHj99dfxxx9/YO3atahfv77qOKaybt06HD58GIsXL4abm5vqOKaxevVqNG/eHG3atFEdhYiug8WbKiS73Y6TJ0/i6NGjf7okJCQAAFxcXNCyZUsMGzYMXbp0QdeuXdG6dWu4u7srTk6qrV27Fp999hkmT57M+cU6s9vteOutt9CsWTOMHj1adRzTSEhIwK5duzB9+nROMyGq4Fi8SSm73Y64uDicPHkSFosFx44dw9GjR2GxWFBUVAQAcHNzQ4sWLdC3b1+0b98eXbt2RceOHVGjRg3F6amiiYqKwlNPPYXOnTvjww8/VB3HdL744gvExMRg48aNXEFGR9999x2EEHjqqadURyGiG+AjI0knhEBqaipiYmJw8uRJxMTEXL6cOXMGxcXFl6975513ok2bNrjvvvvQtm1btG3bFs2aNeNINt1Qeno6hg0bhurVq2Pt2rX8ndFZSkoKZs6ciUGDBmHIkCGq45hGcXExvv32W9x7773c4IuoEmDxpnKRn5+Ps2fP4uzZs4iLi/vbJSsr6/J13d3d0aRJE7Rq1QoPPfQQmjZtimbNmqFZs2bw9fVV+K+gysput+Phhx9GYmIidu7cidq1a6uOZDrTpk1DXl4e5s2bx+kOOtq8eTMSEhLw6aefqo5CRDeBxZtuSAiBtLQ0JCQk/OlyZbFOTk7+0/d4eHigfv36aNCgAXr06HG5XDdt2hT16tXjKhNUrl566SXs3LkTS5YsQffu3VXHMZ0DBw5g4cKFePnll9G8eXPVcUwlNDQUwcHB3DSHqJJg8Ta5oqIiJCcn48KFC5cLdWJi4p8KdmJiImw225++r0qVKqhTpw4aNGiAIUOGoEGDBn+6BAcHc9SLdPHpp5/i66+/xuuvv44nn3xSdRzTsdlsGDduHO68807MmDFDdRxTOXXqFMLCwjBlyhSuIENUSbB4G5DdbkdqaiouXryI5OTk677NzMz82/e7u7ujdu3aqF27Nrp163b5/SsvQUFBHLUm5VavXo1XXnkFDz74IN5//33VcUxp9uzZiIqKQlhYGLy8vFTHMZW5c+fC3d0dkyZNUh2FiG4Si3cFJ4RAdnY20tPTkZaWdvlyrY9TUlIur2v9VzVr1kRQUBCCg4PRsmVL9O/f//LHwcHBl0u1v78/R6tNQgihOsJt2759O8aMGYPevXtzkxxFjhw5gjlz5uCJJ57A/fffrzqOqSQnJ+P777/HU089heDgYCk/ozI/PhBVVCzeOrHb7cjKykJmZuZV32ZlZSEjI+Nygb7y7ZWrflzJ1dUVfn5+8Pf3h7+/P5o3b44+ffpcLtKlpTooKAhBQUGoXr26zv9qqsgq85OryMhIPPjgg2jSpAnWr1/PHUkVKCgowJNPPglfX1/Mnz9fdRzT+eKLL1BUVIRXX3213G+7Mj82EFV0LN5lEBMTg/DwcERERFy3VGdmZiI/P/+6t+Xu7g4fHx/4+/vDz88PzZo1Q69evf5UrEu/Vvr+HXfcwQdIMp24uDgMHjwY3t7e2Lx5M3x8fFRHMqV//etfOHbsGDZt2gQ/Pz/VcUwlPT0dn332GYYPH46mTZuqjkNEt4DFuwy2bdt2eZMOTdPg5eUFb29v+Pj4wNvbG82aNYO3t/efPvfXt6Xve3h4sEQT3UBqairuu+8+2Gw2bNu2jcsGKrJu3TosWLAA//d//4fBgwerjmM6H374IXJzczFz5kzVUYjoFrF4l8Fjjz0GX19fDBo0CHfccQdcXFxURyIyLKvViiFDhiA+Ph7btm1Dy5YtVUcypfj4eIwbNw6dOnXCnDlzVMcxnaSkJHz++ecYPXo0WrdurToOEd0iFu8y8PHxQXBwMLy9vVVHITI0u92OkSNH4o8//sDatWvRs2dP1ZFMqaioCKNGjYLdbsfKlSu5O6gC7777Lux2O9555x3VUYjoNrB4E1GF5nQ68cwzzyA8PBwLFy7EsGHDVEcyrZdeegl79uzB6tWr0aRJE9VxTOfw4cP45ptv8OKLL6JRo0aq4xDRbeDcCCKq0N58800sXboU7733Hp599lnVcUwrNDQUoaGhePPNNzFy5EjVcUzH6XRi0qRJCAgI4NxuokqMI95EVGHNnTsXH330EV544QW8/fbbquOY1u+//44XX3wRgwYNwqxZs1THMaWlS5di7969+P777zm9kagS44g3EVVIy5cvx6uvvoqHH34Yn376KVf9USQhIQEjRoxAvXr1uFGRIpmZmXj99dfRvXt3jBkzRnUcIioDjngTUYXzyy+/YOzYsbj77ruxdOlSlj1FcnJyMHToUOTl5WHbtm1cM12RyZMnIz09HZs3b+bqWUSVHIs3EVUoBw8exEMPPYRWrVph3bp18PDwUB3JlOx2Ox555BEcP34cYWFhaNWqlepIprR+/XosXboU06dPR4cOHVTHIaIyYvEmogrj3LlzGDp0KAICAhAeHg4vLy/VkUxJCIHnnnsOv/zyC7799lvcd999qiOZUnp6OiZMmIB27drxHAcig2DxJqIKITc3F0OHDoXNZsPOnTtRq1Yt1ZFMa/bs2fjuu+8wdepUPPPMM6rjmJIQAuPHj0d6ejq2bNnCNdOJDILFm4iUczgceOyxx2CxWBAeHo4WLVqojmRay5Ytw7Rp0/Dkk0/i3XffVR3HtBYsWICffvoJH374Idq1a6c6DhGVExZvIlLutddeQ1hYGL766isMHDhQdRzT2rFjB8aNG4d+/frh22+/5UoyikREROD//u//cP/99+PVV19VHYeIyhFPjyYipb777jvMmzcPkydPxsSJE1XHMa2oqCgMHz4cTZo0wU8//cSpDYpkZWVh5MiRCAgIwH/+8x+uYkJkMBzxJiJlIiIi8MILL2DAgAH4+OOPVccxrQsXLuD++++Hp6cnNm3axA1aFHE4HHj88cdx9uxZ7NixA/7+/qojEVE5Y/EmIiWysrLw8MMPIyAgACtWrICrKx+OVLBarRgyZAjS09Oxe/du1KtXT3Uk03r77bcRHh6Or7/+Gr1791Ydh4gk4F86ItKdEAJPPfUUzp8/j927dyMgIEB1JFMqLi7GyJEjcfToUWzcuBEdO3ZUHcm0Vq5ciQ8++AATJ07EhAkTVMchIklYvIlId/PmzcOGDRswf/589OjRQ3UcUxJC4Pnnn0d4eDhCQ0MxePBg1ZFMa+/evRg3bhz69OmDTz/9VHUcIpKIZ20Qka6OHTuGt956C//4xz/w0ksvqY5jWv/+97+xcOFCvPXWWxg/frzqOKZlsVgwdOhQ1K5dGz/88ANPaiUyOBZvItJNUVERnnzySXh7e2PhwoVcrk6R5cuXY8qUKXj88ccxa9Ys1XFMKzExEYMGDYKrqyu2bNmCwMBA1ZGISDJONSEi3bzzzjs4cuQI1q9fb5p53UII1RH+ZOfOnXj66afRt29fLFq0qFIvV1fRju2tyMrKwuDBg5GRkYFdu3ahYcOGqiMRkQ5YvIlIFxEREfjggw8wbtw4DBs2THUc6SriaH5sbCxGjBiBxo0bY+3atahatarqSLelIh7bW1FYWIgHH3wQ0dHRCAsL40mtRCbC4k1E0jkcDkycOBH+/v745JNPVMcxpdzcXAwbNgxCCGzcuBE+Pj6qI5mSw+HAE088gV27dmHFihXcqZXIZFi8iUi6hQsX4sCBA1i6dCk3Z1HA6XRizJgxsFgs2LJlCxo1aqQ6kikJITB58mT8+OOPmDt3Lh577DHVkYhIZyzeRCRVcnIy3nrrLfTr1w+jR49WHceUZs6ciXXr1mHevHkYMGCA6jimNWfOHHz55Zf417/+hVdeeUV1HCJSoPKeVUNElcK0adOQl5eHBQsWVPq5uZXR+vXr8e6772Ls2LGYPHmy6jimtWjRIrz99tsYPXo0PvjgA9VxiEgRFm8iksZiseC7777Dc889h+bNm6uOYzpnz57F2LFj0alTJ3z99dd84qPIzz//jPHjx+Pee++t9CvJEFHZ8H8/EUnz1ltvoUaNGpg6darqKKZTVFSEUaNGwel0Ys2aNZV2BZPKbt++fRg5ciTat2/PDXKIiHO8iUiOvXv3Yv369Zg9e7Zp1uyuSN5++23s378fa9as4RrRikRHR2PIkCEICQlBWFgYatasqToSUaUihEBhYSGysrKQnZ2NrKysm768/vrruPvuu1X/E/6GxZuIpHjvvfcQEBDAecUKbNq0CR9//DGee+45PPLII6rjmNLFixf/tCtlUFCQ6khEShUVFSEjIwPp6elIT0+//P7VPnfl+zab7bq36+bmBh8fH3h7e1++1KtXD66uFbPiVsxURFSpHT58GOHh4Xj//fdRvXp11XFMJSMjA88++yxat26NuXPnqo5jSrm5uRgyZAhSU1Oxa9cuLt9IhlRYWIjU1FQkJycjJSXlmpfSAm21Wq95W25ubvD19YWfnx/8/PzQqFEjdO3aFb6+vvD19f1Tqf7rxcPD46rnr+zcuVPiv/72sXgTUbmbPXs2vLy88Pzzz6uOYjqTJ09Gamoqfv75Z3h4eKiOYzp2ux0jR47EkSNHsGHDBnTu3Fl1JKKbVlxcjJSUFCQlJSExMRFJSUlISkq6aqHOycm56m14enoiKCgIgYGBqFWrFlq3bn25UJeW6ytLtq+vL2rUqGGak79ZvImoXMXExOCnn37C1KlT4eXlpTqOqaxbtw7Lli3DjBkzuA25AkIITJw4EZs3b8Y333yD+++/X3UkosuysrJw/vz5PxXq0kvp55KTk+F0Ov/0fS4uLggICEBgYCACAwPRpUuXy+9f7VK9enXTlOjbweJNROVqwYIFcHNzw6RJk1RHMZXs7Gw899xzaNeuHaZMmaI6jim99957WLRoEaZNm4Z//vOfquOQiTidTiQnJ+PcuXM4d+4czp8/f/n90svVRqj9/f0REhKCkJAQtG/f/vL7ISEhuPPOOxESEoLAwEBUqVJFwb/KmFi8iajc5ObmYvHixRg5ciRPJtPZjBkzkJycjI0bN3LJOgXWrFmDGTNm4KmnnsLMmTNVxyEDKioqwtmzZ3H69GmcPn0ap06dwunTp3HmzBmcP3/+bychenl5oV69eqhfvz769u2LevXqoW7duqhduzZCQkJQq1YtLjOqAIs3EZWbZcuWIScnh6PdOjt27Bi++OILTJgwgXOKFTh8+DDGjh2LXr16ITQ0lC+z020TQiAxMRFRUVGwWCw4derU5YJ97ty5P00DqVmzJpo0aYL27dvjwQcfRL169S5f6taty6l+FRSLNxGVCyEEvvzyS3Ts2BHdunVTHcc0hBB44YUX4O3tjdmzZ6uOYzrJycn4xz/+AX9/f/z4448cQaSbIoRAfHw8oqKicOLECZw4ceLy+7m5uZev5+XlhSZNmqBbt2544okn0Lhx48uXgIAAPsmrhFi8iahcHDp0CFFRUfjmm2/4x0BHq1atwq+//oqFCxfC19dXdRxTsdvtGDFiBNLT0/H7779zehVdVXFxMSwWCyIiIi5fIiMjkZ2dffk6QUFBaNmyJcaMGYNWrVqhZcuWaNGiBcu1AbF4E5mYEKLcbmvZsmWoWrUqN2zRUVFREaZOnYp27dph3LhxquOYzpQpU/D7779j5cqV6NChg+o45a48Hx/Mwul04uTJk9i7dy/279+Pw4cP4+jRo5fnX3t6eqJt27Z47LHH0K5du8sl28/PT3Fy0guLN5FJlecoit1ux8qVK/HAAw/A29u73G6Xrm/hwoWIjY3Fpk2b4OLiojqOqYSFhV3eHXTUqFGq45QrjrDevJycHOzfvx/79u3D3r17sW/fPmRmZgK4NE2kU6dOmDRpEjp06IAOHTqgWbNmXCHE5Fi8iajM/ve//yE1NRVPPPGE6iimYbVa8d5776Fv374YNGiQ6jimEh8fjzFjxqB9+/bcHdRk8vPzsWfPHmzfvh3bt2/HoUOH4HA4oGkaWrZsiREjRqBHjx7o3r07mjdvzifE9Dcs3kRUZsuWLYOvry8GDx6sOopphIaGIjk5GT/99BNHKHXkcDjw+OOPo6ioCGvWrOHuoAbndDoRERGBsLAwbNu2Dfv27UNRURFcXV3RtWtXvPXWW+jTpw+6devGVUToprB4E1GZFBYWYsOGDRg9ejTXj9aJzWbD3Llz0a9fP/Ts2VN1HFP55JNP8Ntvv2Hp0qVo0qSJ6jgkQV5eHrZt24aNGzciLCwMFy5cgKZp6NixIyZPnoz+/fujd+/eqFGjhuqoVAmxeBNRmezcuRN5eXn4xz/+oTqKaSxZsgRJSUlYvHix6iimYrFYMH36dDz00EMYPXq06jhUjqxWKzZu3IhVq1Zhy5YtsNlsqFmzJgYNGoShQ4di8ODBCAgIUB2TDIDFm4jKZOPGjahWrRr69++vOoopOBwOfPjhh+jcuTMGDBigOo5pFBcXY+zYsahRowYWLFjA6T0GYLPZsGnTJqxatQobN25EQUEBQkJCMGHCBAwbNgx9+vThq3hU7li8iei2CSGwceNGDBw4kHNddbJp0yacPn0aa9asYfnT0fz583HgwAGsWrWK63VXcidOnMC3336LJUuWID09HQEBAXj66acxatQo9OrViydEklQs3kR0244ePYr4+HjMmDFDdRTT+PLLLxESEoLhw4erjmIaiYmJmDlzJoYOHYqRI0eqjkO3oaioCKtXr8YHH3yAqKgouLm54cEHH8QzzzyDe+65B66urEOkD/6mEdFt27BhAzRNw9ChQ1VHMYXTp09jy5YtmDlzJouCjl5//XXY7XbMnz+frzJUMllZWQgNDcVnn32GpKQk1KlTBx999BHGjBmDwMBA1fHIhPjITUS3bfPmzejcuTNfetfJ119/DVdXV/zzn/9UHcU0du/ejRUrVmDatGlo1KiR6jh0k9LS0vDhhx/iq6++gtVqxT333IPvvvsOVatWRb9+/VTHIxPjRCYiui1WqxUHDhzgCX46sdls+P777zF8+HDUqlVLdRxTcDgcePHFF1G3bl28+eabquPQTcjKysK0adPQoEEDfPzxxxg2bBgiIiKwdetWDBo0iK9YkHIc8Sai2/Lbb7+huLiYq5noJDw8HBkZGRg3bpzqKKaxcuVKHD16FKtWrUK1atVUx6HrsNvt+PLLLzFz5kxkZWVh5MiReOedd9CiRQvV0Yj+RPmIt6ZpwZqmfapp2hlN0wo1TUvWNG2jpmn3qM5GRNe2fft2uLu7cwMXnSxfvhyBgYF8hUEndrsdM2bMQPv27fHII4+ojkPX8b///Q/t2rXDK6+8gm7duiEyMhKrV69m6aYKSWnx1jStLYDjAF4C0BCADYA/gKEA/qdpGl/bI6qgduzYge7du3MkUAfZ2dnYuHEjHn30UZ5UqZNFixYhNjYWs2fP5vJyFVRGRgaeeOIJ3HvvvbDZbNiwYQPCw8PRrl071dGIrknZo4mmaZ4ANgDwAxABoLUQwguAD4BPAGgA5miadq+qjER0dVarFYcPH0bfvn1VRzGFH3/8ETabjbsl6qSoqAizZs1Cz549MXjwYNVx6CrWr1+Pli1bYvXq1Zg2bRqioqLwwAMPcA43VXgqh04mAKgHwArgASFEIgAIIXIA/EvTtEYAHgQwB8AvylIS0d8cPHgQTqcTPXr0UB3FFFasWIHGjRuja9euqqOYwsqVK5GQkIBvvvmGRa6CKSwsxMsvv4zQ0FC0a9cO4eHh6NChg+pYRDdN5etnpUM3K0pL9198VPK2o6ZpzXXKREQ3Ye/evQCAbt26KU5ifFlZWdi1axcefvhhlkAdOJ1OfPTRR2jTpg0GDRqkOg5dIS4uDr169UJoaChef/11HDhwgKWbKh0lI96aptUE0Knkwy3XuNo+ANkAvAD0BxCtQzQiugl79+5Fs2bN4OvrqzqK4W3evBnFxcV44IEHVEcxhU2bNiEqKgpLly7lE50KZMuWLRg1ahSAS9NMhg0bpjgR0e1RNeLdApfmcANA1NWuIIRwAogp+bClHqGI6MaEENi3bx+nmehk48aN8Pf356sLOpk3bx7q1KmDRx99VHUUKrFo0SIMGTIEdevWxR9//MHSTZWaJoTQ/4dq2j8ArCv58A4hRO41rrcWl+Z5/ySEGHG92+zcubM4dOhQ+Qa9gZkbo7DnxHl4e3vr+nNJjaysLEPd16dPn0ZycjJ69ep1S99XWFiA/fsPoGnTpobeyKWs93dKSgosFgu6du0CT8/bW/lFCIE9e/bAz98PzZtxxl2p/39su8LT07PMt1d6XxcUFODAgQNo0KAB6tatWw5JKyebrRD79u1Hs2bNEBwcrDTLuXPncPbsWfj4+KBVq1aoUqVKmW7PaI/jdHUtQ+5A35qpuPvuu3X/2Zqm/SGE6Hytr6s6ubL6Fe8XXOd6+SVva1zti5qmjQcwHgCCgoKwc+fOcgl3sxISbHA4HMjKytL155IaRruvbTYbhBC3/G/Kzs6+/L6RjsdflfX+zs+/9PCVk5MLm63otm7DarWiuLgYnh6ehj7WtyovLw8AkJOTA5vNVubbK72vk5KSoGlAtWrVTH28i4ou/b7m5+crPQ4XL15EcnIyfHx9UKd2HeTmXnWM7pYY7XGcri7BmQNrHbvuvfCmCCF0v+DSiZWi5OJ6nestL7nOlhvdZqdOnYQKO3bsUPJzSX9Gu68nTZokfH19b/n7pkyZIqpUqSIKCgokpKo4ynp/r1ixQgAQ0dHRt30br776qnB3dxe5ubllymI0y5cvFwBETExMudzejh07REFBgfDz8xMPP/xwudxmZXb27FkBQCxatEjJz3c6nWLq1KkCgHj22WeFw+Eot9s22uM4XZuq+xrAIXGdvqpqjrf1ivev9zph6euz1utch4h0FBkZiRYtWsDDw0N1FMPbunUrevfujRo1rvqiH5Wjn376Cenp6ZgwYYLqKKb30UcfYdasWXj22WcRGhrKDYzIUFT9Nidd8X7Ida5X+rULErMQ0S2IjIzkEl46SEtLw5EjR9C/f3/VUUxh6dKlqF+/Po+3Yv/973/xxhtv4NFHfjBy+AAAIABJREFUH2XpJkNS9RsdjUtTSACg1dWuoGmaC4BmJR+e0CMUEV1fSkoKkpKS0L59e9VRDK90biKLoHzZ2dn43//+h1GjRrHoKbR37148+eST6NWrFxYvXsz7ggxJyW+1uLSKSekSJAOvcbVuuLSGNwBskx6KiG4oMjISADjirYPt27ejevXq6Nz5mifHUznZtWsXHA7H5XWiSX8XL17E8OHDUbt2baxbt45T2ciwVD6dXFHydrSmaVdbk+xfJW//EELEXOXrRKSz48ePAwDatm2rOInx7dy5E3369IGbm5vqKIa3fft2tGjRgr/XijgcDjz++OPIycnBunXr4O/vrzoSkTQqi3cogHMAagL4WdO0lsClXS01TfsQwEMl15uiKB8R/YXFYkFAQAD8/PxURzG0jIwMWCwW9OnTR3UUw0tMTMTRo0cxatQo7lSpyKxZs7Bjxw4sWLAArVu3Vh2HSCpV63hDCFFQspHONgAdAURpmpaDS2t2u+DSHPApQohfVGUkoj+Ljo5G8+bcyEW2ffv2AQB3B9XBunXrIITgTpWKRERE4L333sMTTzyBsWPHqo5DJJ3SMxeEEEcAtAbwGYBYAFUBpAMIAzBQCPFvhfGI6C8sFgtatGihOobh7d27Fy4uLujSpYvqKIa3adMm1K5dG82aNbvxlalcFRUV4emnn0ZgYCA+++wz1XGIdKFsxLuUEOIigMklFyKqoNLS0pCens4Rbx3s27cPbdu25frdkhUUFGD79u0YMmSI6iim9NFHH+HIkSNYv349fHx8VMch0gXX6iGim2KxWACAI96SORwO7N+/n9NMdLBz504UFhaiW7duqqOYTkJCAt5//32MGDECw4YNUx2HSDcs3kR0U6KjowGAI96SxcTEIDc3l2VQB5s2bUK1atXQrl071VFM54033oDD4cDHH3+sOgqRrli8ieimREdHw9PTE3Xr1lUdxdAiIiIAAB07dlScxNiEENi0aRP69+8Pd3d31XFMZe/evVixYgVee+011K9fX3UcIl2xeBPRTTl9+jQaNWrE3eQki4iIQNWqVfnKgmSnTp1CbGwsBg8erDqK6UyfPh2BgYF44403VEch0h3/ghLRTYmNjUXDhg1VxzC8iIgItGnThhvnSLZjxw4AwMCB19o8mWTYs2cPtm7ditdee40nD5MpsXgT0Q0JIRAXF4cGDRqojmJoQghERESgffv2qqMY3q5du1CrVi00btxYdRRTmTlzJgICAvDcc8+pjkKkBIs3Ed1QWloa8vLyWLwli4+PR2ZmJjp06KA6iqEJIbBr1y707duXu1Xq6ODBg/jll1/w2muvoXr16qrjECnB4k1ENxQXFwcAnGoiWWRkJABwxFuy2NhYJCUloW/fvqqjmMrnn3+OmjVrYuLEiaqjECnD4k1ENxQbGwsAHPGWLCoqCgDQqlUrxUmMbffu3QCAu+66S3ES80hJScHq1asxduxY1KxZU3UcImVYvIlMTAhxU9crHfHm0l9yWSwWhISEwMvLS3UUQ9u3bx+8vb25GdQN3Ozjw81YuHAhioqK8Pzzz5fbbRJVRizeRCZ1K3Nb4+LiEBAQwFUIJLNYLGjZsqXqGIZ34MABdO3alfO7r6G8j0txcTG++uorDBw4kMtkkumxeBPRDcXFxXF+t2RCCERHR3MUVrL8/HwcO3YMXbt2VR3FNLZt24bExETO7SYCizcR3YSzZ89ymolkCQkJsFqtLN6SHT58GA6Hg8VbR8uXL4e3tzeGDBmiOgqRcizeRHRdQggkJiaidu3aqqMYmsViAQAWb8kOHDgAAOjSpYviJOaQn5+PtWvX4uGHH0bVqlVVxyFSjsWbiK4rKysLBQUFCAkJUR3F0Fi89XHw4EHUqVMHwcHBqqOYwoYNG2C1WjF69GjVUYgqBBZvIrqupKQkAMCdd96pOImxxcTEwNvbG4GBgaqjGNqRI0fQsWNH1TFMY8WKFbjzzju5dCNRCRZvIrquxMREAOCIt2SxsbFo1KgRV9qQqLCwEDExMWjbtq3qKKaQ8//au/cwue76vuOf7+7qtrrYli3L8l3BYBl8waCYQENRbMdtHYIDpaUhbXGSJ86TJiQ0FwolpW4JIdAkhAenSdyEPDVQ4pbUgMFcGicKF99iItmxZVmJdcHaXUm7knbXe9+d+fWPc87oaHZmd3bmfOfM5f16nvPszp5zdn77+83OfOY7v3PO+Li+9rWv6R3veId6eogbgETwBrCMJHhT8fZ18OBBzhzjbN++fSoWiwTvJvnGN76h+fl53XHHHXk3BWgZBG8AS0qmmmzbti3nlnSuQqGgI0eOELydPf3005Kk6667LueWdIcHH3xQ5513nt7whjfk3RSgZRC8ASxpYGBAmzdv1rp16/JuSscaHBzU3Nyctm/fnndTOtrTTz+ttWvX6qqrrsq7KR2vUCjooYce0u23366+vr68mwO0DII3gCUNDg4yzcTZwYMHJYmKt7Onn35a1157rXp7e/NuSsd78sknNTIywrm7gTIEbwBLGhgY4MBKZwTv5njmmWeYZtIkDz/8sCTp1ltvzbklQGsheANY0sDAABVvZwcPHlRPT48uv/zyvJvSsUZHR3X8+HHt2LEj76Z0hYcffljXX3+9tmzZkndTgJZC8AZQ1cLCgo4fP07F29nBgwd1+eWXa9WqVXk3pWMdOHBAknT11Vfn3JLONz09re985zu65ZZb8m4K0HII3gCqGh4eVrFY5Iwmzg4dOsSBlc6S4P2KV7wi55Z0vkceeUSzs7MEb6ACgjeAqoaHhyWJj4udvfjii7rsssvybkZHO3DggHp6evSyl70s76Z0vG9961vq6enRG9/4xrybArQcgjeAqgje/gqFgoaGhnTppZfm3ZSO9vzzz2v79u1avXp13k3peI8++qiuvfZabdq0Ke+mAC2H4A2gqpGREUkE70aEEJZcf+LECRUKBQ5grcNyfZt24MAB5nc3QbFY1GOPPcZFc4AqCN4Aqkoq3hdccEHOLWk/ZlbTdgMDA5JE8F6BWvs2USwWdeDAAeZ3N8G+ffs0Pj6u17/+9Xk3BWhJBG8AVSUV7/PPPz/nlnSuo0ePSiJ4exocHNTU1BTBuwkeeeQRSaLiDVRB8AZQ1fDwsDZv3swlnx0lFW/mePs5fPiwJC5Q1AyPPfaYLrjgAg5iBaogeAOoanh4mGkmzgYGBtTX16cLL7ww76Z0rCR4X3HFFfk2pAvs2bNHr33ta1c8HQjoFgRvAFWNjIxwYKWzgYEBbdu2TT09PB17OXLkiCSCt7e5uTk9++yzevWrX513U4CWxTM9gKqGh4cJ3s4GBgaY3+3syJEjuvDCC7Vu3bq8m9LR9u3bp/n5ed144415NwVoWQRvAFWNjIww1cTZwMAA87udHT58WFdeeWXezeh4e/bskSSCN7AEgjeAikIITDVpgqNHj1LxdnbkyBGmmTTB3r17tX79el111VV5NwVoWQRvABWNjo5qYWGB4O1oYmJCExMT2rZtW95N6VjFYpHg3SR79uzRDTfcwPEKwBL47wBQUXIOb6aa+EkuULR169acW9K5Tpw4odnZWaaaOAsh6JlnntF1112Xd1OAlkbwBlBREgqpePs5ceKEJPrYE6cSbI7h4WGdPn1a11xzTd5NAVoawRtARVS8/fHmxh+nEmyO/fv3SxLBG1gGwRtARaOjo5Kk8847L+eWdK4keHPxHD+Dg4OSxAGszp577jlJ0o4dO3JuCdDaCN4AKhobG5MknXPOOTm3pHMx1cTf0NCQ1qxZwxtIZ/v371d/fz+nxgSWQfAGUFFS8SZ4+xkeHta6deu0fv36vJvSsYaGhnTRRRdxCXNnzz33nHbs2MEZTYBl8B8CoKKxsTH19/dr1apVeTelY3FlUH9DQ0OcrrEJ9u/fzzQToAYEbwAVjY6OUu12duLECeZ3OxscHCR4O5uamtKRI0cI3kANCN4AKhobG9O5556bdzM6GhVvf0NDQ7r44ovzbkZHO3TokCTp5S9/ec4tAVofwRtARVS8/RG8fU1PT2t0dJSKt7ODBw9KkrZv355zS4DWR/AGUBEVb18hBKaaODt27JgkEbydJRVvgjewPII3gIqoePuanJzUzMwMFW9HQ0NDkgje3g4dOqT+/n4ey0ANCN4AKqLi7YurVvojeDfHoUOHtH37dk7ZCNSA4A10sRBC1XVjY2NUvB0lF89hqokfgndjlnp+SEuCN4DlEbyBLrVUdWpmZkazs7MEb0enTp2SJG3evDnnlnSuoaEh9fb28qnCCq2kch1CIHgDK0DwBrBIcrl4ppr4Sa4MyqXM/Rw/flxbtmzhaoqOTp06pZdeeongDdSIZyMAi3C5eH9JH/Pmxs/IyAjVbmec0QRYGYI3gEWoePsjePs7efKkzj///Lyb0dGOHDkiSbryyivzbQjQJgjeABah4u1vdHRUa9as0dq1a/NuSsciePsbHByUJF166aU5twRoDwRvAIskFW+Ct5/R0VHmdzsbGRnRBRdckHczOtrAwIBWrVrFGxygRgRvAIswDcLf6Ogo/esohKBTp04RCJ0NDg7q4osv5hzeQI0I3gAWoeLtj+Dta2xsTIVCgeDtbHBwUJdccknezQDaBsEbwCJjY2Pq6enRhg0b8m5KxyJ4+xoZGZEkppo4GxgY0MUXX5x3M4C2QfAGsMj4+Lg2bNjA+Y8dEbx9nTx5UpKoeDtLppoAqA2vqgAWmZyc1Pr16/NuRkcjePtKgjcVbz8TExMaHx9nqgmwAgRvAIsQvH2FEAjezpKpJlS8/SSnEqTiDdSO4A1gkampKYK3o5mZGc3NzRG8HTHVxB/BG1g5gjeARah4++J0jf5Onjyp3t5ezszjaGBgQJKYagKsAMEbwCIEb18Eb38jIyPavHkzBwg7ouINrBzPSAAWIXj7Inj743Lx/oaGhrR+/Xpt3Lgx76YAbYPgDWARgrcvgre/kydPckYTZyMjI9qyZUvezQDaCsEbwCIEb18Eb38nT57U5s2b825GRxseHubNDbBCBG8Ai0xOTqq/vz/vZnQsgre/sbEx+tcZFW9g5XIJ3ma2xcx+1sz+j5m9YGYzZjZpZs+Z2T1mdlUe7QIQnWOairevJHhzxg0/4+Pj2rRpU97N6GgjIyNUvIEV6svpfgfL7ntC0mpJO+Llp83sp0IIn8ujcUA3m5mZUQiB4O3opZdeUl9fn9auXZt3UzpSCEHj4+O8sXHGVBNg5fKaatIn6ZuS3iVpWwhho6R+ST8oaa+ktZLuM7Prc2of0LUmJyclieDtaHJyUhs2bMi7GR1rampKhUKBirej6elpTU5OMtUEWKG8gvebQghvCiHcF0I4JkkhhEII4TuSbpN0QlE4//c5tQ/oWlNTU5II3p4mJiboX0fj4+OSRPB2lFwZlIo3sDK5BO8QwjeXWDcs6aH45mub0yIACSre/qh4+xobG5PEHHpPw8PDkkTFG1ihVj2rycn4a2+urQC6EMHbHxVvX1S8/Y2MjEii4g2sVKsG7zfFX5/JtRVAFyJ4+6Pi7Yvg7S+peBO8gZVpueBtZndI2hnf/NM82wJ0I4K3Pyrevphq4i+peDPVBFiZvE4nWJGZXSLp3vjml0IIX1tm+7sk3SVJW7du1e7du30bWMHExEQu94vm67SxPnr0qBYWFhb9TU888YQkad++fZqens6hZa2h0fHet2+fpKg/jx07dta64eFhrV+/vqMeT82U9O3jjz+uoaGhRevTj+FTp04t+/s67X+7USdOnJAk7d+/v2q/PPnkk+rp6dHevXvV29s+s0IZ6+7RqmNdc/A2sw9K+mCd9/PREMIHlvn9GyR9QdKFko5I+unlfmkI4V7FQX3nzp1h165ddTavfrt371Ye94vm67SxfuCBB9TX17fobzp8+LAkadeuXdq+fXvzG9YiGh3vJGzfdNNN2rFjx1nrisWitm/f3lGPp2ZKwvbrXvc6XX311YvWP/XUU5Kk2267rabLxnfa/3ajXnzxRUnSjh07qvbL/fffr82bN+uWW25pYssax1h3j1Yd65VUvHtU/8GOS+5nZmslfVHRFJNhSf8khDBS530BaABTTfwxx9tXMtWEOd5+uFw8UJ+ag3cI4W5Jd2fdADNbLenzkm6WNCrpthDC81nfD4DaELz9Mcfb1/j4uPr7+9XX11KzKTvKyMiIzj///LybAbSdXA+uNLM+SZ+T9COKLht/ewhhb55tArpdErzXrVuXc0s60/z8vObm5gjejsbHx6l2Ozt9+nRN03gAnC234G1mPZL+p6S3SZqW9JYQwqN5tQdAZGpqSv39/erpabmTHnWE5I0NU038jI2NEbydjY2NcdYYoA65vLKamSk6KPKdkuYkvS2E8Fd5tAXA2SYnJ6nGOmIqj7/x8XFCoTOCN1CfvEpav6vorCULkv7lcqcNBNA8BG9fExMTkqh4e2Kqia8QAsEbqFPTjzwxs8slvSe+GST9kZn9UbXtQwgXNaVhACQRvL1R8fY3NjamrVu35t2MjjUxMaFisUjwBuqQxyHf6Sr7Kkk8OwItZHJyUv39/Xk3o2NR8fbHVBNfXBkUqF/Tg3cI4bAka/b9AqgNFW9fVLz9MdXEF8EbqB+nLQBwlunpaSrejqh4+yoWiwRvZ0nwPvfcc3NuCdB+CN4AzjIzM6O1a9fm3YyORcXb1+TkpEIIBG9HVLyB+hG8gS4WQlj0M4K3LyrevqampiTxxiYLlZ4fJII30AiCN9ClotPpL0bw9kXF21cSvJkuVb9qzw0JgjdQP4I3gLMQvH1NTEzIzLRu3bq8m9KRCN7+CN5A/QjeAM4yMzNDKHSUnDVmuaoi6kPw9jc2Nqbe3l76GKgDwRvAWah4+5qYmGCaiSOCt7/kqpW8eQRWjuANoGRhYUELCwsEb0eTk5McWOmI4O1vdHSUaSZAnQjeAEpmZmYkieDtiIq3L4K3v7GxMc7hDdSJ4A2ghODtb2pqiuDtiODtL5lqAmDlCN4ASgje/qanpzl41RHB2x/BG6gfwRtACcHbHwev+iJ4+yN4A/UjeAMoIXj7o+Lti+Dtj+AN1I/gDaCE4O2P4O1rampKvb29WrVqVd5N6UghBI2PjxO8gToRvAGUELyzF0I46zZTTbJT3rdSFLz7+/s5x7STmZkZFYtFTokJ1IngDaBkenpaEsE7C9WCHxXvxi0VqpPgDR+Tk5OSmMoD1IvgDaCEire/6elp+tcRwdtXErw5JSZQH4I3gBKCt68QgmZnZ6l4OyJ4+0oOXiV4A/UheAMoIXj7SvqX4O2H4O2LqSZAYwjeAEoI3r7oX38Eb19MNQEaQ/AGUDI3NydJWrNmTc4t6UzJwatUvP1MTU0RCh0x1QRoDMEbQMns7KwkafXq1Tm3pDMRvP1R8fbFVBOgMQRvACVJxZvg7YOpJv4I3r6YagI0huANoISKty8q3v4I3r6SqSb0MVAfgjeAkrm5OfX29qq3tzfvpnQkKt7+CN6+qHgDjSF4AyiZm5vjwEpHVLx9hRAI3s6Y4w00huANoGR2dpZpJo4I3r7m5uZULBYJhY6mpqa0Zs0aPhUD6kTwBlBCxdsXU018Mf/Y3+TkJNNMgAYQvAGUUPH2lRy8ypsbHwRvfwRvoDEEbwAlVLx9Ebx9Ebz9MYceaAzBG0AJFW9fBG9fBG9/VLyBxhC8AZTMzc0RvB0RvH0RvP0RvIHGELwBlMzOzhIKHRG8fRG8/THVBGgMwRtACRVvX3Nzc5K4MqiX5HSNnDXGDxVvoDEEbwAlHFzpa3Z2Vn19ferp4anXA6dr9EfFG2gMz/4ASji40hdTeXwxlccfFW+gMQRvACVUvH0RvH0RvP0RvIHGELwBlFDx9kXw9sVUE18hBKaaAA0ieAMo4eBKXwRvX1S8fc3MzCiEQMUbaADBG0AJwdAXnyj4Inj7mpyclCSCN9AAgjeAEirevphD7yuZakIf++A86UDjCN4ASgiGvvhEwdfs7KxWrVrF6RqdUPEGGsezE4ASpkL4Inj7on99EbyBxhG8gS4WQjjre6aa+CIY+pqZmeGMJhlKPz9IZ64Mum7dujyaA3QEgjfQpczsrNsLCwuSuJy5J4K3L/o3G+XPDYlkDj3BG6gfwRuAJGl+fl6StGrVqpxb0rmYyuNrdnaWircjzpMONI7gDUASwbsZOHjV18zMDP3riLPGAI0jeAOQRPBuBqZC+KJ/fVHxBhpH8AYgieDdDARDXxxc6YvgDTSO4A1A0png3dfXl3NLOhfB2xf96yu5MijBG6gfwRuApDNnNaHi7YfTNfoiePui4g00juANQBJTTZqB4O2LqSa+CN5A4wjeACQRvL2FEDQ/P0//OqLi7WtmZkY9PT1MRwMaQPAGIIng7a1QKCiEQMXbEcHbF58oAI0jeAOQRPD2Rv/6I3j7IngDjSN4A5BEMPRG//rjAkW+uDIo0DiCNwBJBENvc3NzksRUE0dzc3M8fh1R8QYaR/AGIIng7Y3+9cdZY3wRvIHGEbwBSCIYeqN//c3PzxO8Hc3MzDCVB2gQwRuAJIKhN6aa+CoUCioWi/SvIyreQOMI3gAkccl4b7yx8ZW8saF//RC8gcYRvAFI4pLx3gjevvhEwR9nNQEaR/AGIIlg6I1g6Ct5/NK/fqh4A40jeAOQRPD2Rv/64o2NPw6uBBpH8AYgiWDojf71xRxvf1ygCGgcwRuAJIKhNyqyvphq4o/zpAONI3gDkETw9kb/+uKNjT+CN9A4gjcASQRDb1RkfTHVxN/s7CyPX6BBBG8Akgje3giGvqh4+6PiDTSO4A1AEsHbG/3ri08UfIUQOLgSyEDLBG8z6zWzJ80sxMvdebcJ6Cbz8/Pq7e2VmeXdlI5EMPRFxdtXcoEt+hdoTMsEb0nvlvTavBsBdKuFhQUuF++IqSa+6F9fvLEBstESwdvMLpX0IUlHJB3PuTlAV5qfnye0OGKqiS8+UfA1Ozsrif4FGtUSwVvSJyVtkPSLkmZybgvQlQjevgiGvqjI+qJ/gWzkHrzN7C2SfkzSl0MIX8q7PUC3Inj7YiqEL/rXV9K/HFwJNCbX4G1m6yXdI2laUbUbQE4I3r6YauKLiqwv+hfIRt4V7w9JukzSb4YQDuXcFqCrLSwsEAodEbx9MZXHF3O8gWzkFrzN7EZFVe4Dkj6WVzsARBYWFtTb25t3MzrW3Nycenp66GMnVGR90b9ANnI5d5iZ9Uj6I0m9kn4hhDBX5++5S9JdkrR161bt3r07szbWamJiIpf7RfN12lgfPXpUhUKh9DcNDQ1pbm6uo/7GRjQ63s8++6wk6YknntDx48f1wgsvqLe3l/7NwL59+yRFfXvs2LGzfvb4449r48aNK/p9nfa/3ajh4WFJ0vPPP1/ql6R/0z9rR4x192jVsa45eJvZByV9sM77+WgI4QOp2z8v6fsl/e8Qwv+r83cqhHCvpHslaefOnWHXrl31/qq67d69W3ncL5qv08b6i1/8onp7e0t/0z333KONGzd21N/YiEbH+/jx6MyoN910k6655hp96Utf0tq1a+nfDCRh+6abbtKOHTskRSFckm6++WatX79+Rb+v0/63GzUwMCBJuvrqq0v9knxSs3PnzrbuK8a6e7TqWK+k4t2jqEJdj9J+ZnaxpN+QNCHpl+v8fQAyxgV0fM3Pz9O/jpjj7YupJkA2an4VCCHcLenuDO7zI5I2SfpPksbMbEPZ+uR61auTdSGEiQzuF8ASCoUCwdARB6/6SoIhj2EfHFwJZCOPgyuviL9+SNJLFZbL4/XvT/0MgDMOrvTFJwq+ktNhmtnyG2PFqHgD2cj7dIIAWgQVb18Eb1/0ry8uoANko+nBO4SwK4Rg1RZJR+JN/0vqZwAchBBK3xNcfNG/vphDn7308wMVbyAbVLyBLlX+kTxTTXwRvH0xhz47labrELyBbBC8AUhiqok3grcv+tcXB1cC2SB4A5BEcPFG//piqokvKt5ANlruWSqEcGXebQC6EVNNfBG8fTHVxBfBG8gGFW8Akphq4o3g7Yv+9cUFioBsELwBSKLi7SU5MwTBMHvps24w1cRXErzpY6AxBG8AkgiGWat01hj6NxuVzrrBVBNf8/Pz6unpUU8PsQFoBP9BACQx1cQbwdsX/esruTIogMYQvAFIYqqJN4KhL6aa+OITBSAbBG8Akqh4eyN4+yIY+uKNDZANgjcASQRDb/SvL4KhL6aaANkgeAOQxFQTbwRvX/SvL4I3kA2CNwBJTDXxxhsbX0w18UX/AtkgeAOQRMXQG/3ri6kmvuhfIBsEbwCSqMh6I3j7on99MdUEyAbBG4Akppp4Ixj6YiqEL4I3kA2CNwBJVLy9Ebx9MRXCF29sgGwQvAEohEDF2xnB2xfB0BdvbIBsELwBqFgsShIvrI4I3r7oX19MNQGyQfAGoIWFBUliqokjPlHwRUXWF8EbyAbBG4AKhYIkKt6eqMj6YqqJLx6/QDYI3gBKFW9eWP0QXHzRv76oeAPZIHgDYKpJExAMfTHVxBfBG8gGwRsAU02cFYtFhRDoX0dMNfFF8AayQfAGwFQTZ/SvPz5R8EX/AtkgeANgqokzgrevEALB0BkVbyAbBG8ATDVxRvD2lTx+CYZ+CN5ANgjeAKh4OyN4+6J//RG8gWwQvAEQXJzRv77m5+cl0b+emMoDZIPgDYCpJs4I3r6S/qUi64eKN5ANgjcAppo4I3j7ouLtj+ANZIPgDYCKtzOCty/61x/BG8gGwRsAwcUZ/euLqSa+OF0jkB2CNwCmmjgjePvi8euL0zUC2SF4Ayi9sBJcfBAMffHGxlcyh57gDTSO4A2A4O2MYOiLx68vgjeQHYI3ABWLRUkEFy8EQ1/0ry/eOALZIXgDKAWXnh6eEjzwxsYXwdsXwRvIDq+yAAiGznhj44vg7YvgDWSHVwEABENnBENf9K8cXAiVAAAR2UlEQVQvzvMPZIdXWaCLhRAkUfH2RjD0Rf/6SJ4fOCsPkB2CN9ClzKz0PRVvXwRDX0yFyFb6uUHi8QtkiVdZAFS8nRFcfNG/vnhjA2SH4A2Airez5I0N/euD4O2L/gWyw6sAACrezgguvuhfXxxcCWSH4A2AirczgqEv+tcXB1cC2eFVFgAVb2cEQ1/0ry8q3kB2CN4AqHg7Ixj6on99UfEGssOrLAAq3s44uNIXwdsXFW8gO7wKAKDi7Yxg6IvT3fmi4g1kh1dZAFS8nRG8fdG/vqh4A9kheAOg4u2MYOiL/vVFxRvIDq+yAAguzuhfX/SvL/oXyA7BGwBTTZzxiYIvgqEvppoA2eFVAADB0BlvbHwRvH0x1QTIDq+yAAiGzgiGvuhfX1S8gewQvAFQ8XZGMPRF//qi4g1kh1dZAFS8nREMfVGR9UX/AtkheAOg4u2M/vVFRdYX/Qtkh1cBAFS8ndG/vvhEwRcVbyA7BG8AVGSdEQx90b++qHgD2eFVFkCpIkvw9kEw9EX/+qJ/gezwKgtAhUKB0O2I4OKL/vWVVLyZagI0jldaACoWi4QWR0zl8UXw9kX/AtnhVQAAFW9nhUJBZiYzy7spHYn+9cXBlUB2eKUFQMXbGf3rq1AoEAodcXAlkB2CNwAq3s4KhQKhxdHCwgL964iKN5AdXmkBUJF1RvD2Rf/6ouINZIfgDYCKtzP61xfB2xcHBwPZ4b8IABVvZwRDX/Svr2QqDwevAo0jeAOgIuuMNza+CN6+OHgVyA6vtAAIhs4Ihr7oX18cvApkh+ANgIq3M4KhL/rXF/0LZIdXWgBUvJ3xxsYXUyF80b9AdnJ/JTCzC83sw2b2lJmNmdmEmf29mX3OzO7Iu31ANyAY+qJi6IupEL7oXyA7ub6FNbPbJX1G0nnxj6YlFSRdFS9bJH0xn9YB3YOKty+Cty/61xcVbyA7uZW4zOyNkh5QFLrvl3RdCKE/hLBR0vmS3ibpobzaB3QTKt6+eGPji+Dti4o3kJ1c3sKa2VpJfypptaR7Qwg/m14fQjilKJQDaAKCoS+CoS/61xcVbyA7eZW43iHpZZJOS/qVnNoAIEbF2xf964vg7YuKN5CdvF4JfiL++vkQwkRObQAQo+Lti2Doi/71RcUbyE7Tg7dF15x9XXzz22b2GjP7v2Y2bGYzZvYPZvb7ZnZFs9sGdCsqsr4Ihr7oX19UvIHs5PFKu1XSpvj7V0l6TNJbJfVLmlc0BeXfSXoqPgATgJMQgiQq3p5CCPSvk+TxS/D2Qf8C2bPkH6tpd2i2Q9Jz8c0gaUDST0p6OIQQzOwNig68fIWkE5J2hBBOV/ldd0m6K755taTnPdtexQWSRnK4XzQfY91dGO/uwVh3D8a6e+Q11leEELZUW1lz8DazD0r6YJ2N+GgI4QPx73mlpGdT63aFEP667L6uk7RXUUX+vSGE/1bn/bozsydDCDvzbgf8MdbdhfHuHox192Csu0erjvVKppr0SOptYEmkD6Z8qjx0S1II4e8k/UV889YVtBEAAABoSTUH7xDC3SEEq3N5X+pXHVd0dUpp6akhybrLVvg3AQAAAC2n6QdXhhBmJb2Q3KxlF8fmZOHevBuApmGsuwvj3T0Y6+7BWHePlhzrph9cKUlm9t8l/ZyiqSavrrLNNyT9sKSvhhBub2b7AAAAgKzldeLeT8dfbzCzXeUr44Mrb4lvPtSsRgEAAABecgneIYRHJT0Q3/y0md0SX1hHZvZ6SZ+P23ZE0qfyaCMAAACQpTwvVXenpCclXaroDCYTZjYu6RFF5/AeknRHCGEqtxZWYWYXmdknzOyF+Gqbx83sQTO7Zfm90SxmdrmZvScem++Z2ayZvWRmT5nZb5nZtmX2X21m7zWzvWY2YWajZvaomd2VvFFcZv9b4/s+ET9OXogfN1uz+ytRjZltMLMXzSzEy51LbMtYtyEz+z4z+7iZPReP21j8/afM7E1V9mGs24yZ9ZjZT5rZX8RXuZ6Px+1xM/uAmW1cYl/Gu0WY2UYze4uZfcjMvmpmI6nn5x017G/xuD0aj+NLZrbHzH7NzFbXsP9OM/szMxuMx/J7ZvbHZnZVDftuMrPfiJ9fpszspJk9bGZvr/XvLwkh5LZIWiXplxUF8HFJU4rO8f0RSVvybNsSbb5e0QnZQ7yMKTpLS5BUlPS+vNvIEqTobDjF1DglY7WQun1K0g9V2X9T/LhMtp2UNJu6/aCkviXu/wOpbQvxfSe3T0i6Nu8+6vRF0u+Vjf+djHXnLJJ+Kn7NSPp6ouz2HzPW7b8ouqr1wxWey9PP74clfR/j3dqLpB8rG8f0smOZfVdJ+kpq+9my//cnJG1YYv93Kbo6epLVRsueO25eYt9LJR1Mbf9S6ncFSX+won7IeyDaaZG0Lv4HD5L+VtKr4p9vkvTbqUG4Le+2dvsi6cr4n+vLkt4u6bz456sl/bPUP9GYpIsq7H9/vP6kpDdLMkXno3+XpOl43Yer3PftqcfCb0vaGP/8VZL2xD9/QdKavPupUxdJr1H0Juux1FjcWWVbxrrNFkn/SmeC1yeVCl2Stkr615J+irFu/0VRIS4JS++XdE7889Xx4+B0vP4vGe/WXhQF7+OKAvTdkn4m1cfLBe+PxttNx+PXG4/nm+PxDZI+W2Xf6yXNxdt8RnFhV9IVkr4R//y0KhR84/tIXkcOSXpD/PO1kn5NZwqvP1NzP+Q9EO20SHqPzrzbuaTC+gfi9d/Nu63dvkg6R9INS6zfkXri/c9l625MPRm8pcK+vxSvm5J0YYX1e+P1D1RYd2n8+AmS3p13P3XiomgK3d8oCt7psbyzwraMdZstki5U9GlVkPT+FezHWLfhouhYryDpT6qsvzM1rucx3q27SOotu31laoyqBm9JF0maibf7xQrr79CZN2fXV1j/hXj931RowwZJ34vX/06FfZMqfUHSqyus/3i8fkjS6lr6Ic853u3oJ+Kv/yuEMFBhfXJp+9fUMl8JfkIIYyGEp5ZYv1/Ru1hJem3Z6nfGX58PIXypwu73KqqUr5P0tvQKM3uVpBvimx+rcL9HJX0uvvkT5euRiXdL2qno4789y2zLWLefn5N0nqKLrH10Bfsx1u0pmUtd7X/5u6nv+1PfM94tJoRQWH6riv65pDWKxmvRublDCF+UdEBRdfqd6XVmdq6iTy8k6XfL2xBCmJD0h/HNH68w7z8Z378IIeyt0LZktsNFkm6u5Y8heNcoPngjCWhfr7LZY4oeGFKNA4BcnYy/9pb9/Ifir9+otFMIYVrSt+Kb5eOc7Dsm6fEq95s8fm4ysw21NRW1MLNLJH1I0ceZv17DLox1+0leCO8LIRRXsB9j3Z4Ox19vrLI+eV0+Lmkw9XPGu3Mk4/HNEMJMlW2ScS4fyx9UND88vU25ZCy3SbqmbN2usm3OEhdhn61y3xURvGt3jaJ3U9KZTj5L/CKQXOr+lc1oFOpjZn2S/lF885nUz03RNBSpyjjH9sVfy8c5uf3cEqEg2Td9X8jGJyVtlPSrIYSxpTZkrNuPmZ0v6eXxzW+b2c1m9nUzOx2faWCfRWcsuqBsP8a6ff2P+OtPmtn7zOwcqXS2knfozEf9vxriz/4Z746TjEctY3lNWdU62fdYCOGkKtuX+r70WDCzCyUlzyX1PI4qInjXLn3qucGqW51Zt+Sp6pC7n1f00VBR0n2pn2+StD7+vp5x3la2fql9K+2POpnZj0p6q6TdIYTP1LALY91+Xp76/jZFp6K9TWc+tbpG0n+QtNfM0pUrxrp9/Z6k31cUcD8iadTMRhUdo/NnkvYrmsOd/p9nvDvLSsZjQ7zUvG/86cdo2fbl32eW+wjetVuf+n56ie2S847z0VOLMrPrJf1mfPOeEEL6nWyj45zsX8u+lfZHHcxsvaR7FJ3i6edr3I2xbj/npr7/j4qqUK8LIWxS1L+3KzrN2yWS/jz+ZEtirNtWPCf3PZJ+RdEB01J08HySXzZK2lK2G+PdWRoZj1r2Te9fad9a77umxwHBu3bLnmgfrc+ii+Z8QdFBON9VVB07a5PU96Geu6hh33p+L5b2XyVdLunjIYR9y20cY6zbT/o1qyDprSGEJ6Roql8I4auKzu8tRdXvt8bfM9ZtyswukvQdSb8j6bOKDnrcoOjTj/dL+j5JnzKzj6R3S33PeHcOr7Fcbt9696+I4F27idT365bYLjmqemKJbZADM9us6OCK7ZL+XtKPVDhQIz1u/aqu2jhPlK2vJP0umsdJg8zs1YpODfaiogBeK8a6/aT78CshhH8o3yCE8BVFZziQpFsr7MdYt5f7JN2k6HSCd4YQng4hTIYQ/iGE8FuSfjbe7r1mdm38PePdWSbjr7WMpXT2eNQylun1lfat9b5rehwQvGuXnt9z8RLbJeuGHNuCFYoPyPm6pGsVnbPz1hDC8QqbjuvMP3k94zxYtn6pfSvtj5X7hKI5vh9QdFzVhvSS2m5N/LPkSZKxbj/p5+Hnq251Zt1l8VfGug2Z2Ssl/XB88+OVtgkhfFrRGap6FF1MRWK8O81KxmNCZwfgZfc1s3U6M40tPZYuuY/gXbv9OvNRw6sqbWBmPZKujm/W+nE3nMXzfx9SdG7nY4pC9/cqbRsfFf9cfLPiOMeSo5fLxzl9ZHW1/69k3/R9oX5XxF/vU3RRi/Il8Yfx7X0SY92mDurMXMtaPvoNEmPdxtIHyB5aYruD8dcrJca7AyXjUctYPpec3aZs34visyIttW96e4UQhiWNrOC+a8p9BO8ahRBekvRkfPOHq2z2OkUHfUjSw+6NwrLid7IPSnqDoqrIrSGEv19mt7+Kv1YcZzNbK+mN8c3ycU72PUfS91f5/bfFXx8PIUxW2QbNwVi3kfjUbrvjm0udwi0pgBxJ/Yyxbj/pU/ldvsR2yZvv9BttxrtzJOPxxnjcKknGuXwsv63ooHvpzNSzcslYDmnxm6jlHkeX6Ewory33ZXUp0G5YdOaS8eOStlVY/+fx+ifzbitLkKTVkr4aj8lpSa+pcb/0pYbfXGH9u1XbpYY/X2HdxfHjJ0j6hbz7qBuW1FjeyVi3/yLpx+M+nZd0VYX1P1JpTBnr9lsUHY+TjNmiy3nH2/xoapu3M97ts6i+S8Yv6u/UY2C5S8Y/LqmnbN16RW/Qg6TfrrBv+pLxN1RY/zvx+kHVeMn43Du+nRZFB1Uejjv5u5JeGf98o6LLyiYPoNvybmu3L4rm/H5eZ94o/cAK978/3ndE0u2p3/lv4yfqIOnDVfa9PfVY+JikjfHPXxk/boKkFyStybufumFJjcWdjHX7L4o+qX0y7tu/k/T9qZ//U0XTyZIXWWOs23tRdGxOEnw+ojgkKzqzyZ2KPskMiqairC7bl/FusUXRBWmSJf3m6AfK1pUH5I/qzBulfyOpNzVOI/G6z1a5zxskzcXb3Cfpgvjnl6ceX6clbamwrym6KnlQNKXpB+Kfr1F0istCvO5nau6DvAeh3ZZ4AJNBDoouKZt0fFHS+/JuI0uQpH+cGqPp+MW42vI3FfbfpDMv7kHRgTozqdsPSupb4v5/PbXtQvw4SW4PS7o27z7qliXV73dWWc9Yt9ki6dI49CT9nBxMl9zeL+kyxrr9F0UXJdmX6udkvNO3j0m6kfFu/aVs3JZarizbb5Wkr6TWz5T9zz+h+M1Rlft9l6JPyZKsNprad0LSzUvse6mi0J1s/1LqdwVJf7CiPsh7ENpxUfSxxyfiJ/4ZRRds+LKkW/JuG0tpjHat4B/8cJXfsVrxVfDif8wxSY9KuktllbQq+98aPy6G48fJC/HjZmve/dNNS2qc71xiG8a6zRZFFc+7FVW9J+Nx+1tFF9bZwFh3zqLo0+ZfkvTXiircSQj+rqJTiC6qVDLerbms4HX5ygr79ig6feSj8ThOSNoj6b2qYZqHohMs3K9oLvesojOc/YkqTFmrsO8mSR9WNAd8WtIpSX8p6V+stA8s/oUAAAAAHHFWEwAAAKAJCN4AAABAExC8AQAAgCYgeAMAAABNQPAGAAAAmoDgDQAAADQBwRsAAABoAoI3AAAA0AQEbwAAAKAJCN4AAABAExC8AQAAgCb4//PS0uTrSu5tAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import jax.numpy as jnp\n",
    "\n",
    "lm = [1, 2, 3, 4]\n",
    "M = len(lm)\n",
    "D = jnp.array(lm)\n",
    "a = jnp.min(lm)\n",
    "b = jnp.max(lm)\n",
    "t = jnp.linspace(-1, 6, 1000)\n",
    "u = 0.5 * jnp.ones(M)\n",
    "rho = 1\n",
    "def fun(lam):\n",
    "    return 1 + rho * jnp.sum(u**2/(D - lam))\n",
    "res = [fun(lam) for lam in t]\n",
    "plt.figure(figsize=(10,8))\n",
    "plt.plot(res, 'k')\n",
    "plt.plot(jnp.zeros_like(t))\n",
    "plt.ylim([-6, 6])\n",
    "plt.tight_layout()\n",
    "plt.yticks(fontsize=24)\n",
    "plt.xticks(fontsize=24)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- The function has only one root at $[d_i, d_{i+1}]$\n",
    "\n",
    "- We have proved, by the way, the **Cauchy interlacing theorem** (what happens to the eigenvalues under rank-$1$ perturbation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## How to find the root\n",
    "\n",
    "- A Newton method will fail (draw a picture with a tangent line).\n",
    "\n",
    "- Note that Newton method is just approximation of a function $f(\\lambda)$ by a linear function.\n",
    "\n",
    "- Much better approximation is the **hyperbola**:\n",
    "\n",
    "$$f(\\lambda) \\approx c_0 + \\frac{c_1}{d_i - \\lambda} + \\frac{c_2}{d_{i+1} - \\lambda}.$$\n",
    "\n",
    "- To fit the coefficients, we have to evaluate $f(\\lambda)$ and $f'(\\lambda)$ in the particular point.\n",
    "\n",
    "- After that, the approximation can be recovered from solving **quadratic equation**\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Important issues\n",
    "\n",
    "- First, stability: this method was abandoned for a long time due to instability of the computation of the eigenvectors.\n",
    "\n",
    "- In the recursion, we need to compute the eigenvectors of the $D + \\rho uu^*$ matrix.\n",
    "\n",
    "- The exact expression for the eigenvectors is just (let us check!)\n",
    "\n",
    "$$(D - \\alpha_i I)^{-1}u,$$ where $\\alpha_i$ is the computed root.\n",
    "\n",
    "- The reason of instability:\n",
    "    - if $\\alpha_i$ and $\\alpha_{i+1}$ are close, then the corresponding eigenvectors are collinear, but they have to be orthogonal\n",
    "    - if $\\alpha_i$ and $\\alpha_{i+1}$ are close, then they close to $d_i$, therefore matrices $D - \\alpha_i I$ and $D - \\alpha_{i+1} I$ are close to singular"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Loewner theorem\n",
    "\n",
    "- The solution came is to use a strange **Loewner theorem:**\n",
    "\n",
    "If $\\alpha_i$ and $d_i$ satisfy the **interlacing theorem** \n",
    "\n",
    "$$d_n < \\alpha_n < \\ldots < d_{i+1} < \\alpha_{i+1} \\ldots$$\n",
    "\n",
    "Then there exists a vector $\\widehat{u}$ such that $\\alpha_i$ are exact eigenvalues of the matrix\n",
    "\n",
    "$$\\widehat{D} = D + \\widehat{u} \\widehat{u}^*.$$\n",
    "\n",
    "- So, you first compute the eigenvalues, then compute $\\widehat{u}$ and only then the eigenvectors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Divide and conquer and the Fast Multipole Method\n",
    "\n",
    "In the computations of divide and conquer we have to evaluate the sums of the form\n",
    "\n",
    "$$f(\\lambda) = 1 + \\rho \\sum_{i=1}^n \\frac{|u_i|^2}{(d_i - \\lambda)},$$\n",
    "\n",
    "and have to do it at least for $n$ points. \n",
    "\n",
    "- The complexity is then $\\mathcal{O}(n^2)$, as for the QR algorithm.\n",
    "\n",
    "- Can we make it $\\mathcal{O}(n \\log n)$? \n",
    "\n",
    "- The answer is yes, but we have to replace the computations by the approximate ones\n",
    "by the help of **Fast Multipole Method**.\n",
    "\n",
    "Lets explain a little bit..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Few more algorithms\n",
    "\n",
    "- Absolutely different approach is based on the **bisection**.\n",
    "\n",
    "- Given a matrix $A$ its inertia is defined as a triple \n",
    "\n",
    "$$(\\nu, \\zeta, \\pi),$$\n",
    "\n",
    "where $\\nu$ is the number of negative, $\\zeta$ - zero and $\\pi$ - positive eigenvalues. \n",
    "\n",
    "- If $X$ is non-singular, then \n",
    "\n",
    "$$Inertia(A) = Inertia(X^* A X)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bisection via Gaussian elimination\n",
    "\n",
    "- Given $z$ we can do the Gaussian elimination:\n",
    "\n",
    "$$A - zI = L D L^*,$$\n",
    "\n",
    "and inertia of the diagonal matrix is trivial to compute.\n",
    "\n",
    "- Thus, if we want to find all the eigenvalues in the interval $a$, $b$\n",
    "\n",
    "- Using inertia, we can easily count the number of eigenvalues in an interval. \n",
    "- Illustration: if $Inertia(A)=(5,0,2)$ and after shift $Inertia(A-zI)=(4,0,3)$, $z\\in[a,b]$ then we know that $\\lambda(A)\\in[a,z]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Jacobi method\n",
    "\n",
    "- Recall what a Jacobi (Givens rotations) are\n",
    "\n",
    "- In a plane they correspong to a $2 \\times 2$ orthogonal matrix of the form\n",
    "\n",
    "$$\\begin{pmatrix} \\cos \\phi & \\sin \\phi \\\\ -\\sin \\phi & \\cos \\phi \\end{pmatrix},$$\n",
    "\n",
    "and in the $n$-dimensional case we select two variables $i$ and $j$ and rotate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Jacobi method (cont.)\n",
    "\n",
    "- The idea of the Jacobi method is to minimize sum of squares of off-diagonal elements:\n",
    "\n",
    "$$\\Gamma(A) = \\mathrm{off}( U^* A U), \\quad \\mathrm{off}^2(X) = \\sum_{i \\ne j} \\left|X_{ij}\\right|^2 = \\|X \\|^2_F - \\sum\\limits_{i=1}^n x^2_{ii}.$$\n",
    "\n",
    "by applying succesive Jacobi rotations $U$ to zero off-diagonal elements. \n",
    "\n",
    "- When the \"pivot\" is chosen, it is easy to eliminate it. \n",
    "\n",
    "- The main question is then what is the order of **sweeps** we have to make (i.e. in which order to eliminate).\n",
    "\n",
    "- If we always eliminate the largest off-diagonal elements the method has quadratic convergence.\n",
    "\n",
    "- In practice, a cyclic order (i.e., $(1, 2), (1, 3), \\ldots, (2, 3), \\ldots$) is used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Jacobi method: convergence\n",
    "\n",
    "- To show convergence, we firstly show that \n",
    "\n",
    "$$ \\text{off}(B) < \\text{off}(A), $$\n",
    "\n",
    "where $B = U^*AU$.\n",
    "\n",
    "- In this case we use the unitary invariance of Frobenius norm and denote by $p$ and $q$ the indices that is changed after rotation:\n",
    "\n",
    "$ \\Gamma^2(A) = \\text{off}^2(B) =  \\|B\\|^2_F - \\sum\\limits_{i=1}^n b^2_{ii} = \\| A \\|^2_F - \\sum\\limits_{i \\neq p, q} b^2_{ii} - (b^2_{pp} + b^2_{qq}) = \\| A \\|^2_F - \\sum\\limits_{i \\neq p, q} a^2_{ii} - (a^2_{pp} + 2a^2_{pq} + a^2_{qq}) = \\| A \\|^2_F - \\sum\\limits_{i =1}^n a^2_{ii} - 2a^2_{pq} = \\text{off}^2(A) - 2a^2_{pq} < \\text{off}^2(A)$\n",
    "\n",
    "- We show that the ''size'' of off-diagonal elements decreases after Jacobi rotation.\n",
    "\n",
    "- If we always select the largest off-diagonal element $a_{pq} = \\gamma$ to eliminate (pivot), then we have\n",
    "\n",
    "$$ |a_{ij}| \\leq \\gamma, $$\n",
    "\n",
    "thus  \n",
    "\n",
    "$$ \\text{off}(A)^2 \\leq 2 N \\gamma^2, $$\n",
    "\n",
    "where $2N = n(n-1)$ is the number of off-diagonal elements.\n",
    "\n",
    "- Or rewrite this inequality in the form\n",
    "\n",
    "$$2\\gamma^2 \\geq \\frac{\\text{off}^2(A)}{N}.$$\n",
    "\n",
    "Now we use relations $\\Gamma^2(A) = \\text{off}^2(A) - 2\\gamma^2 \\leq \\text{off}^2(A) - \\dfrac{\\text{off}^2(A)}{N}$ and get\n",
    "\n",
    "$$ \\Gamma(A) \\leq \\sqrt{\\left(1 - \\frac{1}{N}\\right)} \\text{off}(A). $$\n",
    "\n",
    "- Aften $N$ steps we have the factor\n",
    "\n",
    "$$\\left(1 - \\frac{1}{N}\\right)^{\\frac{N}{2}} \\approx e^{-\\frac{1}{2}},$$\n",
    "\n",
    "i.e. linear convergence. However, the convergence is locally quadratic (given without proof here)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Jacobi: summary\n",
    "\n",
    "Jacobi method was the first numerical method for the eigenvalues, proposed in 1846.\n",
    "\n",
    "- Large constant\n",
    "- Very accurate (high relative error for small eigenvalues)\n",
    "- Good parallel capabilities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary for this part\n",
    "- Many algorithms for the computation of the SEV solution:\n",
    "    - QR\n",
    "    - Divide-and-conquer\n",
    "    - Bisection\n",
    "    - Jacobi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Next lecture\n",
    "- We start **sparse and/or structured** NLA."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Questions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "    }\n",
       "    div.cell{\n",
       "        /*width:80%;*/\n",
       "        /*margin-left:auto !important;\n",
       "        margin-right:auto;*/\n",
       "    }\n",
       "    h1 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\n",
       "    h2 {\n",
       "        font-family: 'Fenix', serif;\n",
       "    }\n",
       "    h3{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "        margin-top:12px;\n",
       "        margin-bottom: 3px;\n",
       "       }\n",
       "\th4{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "       }\n",
       "    h5 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\t   \n",
       "    div.text_cell_render{\n",
       "        font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "        line-height: 1.2;\n",
       "        font-size: 120%;\n",
       "        /*width:70%;*/\n",
       "        /*margin-left:auto;*/\n",
       "        margin-right:auto;\n",
       "    }\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 50pt;\n",
       "\t\tline-height: 110%;\n",
       "        color:#CD2305;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\t\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 300;\n",
       "        font-size: 16pt;\n",
       "        color: #CD2305;\n",
       "        font-style: italic;\n",
       "        margin-bottom: .5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\n",
       "    \n",
       "    li {\n",
       "        line-height: 110%;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "\n",
       "</style>\n",
       "\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}