inverse_power_method.py

#function [lamda_1, v_1, success]  =  inverse_power_method...
#    (lamda, q0, A, tol, maxitr)
#  Accompanying program for the text
#
#     Classical and Modern Numerical Analysis:
#     Theory, Methods and Practice
#  by Azmy S. Ackleh, Edward J. Allen,
#     R. Baker Kearfott, and Padmanabhan Seshaiyer
#
#     (Taylor and Francis / CRC Press, 2009)
#
#  [lamda_1, v_1, success] = inverse_power_method(lamda, q0, A, tol, maxitr)
#  computes an eigenvalue and eigenvector of the matrix A, according to the
#  inverse power method as described in Section 5.3 (starting on page 303)
#  of the text.  The eigenvalue is approximated using a Rayleigh quotient
#  (described on page 300 of the text).  Iteration stops when either maxitr
#  iterations of (5.19) have been reached, or the magnitude of the
#  difference between successive eigenvalue approximations is less than tol.
#  'success' is set to true if the difference between successive
#  approximations is less than tol upon return, and is set to 'false'
#  otherwise.  In any case, lamda_1 is set to the most recent approximation
#  to the eigenvalue and v_1 is set to the most recent approximation to the
#  eigenvector, normalized to infinity norm 1.

import numpy as np

def inverse_power_method(lamda, q0, A, tol, maxitr):  # lambda is a reserved word in Python
        n = A.shape[0]
        a = np.linalg.inv(A-lamda*np.eye(n))  # inefficient to actually compute the inverse
        q_nu  =  q0
        alam = 2*np.linalg.norm(q_nu,np.inf)
        check = 1
        success  =  False
        for k in range(maxitr):
                # k = k+1  Why was this here?
                alam2 = alam
                q_nu = np.dot(a,q_nu)
                q_nu = q_nu/np.linalg.norm(q_nu,np.inf)
                alam = np.dot( q_nu, np.dot(a,q_nu) ) / np.dot(q_nu,q_nu)
                check = np.linalg.norm(alam-alam2)
                if (check < tol):
                    success  =  True
                    break

        lamda_1 = lamda + 1/alam
        v_1  =  q_nu
        print('{:9.0f} {:15.4f} {:15.5f}'.format(k,lamda,lamda_1) )
        return lamda_1, v_1, success


if __name__ == '__main__':

        A = np.array( [[1,0,0],[0,.375,0],[0,0,-2]] )
        print( inverse_power_method(.3, np.random.rand(3), A, 1.e-7, 100) )

'''
$ py inverse_power_method.py
        4          0.3000         0.37500
(0.37500000000216793, array([ 5.68986953e-06,  1.00000000e+00, -1.61848527e-08]), True)
'''