{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Power Iteration and its Variants" ] }, { "cell_type": "code", "execution_count": 104, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import numpy.linalg as la\n", "np.set_printoptions(precision=3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's prepare a matrix with some random or deliberately chosen eigenvalues:" ] }, { "cell_type": "code", "execution_count": 105, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-2.668 -0.958 -0.33 -0.292 -0.186 -0.144]\n" ] } ], "source": [ "n = 6\n", "\n", "if 1:\n", " np.random.seed(70)\n", " eigvecs = np.random.randn(n, n)\n", " eigvals = np.sort(np.random.randn(n))\n", " # Uncomment for near-duplicate largest-magnitude eigenvalue\n", " # eigvals[1] = eigvals[0] + 1e-3\n", "\n", " A = eigvecs.dot(np.diag(eigvals)).dot(la.inv(eigvecs))\n", " print(eigvals)\n", " \n", "else:\n", " # Complex eigenvalues\n", " np.random.seed(40)\n", " A = np.random.randn(n, n)\n", " print(la.eig(A)[0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's also pick an initial vector:" ] }, { "cell_type": "code", "execution_count": 106, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 2.269, 0.664, 0.899, -0.366, 0.463, 0.08 ])" ] }, "execution_count": 106, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x0 = np.random.randn(n)\n", "x0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Power iteration" ] }, { "cell_type": "code", "execution_count": 107, "metadata": {}, "outputs": [], "source": [ "x = x0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now implement plain power iteration." ] }, { "cell_type": "code", "execution_count": 127, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 1.329e+12 -4.740e+12 -7.728e+11 1.018e+12 -8.466e+11 -2.213e+12]\n", "[-3.545e+12 1.264e+13 2.062e+12 -2.716e+12 2.258e+12 5.903e+12]\n", "[ 9.457e+12 -3.373e+13 -5.500e+12 7.245e+12 -6.025e+12 -1.575e+13]\n", "[-2.523e+13 8.998e+13 1.467e+13 -1.933e+13 1.607e+13 4.201e+13]\n", "[ 6.730e+13 -2.400e+14 -3.914e+13 5.156e+13 -4.287e+13 -1.121e+14]\n", "[-1.795e+14 6.403e+14 1.044e+14 -1.375e+14 1.144e+14 2.989e+14]\n", "[ 4.789e+14 -1.708e+15 -2.785e+14 3.669e+14 -3.051e+14 -7.975e+14]\n", "[-1.278e+15 4.557e+15 7.430e+14 -9.788e+14 8.139e+14 2.127e+15]\n", "[ 3.408e+15 -1.216e+16 -1.982e+15 2.611e+15 -2.171e+15 -5.675e+15]\n", "[-9.092e+15 3.243e+16 5.288e+15 -6.965e+15 5.792e+15 1.514e+16]\n", "[ 2.425e+16 -8.650e+16 -1.411e+16 1.858e+16 -1.545e+16 -4.039e+16]\n", "[-6.470e+16 2.308e+17 3.763e+16 -4.957e+16 4.122e+16 1.077e+17]\n", "[ 1.726e+17 -6.156e+17 -1.004e+17 1.322e+17 -1.100e+17 -2.874e+17]\n", "[-4.604e+17 1.642e+18 2.678e+17 -3.527e+17 2.933e+17 7.667e+17]\n", "[ 1.228e+18 -4.381e+18 -7.143e+17 9.410e+17 -7.825e+17 -2.045e+18]\n", "[-3.277e+18 1.169e+19 1.906e+18 -2.510e+18 2.087e+18 5.456e+18]\n", "[ 8.741e+18 -3.117e+19 -5.083e+18 6.696e+18 -5.569e+18 -1.455e+19]\n", "[-2.332e+19 8.316e+19 1.356e+19 -1.786e+19 1.486e+19 3.883e+19]\n", "[ 6.220e+19 -2.219e+20 -3.618e+19 4.765e+19 -3.963e+19 -1.036e+20]\n", "[-1.659e+20 5.918e+20 9.650e+19 -1.271e+20 1.057e+20 2.763e+20]\n" ] } ], "source": [ "for i in range(20):\n", " x = A @ x\n", " print(x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* What's the problem with this method?\n", "* Does anything useful come of this?\n", "* How do we fix it?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Normalized power iteration" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Back to the beginning: Reset to the initial vector." ] }, { "cell_type": "code", "execution_count": 128, "metadata": {}, "outputs": [], "source": [ "x = x0/la.norm(x0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Implement normalized power iteration." ] }, { "cell_type": "code", "execution_count": 129, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-0.285 0.819 0.106 -0.172 0.149 0.431]\n", "[ 0.243 -0.842 -0.13 0.173 -0.156 -0.402]\n", "[-0.238 0.845 0.135 -0.177 0.153 0.396]\n", "[ 0.237 -0.845 -0.137 0.18 -0.152 -0.395]\n", "[-0.237 0.845 0.137 -0.181 0.151 0.395]\n", "[ 0.237 -0.845 -0.138 0.181 -0.151 -0.394]\n", "[-0.237 0.845 0.138 -0.181 0.151 0.394]\n", "[ 0.237 -0.845 -0.138 0.181 -0.151 -0.394]\n", "[-0.237 0.845 0.138 -0.181 0.151 0.394]\n", "[ 0.237 -0.845 -0.138 0.181 -0.151 -0.394]\n", "2.667662268743778\n" ] } ], "source": [ "for i in range(10):\n", " x = A @ x\n", " nrm = la.norm(x)\n", " x = x/nrm\n", " print(x)\n", "\n", "print(nrm)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* What do you observe about the norm?\n", "* What about the sign?\n", "* What is the vector $x$ now?\n", "\n", "Extensions:\n", "\n", "* Now try the matrix variants above.\n", "* Suggest a better way of estimating the eigenvalue. [Hint](https://en.wikipedia.org/wiki/Rayleigh_quotient)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "------\n", "### Inverse iteration\n", "\n", "What if we want the smallest eigenvalue (by magnitude)?\n", "\n", "Once again, reset to the beginning." ] }, { "cell_type": "code", "execution_count": 143, "metadata": {}, "outputs": [], "source": [ "x = x0/la.norm(x0)" ] }, { "cell_type": "code", "execution_count": 144, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 0.344 -0.55 -0.094 0.23 -0.05 -0.718]\n", "[-0.527 0.236 0.046 -0.123 -0.066 0.803]\n", "[ 0.599 -0.01 -0.037 -0.031 0.127 -0.789]\n", "[-0.617 -0.128 0.043 0.173 -0.153 0.74 ]\n", "[ 0.611 0.213 -0.052 -0.284 0.161 -0.687]\n", "[-0.597 -0.265 0.061 0.365 -0.161 0.64 ]\n", "[ 0.583 0.299 -0.067 -0.424 0.16 -0.601]\n", "[-0.57 -0.322 0.072 0.465 -0.157 0.571]\n", "[ 0.558 0.337 -0.076 -0.496 0.155 -0.547]\n", "[-0.55 -0.348 0.078 0.517 -0.153 0.529]\n", "[ 0.543 0.356 -0.08 -0.533 0.151 -0.515]\n", "[-0.537 -0.362 0.082 0.545 -0.15 0.504]\n", "[ 0.533 0.366 -0.083 -0.554 0.149 -0.496]\n", "[-0.529 -0.369 0.084 0.561 -0.148 0.49 ]\n", "[ 0.527 0.371 -0.084 -0.566 0.147 -0.485]\n", "[-0.525 -0.373 0.085 0.57 -0.147 0.481]\n", "[ 0.523 0.375 -0.085 -0.573 0.146 -0.478]\n", "[-0.522 -0.376 0.086 0.575 -0.146 0.476]\n", "[ 0.521 0.376 -0.086 -0.577 0.146 -0.474]\n", "[-0.521 -0.377 0.086 0.578 -0.146 0.473]\n", "[ 0.52 0.377 -0.086 -0.579 0.146 -0.472]\n", "[-0.52 -0.378 0.086 0.58 -0.145 0.471]\n", "[ 0.519 0.378 -0.086 -0.581 0.145 -0.471]\n", "[-0.519 -0.378 0.086 0.581 -0.145 0.47 ]\n", "[ 0.519 0.379 -0.086 -0.582 0.145 -0.47 ]\n", "[-0.519 -0.379 0.086 0.582 -0.145 0.47 ]\n", "[ 0.519 0.379 -0.086 -0.582 0.145 -0.469]\n", "[-0.518 -0.379 0.086 0.582 -0.145 0.469]\n", "[ 0.518 0.379 -0.086 -0.582 0.145 -0.469]\n", "[-0.518 -0.379 0.086 0.583 -0.145 0.469]\n" ] } ], "source": [ "for i in range(30):\n", " x = la.solve(A, x)\n", " nrm = la.norm(x)\n", " x = x/nrm\n", " print(x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* What's the computational cost per iteration?\n", "* Can we make this method search for a specific eigenvalue?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "------\n", "### Inverse Shift iteration\n", "\n", "What if we want the eigenvalue closest to a give value $\\sigma$?\n", "\n", "Once again, reset to the beginning." ] }, { "cell_type": "code", "execution_count": 148, "metadata": {}, "outputs": [], "source": [ "x = x0/la.norm(x0)" ] }, { "cell_type": "code", "execution_count": 149, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 0.011 -0.809 -0.191 0.223 -0.186 -0.474]\n", "[-0.127 0.775 0.149 -0.232 0.16 0.531]\n", "[ 0.192 -0.728 -0.123 0.246 -0.135 -0.582]\n", "[-0.246 0.674 0.103 -0.261 0.108 0.628]\n", "[ 0.295 -0.616 -0.085 0.271 -0.08 -0.669]\n", "[-0.338 0.556 0.07 -0.277 0.052 0.702]\n", "[ 0.376 -0.497 -0.057 0.277 -0.026 -0.728]\n", "[-0.41 0.441 0.047 -0.273 0.002 0.749]\n", "[ 0.44 -0.387 -0.038 0.264 0.019 -0.765]\n", "[-0.465 0.338 0.031 -0.252 -0.039 0.777]\n", "[ 0.488 -0.291 -0.026 0.237 0.056 -0.786]\n", "[-0.508 0.249 0.022 -0.219 -0.071 0.792]\n", "[ 0.525 -0.209 -0.019 0.199 0.084 -0.796]\n", "[-0.54 0.172 0.017 -0.178 -0.095 0.798]\n", "[ 0.553 -0.138 -0.016 0.155 0.105 -0.799]\n", "[-0.565 0.107 0.016 -0.132 -0.114 0.799]\n", "[ 0.575 -0.078 -0.016 0.108 0.122 -0.798]\n", "[-0.584 0.05 0.016 -0.084 -0.129 0.796]\n", "[ 0.591 -0.025 -0.017 0.06 0.135 -0.792]\n", "[-0.597 0.001 0.018 -0.035 -0.14 0.789]\n", "[ 0.603 0.021 -0.02 0.011 0.144 -0.784]\n", "[-0.607 -0.041 0.022 0.013 -0.148 0.779]\n", "[ 0.611 0.061 -0.024 -0.036 0.151 -0.774]\n", "[-0.614 -0.079 0.025 0.059 -0.154 0.768]\n", "[ 0.616 0.096 -0.027 -0.082 0.156 -0.761]\n", "[-0.618 -0.112 0.03 0.103 -0.158 0.755]\n", "[ 0.619 0.126 -0.032 -0.124 0.16 -0.748]\n", "[-0.619 -0.14 0.034 0.145 -0.161 0.741]\n", "[ 0.619 0.154 -0.036 -0.164 0.162 -0.734]\n", "[-0.619 -0.166 0.038 0.183 -0.163 0.726]\n" ] } ], "source": [ "sigma = 1\n", "A_sigma = A-sigma*np.eye(A.shape[0])\n", "for i in range(30):\n", " x = la.solve(A_sigma, x)\n", " nrm = la.norm(x)\n", " x = x/nrm\n", " print(x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "--------------\n", "### Rayleigh quotient iteration\n", "\n", "Can we feed an estimate of the current approximate eigenvalue back into the calculation? (Hint: Rayleigh quotient)\n", "\n", "Reset once more." ] }, { "cell_type": "code", "execution_count": 153, "metadata": {}, "outputs": [], "source": [ "x = x0/la.norm(x0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Run this cell in-place (Ctrl-Enter) many times." ] }, { "cell_type": "code", "execution_count": 154, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 0.063 -0.792 -0.173 0.23 -0.173 -0.505]\n", "[-0.191 0.726 0.126 -0.248 0.133 0.585]\n", "[ 0.521 -0.254 -0.051 0.121 0.058 -0.802]\n", "[-0.544 -0.35 0.081 0.533 -0.15 0.519]\n", "[ 0.52 0.378 -0.086 -0.58 0.146 -0.472]\n", "[-0.518 -0.379 0.086 0.583 -0.145 0.469]\n", "[ 0.518 0.379 -0.086 -0.583 0.145 -0.469]\n", "[-0.518 -0.379 0.086 0.583 -0.145 0.469]\n", "[-0.518 -0.379 0.086 0.583 -0.145 0.469]\n", "[ 0.518 0.379 -0.086 -0.583 0.145 -0.469]\n" ] } ], "source": [ "for i in range(10):\n", " sigma = np.dot(x, np.dot(A, x))/np.dot(x, x)\n", " x = la.solve(A-sigma*np.eye(n), x)\n", " x = x/la.norm(x)\n", " print(x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* What's a reasonable stopping criterion?\n", "* Computational downside of this iteration?" ] } ], "metadata": { "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.5" } }, "nbformat": 4, "nbformat_minor": 1 }