幂法计算矩阵的特征值与特征向量

本文介绍一种简单有效的计算矩阵特征值与特征向量的近似值的方法——幂法。
主要参考《数值计算方法与算法》，《数值方法：设计、分析和算法实现 》

一、幂法基本方法

在实际问题中, 矩阵的按模最大特征值往往起更重要的作用. 例如: 矩阵的谱半径即矩阵的按模最大特征值的值, 决定了迭代矩阵是否收敛。因此, 矩阵的按模最大的特征值比其余的特征值的地位更加重要. 幂法是计算矩阵按模最大特征值及相应的特征向量的数值方法.

简单地说，任取初始向量 $X^{(0)}$ ，进行迭代计算

$$X^{(k+1)}=A X^{(k)}$$

得到迭代序列 $\left\{X^{(k)}\right\}$, 分析 $X^{(k+1)}$ 与 $X^{(k)}$ 之间的关系, 得到 $A$ 的按模最大特征值及特征向量的近似解.

在幂法中, 假设矩阵 $A$ 有特征值 $\lambda_i, i=1,2, \cdots, n$, 其中

$$\left|\lambda_1\right| \geqslant\left|\lambda_2\right| \geqslant\left|\lambda_3\right| \geqslant \cdots \geqslant\left|\lambda_n\right|$$

并有 $n$ 个线性无关的特征向量 $v_i, A v_i=\lambda_i v_i, i=1,2, \cdots, n$.
任取初始向量 $X^{(0)}, X^{(0)}$ 可由 $A$ 的 $n$ 个线性无关的特征向量 $v_i$ 的线性表示. 设

$$X^{(0)}=\alpha_1 v_1+\alpha_2 v_2+\cdots+\alpha_n v_n$$

那么,

$$X^{(1)}=A X^{(0)}=\alpha_1 \lambda_1 v_1+\alpha_2 \lambda_2 v_2+\cdots+\alpha_n \lambda_n v_n$$

一般地, 有

$$X^{(k)}=A X^{(k-1)}=\alpha_1 \lambda_1^k v_1+\alpha_2 \lambda_2^k v_2+\cdots+\alpha_n \lambda_n^k v_n$$

(1) 按模最大的特征值只有一个, 且是单实根.

设 $\left|\lambda_1\right|>\left|\lambda_2\right| \geqslant\left|\lambda_3\right| \geqslant \cdots\left|\lambda_n\right|$,

由上面的讨论我们知道

$$\begin{aligned} X^{(k)} & =\alpha_1 \lambda_1^k v_1+\alpha_2 \lambda_2^k v_2+\cdots+\alpha_n \lambda_n^k v_n \\ & =\lambda_1^k\left[\alpha_1 v_1+\alpha_2\left(\frac{\lambda_2}{\lambda_1}\right)^k v_2+\cdots+\alpha_n\left(\frac{\lambda_n}{\lambda_1}\right)^k v_n\right] \end{aligned}$$

若 $\alpha_1 \neq 0$, 由于 $\left|\frac{\lambda_i}{\lambda_1}\right|<1, i=2,3, \cdots, n$, 对充分大的 $k$ :

$$\left|\frac{\lambda_i}{\lambda_1}\right|^k \rightarrow 0, \quad i=2,3, \cdots, n$$

故

$$X^{(k)} \approx \lambda_1^k \alpha_1 v_1, \quad X^{(k+1)} \approx \lambda_1^{k+1} \alpha_1 v_1=\lambda_1 \lambda_1^k \alpha_1 v_1=\lambda_1 X^{(k)}$$

得到按模最大的特征值

$$\lambda_1 \approx x_i^{(k+1)} / x_i^{(k)}, \quad i=1,2, \cdots, n$$

相应的特征向量近似地为 $X^{(k)}$.

显然$\left\{x_i^{(k+1)} / x_i^{(k)}\right\}$ 收敛于 $\lambda_1$ 的速度取决于比值 $\left|\frac{\lambda_2}{\lambda_1}\right|$ 的大小

(2) 按模最大的特征值是互为反号的实根.

设 $\lambda_1>0$, 且 $\lambda_1=-\lambda_2$, 即 $\left|\lambda_1\right|=\left|\lambda_2\right|>\left|\lambda_3\right| \geqslant \cdots \geqslant\left|\lambda_n\right|$. 这时有

$$X^{(k)}=\lambda_1^k\left(\alpha_1 v_1+(-1)^k \alpha_2 v_2+\alpha_3\left(\frac{\lambda_3}{\lambda_1}\right)^k v_3+\cdots+\alpha_n\left(\frac{\lambda_n}{\lambda_1}\right)^k v_n\right)$$

当 $k$ 充分大时，有

$$\begin{aligned} & X^{(k)} \approx \lambda_1^k\left(\alpha_1 v_1+(-1)^k \alpha_2 v_2\right) \\ & X^{(k+1)}=\lambda_1^{k+1}\left(\alpha_1 v_1+(-1)^{k+1} \alpha_2 v_2\right) \quad\left(X^{(k)} \text { 与 } X^{(k+1)} \text { 没有比例关系 }\right) \\ & X^{(k+2)} \approx \lambda_1^{k+2}\left(\alpha_1 v_1+(-1)^{k+2} \alpha_2 v_2\right) \approx \lambda_1^2 X^{(k)} \end{aligned}$$

对充分大的 $k, X^{(k)}$ 与 $X^{(k+2)}$ 近似一个常数因子 $\lambda_1^2$. 所以

$$\lambda_1=\sqrt{x_i^{(k+2)} / x_i^{(k)}}$$

再计算相应的特征值

得

$$\begin{aligned} & \left\{\begin{array}{l} X^{(k+1)}=\lambda_1^{k+1}\left(\alpha_1 v_1+(-1)^{k+1} \alpha_2 v_2\right) \\ X^{(k)}=\lambda_1^k\left(\alpha_1 v_1+(-1)^k \alpha_2 v_2\right) \end{array}\right. \\ & \left\{\begin{array}{l} X^{(k+1)}+\lambda_1 X^{(k)} \approx 2 \lambda_1^{k+1} \alpha_1 v_1 \\ X^{(k+1)}-\lambda_1 X^{(k)} \approx(-1)^{(k+1)} 2 \lambda_1 \alpha_2 v_2 \end{array}\right. \end{aligned}$$

按特征向量的性质, 相应的特征向量可以取为

$$\left\{\begin{array}{l} v_1=X^{(k+1)}+\lambda_1 X^{(k)} \\ v_2=X^{(k+1)}-\lambda_1 X^{(k)} \end{array}\right.$$

还有很多更复杂的情况, 可参考有关书籍或选用其他方法. 这里只讨论两种情况. 在计算中若 $x_i^{(k+1)} / x_i^{(k)}$ 的比值趋于一个稳定的值, 则属于第一种情况, 如果 $x_i^{(k+1)} / x_i^{(k)}$ 不能趋于一个稳定的值, 但是 $x_i^{(2 k+2)} / x_i^{(2 k)}$ 和 $x_i^{(2 k+1)} / x_i^{(2 k-1)}$ 的比值分别趋于一个稳定的值, 则属于第二种情况。

二、幂法的规范计算

在幂法迭代计算 $X^{(k+1)}=A X^{(k)}$ 中, 当 $k$ 充分大时, 若 $A$ 的按模最大特征值其绝对值较大时， $X^{(k)}$ 的某些分量迅速增大, 或许会超过计算机实数的值域 (上溢); 而 $A$ 的按模最大特征值其绝对值较小时, $X^{(k)}$ 的分量迅速缩小,当 $k$ 充分大时, 或许会被计算机当零处理 (下溢). 因此, 分量过大或过小都会导致计算失败。在实际计算中, 通常采用规范运算, 对 $X^{(k)}$ 的每个元素除以 $X^{(k)}$ 按模最大的分量 $\left\|X^{(k)}\right\|_{\infty}=\max _{1 \leqslant i \leqslant n}\left|x_i^{(k)}\right|$

规范运算可按下面公式进行

$$\left\{\begin{array}{l} Y^{(k)}=X^{(k)} /\left\|X^{(k)}\right\|_{\infty}, \\ X^{(k+1)}=A Y^{(k)}, \end{array} \quad k=0,1, \cdots\right.$$

规范化运算保证了 $\left\|Y^{(k)}\right\|=1$, 即 $Y^{(k)}$ 按摸最大分量的值保持为 1 或 -1 .下面给出在规范运算中迭代序列的几种情况:

(1) 如果 $\left\{X^{(k)}\right\}$ 收敛, 则 $A$ 的特征值按模最大分量的值仅有一个, 且 $\lambda_1>0$,对充分大的 $k$, 按模最大分量 $x_j^{(k)}$ 不变号, 对应的 $\left|y_j^{(k)}\right|=1, \lambda_1 \approx\left|x_j^{(k+1)}\right|$, 即

$$\lambda_1 \approx \max _{1 \leqslant i \leqslant n}\left|x_i^{(k+1)}\right|=\left|x_j^{(k+1)}\right|$$

相应的特征向量是 $v_1 \approx Y^{(k)}$.

(2) 如果 $\left\{X^{(2 k)}\right\},\left\{X^{(2 k+1)}\right\}$ 分别收敛于互为反号的向量, 则按模最大的特征值也仅有一个单实根, 且 $\lambda_1<0$, 即对充分大的 $k$, 若 $x_j^{(k)}$ 的符号交替变号, 则 $\lambda_1$为负值.

$$\lambda_1 \approx-\max _{1 \leqslant i \leqslant n}\left|x_i^{(k)}\right|=-\left|x_j^{(k)}\right|$$

相应的特征向量是 $v_1 \approx Y^{(k)}$.

(3) 如果 $\left\{X^{(2 k)}\right\},\left\{X^{(2 k+1)}\right\}$ 分别收敛于两个不同的向量（与（2）不同），则按模最大的特征值有两个, 是互为反号的一对实根. 这时, 对充分大的 $k$, 再作一次非规范运算

$$X^{(k+1)}=A X^{(k)}$$

则

$$\lambda_1 \approx \sqrt{x_i^{(k+1)} / y_i^{(k-1)}}, \quad \lambda_2=-\lambda_1$$

而仍有

$$\left\{\begin{array}{l} V_1=X^{(k+1)}+\lambda_1 X^{(k)} \\ V_2=X^{(k+1)}-\lambda_1 X^{(k)} \end{array}\right.$$

(4) 如果 $\left\{X^{(k)}\right\}$ 的趋势无一定的规律, 这时 $A$ 的按模最大的特征值的情况更为复杂，需要另行处理.

三、幂法加速：原点位移法

矩阵 $B=A-p I$ 特征值与矩阵 $A$ 的特征值分别为 $\lambda$ 与 $\lambda-p$. 通过计算矩阵 $A-p I$ 的特征值得到矩阵 $A$ 的特征值的方法称为原点位移法.

若取得的 $p$ 满足

$$\left|\frac{\lambda_2-p}{\lambda_1-p}\right|<\left|\frac{\lambda_2}{\lambda_1}\right|$$

则对矩阵 $B$ 的幂法, 收敛会比矩阵 $A$ 的幂法要快. 只有对矩阵特征值的分布有个大致了解, 才能有效选取 $p$

四、反幂法

反幂法是计算矩阵按模最小的特征值以及相应的特征向量的数值方法，实际上同样使用的是幂法方法

设矩阵 $A$ 可逆， $\lambda$ 和 $v$ 分别为 $A$ 的特征值以及相应的特征向量，对 $A v=\lambda v$两边同乘 $A^{-1}$, 得 $A^{-1} v=\frac{1}{\lambda} v$. 可见 $A$ 和 $A^{-1}$ 的特征值互为倒数, 而且 $v$ 也是 $A^{-1}$ 的特征值 $\frac{1}{\lambda}$ 的特征向量. $A^{-1}$ 的按模最大的特征值正是 $A$ 的按模最小的特征值的倒数. 用幂法计算 $A^{-1}$ 按模最大的特征值从而得到 $A$ 的按模最小的特征值的方法, 称为反幂法. 显然, $A$ 的按模最小特征值越接近于 0 , 收敛越快.

用幂法计算 $A^{-1}$ 按模最大的特征值仍可用规范方法：

任取 $X^{(0)}$ ，规范迭代

$$\left\{\begin{array}{l} Y^{(k)}=X^{(k)} /\left\|X^{(k)}\right\|_{\infty}, \quad k=0,1, \cdots . . . . \\ X^{(k+1)}=A^{-1} Y^{(k)}, \end{array}\right.$$

在实际计算中不是先算出 $A^{-1}$, 再作乘积 $A^{-1} Y^{(k)}$ ，而是解方程 $A X^{(k+1)}=$ $Y^{(k)}$, 求得 $X^{(k+1)}$. 由于需要反复求解, 一般不用 Gauss 消元法, 而用直接分解法解方程. 将 $A$ 分解为 $A=\mathrm{LU}$, 记 $U X^{(k+1)}=Z^{(k+1)}$, 由 $L Z^{(k+1)}=Y^{(k)}$ 解出 $Z^{(k+1)}$,再由 $U X^{(k+1)}=Z^{(k+1)}$ 解出 $X^{(k+1)}$ 。

规范迭代计算公式：

$$\left\{\begin{array}{l} Y^{(k)}=X^{(k)} /\left\|X^{(k)}\right\|_{\infty}, \\ A X^{(k+1)}=Y^{(k)}, \end{array} \quad k=0,1, \cdots\right.$$

若我们想知道最接近 $p$ 的特值值 $\lambda_i$, 可以使用带原点位移的反幂法计算公式

$$\left\{\begin{array}{l} Y^{(k)}=X^{(k)} /\left\|X^{(k)}\right\|_{\infty}, \\ (A-p I) X^{(k+1)}=Y^{(k)}, \end{array} \quad k=0,1, \cdots\right.$$

计算得到特征值 $\mu$,

$$\lambda_i=p+\frac{1}{\mu}$$

五、幂法的程序实现

下面给出一个简单的幂法计算示例，用Python实现

import numpy as np
def power_method(A, num_iterations=1000, tolerance=1e-9):

    # 随机初始化特征向量
    n = A.shape[0]
    b_k = np.random.rand(n)

    # 归一化初始向量
    b_k = b_k / np.linalg.norm(b_k)

    for _ in range(num_iterations):
        b_k1 = np.dot(A, b_k)
        eig_val = np.linalg.norm(b_k1)
        # 归一化特征向量
        b_k1 = b_k1 / eig_val
        # 判断是否收敛
        if np.linalg.norm(b_k1 - b_k) < tolerance:
            return eig_val, b_k1

        b_k = b_k1

    # 如果未收敛，返回最后的结果
    return eig_val, b_k


if __name__ == "__main__":
    A = np.array([
        [1, 1 ,0.5],
        [1, 1 ,0.25],
        [0.5, 0.25, 2]
    ])

    eigenvalue, eigenvector = power_method(A)
    print("最大特征值:", eigenvalue)
    print("对应的特征向量:", eigenvector)

六、小结

总结一下，幂法和反幂法可以帮助我们求解矩阵的主特征值与最小模值特征值，通过原点平移以及归一化方法，可以加速幂法的计算。但这一方法也有很大的局限性，例如并不是所有的矩阵都可以在幂法下收敛，以及很多时候，我们同样需要关心除了最大与最小模特征值以外的特征值，就需要选用其他方法了。