Arnoldi vs Power Iteration

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.

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

Initialization

Set up $Q$ and $H$:

Q = np.zeros((n, n))
H = np.zeros((n, n))

k = 0

Pick a starting vector, normalize it

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:

ritz_values = []
ritz_max = []

Algorithm

Carry out one iteration of Arnoldi iteration.

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

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)

31

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

la.norm(Q[:,:k-1].T @ A @ Q[:,:k-1] - H[:k-1,:k-1])/ la.norm(A)

4.6540588844193635e-11

Check that $AQ-QH$ is fairly small

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

0.044042837433373075

Check that Q is orthogonal:

la.norm((Q.T.conj() @ Q)[:k-1,:k-1] - np.eye(k-1))

4.9488046868226859e-10

Compare max Ritz value to power iteration

#true largest eigenvalue is 50
r = 0
x = x0.copy()
rs = []
for i in range(k):
    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")

49.0975728439 49.9994374247

[<matplotlib.lines.Line2D at 0x7feb5887b668>]

Plot convergence of Ritz values

Enable the Ritz value collection above to make this work.

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