引言
特征值问题 $Ax = \lambda x$ 是科学与工程计算中最基本的代数问题之一：结构动力学中的固有频率与振型、量子力学中的哈密顿量能级、控制系统的稳定性判据、搜索引擎的 PageRank 排序，本质上都归结为求解某类矩阵的特征值与特征向量。由于 Abel–Ruffini 定理指出五次以上多项式方程不存在一般根式解，当矩阵阶数 $n \ge 5$ 时，通过特征多项式 $\det(A - \lambda I) = 0$ 求根的思路在数学上就走不通，一切实用算法都必须是迭代式的。本文讨论两类最具代表性的算法：面向中小规模稠密矩阵的 QR 算法，以及面向大规模稀疏对称矩阵的 Lanczos 算法，并分析其收敛机理与工程选型依据。
理论基础
相似变换与 Schur 分解
特征值算法的基本策略，是通过一系列相似变换 $A \leftarrow Q^{-1}AQ$（特征值在此变换下不变）把矩阵化简为容易读出特征值的形式。Schur 定理保证：任意复方阵 $A$ 都存在酉矩阵 $Q$，使得 $Q^{*}AQ = T$ 为上三角矩阵，$T$ 的对角元即为全部特征值。实矩阵对应实 Schur 形——拟上三角矩阵，对角线上允许出现 $2\times 2$ 小块以承载共轭复特征值。QR 算法的目标正是数值地逼近这一 Schur 形。
从幂法到 QR 迭代
最简单的迭代是幂法：任取 $v_0$，反复执行 $v_{k+1} = Av_k/|Av_k|$。当 $|\lambda_1| > |\lambda_2|$ 时，$v_k$ 以收敛因子 $|\lambda_2/\lambda_1|$ 线性收敛到主特征向量。幂法一次只能给出一个特征对，而 QR 迭代可以理解为”对整个子空间同时做幂法”：
$$A_k = Q_k R_k, \qquad A_{k+1} = R_k Q_k = Q_k^{T} A_k Q_k$$
即对 $A_k$ 做 QR 分解后交换次序相乘。可以证明 $A_k$ 收敛到（拟）上三角矩阵，且对角元的收敛速度由相邻特征值模的比值决定。直接实现每步需 $O(n^3)$ 的 QR 分解，工程上不可接受，因此实用 QR 算法分两步走：先用 Householder 变换把 $A$ 一次性化为上 Hessenberg 形 $H$（代价 $O(n^3)$，且 Hessenberg 形在 QR 迭代下保持不变），此后每步迭代只需 $O(n^2)$；再引入位移策略加速：
$$H_k - \mu_k I = Q_k R_k, \qquad H_{k+1} = R_k Q_k + \mu_k I$$
取 Wilkinson 位移（右下角 $2\times 2$ 子矩阵中靠近 $h_{nn}$ 的那个特征值），对称矩阵可实现接近三次的局部收敛。配合收缩（deflation）技术——当某次对角元低于阈值时把问题分裂为更小的独立子问题——整体复杂度约 $O(n^3)$ 且数值稳定。QR 算法因此被公认为 20 世纪最重要的十大算法之一。
Lanczos 算法与 Krylov 子空间
当阶数达到 $n \sim 10^5$ 甚至更高（例如有限元离散得到的刚度矩阵），$A$ 通常高度稀疏，但 QR 迭代中的 $Q$ 因子是稠密的，仅存储就已失败。此时真正廉价的信息只有矩阵-向量乘积 $v \mapsto Av$。Lanczos 算法在 Krylov 子空间
$$\mathcal{K}_m(A, q_1) = \mathrm{span}{q_1, Aq_1, \dots, A^{m-1}q_1}$$
内寻找特征对的近似。对对称矩阵 $A$，逐次构造标准正交基 ${q_j}$，由对称性，Arnoldi 过程中的长递推退化为三项递推：
$$\beta_{j+1} q_{j+1} = A q_j - \alpha_j q_j - \beta_j q_{j-1}, \qquad \alpha_j = q_j^{T} A q_j,\quad \beta_j = |r_{j-1}|_2$$
于是 $Q_m^{T} A Q_m = T_m$ 为 $m \times m$ 对称三对角矩阵。计算 $T_m$ 的特征值（Ritz 值）作为 $A$ 的近似特征值，相应的 $Q_m y$ 为 Ritz 向量。Kaniel–Paige 收敛理论表明：谱两端（最大、最小）的特征值率先快速收敛，内部特征值收敛缓慢。有限精度下 Lanczos 向量会逐渐丢失正交性，产生”幽灵”重特征值，工程实现需配合完全重正交化或选择性重正交化。每步仅需一次稀疏矩阵乘法与 $O(n)$ 的向量运算，存储 $O(mn)$，这正是它能处理百万阶问题的根本原因。
数值方法与实现
import numpy as np
def shifted_qr_eig(A, max_iter=200, tol=1e-12):
"""带 Wilkinson 位移的 QR 算法（对称矩阵版）"""
H = np.array(A, dtype=float)
n = H.shape[0]
it = 0
while n > 1 and it < max_iter * n:
it += 1
b = H[n-2, n-1]
d = (H[n-2, n-2] - H[n-1, n-1]) / 2.0
s = np.sign(d) if d != 0 else 1.0
mu = H[n-1, n-1] - s * b**2 / (abs(d) + np.hypot(d, b))
Q, R = np.linalg.qr(H[:n, :n] - mu * np.eye(n))
H[:n, :n] = R @ Q + mu * np.eye(n)
if abs(H[n-1, n-2]) < tol * (abs(H[n-1, n-1]) + abs(H[n-2, n-2])):
n -= 1 # 收缩：一个特征值已收敛
return np.sort(np.diag(H))
def lanczos(matvec, n, q1, m, reorth=True):
"""Lanczos 三对角化：matvec(v) 提供矩阵-向量乘积，返回 Ritz 值"""
Q = np.zeros((n, m)); alpha = np.zeros(m); beta = np.zeros(m - 1)
q = q1 / np.linalg.norm(q1)
r = matvec(q)
alpha[0] = q @ r
for j in range(m):
Q[:, j] = q
r = r - alpha[j] * q
if j > 0:
r = r - beta[j-1] * Q[:, j-1]
if reorth: # 完全重正交化
r = r - Q[:, :j+1] @ (Q[:, :j+1].T @ r)
if j == m - 1:
break
beta[j] = np.linalg.norm(r)
q = r / beta[j]
r = matvec(q)
alpha[j+1] = q @ r
T = np.diag(alpha) + np.diag(beta, 1) + np.diag(beta, -1)
return np.linalg.eigvalsh(T)
rng = np.random.default_rng(0)
B = rng.standard_normal((40, 40))
A = B + B.T # 对称稠密矩阵
print(np.allclose(shifted_qr_eig(A), np.linalg.eigvalsh(A), atol=1e-8))
# 稀疏情形：一维离散 Laplacian，只用矩阵-向量乘积
n = 400
diag = 2.0 * np.ones(n); off = -1.0 * np.ones(n - 1)
def matvec(v):
return diag * v + np.concatenate(([0.0], off * v[:-1]))
+ np.concatenate((off * v[1:], [0.0]))
ritz = lanczos(matvec, n, rng.standard_normal(n), m=30)
exact = 2 - 2*np.cos(np.arange(1, 6)*np.pi/(n + 1)) # 最小的 5 个特征值
print(ritz[:5]); print(exact) # 端部特征值率先收敛工程案例分析
桥梁、高层建筑与旋转叶轮的模态分析需求解广义特征值问题 $K\phi = \omega^2 M\phi$，其中 $K$、$M$ 分别为刚度矩阵与质量矩阵。通过对 $M$ 做 Cholesky 分解 $M = LL^{T}$，可将其化为标准对称问题 $\tilde{K}\psi = \lambda\psi$。一座中等规模桥梁的有限元模型自由度数常在 $10^5$ 到 $10^6$ 量级，但工程设计通常只需要最低数十阶模态——这正是 Lanczos 及其块版本、带位移-反演变体的主场，ANSYS、Abaqus 等商业软件的模态求解器内核即构建于此。反之，控制系统极点配置、量子化学中 $n \le 10^4$ 的稠密问题，直接调用基于 QR/分治法的 LAPACK 例程（如 syevd）更为划算。
| 算法 | 适用矩阵 | 产出 | 单步代价 | 典型规模 |
|---|---|---|---|---|
| 幂法/反幂法 | 任意（稀疏友好） | 单个极端特征对 | $O(\mathrm{nnz})$ | 任意大 |
| QR（Hessenberg+位移） | 中小稠密 | 全部特征值 | $O(n^2)$（化简后） | $n less\sim 10^4$ |
| Jacobi 扫描 | 小稠密对称 | 全部特征对 | $O(n^3)$ | $n less\sim 10^3$ |
| Lanczos | 大稀疏对称 | 端部若干特征对 | 一次 SpMV | $n gtr\sim 10^6$ |
| 分治法（syevd） | 稠密对称 | 全部特征对 | $O(n^3)$，常数小 | $n less\sim 10^4$ |
选型时还需关注谱分布：当目标特征值聚集在一起时，Lanczos 收敛显著恶化，宜改用块 Lanczos 或谱变换（shift-invert，通过求 $A - \sigma I$ 的逆把目标区间翻转到谱的端部）。内存受限时，隐式重启动技术（IRLM，即 ARPACK 的核心）用隐式位移 QR 步压缩 Krylov 子空间维数，把存储控制在可接受范围内。
结论与展望
QR 与 Lanczos 代表了特征值计算的两条基本路线：前者以正交相似变换逐步压缩矩阵结构、一次求得全谱；后者以 Krylov 投影充分利用稀疏性、按需提取端部谱。理解二者的共同语言——正交变换、Rayleigh–Ritz 投影与收缩机制——是阅读和选用现代求解器（FEAST、SLEPc、LOBPCG 等）的基础。近年来，随机化算法（randomized SVD/eig）以几乎与精度要求无关的低廉复杂度给出次优但极快的近似，正在数据科学场景中与传统确定性算法形成有力互补。
发表回复