seminars/seminar5/seminar5.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Семинар 5. Степенной метод. QR алгоритм"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Экспериментально посмотрим на сходимость степенного метода\n",
    "\n",
    "$$x_{k+1} = \\frac{Ax_k}{\\|Ax_k\\|} $$\n",
    "\n",
    "- Теретическая скорость сходимости **линейная** с фактором $\\frac{|\\lambda_2|}{|\\lambda_1|}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def power_method(A, x0, max_iter, eps):\n",
    "    x = x0.copy()\n",
    "    conv = [x]\n",
    "    for i in range(max_iter):\n",
    "        x = A @ x\n",
    "        x = x / np.linalg.norm(x)\n",
    "        conv.append(x)\n",
    "        eigval = x @ (A @ x)\n",
    "        res = A @ x - eigval * x\n",
    "        if np.linalg.norm(res) < eps:\n",
    "            break\n",
    "    return x, eigval, conv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Симметричная положительно определённая матрица"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "n = 10\n",
    "A = np.random.randn(n, n)\n",
    "A = A.T @ A\n",
    "\n",
    "true_eigvals, true_eigvec = np.linalg.eigh(A)\n",
    "print(true_eigvals)\n",
    "print(true_eigvec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "x0 = np.random.randn(n)\n",
    "max_eigvec, max_eigval, conv = power_method(A, x0, 50, 1e-6)\n",
    "print(max_eigvec)\n",
    "print(max_eigval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.rc(\"text\", usetex=True)\n",
    "\n",
    "power_method_eigval_conv = np.array([np.linalg.norm(A @ x - x @ (A @ x) * x) for x in conv])\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(power_method_eigval_conv, linewidth=4, label=\"Eigval conv\")\n",
    "power_method_eigvec_conv = np.array([np.linalg.norm(x - x[0] / true_eigvec[0, -1] * true_eigvec[:, -1]) for x in conv])\n",
    "plt.plot(power_method_eigvec_conv, linewidth=4, label=\"Eigvec conv\")\n",
    "plt.legend(fontsize=20)\n",
    "plt.yscale(\"log\")\n",
    "plt.ylabel(\"Residual norm\", fontsize=20)\n",
    "plt.yticks(fontsize=16)\n",
    "plt.xlabel(\"\\# iterations\", fontsize=20)\n",
    "plt.xticks(fontsize=16)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "print(power_method_eigval_conv[1:] / power_method_eigval_conv[:-1])\n",
    "print(power_method_eigvec_conv[1:] / power_method_eigvec_conv[:-1])\n",
    "print(true_eigvals[-2] / true_eigvals[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Симметричная матрица\n",
    "\n",
    "- Покажем, что собственные значений вещественные\n",
    "    - $Ax = \\lambda x$\n",
    "    - Умножим слева на $x^*$: $x^* A x = \\lambda \\|x\\|^2_2$\n",
    "    - Возьмём сопряжение от обеих частей: $x^* A^* x = \\lambda^* \\|x\\|_2^2$\n",
    "    - В силу эрмитовости: $x^* A x = \\lambda^* \\|x\\|_2^2$\n",
    "    - Сравнивая со вторым равенством, получаем $\\lambda^* = \\lambda$\n",
    "    - Значит собственные значения действительные"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "A = np.random.randn(n, n)\n",
    "A = A + A.T\n",
    "true_eigvals, true_eigvec = np.linalg.eig(A)\n",
    "print(true_eigvals)\n",
    "print(true_eigvec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "x0 = np.random.randn(n)\n",
    "max_eigvec, max_eigval, conv = power_method(A, x0, 500, 1e-6)\n",
    "print(max_eigvec)\n",
    "print(max_eigval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "power_method_eigval_conv = np.array([np.linalg.norm(A @ x - x @ (A @ x) * x) for x in conv])\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(power_method_eigval_conv, linewidth=4, label=\"Eigval conv\")\n",
    "\n",
    "idx = np.argmax(np.abs(true_eigvals))\n",
    "\n",
    "power_method_eigvec_conv = np.array([np.linalg.norm(x - x[0] / true_eigvec[0, idx] * true_eigvec[:, idx]) for x in conv])\n",
    "plt.plot(power_method_eigvec_conv, linewidth=4, label=\"Eigvec conv\")\n",
    "plt.legend(fontsize=20)\n",
    "plt.yscale(\"log\")\n",
    "plt.ylabel(\"Residual norm\", fontsize=20)\n",
    "plt.yticks(fontsize=16)\n",
    "plt.xlabel(\"\\# iterations\", fontsize=20)\n",
    "plt.xticks(fontsize=16)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "print(power_method_eigval_conv[1:] / power_method_eigval_conv[:-1])\n",
    "print(power_method_eigvec_conv[1:] / power_method_eigvec_conv[:-1])\n",
    "sorted_abs_eigvals = np.sort(np.abs(true_eigvals))\n",
    "print(sorted_abs_eigvals[-2] / sorted_abs_eigvals[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Несимметричная матрица"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "A = np.random.randn(n, n)\n",
    "\n",
    "true_eigvals, true_eigvec = np.linalg.eig(A)\n",
    "print(true_eigvals)\n",
    "print(np.sort(np.abs(true_eigvals)))\n",
    "print(true_eigvec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "x0 = np.random.randn(n)\n",
    "max_eigvec, max_eigval, conv = power_method(A, x0, 500, 1e-6)\n",
    "print(max_eigvec)\n",
    "print(max_eigval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "power_method_eigval_conv = np.array([np.linalg.norm(A @ x - x @ (A @ x) * x) for x in conv])\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(power_method_eigval_conv, label=\"Eigval conv\")\n",
    "\n",
    "idx = np.argmax(np.abs(true_eigvals))\n",
    "\n",
    "power_method_eigvec_conv = np.array([np.linalg.norm(x - x[0] / true_eigvec[0, idx] * true_eigvec[:, idx]) for x in conv])\n",
    "plt.plot(power_method_eigvec_conv, label=\"Eigvec conv\")\n",
    "plt.legend(fontsize=20)\n",
    "plt.yscale(\"log\")\n",
    "plt.ylabel(\"Residual norm\", fontsize=20)\n",
    "plt.yticks(fontsize=16)\n",
    "plt.xlabel(\"\\# iterations\", fontsize=20)\n",
    "plt.xticks(fontsize=16)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "print(power_method_eigval_conv[1:] / power_method_eigval_conv[:-1])\n",
    "print(power_method_eigvec_conv[1:] / power_method_eigvec_conv[:-1])\n",
    "sorted_abs_eigvals = np.sort(np.abs(true_eigvals))\n",
    "print(sorted_abs_eigvals[-2] / sorted_abs_eigvals[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "print(power_method_eigval_conv[1:] / power_method_eigval_conv[:-1])\n",
    "print(power_method_eigvec_conv[1:] / power_method_eigvec_conv[:-1])\n",
    "sorted_abs_eigvals = np.sort(np.abs(true_eigvals))\n",
    "print(sorted_abs_eigvals[-2] / sorted_abs_eigvals[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Если $\\lambda^*_2 = \\lambda_1$, то степенной метод осциллирует и не сходится! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Другие проблемные случаи\n",
    "\n",
    "1) $\\lambda_1 = \\lambda_2$\n",
    "\n",
    "2) $\\lambda_1 = -\\lambda_2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Если $\\lambda_1 = \\lambda_2$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "A = np.random.randn(n, n)\n",
    "Q, _ = np.linalg.qr(A)\n",
    "A = Q @ np.diagflat([2, 2] + list(np.random.rand(n - 2))) @ Q.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "true_eigvals, true_eigvec = np.linalg.eig(A)\n",
    "print(true_eigvals)\n",
    "print(np.abs(true_eigvals))\n",
    "print(true_eigvec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "x0 = np.random.randn(n)\n",
    "max_eigvec, max_eigval, conv = power_method(A, x0, 500, 1e-6)\n",
    "print(max_eigvec)\n",
    "print(max_eigval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "power_method_eigval_conv = np.array([np.linalg.norm(A @ x - x @ (A @ x) * x) for x in conv])\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(power_method_eigval_conv, linewidth=4, label=\"Eigval conv\")\n",
    "\n",
    "idx = np.argmax(np.abs(true_eigvals))\n",
    "\n",
    "power_method_eigvec_conv = np.array([np.linalg.norm(x - x[0] / true_eigvec[0, idx] * true_eigvec[:, idx]) for x in conv])\n",
    "plt.plot(power_method_eigvec_conv, linewidth=4, label=\"Eigvec conv\")\n",
    "plt.legend(fontsize=20)\n",
    "plt.yscale(\"log\")\n",
    "plt.ylabel(\"Residual norm\", fontsize=20)\n",
    "plt.yticks(fontsize=16)\n",
    "plt.xlabel(\"\\# iterations\", fontsize=20)\n",
    "plt.xticks(fontsize=16)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "alpha = np.linalg.lstsq(true_eigvec[:, n-2:], max_eigvec, rcond=None)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "print(\"Accuracy = {}\".format(np.linalg.norm(true_eigvec[:, n-2:] @ alpha - max_eigvec)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Если $\\lambda_1 = -\\lambda_2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "A = np.random.randn(n, n)\n",
    "Q, _ = np.linalg.qr(A)\n",
    "A = Q @ np.diagflat([2, -2] + list(np.random.rand(n - 2))) @ Q.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "true_eigvals, true_eigvec = np.linalg.eig(A)\n",
    "print(true_eigvals)\n",
    "print(np.abs(true_eigvals))\n",
    "print(true_eigvec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "x0 = np.random.randn(n)\n",
    "max_eigvec, max_eigval, conv = power_method(A, x0, 500, 1e-6)\n",
    "print(max_eigvec)\n",
    "print(max_eigval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "power_method_eigval_conv = np.array([np.linalg.norm(A @ x - x @ (A @ x) * x) for x in conv])\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(power_method_eigval_conv, linewidth=4, label=\"Eigval conv\")\n",
    "\n",
    "idx = np.argmax(np.abs(true_eigvals))\n",
    "\n",
    "power_method_eigvec_conv = np.array([np.linalg.norm(x - x[0] / true_eigvec[0, idx] * true_eigvec[:, idx]) for x in conv])\n",
    "plt.plot(power_method_eigvec_conv, linewidth=4, label=\"Eigvec conv\")\n",
    "plt.legend(fontsize=20)\n",
    "plt.yscale(\"log\")\n",
    "plt.ylabel(\"Residual norm\", fontsize=20)\n",
    "plt.yticks(fontsize=16)\n",
    "plt.xlabel(\"\\# iterations\", fontsize=20)\n",
    "plt.xticks(fontsize=16)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Блочный степенной метод\n",
    "\n",
    "- Ищем $k$ cтарших собственных значения и вектора\n",
    "- Этапы аналогичны степенному методу для максимального по модулю собственного значения\n",
    "    - Умножение матрицы на набор векторов\n",
    "    - Аналог нормализации для матриц?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def block_power_method(A, x0, max_iter, eps):\n",
    "    x = x0.copy()\n",
    "    conv_vec = [x]\n",
    "    conv_val = [np.einsum('ik,ij,jk->k',x, A, x)]\n",
    "    for i in range(max_iter):\n",
    "        x = A @ x\n",
    "        x, _ = np.linalg.qr(x)\n",
    "        conv_vec.append(x)\n",
    "        eigval = np.einsum('ik,ij,jk->k', x, A, x)\n",
    "        conv_val.append(eigval)\n",
    "        res = A @ x - (eigval[:, np.newaxis] * x.T).T\n",
    "        if np.linalg.norm(res) < eps:\n",
    "            break\n",
    "    return x, eigval, conv_vec, conv_val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Немного про einsum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "A = np.random.randn(n, n)\n",
    "print(np.trace(A), np.einsum('ii', A))\n",
    "print(np.diag(A) - np.einsum(\"ii->i\", A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "x = np.random.randn(n)\n",
    "print(A @ x - np.einsum('ij,j', A, x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "B = np.random.randn(n, n)\n",
    "print(np.diag(A @ B) - np.einsum('ij,ji->i', A, B))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Более подробное об этой операции см. [тут](https://obilaniu6266h16.wordpress.com/2016/02/04/einstein-summation-in-numpy/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Вернёмся к блочному степенному методу"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "block_size = 3\n",
    "A = np.random.randn(n, n)\n",
    "A = A.T @ A\n",
    "true_eigvals, true_eigvec = np.linalg.eigh(A)\n",
    "print(true_eigvals)\n",
    "print(true_eigvec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "x0 = np.random.randn(n, block_size)\n",
    "max_eigvec, max_eigval, conv_vec, conv_val = block_power_method(A, x0, 500, 1e-6)\n",
    "print(max_eigvec)\n",
    "print(max_eigval)\n",
    "print(true_eigvec[:, n-block_size:n][:, ::-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "power_method_eigval_conv = np.array([np.linalg.norm(A @ x - lam * x) for x, lam in zip(conv_vec, conv_val)])\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(power_method_eigval_conv, linewidth=4, label=\"Eigval conv\")\n",
    "\n",
    "idx = np.argsort(np.abs(true_eigvals))\n",
    "\n",
    "plt.legend(fontsize=20)\n",
    "plt.yscale(\"log\")\n",
    "plt.ylabel(\"Residual norm\", fontsize=20)\n",
    "plt.yticks(fontsize=16)\n",
    "plt.xlabel(\"\\# iterations\", fontsize=20)\n",
    "plt.xticks(fontsize=16)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(power_method_eigval_conv[1:] / power_method_eigval_conv[:-1])\n",
    "print(true_eigvals[-3] / true_eigvals[-2])\n",
    "print(true_eigvals[-1], true_eigvals[-2], true_eigvals[-3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Несимметричный случай"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "A = np.random.randn(n, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "true_eigvals, true_eigvec = np.linalg.eig(A)\n",
    "block_size = 5\n",
    "print(true_eigvals)\n",
    "print(np.sort(np.abs(true_eigvals)))\n",
    "print(true_eigvec.round(5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "x0 = np.random.randn(n, block_size)\n",
    "max_eigvec, max_eigval, conv_vec, conv_val = block_power_method(A, x0, 500, 1e-6)\n",
    "print(max_eigvec)\n",
    "print(max_eigval)\n",
    "print(true_eigvec[:, n-block_size:n][:, ::-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "power_method_eigval_conv = np.array([np.linalg.norm(A @ x - lam * x) for x, lam in zip(conv_vec, conv_val)])\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(power_method_eigval_conv, linewidth=4, label=\"Eigval conv\")\n",
    "\n",
    "idx = np.argsort(np.abs(true_eigvals))\n",
    "\n",
    "plt.legend(fontsize=20)\n",
    "plt.yscale(\"log\")\n",
    "plt.ylabel(\"Residual norm\", fontsize=20)\n",
    "plt.yticks(fontsize=16)\n",
    "plt.xlabel(\"\\# iterations\", fontsize=20)\n",
    "plt.xticks(fontsize=16)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Какие собственные значения удаётся восстановить, а какие нет?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR алгоритм\n",
    "\n",
    "- Приведение матрицы к форме Шура с помощью унитарных преобразований\n",
    "\n",
    "$$ A = UTU^*, $$\n",
    "\n",
    "$T$ – верхнетреугольная матрица.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def qr_algorithm(A, num_iter, eps):\n",
    "    T = A.copy()\n",
    "    U = np.eye(A.shape[0])\n",
    "    conv = [(T, U)]\n",
    "    for i in range(num_iter):\n",
    "        Q, R = np.linalg.qr(T)\n",
    "        T = R @ Q\n",
    "        U = U @ Q\n",
    "        conv.append((T, U))\n",
    "        if np.sum(np.abs(np.tril(T, k=-1))) < eps:\n",
    "            break\n",
    "    return T, U, conv[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "n = 7\n",
    "A = np.random.randn(n, n)\n",
    "# A = A.T @ A\n",
    "A = A + A.T\n",
    "true_eigvals, true_eigvec = np.linalg.eig(A)\n",
    "print(true_eigvals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "T, U, conv = qr_algorithm(A, 2000, 1e-6)\n",
    "print(np.linalg.norm(A - U @ T @ U.T))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "plt.spy(T, markersize=10, precision=1e-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 7))\n",
    "conv_qr = np.array([np.sum(np.abs(np.tril(T, k=-1))) for T, U in conv])\n",
    "plt.plot(conv_qr)\n",
    "plt.yscale(\"log\")\n",
    "plt.xticks(fontsize=20)\n",
    "plt.yticks(fontsize=20)\n",
    "plt.xlabel(\"\\# iterations\", fontsize=20)\n",
    "plt.ylabel(\"Sum of lower triangular part in modulos\", fontsize=20)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "U.round(4), true_eigvec.round(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Стохастическая матрица и собственное значение равное 1\n",
    "\n",
    "**Утверждение.** Пусть дана матрица $A$, в которой элементы неотрицательны и сумма в каждой строке рнавна 1.\n",
    "Тогда 1 максимальное собственное значение такой матрицы."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Доказательство.** \n",
    "- Возьмём вектор $v$ из всех 1 и вычислим $Av$. \n",
    "- В силу равенства 1 суммы элементов в каждой строке $Av$ также будет вектором из всех 1. \n",
    "- Значит вектор $v$ собственный вектор для собственного значения 1\n",
    "- Теперь воспользуемся теоремой Гершгорина: все собственные значения лежат в объединении кругов с центром в точках $a_{ii}$ и радиусом $1 - a_{ii}$\n",
    "- Значит ни один из шаров не имеет точек расположенных дальше 1 и -1\n",
    "- А значит и все собственные значения по модулю меньше 1 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Резюме\n",
    "\n",
    "- Сходимость степенного метода\n",
    "- Сходимость QR алгоритма\n",
    "- Использование спектрального разложения для вычисления матричных рядов"
   ]
  }
 ],
 "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}