Copyright (C) 2020 Andreas Kloeckner
import numpy as np
import numpy.linalg as la
np.set_printoptions(precision=3, linewidth=120)
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.
x0 = np.random.randn(n)
x = x0
Implement normalized power iteration.
for i in range(10):
x = A @ x
nrm = la.norm(x)
x = x/nrm
print(la.solve(eigvecs, x))
#print(x)
print(nrm)
x = x0
errors = []
coeffs = la.solve(eigvecs, x0)
for i in range(10):
x = A @ x
errors.append(
la.norm(x/eigvals[0]**(i+1) - coeffs[0]*eigvecs[:,0], 2))
print("coefficients:", la.solve(eigvecs, x/la.norm(x,2)))
conv_factor = eigvals[1]/eigvals[0]
errors = np.array(errors)
for i in range(len(errors)-1):
print(f"{i=}: {errors[i]=:.6e}, {errors[i+1]/errors[i]=:.6g}, {conv_factor=:.6g}")
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)