12.4.2 The Power Method.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "    12.4.2 The Power Method"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "With this notebook, we demonstrate how the Power Method can be used to compute the eigenvector associated with the largest eigenvalue (in magnitude).\n",
      "\n",
      "<font color=red> Be sure to make a copy!!!! </font>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We start by creating a matrix with known eigenvalues and eigenvectors\n",
      "\n",
      "How do we do this?  \n",
      "<ul>\n",
      "  <li>\n",
      "    We want a matrix that is not deficient, since otherwise the Power Method may not work. \n",
      "  </li>\n",
      "  <li>\n",
      "    Hence, $ A = V \\Lambda V^{-1} $ for some diagonal matrix $ \\Lambda $ and nonsingular matrix $ V $.  The eigenvalues are then on the diagonal of $ \\Lambda $ and the eigenvectors are the columns of $ V $.\n",
      "    </li>\n",
      "    <li>\n",
      "    So, let's pick the eigenvalues for the diagonal of $ \\Lambda $ and let's pick a random matrix $ V $ (in the hopes that it has linearly independent columns) and then let's see what happens.  \n",
      "    </li>\n",
      "    </ul>\n",
      "\n",
      "<font color=red> Experiment by changing the eigenvalues!  What happens if you make the second entry on the diagonal equal to -4?  Or what if you set 2 to -1? </font>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import laff\n",
      "import flame\n",
      "\n",
      "Lambda = np.matrix( ' 4., 0., 0., 0;\\\n",
      "                      0., 3., 0., 0;\\\n",
      "                      0., 0., 2., 0;\\\n",
      "                      0., 0., 0., 1' )\n",
      "\n",
      "lambda0 = Lambda[ 0,0 ]\n",
      "\n",
      "V = np.matrix( np.random.rand( 4,4 ) )\n",
      "\n",
      "# normalize the columns of V to equal one\n",
      "\n",
      "for j in range( 0, 4 ):\n",
      "    V[ :, j ] = V[ :, j ] / np.sqrt( np.transpose( V[:,j] ) * V[:, j ] )\n",
      "\n",
      "A = V * Lambda * np.linalg.inv( V )\n",
      "\n",
      "print( 'Lambda = ' )\n",
      "print( Lambda)\n",
      "\n",
      "print( 'V = ' )\n",
      "print( V )\n",
      "\n",
      "print( 'A = ' )\n",
      "print( A )\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Pick a random starting vector\n",
      "\n",
      "x = np.matrix( np.random.rand( 4,1 ) )\n",
      "\n",
      "for i in range(0,10):\n",
      "    x = A * x \n",
      "    \n",
      "    # normalize x to length one\n",
      "    x = x / np.sqrt( np.transpose( x ) * x )\n",
      "    \n",
      "    print( 'Rayleigh quotient with vector x:', np.transpose( x ) * A * x / ( np.transpose( x ) * x ))\n",
      "    print( 'inner product of x with v0     :', np.transpose( x ) * V[ :, 0 ] )\n",
      "    print( ' ' )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the above, \n",
      " <ul>\n",
      " <li>\n",
      " The Rayleigh quotient should converge to 4.0 (slowly).\n",
      " </li>\n",
      " <li>\n",
      " The inner product of $ x $ and the first column of $ V $, $ v_0 $, should converge to 1 or -1 since eventually $ x $ should be in the direction of $ v_0 $ (or in the opposite direction).\n",
      " </li>\n",
      " </ul>\n",
      " \n",
      "If you change the \"3\" on the diagonal to \"-4\", then you have two largest eigenvalues (in magnitude), and the vector $ x $ will end up in the space spanned by $ v_0 $ and $ v_1 $. \n",
      " You can check this by looking at $ ( I - V_L ( V_L^T V_L )^{-1} V_L^T ) x $, where $V_L $ equals the matrix with $ v_0 $ and $ v_1 $ as its columns, to see if the vector orthogonal to $ {\\cal C}( V_L ) $ converges to zero. This is seen in the following code block:\n",
      "\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "w = x - V[ :,0:2 ] * np.linalg.inv( np.transpose( V[ :,0:2 ] ) * V[ :,0:2 ] ) * np.transpose( V[ :,0:2 ] ) * x\n",
      "    \n",
      "print( 'Norm of component orthogonal: ', np.linalg.norm( w ) )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}