Copyright (C) 2020 Andreas Kloeckner
import numpy as np
import numpy.linalg as la
np.set_printoptions(precision=3)
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])
[-2.668 -0.958 -0.33 -0.292 -0.186 -0.144]
Let's also pick an initial vector:
x0 = np.random.randn(n)
x0
array([ 2.269, 0.664, 0.899, -0.366, 0.463, 0.08 ])
x = x0
Now implement plain power iteration.
for i in range(20):
x = A @ x
print(x)
[ 1.329e+12 -4.740e+12 -7.728e+11 1.018e+12 -8.466e+11 -2.213e+12] [-3.545e+12 1.264e+13 2.062e+12 -2.716e+12 2.258e+12 5.903e+12] [ 9.457e+12 -3.373e+13 -5.500e+12 7.245e+12 -6.025e+12 -1.575e+13] [-2.523e+13 8.998e+13 1.467e+13 -1.933e+13 1.607e+13 4.201e+13] [ 6.730e+13 -2.400e+14 -3.914e+13 5.156e+13 -4.287e+13 -1.121e+14] [-1.795e+14 6.403e+14 1.044e+14 -1.375e+14 1.144e+14 2.989e+14] [ 4.789e+14 -1.708e+15 -2.785e+14 3.669e+14 -3.051e+14 -7.975e+14] [-1.278e+15 4.557e+15 7.430e+14 -9.788e+14 8.139e+14 2.127e+15] [ 3.408e+15 -1.216e+16 -1.982e+15 2.611e+15 -2.171e+15 -5.675e+15] [-9.092e+15 3.243e+16 5.288e+15 -6.965e+15 5.792e+15 1.514e+16] [ 2.425e+16 -8.650e+16 -1.411e+16 1.858e+16 -1.545e+16 -4.039e+16] [-6.470e+16 2.308e+17 3.763e+16 -4.957e+16 4.122e+16 1.077e+17] [ 1.726e+17 -6.156e+17 -1.004e+17 1.322e+17 -1.100e+17 -2.874e+17] [-4.604e+17 1.642e+18 2.678e+17 -3.527e+17 2.933e+17 7.667e+17] [ 1.228e+18 -4.381e+18 -7.143e+17 9.410e+17 -7.825e+17 -2.045e+18] [-3.277e+18 1.169e+19 1.906e+18 -2.510e+18 2.087e+18 5.456e+18] [ 8.741e+18 -3.117e+19 -5.083e+18 6.696e+18 -5.569e+18 -1.455e+19] [-2.332e+19 8.316e+19 1.356e+19 -1.786e+19 1.486e+19 3.883e+19] [ 6.220e+19 -2.219e+20 -3.618e+19 4.765e+19 -3.963e+19 -1.036e+20] [-1.659e+20 5.918e+20 9.650e+19 -1.271e+20 1.057e+20 2.763e+20]
Back to the beginning: Reset to the initial vector.
x = x0/la.norm(x0)
Implement normalized power iteration.
for i in range(10):
x = A @ x
nrm = la.norm(x)
x = x/nrm
print(x)
print(nrm)
[-0.285 0.819 0.106 -0.172 0.149 0.431] [ 0.243 -0.842 -0.13 0.173 -0.156 -0.402] [-0.238 0.845 0.135 -0.177 0.153 0.396] [ 0.237 -0.845 -0.137 0.18 -0.152 -0.395] [-0.237 0.845 0.137 -0.181 0.151 0.395] [ 0.237 -0.845 -0.138 0.181 -0.151 -0.394] [-0.237 0.845 0.138 -0.181 0.151 0.394] [ 0.237 -0.845 -0.138 0.181 -0.151 -0.394] [-0.237 0.845 0.138 -0.181 0.151 0.394] [ 0.237 -0.845 -0.138 0.181 -0.151 -0.394] 2.667662268743778
Extensions:
What if we want the smallest eigenvalue (by magnitude)?
Once again, reset to the beginning.
x = x0/la.norm(x0)
for i in range(30):
x = la.solve(A, x)
nrm = la.norm(x)
x = x/nrm
print(x)
[ 0.344 -0.55 -0.094 0.23 -0.05 -0.718] [-0.527 0.236 0.046 -0.123 -0.066 0.803] [ 0.599 -0.01 -0.037 -0.031 0.127 -0.789] [-0.617 -0.128 0.043 0.173 -0.153 0.74 ] [ 0.611 0.213 -0.052 -0.284 0.161 -0.687] [-0.597 -0.265 0.061 0.365 -0.161 0.64 ] [ 0.583 0.299 -0.067 -0.424 0.16 -0.601] [-0.57 -0.322 0.072 0.465 -0.157 0.571] [ 0.558 0.337 -0.076 -0.496 0.155 -0.547] [-0.55 -0.348 0.078 0.517 -0.153 0.529] [ 0.543 0.356 -0.08 -0.533 0.151 -0.515] [-0.537 -0.362 0.082 0.545 -0.15 0.504] [ 0.533 0.366 -0.083 -0.554 0.149 -0.496] [-0.529 -0.369 0.084 0.561 -0.148 0.49 ] [ 0.527 0.371 -0.084 -0.566 0.147 -0.485] [-0.525 -0.373 0.085 0.57 -0.147 0.481] [ 0.523 0.375 -0.085 -0.573 0.146 -0.478] [-0.522 -0.376 0.086 0.575 -0.146 0.476] [ 0.521 0.376 -0.086 -0.577 0.146 -0.474] [-0.521 -0.377 0.086 0.578 -0.146 0.473] [ 0.52 0.377 -0.086 -0.579 0.146 -0.472] [-0.52 -0.378 0.086 0.58 -0.145 0.471] [ 0.519 0.378 -0.086 -0.581 0.145 -0.471] [-0.519 -0.378 0.086 0.581 -0.145 0.47 ] [ 0.519 0.379 -0.086 -0.582 0.145 -0.47 ] [-0.519 -0.379 0.086 0.582 -0.145 0.47 ] [ 0.519 0.379 -0.086 -0.582 0.145 -0.469] [-0.518 -0.379 0.086 0.582 -0.145 0.469] [ 0.518 0.379 -0.086 -0.582 0.145 -0.469] [-0.518 -0.379 0.086 0.583 -0.145 0.469]
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)
[ 0.011 -0.809 -0.191 0.223 -0.186 -0.474] [-0.127 0.775 0.149 -0.232 0.16 0.531] [ 0.192 -0.728 -0.123 0.246 -0.135 -0.582] [-0.246 0.674 0.103 -0.261 0.108 0.628] [ 0.295 -0.616 -0.085 0.271 -0.08 -0.669] [-0.338 0.556 0.07 -0.277 0.052 0.702] [ 0.376 -0.497 -0.057 0.277 -0.026 -0.728] [-0.41 0.441 0.047 -0.273 0.002 0.749] [ 0.44 -0.387 -0.038 0.264 0.019 -0.765] [-0.465 0.338 0.031 -0.252 -0.039 0.777] [ 0.488 -0.291 -0.026 0.237 0.056 -0.786] [-0.508 0.249 0.022 -0.219 -0.071 0.792] [ 0.525 -0.209 -0.019 0.199 0.084 -0.796] [-0.54 0.172 0.017 -0.178 -0.095 0.798] [ 0.553 -0.138 -0.016 0.155 0.105 -0.799] [-0.565 0.107 0.016 -0.132 -0.114 0.799] [ 0.575 -0.078 -0.016 0.108 0.122 -0.798] [-0.584 0.05 0.016 -0.084 -0.129 0.796] [ 0.591 -0.025 -0.017 0.06 0.135 -0.792] [-0.597 0.001 0.018 -0.035 -0.14 0.789] [ 0.603 0.021 -0.02 0.011 0.144 -0.784] [-0.607 -0.041 0.022 0.013 -0.148 0.779] [ 0.611 0.061 -0.024 -0.036 0.151 -0.774] [-0.614 -0.079 0.025 0.059 -0.154 0.768] [ 0.616 0.096 -0.027 -0.082 0.156 -0.761] [-0.618 -0.112 0.03 0.103 -0.158 0.755] [ 0.619 0.126 -0.032 -0.124 0.16 -0.748] [-0.619 -0.14 0.034 0.145 -0.161 0.741] [ 0.619 0.154 -0.036 -0.164 0.162 -0.734] [-0.619 -0.166 0.038 0.183 -0.163 0.726]
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)
[ 0.063 -0.792 -0.173 0.23 -0.173 -0.505] [-0.191 0.726 0.126 -0.248 0.133 0.585] [ 0.521 -0.254 -0.051 0.121 0.058 -0.802] [-0.544 -0.35 0.081 0.533 -0.15 0.519] [ 0.52 0.378 -0.086 -0.58 0.146 -0.472] [-0.518 -0.379 0.086 0.583 -0.145 0.469] [ 0.518 0.379 -0.086 -0.583 0.145 -0.469] [-0.518 -0.379 0.086 0.583 -0.145 0.469] [-0.518 -0.379 0.086 0.583 -0.145 0.469] [ 0.518 0.379 -0.086 -0.583 0.145 -0.469]