矩阵的特征值（eigenvalue）和特征向量（eigenvector）在很多应用中都具有重要的数学和物理意义。乘幂法（Power Iteration）是线性代数中一种重要的数值计算方法，用于估计矩阵的最大特征值及其对应特征向量的迭代算法，广泛应用于许多科学和工程领域。
本文将详细介绍带有原点移位的乘幂法，并给出其Python实现。
一、乘幂法
- 给定一个方阵
，如果
是
的最大特征值，而
是对应的特征向量
- 那么对于任意非零向量
，迭代
的结果会趋近于
。当
越来越大时，
的方向会趋近于
的方向，而其模（长度）则会趋近于
的绝对值。
【数值计算方法（黄明游）】矩阵特征值与特征向量的计算（一）：乘幂法【理论到程序】
二、乘幂法的加速
1. 天书
2. 人话
乘幂法（Power Iteration）是一种用于寻找矩阵的最大特征值及其对应特征向量的迭代算法。它通过迭代计算矩阵与向量的乘积，并规范化得到新的向量，最终收敛到矩阵的最大特征值和对应的特征向量。然而，对于某些矩阵，乘幂法的收敛速度可能相对较慢。为了加速乘幂法的收敛，一种常见的做法是通过平移（Shift）矩阵的方式。
a. 基本乘幂法
- 选择初始向量
：通常选择一个非零向量作为初始向量，其选择可能影响到迭代的收敛速度。
- 迭代计算
：对于每一次迭代
，计算
。
- 规范化
：通过规范化
（即将其除以其最大分量），得到新的向量
。
- 判断收敛条件：检查
是否小于一个预定的收敛阈值，满足则停止迭代。
- 更新
：将
赋值给
，重复步骤2。
- 计算特征值：一旦迭代收敛，通过
的比值来估计矩阵
的最大特征值。
b. 加速乘幂法的思想
加速乘幂法的思想是通过平移矩阵
，将其转化为
的形式。这个平移操作使得新的矩阵
的最大特征值与原矩阵
的最大特征值之间的差距变大，从而加速收敛。
c. 具体步骤
- 平移操作： 对给定的矩阵
，选择一个数值
，构造新矩阵
。
- 基本乘幂法： 使用
进行上述乘幂法迭代计算，直至收敛
- 计算特征值： 一旦迭代收敛，通过
来估计矩阵
的最大特征值。
3. 典例
4. 实现
a.代码
代码语言：javascript复制import numpy as np
def normalize_vector(vec):
# 计算向量的最大分量
# max_val = np.max(np.abs(vec))
max_index = np.argmax(np.abs(vec))
max = vec[max_index]
# 对向量进行规范化
return vec / vec[max_index], max
def power_iteration(A, x, max_iter=1000, tol=1e-6):
max_val0 = 0
for i in range(max_iter):
# 计算矩阵与向量的乘积
np.set_printoptions(precision=4, suppress=True)
Ax = np.dot(A, x)
# 对乘积结果进行规范化
print(" x{}：".format(i 1), Ax)
x, max_val = normalize_vector(Ax)
print("第{}次：".format(i 1))
print("特征值max(x{})：".format(i 1), max_val)
print("特征向量 y{}：".format(i 1), x)
# 如果最大分量小于收敛阈值，提前退出迭代
if abs(max_val0 - max_val) < tol:
break
max_val0 = max_val
# 计算对应的特征值
eigenvalue = np.dot(Ax, x) / np.dot(x, x)
return eigenvalue
# 给定矩阵 A
A = np.array([[-4, 14, 0],
[-5, 13, 0],
[-1, 0, 2]])
A0 = np.array([[1, 1, 0.5],
[1, 1, 0.25],
[0.5, 0.25, 2]])
# lambda0 = 0.75
# B = A0-lambda0*np.eye(3)
A3 = np.array([[-11, 11, 1],
[11, 9, -2],
[1, -2, 13]])
lambda0 = 15
B = A3 - lambda0 * np.eye(3)
# 初始向量
x0 = np.array([1, 1, 1])
# 运行乘幂法迭代
# result_vector, eigenvalue = power_iteration(A, x0)
# result_vector, eigenvalue = power_iteration(A0, x0)
# result_vector, eigenvalue = power_iteration(A3, x0)
eigenvalue0 = power_iteration(B, x0)
eigenvalue = lambda0 eigenvalue0
print("估计的特征值:", eigenvalue)
# np.dot(A, result_vector)
# np.dot(eigenvalue, result_vector)b.解析
-
normalize_vector函数：- 输入：一个向量
vec。 - 功能：计算向量的最大分量，并将向量规范化。
- 输出：规范化后的向量和最大分量。
- 输入：一个向量
-
power_iteration函数：- 输入：
A：一个方阵（矩阵）。x：初始向量。max_iter：最大迭代次数，默认为1000。tol：收敛阈值，默认为1e-4。
- 功能：使用乘幂法迭代来估计矩阵的最大特征值及其对应的特征向量。
- 计算矩阵 A 与向量 x 的乘积，得到 Ax。
- 调用
normalize_vector函数对 Ax 进行规范化，得到规范化后的向量和最大分量。 - 打印每次迭代的结果，即特征值、特征向量。
- 判断是否满足收敛条件，如果最大分量变化小于阈值
tol，则提前退出迭代。 - 计算对应的特征值，更新最大分量，并继续迭代。
- 输出：估计得到的特征值。
- 输入：
- 主程序部分：
- 教材例题及课后题的矩阵
A、A1、A2、A3。 - 定义了初始向量
x0、系数
lambda - 构造新矩阵
B = A - lambda I 。
- 调用
power_iteration函数，传入矩阵
B 和初始向量进行乘幂法迭代。
B 的特征值
lambda =
A 的最大特征值。
- 教材例题及课后题的矩阵
c.应用
代码语言：javascript复制 x1： [1.75 1.5 2. ]
第1次：
特征值max(x1)： 2.0
特征向量 y1： [0.875 0.75 1. ]
x2： [1.4688 1.3125 1.875 ]
第2次：
特征值max(x2)： 1.875
特征向量 y2： [0.7833 0.7 1. ]
x3： [1.3958 1.2083 1.8167]
第3次：
特征值max(x3)： 1.8166666666666667
特征向量 y3： [0.7683 0.6651 1. ]
x4： [1.3572 1.1846 1.8005]
第4次：
特征值max(x4)： 1.8004587155963303
特征向量 y4： [0.7538 0.658 1. ]
x5： [1.3464 1.1683 1.7914]
第5次：
特征值max(x5)： 1.7914012738853502
特征向量 y5： [0.7516 0.6522 1. ]
x6： [1.3401 1.1646 1.7888]
第6次：
特征值max(x6)： 1.7888444444444445
特征向量 y6： [0.7491 0.6511 1. ]
x7： [1.3383 1.1619 1.7873]
第7次：
特征值max(x7)： 1.7873301200029814
特征向量 y7： [0.7488 0.6501 1. ]
x8： [1.3373 1.1613 1.7869]
第8次：
特征值max(x8)： 1.7869153405872398
特征向量 y8： [0.7484 0.6499 1. ]
x9： [1.337 1.1608 1.7867]
第9次：
特征值max(x9)： 1.7866588534107755
特征向量 y9： [0.7483 0.6497 1. ]
x10： [1.3368 1.1608 1.7866]
第10次：
特征值max(x10)： 1.786591427357151
特征向量 y10： [0.7482 0.6497 1. ]
x11： [1.3368 1.1607 1.7865]
第11次：
特征值max(x11)： 1.7865478412280433
特征向量 y11： [0.7482 0.6497 1. ]
x12： [1.3367 1.1607 1.7865]
第12次：
特征值max(x12)： 1.7865369093758368
特征向量 y12： [0.7482 0.6497 1. ]
x13： [1.3367 1.1606 1.7865]
第13次：
特征值max(x13)： 1.786529489372963
特征向量 y13： [0.7482 0.6497 1. ]
x14： [1.3367 1.1606 1.7865]
第14次：
特征值max(x14)： 1.7865277239374533
特征向量 y14： [0.7482 0.6497 1. ]
x15： [1.3367 1.1606 1.7865]
第15次：
特征值max(x15)： 1.7865264587497982
特征向量 y15： [0.7482 0.6497 1. ]
x16： [1.3367 1.1606 1.7865]
第16次：
特征值max(x16)： 1.786526175001885
特征向量 y16： [0.7482 0.6497 1. ]
估计的特征值: 2.5365264587497985. 课后题
a. 乘幂法
代码语言：javascript复制 x1： [ 1 18 12]
第1次：
特征值max(x1)： 18
特征向量 y1： [0.0556 1. 0.6667]
x2： [11.0556 8.2778 6.7222]
第2次：
特征值max(x2)： 11.055555555555555
特征向量 y2： [1. 0.7487 0.608 ]
x3： [-2.1558 16.5226 7.407 ]
第3次：
特征值max(x3)： 16.52261306532663
特征向量 y3： [-0.1305 1. 0.4483]
x4： [12.8835 6.6682 3.6974]
第4次：
特征值max(x4)： 12.88351581508516
特征向量 y4： [1. 0.5176 0.287 ]
x5： [-5.0197 15.0842 3.6957]
第5次：
特征值max(x5)： 15.084204811029013
特征向量 y5： [-0.3328 1. 0.245 ]
x6： [14.9056 4.8494 0.8523]
第6次：
特征值max(x6)： 14.90555759004166
特征向量 y6： [1. 0.3253 0.0572]
x7： [-7.364 13.8137 1.0926]
第7次：
特征值max(x7)： 13.813746465256594
特征向量 y7： [-0.5331 1. 0.0791]
x8： [16.9431 2.9778 -1.5049]
第8次：
特征值max(x8)： 16.94313702089086
特征向量 y8： [ 1. 0.1758 -0.0888]
x9： [-9.1556 12.7594 -0.5061]
第9次：
特征值max(x9)： 12.759391141602638
特征向量 y9： [-0.7176 1. -0.0397]
x10： [18.8534 1.1862 -3.2332]
第10次：
特征值max(x10)： 18.853433942481207
特征向量 y10： [ 1. 0.0629 -0.1715]
x11： [-10.4794 11.9093 -1.3552]
第11次：
特征值max(x11)： 11.909254019299365
特征向量 y11： [-0.8799 1. -0.1138]
x12： [20.5655 -0.4517 -4.3593]
第12次：
特征值max(x12)： 20.565503848485033
特征向量 y12： [ 1. -0.022 -0.212]
x13： [-11.4536 11.2263 -1.7117]
第13次：
特征值max(x13)： -11.453578681607251
特征向量 y13： [ 1. -0.9802 0.1494]
x14： [-21.6322 1.8797 4.9031]
第14次：
特征值max(x14)： -21.632241268767345
特征向量 y14： [ 1. -0.0869 -0.2267]
x15： [-12.1825 10.6713 -1.7728]
第15次：
特征值max(x15)： -12.182498255541232
特征向量 y15： [ 1. -0.876 0.1455]
x16： [-20.4899 2.8254 4.6436]
第16次：
特征值max(x16)： -20.48993757175785
特征向量 y16： [ 1. -0.1379 -0.2266]
x17： [-12.7434 10.2122 -1.6704]
第17次：
特征值max(x17)： -12.743448525598359
特征向量 y17： [ 1. -0.8014 0.1311]
x18： [-19.684 3.5255 4.3068]
第18次：
特征值max(x18)： -19.683997941414024
特征向量 y18： [ 1. -0.1791 -0.2188]
x19： [-13.189 9.8256 -1.4861]
第19次：
特征值max(x19)： -13.1889525203576
特征向量 y19： [ 1. -0.745 0.1127]
x20： [-19.0822 4.0697 3.9548]
第20次：
特征值max(x20)： -19.08221542867637
特征向量 y20： [ 1. -0.2133 -0.2073]
x21： [-13.5533 9.495 -1.2677]………………
代码语言：javascript复制第148次：
特征值max(x148)： -15.96947870451083
特征向量 y148： [ 1. -0.4457 -0.0653]
x149： [-15.9683 7.1192 1.0421]
第149次：
特征值max(x149)： -15.96828831423285
特征向量 y149： [ 1. -0.4458 -0.0653]
x150： [-15.9694 7.118 1.0433]
第150次：
特征值max(x150)： -15.969403255924263
特征向量 y150： [ 1. -0.4457 -0.0653]
x151： [-15.9684 7.1191 1.0422]
第151次：
特征值max(x151)： -15.96835897584962
特征向量 y151： [ 1. -0.4458 -0.0653]
x152： [-15.9693 7.1181 1.0432]
第152次：
特征值max(x152)： -15.969337068297888
特征向量 y152： [ 1. -0.4457 -0.0653]
估计的特征向量: [ 1. -0.4457 -0.0653]
估计的特征值: -15.96835897584962b. 带原点移位的乘幂法
代码语言：javascript复制 x1： [-14. 3. -3.]
第1次：
特征值max(x1)： -14.0
特征向量 y1： [ 1. -0.2143 0.2143]
x2： [-28.1429 11.8571 1. ]
第2次：
特征值max(x2)： -28.142857142857142
特征向量 y2： [ 1. -0.4213 -0.0355]
x3： [-30.6701 13.599 1.9137]
第3次：
特征值max(x3)： -30.67005076142132
特征向量 y3： [ 1. -0.4434 -0.0624]
x4： [-30.9398 13.7852 2.0116]
第4次：
特征值max(x4)： -30.93975504799735
特征向量 y4： [ 1. -0.4455 -0.065 ]
x5： [-30.9661 13.8033 2.0211]
第5次：
特征值max(x5)： -30.966052915940043
特征向量 y5： [ 1. -0.4458 -0.0653]
x6： [-30.9686 13.8051 2.0221]
第6次：
特征值max(x6)： -30.968592776449142
特征向量 y6： [ 1. -0.4458 -0.0653]
x7： [-30.9688 13.8052 2.0221]
第7次：
特征值max(x7)： -30.968837849183018
特征向量 y7： [ 1. -0.4458 -0.0653]
x8： [-30.9689 13.8053 2.0221]
第8次：
特征值max(x8)： -30.968861494278652
特征向量 y8： [ 1. -0.4458 -0.0653]
x9： [-30.9689 13.8053 2.0222]
第9次：
特征值max(x9)： -30.968863775583827
特征向量 y9： [ 1. -0.4458 -0.0653]
x10： [-30.9689 13.8053 2.0222]
第10次：
特征值max(x10)： -30.968863995686508
特征向量 y10： [ 1. -0.4458 -0.0653]
估计的特征值: -15.96886377558382
