12.5.3 The Rayleigh Quotient Iteration.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "    12.5.3 The Rayleigh Quotient Iteration"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "With this notebook, we demonstrate how the Inverse Power Method can be accelerated by shifting the matrix, this time by approximating the smallest eigenvalue with the Rayleigh quotient.\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 Rayleigh Quotient Iteration 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": "markdown",
     "metadata": {},
     "source": [
      "The idea is as follows:\n",
      "\n",
      "The eigenvalues of $ A $ are $ \\lambda_0, \\ldots, \\lambda_3 $ with\n",
      "\n",
      "$$\n",
      "\\vert \\lambda_0 \\vert > \\vert \\lambda_1 \\vert > \\vert \\lambda_2 \\vert > \\vert \\lambda_3 \\vert > 0\n",
      "$$\n",
      "\n",
      "and how fast the iteration converges depends on the ratio \n",
      "\n",
      "$$\n",
      "\\left\\vert \\frac{\\lambda_3}{\\lambda_2} \\right\\vert .\n",
      "$$\n",
      "Now, if you pick a value, $ \\mu $ close to $ \\lambda_3 $, and you iterate with $ A - \\mu I $ (which is known as shifting the matrix/spectrum by $ \\mu $) you can greatly improve the ratio\n",
      "$$\n",
      "\\left\\vert \\frac{\\lambda_3-\\mu}{\\lambda_2-\\mu} \\right\\vert .\n",
      "$$\n",
      "\n",
      "Generally we don't know $ \\lambda_3 $ and hence don't know how to choose $ \\mu $.  But we are generating a vector $ x $ that progressively gets closer and closer to an eigenvector.  Thus, we can use the Rayleigh quotient to approximate an eigenvalue.\n",
      "\n",
      "Here we purposely say \"an eigenvalue\" since it could be that the first random vector $ x $ is close to an eigenvector associated with another eigenvalue, and then we may converge to a different eigenvalue."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Pick a random starting vector\n",
      "\n",
      "x = np.matrix( np.random.rand( 4,1 ) )\n",
      "\n",
      "\n",
      "mu = 0.    # Let's start by not shifting, so hopefully we hone in on the smallest eigenvalue\n",
      "\n",
      "for i in range(0,10):\n",
      "    # We should really compute a factorization of A, but let's be lazy, and compute the inverse\n",
      "    # explicitly\n",
      "    Ainv = np.linalg.inv( A - mu * np.eye( 4, 4 ) )\n",
      "    \n",
      "    x = Ainv * x \n",
      "    \n",
      "    # normalize x to length one\n",
      "    x = x / np.sqrt( np.transpose( x ) * x )\n",
      "    \n",
      "    # Notice we compute the Rayleigh quotient with matrix A, not Ainv.  This is because\n",
      "    # the eigenvector of A is an eigenvector of Ainv\n",
      "    \n",
      "    mu = np.transpose( x ) * A * x\n",
      "    \n",
      "    # The above returns a 1 x 1 matrix.  Let's set mu to the scalar\n",
      "    \n",
      "    mu = mu[ 0, 0 ]\n",
      "    \n",
      "    print( 'Rayleigh quotient with vector x:', np.transpose( x ) * A * x / ( np.transpose( x ) * x ))\n",
      "    print( 'inner product of x with v3     :', np.transpose( x ) * V[ :, 3 ] )\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 may converge to 1.0 (but it may converge to another eigenvalue!).\n",
      " </li>\n",
      " <li>\n",
      " The inner product of $ x $ and the last column of $ V $, $ v_{n-1} $, may converge to 1 or -1 since eventually $ x $ may be in the direction of $ v_{n-1} $ (or in the opposite direction).  But not if we start converging to another eigenvalue... If this happens, try rerunning all the code blocks above to get a different $V$ matrix.\n",
      " </li>\n",
      " </ul>\n"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}