{"nbformat": 3, "nbformat_minor": 0, "metadata": {"signature": "sha256:3e4004161fab291bef880f3f60be8c508e302bd24d2b7f820ec8a9adaf1fc667", "name": ""}, "worksheets": [{"metadata": {}, "cells": [{"source": ["# Power Iteration and its Variants"], "metadata": {}, "cell_type": "markdown"}, {"cell_type": "code", "outputs": [], "input": ["import numpy as np\n", "import numpy.linalg as la"], "collapsed": false, "metadata": {}, "language": "python", "prompt_number": 2}, {"source": ["Let's  prepare a matrix with some random or deliberately chosen eigenvalues:"], "metadata": {}, "cell_type": "markdown"}, {"cell_type": "code", "outputs": [{"stream": "stdout", "text": ["[-2.667651   -0.95797093 -0.33019549 -0.29151942 -0.18635343 -0.14418093]\n"], "output_type": "stream"}], "input": ["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])"], "collapsed": false, "metadata": {}, "language": "python", "prompt_number": 6}, {"source": ["Let's also pick an initial vector:"], "metadata": {}, "cell_type": "markdown"}, {"cell_type": "code", "outputs": [{"output_type": "pyout", "text": ["array([ 2.26930477,  0.66356156,  0.8991019 , -0.36580094,  0.46269004,\n", "        0.079874  ])"], "metadata": {}, "prompt_number": 4}], "input": ["x0 = np.random.randn(n)\n", "x0"], "collapsed": false, "metadata": {}, "language": "python", "prompt_number": 4}, {"source": ["### Power iteration"], "metadata": {}, "cell_type": "markdown"}, {"cell_type": "code", "outputs": [], "input": ["x = x0"], "collapsed": false, "metadata": {}, "language": "python", "prompt_number": 5}, {"source": ["Now implement plain power iteration.\n", "\n", "Run the below cell in-place (Ctrl-Enter) many times."], "metadata": {}, "cell_type": "markdown"}, {"cell_type": "code", "outputs": [{"output_type": "pyout", "text": ["array([-6.40496816,  7.1851682 ,  1.97149585,  4.77787616,  3.09099127,\n", "        4.33054803])"], "metadata": {}, "prompt_number": 6}], "input": ["x = np.dot(A, x)\n", "x"], "collapsed": false, "metadata": {}, "language": "python", "prompt_number": 6}, {"source": ["* What's the problem with this method?\n", "* Does anything useful come of this?\n", "* How do we fix it?"], "metadata": {}, "cell_type": "markdown"}, {"source": ["### Normalized power iteration"], "metadata": {}, "cell_type": "markdown"}, {"source": ["Back to the beginning: Reset to the initial vector."], "metadata": {}, "cell_type": "markdown"}, {"cell_type": "code", "outputs": [], "input": ["x = x0/la.norm(x0)"], "collapsed": false, "metadata": {}, "language": "python", "prompt_number": 7}, {"source": ["Implement normalized power iteration.\n", "\n", "Run this cell in-place (Ctrl-Enter) many times."], "metadata": {}, "cell_type": "markdown"}, {"cell_type": "code", "outputs": [{"stream": "stdout", "text": ["4.67639723405\n", "[-0.52706768  0.59127069  0.16223527  0.39317355  0.25435905  0.35636272]\n"], "output_type": "stream"}], "input": ["x = np.dot(A, x)\n", "nrm = la.norm(x)\n", "x = x/nrm\n", "\n", "print(nrm)\n", "print(x)"], "collapsed": false, "metadata": {}, "language": "python", "prompt_number": 8}, {"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)"], "metadata": {}, "cell_type": "markdown"}, {"source": ["------\n", "What if we want the smallest eigenvalue (by magnitude)?\n", "\n", "Once again, reset to the beginning."], "metadata": {}, "cell_type": "markdown"}, {"cell_type": "code", "outputs": [], "input": ["x = x0/la.norm(x0)"], "collapsed": false, "metadata": {}, "language": "python", "prompt_number": 9}, {"source": ["Run the cell below in-place many times."], "metadata": {}, "cell_type": "markdown"}, {"cell_type": "code", "outputs": [{"stream": "stdout", "text": ["0.0464410512401\n", "[ 0.07260214 -0.82392938  0.17918176 -0.52684691  0.07819739  0.00898414]\n"], "output_type": "stream"}], "input": ["x = la.solve(A, x)\n", "nrm = la.norm(x)\n", "x = x/nrm\n", "\n", "print(1/nrm)\n", "print(x)"], "collapsed": false, "metadata": {}, "language": "python", "prompt_number": 10}, {"source": ["* What's the computational cost per iteration?\n", "* Can we make this method search for a specific eigenvalue?\n", "* What is this [method](https://en.wikipedia.org/wiki/Inverse_iteration) called?"], "metadata": {}, "cell_type": "markdown"}, {"source": ["--------------\n", "Can we feed an estimate of the current approximate eigenvalue back into the calculation? (Hint: Rayleigh quotient)\n", "\n", "Reset once more."], "metadata": {}, "cell_type": "markdown"}, {"cell_type": "code", "outputs": [], "input": ["x = x0/la.norm(x0)"], "collapsed": false, "metadata": {}, "language": "python", "prompt_number": 11}, {"source": ["Run this cell in-place (Ctrl-Enter) many times."], "metadata": {}, "cell_type": "markdown"}, {"cell_type": "code", "outputs": [{"stream": "stdout", "text": ["-1.17969302281\n", "[-0.06324034 -0.54447197  0.4628917  -0.37622458  0.41482856  0.41431213]\n"], "output_type": "stream"}], "input": ["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", "\n", "print(sigma)\n", "print(x)"], "collapsed": false, "metadata": {}, "language": "python", "prompt_number": 12}, {"source": ["* What's this [method](https://en.wikipedia.org/wiki/Rayleigh_quotient_iteration) called?\n", "* What's a reasonable stopping criterion?\n", "* Computational downside of this iteration?"], "metadata": {}, "cell_type": "markdown"}]}]}