Logo

特征值的迭代估计

1. 特征值近似的背景

在科学计算和工程应用中，精确地计算出所有特征值往往既有难度且没有必要。例如，对于 n×n ~n\times n~的矩阵 A ~\mathbf{A}~，当 n>5 ~n \gt 5~时，特征方程 det⁡(A−λI)=0 ~\det(\mathbf{A} - \lambda\mathbf{I}) = 0~没有通用的解析解。再例如前面介绍过的动力系统，我们通常只需要确定最大的特征值（即主特征值）即可，而无需求解所有特征值。基于这些原因，获取主特征值的近似值往往已经足够满足需求。为了高效地估算最大特征值，下面介绍一种经典的迭代方法：幂法(Power Method)(\textbf{Power Method})。

2. 幂法：主特征值估算

如果一个矩阵 A ~\mathbf{A}~反复对任意一个初始向量 x0 ~\mathbf{x}_0~做乘法，那么经过足够多的次数后，最终得到的向量yk=Akx0 ( k=1,2,3⋯ )\mathbf{y}_k = \mathbf{A}^k\mathbf{x}_0\,(\,k=1,2,3\cdots)的方向会无限趋近于主特征向量的方向。我们设计这样一个实验，取矩阵 A ~\mathbf{A}~和初始向量 x0 ~\mathbf{x}_0~如下：
A=[1.80.80.21.2],x0=[−0.51]\mathbf{A} = \begin{bmatrix}1.8 & 0.8 \\ 0.2 & 1.2\end{bmatrix},\quad \mathbf{x}_0 = \begin{bmatrix}-0.5 \\ 1\end{bmatrix}
矩阵 A ~\mathbf{A}~的特征值 λ1=2~\lambda_1=2（主特征值）、 λ2=1  ~\lambda_2=1~~，对应的特征向量为v1=[41]T, v2=[−11]T\mathbf{v}_1 = \begin{bmatrix}4 & 1\end{bmatrix}^T,\,\mathbf{v}_2=\begin{bmatrix}-1 & 1\end{bmatrix}^T。下面绘制出每次迭代后向量的方向， 请观察这个过程：

开通会员解锁全部动画

随着迭代次数的增加，新的向量yk\mathbf{y}_k与主特征向量 v1 ~\mathbf{v}_1~之间的夹角 θ ~\theta~开始趋近于 0 ~0~。 我们在迭代到第 10 ~10~次后停止，此时向量 y10=[408.7103.3]T ~\mathbf{y}_{10} = \begin{bmatrix}408.7 & 103.3\end{bmatrix}^T~，可近似为特征向量。特征值可以通过向量 y10 ~\mathbf{y}_{10}~和 y9 ~\mathbf{y}_9~中最大分量的比值来近似：
λk=(y10)1(y9)1=408.7203.9≈λ1=2\lambda_k = \frac{(\mathbf{y}_{10})_{1}}{(\mathbf{y}_{9})_{1}} = \frac{408.7}{203.9} \approx \lambda_1 = 2
幂法正是基于这一现象，用于估算矩阵 A ~\mathbf{A}~的主特征值 λ1 ~\lambda_1~及其对应的特征向量 v1 ~\mathbf{v}_1~。其基本思想是：从一个初始向量 x0 ~\mathbf{x}_0~出发，不断地对其进行矩阵幂乘操作 Akx0 ~\mathbf{A}^k\mathbf{x}_0~，利用主特征值的绝对值显著大于其他特征值的特性，使向量逐步沿着主导特征向量的方向收敛。

3. 幂法的数学基础

假设 A ~\mathbf{A}~是一个 n×n ~n\times n~矩阵，其特征分解为：
A=VΛV−1\mathbf{A} = \mathbf{V} \Lambda \mathbf{V}^{-1}
其中 V\mathbf{V} 是特征向量矩阵，Λ\Lambda 是包含特征值λ1,λ2,…,λn\lambda_1, \lambda_2, \dots, \lambda_n的对角矩阵，并满足∣λ1∣>∣λ2∣≥⋯≥∣λn∣|\lambda_1| > |\lambda_2| \geq \cdots \geq |\lambda_n|。 对于一个初始向量 x0 ~\mathbf{x}_0~，可以表示为特征向量的线性组合：
x0=c1v1+c2v2+⋯+cnvn\mathbf{x}_0 = c_1 \mathbf{v}_1 + c_2 \mathbf{v}_2 + \cdots + c_n \mathbf{v}_n
对 x0 ~\mathbf{x}_0~反复施加矩阵 A ~\mathbf{A}~后：
Akx0=c1λ1kv1+c2λ2kv2+⋯+cnλnkvn\mathbf{A}^k \mathbf{x}_0 = c_1 \lambda_1^k \mathbf{v}_1 + c_2 \lambda_2^k \mathbf{v}_2 + \cdots + c_n \lambda_n^k \mathbf{v}_n
由于∣λ1∣≫∣λ2∣,∣λ3∣,…|\lambda_1| \gg |\lambda_2|, |\lambda_3|, \dots，当 k→∞ ~k \to \infty~时：
Akx0≈c1λ1kv1\mathbf{A}^k \mathbf{x}_0 \approx c_1 \lambda_1^k \mathbf{v}_1
即， Akx0 ~\mathbf{A}^k \mathbf{x}_0~ 的方向趋于主导特征向量 v1 ~\mathbf{v}_1~，而其模长以 ∣λ1∣k|\lambda_1|^k的速率增长。

4. 算法步骤

在前面的示例中，第 10 ~10~次迭代停止后， 向量的最大分量 (y10)1=408.7 ~(\mathbf{y}_{10})_1 = 408.7~和 (y9)1=203.9 ~(\mathbf{y}_9)_1=203.9~已经很大了，如果我们再继续迭代下去，向量的分量值就会迅速膨胀、甚至溢出。其实，我们只需要关心相邻两个向量之间的比值就可以近似得到主特征值 λ ~\lambda~以及主特征向量 v ~\mathbf{v}~。所以，幂法在每次迭代得到新的向量 yk ~\mathbf{y}_k~后，会通过对新向量 yk ~\mathbf{y}_k~进行归一化处理来防止数值溢出。具体步骤如下：

初始向量 x0 ~\mathbf{x}_0~可以是任意的非零向量，但为了方便计算， x0 ~\mathbf{x}_0~的最大分量建议不超过 1 ~1~。

对于 k=0,1,2,⋯ ~k=0,1,2,\cdots~：
  1. 计算矩阵乘积：
    yk=Axk\mathbf{y}_k = \mathbf{A} \mathbf{x}_k
  2. 找到 yk ~\mathbf{y}_k~的最大分量（绝对值最大元素）：
    λk=max⁡(∣yk∣)\lambda_k = \max(|\mathbf{y}_k|)
  3. 将 yk ~\mathbf{y}_k~归一化：
    xk+1=ykλk\mathbf{x}_{k+1} = \frac{\mathbf{y}_k}{\lambda_k}

  • 当两次迭代之间的特征值估算 ∣λk+1−λk∣∣λk∣<ϵ ~ \frac{|\lambda_{k+1} - \lambda_{k}|}{|\lambda_{k}|} < \epsilon~的相对变化量足够小时（也可以是连续两个向量间的夹角 cos⁡−1(vk⋅vk−1∥vk∥∥vk−1∥)<ϵ ~\cos^{-1} \left( \frac{\mathbf{v}_k \cdot \mathbf{v}_{k-1}}{\|\mathbf{v}_k\| \|\mathbf{v}_{k-1}\|} \right) < \epsilon~），算法终止。

  • 输出 λk ~\lambda_k~和对应的特征向量 xk ~\mathbf{x}_k~。

我们这次用幂法的步骤重新迭代计算 A ~\mathbf{A}~的主特征值，初始条件以及终止条件如下：
A=[1.80.80.21.2],x0=[−0.51],ϵ=10−3\mathbf{A}=\begin{bmatrix}1.8 & 0.8 \\ 0.2 & 1.2\end{bmatrix},\quad\mathbf{x}_0=\begin{bmatrix}-0.5\\ 1\end{bmatrix},\quad \epsilon=10^{-3}
借助 python ~\text{python}~代码，可得幂法迭代的结果如下：
迭代次数 k vk~\mathbf{v}_k归一化后 vk~\mathbf{v}_kε（弧度）
0(-0.5, 1)(-0.5, 1)N/A
1(-0.1, 1.1)(-0.09, 1)0.373
2(0.64, 1.18)(0.54, 1)0.5846
3(1.77, 1.31)(1, 0.74)0.4403
4(2.39, 1.09)(1, 0.45)0.2099
5(2.16, 0.75)(1, 0.34)0.0948
6(2.08, 0.61)(1, 0.3)0.0444
7(2.04, 0.55)(1, 0.27)0.0215
8(2.02, 0.53)(1, 0.26)0.0105
9(2.01, 0.51)(1, 0.26)0.0052
10(2, 0.51)(1, 0.25)0.0026
11(2, 0.5)(1, 0.25)0.0013
12(2, 0.5)(1, 0.25)0.0006

主特征值收敛于 2 ~2~，主特征向量收敛于[10.25]T\begin{bmatrix}1 & 0.25\end{bmatrix}^T，其渐进轨迹如下：

开通会员解锁全部动画

5. 幂法的收敛性

幂法的收敛速度决定了算法的运算效率，而收敛速度主要和下面三个因素有关：
  1. 初始向量 x0 ~\mathbf{x}_0~的选择：如果初始向量在主特征向量方向上的分量 c1 ~c_1~为 0 ~0~，理论上算法将无法收敛到主特征向量。例如用前面的示例中的矩阵 A ~\mathbf{A}~，以及初始向量 x0′ ~\mathbf{x}'_0~：
    A=[1.80.80.21.2],x0′=[−11]\mathbf{A} = \begin{bmatrix}1.8 & 0.8 \\ 0.2 & 1.2\end{bmatrix},\quad \mathbf{x}'_0 = \begin{bmatrix}-1 \\ 1\end{bmatrix}
    由于这里特意选择了 A ~\mathbf{A}~的非主特征向量 v2 ~\mathbf{v}_2~作为初始向量，幂法的迭代过程满足：Akx0=λ2kx0 \mathbf{A}^k\mathbf{x}_0 = \lambda_2^k \mathbf{x}_0~。此时，算法将无法收敛到主特征向量 v1 ~\mathbf{v}_1~。相反，如果初始向量在主特征向量方向上的分量 c1 ~c_1~越大，幂法的收敛速度会越快。

  2. 特征值的分布：主特征值 λ1 ~\lambda_1~与次特征值 λ2 ~\lambda_2~的比值∣λ1λ2∣\left| \frac { \lambda _ { 1 } } { \lambda _ { 2 } } \right|越大，算法收敛越快。反之，如果该比值越接近 1 ~1~（也即 λ1≈λ2 ~\lambda_1 \approx \lambda_2~），收敛速度越慢。

  3. 矩阵的性质：对于正定矩阵（对称矩阵A=AT\mathbf{A} = \mathbf{A}^T，且特征值为正 ），幂法可以有效地分离主特征向量。而对于具有复特征值的矩阵，由于复特征值可能导致迭代过程中的振荡行为，使得算法无法收敛到实数特征向量，因此幂法可能会失效。

6. 反幂法：任意特征值估算

反幂法 (Inverse Power Method) ~(\textbf{Inverse Power Method})~最典型的应用场景是当我们需要计算矩阵的某个特定特征值，并且已经有一个接近该特征值的初始估计 α ~\alpha~。反幂法是基于幂法的思想，通过构造特定的矩阵 B ~\mathbf{B}~来实现对任意特征值的逼近。下面做一个简单推导，假设矩阵 A ~\mathbf{A}~的特征值为 λ1,λ2,⋯ ,λn ~\lambda_1,\lambda_2,\cdots,\lambda_n~，那么就有如下结论：
A∼[λ10⋯00λ2⋯0⋮⋮⋱⋮00⋯λn],A−αI∼[λ1−α0⋯00λ2−α⋯0⋮⋮⋱⋮00⋯λn−α]\mathbf{A} \sim \begin{bmatrix} \lambda_1 & 0 & \cdots & 0 \\ 0 & \lambda_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda_n \end{bmatrix},\quad \mathbf{A} - \alpha \mathbf{I} \sim \begin{bmatrix} \lambda_1 - \alpha & 0 & \cdots & 0 \\ 0 & \lambda_2 - \alpha & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda_n - \alpha \end{bmatrix}
我们构造矩阵 B=A−αI ~\mathbf{B}=\mathbf{A} - \alpha\mathbf{I}~，那么 B ~\mathbf{B}~的特征值就是 λ1−α, λ2−α, ⋯ , λn−α ~\lambda_1-\alpha,~\lambda_2-\alpha,~\cdots,~\lambda_n-\alpha~。因此， B−1=(A−λI)−1 ~\mathbf{B}^{-1}=(\mathbf{A} - \lambda\mathbf{I})^{-1}~的特征值如下：
λ1′=1λ1−α,λ2′=1λ2−α,⋯ ,λn′=1λn−α\lambda_1' = \frac{1}{\lambda_1-\alpha},\quad \lambda_2' = \frac{1}{\lambda_2-\alpha},\quad \cdots, \quad \lambda_n' = \frac{1}{\lambda_n-\alpha}
如果 α ~\alpha~接近某个特征值 λi ~\lambda_i~，则 (λi−α)−1 ~(\lambda_i-\alpha)^{-1}~会成为矩阵 B−1 ~\mathbf{B}^{-1}~的主特征值。这样，问题就被转化为了用幂法来近似计算 B−1 ~\mathbf{B}^{-1}~的主特征值。接下来就是选择一个初始向量 x0 ~\mathbf{x}_0~，然后开始迭代计算过程：
yk=(A−αI)−1xk ,(k=0,1,2,⋯ )\mathbf{y}_k = (\mathbf{A} - \alpha\mathbf{I})^{-1}\mathbf{x}_k\,,\quad (k = 0, 1, 2,\cdots)
要注意的是，在实际计算中并不会直接计算出逆矩阵 (A−αI)−1 ~(\mathbf{A} - \alpha \mathbf{I})^{-1}~后再和 xk ~\mathbf{x}_k~相乘计算 yk ~\mathbf{y}_k~，而是通过解线性方程组的方式来计算 yk ~\mathbf{y}_k~：
(A−αI) yk=xk(\mathbf{A} - \alpha \mathbf{I})\,\mathbf{y}_k = \mathbf{x}_k
后面迭代的步骤和幂法步骤完全一样，不过每次得到估计 λk′ ~\lambda_k'~后，需要转化为原特征值的估计：
λk=α+1λk′\lambda_k = \alpha + \frac{1}{\lambda_k'}
下面我们通过一个具体示例来演示反幂法的具体步骤。

7. 反幂法实例分析

给定矩阵 A ~\mathbf{A}~，初始估计 α ~\alpha~，初始向量 x0 ~\mathbf{x}_0~如下：
A=[10−8−4−8134−454],α=1.9 ,x0=[111]\mathbf{A}=\begin{bmatrix} 10 & -8 & -4 \\ -8 & 13 & 4 \\ -4 & 5 & 4 \end{bmatrix},\quad \alpha = 1.9\, ,\quad \mathbf{x}_0 = \begin{bmatrix}1 \\ 1 \\ 1\end{bmatrix}
迭代结束条件设置为相邻两个特征估计值的相对变化量 ϵ<10−6 ~\epsilon < 10^{-6}~。下面是 1 ~1~次迭代过程：

通过解线性方程组来计算 y0 ~\mathbf{y}_0~：
(A−1.9I)y0=x0(\mathbf{A} - 1.9\mathbf{I})\mathbf{y}_0 = \mathbf{x}_0
将初始条件代入，可得：
A−1.9I=[10−1.9−8−4−813−1.94−454−1.9]=[8.1−8−4−811.14−452.1]\mathbf{A} - 1.9\mathbf{I} = \begin{bmatrix} 10 - 1.9 & -8 & -4 \\ -8 & 13 - 1.9 & 4 \\ -4 & 5 & 4 - 1.9 \end{bmatrix} = \begin{bmatrix} 8.1 & -8 & -4 \\ -8 & 11.1 & 4 \\ -4 & 5 & 2.1 \end{bmatrix}
解方程组：
[8.1−8−4−811.14−452.1]y0=[111]\begin{bmatrix} 8.1 & -8 & -4 \\ -8 & 11.1 & 4 \\ -4 & 5 & 2.1 \end{bmatrix} y_0 = \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix}
解得：
y0≈[4.4500.5027.759]\mathbf{y}_0 \approx \begin{bmatrix} 4.450 \\ 0.502 \\ 7.759 \end{bmatrix}

用 y0 ~\mathbf{y}_0~的最大分量作为缩放因子进行归一化：
μ1=max⁡(∣4.450∣,∣0.502∣,∣7.759∣)=7.759λ1′=7.759\mu_1 = \max(|4.450|, |0.502|, |7.759|) = 7.759 \quad \lambda_1' = 7.759
更新向量：
x1=17.759y0=[0.5740.0651]\mathbf{x}_1 = \frac{1}{7.759}\mathbf{y}_0 = \begin{bmatrix}0.574 \\0.065 \\1\end{bmatrix}

λ1=α+1λ1′=1.9+17.759=2.029\lambda_1 = \alpha + \frac{1}{\lambda_1'} = 1.9 + \frac{1}{7.759} = 2.029

重复上面的过程，直到满足结束条件，由于初始预估值 α ~\alpha~和实际特征值很接近，所以算法收敛速度很快。下面是迭代的结果：

kk xk~\mathbf{x}_k yk~\mathbf{y}_k λk′ ~\lambda_k'~ λk ~\lambda_k~
0(0.57359 ,0.064649,1)(4.450374,0.501601,7.758805)7.7588052.028886
1(0.505364,0.004453,1)(5.013058,0.044172,9.9197)9.9197002.000809
2(0.500378,0.000313,1)(5.001247,0.003126,9.994931)9.9949312.000051
3(0.500027,0.000022,1)(5.000089,0.00022 ,9.999646)9.9996462.000004
4(0.500002,0.000002,1)(5.000006,0.000015,9.999975)9.9999752.000000
5(0.5, 1.09e-7, 1)(5, 1.09e-6, 9.999998)9.9999982.000000

5.7 在微分方程中的应用

5.9 在马尔科夫链中的应用

bei an浙ICP备2024102481号