lectures/lecture7/lecture7.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Лекция 7. Симметричная задача на собственные значения и вычисление сингулярного разложения"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## На прошлой лекции\n",
    "\n",
    "- Собственные векторы и собственные значения\n",
    "- Степенной метод и его сходимость\n",
    "- Разложение Шура и QR алгоритм"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## План на сегодня\n",
    "\n",
    "- Ускорение QR алгоритма\n",
    "- Концепция сдвигов\n",
    "- Подход \"разделяй и властвуй\" к вычислению сингулярных чисел\n",
    "- Бисекция\n",
    "- Метод Якоби"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Сложность вычисления разложения Шура\n",
    "\n",
    "- Напоминаем, что мы хотим избежать сложности $\\mathcal{O}(n^4)$, которую получили при наивной реализации QR алгоритма\n",
    "- Это можно сделать, приведя матрицу к более простой форме\n",
    "- Эта форма призвана удешевить одну итерацию QR алгоритма"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Верхне-гессенбергова форма матрицы\n",
    "\n",
    "Матрица $A$ имеет верхне-гессенбергову форму, если\n",
    "\n",
    "$$a_{ij} = 0, \\; \\mbox{при  } i \\geq j+2.$$\n",
    "\n",
    "$$H = \\begin{bmatrix} * & * & * & * & * \\\\ * & * & * & * & * \\\\ 0 & * & * & * & *\\\\ 0 & 0 & * & * & *\\\\ 0 & 0 & 0 & * & * \\\\ \\end{bmatrix}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Приведение произвольной матрицы к верхне-гессенберговой форме\n",
    "\n",
    "С помощью отражений Хаусхолдера можно привести любую матрицу к верхне-гессенберговой форме:\n",
    "\n",
    "$$U^* A U = H.$$\n",
    "\n",
    "- Единственное отличие от вычисления разложения Шура заключается в занулении последних $n-2, n-3, ...$ элементов в первом, втором и так далее столбцах\n",
    "\n",
    "- Сложность такого приведения $\\mathcal{O}(n^3)$ операций\n",
    "\n",
    "- Если матрица приведена к верхне-гессенберговой форме, то одна итерация QR алгоритма имеет сложность $\\mathcal{O}(n^2)$ операций (например, используя вращения Гивенса) (почему?)\n",
    "\n",
    "- Также верхне-гессенбергова форма матрицы сохраняется после выполнения одной итерации QR алгоритма (проверьте!)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Случай симметричной (эрмитовой) матрицы\n",
    "\n",
    "- Если матрица $A$ симметричная (эрмитова), то $A = A^*$, тогда $H = H^*$ и верхне-гессенбергова форма оказывается трёхдиагональной матрицей\n",
    "\n",
    "- Далее мы будем говорить только о симметричном трёхдиагональном виде верхне-гессенберговой формы\n",
    "\n",
    "- Любая эрмитова матрица может быть приведена к трёхдиагональной форме с помощью отражений Хаусхолдера\n",
    "\n",
    "- **Основная идея:** трёхдиагональная форма сохраняется при выполнении QR алгоритма, и сложность одной итерации может быть сокращена до $\\mathcal{O}(n)$!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Одна итерация QR алгоритма\n",
    "\n",
    "- Одна итерация QR алгоритма имеет следующий вид:\n",
    "\n",
    "$$A_k = Q_k R_k, \\quad A_{k+1} = R_k Q_k.$$\n",
    "\n",
    "- Если $A_0 = A$ <font color='red'> симметричная трёхдиагональная матрица </font>, то эта форма сохраняется.\n",
    "\n",
    "- Давайте это проверим!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.25844981, -1.86142358,  0.        ,  0.        ,  0.        ,\n",
       "        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,\n",
       "        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,\n",
       "        0.        ,  0.        ,  0.        ,  0.        ,  0.        ])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAEQCAYAAAB2hgJTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAPnUlEQVR4nO3df6xkZX3H8fcHdtHyq/LjLiW2sGgR0RY37HXbBIJLjJZiK1BbBG1LK3ZpCmmrQUtoG9b+YWnVoG0qZtXNYlKgpEKBQCgUtdsosVwN0RW6soWLRZbdu4DugqDs8u0f57k4DDP3zp05c8/c+X5eyWTuPXPmPN8z937mOT+eOaOIwMzG335NF2Bmi8NhN0vCYTdLwmE3S8JhN0vCYTdLwmFfRJI2Ser7XKeklZJC0voay5qrvfWlvZWL0d5iGdf1mo/DPgdJh0r6a0nflLRH0o8k3S/pY5KOarq+zCS9qoR2bdO1LBXyoJrOJL0O+HfgWOBG4MvA88CvAr8L7AZ+MyLuWcAylwP7R8RzfdYk4BXA3ojY288yFtjeeuAK4LiImB52ewtReuWHgY9ExPoFPncZsAz4cSQKwLKmCxhFkg4EbgVeTRXo21oe3iDp08B/ADdL+uWI2DHHsgQcFBFPR8TzVG8YfSn/mH29URhIOiQi9pQ3yqG/WY4ab8Z3diHwOuCqtqADEBFTwOXABPCh2emS1pZ9wT+QdLGk+6nCeWl5vOM+u6S3SLpH0rOSHpf0KUlvbN8/77TP3jpN0m9IulfSc5K2l92NZW1trSl1fLfsluyR9FVJ5wzygkn6iqTpUs9Nkn4g6anS1sGS9pN0uaSHS33flHRK2zL2k/SXkjaX1+Enkr4n6WpJR7S+zlS9OsAVZf1D0nSH1+Tdkr4h6VngH8vjL9lnL7tr2yQ9JmlFW00fLfO+b5DXZxS4Z+/st8v9Z+eYZxPwSeBdlDC3+HPgiPL8x4H/67YQSacCdwJPAVcCPwDOBU7p9pwuzgT+BPgMsBE4q9T1FPDRlvnOAV4P3AA8Uuq8ALhR0nsj4toFttvqIOBLwGbgMuDNwPuAVwJPAL9CFbjlpbZbJR0bEXvK8w+gevP8InAz8ExZxoXAqZJWR8RPgAeADwBXATdR7WYBPN1Wz9nAnwJXl9dld6eiI2K3pPOArwLXSDozIkLSW4G/AK6PiI39vigjIyJ8a7tR/WPu7mG+bwMBHFx+X1t+fxJY0WH+TZSt8ZZp/03V+7+mZdpyqn+8ANa3TF85x7RngJUt0wVsAba3tXdQh7oOBLYC97dNX1+WvbLT+rfN+5Uy74fapt8IvABMActbpr+zzH9RW80/02HZF5Z5z53rtejw2PPAiR0e77hewAfL9Eupttq2Aw8Bhzb9P1nHzZvxnR0K/LCH+Wbn+dm26V+IiJ3zPbkc0X8zcHNEPDQ7Pap9+0/1WOusf4uWg2hR/fd+Gfg5SQe3TH+mpf0Dy+bxgVQ98omSDl1gu632UTaVW/wXVYg/U9ardTrA8a01R8Szpbb9yxH3I0ttUG0ZLMRtEfHAAua/CridakvoNqqtnvMjouMWwVLjsHe2myrw85mdp/2N4bs9tnNcud/a4bFO0+byUIdpT5T71v3dFZI2SNpBtTWwC5gB/rjM8qoFtttqe7z8TMNT5f7h1okRMTv9iNbpks6V9HXg2fLcGX66boctsJ5e/w6zNQXVLs3s7sP6iPj6AtscWd5n72wLcJqkX4yIbZ1mKEfsTwCmI6J9X/FHPbajAWpst2++dsqZgTuBE4F/AO6leqPaB/wh8B4G6wDmqqHbYy++BpJ+C/gXql2bP6M61vEcsD9wRx+19fp3aHUaP33DW9XH80eWw97ZjVR/9PdTHWjq5PepDijd2OXxXsz2WCd0eKzTtEGdBLwJ+JuIuKL1AUnvH0J7C/V7VOE+PSJeDKqk13eYt/bz45KOAT5H9WZ/J/BBSX8UEXMdqF0yvBnf2eeAbcAHJJ3R/qCkk4G/pdrE/Fi/jUR1fn4KOEvSa1qWv5yqZ6vbbO/6ki0KSb9EdZS+afuoQvzi/2XZGvmrDvPObk0dXkfDkvYHrqU6c/BuqqPw9wCflHRiHW00zT17BxHxjKR3Um063ibpi1RHm/cCa6h6oKeBsyPi8QGbuxS4C/haGazzQ6pTbwfMljPg8ls9AHwH+HDZDdlKNZ7gIqre7OQa2+rHv1KdyvySpC9QnZU4m+oA4ktExBOStgHnSfpfYAfwTETc2mfb66lOd66LiPsBJL0HuA+4XtKaiPhxn8seCe7ZuyhHcU8CPkIViL+nOkJ+KtUR5zdGxNdqaOc/gTOAaaqBOpdT9faXlFmeHbSNlrb2Ae+gGh14AdX6vKX83G9IahMR1wPrgIOBjwMfpnpD+rUuT3kv8CDV0fPrePmZgJ6UQTqXAze0brKXsxvrqP4PPt7PskeJx8aPKEnvourpzi8hMBuIe/aGqfLKtmnLqQZ47KXafTAbmPfZm/cK4BFJ/0y1yXoE1QGik4C/q+GYgBngsI+C56lGa50FHE11pHwrcHFEfLrJwmy8eJ/dLAnvs5sl0WjYJZ0haWv5LHG3kWpLTvlc97cl3Sdpqul6BiFpo6Sdkra0TDtc0l2SHiz3Cx2z3rgu67Ve0vfL3+0+SWc2WWPdGgt7GbH0T8CvA28Azpf0hqbqGYLTI2JVREw2XciANlGNA2h1GXB3RBwP3E33IcWjbBMvXy+oLliyqtxuX+SahqrJnn0NsC0iHorqggTXUx2kshESEZupPp/f6izgmvLzNVSj3JaULus11poM+6t56RVcHi3TxkEAd5bLIa1rupghOCoitgOU+xXzzL+UXCLpW2Uzf8ntnsylybB3+njnuJwaOCUiTqbaRblY0mlNF2Q9uRp4LdVHW7cDn2i0mpo1GfZHgV9o+f3ngccaqqVWEfFYud9JdY20Nc1WVLsdko4GKPfzXpVnKYiIHRGxLyJeoLp+4Fj93ZoM+73A8ZKOk3QAcB5wS4P11ELSQZIOmf0ZeDvVJ8rGyS1UH56h3N/cYC21mX0DK85hzP5ujY2gi4i9ki6h+iKG/YGNEfGdpuqp0VHATdXHsFkGXBsRdzRbUv8kXUd1Ic0jJT1K9aURVwI3SLoQ+B7wO81V2J8u67VW0iqq3clpqo/+jg2PoDNLwiPozJJw2M2ScNjNknDYzZJw2M2SaDzsYzqcFBjfdfN6LU2Nh53q6p3jalzXzeu1BI1C2M1sESzqoBpJQ2ls9erVw1jswGZmZpiYmGi6jNp5vUbX9PQ0u3bt6vgdgmNxwcmpqSV9MRiz2kxOdr9WykCb8eN6WSmzcdR32BNcVspsrAzSs/uyUmZLyCBh7+myUpLWSZpa6ldZNVvqBjlA19NlpSJiA7ABhnc03szmN0jPPraXlTIbR4OEfSwvK2U2rvrejB/jy0qZjaWBBtWUb8xo/FszyvXeeuLLcFlWHhtvloTDbpaEw26WhMNuloTDbpaEw26WhMNuloTDbpaEw26WhMNulsSihn316tVERE+3YZHU881snLhnN0vCYTdLwmE3S8JhN0vCYTdLwmE3S8JhN0vCYTdLwmE3S8JhN0tiZL+yeSFDZoc1tLXX5fqKtbYUuGc3S8JhN0vCYTdLwmE3S8JhN0vCYTdLwmE3S8JhN0vCYTdLwmE3S2Jkh8suRNNDaxeyTA+ttaa4ZzdLwmE3S8JhN0vCYTdLwmE3S8JhN0vCYTdLwmE3S8JhN0tiLEbQLYRH21lW7tnNknDYzZIYaDNe0jSwB9gH7I2IyTqKMrP61bHPfnpE7KphOWY2RN6MN0ti0LAHcKekb0haV0dBZjYcg27GnxIRj0laAdwl6X8iYnPrDOVNYB3AMcccM2BzZtavgXr2iHis3O8EbgLWdJhnQ0RMRsTkxMTEIM2Z2QD6DrukgyQdMvsz8HZgS12FmVm9BtmMPwq4qYwIWwZcGxF31FKVmdWu77BHxEPAm2qsZeT0Oly16e+HBw+ttfn51JtZEg67WRIOu1kSDrtZEg67WRIOu1kSDrtZEg67WRIOu1kSDrtZEumuLjsMTV+xdqHL9dDanNyzmyXhsJsl4bCbJeGwmyXhsJsl4bCbJeGwmyXhsJsl4bCbJeGwmyXh4bKLzENrrSnu2c2ScNjNknDYzZJw2M2ScNjNknDYzZJw2M2ScNjNknDYzZJw2M2S8HDZEeahtVYn9+xmSTjsZkk47GZJOOxmSTjsZkk47GZJOOxmSTjsZkk47GZJOOxmSXi47Jjw0Fqbj3t2syTmDbukjZJ2StrSMu1wSXdJerDcHzbcMs1sUL307JuAM9qmXQbcHRHHA3eX381shM0b9ojYDDzZNvks4Jry8zXA2fWWZWZ163ef/aiI2A5Q7lfUV5KZDcPQD9BJWidpStLUzMzMsJszsy76DfsOSUcDlPud3WaMiA0RMRkRkxMTE302Z2aD6jfstwAXlJ8vAG6upxwzG5ZeTr1dB9wDnCDpUUkXAlcCb5P0IPC28ruZjbB5R9BFxPldHnprzbWY2RB5uGxCHlqbk4fLmiXhsJsl4bCbJeGwmyXhsJsl4bCbJeGwmyXhsJsl4bCbJeGwmyXh4bI2Jw+tHR/u2c2ScNjNknDYzZJw2M2ScNjNknDYzZJw2M2ScNjNknDYzZJw2M2S8HBZq81SGlqbcVite3azJBx2syQcdrMkHHazJBx2syQcdrMkHHazJBx2syQcdrMkHHazJDxc1hrR9NDajFesdc9uloTDbpaEw26WhMNuloTDbpaEw26WhMNuloTDbpaEw26WhEfQ2cjzaLt6uGc3S8JhN0ti3rBL2ihpp6QtLdPWS/q+pPvK7czhlmlmg+qlZ98EnNFh+lURsarcbq+3LDOr27xhj4jNwJOLUIuZDdEg++yXSPpW2cw/rNtMktZJmpI0NTMzM0BzZjaIfsN+NfBaYBWwHfhEtxkjYkNETEbE5MTERJ/Nmdmg+gp7ROyIiH0R8QLwWWBNvWWZWd36Cruko1t+PQfY0m1eMxsN846gk3QdsBY4UtKjwBXAWkmrgACmgYuGV6KZ1WHesEfE+R0mf34ItZgNrNfhqk1/Pzws/tBaj6AzS8JhN0vCYTdLwmE3S8JhN0vCYTdLwmE3S8JhN0vCYTdLwmE3S8JXl7WUmr5i7UKXW8fQWvfsZkk47GZJOOxmSTjsZkk47GZJOOxmSTjsZkk47GZJOOxmSTjsZkl4uKzZPJba0Npu3LObJeGwmyXhsJsl4bCbJeGwmyXhsJsl4bCbJeGwmyXhsJsl4bCbJeHhsmY1GoWhtd24ZzdLwmE3S8JhN0vCYTdLwmE3S8JhN0vCYTdLwmE3S8JhN0vCYTdLQnV8yXvPjUkzwCNtk48Edi1aEYtrXNfN6zW6jo2IiU4PLGrYOxYgTUXEZKNFDMm4rpvXa2nyZrxZEg67WRKjEPYNTRcwROO6bl6vJajxfXYzWxyj0LOb2SJw2M2ScNjNknDYzZJw2M2S+H8gemCd25mfgAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUcAAAEQCAYAAAAjwrYkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAWKElEQVR4nO3debRdZXnH8e8DEVFERRMQrRIESp1qlIiu5RRUFGktYmtbR7QgTrS1da6txiXOWodVtQTFIBZEZVSpBYcspOIQJ8QighhqhJJEHABFBJ7+8b6XnJy8595zxzPc72ets869++yz97un33738O4TmYkkaVs7DLoAkjSMDEdJajAcJanBcJSkBsNRkhoMR0lqmJdwjIiDIuJrEXFdRGREPG8+xjMu6jxaO+hyzFRELK/TsHrQZelXROwQEasj4oqIuDkixuaetoh4Xl0eq0atLMNU9mmFY0TsFhE31sI/u1c/wOnALsDLgecA59eJftlsC7xY1A33qYMux4QagKsjYsWgyzJHjgDeAHwZOJKynk4pIg6LiM9FxKaIuCkiroqIT0fEo3v0v65uLxOv39fvnBoRD5y7yRkdEbGirkvLB12WySyZZv/PAnYCfkJZoT7e6OdhwF2BIzPz9ImOEXECsBx47wzKOe7uANzS1e0NwInAmQtemrbllDJtAL7b9dmVlGm4eUFLNDsHA78Cjso+WkJExI7ARykh+j/A+4Crgb2A51IqAMdm5r80vv474Kj69x2AA4DnA4dGxMrMvHS2EzPETgI+AdzU0W0FZV1aR1mfhtJ0w/FIyp72LOC9EbFPZv64q5971PdrZ1u4fkXEHYDfZ+YobZy3ycwbF3J8EbFrZl43V8Or4bKg0zAH7gH8sp9grFZTgvFjlB3/betaRLydsk38c0Rcnpkndn335szsrEgcHxETAXsM8LcznIahNbGOZeYtbL/jHw2Z2dcLeCiQlL3k3Sl7w2O7+tlQ++l+9eq+quO7+1H2MldT9jIbgHcCu3SNY2397jLgBOAa4FZg+RTlXwp8APhpHf5P6/937+rveXX4jwNeAfy4TuuPgCN6DPsJwLnALykhcRHwomnM2wTW1r+X95hXOZNx1vm4DngI8F+U2tJP6me7AscCXwe21Om8HHgbcMfGPOl+resq8+qucS8BXk2pad0I/Bw4A3hQV3+3fR/4U+Cbtf+r6zqwZBrz8ijg28Bv67SeCzyq4/NVPaZl7STD3L0O70pg50n6uR74GbBTR/d1wPWN/h9Qx/v5Pqfrj4APAj8ArgN+A3wLeEGj34nltaqr+3LgNODXdd6cBew9sY5Md152r7/A44EL6nxY1ypLXcY95z/T3P46xv044MI6XzYCr66f7wZ8BNhUP/sscM9+5vl0ao5HAjcAp2XmDRHxOeCIiHh9Zt5a+3kZ8GTgaOAtwCW1+/XAWykB9Q8dw7wEICIOAL5E2dCPo6xgDwb+DnhkRDw2M3/fVZ7zgP8D3kQ5v3l9r4JHxF2ArwL7UgL125SweDHwuIg4MLevSb2Fcgh0HGXhvBhYW2sG/90x7KOBfwe+Bry5zqODgQ/VmvUre5Wrh82UGspJwFeANY3pme4470OZv5+ibBx3qt3vRdkATgNOphwWPxZ4FWX+PKn2d36dH/9Uy/OV2v2aKablP4C/pCyrD1Fqay8FLoyIR2fmd7r6PxR4SZ22E4DDKBvIL+r4J1VrcK8CvlHLuitlXfxyRByWmedQ1rnnAK9j2/Wx+wio058AOwMnZY9afmZuioizgGcCj6DMs8nsU9/7PcJaBTyGsnH/hLLOPx1YExFLM/Otk305Iu5OWW57UObvJcCjKUeCuzT672dedloJ/DlwPOV0UC+nA3uyfUZ0z/++tr/qIcBTKOvmxyjr3Nsi4kbKueUNlFDel5IpH6NULibX515rZ8pCXNvR7TBKaj+5z73WOmBDj+F/D/ghsGtX98PrsJ7X0W1t7fbxadQm3ly/85Ku7i+t3d/UKP932LYGcK+6kE7p6LYnpYZzcmOc76McTuzTR/m2q7m0us1knGyttR/V6H8n4HaN7m+q3zmwo9uq7mXRVSPZpuZICesETgWio/sfU0L4K43v30DHEQAQwMXA1X3Mw/0pRxAXdC23e1J2uhuAHftZHxvDfnct39Om6O/ltb9jusZzPSWIlwL3Bp7asVwO7bMMuzS67VCH/6vO5UhjGwTeUbs9q2sYE93XzWJeTtT+ntAoY6ss23WbyfbXMe5bgYd3rddX1+7v7+r/X+t39p9qnvd7tfpplOpp5x7hc5Sq6t/0OYymiHgQZYM5Gbh9RCydeFEWzg3AExtffdc0RnM4pUbWXQs7jnI4eXjjOx/MzNtOImfmzyhV+/06+vkL4PbARzrLXcv+GcrK+/hplLMfMxnntZSLCdvIzJuy1sgjYkm9G2Ep8IXay8NnUc6JefrmrGtlHedFlNrPoyJiWdd3zszMDR39JqVmc4+IuBOTO4wSpu/oWm5XUXaoe1FqGDNx5/r+qyn6m/h8167uu1DWv83A/1JOLexEOUzsroE1ZeYNE39HxM61Jng3yqHunSmH3ZN5CiUwTunq3tqOZjIvv5eZX2Du9LP9TbgwM7/e0e9NlBpvAO/v6nfiqKc1nG30e1h9JGXBboyIfTu6nwc8vVbrt/Q5rG73q+9vrK+WPRrdfjSNcewNrM+uCzaZeXNEXEo5n9rtika3n1NWjAkTZZ9spWiVfTZmMs4fZzkxvp2IeAnwIso5sO6d5W4zKmGxN2XPfUnjs4spG+DelPVqQq95DuU8d89TJ3VYUM7JtcYHcF9g/STD6OXX9f0uU/Q3EaLdpxtupIQTlEB7LqVm3fetdHXnsJpyyHjvRi9TLau9gW/k1lNgwG2nA37Z6BemNy+nsz32o5/tb7J+f1Hff9Kj+92nKsCU4RgRewMHUVK41wx4NjO/RSfq+7uBz/fo5xfdHTLzNzMcX796XWGLxt/PpeyVW1oLbjZmMs7mvIqIf6TM93Mpe9irKBer7kWpIcymkUBM3ct2JruqOdXwZjK+fk0EwkMp58x6mdjJXt7V/ZbOWlVEfJpSe14TEd+utempnEy5WLWGcj7zWsrpiUMp503nskHHTOblXG+P/Wx/U/VLr0pBj+Fso5+a4/PrgF5AOd/Q7VhKzfK9Uwwne3S/rL7fMsfV8k5XAPtHxJLc9haMJcAfMvMAmyj7lnks+3yO8zmU80dP7qxRRMQhjX57Lb9efky5oHM/ypX0Tvev79179dmYOKH/ALY/uT8xvpku589Ran/Pjog3ZebvunuopwgOo1zRvmCygWXmrRHx95Sr+O+ifdqoc9h3pQTjSZn5oq7Ppr6wUGwA9o2IHbqW9e6U+5I7zee8hOmvSwMx6d4mInagnCD9fmZ+ODM/3f2inMN4YEQ8bIpxXQ/sFhHdif0dyp75RRFx30YZlkTE3fqdoB7OpNz6c1RX9xfU7mfMcLifpJwkfmO913IbEXGXiLj9DId9PeUQbD7HeQtlRb1tmdQdxmt6lIceZWo5s76/tnOZ11YhfwZckJmbW1+cobMp0/LKiLhdx/j2pOzgr6Ssa9OWmZsoNey9KHcE7Nj5eV0OJ1HuAnh996Frj2FeRqkNHhwRj5qi94nazzbbTp227nW6l89QLuY9o6v7Kxr9ztu8rKa7Lg3EVDXHJ1LOb3xkkn5Oo5wLOZJyf1ovX6Ps/f4tIr5KWeBfquc8nkO51eSi2pLmB8AdKZfenwa8lnKYN1PvoNz28IGIeChlwT6klvnS+vm0ZebGiHgx8GHgkog4ibLiLAMeRLkqeX9m1grga8ATIuLVlJP4mZmfmONxfppyi9V/RsTplHNmzwS6b5uCUsu5DnhJRPyGchSxKTO/1BpwZp4XEZ8E/pqyU/wsW2/luZFyS8WcycxLI+KdlNtPzo+IU9l6+8mdKFdpZ3Mz8hso4fh84GERcQrbtpDZi3K1/mPTGOZbKKek3sgkF+4y87qIOJdSc/0tZTvbC3ghpfY95fkz4O2UZfvRiDiQcnfIo4BHUi5Kdl40m+95+U3K+ejXRWlufAPl3tuvT/61BTbF7QOfosy0B03R36WUjeUO9L6VZxdKyF7D1hpL5+X9vSj3X22gnPf6OeUm17cC9+7oby1dN0T386KExwcpN4j+vr5/AFg61a0HHZ+to3H7B2UFO4Ny9f4myrm7L1Nu7WjeNNz1/datPPtRzgX+mvZN4H2Nkx43+NbPdqTseC6n1EavpOwo7kf7pu5DKfeI3sj0bgK/pA7/WkqNsudN4I0yrq6fLe9zOb+AsvO7sc6784BH97ss+xj+U4FzKBeSJtbjW4En9uh/HY2bwDs+P6UO47FTjHcpZYd4VZ2279dp3W597bUOUy60nE7Zyf2arTeBbwHOmcW83G797aMsR1B2uDd1fr9X/72WWa9x0yMnmOSWtO5X1C9ImoGI+CvKze5fB56UmZNdUR869ZagLcBx2XU+c7HzeY7SLGTmqZR7fR8BnN06DzwsepTt1fX9vIUsyyiw5igtEhGxjnLqZD3llMrjKdcBvgo8Jmd3HnHsGI7SIhERL6dcPFpOuT6wkXIO8o05h09pGheGoyQ1eM5RkhoWVThGxCERcWlEXB4RrRudR1JEbIiI70fEdyNiJm2Hh0ZEnBDlJwgu7uh2t4g4LyIuq++zafM9ED2ma3VE/Kwut+9GxKGDLKO2tWjCsbZq+ADleZP3B54REfef/Fsj5aDMXJGZKwddkFlaC3Q3X3wN8MXM3A/4Iu0WPMNuLdtPF8B76nJbkX0+oUcLY9GEI3AgcHlmXpHlkUafoLSF1RDJzImHKnQ6jK2PyzuRciP2SOkxXRpiiykc70X5aYQJG2u3cZDAuRHxrShPCR83e2Tm1QD1ffcBl2cuHRMRF9XD7pE7XTDOFlM4th5RNC6X6h+ZmQ+lnDJ4aUQ8ZtAFUl8+RPm5hBWUdtrvHmhptI3FFI4b2fYhoX9Aaac68rI8oZksT485g3IKYZxcU58IM/FkmE0DLs+cyMxrMvOWLE/xOZ7xW24jbTGF4zeB/SJi74jYifK0mLMHXKZZi4hdImLXib8pT1K6ePJvjZyzKQ8qoL6fNcCyzJmJwK8OZ/yW20ib7u9Wj6wsP4lwDOXnSXcETsjM1mPgR80ewBn1kYlLKD+81euJ6kOvPgpsFbA0IjZSHhX2NuCTEXEk5fFtTx9cCWemx3StiogVcNvPF79wUOXT9mwhI0kNi+mwWpL6ZjhKUoPhKEkNhqMkNRiOktSw6MJxTJvXAeM7bU6XBmHRhSPl5yXH1bhOm9OlBbcYw1GSpjTWN4FHxLxM3AEHHDAfg521zZs3s2zZskEXY845XcNrw4YNbNmypfVQl5G3aJoPzqX160f6YdvSnFm5ctSfrdzbSB1Wj+vPHEgaPiMTjovgZw4kDZGRCUf8mQNJC2iUwrGvnzmIiKMjYv2o/wqfpMEapQsyff3MQWauAdbA/F2tljT+RqnmOLY/cyBp+IxSOI7lzxxIGk4jc1g9xj9zIGkIjUw4AmTmOcA5gy5H/b2WvoxzCyRpnI3SYbUkLRjDUZIaDEdJajAcJanBcJSkBsNRkhoMR0lqMBwlqcFwlKQGw1GSGsY6HA844AAys6/XfImIvl+ShsdYh6MkzZThKEkNhqMkNRiOktRgOEpSg+EoSQ2GoyQ1GI6S1GA4SlKD4ShJDSP164PzaTpNCOerqV+/w/UXDaX5Z81RkhoMR0lqMBwlqcFwlKQGw1GSGgxHSWowHCWpwXCUpAbDUZIaDEdJarD54AwMuqnhdIZpU0NpZqw5SlKD4ShJDYajJDUYjpLUYDhKUoPhKEkNhqMkNRiOktRgOEpSgy1k5pmtaaTRZM1RkhoMR0lqGKnD6ojYAFwH3ALcnJkrB1siSeNqpMKxOigztwy6EJLGm4fVktQwauGYwLkR8a2IOHrQhZE0vkbtsPqRmXlVROwOnBcRP8zM8zt7qKF5NMB97nOfQZRR0hgYqZpjZl5V3zcBZwAHNvpZk5krM3PlsmXLFrqIksbEyIRjROwSEbtO/A08Ebh4sKWSNK5G6bB6D+CM2uJjCXByZn5+sEWSNK5GJhwz8wrgwYMux3zqt/nefDQznO5wbWqocTcyh9WStJAMR0lqMBwlqcFwlKQGw1GSGgxHSWowHCWpwXCUpAbDUZIaDEdJahiZ5oPaatC/aDjd4drUUKPImqMkNRiOktRgOEpSg+EoSQ2GoyQ1GI6S1GA4SlKD4ShJDYajJDUYjpLUYPPBMWdTQ2lmrDlKUoPhKEkNhqMkNRiOktRgOEpSg+EoSQ2GoyQ1GI6S1GA4SlKD4ShJDTYf1G1saihtZc1RkhoMR0lqMBwlqcFwlKQGw1GSGgxHSWowHCWpwXCUpAbDUZIaDEdJarD5oGbEpoYad9YcJalh6MIxIk6IiE0RcXFHt7tFxHkRcVl9322QZZQ0/oYuHIG1wCFd3V4DfDEz9wO+WP+XpHkzdOGYmecD13Z1Pgw4sf59IvDUhSyTpMVn6MKxhz0y82qA+r77gMsjacyNSjj2LSKOjoj1EbF+8+bNgy6OpBE1KuF4TUTsCVDfN/XqMTPXZObKzFy5bNmyBSugpPEyKuF4NnBE/fsI4KwBlkXSIjB04RgRpwAXAvtHxMaIOBJ4G3BwRFwGHFz/l6R5M3QtZDLzGT0+evyCFkTSojZ04ajxY1NDjaKhO6yWpGFgOEpSg+EoSQ2GoyQ1GI6S1GA4SlKD4ShJDYajJDUYjpLUYDhKUoPNBzVUbGqoYWHNUZIaDEdJajAcJanBcJSkBsNRkhoMR0lqMBwlqcFwlKQGw1GSGgxHSWqw+aBG1ig1NbSZ4eix5ihJDYajJDUYjpLUYDhKUoPhKEkNhqMkNRiOktRgOEpSg+EoSQ2GoyQ12HxQi8Kgmxr6i4ajx5qjJDUYjpLUYDhKUoPhKEkNhqMkNRiOktRgOEpSg+EoSQ2GoyQ12EJG6mJrGoE1R0lqMhwlqWHowjEiToiITRFxcUe31RHxs4j4bn0dOsgyShp/QxeOwFrgkEb392Tmivo6Z4HLJGmRGbpwzMzzgWsHXQ5Ji9vQheMkjomIi+ph9269eoqIoyNifUSs37x580KWT9IYGZVw/BCwD7ACuBp4d68eM3NNZq7MzJXLli1boOJJGjcjEY6ZeU1m3pKZtwLHAwcOukySxttIhGNE7Nnx7+HAxb36laS5MHQtZCLiFGAVsDQiNgJvAFZFxAoggQ3ACwdVPkmLw9CFY2Y+o9H5IwteEKkP/Tbfm49mhtMdrk0Np2ckDqslaaEZjpLUYDhKUoPhKEkNhqMkNRiOktRgOEpSg+EoSQ2GoyQ1GI6S1DB0zQelcTToXzSc7nBtamjNUZKaDEdJajAcJanBcJSkBsNRkhoMR0lqMBwlqcFwlKQGw1GSGgxHSWqw+aA0ZEatqeG4suYoSQ2GoyQ1GI6S1GA4SlKD4ShJDYajJDUYjpLUYDhKUoPhKEkNhqMkNdh8UBphw9DUcFxZc5SkBsNRkhoMR0lqMBwlqcFwlKQGw1GSGgxHSWowHCWpwXCUpAbDUZIaYjrNj0ZNRGwGruzqvBTYMoDiLIRxnTana3jtlZnLBl2I+TDW4dgSEeszc+WgyzEfxnXanC4NgofVktRgOEpSw2IMxzWDLsA8Gtdpc7q04BbdOUdJ6sdirDlK0pQMR0lqMBwlqcFwlKQGw1GSGv4fbxUl7GfP9psAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#Generate a random tridiagonal matrix\n",
    "n = 20\n",
    "d = np.random.randn(n)\n",
    "sub_diag = np.random.randn(n-1)\n",
    "\n",
    "mat = np.diag(d) + np.diag(sub_diag, -1) + np.diag(sub_diag, 1)\n",
    "plt.spy(mat)\n",
    "plt.title(\"Original matrix\", fontsize=18)\n",
    "q, r = np.linalg.qr(mat)\n",
    "plt.figure()\n",
    "b = r.dot(q)\n",
    "b[abs(b) <= 1e-12] = 0\n",
    "plt.spy(b)\n",
    "plt.title(\"After one iteration of QR algorithm\", fontsize=18)\n",
    "#plt.figure()\n",
    "#plt.imshow(np.abs(r.dot(q)))\n",
    "b[0, :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Трёхдиагональная форма\n",
    "\n",
    "- Работая с трёхдиагональной формой, вам не нужно вычислять матрицу $Q$: нужно лишь вычислить трёхдиагональную часть, которая получается после итерации \n",
    "\n",
    "    $$A_k = Q_k R_k, \\quad A_{k+1}  = R_k Q_k,$$\n",
    "\n",
    "    в случае $A_k = A^*_k$.\n",
    "\n",
    "- Такая матрица определяется $\\mathcal{O}(n)$ параметрами\n",
    "- Вычисление QR разложения более сложное, но возможно вычислить $A_{k+1}$ напрямую без вычисления $Q_k$.\n",
    "- Это называется **неявный QR-шаг**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Теорема о неявном QR алгоритме\n",
    "\n",
    "\n",
    "- Все реализации неявного QR алгоритма основаны на следующей теореме \n",
    "\n",
    "**Теорема.** Пусть \n",
    "\n",
    "$$Q^* A Q = H$$ \n",
    "\n",
    "верхне-гессенбергова форма матрицы.\n",
    "Тогда первый столбец матрицы $Q$ определяет все остальные её столбцы. Он может быть найден из следующего уравнения\n",
    "\n",
    "$$A Q = Q H. $$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Сходимость QR алгоритма\n",
    "\n",
    "- Сходимость QR алгоритма – непростой вопрос (см. Е.Е. Тыртышников \"Краткий курс численного анализа\")! \n",
    "\n",
    "\n",
    "**Итог.** Если у нас есть разложение вида\n",
    "\n",
    "$$A = X \\Lambda X^{-1}, \\quad A = \\begin{bmatrix}A_{11} & A_{12} \\\\ A_{21} & A_{22}\\end{bmatrix}$$\n",
    "\n",
    "и\n",
    "\n",
    "$$\n",
    "\\Lambda = \\begin{bmatrix} \\Lambda_1 & 0 \\\\ \n",
    "0 & \\Lambda_2 \\end{bmatrix}, \\quad \\lambda(\\Lambda_1)=\\{\\lambda_1,\\dots,\\lambda_m\\}, \\ \\lambda(\\Lambda_2)=\\{\\lambda_{m+1},\\dots,\\lambda_r\\},\n",
    "$$\n",
    "\n",
    "а также есть **зазор** между собственными значениями в матрице $\\Lambda_1$ и $\\Lambda_2$ ($|\\lambda_1|\\geq \\dots \\geq |\\lambda_m| > |\\lambda_{m+1}| \\geq\\dots \\geq |\\lambda_r| >0$), тогда блок $A^{(k)}_{21}$  матрицы $A_k$\n",
    "сходится к нулевому в процессе работы QR алгоритма со скоростью\n",
    "\n",
    "$$\\Vert A^{(k)}_{21} \\Vert \\leq  C q^k, \\quad q = \\left| \\frac{\\lambda_{m+1}}{\\lambda_{m}}  \\right |,$$\n",
    "\n",
    "где $m$ размер матрицы $\\Lambda_1$.\n",
    "\n",
    "Таким образом, нам нужно увеличить зазор между $\\Lambda_1$ и $\\Lambda_2$. Это можно сделать с помощью **QR алгоритма со сдвигами**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR алгоритм со сдвигами\n",
    "\n",
    "$$A_{k} - s_k I = Q_k R_k, \\quad A_{k+1} = R_k Q_k + s_k I$$\n",
    "\n",
    "\n",
    "- Сходимость такого алгоритма линейная с фактором\n",
    "\n",
    "$$\\left| \\frac{\\lambda_{m+1} - s_k}{\\lambda_{m} - s_k}  \\right |,$$\n",
    "\n",
    "где $\\lambda_m$ – $m$-ое большее по модулю собственное значение. Если сдвиг близок к собственному вектору, сходимость более быстрая.\n",
    "\n",
    "- Существуют различные стратегии выбора сдвигов. \n",
    "\n",
    "- Использование сдвигов – это общий подход к ускорению сходимости итерационных методов вычисления собственных значений. Далее покажем, как выбирать сдвиги для более простых алгоритмов"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Сдвиги и степенной метод\n",
    "\n",
    "Напоминаем, что такое степенной метод для вычисления собственых значений.\n",
    "\n",
    "$$x_{k+1} := A x_k, \\quad x_{k+1} := \\frac{x_{k+1}}{\\Vert x_{k+1} \\Vert}.$$\n",
    "\n",
    "- Он сходится к собственному вектору, соответствующему максимальному по модулю собственному значению.\n",
    "\n",
    "- Сходимость может быть очень медленной.\n",
    "\n",
    "- Используем сдвиги: если мы преобразуем матрицу\n",
    "\n",
    "$$  A := A - \\lambda_k I,$$\n",
    "\n",
    "то соответствующее собственное значение уменьшится (а нам нужно максимальное по модулю). Это не то, что нам было нужно!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Обратная итерация и итерация Релея\n",
    "\n",
    "- Для того чтобы из малого собственного значения сделать большое, нам нужно **обратить матрицу**, и это приводит нас к **методу обратной итерации**\n",
    "\n",
    "$$x_{k+1} = (A - \\lambda I)^{-1} x_k,$$\n",
    "\n",
    "где $\\lambda$ – сдвиг, который близок к собственному значению, которое мы хотим найти. Аналогично степенному методу сходимость линейная.\n",
    "\n",
    "- Для ускорения сходимости можно использовать **итерацию Релея**, которая задаётся с помощью **адаптивного выбора** параметра сдвига:\n",
    "\n",
    "$$x_{k+1} = (A - \\lambda_k I)^{-1} x_k,$$\n",
    "\n",
    "$$\\lambda_k = \\frac{(Ax_k, x_k)}{(x_k, x_k)}$$\n",
    "\n",
    "В симметричном случае $A = A^*$ сходимость **локально кубическая**, и **локально квадратичная** иначе."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Сингулярные значения и собственные значения (1)\n",
    "\n",
    "\n",
    "- Сингулярное разложение имеет вид\n",
    "\n",
    "$$A = U \\Sigma V^*$$\n",
    "\n",
    "и существует для любой матрицы.\n",
    "\n",
    "- Его также можно считать способом приведения данной матрицы к диагональному виду с помощью двух унитарных преобразований:\n",
    "\n",
    "$$\\Sigma = U^* A V.$$\n",
    "\n",
    "- С помощью двусторонних преобразований Хаусхолдера мы можем привести любую матрицу к **бидиагональной форме** $B$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Сингулярные значения и собственные значения (2)\n",
    "\n",
    "**Неявный QR алгоритм** (со сдвигами) вычисляет собственные значения (разложение Шура).\n",
    "Но мы не можем применить его напрямую к бидиагональной матрице, поскольку она может быть недиагонализуема в общем случае.\n",
    "\n",
    "Однако задачу вычисления сингулярного разложения можно свести к **симметричной задаче на собственные значения** двумя способами:\n",
    "\n",
    "1. Работать с трёхдиагональной матрицей\n",
    "\n",
    "$$T = B^* B$$ \n",
    "\n",
    "2. Работать с расширенной матрицей \n",
    "\n",
    "$$T = \\begin{bmatrix} 0 & B \\\\ B^* & 0 \\end{bmatrix}$$\n",
    "\n",
    "\n",
    "Случай 1 практически реализуем, если **не формировать матрицу $T$ явно**!\n",
    "\n",
    "Таким образом, задача вычисления сингулярных чисел может быть сведена к задаче вычисления собственных чисел симметричной трёхдиагональной матрицы."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Алгоритмы для симметричной задачи на собственные значения\n",
    "\n",
    "Рассмотрены:\n",
    "\n",
    "- QR алгоритм: \"золотой стандарт\" для вычисления собственных значений\n",
    "- Метод итерации Релея\n",
    "\n",
    "Далее:\n",
    "\n",
    "- Алгоритм на основе стратегии \"разделяй и властвуй\"\n",
    "- Метод бисекции\n",
    "- Метод Якоби"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Метод разделяй и властвуй\n",
    "\n",
    "Пусть у нас есть трёхдиагональная матрица и мы разделили её на блоки:\n",
    "\n",
    "\n",
    "$$T = \\begin{bmatrix} T'_1 & B \\\\ B^{\\top} & T'_2 \\end{bmatrix}$$\n",
    "\n",
    "Можем записать матрицу $T$ в виде\n",
    "\n",
    "$$T = \\begin{bmatrix} T_1 & 0 \\\\ 0 & T_2 \\end{bmatrix} + b_m v v^*$$\n",
    "\n",
    "где $vv^*$ – матрица ранга 1, $v = (0,\\dots,0,1,1,0,\\dots,0)^T$.\n",
    "\n",
    "Пусть мы уже разложили матрицы $T_1$ и $T_2$:\n",
    "\n",
    "$$T_1 = Q_1 \\Lambda_1 Q^*_1, \\quad T_2 = Q_2 \\Lambda_2 Q^*_2$$\n",
    "\n",
    "Тогда (проверьте!),\n",
    "\n",
    "$$\\begin{bmatrix} Q^*_1 & 0 \\\\ 0 & Q^*_2 \\end{bmatrix} T\\begin{bmatrix} Q_1 & 0 \\\\ 0 & Q_2 \\end{bmatrix} = D + \\rho u u^{*}, \\quad D = \\begin{bmatrix} \\Lambda_1 & 0 \\\\ 0 & \\Lambda_2\\end{bmatrix},$$\n",
    "\n",
    "то есть мы свели задачу к задаче вычисления собственных значений у матрицы вида\n",
    "\n",
    "\n",
    "<font color='red'> диагональная матрица плюс матрица малого ранга. </font>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Матрица вида диагональная матрица плюс матрица малого ранга\n",
    "\n",
    "- Собственные значения матрицы  вида\n",
    "\n",
    "$$D + \\rho u u^* $$\n",
    "\n",
    "вычислить не так просто!\n",
    "\n",
    "- Характеристический многочлен имеет вид\n",
    "\n",
    "$$\\det(D + \\rho uu^* - \\lambda I) = \\det(D - \\lambda I)\\det(I + \\rho (D - \\lambda I)^{-1} uu^*) = 0.$$\n",
    "\n",
    "Тогда (докажите!!) \n",
    "\n",
    "$$\\det(I + \\rho (D - \\lambda I)^{-1} uu^*) = 1 + \\rho \\sum_{i=1}^n \\frac{|u_i|^2}{d_i - \\lambda} = 0$$\n",
    "\n",
    "Подсказка: найдите $\\det(I + w u^*)$ используя два факта: \n",
    "- $\\text{det}(C) = \\prod_{i=1}^n\\lambda_i(C)$\n",
    "- $\\text{trace}(C) = \\sum_{i=1}^n \\lambda_i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Характеристическое уравнение\n",
    "\n",
    "$$1 + \\rho \\sum_{i=1}^n \\frac{|u_i|^2}{d_i - \\lambda} = 0$$\n",
    "\n",
    "Как найти корни?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEdCAYAAAALugwIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzq0lEQVR4nO3deXxU9b3/8dcnQNgJOwQQIi6AuLBEtNpaRFGp2vorilTRVlQqFa3W661a61KpF+t6rdILbhUUr1aRIm1RuGBVqhJQFIHIEhBIgCyE7CHb9/fHzFgICZwsZ858mM/z8ZhHyEwy521M3vOd7znne8Q5hzHGmPiREHQAY4wx0WXFb4wxccaK3xhj4owVvzHGxBkrfmOMiTNW/MYYE2es+I0xJs5ErfhFpKuIPCYim0WkXERyRGS5iHwvWhmMMcZAy2hsREQGAO8DHYAXgI1AEnAq0DcaGYwxxoREpfiBV8LbOtU5tytK2zTGGFMH34tfRM4Bvgvc6pzbJSKtgFbOuVK/t22MMeZQ0Zjj/0H443YReQcoA0pEZKOITIrC9o0xxhwgGlM9g8IfnwM2AT8FWgO/AuaKSCvn3Eu1v0lEpgBTANq3bz9y8ODBUYhqDqe8vJx169YxcOBAcnNzqaqqYsiQIUHHqldWVha7du1i2LBhrFmzhn79+tGrV6+gY9VrzZo1dOvWjezsbPr06UNycnLQkeqVn59PRkYGvXv3Zvfu3Zx88sm0bt066FimltWrV+c653oc8oBzztcbsBRwwBYg8YD7uwD5wC4g4XDPMXLkSGeC99lnnznAvf32227cuHHu9NNPDzrSYd11110uMTHRFRQUOMA9/vjjQUc6rLZt27o77rjDAe7BBx8MOs5hzZ071wHuvvvuc4DbtGlT0JFMHYBVro5OjcZUT1n442vOuYoDXnDygYVAb/79rsDEsP379wOoGdmVl5fTpk2boGN4VlFRQatWrYKO4UllZSWAmrzmYNEo/p3hj7vreCxyhE+XKOQwTVRREXrdTkxMDDiJN5qKv6amhurqajU/Wyt+3aJR/CvDH/vV8Vjkvuwo5DBNZCN+/0SK1IrfREM0in8BUARMEpEOkTtFJBm4DNjknNschRymiTQWv5askXdTWorUil8334/qcc7li8h/ALOAT0TkRSARmBr+OM3vDKZ5aCx+G/H7w4pft6icueucmy0iucB/Ag8BNcDHwFXOuRXRyGCarqSkBID27dsHnMSb/fv3qyl+bftPtL1QmYNFa8kGnHPzgfnR2p5pfqWloZOt27VrF3ASbzSN+LVN9WjLaw5myzIbz7SN+DUVv7YRdGVlJS1atEBEgo5iGsGK33hmI37/aJzqsdG+Xlb8xrOSkhISExNp2TJqM4RNYkf1+MeKXzcrfuNZSUmJmmke0DXiLysLneDetm3bgJN4Y8WvmxW/8ay0tFTNNA/oOqonMo1mxW+iwYrfeGYjfv9ERvxaXlgrKirU7I8wh7LiN56VlpZa8ftE21RPWVmZmqzmUFb8xrOSkhI1I1LnnKqdu1b8Jpqs+I1nmkb8VVVV1NTUqBnxa5vjt+LXzYrfeKZpxF9eXg6gpvi1zfFrmkYzh7LiN55pGvFHzjLu0KHDEb4yNthUj4kmK37jmaajeoqKigBdxS8iqvZJWPHrZcVvPCsqKlJTpMXFxYCe4i8tLaVNmzZq1r4pLy+34lfMit94UlVVRUlJCZ07dw46iifair+srEzN/D6E8tocv15W/MaTgoICAJKSkgJO4o3G4tc0gtaW1xzMit94Eil+G/H7o7S0VFWR2lSPblb8xpN9+/YBNuL3i6YRtHNOVV5zKCt+44m2EX/kqJ6OHTsGnMQbTXP8kWsv2xy/Xlb8xhMb8ftL0wha2zkH5lBW/MaTSPFrGfEXFxfTsmVLNStIajwr2opfLyt+44nGo3o6dOig5rj4oqIiOnXqFHQMT2zEr58Vv/EkMuLXUk6R4teisLBQzc82Uvw2x69X1ItfRNqJyFYRcSLyTLS3bxpn7969JCUl0aJFi6CjeKLpLGPQVfw21aNfECP+3wHdA9iuaYKcnBx69OgRdAzP8vPz6dKlS9AxPKmsrKS8vFxN8UcWwNOyT8IcKqrFLyIjgNuA+6O5XdN02op/7969dO3aNegYnkQOPdVS/NoOlTWHilrxi0gL4DlgMTA/Wts1zSM7O5uePXsGHcMzTSP+wsJCQE+RanuhMoeK5oj/dmAwMC2K2zTNRNuIPz8/X82IP1L8WopU2wuVOVRUil9EjgUeBH7nnNsWjW2a5lNTU0Nubu4hI37nXECJDq+6upqCggJ1I34txV/XVE+s/i6YukVrxP8nYCvwhNdvEJEpIrJKRFbl5OT4l8wc0b59+6iqqjpoxB/Lx8dHDj21Eb8/DrzITSz/Hpj6+V78IjIJuAC4yTlX6fX7nHOznXOpzrlUTVMMR6Ps7GwANXP8+fn5ADbi90nkUNmEBDsNSKuWfj65iLQmNMr/O7BbRI4PP9Q3/DEpfF+uc26fn1lM4+3cuROAfv36BZzEm7179wI24vdLYWGhze8r5/dLdlugB3AxsOmA2/vhxyeFP7/B5xymCXbs2AHoK34tI/68vDxAzwtVUVGRFb9yvo74gRLgijru7wHMJHRo5wvAlz7nME0QGfH37dv3CF8ZG/bs2QNAr169Ak7iTW5uLu3ataNt27bU1NQEHeeIrPj187X4w3P6b9a+X0RSwv/c4pw75HETW3bu3EnPnj1p3bp10FE82b17N6Cn+PPy8ujWrVvQMTzTtKCcqZvtnTFHtGPHDjXTPBAa8bdv317NWj25ubl0765nFROb49cvkOJ3zm1zzolzzk7mUiAjI4OUlJSgY3i2e/duevfuHXQMz7SN+Pft26fmugymbjbiN4dVVVXFli1bOPHEE4OO4tmePXvUTPOAvhG/thcqcygrfnNY33zzDVVVVaqK30b8/qmoqKC4uFjNEUimblb85rA2btwIoKr4MzMzSU5ODjqGJ1VVVeTn56sp/sihp1rymrpZ8ZvD0lb8BQUFFBQUMGDAgKCjeBI59LRPnz4BJ/Emco6EFb9uVvzmsDZu3Ejnzp3VzEFHTjbr379/wEm8ycrKAvQUv7aTzUzdrPjNYX311VcMHjxYzWJc33zzDYCaEb/W4rcRv25W/KZeNTU1fPbZZ6SmpgYdxbPt27cDekb8u3btAqz4TXRZ8Zt6bdy4keLiYkaOHBl0FM8yMjJo3bq1mqN6srKySEhIULPyaW5uLmDFr50Vv6nXqlWrAFSN+NPT0xk0aJCaJYMzMzPp1asXLVq0CDqKJ7t27aJjx460b98+6CimCXT8dZhApKWl0bZtWwYPHhx0FM/S09PrzRuLV4natm2bqrOid+3apWZaytTPit/Ua/ny5Zx11lm0bOn3Iq7No7y8nIyMjEOKP5Z3TGdkZHDssccGHcOzXbt2qTlHwtTPit/UKTs7m7Vr13LeeecFHcWzzZs3U1NTo+YdSlVVFTt27LDiN1FnxW/qtHz5cgBVxZ+eng6gpvh37NhBdXW1muJ3zlnxHyWs+E2dFi1aRNeuXRkxYkTQUTxbs2YNLVq0YNCgQUFH8WTz5s0AHHfccQEn8aawsJDS0lIr/qOAFb85xP79+1m4cCGXXXaZmvl9CO2MPvnkk2nXrl3QUTxZt24dACeddFLASbzZtm0bgKqd0aZuVvzmEO+99x6FhYVcfvnlQUfxzDlHWloap59+etBRPFu3bh3du3dXcwy/Ff/Rw4rfHOKFF16gR48equb3t2zZQn5+PqNGjQo6imfr1q1j6NChQcfwzIr/6GHFbw6SmZnJokWLmDx5MomJiUHH8ezjjz8GUFP8zjl1xb9161Y6dOhgZ+0eBaz4zUGeeuopnHNMmTIl6CgNsnTpUrp168Ypp5wSdBRPdu7cSWFhoari37JlC8cee2xMnxdhvLHiN9/Kzs7m2Wef5eqrr2bgwIFBx/HMOcfSpUs577zz1CzVsHLlSgBV6yCtX79ezY5oc3g6/kpMVDz00EPs37+fe++9N+goDbJhwwaysrI4//zzg47i2ccff0zr1q0ZPnx40FE8KS0tZevWrVb8RwkrfgPA6tWrmTlzJjfddJOaq21FLFiwAIBx48YFG6QBPvnkE0aMGKFmP0p6ejrOOSv+o4QVv6GsrIzJkyfTs2dPHn744aDjNNibb77JWWedRb9+/YKO4klFRQWrV6/mO9/5TtBRPFu/fj2Aqn0Spn6+F7+InCgivxORT0QkR0SKRGSNiPxGRGxt14A557jtttv48ssveeGFF0hKSgo6UoN8/fXXfP7554wfPz7oKJ6tWLGC8vJyzjnnnKCjeLZu3TpatmzJ8ccfH3QU0wyiMeKfDNwObAF+B9wJfA1MB/4lIm2jkMHU49FHH2X27Nncdddd/OAHPwg6ToPNnj2bli1bcvXVVwcdxbPFixfTqlUrxowZE3QUz9LS0jj11FNp1apV0FFMM4hG8b8J9HPOXe2c+6Nz7n+cc1cCvwdOBa6PQgZTh5kzZ/LrX/+aiRMn8vvf/z7oOA1WVlbGn//8Z3784x/Tq1evoON4tnjxYs4++2w6duwYdBRPqqurWblyJWeeeWbQUUwz8b34nXOrnHMFdTz0evjjyX5nMAdzznHvvfdy8803c8kll/Dyyy+rOQzyQPPmzWPv3r38/Oc/DzqKZxkZGXz55Zeq3l1t2LCBoqIiK/6jSJArcEX2xO0JMEPcycnJYfLkySxatIgbb7yRmTNnNnohtiCvaFVZWcnDDz/MyJEjOffccwPL0VCvvx4a70yYMOGwXxdLVwuLnBV9uOKPpbzmyAIpfhFpAdwHVAHz6vmaKcAUgP79+0cv3FHKOccbb7zB7bffTl5eHv/93//NLbfc0uizMIM+e3POnDlkZGTwzjvvBJ7FK+cc8+bN46yzzmLAgAF1fk0s/rcsW7aMXr161bljNxbzmiML6v39U8CZwH3Oua/r+gLn3GznXKpzLrVHjx5RDXe0WblyJeeffz4TJ06kd+/efPrpp9x6661q/2gLCgr4zW9+wxlnnMHFF18cdBzPVqxYwVdffcVPf/rToKN4Vl1dzZIlS7jgggvU/r6YQ0W9+EXkIWAaMNs591/R3n68cM6xbNkyLrzwQs444wzWrFnDzJkzSUtLY9iwYUHHa5L777+f7OxsnnnmGVVl9Mwzz9C5c2dVRyB99tln5OXlcdFFFwUdxTSjqE71iMgDwL3AS8BN0dx2vNi9ezdz5szh+eefZ9OmTfTo0YNHHnmEqVOnqjmK5HCWL1/O008/zdSpU0lNTQ06jmfbt2/nrbfe4tZbb6V9ez2nryxcuJCEhATGjh0bdBTTjKJW/CJyP3A/MAe4wdneoGazdetWFixYwPz581mxYgXOOb73ve/x29/+lssvv5y2bY+OUyXy8/O59tprOeGEE/jDH/4QdJwGmT59OgkJCfzyl78MOopnzjlee+01xowZg023Hl2iUvwich/wADAXuM45VxON7R6tdu/ezfvvv8+yZctYtmwZW7ZsAeC0007jgQceYMKECWouOO5VVVUVEyZMYM+ePaxYsULVqHnz5s28+OKLTJ06VdWBCqtWrWLLli3cfffdQUcxzcz34heRm4EHge3AUuCqWvOye5xzS/zOoVV+fj6rV69m1apVpKWlkZaWxo4dOwBISkpi9OjR3HLLLVxyySVqLtrdGLfffjtLly7lxRdfVHV5Recc06ZNo23bttxzzz1Bx2mQWbNm0a5dO1XLYRhvojHij/yV9gderuPxfwJxXfzOOfbs2UN6ejobNmxg/fr1395279797dcdd9xxnH322aSmpnLOOecwYsQIWrRoEWBy/znnuOeee3jmmWe44447uO6664KO1CBvvPEG7777Lk899RTJyclBx/EsLy+PV199lWuvvZbOnTsHHcc0M9+L3zn3M+Bnfm8nllVVVZGdnU1WVhZbt25l27ZtbNu27aB/l5WVffv1HTt2ZMiQIVx00UUMGTKEYcOGkZqaSteuXQP8r4i+mpoabr31Vp599lluuukmdfP6u3bt4pZbbmHEiBHcfPPNQcdpkNmzZ1NeXs60adOCjmJ8EOSZu6pVVVWxd+9e8vLyyMvLIycnh127dh10y8rKYteuXeTk5FBTc/BujS5dupCSksKQIUMYN24cKSkpDBo0iJNOOom+ffuqOkzRD1VVVVx//fXMmTOHO++8k0ceeUTVz6S6uppJkyZRXFzMnDlzGn12dBAKCgp47LHHGDdunJpLWZqG0fPb6IOKigoKCwspLCykoKDgoH8XFBQcVOwH3nJzcykoqGv5IUhISKBXr14kJyfTt29fUlNTSU5OJjk5mT59+pCSkkJKSoq65Y+jaf/+/Vx11VXMnz+fhx56iN/85jeqSh/g3nvvZdmyZTz//PPq1rB//PHH2bt3L9OnTw86ivHJUV38f/zjH0lLS/u20GsXfHl5+RGfo2PHjnTr1o3u3bvTrVs3jj/+eLp163bQrXv37nTv3p3k5GR69Ohx1M+7+6m0tJTx48ezePFinnzySW677bagIzXYc889x4wZM/j5z3/O5MmTg47TIJs2beIPf/gDEyZMYMSIEUHHMT45qot/5cqVfPTRR3Tq1ImkpCR69+7NiSeeSFJSEp06dfr2duDnB/67S5cuai6NdzQoLCzk0ksv5cMPP+T555/n+uv1rdj9t7/9jalTp3LRRRepO7PYOcdNN91E69atefLJJ4OOY3x0VBf/3Llzg45gPMrLy2PcuHF8/vnnvPbaa1x55ZVBR2qwt99+myuvvJLTTjuN119/XdW8PoSuz7Bs2TJmzpxJnz59go5jfKRvEXZz1Nm9ezejR4/myy+/ZP78+SpL/3//93+54oorGDFiBEuXLqVTp05BR2qQVatW8atf/Yof/OAHqq5vYBrHit8EaufOnZxzzjls3bqVv/3tb1x66aVBR2qwl156iauuuoqzzz6bJUuW0KVLl6AjNUhWVta3VzGbM2eOyovymIbR9V7UHFV27NjBueeeS3Z2Nu+99x5nnXVW0JEabObMmdx8882MHTuWBQsW0K5du6AjNUhRUREXX3wx+fn5fPDBB3Tr1i3oSCYK7KXdBGLnzp2MHj2anJwclixZorL0n3jiCW6++WYuvfRSFi5cqK70CwsLueSSS1i7di1vvvkmw4cPDzqSiRIrfhN1+/bt46KLLiI3N5f33nuPM844I+hIDTZ9+nTuuOMOrrjiCt58803atGkTdKQGyc3NZcyYMfzrX//ilVde4cILLww6kokim+oxUVVRUcGPf/xjNm7cyOLFi9WVvnOO++67j+nTp3PNNdfw4osvqjt6JzMzk7Fjx367nLemq5iZ5qHrN9aod+edd7J8+XLmzp3LmDFjgo7TIM457r33Xh5++GFuuOEGZs2apW5H6ObNmxk7dix5eXm8++67nHPOOUFHMgGw4jdR89Zbb/H0009z2223MWnSpKDjNEhkldAZM2YwZcoU/vSnP6kr/bVr1zJ27Fiqq6tZvnw5I0eODDqSCYiu31yjVnZ2NlOmTGHUqFE88sgjQcdpsBkzZjBjxgxuuukmlaX/ySef8P3vf5+WLVvywQcfWOkr4ZzDj4sV2ojfRMXtt99OUVERf/7zn9UtgzF37lzuuecerr76ap599ll1pb906VIuu+wykpOTWbJkCSkpKUFHigs1NTUUFRWxb9++I94KCgooKir69lZYWEhRURHFxcVs2LCBE044oVmzWfEb3y1btox58+Zx//33M2TIkKDjNMgHH3zA5MmTGTNmDC+++KK60l+0aBHjx49n0KBBvPfee/Tu3TvoSCo55ygqKiI7O5ucnJxvb7m5uQd9HrkvUuZHGq137NiRzp0706lTJzp27EinTp3o27fvt59H7mtuVvzGV8457r77bo455hjuuuuuoOM0SHZ2NhMnTmTgwIHMnz9f3TuVN998k5/85CcMGzaMd999N+4u5OOVc47du3ezbds2MjMzD7plZWV9++/S0tI6v79t27b06NHj29vgwYPp0qULnTt3PuytU6dOgR0RZsVvfPXOO++wcuVKXnjhBVXHutfU1DBp0iTy8/NZvHixuusnvPbaa0yaNIkzzzyTv//97+ryN7eamhq2bdvGhg0b2LJlC1u2bCEjI4OMjAy2bt160BXwAFq3bk2fPn3o27cvI0aM4JJLLqFPnz707NnzoJLv0aOHuhP3wIrf+Oyxxx7j2GOP5dprrw06SoP8z//8D0uWLGHWrFmceuqpQcdpkMWLF3Pttdfyve99j0WLFtGhQ4egI0VVfn4+aWlpfPHFF3z11VesW7eODRs2HDRi79ChAwMHDmTQoEGMGzeOgQMHkpKSQr9+/ejTpw/dunVTtaR2Q1nxG9988cUXfPjhhzz++OOqTnLatWsXd999N+effz433nhj0HEa5NNPP2X8+PGcfPLJLFy48Kgv/ZqaGr766ivef/99Vq5cycqVK9m0adO3jycnJzN06FCmTJnC0KFDOemkkzjhhBPo3r37UV3sR6Lnr9Go8+yzz9KuXTuuu+66oKM0yJ133sn+/fuZOXOmqnLIysrihz/8Ib169eIf//iHuqWhvdqzZw/vvPMOS5cuZdmyZeTk5AChkj/jjDO47rrrGDVqFMOHD7f9GvWw4je+qKio4I033uDyyy9XtUzx2rVrmTdvHr/+9a+b/RA6P1VVVTFx4kSKi4tZtmzZUXf0zu7du3njjTd46623+PDDD3HOkZyczIUXXsh5553HmDFj6N+/f9Ax1fC9+EUkAfgl8HMgBcgB3gDuc86V+L19E4wlS5ZQUFDAhAkTgo7SIPfffz8dO3bkzjvvDDpKgzzwwAN8+OGHzJ07V93F3evjnGPp0qXMmjWLv/71r1RVVTF06FB++9vfMn78eE455RRV78hiSTRG/E8CtwJvA48DQ8KfDxeR851zNVHIYKLsL3/5C0lJSYwdOzboKJ5t2rSJt99+m/vuu0/VFMEXX3zBjBkz+OlPf6puKYy61NTUMH/+fKZPn84XX3xBt27duPXWW7n++us56aSTgo53VPC1+EVkKHALMN85N/6A+7cCTwMTgXl+ZjDRV1VVxV//+ld+9KMfqTr2febMmbRq1YqpU6cGHcWz6upqbrzxRrp168YTTzwRdJwmW7FiBdOmTWPNmjWceOKJvPTSS0ycOFHVocAa+H0a4k8AAZ6qdf9zQCmgf3hiDvH555+zb98+xo0bF3QUz0pKSnjppZe4/PLLVc2Pv/rqq6SlpfHkk0+qepdSW2FhIZMnT+a73/0ueXl5vPLKK6xfv56f/exnVvo+ED8WAPr2yUXeBc4H2jnn9td6bAVwonOux5GeJzU11a1atarB23/wnXWszyps8PeZI1v71Voq9lfUudjXjh07yMjI4KyzvkOrVsGP+Kurq/noo48YOHAgxxxzTJ1fk52dzYYNGxg27DSSkjpHN2At//znP0lJSWHAgAGH/TrnHCtXrqRlq5aMHBHMomuRn9uoUafTtm3jTmQqLi5i3fr1lJeV079/fwYM6E9CQotmTqrTSX06cf+ljd9nIyKrnXOpte/3e8TfB8itXfphmUB3EamzGURkioisEpFVkcO1TOwQ6t+plp+fT7v27WKi9L3KyckhMTEx8NJviN27d1NeXs6xKccGHaXR8vLy+PzzNbgax7Bhwzj22GOt9KMhsuynHzdgC7C9nsfmAA7ofKTnGTlypDOx5ZJLLnEjRow45P7q6mrXoUMHN3Xq1ABS1a2oqMgB7tFHH6338TZt2rhp06ZFOdmhampqHOAeeOCBI37d0KFD3fDhw11NTU2U0h1q3rx5DnDp6ekN/t6XX37ZJSQkuNNPP93t2bPHh3QGWOXq6FS/j+opBXrW81ibA77GHCU2btxIcXExp59+etBRPPvHP/5BeXk5l19+edBRPFu+fDnr1q3jpZdeUnlI48KFC7nuuusYM2YMCxYsoH379kFHiit+T/VkEZrOaV3HY30JTQNV+JzBRNHq1asBVF3oY8mSJSQlJfHd73436CiezZo1i+7duzNx4sSgozTY2rVrufLKK0lNTeXtt9+20g+A38WfFt7GqAPvFJE2wDCg4XtsTUxbvXo1bdq0UXW89f/93/8xevRoWrTQMbdcUlLCO++8w4QJE9Qd8VJSUsKVV15J586d42ItoVjld/G/Tmge/7Za998ItANe9Xn7JspWr17NaaedpmZRtm3btpGRkaHqwu+LFi2irKxM3VnRAHfddRfp6enMnTuXXr16BR0nbvla/M65tcCzwI9FZL6I3CAijwNPAP/ETt466qSnp3PyyScHHcOzjz76CIDRo0cHG6QBXn/9dZKTk1VNTUHoDOOZM2fyi1/8gvPPPz/oOHEtGsOy24BtwBTgYiAX+COhtXpsuYajyL59+8jOzubEE08MOopnq1atol27dmqmpiorK3nvvfe45ppr1ExNRfzqV7+ia9euPPTQQ0FHiXu+F79zrprQGj2P+70tE6zIOujain/48OFqpqY+/fRTSkpKuOCCC4KO0iAff/wxy5Yt48knn1S1WuvRSteVo01M27hxI6Cn+Kuqqvj8889JTT3kxMaYtXTpUhISElRNTQHMmDGDrl27qruwzdHKit80m40bN5KQkMBxxx0XdBRPvv76a0pLS1Uderp8+XJGjBihatS8detWFi5cyLRp0+zQzRhhxW+azcaNGxkwYACtW9d12kbsWb9+PQCnnHJKwEm8qaqqYtWqVXznO98JOkqDvPLKKwBcf/31AScxEVb8ptls376dlJSUoGN49vXXXwOoudJW5ILho0aNOvIXxwjnHHPnzmX06NF2hawYYsVvmk1mZiZ9+/YNOoZnX3/9Ncccc4ya6Ye0tDQAVcthrF69mk2bNnHNNdcEHcUcwIrfNAvnHFlZWeqKf9CgQUHH8CwtLY1OnTqpeYcCoZPNEhIS+NGPfhR0FHMAK37TLHJzc6msrKRPnz5BR/HEOaeu+NetW8epp55KQoKeP9vFixczatQounXrFnQUcwA9v0EmpmVmZgKoGfFnZ2dTWFioqvjT09MZPHhw0DE8y8vLY+XKlVx00UVBRzG1WPGbZqGt+L/55hsANTuj8/LyyMnJUVX8K1aswDmnah2keGHFb5pFVlYWgJqpnsgLVb9+/QJO4k3kCCRNxf/JJ5/QsmVLVSfIxQsrftMsIkWanJwccBJvdu7cCeh5h5Keng7oK/5hw4bRtm3boKOYWqz4TbPIysqiZ8+etGrVKugonmRmZpKYmEj37t2DjuJJeno6iYmJaqamqqurWblyJWeeeWbQUUwdrPhNs8jNzaVHjx5Bx/AsMzOTPn36qDlCZvv27fTv31/NipybN2+mpKTEpnlilI7fehPz8vLyVB2yt3PnTjXTPAA7duzgmGOOCTqGZ+vWrQNg6NChAScxdbHiN81i7969dO3aNegYnmVmZqrZsQuhFypNeSPrIGnaJxFPrPhNs9A04nfOqRrxV1dXk5mZqWrEv379elJSUuyaujHKit80mXNOVfEXFhZSVlam5tDTPXv2UF1drW7Er+WqZvHIit80WWlpKRUVFWqKPzc3F0DNET07duwAUDPir6mpIT09nSFDhgQdxdTDit80WV5eHoCaOf5IXi3FHznnQMuIf8+ePezfv5+BAwcGHcXUw4rfNFmkSLWN+LXk1Tbi3759O4Ctvx/DrPhNkxUUFADQuXPnYIN4pG3Ev2fPHlq1aqXmHZUVf+yz4jdNVlhYCECnTp0CTuKNthF/bm4u3bt3R0SCjuKJFX/s87X4RaSviNwtIv8UkV0iUiIi60TkURHR8VdnjqioqAiAjh07BpzEm7y8PFq0aEFSUlLQUTyJFL8W27dvp2PHjmp+vvHI7xH/pcADQB7wKHAb8K/wx89FpLfP2zdRoK34c3Nz6dq1q5rlGnJyctQV/zHHHKPmHUo8aunz838IDHDO7T7gvudE5FPgOeA/wjejmLbi37t3r5ppHgi9UJ122mlBx/Assq6QiV2+Dnmcc+tqlX7E6+GPJ/u5fRMdRUVFiIiai5YXFhaqmobQNtWjbV2heBTUe93IAcl7Atq+aUaFhYV07NhRzVv7SF4Nqqur2bt3r5rir66uJjc3l169egUdxRxGUMX/YPjjy/V9gYhMEZFVIrIqJycnSrFMQzjngNCIX0uRQqj4Y/0IpMjPNj8/H+dczE9NRfLu3bsX5xw9e/YMOJE5HE9z/CLSmdAOWa+eds7tree57gCuAGY755bV9wTOudnAbIDU1FTXgG2bKDhwdG/F33xqv2vat28fAF26dAkgzZHVzpudnQ2g6toM8cjrzt3OwP0NeN5XgEOKX0RuIHR0z9+AaQ14PhPDioqKYrZI6xLLxV9b5OQ4LXkj785txB/bPBW/c24b0KQJXBGZTGgE/x4w3jlX2ZTnM7FD05y5c07VC1Xk5DgtO6MjxW8j/tgWlTl+EbmO0OGbS4HLnHP7o7FdEx2apnpKS0upqalRU/yREb+W4o9M9diIP7b5Xvwi8jPgeWA58CPnXLnf2zTRVVpaSrt27YKO4UlkBK3lhUrrVE+s74yOd76ewCUiPwReAAoJHbs/vtbOoGLn3AI/Mxj/lZeX07Zt26BjeKJtXSFtI/78/HySkpJo2dLvc0NNU/j9f2cEoXcVnQkfoVPLN8ACnzMYn5WVlVnx+0Rb3oKCAjUvUvHM7zN3H3DOyWFuKX5u30RHWVkZbdq0CTqGJxqnetq0aUNiYmLQUTzZt2+fmuW545mOVapMzHLOqZrqKS0tBVCzT0LbCFpb3nhlxW+apKKiAuecmuIvKysDrPj9YiN+Haz4TZNEilTLVE8kr5YXquLiYjXTUqDvhSpeWfGbJikvDx2dq6VII1M9mvJqWfUUbMSvhRW/aRKtI34tUz2lpaVqXqScczbiV8KK3zSJtqkTbXk1nRxXXFxMTU2NjfgVsOI3TaJtqqesrIyEhARatWoVdBRPNBW/tpPN4pkVv2kSbVM9kakTLReN0VT8xcXFgJ5zJOKZFb9pEm1TJ2VlZWqKFHQVf0lJCaBn/0k8s+I3TaKp+J1zqpaXAF3FHzliStNRSPHKit80SWSOP5aneg6c1tFU/JWVlVRVVakp/siI34o/9lnxmybRNOIHXYdHajz0FPTkjWdW/KZJtBW/pjl+bSeb2YhfDyt+0yQapnoOpGmqR9sIWlveeGbFb5qkoqICQM2ywVb8/rERvx5W/KZJKisrAVSdEGXF7w87nFMPK37TJNqK364d4J/S0lISExPtsosKWPGbJqmsrCQhIYGEBB2/ShUVFaqmpUBP8ZeUlKjJGu90/LWamFVZWalmtA+hvFqKP7LjvHXr1gEn8UbbEtLxzIrfNIm24q+oqFCTN7LjXEvxl5SUWPErYcVvmkRTkYKuEb+2I6Y0LS8R76z4TZNoKlLQ9UJlxW/8EtXiF5EEEflYRJyILIrmto0/NE311NTUUF1draZItRV/eXm5mhP54l20R/y/AE6J8jaNjzQVv7ZDT7UVv6YjpuJd1IpfRPoBDwP3RWubxn+ail9jkYKevPv371ezIzreRXPE/wyQAfx3FLdpfKap+G3E7y8b8esRlVPsRORy4IfAWc65ai2XvTNHpqn4NRYpoOZMWBvx6+H7iF9EkoCngVnOuU/83p6JHuecquLXOOJPTExUcX1g55yN+BXxNJQQkc7AbQ143qedc3vD//4DoReYuxsSTESmAFMA+vfv35BvNVEQKSNNxa9xxB/rWQ98UbIRvx5e30N2Bu5vwPO+AuwVke8CNwLXOOf2NSSYc242MBsgNTXVNeR7TfRoKn6tI34ttOWNZ56K3zm3DWjM+81ngS+AT0Xk+FqPtQvft885l9uI5zYxoLKyUs1p+jbi95eN+PXwe6/RACAJ2FTHY+eG738WmOZzDuMTjSN+LWWqrfi15Y1nfhf/tUBdvwl/AVYDM4DNPmcwPtJU/JERv6a8Woq0qqqKmpoaG/Er4WvxO+cW1nV/eIfQbufcm35u3/hP09o3NuL3j7aVROOdLdJmmkTTIm024vfP/v37AT0vqvEukDNDnHOxf2Cy8UTTVI+N+P1jI35dbMRvmkRT8duI3z824tfFit80iabitxG/f2zEr4sVv2kSTcVvI37/2IhfFyt+0ySait9G/P6xEb8uVvymSTQVv434/WMjfl2s+E2TaCp+G/H7x0b8uljxm0arrq7GOaem+G3E7x8b8etixW8aTdtqlzbi94+N+HWx4jeNpq34tY34KysrVV19C/S8qMY7K37TaDU1NQC0aNEi4CTeVFdXA3ryVlVVqcmqbcnreGfFbxotUqQJCTp+jTTm1TLir6qqAvRcHzje6fgLMDFJ44g/ISFBxTVsIZRXy8/Wil8XK37TaJHi1zKCrqmpUVOkoKv4tU2jxTsdf7EmJmmcOtGSFXQWv434ddDzV2BijrapHuecmqygq/gjUz1a8sY7K37TaNpG/KCnmLS9qNpUjy56/mJNzNFWTqAnq7YitZ27uljxm0bTtnMX9GTVVvza8sY7HX8FJibZVI9/tBWpjfh10fMXa2KOTfX4R1vxa8sb76z4TaPZiN8/2orUjurRRc9frIk5Gkf8Wl6ktBV/dXU1IqLm5xvv7P+SaTSNO3c1FSnoyatpQTkTpeIXkXYicp+IrBORMhHZKyIfi8j/i8b2jT+s+P2jrfg1LShnwPf/UyLSBfg/4ATgJeAJoD0wBBjg9/aN/7SUE+jJqq34bcSvSzReop8GjgfOcM6tj8L2TJRpGvFryaqt+G3Er4uvfwUikgJcBTznnFsvIi1EpIOf2zTRp6WcQE9WbcVvI35d/B7+XBTexnoRmQuUAkUislNEbvd52yZKtIyiQU+Rait+TQvKGRDnnH9PLvIkcBuQA+QCjwIVwE3Ad4HfOefur+d7pwBTwp8OAr5uZIzu4W1roSmvZfWPpryasoKuvE3NOsA516P2nZ6KX0Q6Eypwr552zu0VkeeB64EC4DjnXF74+VoB64F+QB/nXH4DnrtBRGSVcy7Vr+dvbpryWlb/aMqrKSvoyutXVq97YzoDdY7M6/EKsBcoC3++KFL6AM65ShGZB9wHnAn8owHPbYwxpgk8Fb9zbhvQmAuV7gx/3F3HY7vCH7s04nmNMcY0kt975VaGP/ar47HIfdk+Z5jt8/M3N015Lat/NOXVlBV05fUlq987d1sAW4AewInOuczw/e2BDUBHoJ9zrsS3EMYYYw7i6xkXzrlqEZkKvAN8LCIzCR3Vcx1wDHC9lb4xxkSXryP+bzci8n1CO4dPB1oAnwMznHPv+L5xY4wxB4lK8RtjjIkdek65bCYi8nMReVVE0kWkWkQCfeUTkQQRuT2cp1xEdojI4+H9IDFFRO4Wkb+ISIaIOBHZFnSm+ojIiSLyOxH5RERyRKRIRNaIyG9i9Gc7KPx7uUFECkSkNPw78YSIJAed70jCK/BuDf9ePBN0ntrCueq6FQedrS4i0lVEHhORzeFeyBGR5SLyveZ4/nhcVeluoBuh6ab21H3EUTQ9CdwKvA08TmjV0luB4SJyvnOuJshwtTxM6PyMzwid2xHLJgM3AwuBV4FK4FxgOjBBRM50zpUd5vujrR+QTOj3YCdQBZxC6Oz1iSIyzDnn9xFwTfE7QmeZxrIPOfQomcogghyOiAwA3gc6AC8AG4Ek4FSgb7NsxDkXVzcgBUgI/3tR6EcQWJahQA3wVq37bwEccFXQP69auQYe8O+vgG1BZzpM1lQgqY77p4d/ttOCzujxv+OKcN7/DDrLYTKOIPRC9atw1meCzlRHRgf8OegcHrN+COwAkv3aRtxN9TjntrnYGUX/hNCJcU/Vuv85QgvaTYp2oMNxzmUEncEr59wq51xBHQ+9Hv54cjTzNME34Y8xeaJj+JDt54DFwPyA4xyRiCTG8grBInIOoXXM/uCc2yUirUSkXXNvJ+6KP8acTmjEv/LAO51z5cCa8OOmeUWm9vYEmqIeItJGRLqLSD8RuQCYFX7o70HmOozbgcHAtKCDeHA5/14hOFtE/igiSUGHquUH4Y/bReQdQsvelIjIRhFptoGgFX+w+gC5zrn9dTyWCXQXkcQoZzpqhUen9xGalpgXcJz63EBoNdsdwLuE9qVMcs59GGSouojIscCDhFbZ3RZwnCNZCTxAqPx/Ciwj9GL1YYy9AxgU/vgc0JVQ1usJnf80V0Sua46NqNy529jVQv1J0yTtgLpKH6D8gK+piE6co95ThBYFvMc519hlvv22AEgntGNvOPBDQme+x6I/AVsJXU41pjnnzqh11xwR+RL4PfDL8MdY0DH8sQg41zlXASAibwMZwMMi8nJTp6tVFj+NXy001pQCPet5rM0BX2OaSEQeIjTCm+2c+6+g89THObeTfy9uuEBE3gLSRKRtLOUOTztcAJzjnIu5I2M8epRQj1xM7BR/5Eiz1yKlD+CcyxeRhcC1hN4VbGjKRlRO9YR30EoDbpuDzlyPLELTOa3reKwvoWkgG+03kYg8ANwLvEToIkBqOOe+JHTo8S+CzhIR/n19gtB+h90icryIHA8MCH9JUvi+zkFl9CL8gpVFbB2GGpUVjVUW/1EkjdD/g1EH3ikibYBhwKoAMh1VROR+QqO6OcANLny8nDJtCc33xoq2hKafLgY2HXB7P/z4pPDnNwQRzqvw31k/YmtHf1RWNLbiD9brhI4vvq3W/TcSmtt/NdqBjiYich+hHXpzgeti6DDeQ4hI73ruP5fQoaefRDfRYZUQOr+g9i3yrmRx+POFgaSrRUS61fPQQ4Smu2NpzbAFhOb3Jx240zl89vZlwKbmmMHQOsffaCJyKXBa+NPjw/fdG/58n3MuaqebO+fWisizwDQRmU/orXPkzN1/EmNHnojINfz77XwPIPGAn903zrm5wSQ7lIjcTOiIk+3AUuAqkYOuJbTHObckiGz1+FP4j3sZoWP32wAjgYmEiuCOALMdJDxF8mbt+0UkJfzPLc65Qx4P0L0iciawnNDvQwdCh02eC3wK/DHAbAcJz+X/B6HDeD8RkReBRGBq+GPzHDYb9Flq0b4BfyY0yq7rti2APC0I/VF/TegIn0xC86cdgv5Z1ZH1/cP87N4POl8D/j/HYt4JwN8IHcZZTmgnXzqhUuofdD6P/w0pxOCZu8CPCB0amxn+2ZYQOk/mHqBN0PnqyfxjQu/ySgi98L8HnN1cz2+rcxpjTJyxOX5jjIkzVvzGGBNnrPiNMSbOWPEbY0ycseI3xpg4Y8VvjDFxxorfGGPijBW/MY0gIv8Skf3hi7mnBJ3HmIaw4jemcR4ntPDbGcB/BJzFmAaxM3eNaSQRaQnkA185574TdB5jvLIRvzGN5JyrAr4CTpZaK8AZE8us+I1ppHDZJxJa7TEl2DTGeGfFb0zjTQVGhP99SpBBjGkIK35jGkFE+gD/xb8vkWfFb9Sw4jemcZ4BWgHjw59b8Rs14u4KXMY0lYhcBvw/4D+dc/8SkWxCl0c0RgU7nNOYBhCRTsB6QlM8ZzjnqkVkCTAaaO+cqwgynzFe2FSPMQ3zX0Av4AbnXHX4vi8JvXseHFgqYxrAit8Yj8IX7L4JeMw5t+aAh74Mf7R5fqOCFb8xHohIK+A5YAvwYK2HrfiNKrZz1xhv/hMYCpzrnCuv9dh6oAorfqOE7dw1xpg4Y1M9xhgTZ6z4jTEmzljxG2NMnLHiN8aYOGPFb4wxccaK3xhj4owVvzHGxBkrfmOMiTNW/MYYE2f+P7OseJdr7RWMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "lm = [1, 2, 3, 4]\n",
    "M = len(lm)\n",
    "D = np.array(lm)\n",
    "a = np.min(lm)\n",
    "b = np.max(lm)\n",
    "t = np.linspace(-1, 6, 1000)\n",
    "u = 0.5 * np.ones(M)\n",
    "rho = 1\n",
    "def fun(lam):\n",
    "    return 1 + rho * np.sum(u**2/(D - lam))\n",
    "res = [fun(lam) for lam in t]\n",
    "plt.plot(t, res, 'k')\n",
    "plt.plot(t, np.zeros_like(t))\n",
    "plt.ylim([-6, 6])\n",
    "plt.xticks(fontsize=18)\n",
    "plt.yticks(fontsize=18)\n",
    "_ = plt.xlabel(\"$\\lambda$\", fontsize=18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- Эта функция имеет только один корень на каждом отрезке $[d_i, d_{i+1}]$\n",
    "\n",
    "- Мы показали справедливость теоремы Коши о чередовании (что происходит с собственными числами после внесения возмущения ранга 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Как найти корни\n",
    "\n",
    "- Метод Ньютона не сработает (изобразите касательные к графику).\n",
    "\n",
    "- Заметим, что метод Ньютона – это по сути линейная аппроксимация функции $f(\\lambda)$ на каждой итерации.\n",
    "\n",
    "- Лучше аппроксимировать с помощью **гиперболы** вида:\n",
    "\n",
    "$$f(\\lambda) \\approx c_0 + \\frac{c_1}{d_i - \\lambda} + \\frac{c_2}{d_{i+1} - \\lambda}.$$\n",
    "\n",
    "- Для вычисления коэффициентов нам нужно вычислить $f(\\lambda)$ и $f'(\\lambda)$  в некоторых точках.\n",
    "\n",
    "- После чего получить апроксимацию из решения **квадратного уравнения**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Важный недостаток метода\n",
    "\n",
    "- Устойчивость: этот метод игнорировали долгое время из-за неустойчивого вычисления собственных векторов.\n",
    "\n",
    "- Нам нужно вычислить собственные векторы матрицы $D + \\rho uu^*$.\n",
    "\n",
    "- Точное выражение для собственных векторов $v_i$, для которых найдены собственные значения $\\alpha_i$:\n",
    "\n",
    "$$(D + \\rho uu^*)v_i = \\alpha_i v_i \\qquad (D - \\alpha_i I)v_i + \\rho (u^*v_i) u  = 0$$\n",
    "\n",
    "То есть $v_i \\in \\mathrm{span}(\\{(D - \\alpha_i I)^{-1}u\\})$.\n",
    "\n",
    "- Причины неустойчивости:\n",
    "    - если есть два близких собственных числа $\\alpha_i$ и $\\alpha_{i+1}$, то соответствующие векторы $(D - \\alpha_i I)^{-1}u$ и $(D - \\alpha_{i+1} I)^{-1}u$ будут близки, хотя долждны быть ортогональны\n",
    "    - если $\\alpha_i$ и $\\alpha_{i+1}$ очень близки, то они близки к числу $d_i$ между ними, то есть матрицы $D - \\alpha_i I$ и $D - \\alpha_{i+1} I$ близки к вырожденным"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Теорема Лёвнера (Charles Loewner)\n",
    "\n",
    "- Решение проблемы неустойчивости можно получить с помощью теоремы Лёвнера:\n",
    "\n",
    "Если $\\alpha_i$ и $d_i$ удовлетворяют **теореме о чередовании** \n",
    "\n",
    "\n",
    "$$d_n < \\alpha_n < \\ldots < d_{i+1} < \\alpha_{i+1} \\ldots$$\n",
    "\n",
    "Тогда существует вектор $\\widehat{u}$ такой что $\\alpha_i$ – точное собственное значение матрицы\n",
    "\n",
    "$$\\widehat{D} = D + \\widehat{u} \\widehat{u}^*$$\n",
    "\n",
    "и \n",
    "\n",
    "$$ \\widehat{u}^2_i = \\frac{\\prod_{j=1}^n (\\alpha_j - d_i)}{\\prod_{j=1, j \\neq i}^n (d_j - d_i)}  $$\n",
    "\n",
    "- Использование вектора $\\widehat{u}$ вместо $u$ даёт устойчивое вычисление собственного вектора!\n",
    "- Таким образом, сначала вычисляются собственные значения, затем $\\widehat{u}$ и только потом собственные векторы."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Ещё несколько алгоритмов\n",
    "\n",
    "- Абсолютно другой подход основан на идеи **бисекции**\n",
    "\n",
    "- Дана матрица $A$, инерция которой определяется как триплет $(\\nu, \\zeta, \\pi)$, где $\\nu$ – число отрицательных, $\\zeta$ - число нулевых, и $\\pi$ - число положительных собственных значений. \n",
    "\n",
    "- Если $X$ невырождена, тогда\n",
    "\n",
    "$$Inertia(A) = Inertia(X^* A X)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Бисекция с помощью метода Гаусса\n",
    "\n",
    "- Для данного $z$ мы можем запустить метод Гаусса и получить разложение:\n",
    "\n",
    "$$A - zI = L D L^*,$$\n",
    "\n",
    "а инерция для диагональной матрицы вычисляется просто.\n",
    "\n",
    "- С помощью инерции мы можем легко посчитать число собственных значений в заданном интервале.\n",
    "- Пример: если $Inertia(A)=(5,0,2)$ и после сдвига $Inertia(A-zI)=(4,0,3)$, $z\\in[a,b]$ тогда это значит, что $\\lambda(A)\\in[a,z]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Метод Якоби\n",
    "\n",
    "- Вспомним что такое вращения Гивенса (Якоби):\n",
    "\n",
    "на плоскости им соответствуют ортогональные матрицы $2 \\times 2$ вида\n",
    "\n",
    "$$\\begin{pmatrix} \\cos \\phi & \\sin \\phi \\\\ -\\sin \\phi & \\cos \\phi \\end{pmatrix},$$\n",
    "\n",
    "а в $n$-мерном пространстве мы выбираем два индекса $i$ и $j$ и вращаем относительно соответствующих элементов $n$-мерного вектора."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Метод Якоби: продолжение\n",
    "\n",
    "- Идея метода Якоби состоит в минимизации суммы квадратов недиагональных элементов:\n",
    "\n",
    "$$\\Gamma(A) = \\mathrm{off}( U^* A U), \\quad \\mathrm{off}^2(X) = \\sum_{i \\ne j} \\left|X_{ij}\\right|^2 = \\|X \\|^2_F - \\sum\\limits_{i=1}^n x^2_{ii}$$\n",
    "\n",
    "с помощью последовательных вращений Якоби для их зануления. \n",
    "\n",
    "- Когда элементы выбраны, их легко занулить.\n",
    "\n",
    "- **Главный вопрос:** в каком порядке нужно проводить зануление?\n",
    "\n",
    "- Если мы всегда зануляем максимальный недиагональный элемент, метод имеет глобально линейную сходимость и локально квадратичную.\n",
    "\n",
    "- На практике используется циклический порядок (то есть, $(1, 2), (1, 3), \\ldots, (2, 3), \\ldots$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Главное про метод Якоби\n",
    "\n",
    "Метод Якоби был первым численным методом для вычисления собственных значений, предложен в 1846.\n",
    "\n",
    "- Большая константа в оценке сложности\n",
    "- Очень точный (высокая относительная точность для малых собственных значений по сравнению с другими методами)\n",
    "- Хорошая возможность параллелизации"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Главное по сегодняшней лекции\n",
    "\n",
    "- Методы решения симметричной задачи на собстенные значения\n",
    "- Использование сдвигов для ускорения сходимости\n",
    "- Способы вычисления сингулярного разложения"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "    }\n",
       "    div.cell{\n",
       "        /*width:80%;*/\n",
       "        /*margin-left:auto !important;\n",
       "        margin-right:auto;*/\n",
       "    }\n",
       "    h1 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\n",
       "    h2 {\n",
       "        font-family: 'Fenix', serif;\n",
       "    }\n",
       "    h3{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "        margin-top:12px;\n",
       "        margin-bottom: 3px;\n",
       "       }\n",
       "\th4{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "       }\n",
       "    h5 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\t   \n",
       "    div.text_cell_render{\n",
       "        font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "        line-height: 1.2;\n",
       "        font-size: 120%;\n",
       "        /*width:70%;*/\n",
       "        /*margin-left:auto;*/\n",
       "        margin-right:auto;\n",
       "    }\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 50pt;\n",
       "\t\tline-height: 110%;\n",
       "        color:#CD2305;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\t\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 300;\n",
       "        font-size: 16pt;\n",
       "        color: #CD2305;\n",
       "        font-style: italic;\n",
       "        margin-bottom: .5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\n",
       "    \n",
       "    li {\n",
       "        line-height: 110%;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "\n",
       "</style>\n",
       "\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}