{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "    12.5.2 Shifting the Inverse Power Method"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "With this notebook, we demonstrate how the Inverse Power Method can be accelerated by shifting the matrix.\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 Shifted Inverse 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": "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",
      "Try different values of $ \\mu$.  What if you pick $ \\mu \\approx 2 $?  \n",
      "What if you pick $ \\mu = 0.8 $?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Pick a random starting vector\n",
      "\n",
      "x = np.matrix( np.random.rand( 4,1 ) )\n",
      "\n",
      "# We should really compute a factorization of A, but let's be lazy, and compute the inverse\n",
      "# explicitly\n",
      "\n",
      "mu = 0.8\n",
      "\n",
      "Ainv = np.linalg.inv( A - mu * np.eye( 4, 4 ) )\n",
      "\n",
      "for i in range(0,10):\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",
      "    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 should converge to 1.0 (quickly if $ \\mu \\approx 1 $).\n",
      " </li>\n",
      " <li>\n",
      " The inner product of $ x $ and the last column of $ V $, $ v_{n-1} $, should converge to 1 or -1 since eventually $ x $ should be in the direction of $ v_{n-1} $ (or in the opposite direction).\n",
      " </li>\n",
      " </ul>\n",
      " \n",
      " This time, if you change the \"2\" on the diagonal to \"-1\", you still converge to $ v_{n-1} $ because for the matrix $ A - \\mu I $, $ -1 - \\mu $ is not as small as $ 1 - \\mu $ (in magnitude).\n",
      "\n",
      " You can check this by looking at $ ( I - V_R ( V_R^T V_R )^{-1} V_R^T ) x $, where $V_R $ equals the matrix with $ v_2 $ and $ v_3 $ as its columns, to see if the vector orthogonal to $ {\\cal C}( V_R ) $ converges to zero. This is seen in the following code block:\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "w = x - V[ :,2:4 ] * np.linalg.inv( np.transpose( V[ :,2:4 ] ) * V[ :,2:4 ] ) * np.transpose( V[ :,2:4 ] ) * 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": {}
  }
 ]
}