import numpy as np
import numpy.linalg as la
np.set_printoptions(precision=3)
Let's prepare a matrix with some random or deliberately chosen eigenvalues:
n = 6
if 1:
np.random.seed(70)
eigvecs = np.random.randn(n, n)
eigvals = np.sort(np.random.randn(n))
# Uncomment for near-duplicate largest-magnitude eigenvalue
# eigvals[1] = eigvals[0] + 1e-3
A = eigvecs.dot(np.diag(eigvals)).dot(la.inv(eigvecs))
print(eigvals)
else:
# Complex eigenvalues
np.random.seed(40)
A = np.random.randn(n, n)
print(la.eig(A)[0])
Let's also pick an initial vector:
x0 = np.random.randn(n)
x0
x = x0
Now implement plain power iteration.
for i in range(20):
x = A @ x
print(x)
Back to the beginning: Reset to the initial vector.
x = x0/la.norm(x0)
Implement normalized power iteration.
for i in range(10):
x = A @ x
nrm = la.norm(x)
x = x/nrm
print(x)
print(nrm)
Extensions:
What if we want the smallest eigenvalue (by magnitude)?
Once again, reset to the beginning.
x = x0/la.norm(x0)
for i in range(30):
x = la.solve(A, x)
nrm = la.norm(x)
x = x/nrm
print(x)
What if we want the eigenvalue closest to a give value $\sigma$?
Once again, reset to the beginning.
x = x0/la.norm(x0)
sigma = 1
A_sigma = A-sigma*np.eye(A.shape[0])
for i in range(30):
x = la.solve(A_sigma, x)
nrm = la.norm(x)
x = x/nrm
print(x)
Can we feed an estimate of the current approximate eigenvalue back into the calculation? (Hint: Rayleigh quotient)
Reset once more.
x = x0/la.norm(x0)
Run this cell in-place (Ctrl-Enter) many times.
for i in range(10):
sigma = np.dot(x, np.dot(A, x))/np.dot(x, x)
x = la.solve(A-sigma*np.eye(n), x)
x = x/la.norm(x)
print(x)