幂法

幂法（power method）反复用矩阵乘一个向量并归一化，以提取模最大的特征值方向。它不需要完整矩阵分解，适合能够快速计算矩阵向量积的大型问题。成功的关键是主特征值的模严格占优，以及初始向量确实含有该方向的分量。

同一矩阵连续放大一个方向

取 $A=\left(\begin{array}{cc}2& 1\\ 1& 2\end{array}\right)$、$x_0=(1,0{)}^T$。乘一次得到 $(2,1{)}^T$，用最大分量归一化为 $x_1=(1,1/2{)}^T$；下一次得到 $(5/2,2{)}^T$，归一化为 $x_2=(1,4/5{)}^T$；再一次为 $x_3=(1,13/14{)}^T$。方向逐渐接近 $(1,1{)}^T$。

这个矩阵沿 $(1,1{)}^T$ 放大三倍，沿 $(1,-1{)}^T$ 只放大一倍。初始向量恰好是两个方向各一半，因此 $$A^k{x}_0=\frac{1}{2}\left(\begin{array}{c}{3}^k+1\\ 3^k-1\end{array}\right),x_k=\left(\begin{array}{c}1\\ (3^k-1)/(3^k+1)\end{array}\right).$$ 归一化避免数值随 $3^k$ 增长，却保留方向的相对比例。

按最大分量归一化后第一坐标恒为一，第二坐标逐步趋于一。图展示方向变化，而不是向量长度的指数增长。

模占优条件的收敛证明

设矩阵可对角化，特征向量为 $v_1,\dots ,v_n$，特征值按模排序满足 $|{\lambda }_1|>|{\lambda }_2|\ge \cdots$。将初始向量写成 $x_0=\sum _i{c}_i{v}_i$，并要求 $c_1\ne 0$。则 $$A^k{x}_0={\lambda }_1^k{c}_1[v_1+\sum _{i=2}^n\frac{c_i}{c_1}{\left(\frac{{\lambda }_i}{{\lambda }_1}\right)}^k{v}_i].$$ 所有括号内的次要分量都以模小于一的比值衰减，因此归一化方向收敛到 $v_1$，渐近方向误差主要由 $|{\lambda }_2/{\lambda }_1{|}^k$ 控制。主特征值为负时向量可能每轮变号，所以应比较方向或残差，不能要求带符号向量本身收敛。

可对角化假设使上式直接成立。含 Jordan 块的矩阵可能多出多项式因子，非正规矩阵也可能出现很长的暂态；不能只按两个特征值的比值预报从第一轮起的误差。幂法、逆迭代与正交迭代的关系可参见 FNC 的特征值迭代章节。

特征值估计与残差验收

非零向量的 Rayleigh 商为 $\mu (x)=x^{T}Ax/(x^{T}x)$。在实对称情形，将单位向量写成 $x=\cos \theta v_1+\sin \theta w$，其中 $w\perp v_1$、$\Vert w\Vert =1$，交叉项消失，得到 $$\mu (x)-{\lambda }_1={\sin }^2\theta (w^{T}Aw-{\lambda }_1).$$ 因此特征值误差对小方向夹角为二阶；这个结论依赖对称性。主例有精确公式 $$\mu (x_k)=3-\frac{2}{9^k+1},$$ 所以首两轮估计分别为 $14/5=2.8$、$122/41\approx 2.97560976$。

验收还应计算 $r=Ax-\mu x$。对单位向量和实对称矩阵，若正交特征展开为 $x=\sum _j{c}_j{v}_j$，则 $\Vert r{\Vert }^2=\sum _j{c}_j^2({\lambda }_j-\mu {)}^2\ge \underset{j}{min}|{\lambda }_j-\mu {|}^2$。所以至少存在一个特征值距 $\mu$ 不超过残差范数；这不自动证明该值就是模最大的一个，也不在缺少谱间隔时单独保证特征向量准确。

有解变式把初始向量改成 $(1,-1{)}^T$。它与主方向正交，所有迭代都停留在特征值一的方向，得到的是正确的次特征对。另一个反例是平面九十度旋转矩阵，两个复特征值模相同，实向量不断旋转，不会选择一个实主方向。二者分别说明初始分量和严格模间隔的重要性。

可运行的归一化迭代

代码接受实矩阵和非零初始向量，返回 Rayleigh 商、单位向量和迭代数，不修改输入。它以实际特征对残差停止，并不验证模占优条件；调用者应根据上述条件解释输出。稠密每轮 $O(n^2)$，稀疏算子每轮可为 $O(\mathrm{nnz}(A))$，幂迭代的向量工作区为 $O(n)$；下面的教学实现另行复制并缩放稠密矩阵，总额外存储为 $O(n^2)$。

Python 3

import numpy as np


def power_iteration(A, start, tol=1e-12, maxiter=1000):
    A, x = np.array(A, dtype=float, copy=True), np.array(start, dtype=float, copy=True)
    if x.ndim != 1 or not x.size or A.shape != (x.size, x.size):
        raise ValueError("square matrix matching vector required")
    if not np.isfinite(A).all() or not np.isfinite(x).all() or np.max(np.abs(x)) == 0:
        raise ValueError("finite data and nonzero start required")
    if not np.isfinite(tol) or tol <= 0 or not isinstance(maxiter, int) or maxiter < 0:
        raise ValueError("positive tolerance and nonnegative integer limit required")
    x /= np.max(np.abs(x)); x /= np.linalg.norm(x)
    matrix_scale = float(np.max(np.abs(A)))
    C = A/matrix_scale if matrix_scale > 0 else A.copy()
    scale = max(1.0, np.linalg.norm(C))
    threshold = tol*scale
    if not np.isfinite(threshold):
        raise ArithmeticError("nonfinite stopping threshold")
    for k in range(maxiter+1):
        w = C@x; mu_scaled = float(x@w)
        mu = mu_scaled*matrix_scale
        if not np.isfinite(w).all() or not np.isfinite(mu):
            raise ArithmeticError("nonfinite matrix-vector product")
        if np.linalg.norm(w-mu_scaled*x) <= threshold:
            return mu, x, k
        if k == maxiter:
            break
        length = np.linalg.norm(w)
        if length == 0 or not np.isfinite(length):
            raise ArithmeticError("cannot normalize next iterate")
        x = w/length
    raise RuntimeError("eigenpair residual did not reach tolerance")


print(power_iteration([[2, 1], [1, 2]], [1, 0])[:1])  # near (3.0,)
print(power_iteration([[2, 1], [1, 2]], [1, -1])[:1]) # (1.0,), a different invariant direction

若要查找靠近指定数值的特征值，可以使用逆幂法；同时追踪多个方向则需要子空间迭代或 Arnoldi、Lanczos 方法。它们改变搜索空间或算子，不能仅靠多运行几次同一随机初值幂法替代完整的特征值算法。