lectures/lecture16/lecture-16.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 16:  Iterative methods for large scale eigenvalue problems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Previous lecture\n",
    "\n",
    "- Matrix functions: matrix exponential, methods to compute matrix exponential, other matrix functions\n",
    "- Krylov methods for computing action of matrix functions\n",
    "- Matrix equations (Lyapunov equations, Sylvester equation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Partial eigenvalue problem\n",
    "\n",
    "- Recall that to find eigenvalues of matrix of size $N\\times N$ one can use, e.g. the QR algorithm.\n",
    "\n",
    "- However, in some applications matrix is so large, that we even can not store it exactly.\n",
    "\n",
    "- Typically such matrices are given as a **black-box** that is able only to multiply matrix by vector (sometimes even without access to matrix elements). This is what we assume today.\n",
    "\n",
    "- In this case the best we can do is to solve partial eigenvalue problem, e.g.\n",
    "\n",
    "    - Find $k\\ll N$ smallest or largest eigenvalues (and eigenvectors if needed)\n",
    "    - Find $k\\ll N$ eigenvalues closest to a given number $\\sigma$\n",
    "\n",
    "- For simplicity we will consider the case when matrix is normal and thus has orthonormal basis of eigenvectors. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Power method and related methods\n",
    "\n",
    "### Power method\n",
    "\n",
    "Recall that the simplest method to find the largest eigenvalue is the **power method**\n",
    "\n",
    "$$\n",
    "    x_{i+1} = \\frac{Ax_{i}}{\\|Ax_{i}\\|}\n",
    "$$\n",
    "\n",
    "The convergence is linear with rate $q = \\left|\\frac{\\lambda_1}{\\lambda_2}\\right|$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Inverse iteration\n",
    "\n",
    "To find the smallest eigenvalue one may run power method for $A^{-1}$:\n",
    "\n",
    "$$x_{i+1} = \\frac{A^{-1}x_{i}}{\\|A^{-1}x_{i}\\|}.$$\n",
    "\n",
    "To accelerate convergence <font color='red'>shift-and-invert</font> strategy can be used:\n",
    "\n",
    "$$x_{i+1} = \\frac{(A-\\sigma I)^{-1}x_{i}}{\\|(A-\\sigma I)^{-1}x_{i}\\|},$$\n",
    "\n",
    "where $\\sigma$ should be close to the eigenvalue we want to find."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Rayleigh quotient (RQ) iteration\n",
    "\n",
    "In order to get superlinear convergence one may use adaptive shifts:\n",
    "\n",
    "$$x_{i+1} = \\frac{(A-R(x_i) I)^{-1}x_{i}}{\\|(A-R(x_i) I)^{-1}x_{i}\\|},$$\n",
    "\n",
    "where $R(x_k) = \\frac{(x_i, Ax_i)}{(x_i, x_i)}$ is Rayleigh quotient. \n",
    "\n",
    "The method converges **cubically for Hermitian matrices** and quadratically for non-Hermitian case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Inexact inverse iteration framework\n",
    "\n",
    "- Matrices $(A- \\sigma I)$ as well as $(A-R(x_i) I)$ are ill-conditioned if $\\sigma$ or $R(x_i)$ are close to eigenvalues.\n",
    "\n",
    "- Thus, if you are not given e.g. LU factorization of such matrix you might face a problem.\n",
    "\n",
    "- In practice you can solve systems only with some accuracy. Recall also that condition number is an upper bound and is overestimated for cosistent rhs. So, even in RQ iteration  letting\n",
    "the shift tend to the eigenvalue [does not harm](http://www.sciencedirect.com/science/article/pii/S0024379505005756) significantly\n",
    "the performance of the iterative methods.\n",
    "\n",
    "- If accuracy of solution of systems increases from iteration to iteration, superlinear convergence for RQ iteration can still be present, see [Theorem 2.1](http://www.sciencedirect.com/science/article/pii/S0024379505005756).\n",
    "Otherwise, you will get linear convergence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Block power method\n",
    "\n",
    "The block power method (also known as subspace iteration method or simultaneous vector iteration) is a natural generalization of the power method for several largest eigenvalues.<br>\n",
    "It looks as:\n",
    "\n",
    "1. $Y_0$ is $N\\times k$ matrix of rank $k$, $Y_0 = X_0 R_0$ (QR-decomposition)\n",
    "2. $Y_i = AX_{i-1}$ \n",
    "3. $Y_i = X_i R_i$ (QR-decomposition)\n",
    "\n",
    "QR-decomposition plays role of normalization in the standard power method. \n",
    "\n",
    "Moreover, orthogonalization prevents the columns of the $X_i$ from converging all to the eigenvector corresponding to the largest modulus eigenvalue."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Accelerating convergence of the block power method\n",
    "\n",
    "* For Hermitian matrices convergence of the $j$-column is **linear** as for the power method with $q=\\frac{|\\lambda_{j}|}{|\\lambda_{j+1}|}$. \n",
    "\n",
    "\n",
    "* Hence, applying the block power method to the matrix $(A-\\sigma I)^{-1}$ will accelerate convergence (<font color='red'>shift-and-invert</font> strategy).\n",
    "\n",
    "\n",
    "* You can also accelerate the convergence by applying the **Rayleigh-Ritz procedure** discussed below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7ff0b30b3e10>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAA4b0lEQVR4nO3deXxU5fn//9eVBRICJEBYsxAgAYQAAcK+hd0iYBUV0Nb1U0TL1n7auvTX6ret1X60Le6KiqilAsUVd0F2qWwisieBAGHfErKQde7fH2cmGUI2yEwmmbmej8c85pwzZ865TgLzzn3uM/cRYwxKKaVURfw8XYBSSqm6TYNCKaVUpTQolFJKVUqDQimlVKU0KJRSSlUqwNMFuEN4eLiJiYnxdBlKKVWvbNu27awxpmXZ5V4ZFDExMWzdutXTZSilVL0iIofLW+5Vp55EZJKILMjMzPR0KUop5TW8KiiMMSuMMTNCQ0M9XYpSSnkNrwoKpZRSrueVfRRKqeorLCwkPT2dvLw8T5eiaklQUBCRkZEEBgZWa32vCgoRmQRMio2N9XQpStUb6enpNGnShJiYGETE0+UoNzPGcO7cOdLT0+nQoUO13uNVp560j0Kpq5eXl0eLFi00JHyEiNCiRYurakF6VVAopa6NhoRvudrftwaFk5V7TvHauoOeLkMppeoUDQonaw+c4Z8rD1BQZPN0KUr5lCeeeILu3bvTs2dPEhIS+O677zxSR0ZGBi+99NJVv+/xxx/nmWeeqfb6H3/8MU899dRV76cyaWlpxMfHu3SbDhoUTobGhZNbUMz3Ry54uhSlfMamTZv45JNP2L59Ozt37mTlypVERUW5bX9FRUUVvnatQXG1Jk+ezMMPP+z2/biKBoWTQZ1a4CewIeWsp0tRymecOHGC8PBwGjZsCEB4eDjt2rUDYNu2bYwYMYK+ffsyfvx4Tpw4AUBSUhLz5s1j8ODBxMfHs3nzZgA2b97M4MGD6d27N4MHD2b//v0ALFq0iFtvvZVJkyYxbtw4srOzGT16NH369KFHjx589NFHADz88MOkpqaSkJDAb3/7WwCefvpp+vXrR8+ePXnsscdK6n7iiSfo0qULY8aMKdlPWWfOnGHKlCn069ePfv36sXHjxpJ6Zs2aBUBqaioDBw6kX79+/PGPf6Rx48YATJ06lc8++6xkW3fffTfvvfceaWlpDBs2jD59+tCnTx++/fbbK/a7e/du+vfvT0JCAj179iQ5OflafjUl9PJYJ02DAukVFcaGlLP877guri1OqXrg/63YzZ7jF126zW7tmvLYpO4Vvj5u3Dj+9Kc/0blzZ8aMGcPUqVMZMWIEhYWFzJ49m48++oiWLVuydOlSfv/737Nw4UIAcnJy+Pbbb1m3bh333nsvu3btomvXrqxbt46AgABWrlzJo48+ynvvvQdYLZedO3fSvHlzioqK+OCDD2jatClnz55l4MCBTJ48maeeeopdu3axY8cOAL766iuSk5PZvHkzxhgmT57MunXrCAkJYcmSJXz//fcUFRXRp08f+vbte8WxzZ07l1/96lcMHTqUI0eOMH78ePbu3XvFOnPnzmX69Om88sorJcunTZvG0qVLmTBhAgUFBaxatYqXX34ZYwxff/01QUFBJCcnM3369CvGtnvllVeYO3cud9xxBwUFBRQXF1/T787Bq4LCGLMCWJGYmPiLa93GsNhwXlidQualQkKDq/dlFKXUtWvcuDHbtm1j/fr1rF69mqlTp/LUU0+RmJjIrl27GDt2LADFxcW0bdu25H3Tp08HYPjw4Vy8eJGMjAyysrK46667SE5ORkQoLCwsWX/s2LE0b94csL5L8Oijj7Ju3Tr8/Pw4duwYp06duqK2r776iq+++orevXsDkJ2dTXJyMllZWdx00000atQIsE4llWflypXs2bOnZP7ixYtkZWVdts6mTZv48MMPAbj99tv5zW9+A8BPfvIT5syZQ35+Pl988QXDhw8nODiYzMxMZs2axY4dO/D39+fAgQNX7HfQoEE88cQTpKenc/PNNxMXF1fJb6BqXhUUrjA0riXPfZPCptRzXB/fxtPlKFWrKvvL3538/f1JSkoiKSmJHj168NZbb9G3b1+6d+/Opk2byn1P2Us8RYQ//OEPjBw5kg8++IC0tDSSkpJKXg8JCSmZXrx4MWfOnGHbtm0EBgYSExNT7vcKjDE88sgj3H///Zctnz9/frUuMbXZbGzatIng4OAq1y0rKCiIpKQkvvzyS5YuXVoSjP/85z9p3bo1P/zwAzabjaCgoCvee/vttzNgwAA+/fRTxo8fz+uvv86oUaOuugYH7aMoIyEqjEYN/Nmo/RRK1Yr9+/dfdg59x44dtG/fni5dunDmzJmSoCgsLGT37t0l6y1duhSADRs2EBoaSmhoKJmZmURERABWP0BFMjMzadWqFYGBgaxevZrDh63RtZs0aXLZX/zjx49n4cKFZGdnA3Ds2DFOnz7N8OHD+eCDD7h06RJZWVmsWLGi3P2MGzeOF1544bJjK2vgwIElp8eWLFly2WvTpk3jzTffZP369YwfP76k9rZt2+Ln58c777xT7mmlgwcP0rFjR+bMmcPkyZPZuXNnhT+L6tCgKKNBgB8DO7bQDm2lakl2djZ33XUX3bp1o2fPnuzZs4fHH3+cBg0asHz5ch566CF69epFQkLCZR23zZo1Y/DgwcycOZM33ngDgN/97nc88sgjDBkypNLz8nfccQdbt24lMTGRxYsX07VrVwBatGjBkCFDiI+P57e//S3jxo3j9ttvZ9CgQfTo0YNbbrmFrKws+vTpw9SpU0lISGDKlCkMGzas3P0899xzbN26lZ49e9KtW7fL+iAc5s+fzz/+8Q/69+/PiRMncB5ZYty4caxbt44xY8bQoEEDAB588EHeeustBg4cyIEDBy5rKTksXbqU+Ph4EhIS2LdvH3feeWc1fhMVE2NMjTZQFyUmJpqa3Lho4YZD/OmTPWx4aCSRzRq5sDKl6p69e/dy3XXXebqMq5KUlMQzzzxDYmKip0upsdzcXIKDgxERlixZwrvvvltyFZY7lfd7F5FtxpgrfqjaR1GOYXHhAGxMOcvUftEerkYp5c22bdvGrFmzMMYQFhZWclVXXaJBUY7YVo1p3bQh65M1KJSqi9asWePpElxm2LBh/PDDD54uo1J1PihEJAR4CSgA1hhjFtfCPhkSG86a/Wew2Qx+fjpgmlLKd3mkM1tEForIaRHZVWb59SKyX0RSRMTx/fabgeXGmF8A5V+s7AbD4sI5n1PAnhOu/fKRUkrVN5666mkRcL3zAhHxB14EfgJ0A6aLSDcgEjhqX61mXy+8CkNirX4KvfpJKeXrPBIUxph1wPkyi/sDKcaYg8aYAmAJcCOQjhUWUEm9IjJDRLaKyNYzZ87UuMZWTYLo0roJG5I1KJRSvq0ufY8igtKWA1gBEQG8D0wRkZeB8r/VAhhjFhhjEo0xiS1btnRJQUPjwtmcdp68wlpryCjlk06ePMm0adPo1KkT3bp1Y8KECRw4cKBGQ2cvWrSI48ePu7hSy7XUNXjwYJfXkZSUdMU4T+5Ql4KivB5jY4zJMcbcY4x5oKqObBGZJCILMjMzXVLQ0LhwCopsbEkr2/hRSrmKMYabbrqJpKQkUlNT2bNnD3/961/LHXvpalxLUFQ2BHlNlTfKa31Rl4IiHXAehD4ScM+fA9U0oENzGvj76eknpdxo9erVBAYGMnPmzJJlCQkJV3zb2XloboCJEyeyZs0aiouLufvuu4mPj6dHjx7885//ZPny5WzdupU77riDhIQELl26VOmQ5Y8++igjRozg2WefrXC9bdu20atXLwYNGsSLL75Y4fFUNCy5Y/hwm83Ggw8+SPfu3Zk4cSITJkxg+fLlfP7559x2220l669Zs4ZJkyYB8MADD5CYmEj37t0v26ZDeT8DV6pLl8duAeJEpANwDJgG3H41G3DF6LHOGjUIoE/7MNYnn+URV2xQqbru84fh5I+u3WabHvCTiu/mtmvXrnKH6K6uHTt2cOzYMXbtsi6izMjIICwsjBdeeKHk29tVDVmekZHB2rVrKSwsZMSIEeWud8899/D8888zYsSIkntVlFXRsOTDhw8vWef9998nLS2NH3/8kdOnT3Pddddx7733MnbsWO6//35ycnIICQlh6dKlTJ06FbDufdG8eXOKi4sZPXo0O3fupGfPnpX+DFzJU5fHvgtsArqISLqI3GeMKQJmAV8Ce4FlxpjdlW2nNgyLa8meExc5m53v6VKUUuXo2LEjBw8eZPbs2XzxxRc0bdr0inX2799fMmR5QkICf/nLX0hPTy953fGBXNF6mZmZZGRkMGLECAB+/vOfl1uL87Dkffr0Yd++fVfcNGjDhg3ceuut+Pn50aZNG0aOHAlAQEAA119/PStWrKCoqIhPP/2UG2+8EYBly5bRp08fevfuze7duy8bury6P4Oa8EiLwhgzvYLlnwGflfdaddT0xkXlGRIbztNf7mdjylluTIhw2XaVqpMq+cvfXbp3787y5curXC8gIACbrfR+9o5hwZs1a8YPP/zAl19+yYsvvsiyZcuuGAbDGFPpkOWOgfUqWi8jI6Naw4pXNCx52XUqMnXqVF588UWaN29Ov379aNKkCYcOHeKZZ55hy5YtNGvWjLvvvvuKIdGr8zOoibrUR1FjxpgVxpgZzqMv1lSPiFBCgwO1n0IpNxk1ahT5+fm89tprJcu2bNnC2rVrL1svJiaGHTt2YLPZOHr0aMntT8+ePYvNZmPKlCn8+c9/Zvv27cDlQ4ZXNWS5Q0XrhYWFERoayoYNGwDrfhblqWhYcmdDhw7lvffew2azcerUqcuGI0lKSmL79u289tprJa2cixcvEhISQmhoKKdOneLzzz+/Yr8V/QxcpS71UdSYO1oU/n7C4E7WsOPGmGr9VaGUqj4R4YMPPmDevHk89dRTBAUFERMTw/z58y9bb8iQIXTo0IEePXoQHx9Pnz59AOvD+J577ilpbTz55JOAdY/pmTNnEhwczKZNm1i+fDlz5swhMzOToqIi5s2bR/ful9+oyTG0eXnrvfnmm9x77700atSo5N4QZY0bN469e/cyaNAgwOrA/te//kWrVq1K1pkyZQqrVq0iPj6ezp07M2DAgJKhxf39/Zk4cSKLFi3irbfeAqBXr1707t2b7t2707FjR4YMGXLFfiv6GbiKDjNeDcu2HuV3y3eyfOYgEmOau2y7StUF9XGY8fouOzubxo0bc+7cOfr378/GjRtp06Z276h5NcOMe9WpJ3e5oUdbQhr4s2TL0apXVkqpKkycOLHkEuA//OEPtR4SV0tPPVVDSMMAJie044Pvj/HHSd1oGhTo0u0rpXxLfRsm3ataFO7ozHaY1i+avEIbH+3w6HcAlXILbzwFrSp2tb9vrwoKd+oZGcp1bZuyZPMRT5eilEsFBQVx7tw5DQsfYYzh3LlzBAUFVfs9XnXqyZ1EhOn9o/jjR7v5MT2THpGub7Uo5QmRkZGkp6fjilGXVf0QFBREZGRk1SvaeVVQ1LiP4rtX4dQumPx8uS/f2CuCJz7dy5ItR+gR2ePaC1WqDgkMDKRDhw6eLkPVYV516qnGfRSZR+GHJVCYV+7LoY0CuaFHWz7acZzsfPeNMqmUUnWJVwVFjUUNhOICOLGjwlXuHhJDdn4RizYeqr26lFLKgzQonEUNsJ6P/LfCVXpGhjHmulYsWHeQi3mFtVSYUkp5jlcFRY1vXNS4JTTvBEc3V7ravDGduZhXxMIN2qpQSnk/rwoKl3yPInogHP0OKrlUMD4ilOu7t+GN9YfIzNVWhVLKu3lVULhE1ADIPQvnUitdbd7YOLLyi3ht/cFaKkwppTxDg6Ks6IHW89GK+ykAurZpyg092/LmxkOc05saKaW8mAZFWS3iILhZpR3aDr8a05m8IhvPrUqucl2llKqvNCjK8vODyP5WP0UVYls1Znr/KBZ/d4TUM9m1UJxSStU+DYryRA+Aswcg93yVq84b05mgQH/+9vm+WihMKaVqn1cFRY0vj3WIcvRTVN2qCG/ckAeSOvHVnlN8d/BczfarlFJ1kFcFhcuGGY/oA36B1eqnALhvaAfahgbxxGd7sdl0BE6llHfxqqBwmcBgaNurWi0KgKBAf347vgs70zP54Ptjbi5OKaVqlwZFRToMs76hnb6tWqv/NCGChKgwnvx8nw7toZTyKhoUFRkyD5pGwHv3Qt7FKlf38xP+8tN4zuXk88+vD7i/PqWUqiUaFBUJDoMpr0PGUfj015UO6eEQHxHKHQOieevbNPaeqDpclFKqPtCgqEz0AEh6BH78j3Wfimr4zbguhAYH8sePdumtJZVSXkGDoirDfg3th1qtiqNbqlw9rFEDHrq+K1vSLvCfrem1UKBSSrlXnQ8KEekoIm+IyHKPFODnD7cshCZtYPEUOLGzyrfclhjFgA7N+fMneziecakWilRKKfdxa1CIyEIROS0iu8osv15E9otIiog8XNk2jDEHjTH3ubPOKjVpDXd+BA2awDs3wZnKO6v9/ISnb+lFsTE89N5OPQWllKrX3N2iWARc77xARPyBF4GfAN2A6SLSTUR6iMgnZR6t3Fxf9YVFW2EhfvD2ZDj5Y6WrR7doxKMTrmN98ln+vflILRWplFKu59agMMasA8oOmNQfSLG3FAqAJcCNxpgfjTETyzxOV3dfIjJDRLaKyNYzZ8648CichMfCnR8CAm+Mh70rKl39jgHRDI0N54lP93LkXK57alJKKTfzRB9FBHDUaT7dvqxcItJCRF4BeovIIxWtZ4xZYIxJNMYktmzZ0nXVltW6O8xYDa2ug6U/g7VPg81W7qoiwt9u6UmAn3D/v7ZxqaDYfXUppZSbeCIopJxlFZ7EN8acM8bMNMZ0MsY8WemGXTUoYFWatIG7P4WeU2H1X6xTURfSyl01IiyY56b3Zt/Ji/xO+yuUUvWQJ4IiHYhymo8EjnugjpoJDIKbXoXJL8DxHfDSYNjyRrmti6Qurfjd+K6s+OE4C9bprVOVUvWLJ4JiCxAnIh1EpAEwDfjYFRt22eix1SUCfX4OD26CqH7Wdy1eH13uqLMzR3Tkhp5t+dsX+1i551Tt1KeUUi7g7stj3wU2AV1EJF1E7jPGFAGzgC+BvcAyY8xud9bhdmFR8PMPrRZG1klYOB6W3QVnS2+RKiI8fUtP4iNCeXDxdlbvq3Y/vVJKeZR40zlzEZkETIqNjf1FcrKH7mNdkAPfPg8bn4WiPOh+Ewz/rdX5DWTmFnLHG//lwMlsXv15X0Z2rTtXACulfJuIbDPGJF6x3JuCwiExMdFs3brVs0Vkn4FNL8Dm16AwB+LGwYD7oeMoMvOKS8Li+dt7M757G8/WqpRS+EhQ1IkWRVm552HzAqujO+c0tIiFvneTGXcTdy5N44f0TH49tjOzRsbi51feBWFKKVU7fCIoHOpEi6KsogLY86EVGulbwC+A4thxvJMzgKdSoxneLZp/TE2gccMAT1eqlPJRGhR1yel9sGMx7FwK2aco9A/mi4IENgcPY+LNP2NA1/aerlAp5YN8Iijq5KmnytiKIW0D7H6fwl0fEZh/gXwTwOGmiUQPupmgrmOheUdPV6mU8hE+ERQOdb5FUZ7iIvIPfsv2lf+m7clviBHruxamWQek0yjrHt4xwyAk3MOFKqW8lQZFPbL9yAVe/eBrWp/ewKSQPfSx7cG/KMd6sVU3iB4I0YMgaoA1qq1oJ7hSquY0KOoZm83w/vfH+L8v9nE+K4e7YzL4RVQ6rc9vg6OboSDLWrFxa4jsBxF9IaIPtE2w7vetlFJXySeCot71UVRDbkERb25M49W1qVzMK2JCjzbMTurIdX5H4eh3kL7VuorqfGrpm5p3hDY9oW1PaNML2sRbgaItD6VUJXwiKBy8oUVRVualQl5ff5A3N6aRnV/EuG6tmTM6jvgI+7hWuefhxA44tt16PrETMg6XbqBRuDVEeqtu1rfEW3WDll0gqKknDkcpVQdpUHiJzNxCFm48xMKNh8jKK2JE55b8cmQs/Ts0v3LlSxnWnfhO7YZTu6znM/ug0OkmSk3aQcvOEN4FwuMgvLP13KSttkCU8jEaFF7mYl4h72w6zMINhziXU0C/mGY8mBRLUpeWSGUf8Dab1dI4vQfO7IezB6zwOJsMBdml6wWGQItO1qO547mj9QhpqSGilBfyiaDwxj6KqlwqKGbJliO8tu4gxzPzuK5tUx5I6sSE+DYE+F/F4MDGQNYJKzzOpcC5VOv5fCpcOAzG6e58gSHQvAM0iyl9hLWHZu2tq7ACg118lEqp2uATQeHgCy2KsgqKbHy04xivrE0l9UwO7Vs0YsbwjkzpE0lQoH/NNl5cCBlH4PxBOH/Ier6QBhcOWc9FeZev37i1FRhh0RAaZQ3DHhptf46Ehk1qVo9Syi00KHyEzWb4as8pXl6Twg/pmbRs0pD7hnbgjgHRNAkKdMcOrcEOLxy2QiPjiHVqK+OI9chMB1vh5e8JCrUCpGmEFRyhEdDU8dzO6jcJDHJ9rUqpSmlQ+BhjDJtSz/HSmlQ2pJylSVAAdw5qzz1DOhDeuGHtFWKzQfZJyDgKF9Ot4Mg4ChePQeZRazov48r3NWphBUZT+6NJW2ja1r6srTUf3Ez7SpRyIQ0KH7YzPYOX16Tyxe6TNAzwY1q/aGYM70i7sDrSl1CQAxdPWEFy8bg9RI5ZfSYXj1nLcs9d+T7/htCktRUaTdpA4zal841bW48mbSC4Ofh54q6/StUvGhSKlNPZvLo2lQ++PwbATb0jmJnUiU4tG3u4smooyrduM5t1wv44aQVI9il7oJywpvMvXvlevwAIaQWNW9kDxGk6pGXpspCW1mkxbaUoH6VBoUocy7jEa+sOsmTLEfKLbEyIb8sDSZ1Kv7xXnxXk2MPjlHXKK+uUNZ992prPPm09ck6DsV35fv8GVmA4Ho1bWQMxliyzTzcKt6YDavE0nlJu5hNB4YuXx9bE2ex8Fm44xDubDpOVX8TILi2ZNSqWvu3L+fKet7EVW99mz3EEx5nSAMk+Y807pnPPQnFB+dtp2NTqTwkJt4dHC+u5ZJljvrk13bCJtlhUneUTQeGgLYqrk3mpkHc2pfHGhkNcyC1kYMfmzB4Vx+BOLSr/8p6vMMY6peUIjRxHkJyzz5+1PzvNl73Sy8Ev0B4eLezh0dzqQ3EESXCZZcHNIChM+1hUrdCgUFXKLSji398dYcG6g5zOyichKozZo2IZ1bWVBsbVMAbys6zQyD1vdcTnnIVLztMXrOncc9Y6ly5c/qXGy4g1InBwMytAgpuVhkjZR5BjvTBr2l9vrauqT4NCVVteYTHLt6Xz8ppUjmVcolvbpsweFcv47m3w89PAcAubDfIzS0Mj97wVLGWnS+YvWGN55WdWvt0GTUpDIzisdDoo1Gl5M2vesdzxmva/+BwNCnXVCottfLTjOC+tTuHg2RziWjVm1qhYJvZsh78GRt1QXAR5mfbgOG+FR15GaZBcMe30XHSp8m0HBJUGR0WPhk3LmbcvC2yk/TH1jAaFumbFNsOnP57ghW+SOXAqmw7hITyQ1ImbekcQeDXjSam6pTDP6ntxhEdeptO003z+xTLT9vmK+mEcxN8KDUd4NAwtM9/U6twvmbbPlyxrYrWI9PRZrdGgUDXmGB7k+W+S2X38IpHNgnkgqRO39I2kYUANx5NS9YsxUHipNETyMu0BklEaJhU+29fNz6qkX8ZJYKPSACl5NL18vkHjcuYbW+s5pgND9KKAKmhQKJcxxrB6/2meW5XCjqMZtA0NYuaITkztF1XzAQiV7zDGujeKIzTys6wQcUyXLLeHTH6207zTdEE22IqqsUMpDY3LnpuUMx9y+bKy0w1CrIeXnVqr10EhIj8FbgBaAS8aY76qbH0NitphjGFDylmeW5XMlrQLtGrSkPtHdOL2/tEEN9DAULXEGGsE45KwsT8Ksu2BcrF0uiDb6bUs+7Ic6x70jtfLjoZcISkNjJKHc5A0Luf1kHJea2y1mho0sqb93TB4Z3WPyFNBISILgYnAaWNMvNPy64FnAX/gdWPMU9XYVjPgGWPMfZWtp0FRu4wx/PfgeZ5blcymg+cIb9yAGcM7cseA9oQ01PPLqp4pLnQKlhxrusA+7QgTx+uFuaXBU5DrtG5uaQAV5FSzxWPnF+gUJo2s6cAQazqwUWnAlDcd2AiiB1lf/LwGngyK4UA28LYjKETEHzgAjAXSgS3AdKzQeLLMJu41xpy2v+/vwGJjzPbK9qlB4TmbD53n+W+SWZ98luYhDfifYR24c1AMjTUwlC8rKigNG8ej0Gn6snCxB5AjeK6Ydrw/t/wRA+7+DGKGXFOZHj31JCIxwCdOQTEIeNwYM94+/wiAMaZsSDjeL8BTwNfGmJUVrDMDmAEQHR3d9/Dhw64+DHUVth+5wHOrklmz/wxhjQK5b0gH7hoSQ1N33BNDKV9VXGQPnNzS8GneyepPuQYVBUWVlwCIiJ+IDL6mvVYsAjjqNJ9uX1aR2cAY4BYRmVneCsaYBcaYRGNMYsuWLV1XqbomfaKbseie/nz0yyH0jW7G378+wNCnvmH+ygNk5lZxWaVSqnr8A6zvrDRtC+Gx0LbXNYdEZao8H2CMsdlP+Qxy4X7Lu1SgwqaNMeY54LkqN1o6KGANSlOu1CsqjDfu7seuY5k8tyqZ+SuTeWP9Ie4eEsN9QzsQ1qiBp0tUSlWhuhcVfyUiU8R1A/6kA1FO85HA8Zpu1BizwhgzIzTUC4bL9jLxEaEsuDORz+YMY2hcOM9/k8LQv63m6S/3cSGngpFZlVJ1QrX6KEQkCwgBioFLWC0CY4xpWq2dXNlHEYDVmT0aOIbVmX27MWb3NRyD8350mPF6Yt/Jizz/TQqf/XiCRoH+/HxQDL8Y1oEWtXmbVqXUZTx51dO7QBIQDpwCHjPGvCEiE4D5WFc6LTTGPOGqfepVT/XHgVNZPP9NCp/sPE5QgD8/H9SeGcM71u59vZVSgAuCQkQmA8Pts2uMMZ+4sD6X0BZF/ZVy2gqMFT8cp2GAPz8bGM0vhnekVZMgT5emlM+oUVCIyFNAP2CxfdF0YJsx5mGXVuki2qKov1LPZPPiNyl8uOMYDQL8uGNAe+4f3pFWTTUwlHK3mgbFTiDBGOsmw/YvzH1vjOnp8kprQFsU3uPQ2RxesAdGgJ8wvX80DyR1orUGhlJu44qgSDLGnLfPN8c6/VSngsJBWxTeI+1sDi+uTuH974/h7yfc3j+amSM60SZUA0MpV6tpUEwD/gasxrriaTjwiDFmiasLdQUNCu9z5FwuL65O4b3t6fiJMK1/FA8kdaJtaLCnS1PKa1xzUIiIH3ALsB6rn0KA74wxJ91RaE3oqSfvd/R8Li+tSeE/W63AuK1fJA8mxdIuTANDqZqqaYtinTFmeJUr1hHaovB+VmCksnybNRLMbYlRPDgylggNDKWuWU2D4g9YX7RbCuQ4ljv6LOoaDQrfkX4hl5fXpLJsqxUYtyZG8WBSJyKbNfJwZUrVPzUNikPlLDbGmI6uKM7VNCh8z7GMS7y8JoWlW6zAuKVvFL8cqYGh1NWoaR/FrcaYpe4qzlW0j0Idz7jEy2tSWbrlKDZjuDXR6sOIaq6BoVRVtI9C+ZQTmVZgLNlsBcYtfSP55UgNDKUqo30UyiedzMzjlbWp/HvzEWw2w819Ipg1Mo7oFhoYSpWlfRTKp526mMfLa6zAKLYZbu4dwaxRsbRvEeLp0pSqMzx6K9Taon0UqiqnLtpbGN8dochmuKl3BLNGxhITroGh1DUFhYj8zhjzf/bpW40x/3F67a/GmEfdUm0NaYtCVeX0xTxeWXuQxd8dpshm+GlCBLNHaWAo33atQbHdGNOn7HR583WJBoWqrtMX83h13UH+9d/DFBbb+GnvCGaPiqODBobyQRUFRVW3QpUKpsubV6readU0iD9M7Mb6h0Zy75AOfPbjCUb/fQ2/XrqDg2eyPV2eUnVCVUFhKpgub16peqtVkyD+v4ndWP+7Udw3tAOf7TrBmH+s5VdLd5CqgaF8XFWnnoqxLocVIBjIdbwEBBljAt1e4TXQU0+qps5k5fPa+oO8vSmNgiIbk3u1Y9aoOGJbNfZ0aUq5jU9c9eSgQaFc5Wx2Pq+tO8jbmw6TV1TM5F7tmK2BobyUTwSFXh6r3KVsYEzq2Y45o2OJbdXE06Up5TI+ERQO2qJQ7nIuO5/X1h/i7U1pXCosZmLPdswZFUtcaw0MVf9pUCjlQudzCnht/UHe+tYKjBt6tGXO6Dg6a2CoekyDQik3cATG29+mkVtYzIQebZkzKo4ubTQwVP2jQaGUG13IKeD1DQdZtDGNnILSFoYGhqpPNCiUqgUXcgp4Y8MhFn2bRnZ+ERN6tGHO6Di6tmnq6dKUqpIGhVK1KCO3gNfXa2Co+kWDQikPyMgtYOGGQ7y5MY2s/CJ+Em8FxnVtNTBU3VNvg0JErgPmAuHAKmPMy1W9R4NC1TWZuYW8sfEQb244RFZ+Edd3twKjWzsNDFV3XOuggDXd6UIROS0iu8osv15E9otIiog8XNk2jDF7jTEzgduAKw5AqfogtFEgvx7bmQ0PjWLO6Dg2ppxlwnPruf+drew+nunp8pSqlFtbFCIyHMgG3jbGxNuX+QMHgLFAOrAFmA74A0+W2cS9xpjTIjIZeBh4wRjz76r2qy0KVddl5haycOMhFtpbGOO6tWbO6DjiI0I9XZryYR479SQiMcAnTkExCHjcGDPePv8IgDGmbEiUt61PjTE3VPDaDGAGQHR0dN/Dhw+75gCUcqPMS4W8ufEQb2w4RFZeEWO7tWauBobykIqCIsADtUQAR53m04EBFa0sIknAzUBD4LOK1jPGLAAWgNWicEGdSrldaHAg88Z05p4hHUoC4+s9pxhzXWvmjdHAUHWDJ4KivBseVfjBboxZA6yp1oZLBwW8psKU8hTnwFi0MY03Nhxk4vMaGKpucGtndgXSgSin+UjguCs2bIxZYYyZERqq/6lU/RQaHMjcMXFseHgU/zu2M1vSzjPx+Q38z1tb+DFdO72VZ3giKLYAcSLSQUQaANOAj12xYRGZJCILMjP1P5Sq35oGBTJ7dBwbHhrJb8Z1ZkvaBSa9sIH7Fm1hZ3qGp8tTPsbdVz29CyRhfQfiFPCYMeYNEZkAzMe60mmhMeYJV+5Xr3pS3iYrr5C3vk3jtfWHyLxUyKiurZg7Oo5eUWGeLk15kXr7hburoTcuUt4uK6+Qtzcd5rX1B8nILWRkl5bMHdOZBA0M5QI+ERQO2qJQ3i47v8jewrACI6lLS+aOjqN3dDNPl6bqMZ8ICm1RKF+TnV/E25vSeG3dQS5oYKga8omgcNAWhfI1ZQNjROeWzB0TRx8NDHUVNCiU8gE5+UW889/DLFh3kPM5BQzvbLUw+rbXwFBV84mg0FNPSlly8ov4138P86o9MIbFhTNvTBx92zf3dGmqDvOJoHDQFoVSltyCIt7ZZLUwztkDY+7oOBJjNDDUlTQolPJhuQX2FsZaKzCGxlotDA0M5UyDQilFbkERi/97hFfXpXI22wqMuWPi6KeBofCRoNA+CqWq51JBsb0PwwqMIbEtmDemswaGj/OJoHDQFoVS1XOpoJjF3x3mlbUHOZudz+BOVmD076CB4Ys0KJRSFbpUUMy/Nx/hlbWpnMnKZ1DHFswbE8eAji08XZqqRRoUSqkq5RUWs/i70sAY2LE588Z0ZqAGhk/wiaDQPgqlXCOvsJh/f3eEl+2BMaCDFRiDOmlgeDOfCAoHbVEo5Rp5hcW8u/kIL69J5bQGhtfToFBKXbOygdG/Q3PmjYljUMcWiJR3d2NVH2lQKKVqLK+wmCWbrVNSpy7aA2N0HIM6aWB4Aw0KpZTL5BUWs2zrUV5ancrJi3n0i2nGvDGdGayBUa9pUCilXK5sYCS2twJjSKwGRn2kQaGUcpv8omKWbTnKi06BMXdMHENjwzUw6hGfCAq9PFYpz8ovKmbZ1nReWp3Cicw8+rZvxtzRcQyL08CoD3wiKBy0RaGUZzkC4+XVKRzPzKN3dBjzxnRmuAZGnaZBoZSqdflFxfzH3sLQwKj7NCiUUh6TX1TM8m3pvLQ6lWMZl0iICmPemDhGdG6pgVGHaFAopTyuoMjG8m3pvLg6hWMZl+hlD4wkDYw6QYNCKVVnFBTZeG97Oi98Yw+MyFDmjelMUhcNDE/SoFBK1TnlBcbcMXGM7NJKA8MDNCiUUnVWQZGN97en88LqFNIvXKJnZChzR8cxqqsGRm2qKCj8PFHM1RKREBHZJiITPV2LUsr1GgT4Ma1/NKt/k8TfpvTgQm4B9721lckvbGTV3lN44x+09Ylbg0JEForIaRHZVWb59SKyX0RSROThamzqIWCZe6pUStUVgf5+TO0XzTf/m8T/TelJ5qXCksBYuUcDw1PceupJRIYD2cDbxph4+zJ/4AAwFkgHtgDTAX/gyTKbuBfoCYQDQcBZY8wnVe1XTz0p5R0Ki2188P0xXvgmhSPnc4mPaMrc0Z0Zc52eknIHj/VRiEgM8IlTUAwCHjfGjLfPPwJgjCkbEo73PwGEAN2AS8BNxhhbOevNAGYAREdH9z18+LDrD0Yp5RGFxTY+/P4Yz9sDo3u7pswdHcfYbq01MFyooqAI8EAtEcBRp/l0YEBFKxtjfg8gIndjtSiuCAn7eguABWC1KFxVrFLK8wL9/bg1MYqf9o7gw++P8cLqFGa8s41ubZsyd0wc4zQw3MoTndnl/Tar/GA3xiyq6rSTiEwSkQWZmZnXXJxSqu5yBMaqX4/gmVt7kVtQxP3vbOOG5zbwxa6T2Gz6N6I7eCIo0oEop/lI4LgrNmyMWWGMmREaGuqKzSml6qgAfz9u6RvJyl+P4O+39uJSYTEz/7WNG57XwHAHTwTFFiBORDqISANgGvCxKzasLQqlfEuAvx9T+kby9a+G84/bepFvD4wJz63ni10nNDBcxN1XPb0LJGFdtXQKeMwY84aITADmY13ptNAY84Qr96tXPSnlm4pthhU/HOe5VckcPJtD1zZNmDs6jvHd2+Dnp30YVfGJb2brjYuUUlB+YMwZHcf1GhiV8omgcNAWhVIKnALjm2QOnsmhS2srMH4Sr4FRHp8ICm1RKKXKU2wzfLLTamGknsmhc+vGzBkdx4T4thoYTnwiKBy0RaGUKk+xzfDpjyd4blUyKaeziWvVmLljNDAcNCiUUsquvMCYMzqOCT3a4u/DgeETQaGnnpRSV6PYZvjMHhjJ9sCYPTqOG3w0MHwiKBy0RaGUuho2m+GzXSd4dqUVGLGtGjN7VCwTe7bzqcDQoFBKqSrYbIbPd53k2VUHOHAqm04tQ5gzOs5nAkODQimlqskRGM+tSmb/qSw6tgxhzqg4JvXy7sDwiaDQPgqllCvZbIYvdluBse+k9weGTwSFg7YolFKuZLMZvtx9kmcdgREewuzRsUzq2Y4A/3pxR+lq0aBQSqkastkMX+05yfyVpYExa1Qsk3t5R2BoUCillItYgXGKZ1cls/fERTqEhzBrZCw3JtTvwPCJoNA+CqVUbbLZDF/vPcX8lVZgxLRoxKxRcfy0ngaGTwSFg7YolFK1yREYz65MZs+Ji7Rv0YhZI2O5qXdEvQoMDQqllHIzYwxf209J7T5e/wJDg0IppWqJMYaVe0/z7KoD7Dp2kejmjZg1ygqMwDocGBoUSilVy4wxrNp7mmdXJfPjscw6HxgaFEop5SHGGL7Zd5r5K63AiGoezOyRcdzUp24FhgaFUkp5mDGG1futwNiZbgXGrJGx3Nwnsk4Ehk8EhV4eq5SqDxyB8ezKZH5IzySyWWlgNAjwXGD4RFA4aItCKVUfGGNYs/8M81cl88PRDCLCgpk1KpYpHgoMDQqllKqjjDGsOXCG+StLA+OXI2O5pW/tBoYGhVJK1XHGGNbaA2OHPTAeHNmJW/tG1UpgaFAopVQ9YYxhXfJZ5q88wPdHMmgXGsSDI2O5NTGShgH+btuvBoVSStUzxhjWJ5/ln06B8cDIWG5zU2BoUCilVD3lCIz5Kw+w/UgGbUODeDCpE7f1i3JpYGhQKKVUPWeMYUPKWeavTGbb4Qu0DQ3igaRO3JYYRVBgzQOjoqDw/Dc8qiAiSSKyXkReEZEkT9ejlFKeIiIMi2vJ8pmD+Nd9A4gIC+aPH+0m6ek1vPVtGnmFxW7Zr1uDQkQWishpEdlVZvn1IrJfRFJE5OEqNmOAbCAISHdXrUopVV+ICEPjwvmPPTCimgfz2Me7GfH0aramnXf5/gJcvsXLLQJeAN52LBARf+BFYCzWB/8WEfkY8AeeLPP+e4H1xpi1ItIa+Adwh5trVkqpesERGENiW7Ap9Rwvr00lJjzE5ftxa1AYY9aJSEyZxf2BFGPMQQARWQLcaIx5EphYyeYuAA0relFEZgAzAKKjo2tStlJK1SsiwuDYcAbHhrtl+57oo4gAjjrNp9uXlUtEbhaRV4F3sFon5TLGLDDGJBpjElu2bOmyYpVSyte5+9RTeaScZRVeemWMeR94v1obLh0U8BpLU0opVZYnWhTpQJTTfCRw3BUbNsasMMbMCA0NdcXmlFJK4Zmg2ALEiUgHEWkATAM+dsWGRWSSiCzIzMx0xeaUUkrh/stj3wU2AV1EJF1E7jPGFAGzgC+BvcAyY8xuV+xPWxRKKeV67r7qaXoFyz8DPnP1/rSPQimlXK/OfzP7amiLQimlXM+rgkL7KJRSyvW8clBAETkDHL7Gt4cDZ11YTn3hi8fti8cMvnncvnjMcPXH3d4Yc8UX0bwyKGpCRLaWN3qit/PF4/bFYwbfPG5fPGZw3XF71aknpZRSrqdBoZRSqlIaFFda4OkCPMQXj9sXjxl887h98ZjBRcetfRRKKaUqpS0KpZRSldKgUEopVSkNCidXeYvWeklEokRktYjsFZHdIjLXvry5iHwtIsn252aertXVRMRfRL4XkU/s875wzGEislxE9tl/54O8/bhF5Ff2f9u7RORdEQnyxmMu71bTlR2niDxi/2zbLyLjr2ZfGhR2Trdo/QnQDZguIt08W5VbFAH/a4y5DhgI/NJ+nA8Dq4wxccAq+7y3mYs1EKWDLxzzs8AXxpiuQC+s4/fa4xaRCGAOkGiMice6xfI0vPOYFwHXl1lW7nHa/49PA7rb3/OS/TOvWjQoSpXcotUYUwAsAW70cE0uZ4w5YYzZbp/OwvrgiMA61rfsq70F/NQjBbqJiEQCNwCvOy329mNuCgwH3gAwxhQYYzLw8uPGGuw0WEQCgEZY97vxumM2xqwDzpdZXNFx3ggsMcbkG2MOASlYn3nVokFR6qpu0eoN7Pcz7w18B7Q2xpwAK0yAVh4szR3mA78DbE7LvP2YOwJngDftp9xeF5EQvPi4jTHHgGeAI8AJINMY8xVefMxlVHScNfp806AodVW3aK3vRKQx8B4wzxhz0dP1uJOITAROG2O2ebqWWhYA9AFeNsb0BnLwjlMuFbKfk78R6AC0A0JE5GeerapOqNHnmwZFKbfdorWuEZFArJBYbL8nOcApEWlrf70tcNpT9bnBEGCyiKRhnVIcJSL/wruPGax/0+nGmO/s88uxgsObj3sMcMgYc8YYUwi8DwzGu4/ZWUXHWaPPNw2KUm67RWtdIiKCdc56rzHmH04vfQzcZZ++C/iotmtzF2PMI8aYSGNMDNbv9RtjzM/w4mMGMMacBI6KSBf7otHAHrz7uI8AA0Wkkf3f+misfjhvPmZnFR3nx8A0EWkoIh2AOGBzdTeq38x2IiITsM5l+wMLjTFPeLYi1xORocB64EdKz9c/itVPsQyIxvrPdqsxpmxHWb0nIknAb4wxE0WkBV5+zCKSgNWB3wA4CNyD9Qei1x63iPw/YCrWFX7fA/8DNMbLjtl+q+kkrKHETwGPAR9SwXGKyO+Be7F+LvOMMZ9Xe18aFEoppSqjp56UUkpVSoNCKaVUpTQolFJKVUqDQimlVKU0KJRSSlVKg0LVCyJiROTvTvO/EZHHXbTtRSJyiyu2VcV+brWP4Lq6zPJ2IrLcPp1gv0zbVfsME5EHy9uXUtWlQaHqi3zgZhEJ93Qhzq5mBE7gPuBBY8xI54XGmOPGGEdQJQBXFRT2we8qEgaUBEWZfSlVLRoUqr4owrr/76/KvlC2RSAi2fbnJBFZKyLLROSAiDwlIneIyGYR+VFEOjltZoyIrLevN9H+fn8ReVpEtojIThG532m7q0Xk31hfXCxbz3T79neJyN/sy/4IDAVeEZGny6wfY1+3AfAnYKqI7BCRqSISYr/vwBb7wH432t9zt4j8R0RWAF+JSGMRWSUi2+37dox8/BTQyb69px37sm8jSETetK//vYiMdNr2+yLyhVj3Nfg/p5/HInutP4rIFb8L5Z0q+0tEqbrmRWCn44OrmnoB12ENx3wQeN0Y01+sGzbNBubZ14sBRgCdgNUiEgvciTX6aD8RaQhsFJGv7Ov3B+LtQzaXEJF2wN+AvsAFrA/xnxpj/iQio7C+Fb61vEKNMQX2QEk0xsyyb++vWEOO3CsiYcBmEVlpf8sgoKcx5ry9VXGTMeaivdX1XxH5GGsQwHhjTIJ9ezFOu/ylfb89RKSrvdbO9tcSsEYWzgf2i8jzWCORRtjv84C9HuUDtEWh6g37KLdvY92Yprq22O/BkQ+kAo4P+h+xwsFhmTHGZoxJxgqUrsA44E4R2YE1xEkLrDFyADaXDQm7fsAa+6B0RcBirHtCXKtxwMP2GtYAQVjDMwB87TQMhQB/FZGdwEqsIaRbV7HtocA7AMaYfcBhwBEUq4wxmcaYPKzxodpj/Vw6isjzInI94NWjDqtS2qJQ9c18YDvwptOyIux/9NgHgmvg9Fq+07TNad7G5f/+y45lY7A+fGcbY750fsE+XlROBfWVN5xzTQgwxRizv0wNA8rUcAfQEuhrjCkUa6TcoGpsuyLOP7diIMAYc0FEegHjsVojt2GNHaS8nLYoVL1i/wt6GVbHsEMa1qkesO5FEHgNm75VRPzs/RYdgf3Al8ADYg3Ljoh0FuvGP5X5DhghIuH2ju7pwNqrqCMLaOI0/yUw2x6AiEjvCt4XinXPjUJ7X0P7CrbnbB1WwGA/5RSNddzlsp/S8jPGvAf8AWvIcuUDNChUffR3rBEzHV7D+nDeDJT9S7u69mN9oH8OzLSfcnkd67TLdnsH8KtU0Qq331XsEWA18AOw3RhzNUNarwa6OTqzgT9jBd9Oew1/ruB9i4FEEdmK9eG/z17POay+lV1lO9GBlwB/EfkRWArcbT9FV5EIYI39NNgi+3EqH6CjxyqllKqUtiiUUkpVSoNCKaVUpTQolFJKVUqDQimlVKU0KJRSSlVKg0IppVSlNCiUUkpV6v8Hdh35by8Rr5YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import copy\n",
    "%matplotlib inline\n",
    "\n",
    "n = 100\n",
    "k = 10\n",
    "A = np.diag(1./(1. + np.arange(n))) # diagonal matrix with well-separated maximum eigenvalues\n",
    "A_clustered = np.diag(1 - 1./(1. + np.arange(n))) # diagonal matrix with clustered maximum eigenvalues\n",
    "\n",
    "def subspace_iter(A, Y0, num_iter=100):\n",
    "    Y0, _ = np.linalg.qr(Y0)\n",
    "    Y = Y0.copy()\n",
    "    Y_old = Y0.copy()\n",
    "    err = []\n",
    "    for i in range(num_iter):\n",
    "        X = A.dot(Y)\n",
    "        Y, _ = np.linalg.qr(X)\n",
    "        err.append(np.linalg.norm(Y_old - Y.dot(Y.T.dot(Y_old))))\n",
    "        Y_old = Y.copy()\n",
    "    return Y, err\n",
    "\n",
    "Y0 = np.random.random((n, k))\n",
    "Y, err = subspace_iter(A, Y0, num_iter=100)\n",
    "Y, err_clustered = subspace_iter(A_clustered, Y0, num_iter=100) #np.diag((diagonal - sigma)**(-1))\n",
    "plt.semilogy(err, label='Separated eigvals')\n",
    "plt.semilogy(err_clustered, label='Clustered eigvals')\n",
    "plt.xlabel('Number of iterations')\n",
    "plt.ylabel('Error')\n",
    "plt.legend(loc='best')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Before we go to advanced methods let us discuss the important concept of **Ritz approximation**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Ritz approximation\n",
    "\n",
    "Given subspace spanned by columns of unitary matrix $Q_k$ of size $N\\times k$ we consider the projected matrix $Q_k^* A Q_k$.\n",
    "\n",
    "Let $\\Theta_k=\\mathrm{diag}(\\theta_1,\\dots,\\theta_k)$ and $S_k=\\begin{bmatrix}s_1 & \\dots & s_k \\end{bmatrix}$ be matrices of eigenvalues and eigenvectors of $(Q_k^* A Q_k)$: \n",
    "\n",
    "$$(Q_k^* A Q_k)S_k = S_k \\Theta_k$$\n",
    "\n",
    "then $\\{\\theta_i\\}$ are called **Ritz values** and $y_i = Q_k s_i$ - **Ritz vectors**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Properties of the Ritz approximation\n",
    "\n",
    "- Note that they are not eigenpairs of the initial matrix $AY_k\\not= Y_k \\Theta_k$, but the following equality holds:\n",
    "\n",
    "    $$Q_k^* (AY_k - Y_k \\Theta_k) = Q_k^* (AQ_k S_k - Q_k S_k \\Theta_k) = 0,$$\n",
    "\n",
    "   so the residual for the Ritz approximation is **orthogonal** to the subspace spanned by columns of $Q_k$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- $\\lambda_\\min(A) \\leq \\theta_\\min \\leq \\theta_\\max \\leq \\lambda_\\max(A)$. Indeed, using Rayleigh quotient:\n",
    "\n",
    "    $$\\theta_\\min = \\lambda_\\min (Q_k^* A Q_k) = \\min_{x\\not=0} \\frac{x^* (Q_k^* A Q_k) x}{x^* x} = \\min_{y\\not=0:y=Q_k x} \\frac{y^*  A y}{y^* y}\\geq \\min_{y\\not= 0} \\frac{y^*  A y}{y^* y} = \\lambda_\\min(A).$$\n",
    "\n",
    "    Obviously, $\\lambda_\\min (Q_k^* A Q_k) = \\lambda_\\min(A)$ if $k=N$, but we want to construct a basis $k\\ll N$ such that $\\lambda_\\min (Q_k^* A Q_k) \\approx \\lambda_\\min(A)$.\n",
    "\n",
    "    Similarly, $\\theta_\\max \\leq \\lambda_\\max(A)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### <font color='red'>Rayleigh-Ritz method</font>\n",
    "\n",
    "Thus, if a subspace $V$ approximates first $k$ eigenvectors, then one can use the **Rayleigh-Ritz method**:\n",
    "\n",
    "1. Find orthonormal basis $Q_k$ in $V$ (e.g. by using QR decomposition)\n",
    "2. Calculate $Q_k^*AQ_k$\n",
    "3. Compute Ritz values and vectors\n",
    "4. Note that alternatevly one could use $V$ with no orthogonalization, but then generalized eigenvalue problem $(V^*AV)s_i = \\theta_i (V^*V)s_i$ has to be solved.\n",
    "\n",
    "The question is how to find a good subspace $V$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Lanczos and Arnoldi methods\n",
    "\n",
    "The good choice for $V$ is the Krylov subspace.\n",
    "\n",
    "Recall that in the power method we used only one Krylov vector \n",
    "\n",
    "$$x_k = \\frac{A^k x_0}{\\|A^k x_0\\|}.$$\n",
    "\n",
    "In this case $\\theta_k = \\frac{x_k^* A x_k}{x_k^* x_k}$ is nothing but a Ritz value. Natural idea is to use a bigger Krylov subspace.\n",
    "\n",
    "As a result we can find more eigenvalues (power method only gives $\\lambda_\\max$). Moreover,convergence of the eigenvalue corresponding to $\\lambda_\\max$ will be faster."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "For Hermitian matrices from the Arnoldi relation we have\n",
    "\n",
    "$$\n",
    "Q_k^*AQ_k = T_k,\n",
    "$$\n",
    "\n",
    "where $Q_k$ is orthogonal basis in the Krylov subspace generated by the Lanczos procedure and $T_k$ is triangular matrix.\n",
    "\n",
    "According to the Rayleigh-Ritz method we expect that eigenvalues of $T_k$ approximate eigenvalues of $A$. This method is called the **Lanczos method**. For nonsymmetric matrices it is called the **Arnoldi method** and instead of tridiagonal $T_k$ we would get upper=Hessenberg matrix.\n",
    "\n",
    "Let us show that  $\\theta_\\max \\approx\\lambda_\\max$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Why is $\\theta_\\max \\approx \\lambda_\\max$?\n",
    "\n",
    "Let us denote $\\theta_1 \\equiv \\theta_\\max$ and $\\lambda_1 \\equiv \\lambda_\\max$. Then\n",
    "\n",
    "$$\n",
    "    \\theta_1 = \\max_{y\\in \\mathcal{K}_i, y\\not=0}\\frac{(y,Ay)}{(y,y)} = \\max_{p_{i-1}} \\frac{(p_{i-1}(A)x_0, A p_{i-1}(A)x_0)}{(p_{i-1}(A)x_0, p_{i-1}(A)x_0)},\n",
    "$$\n",
    "\n",
    "where $p_{i-1}$ is a polynomial of degree not greater than $i-1$ such that $p_{i-1}(A)x_0\\not=0$.\n",
    "\n",
    "Expand $x_0 = \\sum_{j=1}^N c_j v_j$, where $v_j$ are eigenvectors of $A$ (form orthonormal basis).\n",
    "\n",
    "Since $\\theta_1 \\leq \\lambda_1$ we get\n",
    "$$\n",
    "    \\lambda_1 - \\theta_1 \\leq \\lambda_1 - \\frac{(p_{i-1}(A)x_0, A p_{i-1}(A)x_0)}{(p_{i-1}(A)x_0, p_{i-1}(A)x_0)}\n",
    "$$\n",
    "for any polynomial $p_{i-1}$. Hence\n",
    "$$\n",
    "\\lambda_1 - \\theta_1 \\leq \\lambda_1 - \\frac{\\sum_{k=1}^N \\lambda_k |p_{i-1}(\\lambda_k)|^2 |c_k|^2}{\\sum_{k=1}^N |p_{i-1}(\\lambda_k)|^2 |c_k|^2} =\n",
    "$$\n",
    "$$\n",
    "= \\frac{\\sum_{k=2}^N (\\lambda_1 - \\lambda_k) |p_{i-1}(\\lambda_k)|^2 |c_k|^2}{|p_{i-1}(\\lambda_1)|^2 |c_1|^2 + \\sum_{k=2}^N |p_{i-1}(\\lambda_k)|^2 |c_k|^2} \\leq \n",
    "(\\lambda_1 - \\lambda_n) \\frac{\\max_{2\\leq k \\leq N}|p_{i-1}(\\lambda_k)|^2}{|p_{i-1}(\\lambda_1)|^2 }\\gamma, \\quad \\gamma = \\frac{\\sum_{k=2}^N|c_k|^2}{|c_1|^2}\n",
    "$$\n",
    "\n",
    "Since the inequality holds for any polynomial $p_{i-1}$ we will choose a polynomial: \n",
    "\n",
    "$$|p_{i-1}(\\lambda_1)| \\gg \\max_{2\\leq k \\leq N}|p_{i-1}(\\lambda_k)|.$$\n",
    "\n",
    "This holds, e.g. for the Chebyshev polynomial on $[\\lambda_n,\\lambda_2]$. Thus, $\\theta_1 \\approx \\lambda_1$ or more precisely (Paige-Kaniel error bound, check it!):\n",
    "$$\n",
    "    \\lambda_1 - \\theta_1 \\leq \\frac{\\lambda_1 - \\lambda_n}{T_{i-1}^2(1 + 2\\mu)}\\gamma, \\quad \\mu = \\frac{\\lambda_1 - \\lambda_2}{\\lambda_2 - \\lambda_n},\n",
    "$$\n",
    "where $T_{i-1}$ is a Chebyshev polynomial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Demo: approximating largest eigenvalue with Lanczos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k=10, err = 0.1681667143962713\n",
      "k=20, err = 0.022651261043122872\n",
      "k=100, err = 6.365130644780947e-10\n"
     ]
    }
   ],
   "source": [
    "import scipy as sp\n",
    "import scipy.sparse\n",
    "from scipy.sparse import csc_matrix, csr_matrix\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.linalg\n",
    "import scipy.sparse.linalg\n",
    "import copy\n",
    "n = 40\n",
    "ex = np.ones(n)\n",
    "lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr')\n",
    "e = sp.sparse.eye(n)\n",
    "A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)\n",
    "\n",
    "def lanczos(A, m):\n",
    "    n = A.shape[0]\n",
    "    v = np.random.random((n, 1))\n",
    "    v = v / np.linalg.norm(v)\n",
    "    v_old = np.zeros((n, 1))\n",
    "    beta = np.zeros(m)\n",
    "    alpha = np.zeros(m)\n",
    "    for j in range(m-1):\n",
    "        w = A.dot(v)\n",
    "        alpha[j] = w.T.dot(v)\n",
    "        w = w - alpha[j] * v - beta[j] * v_old\n",
    "        beta[j+1] = np.linalg.norm(w)\n",
    "        v_old = v.copy()\n",
    "        v = w / beta[j+1]\n",
    "    w = A.dot(v)\n",
    "    alpha[m-1] = w.T.dot(v)\n",
    "    A = np.diag(beta[1:], k=-1) + np.diag(beta[1:], k=1) + np.diag(alpha[:], k=0)\n",
    "    l, _ = np.linalg.eigh(A)\n",
    "    return l\n",
    "\n",
    "# Approximation of the largest eigenvalue for different k\n",
    "l_large_exact = sp.sparse.linalg.eigsh(A, k=99, which='LM')[0][0]\n",
    "print('k=10, err = {}'.format(np.abs(l_large_exact - lanczos(A, 10)[0])))\n",
    "print('k=20, err = {}'.format(np.abs(l_large_exact - lanczos(A, 20)[0])))\n",
    "print('k=100, err = {}'.format(np.abs(l_large_exact - lanczos(A, 100)[0])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Practical issues and stability\n",
    "\n",
    "- The Lanczos vectors may loose orthogonality during the process due to floating-point errors, thus all practical implementations of it use **restarts**.\n",
    "\n",
    "- A very good introduction to the topic is given in the book of **Golub and Van-Loan (Matrix Computations)**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### More problems with the Lanczos method\n",
    "\n",
    "- Applying Lanczos directly to $A$ may result into a very slow convergence if $\\lambda_i\\approx \\lambda_{i+1}$ <br> (typically holds for smallest eigenvalues that are not well-separated)\n",
    "\n",
    "\n",
    "- To accelerate the convergence one may apply Lanczos to $(A-\\sigma I)^{-1}$, but in this case systems have to be solved **very accurately**. <br>\n",
    "Otherwise the Arnoldi relation does not hold anymore.\n",
    "\n",
    "An alternative to this approach are the so-called preconditioned iterative methods that include:\n",
    "1. PINVIT (Preconditioned Inverse Iteration)\n",
    "2. LOBCPG (Locally optimal block preconditioned CG)\n",
    "3. Jacobi-Davidson method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## PINVIT (preconditioned inverse iteration)\n",
    "\n",
    "### Derivation\n",
    "\n",
    "Consider Rayleigh quotient $R(x) = \\frac{(x,Ax)}{(x,x)}$. Then\n",
    "$$\n",
    "\\nabla R(x) = \\frac{2}{(x,x)} (Ax - R(x) x),\n",
    "$$\n",
    "\n",
    "so the simplest gradient descent method with a preconditioner $B$ reads\n",
    "\n",
    "$$\n",
    "    x_{i+1} = x_{i} - \\tau_i B^{-1} (Ax_i - R(x_i) x_i),\n",
    "$$\n",
    "\n",
    "$$\n",
    "    x_{i+1} = \\frac{x_{i+1}}{\\|x_{i+1}\\|}.\n",
    "$$\n",
    "\n",
    "Typically $B\\approx (A-\\sigma I)$, where $\\sigma$ is called shift.\n",
    "\n",
    "The closer $\\sigma$ to the required eigenvalue is, the faster the convergence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- Parameter $\\tau_k$ is chosen to minimize the $R(x_{i+1})$ over $\\tau_k$ (steepest descent method).\n",
    "\n",
    "- One can think of this minimization procedure as minimization in basis $V = [x_i, r_i]$, where $r_{i}=B^{-1} (Ax_i - R(x_i) x_i)$.\n",
    "\n",
    "- This results into the generalized eigenvalue problem $(V^*AV)\\begin{bmatrix}1 \\\\ -\\tau_i \\end{bmatrix} = \\theta (V^*V) \\begin{bmatrix}1 \\\\ -\\tau_i \\end{bmatrix}$ (Rayleigh-Ritz procedure with no orthogonalization of $V$). Here $\\theta$ is the closest to the required eigenvalue."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Convergence theory\n",
    "\n",
    "**Theorem** ([Knyazev and Neymeyr](http://www.sciencedirect.com/science/article/pii/S002437950100461X)) \n",
    "\n",
    "Let \n",
    "- $R(x_{i})\\in [\\lambda_j,\\lambda_{j+1}]$\n",
    "- $R(x_{i+1})\\in [R(x_{i}),\\lambda_{j+1}]$ (case $R(x_{i+1})\\in [\\lambda_{j}, R(x_{i})]$ is similar)\n",
    "- $\\|I - B^{-1} A\\|_A \\leq \\gamma < 1$\n",
    "\n",
    "then\n",
    "\n",
    "$$\n",
    "\\left|\\frac{R(x_{i+1}) - \\lambda_j}{R(x_{i+1}) - \\lambda_{j+1}}\\right| < \\left[ 1 - (1-\\gamma)\\left(1 - \\frac{\\lambda_j}{\\lambda_{j+1}}\\right) \\right]^2 \\cdot \\left|\\frac{R(x_{i}) - \\lambda_j}{R(x_{i}) - \\lambda_{j+1}}\\right|\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Block case\n",
    "\n",
    "To find, e.g. $k$ eigenvalues one can do a one step of PINVIT for each vector:\n",
    "\n",
    "\n",
    "$$\n",
    "    x^{(j)}_{i+1} = x^{(j)}_{i} - \\tau^{(j)}_i B^{-1} (Ax^{(j)}_i - R(x^{(j)}_i) x^{(j)}_i), \\quad j=1,\\dots,k\n",
    "$$\n",
    "\n",
    "$$\n",
    "    x^{(j)}_{i+1} = \\frac{x^{(j)}_{i+1}}{\\|x^{(j)}_{i+1}\\|}.\n",
    "$$\n",
    "\n",
    "And then orthogonalize them using the QR-decomposition. However, it is better to use the Rayleigh-Ritz procedure:\n",
    "\n",
    "- Set $X^{i}_k = [x^{(1)}_{i},\\dots, x^{(k)}_{i}]$ and $R^{i}_k = [B^{-1}r^{(1)}_{i},\\dots, B^{-1}r^{(k)}_{i}]$, where $r^{(j)}_{i} = Ax^{(j)}_i - R(x^{(j)}_i) x^{(j)}_i$\n",
    "\n",
    "\n",
    "- Set $V = [X^{i}_k, R^{i}_k]$, use Rayleigh-Ritz procedure for $V$ to find new $X^{i+1}_k$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## LOBPCG (Locally Optimal Block Preconditioned CG)\n",
    "\n",
    "### Locally optimal PCG (not \"Block\" so far :))\n",
    "LOPCG method\n",
    "\n",
    "$$\n",
    "    x_{i+1} = x_{i} - \\alpha_i B^{-1} (Ax_i - R(x_i) x_i) + \\beta_i x_{i-1} ,\n",
    "$$\n",
    "\n",
    "$$\n",
    "    x_{i+1} = \\frac{x_{i+1}}{\\|x_{i+1}\\|}.\n",
    "$$\n",
    "\n",
    "\n",
    "is a superior to PINVIT method as it adds to basis not only $x_i$ and $r_i$, but also $x_{i-1}$.\n",
    "\n",
    "However, this interpretation leads to an unstable algorithm as $x_{i}$ is becoming colinear to $x_{i-1}$ as the procedure converges."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### LOPCG (stable version)\n",
    "\n",
    "Knyazev suggested an equivalent stable version, which introduces new vectors $p_i$ (conjugate gradients)\n",
    "\n",
    "$$\n",
    "p_{i+1} = r_{i} + \\beta_i p_{i},\n",
    "$$\n",
    "\n",
    "$$\n",
    "x_{i+1} = x_{i} + \\alpha_i p_{i+1}.\n",
    "$$\n",
    "\n",
    "One can check that $\\mathcal{L}(x_{i},x_{i-1},r_{i})=\\mathcal{L}(x_{i},p_{i},r_{i})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The stable version explains name of the method:\n",
    "\n",
    "In standard CG method we would minimze Rayleigh quotient $R$ in the conjugate gradient direction $p_{i+1}$: \n",
    "\n",
    "$$\\alpha_i = \\arg\\min_{\\alpha_i} R(x_i + \\alpha_i p_{i+1}).$$\n",
    "\n",
    "In the locally-optimal CG we minimize over two parameters:  \n",
    "\n",
    "$$\\alpha_i, \\beta_i = \\arg\\min_{\\alpha_i,\\beta_i} R\\left(x_i + \\alpha_i p_{i+1}\\right) = \\arg\\min_{\\alpha_i,\\beta_i} R\\left(x_i + \\alpha_i (r_{i} + \\beta_i p_{i})\\right)$$\n",
    "\n",
    "and we locally obtain more optimal solution. That is why the method is called **locally optimal**.\n",
    "\n",
    "As for PINVIT coefficients $\\alpha_i,\\beta_i$ can be found by the Rayleigh-Ritz procedure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Locally optimal <font color='red'> block </font> PCG\n",
    "\n",
    "In the block version similarly to PINVIT on each iteration we are given basis $V=[X^{(i)}_k,B^{-1}R^{(i)}_k, P^{(i)}_k]$ and use Rayleigh-Ritz procedure.\n",
    "\n",
    "The overall algorithm:\n",
    "\n",
    "1. Find $\\tilde A = V^* A V$\n",
    "2. Find $\\tilde M = V^*V$\n",
    "3. Solve generalized eigenvalue problem $\\tilde A S_k = \\tilde M S_k \\Theta_k$\n",
    "4. $P^{(i+1)}_{k} = [B^{-1}R^{(i)}_k, P^{(i)}_k]S_k[:,k:]$\n",
    "5. $X^{(i+1)}_{k} = X^{(i)}_k S_k[:,:k] + P^{(i+1)}_{k}$ (equivalent to $X^{(i+1)}_{k} = VS_k$)\n",
    "6. Calculate new $B^{-1}R^{(i+1)}_k$\n",
    "7. Set $V=[X^{(i+1)}_k,B^{-1}R^{(i+1)}_k, P^{(i+1)}_k]$, goto 1.\n",
    "\n",
    "**Deflation technique** which stops iterating converged eigestates can also be applied here.\n",
    "\n",
    "The method also converges linearly, but faster than PINVIT."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### LOBPCG summary\n",
    "\n",
    "- Locally optimal preconditioned solver\n",
    "\n",
    "- Linear convergence\n",
    "\n",
    "- Preconditioner $(A-\\sigma I)$ is not always good for eigenvalue problems\n",
    "\n",
    "The next method (Jacobi-Davidson) has smart preconditioning and superlinear convergence (if systems are solved accurately)!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Jacobi-Davidson (JD) method\n",
    "\n",
    "Jacobi-Davidson method is a very popular technique for solving eigvalue problems (not only symmetric!).\n",
    "\n",
    "It consits of two **key ingredients**:\n",
    "\n",
    "- Given a preconditioner for $A-R(x_j) I$ it automatically constructs a good preconditioner for the eigevalue problem:\n",
    "$$\n",
    "    B = (I - x_j x^*_j) (A - R(x_j) I) (I - x_j x^*_j),\n",
    "$$\n",
    "where $x_j$ - is approximation to the eigenvector on the $j$-th iteration.<br> **Note** that sometimes approximation to $(A-R(x_j) I)^{-1}$ is not a good preconditioner.\n",
    "\n",
    "\n",
    "- It additionally adds to a subspace $V$ solutions from previous iterations (**subspace acceleration**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### JD derivation\n",
    "\n",
    "- Jacobi-Davidson method has a nice manifold optimization interpretation. \n",
    "- It is a **Riemannian Newton** method on a sphere and $P = I - x_j x^*_j$ is a projection on a tanget space of a sphere at $x_j$.\n",
    "\n",
    "But we will derive it similarly to the original paper."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Jacobi correction equation\n",
    "\n",
    "Jacobi not only presents the way to solve the eigenvalue problem by Jacobi rotations, but also proposed an iterative procedure. Let $x_j$ be the current approximation, and $t$ the correction:\n",
    "\n",
    "$$A(x_j + t) = \\lambda (x_j + t),$$\n",
    "\n",
    "and we look for the correction $t \\perp x_j$ (new orthogonal vector).\n",
    "\n",
    "Then, the parallel part has the form\n",
    "\n",
    "$$x_j x^*_j A (x_j + t) = \\lambda x_j,$$\n",
    "\n",
    "which simplifies to \n",
    "\n",
    "$$R(x_j) + x^* _j A t = \\lambda.$$\n",
    "\n",
    "The orthogonal component is \n",
    "\n",
    "$$( I - x_j x^*_j) A (x_j + t) = (I - x_j x^*_j) \\lambda (x_j + t),$$\n",
    "\n",
    "which is equivalent to \n",
    "\n",
    "$$\n",
    "  (I - x_j x^*_j) (A - \\lambda I) t = (I - x_j x^*_j) (- A x_j + \\lambda x_j) = - (I - x_j x^*_j) A x_j = - (A - R(x_j) I) x_j = -r_j.\n",
    "$$\n",
    "\n",
    "$r_j$ is the residual.\n",
    "\n",
    "Since $(I - x_j x^*_j) t  = t$, we can rewrite this equation in the symmetric form"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "$$ (I - x_j x^*_j) (A - \\lambda I) (I - x_j x^*_j) t = -r_j.$$\n",
    "\n",
    "Now we replace $\\lambda$ by $R(x_j)$, and get the **Jacobi correction equation**:\n",
    "\n",
    "$$\n",
    " (I - x_j x^*_j) (A - R(x_j) I) (I - x_j x^*_j) t = -r_j.\n",
    "$$\n",
    "\n",
    "Since $r_j \\perp x_j$ this equation is consistent, if $(A - R(x_j) I)$ is non-singular."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Solving Jacobi correction equation\n",
    "\n",
    "Typically Jacobi equation is solved inexactly by the appropriate Krylov method.\n",
    "\n",
    "Even inexact solution of Jacobi equation ensures (why?) that the correction $t$ is orthogonal to $x_j$, which is good for computations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Connection to the Rayleigh quotient iteration\n",
    "\n",
    "If this equation is solved exactly, we will get Rayleigh quotient iteration! Let us show that.\n",
    "\n",
    "$$ (I - x_j x^*_j) (A - R(x_j) I) (I - x_j x^*_j) t = -r_j.$$\n",
    "\n",
    "$$ (I - x_j x^*_j) (A - R(x_j) I) t = -r_j.$$\n",
    "\n",
    "$$  (A - R(x_j) I) t - \\alpha x_j = -r_j, \\quad \\alpha = x^*_j (A - R(x_j) I) x_j$$\n",
    "\n",
    "$$   t = \\alpha (A - R(x_j) I)^{-1}x_j  - (A - R(x_j) I)^{-1}r_j,$$\n",
    "\n",
    "Thus, since $(A - R(x_j) I)^{-1}r_j = (A - R(x_j) I)^{-1}(A - R(x_j) I)x_j = x_j$ we get\n",
    "\n",
    "$$x_{j+1} = x_j + t = \\alpha (A - R(x_j) I)^{-1}x_j,$$\n",
    "\n",
    "which is Rayleigh quotient iteration up to normalization."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Preconditioning of Jacobi equation\n",
    "\n",
    "A popular preconditioner for solving Jacobi equation by Krylov method has the form\n",
    "\n",
    "$$\n",
    "\\widetilde K = (I - x_j x^*_j) K (I - x_j x^*_j)\n",
    "$$\n",
    "\n",
    "where $K$ is easily-inverted approximation of $(A - R(x_j) I)$.\n",
    "\n",
    "We need to derive how to solve a system with $\\widetilde K$ in terms of solving a system with $K$.\n",
    "\n",
    "We already showed that equation\n",
    "\n",
    "$$ (I - x_j x^*_j) K (I - x_j x^*_j) \\tilde t = f $$\n",
    "\n",
    "is equavelnt to \n",
    "\n",
    "$$  \\tilde t = \\alpha K^{-1}x_j  + K^{-1}f $$\n",
    "\n",
    "The trick now is to forget about the value of $\\alpha$ and find it from $\\tilde t\\perp x_j$ to maintain orthogonality:\n",
    "\n",
    "$$\n",
    "    \\alpha = \\frac{x_j^*K^{-1}f}{x_j^* K^{-1}x_j}\n",
    "$$\n",
    "Thus for each iteration of the Jacobi equation we calculate $K^{-1}x_j$ and then update only $K^{-1}f$ on each internal Krylov iteration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Subspace acceleration in JD\n",
    "\n",
    "On each iteration of the method we expand a basis with new $t$.\n",
    "\n",
    "Namely, $V_j = [v_1,\\dots,v_{j-1},v_j]$, where $v_j$ is vector $t$ orthogonalized to $V_{j-1}$.\n",
    "\n",
    "Then standard Rayleigh-Ritz procedure is used.\n",
    "\n",
    "**Historal fact:** Initially subspace acceleration was used in the **Davidson method**. <br> \n",
    "\n",
    "However, instead of the Jacobi equation, equation $(\\mathrm{diag}(A) - R(x_j)I)t = -r_j$ was used. <br>\n",
    "Davidson method was very popular in quantum chemistry computations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### The block case of JD\n",
    "\n",
    "If we want many eigenvectors, we just compute **partial Schur decomposition:**\n",
    "\n",
    "$$A Q_k = Q_k T_k, $$\n",
    "\n",
    "and then want to update $Q_k$ by one vector added to $Q_k$. We just use instead of $A$ the matrix $(I - Q_k Q^*_k) A (I - Q_k Q^*_k)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Jacobi-Davidson: summary\n",
    "\n",
    "The correction equation can be solved only roughly, and JD method is often the fastest."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Software\n",
    "\n",
    "- The [ARPack](http://www.caam.rice.edu/software/ARPACK/) is the most widely used (it powers scipy sparse eigensolver). Includes versions of Lanczos and Arnoldi algorithms.\n",
    "- The [PRIMME](https://github.com/primme/primme) is the best from my experience (it employs dynamic switching between different methods including LOBPCG and JD)\n",
    "- [PROPACK](http://sun.stanford.edu/~rmunk/PROPACK/) works well for the SVD."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Comparison of Lanczos, LOBPCG and PRIMME "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of nnz 14998\n",
      "[3.99999961 3.99999842 3.99999645 3.99999369 3.99999013 3.99998579\n",
      " 3.99998066 3.99997474 3.99996804 3.99996054]\n",
      "[3.99999961 3.99999842 3.99999645 3.99999369 3.99999013 3.99998579\n",
      " 3.99998066 3.99997474 3.99996804 3.99996054]\n",
      "[3.99999821 3.9999956  3.99999403 3.99999153 3.99998655 3.99998388\n",
      " 3.99997659 3.99996981 3.99996498 3.99995576]\n",
      "16.7 s ± 667 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "2.66 s ± 74.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "1.57 s ± 93.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQgAAAD8CAYAAACLgjpEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQf0lEQVR4nO3cf2jc933H8ecrTuJmTb3YqeJ5koO9oIw63ZbGh9HIWEOoiR2XOn960EWQgiFkkLBBsSkM+sfA2x/FhC0Bk4YptKsJtBCTzRnGTdZ/vLhykzRxXNfy7GERJ1JW2roZeIv73h/3WfONdCedpLvvz9cDjvve5+57el9Rnv7e905VRGBm1sl1RQ9gZuXlQJhZVw6EmXXlQJhZVw6EmXXlQJhZV6UNhKQdks5KmpK0r8A5npU0I+mtzNo6SccknUvXazP37U8zn5X0QGZ9q6Q3031PStIAZt0o6WVJZySdlvR4WeeV9AlJJyW9kWb9ellnzfycVZJek/Ri2Wftm4go3QVYBZwHfg+4EXgD2FLQLH8K3AO8lVn7O2Bf2t4H/G3a3pJmXQ1sTq9hVbrvJPDHgICjwM4BzLoBuCdtfwr4aZqpdPOm5705bd8AvAqMlXHWzMx/CfwT8GKZfw/6eSnrEcQ2YCoi/iMi/gc4DOwuYpCI+AHwsznLu4GJtD0BPJRZPxwRVyPiAjAFbJO0AVgTESei/VvyXGaffs56OSJ+lLavAGeA4TLOG22/SjdvSJco46wAkkaAXcAzmeVSztpPZQ3EMHApc3s6rZXF+oi4DO3/KIHb0nq3uYfT9tz1gZG0Cfgc7X+ZSzlvOmR/HZgBjkVEaWcFDgJfBX6dWSvrrH1T1kB0el9Whe+Ed5s719cj6Wbgu8ATEfHLhR7aYS23eSPiWkTcDYzQ/hf2sws8vLBZJX0RmImIU73u0mEt99+DfihrIKaBjZnbI8A7Bc3SyXvpcJF0PZPWu809nbbnrvedpBtox+HbEfG9ss8LEBE/B14BdpR01nuBL0m6SPvt7v2SvlXSWfuqrIH4ITAqabOkG4E9wJGCZ8o6Aoyn7XHghcz6HkmrJW0GRoGT6fDziqSxdNb64cw+fZOe+5vAmYj4RpnnlTQk6Za0fRPwBeAnZZw1IvZHxEhEbKL9u/j9iPhyGWftu6LPkna7AA/SPgt/HvhagXN8B7gM/C/tfwG+AtwKHAfOpet1mcd/Lc18lswZaqAFvJXu+3tAA5j1T2gfsv4YeD1dHizjvMAfAq+lWd8C/jqtl27WOXPfx0efYpR61n5clIY2M5unrG8xzKwEHAgz68qBMLOuHAgz6yr3QKgkf4RlZovLNRCSVgH/AOyk/QctfyZpyyL77M1jtn6o0qxQrXmrNCtUb95u8j6CWM4fYVXpf+gqzQrVmrdKs0L15u0o70CU/Y+wzCzj+px/Xk9/rJIOz/YCXHfTmq2rN4wGwB8M//Zgp1uh22+/nVarVZlvnlVp3irNCtWa99SpU+9HxFCn+/IORE9/hBURh4BDAKs3jMaG8YMAvA9cPLBr4EOaNYmk/+x2X95vMVb8R1ib9v3zQAYzs/lyDUREfAj8BfCvtP/fjp6PiNML7fOp1fMPchwJs3zk/j2IiPiXiLgzIu6IiL9Z7PGbPv1JPn/nrfPXHQmzgavENyknHhlzJMwKUIlAgCNhVoTKBAIcCbO8VSoQ4EiY5alygQBHwiwvlQwEOBJmeahsIMCRMBu0SgcCHAmzQap8IMCRMBuUWgQCHAmzQahNIMCRMOu3WgUCHAmzfqpdIMCRMOuXWgYCHAmzfqhtIMCRMFupWgcCHAmzlah9IMCRMFuuRgQCHAmz5WhMIMCRMFuqRgUCHAmzpWhcIMCRMOtVIwMBjoRZLxobCHAkzBbT6ECAI2G2kMYHAhwJs24ciMSRMJvPgchwJMw+zoGYw5Ew+4gD0YEjYdbmQHThSJg5EAtyJKzpHIhFOBLWZA5EDxwJayoHokeOhDWRA7EEjoQ1zaKBkPSspBlJb2XW1kk6Julcul6buW+/pClJZyU9kFnfKunNdN+TktT/lzN4joQ1SS9HEP8I7Jiztg84HhGjwPF0G0lbgD3AXWmfpyStSvs8DewFRtNl7nNWhiNhTbFoICLiB8DP5izvBibS9gTwUGb9cERcjYgLwBSwTdIGYE1EnIiIAJ7L7FNJjoQ1wXLPQayPiMsA6fq2tD4MXMo8bjqtDaftuesdSdoraVLS5Ozs7DJHHDxHwuqu3ycpO51XiAXWO4qIQxHRiojW0NBQ34YbBEfC6my5gXgvvW0gXc+k9WlgY+ZxI8A7aX2kw3otTDwyxmd+55Pz1h0Jq7rlBuIIMJ62x4EXMut7JK2WtJn2yciT6W3IFUlj6dOLhzP71MLRJ+5zJKx2evmY8zvACeD3JU1L+gpwANgu6RywPd0mIk4DzwNvAy8Bj0XEtfRUjwLP0D5xeR442ufXUjhHwupG7Q8VyqvVasXk5GTRYyzJzoOvcObdD+atXzywq4BpzBYm6VREtDrd529SDoCPJKwuHIgBcSSsDhyIAXIkrOociAFzJKzKHIgcOBJWVQ5EThwJqyIHIkeOhFWNA5EzR8KqxIEogCNhVeFAFMSRsCpwIArkSFjZORAFcySszByIEnAkrKwciJJwJKyMHIgScSSsbByIknEkrEwciBJyJKwsHIiSciSsDByIEnMkrGgORMk5ElYkB6ICHAkrigNREY6EFcGBqBBHwvLmQFSMI2F5ciAqyJGwvDgQFeVIWB4ciApzJGzQHIiKcyRskByIGnAkbFAciJpwJGwQHIgacSSs3xyImnEkrJ8ciBpyJKxfHIiaciSsHxYNhKSNkl6WdEbSaUmPp/V1ko5JOpeu12b22S9pStJZSQ9k1rdKejPd96QkDeZlGTgStnK9HEF8CPxVRHwGGAMek7QF2Accj4hR4Hi6TbpvD3AXsAN4StKq9FxPA3uB0XTZ0cfXYh04ErYSiwYiIi5HxI/S9hXgDDAM7AYm0sMmgIfS9m7gcERcjYgLwBSwTdIGYE1EnIiIAJ7L7GMD5EjYci3pHISkTcDngFeB9RFxGdoRAW5LDxsGLmV2m05rw2l77rrlwJGw5eg5EJJuBr4LPBERv1zooR3WYoH1Tj9rr6RJSZOzs7O9jmiLcCRsqXoKhKQbaMfh2xHxvbT8XnrbQLqeSevTwMbM7iPAO2l9pMP6PBFxKCJaEdEaGhrq9bVYDxwJW4pePsUQ8E3gTER8I3PXEWA8bY8DL2TW90haLWkz7ZORJ9PbkCuSxtJzPpzZx3LkSFivejmCuBf4c+B+Sa+ny4PAAWC7pHPA9nSbiDgNPA+8DbwEPBYR19JzPQo8Q/vE5XngaD9fjPXOkbBeqP2BQnm1Wq2YnJwseoza2nnwFc68+8G89YsHdhUwjRVB0qmIaHW6z9+kbDgfSdhCHAhzJKwrB8IAR8I6cyDsNxwJm8uBsI9xJCzLgbB5HAn7fw6EdeRIGDgQtgBHwhwIW5Aj0WwOhC3KkWguB8J64kg0kwNhPXMkmseBsCVxJJrFgbAlcySaw4GwZXEkmsGBsGVzJOrPgbAVcSTqzYGwFXMk6suBsL5wJOrJgbC+cSTqx4GwvnIk6sWBsL5zJOrDgbCBcCTqwYGwgXEkqs+BsIFyJKrNgbCBcySqy4GwXDgS1eRAWG4ciepxICxXjkS1OBCWO0eiOhwIK4QjUQ0OhBXGkSg/B8IK5UiUmwNhhXMkysuBsFJwJMpp0UBI+oSkk5LekHRa0tfT+jpJxySdS9drM/vslzQl6aykBzLrWyW9me57UpIG87KsihyJ8unlCOIqcH9E/BFwN7BD0hiwDzgeEaPA8XQbSVuAPcBdwA7gKUmr0nM9DewFRtNlR/9eitWBI1EuiwYi2n6Vbt6QLgHsBibS+gTwUNreDRyOiKsRcQGYArZJ2gCsiYgTERHAc5l9zH7DkSiPns5BSFol6XVgBjgWEa8C6yPiMkC6vi09fBi4lNl9Oq0Np+25651+3l5Jk5ImZ2dnl/ByrC4ciXLoKRARcS0i7gZGaB8NfHaBh3c6rxALrHf6eYciohURraGhoV5GtBpyJIq3pE8xIuLnwCu0zx28l942kK5n0sOmgY2Z3UaAd9L6SId1s64ciWL18inGkKRb0vZNwBeAnwBHgPH0sHHghbR9BNgjabWkzbRPRp5Mb0OuSBpLn148nNnHrCtHoji9HEFsAF6W9GPgh7TPQbwIHAC2SzoHbE+3iYjTwPPA28BLwGMRcS0916PAM7RPXJ4HjvbxtViNORLFUPsDhfJqtVoxOTlZ9BhWEjsPvsKZdz+Yt37xwK4CpqkHSaciotXpPn+T0irFRxL5ciCschyJ/DgQVkmORD4cCKssR2LwHAirNEdisBwIqzxHYnAcCKsFR2IwHAirDUei/xwIqxVHor8cCKsdR6J/HAirJUeiPxwIqy1HYuUcCKs1R2JlHAirPUdi+RwIawRHYnkcCGsMR2LpHAhrFEdiaRwIaxxHoncOhDWSI9EbB8Iay5FYnANhjeZILMyBsMZzJLpzIMxwJLpxIMwSR2I+B8Isw5H4OAfCbA5H4iMOhFkHjkSbA2HWhSPhQJgtqOmRcCDMFtHkSDgQZj1oaiQcCLMeNTESDoTZEjQtEg6E2RI1KRI9B0LSKkmvSXox3V4n6Zikc+l6beax+yVNSTor6YHM+lZJb6b7npSk/r4cs3w0JRJLOYJ4HDiTub0POB4Ro8DxdBtJW4A9wF3ADuApSavSPk8De4HRdNmxounNCtSESPQUCEkjwC7gmczybmAibU8AD2XWD0fE1Yi4AEwB2yRtANZExImICOC5zD5mlVT3SPR6BHEQ+Crw68za+oi4DJCub0vrw8ClzOOm09pw2p67blZpdY7EooGQ9EVgJiJO9ficnc4rxALrnX7mXkmTkiZnZ2d7/LFmxalrJHo5grgX+JKki8Bh4H5J3wLeS28bSNcz6fHTwMbM/iPAO2l9pMP6PBFxKCJaEdEaGhpawssxK04dI7FoICJif0SMRMQm2icfvx8RXwaOAOPpYePAC2n7CLBH0mpJm2mfjDyZ3oZckTSWPr14OLOPWS0cfeI+Pn/nrfPWqxqJlXwP4gCwXdI5YHu6TUScBp4H3gZeAh6LiGtpn0dpn+icAs4DR1fw881KaeKRsdpEQu0PFMqr1WrF5ORk0WOYLdn4s//Ov/30v+atXzywq4BpupN0KiJane7zNynNBqQORxIOhNkAVT0SDoTZgFU5Eg6EWQ6qGgkHwiwnVYyEA2GWo6pFwoEwy1mVIuFAmBWgKpFwIMwKUoVIOBBmBSp7JBwIs4KVORIOhFkJlDUSDoRZSZQxEg6EWYmULRIOhFnJlCkSDoRZCZUlEg6EWUmVIRIOhFmJFR0JB8Ks5IqMhANhVgFFRcKBMKuIIiLhQJhVSN6RcCDMKibPSDgQZhWUVyQcCLOKyiMSDoRZhQ06Eg6EWcUNMhIOhFkNDCoSDoRZTQwiEg6EWY30OxIOhFnN9DMSDoRZDfUrEg6EWU31IxIOhFmNrTQSDoRZza0kEg6EWQMsNxIOhFlDdIvEjevvuKfbPg6EWYN0jISkbo9XRAx6phWRdAU4W/QcPfo08H7RQyxBleat0qxQ8nmvX/u7d1y3+rduAfjwFzNc++9fdIzE9blOtTxnI6JV9BC9kDRZlVmhWvNWaVao3rzd+C2GmXXlQJhZV1UIxKGiB1iCKs0K1Zq3SrNC9ebtqPQnKc2sOFU4gjCzgjgQZtaVA2FmXTkQZtaVA2FmXf0fVMQweXhgfRMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# !pip install primme\n",
    "import primme\n",
    "import scipy.sparse.linalg as spsplin\n",
    "import scipy.sparse as spsp\n",
    "\n",
    "# A = spsp.rand(10000, 10000, density=1e-5, random_state=10, format=\"csr\")  * 10\n",
    "# A = A + A.T\n",
    "\n",
    "n = 5000\n",
    "ex = -np.ones(n)\n",
    "A = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr')\n",
    "# e = sp.sparse.eye(n)\n",
    "# A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)\n",
    "\n",
    "print(\"Number of nnz\", A.nnz)\n",
    "plt.spy(A, markersize=1)\n",
    "k = 10\n",
    "print(spsplin.eigsh(A, k=k)[0][::-1])\n",
    "print(primme.eigsh(A, k=k)[0])\n",
    "print(spsplin.lobpcg(A, np.random.randn(A.shape[0], k), maxiter=700)[0])\n",
    "\n",
    "%timeit spsplin.eigsh(A, k=k)\n",
    "%timeit primme.eigsh(A, k=k)\n",
    "%timeit spsplin.lobpcg(A, np.random.randn(A.shape[0], k), maxiter=700)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Take-home message about large scale partial eigenvalue problem\n",
    "\n",
    "- Arnoldi and Lanczos methods. Shift-and-invert strategy is very expensive since inversion must be done very accurately.\n",
    "- Preconditioned iterative methods (PINVIT, LOBPCG, JD). Good for inexact inversions. \n",
    "- There is a software for using them\n",
    "- There is a lot of technical issues hidden (restarts, spurious eigenvalues, stability)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Intro to streaming algorithms\n",
    "\n",
    "- Streaming algorithms process data samples one-by-one and do not store entire dataset\n",
    "- Sketching approach\n",
    "- Typical applications are solving linear systems, PCA, graph sparsification "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Oja's method\n",
    "\n",
    "- Embed the stream data in the orthonormal low dimensional subspace linearly such that the variance is maximized\n",
    "- Optimization problem\n",
    "\n",
    "$$ \\max_{W \\in \\mathbb{R}^{n \\times k}, W^{\\top}W = I} \\mathrm{trace}(W^{\\top}\\hat{\\Sigma}W), $$\n",
    "\n",
    "where $\\hat{\\Sigma}$ is sampled covariance matrix.\n",
    "\n",
    "- The solution is the $k$ leading eigenvector\n",
    "- But the samples are coming one-by-one, so you can not store them all in memory\n",
    "- Stochastic gradient method helps in this case\n",
    "- The result Oja's method\n",
    "    - Update: $W_{k+1} = W_k + \\eta_k x_t x_t^{\\top}W_k$, where $\\eta_k$ is a learning rate\n",
    "    - Project: compute $Q$ from QR decompisition of $W_{k+1}$ and replace $W_{k+1} = Q$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Sketching\n",
    "\n",
    "- Replace matrix $A$ with another matrix $SA$, where matrix $SA$ has significantly smaller number of rows but preserves some important properties of matrix $A$\n",
    "- Possible choices:\n",
    "    - random projection\n",
    "    - random row selection\n",
    "- Example: linear least squares problem $\\|Ax - b\\|_2^2 \\to \\min_x$ transforms to $\\| (SA)y - Sb \\|_2^2 \\to \\min_y$ and we expect that $x \\approx y$\n",
    "- [Blendenpick](https://pdos.csail.mit.edu/~petar/papers/blendenpik-v1.pdf) solver is based on that idea and outperforms LAPACK routine \n",
    "- More details see in [Sketching as a Tool for Numerical Linear Algebra](https://researcher.watson.ibm.com/researcher/files/us-dpwoodru/wNow3.pdf) by D. Woodruff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Next lecture\n",
    "\n",
    "\n",
    "- Final lecture: we will talk about multilinear algebra and tensors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "metadata": {
  "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.7.7"
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": true,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}