lectures/lecture12/lecture12.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Лекция 12\n",
    "\n",
    "## Итерационные методы для поиска собственных значений"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## На прошлой лекции\n",
    "\n",
    "- Методы MINRES, BiCG и BiCGStab \n",
    "- Предобуславливатели: метод Якоби, Гаусса-Зейделя и SOR($\\omega$)\n",
    "- Неполное LU разложение и его варианты"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Частичная задача на собственные значения\n",
    "\n",
    "- Напомним, что для поиска собственных значений матрицы $N\\times N$ можно использовать например QR алгоритм.\n",
    "\n",
    "- Однако в некоторых приложениях матрицы настолько большие, что мы даже не можем их хранить в памяти \n",
    "\n",
    "- Обычно такие матрицы даны в виде **чёрного ящика**, который позволяет только умножить эту матрицу на вектор. Далее рассмотрим именно этот случай.\n",
    "\n",
    "- Лучшее, что мы можем сделать это решить частичную задачу на собственные значения, то есть\n",
    "\n",
    "    - Найти $k\\ll N$ наименьших или наибольших собственных значений (и собственных векторов, если необходимо)\n",
    "    - Найти $k\\ll N$ собственных значений, ближайших к заданному числу $\\sigma$\n",
    "\n",
    "- Для простоты рассмотрим случай нормальной матрицы, которая имеет ортонормированный базис из собственных векторов. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Степенной метод и его аналоги\n",
    "\n",
    "#### Степенной метод: напоминание\n",
    "\n",
    "- Простейший метод для поиска максимального по модулю собственного значения – это **степенной метод**\n",
    "\n",
    "$$ x_{i+1} = \\frac{Ax_{i}}{\\|Ax_{i}\\|_2}. $$\n",
    "\n",
    "- Сходимость линейная с коэффициентом $q = \\left|\\frac{\\lambda_1}{\\lambda_2}\\right|$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Метод обратной итерации: напоминание\n",
    "\n",
    "- Для поиска наименьшего собственного значения можно запустить степенной метод для матрицы $A^{-1}$:\n",
    "\n",
    "$$x_{i+1} = \\frac{A^{-1}x_{i}}{\\|A^{-1}x_{i}\\|}.$$\n",
    "\n",
    "- Для ускорения сходимости может быть использована стратегия <font color='red'>shift-and-invert</font>:\n",
    "\n",
    "$$x_{i+1} = \\frac{(A-\\sigma I)^{-1}x_{i}}{\\|(A-\\sigma I)^{-1}x_{i}\\|},$$\n",
    "\n",
    "где $\\sigma$ должна лежать близко от целевого собственного значения."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Метод Релея\n",
    "\n",
    "- Для получения сверхлинейной сходимости можно воспользоваться адаптивными сдвигами:\n",
    "\n",
    "$$x_{i+1} = \\frac{(A-R(x_i) I)^{-1}x_{i}}{\\|(A-R(x_i) I)^{-1}x_{i}\\|},$$\n",
    "\n",
    "где $R(x_k) = \\frac{(x_i, Ax_i)}{(x_i, x_i)}$ соотношение Релея. \n",
    "\n",
    "- Метод сходится **кубически для эрмитовых матриц** и квадратично в противном случае."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Неточный метод обратной итерации\n",
    "\n",
    "- Матрицы $(A- \\sigma I)$ также как $(A-R(x_i) I)$ плохо обусловлены, если $\\sigma$ или $R(x_i)$ близки к собственным значениям.\n",
    "\n",
    "- Поэтому если у вас нет LU разложения этой матрицы, то могут возникнуть проблемы при решении систем на каждой итерации\n",
    "\n",
    "- На практике вы можете решать системы **только с некоторой точностью**. Напомним также, что число обусловленности даёт оценку сверху, которая может быть завышенной для подходящих правых частей. Поэтому даже в методе Релея близость сдвига к собственному значению существенно [не ухудшает](http://www.sciencedirect.com/science/article/pii/S0024379505005756) сходимость итерационного метода.\n",
    "\n",
    "- Если точность решения системы возрастает от итерации к итерации, сверхлинейная сходимость для метода Релея также присутствует, см [Theorem 2.1](http://www.sciencedirect.com/science/article/pii/S0024379505005756). Иначе вы получите только линейную сходимость."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Перед тем как мы перейдём к продвинутым методам, обсудим важную концепцию **аппроксимации Ритца**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Аппроксимация Ритца\n",
    "\n",
    "- Для данного подпространства, натянутого на столбцы унитарной матрицы $Q_k$ размера $N\\times k$, рассмотрим спроецированную матрицу $Q_k^* A Q_k$.\n",
    "\n",
    "- Пусть $\\Theta_k=\\mathrm{diag}(\\theta_1,\\dots,\\theta_k)$ и $S_k=\\begin{bmatrix}s_1 & \\dots & s_k \\end{bmatrix}$ матрицы собственных значений и собственных векторов для матрицы $Q_k^* A Q_k$: \n",
    "\n",
    "$$(Q_k^* A Q_k)S_k = S_k \\Theta_k$$\n",
    "\n",
    "тогда $\\{\\theta_i\\}$ называются **числами Ритца** и $y_i = Q_k s_i$ - **векторы Ритца**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Свойства аппроксимации Ритца\n",
    "\n",
    "- Заметим, что числа и векторы Ритца не являются собственными значениями и собственными векторами исходной матрицы $AY_k\\not= Y_k \\Theta_k$, но выполнено следующее равенство:\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",
    "   таким образом невязка для аппроксимации Ритца **ортогональна** подпространству, натянутому на столбцы $Q_k$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- $\\lambda_\\min(A) \\leq \\theta_\\min \\leq \\theta_\\max \\leq \\lambda_\\max(A)$. Действительно, используя отношение Релея:\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",
    "- Очевидно, что $\\lambda_\\min (Q_k^* A Q_k) = \\lambda_\\min(A)$, если $k=N$, но мы хотим построить базис размера $k\\ll N$ такой что $\\lambda_\\min (Q_k^* A Q_k) \\approx \\lambda_\\min(A)$.\n",
    "\n",
    "- Таким же образом можно показать, что $\\theta_\\max \\leq \\lambda_\\max(A)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### <font color='red'>Метод Релея-Ритца</font>\n",
    "\n",
    "Таким образом, если подпространство $V$ приближает первые $k$ собственных векторов, тогда вы можете использовать **метод Релея-Ритца**:\n",
    "\n",
    "1. Найти ортонормированный базис $Q_k$ в $V$ (например с помощью QR разложения)\n",
    "2. Вычислить $Q_k^*AQ_k$\n",
    "3. Вычислить вектора и числа Ритца\n",
    "4. Заметим, что также можно использовать $V$ без ортогонализации, но в этом случае нужно будет решать обобщённую задачу на собственные векторы $(V^*AV)s_i = \\theta_i (V^*V)s_i$.\n",
    "\n",
    "Вопрос в том, как найти хорошее подпространство $V$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Какое подпространство мы будем использовать?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Метод Ланцоша и Арнольди\n",
    "\n",
    "- Хорошим выбором $V$ будет **Крыловское подпространство**.\n",
    "\n",
    "- Напомним, что в степенном методе мы использовали только один Крыловский вектор\n",
    "\n",
    "$$x_k = \\frac{A^k x_0}{\\|A^k x_0\\|}.$$\n",
    "\n",
    "- В этом случае $\\theta_k = \\frac{x_k^* A x_k}{x_k^* x_k}$ есть ни что иное как число Ритца. Естественная идея – использовать Крыловское пространство большей размерности.\n",
    "\n",
    "- В результате мы найдём больше собственных значений, а также сходимость к собственному вектору для $\\lambda_\\max$ будет быстрее, чем в степенном методе."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- Для эрмитовой матрицы из соотношения Арнольди следует, что \n",
    "\n",
    "$$ Q_k^*AQ_k = T_k, $$\n",
    "\n",
    "где $Q_k$ ортогональный базис в Крыловском подпространстве, сгенерированный процессом Ланцоша и $T_k$ трёхдиагональная матрица.\n",
    "\n",
    "- В соответствии с методом Релея-Ритца мы ожидаем, что собственные значения $T_k$ приближают собственные значения матрицы $A$. \n",
    "- Этот метод называется **методом Ланцоша**. \n",
    "- Для несимметричных матриц он называется **методом Арнольди** и вместо трёхдиагональной матрицы мы получим верхне-гессенбергову матрицу.\n",
    "\n",
    "Давайте теперь покажем, что  $\\theta_\\max \\approx\\lambda_\\max$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Почему $\\theta_\\max \\approx \\lambda_\\max$?\n",
    "\n",
    "- Обозначим $\\theta_1 \\equiv \\theta_\\max$ и $\\lambda_1 \\equiv \\lambda_\\max$. Тогда\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",
    "где $p_{i-1}$ полином степени не выше $i-1$ такой что $p_{i-1}(A)x_0\\not=0$.\n",
    "\n",
    "- Разложим $x_0 = \\sum_{j=1}^N c_j v_j$, где $v_j$ собственные векторы $A$ (которые образуют ортонормированный базис).\n",
    "\n",
    "- Поскольку $\\theta_1 \\leq \\lambda_1$ получим\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",
    "для любого полинома $p_{i-1}$. Таким образом\n",
    "\n",
    "\\begin{align*}\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} = \\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} \\\\\n",
    "& \\leq (\\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",
    "\\end{align*}\n",
    "\n",
    "- Так как неравенство выполнено для любого полинома $p_{i-1}$, мы выберем полином: \n",
    "\n",
    "$$|p_{i-1}(\\lambda_1)| \\gg \\max_{2\\leq k \\leq N}|p_{i-1}(\\lambda_k)|.$$\n",
    "\n",
    "- Это неравенство выполнено например для полиномов Чебышёва на $[\\lambda_n,\\lambda_2]$. \n",
    "- В итоге $\\theta_1 \\approx \\lambda_1$ или более точно:\n",
    "\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",
    "\n",
    "где $T_{i-1}$ – полином Чебышёва."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Демо: аппроксимация максимального собственного значения с помощью метода Ланцоша"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k=10, err = 0.14693010320059852\n",
      "k=20, err = 0.028873846136214354\n",
      "k=100, err = 1.503028812521734e-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",
    "import numpy as np\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": [
    "### Устойчивость\n",
    "\n",
    "- Векторы Ланцоша могут терять ортогональность в процессе вычисления из-за арифметики с плавающей точкой, поэтому все практически полезные реализации используют **рестарты**.\n",
    "\n",
    "- Очень хорошее введение в тему содержится в книге **Matrix Computations** авторов G. Golub и Ch. Van-Loan."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Другие недостатки метода Ланцоша\n",
    "\n",
    "- Применение метода Ланцоша напрямую к матрице $A$ может привести к очень медленной сходимости, если  $\\lambda_i\\approx \\lambda_{i+1}$ (обычно это происходит для наименьших собственных значений, которые плохо разделены)\n",
    "\n",
    "- Для ускорения сходимости можно применить метод Ланцоша к матрице $(A-\\sigma I)^{-1}$, но в этом случае системы должны решаться **очень точно**. В противном случае соотношение Арнольди перестанет выполняться.\n",
    "\n",
    "Альтернативой этому подходу являются так называемые предобусловленные итерационные методы, например:\n",
    "1. PINVIT (Предобусловленный метод обратной итерации)\n",
    "2. LOBPCG (локально оптимальный блочный предобусловленный метод сопряжённых градиентов)\n",
    "3. Метод Якоби-Дэвидсона (Jacobi-Davidson method)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### PINVIT (предобусловленный метод обратной итерации)\n",
    "\n",
    "#### Получение метода\n",
    "\n",
    "- Рассмотрим отношение Релея $R(x) = \\frac{(x,Ax)}{(x,x)}$. Тогда\n",
    "\n",
    "$$ \\nabla R(x) = \\frac{2}{(x,x)} (Ax - R(x) x), $$\n",
    "\n",
    "и простейший метод градиентного спуска с предобуславливателем $B$ записывается в виде\n",
    "\n",
    "$$ x_{i+1} = x_{i} - \\tau_i B^{-1} (Ax_i - R(x_i) x_i), $$\n",
    "\n",
    "$$ x_{i+1} = \\frac{x_{i+1}}{\\|x_{i+1}\\|}. $$\n",
    "\n",
    "- Обычно $B\\approx A-\\sigma I$, где $\\sigma$ обозначает сдвиг.\n",
    "\n",
    "- Чем $\\sigma$ ближе к необходимому собственному значению, тем быстрее сходимость."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- Параметр $\\tau_k$ выбирается так, чтобы минимизировать $R(x_{i+1})$ по $\\tau_k$ (метод наискорейшего спуска).\n",
    "\n",
    "- Эта процедура минимизации может быть рассмотрена как минимизация в базисе $V = [x_i, r_i]$, где $r_{i}=B^{-1} (Ax_i - R(x_i) x_i)$.\n",
    "\n",
    "- Это приводит к обобщённой задаче на собственные значения $(V^*AV)\\begin{bmatrix}1 \\\\ -\\tau_i \\end{bmatrix} = \\theta (V^*V) \\begin{bmatrix}1 \\\\ -\\tau_i \\end{bmatrix}$ (процедура Релея-Ритца без ортогонализации $V$). Здесь $\\theta$ – ближайшее число к требуемому собственному значению."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Сходимость\n",
    "\n",
    "**Теорема** ([Knyazev и Neymeyr](http://www.sciencedirect.com/science/article/pii/S002437950100461X)) \n",
    "\n",
    "Пусть \n",
    "- $R(x_{i})\\in [\\lambda_j,\\lambda_{j+1}]$\n",
    "- $R(x_{i+1})\\in [R(x_{i}),\\lambda_{j+1}]$ (случай $R(x_{i+1})\\in [\\lambda_{j}, R(x_{i})]$ аналогичен)\n",
    "- $\\|I - B^{-1} A\\|_A \\leq \\gamma < 1$\n",
    "\n",
    "тогда\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": [
    "#### Блочный случай\n",
    "\n",
    "- Для поиска $k$ собственных значений можно делать один шаг метода PINVIT для каждого вектора:\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",
    "$$ x^{(j)}_{i+1} = \\frac{x^{(j)}_{i+1}}{\\|x^{(j)}_{i+1}\\|}. $$\n",
    "\n",
    "- После чего ортогонализовать их с помощью QR разложения. Однако лучше использовать процедуру Релея-Ритца:\n",
    "\n",
    "    - Пусть $X^{i}_k = [x^{(1)}_{i},\\dots, x^{(k)}_{i}]$ и $R^{i}_k = [B^{-1}r^{(1)}_{i},\\dots, B^{-1}r^{(k)}_{i}]$, где $r^{(j)}_{i} = Ax^{(j)}_i - R(x^{(j)}_i) x^{(j)}_i$\n",
    "    - $V = [X^{i}_k, R^{i}_k]$, используем процедуру Релея-Ритца для $V$, чтобы найти новый $X^{i+1}_k$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## LOBPCG (локально оптимальный блочный предобусловленный метод CG)\n",
    "\n",
    "### Локально оптимальный предобусловленный метод  СG (пока ещё не блочный)\n",
    "\n",
    "LOPCG метод\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",
    "$$ x_{i+1} = \\frac{x_{i+1}}{\\|x_{i+1}\\|} $$\n",
    "\n",
    "превосходит метод PINVIT, поскольку он добавляет в базис не только $x_i$ и $r_i$, но также $x_{i-1}$.\n",
    "\n",
    "Однако такая интерпретация ведёт к неустойчивому алгоритму, так как $x_{i}$ становится коллинеарен $x_{i-1}$ вместе со сходимостью метода."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### LOPCG (устойчивая версия)\n",
    "\n",
    "- A. Knyazev предложил эквивалентную устойчивую версию, которая вводит новые векторы $p_i$ (сопряжённые градиенты)\n",
    "\n",
    "$$ p_{i+1} = r_{i} + \\beta_i p_{i}, $$\n",
    "\n",
    "$$ x_{i+1} = x_{i} + \\alpha_i p_{i+1}. $$\n",
    "\n",
    "- Можно показать, что $\\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": [
    "Устойчивая версия объясняет название метода:\n",
    "\n",
    "- В стандартном методе CG мы бы минимизировали отношение Релея $R$ по направлению сопряжённого направления $p_{i+1}$: \n",
    "\n",
    "$$\\alpha_i = \\arg\\min_{\\alpha_i} R(x_i + \\alpha_i p_{i+1}).$$\n",
    "\n",
    "- В локально оптимальном CG мы минимизируем по отношению к двум параметрам:  \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",
    "и получаем локально более оптимальное решение. Поэтому метод называется **локально оптимальным**.\n",
    "\n",
    "- По аналогии с методом PINVIT коэффициенты $\\alpha_i,\\beta_i$ можно найти из процедуру Релея-Ритца."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Локально оптимальный <font color='red'> блочный </font> предобусловленный метод CG\n",
    "\n",
    "- В блочной версии по аналогии с методом PINVIT на каждой итерации дан базис $V=[X^{(i)}_k,B^{-1}R^{(i)}_k, P^{(i)}_k]$ и используется процедура Релея-Ритца.\n",
    "\n",
    "Общая схема алгоритма:\n",
    "\n",
    "1. Найти $\\tilde A = V^* A V$\n",
    "2. Найти $\\tilde M = V^*V$\n",
    "3. Решить обобщённую задачу на собственные значения $\\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}$ (аналогично для $X^{(i+1)}_{k} = VS_k$)\n",
    "6. Вычислить новое значение $B^{-1}R^{(i+1)}_k$\n",
    "7. Задать $V=[X^{(i+1)}_k,B^{-1}R^{(i+1)}_k, P^{(i+1)}_k]$, goto 1.\n",
    "\n",
    "\n",
    "- Метод также сходится линейно, но быстрее чем PINVIT."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### LOBPCG: резюме\n",
    "\n",
    "- Локально оптимальный предобусловленный метод\n",
    "\n",
    "- Линейная сходимость\n",
    "\n",
    "- Предобуславливатель $A-\\sigma I$ не всегда хорош для задачи на собственные значения\n",
    "\n",
    "Следующий метод (Якоби-Дэвидсона) оснащён более \"умным\" предобуславливаталем и даёт сверхлинейную сходимость (если системы решаются точно)!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Метод Якоби-Дэвидсона (JD)\n",
    "\n",
    "Метод Якоби-Дэвидсона очень популярен для решения задачи на собственные значения (не только симметричной!).\n",
    "\n",
    "Он состоит из двух **основных ингредиентов**:\n",
    "\n",
    "- Для данного предобуславливателя к $A-R(x_j) I$ он автоматически строит хороший предобуславливатель для задачи на собственные значения:\n",
    "\n",
    "$$ B = (I - x_j x^*_j) (A - R(x_j) I) (I - x_j x^*_j), $$\n",
    "\n",
    "где $x_j$ аппроксимация собственного вектора на $j$-ой итерации\n",
    "\n",
    "- Заметим, что приближение $(A-R(x_j) I)^{-1}$ иногда не является хорошим предобуславливателем.\n",
    "\n",
    "- Дополнительно он прибавляет к подпространству $V$ решения с предыдущих итераций (**ускорение на подпространстве**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Вывод метода JD\n",
    "\n",
    "- Метод Якоби-Дэвидсона имеет красивую интерпретацию через оптимизацию на многообразии. \n",
    "- Он является **римановым методом Ньютона** на сфере и $P = I - x_j x^*_j$ проекция на касательное пространство к сфере в точке $x_j$.\n",
    "\n",
    "Но мы получим этот метод по аналогии с оригинальной работой."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Уравнение коррекции Якоби\n",
    "\n",
    "- Якоби предложил не только метод решения задачи на собственные значения через вращения, но также итерационный метод для решения этой задачи. \n",
    "- Пусть $x_j$ текущая аппроксимация собственного вектора, а $t$ коррекция к нему:\n",
    "\n",
    "$$A(x_j + t) = \\lambda (x_j + t),$$\n",
    "\n",
    "и мы ищем коррекцию $t \\perp x_j$ (новый ортогональный вектор).\n",
    "\n",
    "- Тогда его параллельная компонента имеет вид\n",
    "\n",
    "$$x_j x^*_j A (x_j + t) = \\lambda x_j,$$\n",
    "\n",
    "что приводит к выражению\n",
    "\n",
    "$$R(x_j) + x^* _j A t = \\lambda.$$\n",
    "\n",
    "- Ортогональная часть  \n",
    "\n",
    "$$( I - x_j x^*_j) A (x_j + t) = (I - x_j x^*_j) \\lambda (x_j + t),$$\n",
    "\n",
    "что эквивалентно \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$ – вектор невязки.\n",
    "\n",
    "- Так как $(I - x_j x^*_j) t  = t$, мы можем переписать это уравнение в симметричной форме"
   ]
  },
  {
   "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",
    "- Теперь заменим $\\lambda$ на $R(x_j)$, и получим **уравнение коррекции Якоби**:\n",
    "\n",
    "$$\n",
    " (I - x_j x^*_j) (A - R(x_j) I) (I - x_j x^*_j) t = -r_j.\n",
    "$$\n",
    "\n",
    "Так как $r_j \\perp x_j$, это уравнение совместно, если $(A - R(x_j) I)$ невырождена."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Решение уравнения коррекции Якоби\n",
    "\n",
    "- Обычно уравнение Якоби решается неточно с помощью подходящего Крыловского метода.\n",
    "\n",
    "- Даже неточное решение уравнения Якоби обеспечивает ортогональность $t$ к вектору $x_j$ (почему?), что хорошо для вычислений."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Связь с методом Релея\n",
    "\n",
    "- Если уравнение решается точно, мы получим метод Релея! Давайте это покажем.\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",
    "- Таким образом, так как $(A - R(x_j) I)^{-1}r_j = (A - R(x_j) I)^{-1}(A - R(x_j) I)x_j = x_j$ мы получим\n",
    "\n",
    "$$x_{j+1} = x_j + t = \\alpha (A - R(x_j) I)^{-1}x_j,$$\n",
    "\n",
    "что совпадает с методом Релея с точностью до нормировочной постоянной."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Предобуславливание уравнения Якоби\n",
    "\n",
    "Популярный предобуславливатель для решения уравнения Якоби имеет вид \n",
    "\n",
    "$$\n",
    "\\widetilde K = (I - x_j x^*_j) K (I - x_j x^*_j)\n",
    "$$\n",
    "\n",
    "где $K$ легко обратимая аппроксимация $A - R(x_j) I $.\n",
    "\n",
    "- Нам нужно получить метод решения системы с $\\widetilde K$ в терминах решения системы с $K$.\n",
    "\n",
    "- Мы уже видели, что уравнение\n",
    "\n",
    "$$ (I - x_j x^*_j) K (I - x_j x^*_j) \\tilde t = f $$\n",
    "\n",
    "эквивалентно\n",
    "\n",
    "$$  \\tilde t = \\alpha K^{-1}x_j  + K^{-1}f $$\n",
    "\n",
    "- А сейчас мы забудем о значении $\\alpha$ и будем искать его из требования $\\tilde t\\perp x_j$ для поддержки ортогональности:\n",
    "\n",
    "$$ \\alpha = \\frac{x_j^*K^{-1}f}{x_j^* K^{-1}x_j} $$\n",
    "\n",
    "- Таким образом, для каждой итерации решения уравнения Якоби мы вычислим $K^{-1}x_j$ после чего обновим только $K^{-1}f$ на каждой внутренней итерации Крыловского метода"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Ускорение на подпространстве в JD\n",
    "\n",
    "- На каждой итерации метода мы расширим базис с помощью нового $t$.\n",
    "\n",
    "- А именно $V_j = [v_1,\\dots,v_{j-1},v_j]$, где $v_j$ – вектор $t$, ортогонализованный к $V_{j-1}$.\n",
    "\n",
    "- Затем используется стандартная процедура Релея-Ритца\n",
    "\n",
    "**Исторический факт:** сначала ускорение на подпространстве использовалось в методе Дэвидсона\n",
    "\n",
    "- В этом методе вместо уравнения Якоби решалось уравнение $(\\mathrm{diag}(A) - R(x_j)I)t = -r_j$\n",
    "- Метод Дэвидсона был очень популярен при решении задач квантовой химии."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Блочный случай в методе JD\n",
    "\n",
    "Если мы хотим найти несколько собственных векторов, мы вычислим **частичное разложение Шура:**\n",
    "\n",
    "$$A Q_k = Q_k T_k, $$\n",
    "\n",
    "и тогда хотим обновить $Q_k$ с помощью одного вектора, добавленного к $Q_k$. Мы будем использовать вместо матрицы $A$ матрицу $(I - Q_k Q^*_k) A (I - Q_k Q^*_k)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### JD: итог\n",
    "\n",
    "- Уравнение коррекции может быть решено неточно и метод JD часто самый быстрый."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Реализации\n",
    "\n",
    "- [ARPack](http://www.caam.rice.edu/software/ARPACK/) наиболее широко используемый пакет для решения частичной задачи на собственные значения (он также запускается внутри SciPy). Содержит варианты алгоритмов Ланцоша и Арнольди\n",
    "- [PRIMME](https://github.com/primme/primme) – лучший по нашему опыту (использует динамическое переключение между различными методами, включая LOBPCG и JD)\n",
    "- [PROPACK](http://sun.stanford.edu/~rmunk/PROPACK/) работает хорошо для вычисления SVD."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Резюме\n",
    "\n",
    "- Методы Арнольди и Ланцоша. Shift-and-invert стратегия очень затратна, поскольку надо очень точно решать систему на каждой итерации.\n",
    "- Итерационные методы с предобуславливателем (PINVIT, LOBPCG, JD) подходят для случая неточного решения системы. \n",
    "- Есть готовые пакеты для их использования\n",
    "- Большое количество технических сложностей осталось за кадром (рестарты, устойчивость)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "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": 41,
     "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": {
  "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}