幂法与逆幂法

在上一节中，我们定义了矩阵的特征值和特征向量，并介绍了一种通用的求特征值的方法——QR 算法。要理解该算法的收敛原理，多少有点麻烦。本文中我们先介绍两个迭代求特征值的方法，做些铺垫，再来研究QR算法的收敛原理。

假设 $\lambda, \boldsymbol{\alpha}$ 是矩阵 $\boldsymbol{A}$ 的特征值和特征向量，满足关系 $\boldsymbol{A}\boldsymbol{\alpha} = \lambda \boldsymbol{\alpha}$。在等式两侧分别左乘行向量 $\boldsymbol{\alpha}^T$，有 $\boldsymbol{\alpha}^T\boldsymbol{A}\boldsymbol{\alpha} = \lambda\boldsymbol{\alpha}^T\boldsymbol{\alpha}$。因此，根据特征向量 $\boldsymbol{\alpha}$ 我们可以直接写出对应的特征值 $\lambda = \frac{\boldsymbol{\alpha}^T\boldsymbol{A}\boldsymbol{\alpha}}{\boldsymbol{\alpha}^T\boldsymbol{\alpha}}$。对于任意的向量 $\boldsymbol{x}$，都可以像这样计算出一个值 $r(\boldsymbol{x})$，称之为瑞利商Rayleigh Quotient: $$ r(\boldsymbol{x}) = \frac{\boldsymbol{x}^T\boldsymbol{A}\boldsymbol{x}}{\boldsymbol{x}^T\boldsymbol{x}} $$ $|r(\boldsymbol{x}) - \lambda|$常被拿来衡量向量 $\boldsymbol{x}$ 与 $\lambda$ 对应的特征向量有多接近。幂法与逆幂法的迭代过程，实际上就是瑞利商不断收敛到特征值的过程。幂法会收敛到绝对值最大的那个特征值上，逆幂法则收敛到绝对值最小的那个特征值上。

1. 幂法(Power Method)

设一个 $n \times n$ 的矩阵 $\boldsymbol{A}$ 有 $n$ 个线性无关的特征向量 $\{ \boldsymbol{\alpha}_1, \cdots, \boldsymbol{\alpha}_n \}$，对应的特征值为 $\{\lambda_1, \cdots, \lambda_n\}$，并且满足 $|\lambda_1| > |\lambda_2| ≥ \cdots ≥ |\lambda_n|$。其中 $\lambda_1$ 是绝对值最大的那个特征值，有些资料里面称之为主特征值。这 $n$ 个特征向量构成了一组基，向量 $\boldsymbol{x}$ 在这组基下的坐标为 $\{ x_1, \cdots, x_n\}$，我们任选一个 $x_1 \neq 0$ 的向量 $\boldsymbol{\nu}$ 作为迭代的初始向量，对其展开如下: $$ \boldsymbol{\nu} = x_1 \boldsymbol{\alpha}_1 + x_2 \boldsymbol{\alpha}_2 + \cdots + x_n \boldsymbol{\alpha}_n $$ 我们将矩阵 $\boldsymbol{A}$ 分别左乘到上面等式的左右两侧，根据特征值与特征向量的定义有 $$ \boldsymbol{A}\boldsymbol{\nu} = x_1 \lambda_1\boldsymbol{\alpha}_1 + x_2 \lambda_2\boldsymbol{\alpha}_2 + \cdots + x_n \lambda_n \boldsymbol{\alpha}_n $$ 左乘 $k$ 次之后，有 $$ \boldsymbol{A}^k\boldsymbol{\nu} = x_1 \lambda_1^k\boldsymbol{\alpha}_1 + x_2 \lambda_2^k\boldsymbol{\alpha}_2 + \cdots + x_n \lambda_n^k \boldsymbol{\alpha}_n \\ $$ $$ \boldsymbol{A}^k\boldsymbol{\nu} = x_1 \lambda_1^k\left( \boldsymbol{\alpha}_1 + \frac{x_2 \lambda_2^k}{x_1 \lambda_1^k}\boldsymbol{\alpha}_2 + \cdots + \frac{x_n \lambda_n^k}{x_1 \lambda_1^k}\boldsymbol{\alpha}_n \right) $$

        double PowerIterate(MatrixA const & A, VectorV & v) {
            double re = 0; v = v / v.Norm();
            for (int i = 0; i < max_iter; i++) {
                auto w = A * v;
                re = v.Dot(w);
                v = w / w.Norm();
            }
            return re;
        }

令 $\boldsymbol{\delta}_k = \frac{x_2 \lambda_2^k}{x_1 \lambda_1^k}\boldsymbol{\alpha}_2 + \cdots + \frac{x_n \lambda_n^k}{x_1 \lambda_1^k}\boldsymbol{\alpha}_n$，有 $\boldsymbol{A}^k\boldsymbol{\nu} = x_1 \lambda_1^k\left(\boldsymbol{\alpha}_1 + \boldsymbol{\delta}_k\right)$，对其归一化，有

$$ \boldsymbol{\kappa} = \frac{\boldsymbol{A}^k\boldsymbol{\nu}}{\|\boldsymbol{A}^k\boldsymbol{\nu}\|} = \frac{x_1 \lambda_1^k}{|x_1 \lambda_1^k|}\frac{\boldsymbol{\alpha}_1 + \boldsymbol{\delta}_k}{\|\boldsymbol{\alpha}_1 + \boldsymbol{\delta}_k\|} $$

因为 $|\lambda_1| > |\lambda_2| ≥ \cdots ≥ |\lambda_n|$ 所以当 $k \to \infty$ 时，有 $\boldsymbol{\delta}_k \to 0, \boldsymbol{\kappa} \to \boldsymbol{\alpha}_1$，所以有瑞利商 $r(\boldsymbol{\kappa}) \to \lambda_1$: $$ r(\boldsymbol{\kappa}) = \frac{\boldsymbol{\kappa}^T\boldsymbol{A}\boldsymbol{\kappa}}{\boldsymbol{\kappa}^T\boldsymbol{\kappa}} \approx \lambda_1 $$

所以幂法的迭代过程实际就只有三步: ① 对迭代向量归一化 ② 左乘矩阵 $\boldsymbol{A}$ ③ 计算瑞利商。如左侧的伪代码所示，一直重复这个过程直到瑞利商收敛，就可以得到绝对值最大的特征值，及其对应的特征向量。具体的实现参见PowerIterate。幂法收敛有两个条件，其一需要主特征值唯一，即 $|\lambda_1| > |\lambda_2|$。如果$\lambda_1, \lambda_2$的绝对值比较接近，会降低收敛速度。其二初始向量需要在主特征方向上有分量，一般我们会选择一个全为 1 的向量作为初始向量。

2. 逆幂法(Inverse Power Method)

逆幂法与幂法的区别就在于，它每次迭代的时候都是左乘的逆矩阵 $\boldsymbol{A}^{-1}$，一般用来求解矩阵的绝对值最小的那个特征值。如果矩阵 $\boldsymbol{A}$ 可逆，并且其特征值为 $\lambda_1, \cdots, \lambda_n$，并且满足 $|\lambda_1| ≥ |\lambda_2| ≥ \cdots ≥ |\lambda_n|$，那么可以证明它的逆矩阵 $\boldsymbol{A}^{-1}$ 的特征值为 $1/\lambda_1, \cdots, 1/\lambda_n$，并且保持对应的特征向量不变。

        double InversePowerIterate(MatrixA const & A, VectorV & v) {
            double re = 0; LU lu(A);
            for (int i = 0; i < max_iter; i++) {
                lu.Solve(v, v);
                v = v / v.Norm();
                re = v.Dot(A * v);
            }
            return re;
        }

因为 $\boldsymbol{A}$ 可逆，所以不存在某个特征值为 0。根据特征值的定义，有 $\boldsymbol{A}\boldsymbol{x} = \lambda \boldsymbol{x}$。在等式的左右两侧分别左乘 $\boldsymbol{A}$ 的逆，有 $\boldsymbol{A}^{-1}\boldsymbol{A}\boldsymbol{x} = \lambda \boldsymbol{A}^{-1}\boldsymbol{x}$，整理一下并在等式的两侧分别除以 $\lambda$，有 $\boldsymbol{A}^{-1}\boldsymbol{x} = \frac{1}{\lambda}\boldsymbol{x}$。所以$\boldsymbol{A}$和$\boldsymbol{A}^{-1}$具有相同的特征向量，并且特征值成倒数关系。如果我们对矩阵 $\boldsymbol{A}^{-1}$ 进行幂法迭代，那么将收敛到$\boldsymbol{A}^{-1}$ 的绝对值最大的特征值 $1 / \lambda_n$ 上，对应就是 $\boldsymbol{A}$ 的最小特征值。

算法的实现也很简单，只需要对幂法做些简单改造。如右侧的伪代码所示，我们先对目标矩阵 $\boldsymbol{A}$ 进行 LU 分解，这样可以在迭代的过程中通过求解线性方程组的方式替代左乘逆矩阵 $\boldsymbol{A}^{-1}$。然后就是对迭代之后的向量 v 归一化，并计算其关于原始矩阵 $\boldsymbol{A}$ 的瑞利商，得到特征值的估计。详细的实现参见 InversePowerIterate。

3. 带偏移量的逆幂法

幂法和逆幂法一个只能求最大特征值，另一个只能求最小特征值，应用场景有限。如果我们已知特征值接近某个特定值，我们通常会选择带偏移量的逆幂法来求出精确的结果。

对于矩阵 $\boldsymbol{A}$ 假设我们希望在 $u$ 附近找到一个特征值。这个 $u$ 就是所谓的偏移量，我们可以构建偏移矩阵 $\boldsymbol{B} = \boldsymbol{A} - u\boldsymbol{I}$。如果 $\lambda, \boldsymbol{x}$ 是矩阵$\boldsymbol{A}$ 的特征值和特征向量，$\boldsymbol{A}\boldsymbol{x} = \lambda\boldsymbol{x}$，那么很容易证明，矩阵 $\boldsymbol{B}$ 的特征值为 $\lambda - u$，并且特征向量与 $\boldsymbol{A}$ 相同，即$\boldsymbol{B}\boldsymbol{x} = (\lambda - u)\boldsymbol{x}$。此时我们对矩阵 $\boldsymbol{B}$ 实施逆幂法，就可以直接求出最接近偏移量 $u$ 的那个特征值。

我们只需要将逆幂法中对矩阵 $\boldsymbol{A}$ 进行 LU 分解，修改为对 $\boldsymbol{A} - u \boldsymbol{I}$ 进行 LU 分解就可以了。详细的实现，参见OffInvPowerIterate。

幂法与逆幂法

1. 幂法(Power Method)

2. 逆幂法(Inverse Power Method)

3. 带偏移量的逆幂法

4. 完