Power Method for Eigenvalues
幂法：矩阵本征值迭代求解的基本方法¶
基本的幂法：求矩阵的最大本征值¶
$$ Ax^{(k)} = x^{(k+1)} $$
只要把任何一个初始随机向量反复地被一个矩阵作用，同时除以一个足够大的数保证不超过数值精度，那么最终的收敛结果就是矩阵的最大本征矢。
这种方法可行的原因是，矩阵可以看做是一个线性变换，对任何一个向量连续地作用同一个变换的结果，总是会让这个向量趋近于变换的主向量，也就是矩阵的最大本征值对应的本征矢。
In [1]:
import numpy as np
In [2]:
def power_iteration(A, num_simulations: int, v0=None):
if v0 is None:
v = np.random.rand(A.shape[1])
else:
v = v0
for _ in range(num_simulations):
v = A @ v / np.max(v)
return v
In [3]:
a = np.random.rand(4,4)
e2,v2=np.linalg.eig(a)
ei2 = e2.argmax()
e2[ei2], v2[:,ei2] / np.linalg.norm(v2[:,ei2])
Out[3]:
In [4]:
v1 = power_iteration(a, 1000)
(a @ v1 / v1)[0], v1/ np.linalg.norm(v1)
Out[4]:
反幂法：求矩阵的最小本征矢¶
与幂法完全相同，只不过迭代的是矩阵的逆，显然矩阵的逆的最大本征值就是原矩阵的模最小的本征值。 $$ A^{-1}x^{(k)} = x^{(k+1)} $$ 不过由于矩阵求逆计算复杂度很大，可以把它转化成求线性方程组： $$ A x^{(k+1)} = x^{(k)} $$
In [5]:
def inverse_power_iteration(A, num_simulations: int, v0=None):
if v0 is None:
v = np.random.rand(A.shape[1])
else:
v = v0
for _ in range(num_simulations):
v = np.linalg.solve(A, v)
# v = np.linalg.inv(A) @ v
v = v / np.max(v)
return v
In [6]:
a = np.random.rand(4,4)
e2,v2=np.linalg.eig(a)
ei2 = np.abs(e2).argmin()
e2[ei2], v2[:,ei2] / np.linalg.norm(v2[:,ei2])
Out[6]:
In [7]:
v1 = inverse_power_iteration(a, 1000)
1/ (np.linalg.inv(a) @ v1 / v1)[0], v1/ np.linalg.norm(v1)
Out[7]:
In [ ]:
评论
Comments powered by Disqus