【数值计算方法（黄明游）】矩阵特征值与特征向量的计算（一）：乘幂法【理论到程序】

A_1

A_2

矩阵的特征值（eigenvalue）和特征向量（eigenvector）在很多应用中都具有重要的数学和物理意义。乘幂法（Power Iteration）是线性代数中一种重要的数值计算方法，用于估计矩阵的最大特征值及其对应特征向量的迭代算法，广泛应用于许多科学和工程领域。

本文将详细介绍乘幂法的基本原理和步骤，并给出其Python实现。

一、乘幂法

1. 天书

a. 乘幂法

本文仅考虑有唯一的主特征值情况，的主特征值不唯一情况不做介绍

b. 理论证明

c. 规范化的乘幂法

2. 人话

a. 数学原理

乘幂法（Power Iteration）是一种用于估计矩阵的最大特征值及其对应特征向量的迭代算法，基于以下的数学原理：

给定一个方阵

，如果

\lambda

是

的最大特征值，而

是对应的特征向量

那么对于任意非零向量

，迭代

A^n v

的结果会趋近于

\lambda^n x

。当

越来越大时，

A^n v

的方向会趋近于

的方向，而其模（长度）则会趋近于

\lambda^n

的绝对值。

理论证明详见上述天书部分

b. 基本步骤

选择初始向量

x_0

：通常选择一个非零向量作为初始向量，其选择可能影响到迭代的收敛速度。

迭代计算

A x_k

：对于每一次迭代

，计算

A x_k

。

规范化

A x_k

：通过规范化

A x_k

（即将其除以其最大分量），得到新的向量

x_{k+1}

。

判断收敛条件：检查

\frac{\|x_{k+1} - x_k\|}{\|x_{k+1}\|}

是否小于一个预定的收敛阈值，满足则停止迭代。

更新

x_k

：将

x_{k+1}

赋值给

x_k

，重复步骤2。

计算特征值：一旦迭代收敛，通过

\frac{A x_k}{x_k}

的比值来估计矩阵

的最大特征值。

乘幂法的优点是它的简单性和易实现性。然而，它只能找到最大特征值，且其收敛速度可能相对较慢。对于一些特殊的矩阵，可能需要使用其他的迭代方法。

c. 注意事项

收敛性：
- 乘幂法只能估计最大特征值，并且其收敛速度取决于初始向量的选择以及特征值之间的差异。
- 对于对称正定矩阵，收敛是保证的。
复杂性：
- 乘幂法是一种简单且易于实现的方法，但对于某些情况下的矩阵，收敛速度可能较慢。
- 在某些情况下，可能需要使用其他迭代方法。
对称矩阵： 乘幂法在处理对称矩阵时效果更好，因为对称矩阵的特征向量是正交的。
扩展： 乘幂法的扩展形式包括反幂法、带有原点移位的乘幂法等。

3. 典例

4. 实现

代码语言：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-4):

max_val0 = 0

for i in range(max_iter):

# 计算矩阵与向量的乘积

np.set_printoptions(precision=3, 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 x, eigenvalue

给定矩阵 A
A = np.array([[-4, 14, 0],

[-5, 13, 0],

[-1, 0, 2]])
A1 = np.array([[4, 2, 2],

[2, 5, 1],

[2, 1, 6]])
A2 = np.array([[7, 3, -2],

[3, 4, -1],

[-2, -1, 3]])

A3 = np.array([[-11, 11, 1],

[11, 9, -2],

[1, -2, 13]])
初始向量
x0 = np.array([1, 1, 1])
运行乘幂法迭代
result_vector, eigenvalue = power_iteration(A, x0)
result_vector, eigenvalue = power_iteration(A1, x0)
result_vector, eigenvalue = power_iteration(A2, x0)
result_vector, eigenvalue = power_iteration(A3, x0)
print("估计的特征向量:", result_vector)

print("估计的特征值:", eigenvalue)

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。
- 调用 power_iteration 函数，分别传入不同的矩阵和初始向量进行乘幂法迭代。
- 打印估计得到的特征向量和特征值。

代码语言：javascript

复制

          x1： [10  8  1]
第1次：
特征值max(x1)： 10
特征向量   y1： [1.  0.8 0.1]
          x2： [ 7.2  5.4 -0.8]
第2次：
特征值max(x2)： 7.200000000000001
特征向量   y2： [ 1.      0.75   -0.1111]
          x3： [ 6.5     4.75   -1.2222]
第3次：
特征值max(x3)： 6.499999999999998
特征向量   y3： [ 1.      0.7308 -0.188 ]
          x4： [ 6.2308  4.5    -1.3761]
第4次：
特征值max(x4)： 6.23076923076923
特征向量   y4： [ 1.      0.7222 -0.2209]
          x5： [ 6.1111  4.3889 -1.4417]
第5次：
特征值max(x5)： 6.1111111111111125
特征向量   y5： [ 1.      0.7182 -0.2359]
          x6： [ 6.0545  4.3364 -1.4718]
第6次：
特征值max(x6)： 6.054545454545453
特征向量   y6： [ 1.      0.7162 -0.2431]
          x7： [ 6.027   4.3108 -1.4862]
第7次：
特征值max(x7)： 6.0270270270270245
特征向量   y7： [ 1.      0.7152 -0.2466]
          x8： [ 6.0135  4.2982 -1.4932]
第8次：
特征值max(x8)： 6.013452914798206
特征向量   y8： [ 1.      0.7148 -0.2483]
          x9： [ 6.0067  4.2919 -1.4966]
第9次：
特征值max(x9)： 6.006711409395972
特征向量   y9： [ 1.      0.7145 -0.2492]
          x10： [ 6.0034  4.2888 -1.4983]
第10次：
特征值max(x10)： 6.003351955307263
特征向量   y10： [ 1.      0.7144 -0.2496]
          x11： [ 6.0017  4.2873 -1.4992]
第11次：
特征值max(x11)： 6.001675041876046
特征向量   y11： [ 1.      0.7143 -0.2498]
          x12： [ 6.0008  4.2865 -1.4996]
第12次：
特征值max(x12)： 6.00083728718951
特征向量   y12： [ 1.      0.7143 -0.2499]
估计的特征向量: [ 1.      0.7143 -0.2499]
估计的特征值: 6.001675041876047

5. 课后题

黄明游~《数值计算方法》 ~ P91

A_1

代码语言：javascript

复制

          x1： [8 8 9]
第1次：
特征值max(x1)： 9
特征向量   y1： [0.8889 0.8889 1.    ]
          x2： [7.3333 7.2222 8.6667]
第2次：
特征值max(x2)： 8.666666666666666
特征向量   y2： [0.8462 0.8333 1.    ]
          x3： [7.0513 6.859  8.5256]
第3次：
特征值max(x3)： 8.525641025641026
特征向量   y3： [0.8271 0.8045 1.    ]
          x4： [6.9173 6.6767 8.4586]
第4次：
特征值max(x4)： 8.458646616541353
特征向量   y4： [0.8178 0.7893 1.    ]
          x5： [6.8498 6.5822 8.4249]
第5次：
特征值max(x5)： 8.424888888888889
特征向量   y5： [0.813  0.7813 1.    ]
          x6： [6.8147 6.5325 8.4074]
第6次：
特征值max(x6)： 8.407364422874025
特征向量   y6： [0.8106 0.777  1.    ]
          x7： [6.7963 6.5061 8.3981]
第7次：
特征值max(x7)： 8.398130137416075
特征向量   y7： [0.8093 0.7747 1.    ]
          x8： [6.7865 6.4921 8.3932]
第8次：
特征值max(x8)： 8.39322778520782
特征向量   y8： [0.8086 0.7735 1.    ]
          x9： [6.7812 6.4846 8.3906]
第9次：
特征值max(x9)： 8.390615458295574
特征向量   y9： [0.8082 0.7728 1.    ]
          x10： [6.7784 6.4806 8.3892]
第10次：
特征值max(x10)： 8.389220813597767
特征向量   y10： [0.808  0.7725 1.    ]
          x11： [6.777  6.4784 8.3885]
第11次：
特征值max(x11)： 8.388475552357708
特征向量   y11： [0.8079 0.7723 1.    ]
估计的特征向量: [0.8079 0.7723 1.    ]
估计的特征值: 8.389220813597769

A_2

代码语言：javascript

复制

          x1： [8 6 0]
第1次：
特征值max(x1)： 8
特征向量   y1： [1.   0.75 0.  ]
          x2： [ 9.25  6.   -2.75]
第2次：
特征值max(x2)： 9.25
特征向量   y2： [ 1.      0.6486 -0.2973]
          x3： [ 9.5405  5.8919 -3.5405]
第3次：
特征值max(x3)： 9.54054054054054
特征向量   y3： [ 1.      0.6176 -0.3711]
          x4： [ 9.5949  5.8414 -3.7309]
第4次：
特征值max(x4)： 9.594900849858357
特征向量   y4： [ 1.      0.6088 -0.3888]
          x5： [ 9.6041  5.824  -3.7753]
第5次：
特征值max(x5)： 9.604074402125775
特征向量   y5： [ 1.      0.6064 -0.3931]
          x6： [ 9.6054  5.8187 -3.7857]
第6次：
特征值max(x6)： 9.605429001813766
特征向量   y6： [ 1.      0.6058 -0.3941]
          x7： [ 9.6056  5.8172 -3.7881]
第7次：
特征值max(x7)： 9.60557200236834
特征向量   y7： [ 1.      0.6056 -0.3944]
估计的特征向量: [ 1.      0.6056 -0.3944]
估计的特征值: 9.605429001813766