{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from __future__ import division\n",
      "\n",
      "import numpy as np\n",
      "import numpy.linalg as la"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "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 = np.dot(la.solve(eigvecs, np.diag(eigvals)), 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]\n",
      "\n",
      "x0 = np.random.randn(n)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[-2.667651   -0.95797093 -0.33019549 -0.29151942 -0.18635343 -0.14418093]\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Power iteration"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = x0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Run this cell in-place (Ctrl-Enter) many times.\n",
      "x = np.dot(A, x)\n",
      "x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 48,
       "text": [
        "array([ -4.28616666e+19,  -4.53997702e+19,   8.19441238e+19,\n",
        "        -3.31918771e+19,   7.97142990e+19,   8.57108316e+19])"
       ]
      }
     ],
     "prompt_number": 48
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Normalized power iteration"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = x0/la.norm(x0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 49
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Run this cell in-place (Ctrl-Enter) many times.\n",
      "\n",
      "# x is unit-length here\n",
      "\n",
      "x = np.dot(A, x)\n",
      "nrm = la.norm(x)\n",
      "x = x/nrm\n",
      "\n",
      "print nrm\n",
      "print x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "2.66765099561\n",
        "[ 0.2688559   0.28477652 -0.51400617  0.20820077 -0.50001928 -0.53763338]\n"
       ]
      }
     ],
     "prompt_number": 115
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "What if we want smallest, not largest?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = x0/la.norm(x0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 117
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Run this cell in-place (Ctrl-Enter) many times.\n",
      "\n",
      "x = ...\n",
      "nrm = la.norm(x)\n",
      "x = x/nrm\n",
      "\n",
      "print 1/nrm\n",
      "print x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.144180924949\n",
        "[-0.26610396 -0.17430691  0.45852284 -0.11352085  0.51716735  0.63891591]\n"
       ]
      }
     ],
     "prompt_number": 191
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<!--\n",
      "x = la.solve(A, x)\n",
      "-->\n",
      "(Edit this cell for solution.)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$\\uparrow$ Inverse iteration"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "How do we use the Rayleigh quotient?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = x0/la.norm(x0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 197
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Run this cell in-place (Ctrl-Enter) many times.\n",
      "\n",
      "sigma = np.dot(x, np.dot(A, x))/np.dot(x, x)\n",
      "x = ...\n",
      "x = x/la.norm(x)\n",
      "\n",
      "print sigma\n",
      "print x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "-0.957970929184\n",
        "[ 0.0539665   0.58040235 -0.44257009  0.39097236 -0.39553142 -0.39376129]\n"
       ]
      }
     ],
     "prompt_number": 208
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<!--\n",
      "x = la.solve(A-sigma*np.eye(n), x)\n",
      "-->\n",
      "(Edit this cell for solution.)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$\\uparrow$ Rayleigh quotient iteration"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "* What's a reasonable stopping criterion?\n",
      "* Computational downside of Rayleigh quotient iteration?"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}