{ "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", " Be sure to make a copy!!!! " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We start by creating a matrix with known eigenvalues and eigenvectors\n", "\n", "How do we do this? \n", "\n", "\n", " 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? " ] }, { "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", " \n" ] } ], "metadata": {} } ] }