逆幂法

逆幂法（inverse iteration）反复求解同一个线性系统并归一化，以提取接近指定移位的特征值方向。它利用逆算子将靠近移位的特征值放大为占优模式。实际每轮解方程，不应显式形成矩阵的逆；移位的选择既影响收敛速度，也影响线性系统的条件性。

靠近一的特征方向被放大

取 $A=\left(\begin{array}{cc}2& 1\\ 1& 2\end{array}\right)$，其特征值为一和三。选择移位 $\sigma =0.8$，从 $x_0=(1,0{)}^T$ 开始，先解 $$(A-0.8I)w=x_0.$$ 这里 $A-0.8I=\left(\begin{array}{cc}1.2& 1\\ 1& 1.2\end{array}\right)$，行列式为 $0.44$，解为 $w=(30/11,-25/11{)}^T$。按最大分量归一化后 $x_1=(1,-5/6{)}^T$，下一轮为 $x_2=(1,-60/61{)}^T$，迅速接近特征值一的方向 $(1,-1{)}^T$。

第一坐标归一化为一，第二坐标依次为0、−5/6、−60/61，目标方向为(1,−1)。

逆算子下的谱变换

若 $A{v}_j={\lambda }_j{v}_j$ 且 $\sigma$ 不等于任何特征值，便有 $$(A-\sigma I{)}^{-1}{v}_j=\frac{1}{{\lambda }_j-\sigma }{v}_j.$$ 特征向量不变，特征值距离的倒数成为新的放大倍数。假设 $A$ 可对角化，有唯一特征值 ${\lambda }_{\ast }$ 距 $\sigma$ 最近，并且初始向量在其方向上的展开系数非零。连续求解得到 $$(A-\sigma I{)}^{-k}{x}_0=({\lambda }_{\ast }-\sigma {)}^{-k}{c}_{\ast }[v_{\ast }+\sum _{j\ne \ast }\frac{c_j}{c_{\ast }}{\left(\frac{{\lambda }_{\ast }-\sigma }{{\lambda }_j-\sigma }\right)}^k{v}_j].$$ 括号中的其他系数几何衰减，归一化后方向趋向 $v_{\ast }$。这正是幂法作用在移位逆算子上的收敛证明。主例两个放大倍数为五和 $5/11$，方向误差比值为 $1/11$，比直接幂法的三分之一更小。

特征值估计与主例的精确公式

获得方向后，仍可用原矩阵的 Rayleigh 商 $\mu =x^{T}Ax/(x^{T}x)$ 估计原特征值，而非把逆算子的放大倍数当成答案。主例按上述归一化有 $$x_k=\left(\begin{array}{c}1\\ -(11^k-1)/(11^k+1)\end{array}\right),\mu (x_k)=1+\frac{2}{121^k+1}.$$ 第一轮商为 $62/61\approx 1.01639344$，第二轮为 $1+2/14642\approx 1.00013659$。实际验收还应计算原残差 $Ax-\mu x$。近似解线性系统时的残差与特征对残差是两个不同量，均需要与容限匹配。

移位过近、等距与改变移位

移位越接近简单目标特征值，逆迭代的方向分离通常越快，但 $A-\sigma I$ 越接近奇异。若移位正好等于特征值，普通求解式无法定义；有些专门特征向量算法会利用近奇异系统，但不能把它们当成本篇的普通逆迭代。

有解变式取 $\sigma =2$。它到一和三的距离都为一，逆矩阵为 $\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)$；从 $(1,0)$ 出发在两个坐标轴之间交替，不能选择唯一方向。若初始向量本来就是目标特征向量，任意合法移位的一轮求解都只改变长度与符号，归一化后方向立即正确。

每轮用新 Rayleigh 商作为移位得到 Rayleigh 商迭代。它对实对称矩阵在简单特征向量附近可有很快的局部收敛，但移位改变后通常需要重新分解矩阵；不能同时宣称始终复用同一分解。对大稀疏系统，移位线性求解成本常是决定性因素。

可运行的固定移位逆迭代

代码输入实矩阵、非零初向量和固定移位，返回原特征值估计、单位方向及轮数，不修改输入。为简洁每轮调用 NumPy 通用求解，所以稠密 $K$ 轮成本约 $O(K{n}^3)$、存储 $O(n^2)$。生产实现应预先分解 $A-\sigma I$，将后续三角求解降为每轮 $O(n^2)$，稀疏填充则另计。

Python 3

import numpy as np


def inverse_iteration(A, start, shift, tol=1e-12, maxiter=100):
    A, x = np.array(A, dtype=float, copy=True), np.array(start, dtype=float, copy=True)
    shift = float(shift)
    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 not np.isfinite(shift) 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")
    B = A-shift*np.eye(x.size)
    if not np.isfinite(B).all():
        raise ArithmeticError("shifted matrix outside floating-point range")
    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):
        mu_scaled = float(x@(C@x)); mu = mu_scaled*matrix_scale
        residual = C@x-mu_scaled*x
        if not np.isfinite(mu) or not np.isfinite(residual).all():
            raise ArithmeticError("nonfinite eigenvalue or residual")
        if np.linalg.norm(residual) <= threshold:
            return mu, x, k
        if k == maxiter:
            break
        w = np.linalg.solve(B, x)
        magnitude = np.max(np.abs(w))
        if magnitude == 0 or not np.isfinite(magnitude) or not np.isfinite(w).all():
            raise ArithmeticError("nonfinite or zero inverse iterate")
        w = w/magnitude
        x = w/np.linalg.norm(w)
    raise RuntimeError("original eigenpair residual did not reach tolerance")


print(inverse_iteration([[2, 1], [1, 2]], [1, 0], 0.8)[0])  # near 1
print(inverse_iteration([[2, 1], [1, 2]], [1, 0], 3.2)[0])  # near 3

对非正规矩阵或近重特征值，方向误差与残差之间的关系可能较差；严格的误差解释需考虑条件数和谱间隔。求得某一个特征对也不等于找全靠近移位的全部特征值，多维子空间方法更适合成簇目标。