Copyright (C) 2020 Andreas Kloeckner
import numpy as np
import numpy.linalg as la
Let's make a matrix with given eigenvalues:
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:
X = np.random.randn(n, n)
Next, implement orthogonal iteration:
Run this cell in-place (Ctrl-Enter) many times.
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:
la.norm(
Q @ R @ Q.T
- A)
2.4553239646805687e-07
Do the eigenvalues match?
R
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?