Orthogonal Iteration

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

Let's make a matrix with given eigenvalues:

In [47]:
n = 5

np.random.seed(70)
eigvecs = np.random.randn(n, n)
eigvals = np.sort(np.random.randn(n))

A = np.dot(la.solve(eigvecs, np.diag(eigvals)), eigvecs)
print(eigvals)
[-1.3657822  -0.78460489 -0.08829521  0.30824369  0.52110266]

Let's make an array of iteration vectors:

In [48]:
X = np.random.randn(n, n)

Next, implement orthogonal iteration:

  • Orthogonalize.
  • Apply A
  • Repeat

Run this cell in-place (Ctrl-Enter) many times.

In [90]:
Q, R = la.qr(X)
X = A @ Q
print(Q)
[[-0.36366705 -0.2446216   0.39052837  0.06435026 -0.8070026 ]
 [-0.42943052 -0.63956476  0.08336453 -0.49910706  0.3879289 ]
 [-0.040959   -0.42224059 -0.83670034  0.25124759 -0.23841653]
 [ 0.80079022 -0.35639657  0.08781548 -0.40630442 -0.24273785]
 [-0.20098031  0.47519634 -0.364361   -0.72009898 -0.28721745]]

Now check that the (hopefully) converged $Q$ actually led to Schur form:

In [91]:
la.norm(
    Q @ R @ Q.T
    - A)
Out[91]:
2.4553239646805687e-07

Do the eigenvalues match?

In [92]:
R
Out[92]:
array([[-1.3657822 ,  0.08740833,  3.037871  , -0.26244477, -0.29760524],
       [ 0.        , -0.78460491, -0.45833236, -1.38868012,  0.24515348],
       [ 0.        ,  0.        ,  0.52110264, -0.1520245 ,  0.07378813],
       [ 0.        ,  0.        ,  0.        ,  0.30824369,  0.17685701],
       [ 0.        ,  0.        ,  0.        ,  0.        , -0.08829521]])

What are possible flaws in this plan?

  • Will this always converge?
  • What about complex eigenvalues?