Arnoldi vs Power Iteration

In [1]:
import numpy as np
import numpy.linalg as la

import matplotlib.pyplot as pt

Let us make a matrix with a defined set of eigenvalues and eigenvectors, given by eigvals and eigvecs.

In [2]:
np.random.seed(40)

# Generate matrix with eigenvalues 1...50
n = 50
eigvals = np.linspace(1., n, n)
eigvecs = np.random.randn(n, n)
#To work with symmetric matrix, orthogonalize eigvecs
eigvecs, R = la.qr(eigvecs)

print(eigvals)

A = la.solve(eigvecs, np.dot(np.diag(eigvals), eigvecs))
print(la.eig(A)[0])

[ 1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11. 12. 13. 14. 15. 16. 17. 18.
 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36.
 37. 38. 39. 40. 41. 42. 43. 44. 45. 46. 47. 48. 49. 50.]
[50.  1. 49.  2. 48.  3.  4. 47.  5. 46. 45. 44.  6.  7.  8. 43.  9. 42.
 10. 41. 40. 11. 39. 12. 13. 14. 38. 37. 15. 36. 16. 35. 34. 17. 33. 32.
 18. 31. 30. 29. 19. 28. 20. 21. 27. 22. 26. 25. 24. 23.]

Initialization

Set up $Q$ and $H$:

In [3]:
Q = np.zeros((n, n))
H = np.zeros((n, n))

k = 0

Pick a starting vector, normalize it

In [4]:
x0 = np.random.randn(n)
x0 = x0/la.norm(x0)

# Poke it into the first column of Q
Q[:, k] = x0.copy()

Make a list to save arrays of Ritz values:

In [5]:
ritz_values = []
ritz_max = []

Algorithm

Carry out one iteration of Arnoldi iteration.

Run this cell in-place (Ctrl-Enter) until H is filled.

In [26]:
print(k)

u = A @ Q[:, k]

# Carry out Gram-Schmidt on u against Q
# to do Lanczos change range start to k-1
for j in range(0,k+1):
    qj = Q[:, j]
    H[j,k] = qj @ u
    u = u - H[j,k]*qj

if k+1 < n:
    H[k+1, k] = la.norm(u)
    Q[:, k+1] = u/H[k+1, k]

k += 1

pt.spy(H)


if k>1:
    D = la.eig(H)[0]
    max_ritz = D[np.argmax(np.abs(D))]
    ritz_vals = np.zeros(k)
    for i in range(k):
        ritz_vals[i] = D[np.argmax(np.abs(D))]
        D[np.argmax(np.abs(D))] = 0
    ritz_max.append(max_ritz)
    ritz_values.append(ritz_vals)

19

Check that $Q^T A Q =H$:

In [27]:
la.norm(Q[:,:k-1].T @ A @ Q[:,:k-1] - H[:k-1,:k-1])/ la.norm(A)
Out[27]:
2.748225022299607e-14

Check that $AQ-QH$ is fairly small

In [28]:
la.norm(A @ Q[:,:k-1] - Q[:,:k-1]@H[:k-1,:k-1])/ la.norm(A)
Out[28]:
0.05635991594049524

Check that Q is orthogonal:

In [29]:
la.norm((Q.T.conj() @ Q)[:k-1,:k-1] - np.eye(k-1))
Out[29]:
2.437715880278286e-13

Compare Max Ritz Value to Power Iteration

In [30]:
#true largest eigenvalue is 50
r = 0
x = A @ x0.copy()
x = x / la.norm(x)
rs = []
for i in range(k-1):
    y = A @ x
    r = x @ y
    x = y / la.norm(y)
    rs.append(r)
print(r,max_ritz)
pt.plot(rs, "x")
pt.plot(ritz_max, "o")

48.97352938503617 49.99902032411266
Out[30]:
[<matplotlib.lines.Line2D at 0x7f140bf9fbe0>]

Plot convergence of Ritz values

Enable the Ritz value collection above to make this work.

In [31]:
for i, rv in enumerate(ritz_values):
    pt.plot([i] * len(rv), rv, "x")