In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [54]:
def powermethod(A, x, lambdas=None, maxiter=10000):
    for i in range(maxiter):
        y = A.dot(x)
        thislambda = y[0] / x[0]
        x = y / np.abs(y).max()
        
        if lambdas is not None:
            lambdas.append(thislambda)
    return thislambda, x
In [55]:
A = np.array([[0.999,0],[0,1.0]])
Q, _ = np.linalg.qr(np.random.rand(2,2))
A = Q.dot(A).dot(Q.T)
print(np.linalg.eig(A))
lambdas = []
l, x = powermethod(A, x=np.random.rand(2), lambdas=lambdas)
(array([ 1.   ,  0.999]), array([[ 0.87675537,  0.48093662],
       [-0.48093662,  0.87675537]]))
In [57]:
A = np.array([[0.99,0],[0,1.0]])
Q, _ = np.linalg.qr(np.random.rand(2,2))
A = Q.dot(A).dot(Q.T)
print(np.linalg.eig(A))
lambdas2 = []
l, x = powermethod(A, x=np.random.rand(2), lambdas=lambdas2)
(array([ 1.  ,  0.99]), array([[ 0.99707147,  0.07647534],
       [-0.07647534,  0.99707147]]))
In [58]:
plt.semilogy(np.abs(1.0-np.array(lambdas)))
plt.semilogy(np.abs(1.0-np.array(lambdas2)))
Out[58]:
[<matplotlib.lines.Line2D at 0x105315c50>]
In [ ]:
lambdas =