引言
科学与工程计算的核心问题之一是求解大型稀疏线性方程组 $Ax = b$。有限元、有限差分离散 PDE 产生的矩阵规模常达 $10^6 \sim 10^9$ 阶，直接法（LU 分解）的 $\mathcal{O}(n^3)$ 运算量与 $\mathcal{O}(n^2)$ 存储完全不可行，而矩阵-向量乘 $Av$ 却能利用稀疏性以 $\mathcal{O}(n)$ 代价完成。Krylov 子空间方法正是只依赖矩阵-向量乘、把 $n$ 维问题投影到低维子空间求解的一类迭代法，配合预处理技术已成为大规模计算的标准引擎。本文介绍 Krylov 子空间方法的投影原理、共轭梯度法（CG）与 GMRES 的构造和收敛理论，并讨论预处理子的设计原则。
理论基础
Krylov 子空间与投影原理
给定初始残差 $r_0 = b – Ax_0$，$m$ 阶 Krylov 子空间定义为
$$\mathcal{K}_m(A, r_0) = \mathrm{span}{r_0, Ar_0, A^2 r_0, \dots, A^{m-1}r_0}$$
Krylov 方法在第 $m$ 步从仿射空间 $x_0 + \mathcal{K}_m$ 中选取近似解 $x_m$，并施加 Petrov-Galerkin 正交条件 $r_m = b - Ax_m \perp \mathcal{L}_m$。不同的约束空间 $\mathcal{L}_m$ 导出不同方法：CG 取 $\mathcal{L}_m = \mathcal{K}_m$（要求 $A$ 对称正定），GMRES 取 $\mathcal{L}_m = A\mathcal{K}_m$（适用于一般非奇异矩阵）。
Arnoldi 与 Lanczos 过程
直接以 ${A^j r_0}$ 为基在数值上是病态的（向量迅速趋向主特征方向），需要正交化。Arnoldi 过程用 Gram-Schmidt 逐步构造 $\mathcal{K}_m$ 的标准正交基 $V_m$，得到关系 $AV_m = V_{m+1}\bar{H}_m$，其中 $\bar{H}_m$ 为上 Hessenberg 矩阵。GMRES 极小化残差范数等价于求解最小二乘问题 $\min_y |\beta e_1 - \bar{H}_m y|_2$。当 $A$ 对称时，Hessenberg 矩阵退化为三对角矩阵，Arnoldi 退化为短递推的 Lanczos 过程——这是 CG 每步只需常数工作量的根本原因。
共轭梯度法及其收敛理论
对对称正定矩阵，CG 等价于在能量范数 $|e|_A = \sqrt{e^T A e}$ 下求最优近似，其误差满足经典界
$$\frac{|e_m|_A}{|e_0|_A} \le 2\left(\frac{\sqrt{\kappa} - 1}{\sqrt{\kappa} + 1}\right)^m, \qquad \kappa = \kappa_2(A) = \frac{\lambda_{\max}}{\lambda_{\min}}$$
更精细的分析表明收敛速度由特征值分布决定：特征值聚集在少数几个区间时收敛远快于上述均匀界，且 $A$ 仅有 $k$ 个互异特征值时 CG 至多 $k$ 步精确收敛。GMRES 的收敛界同样依赖特征值分布（正规矩阵情形），对高度非正规矩阵还涉及伪谱与极小多项式。
预处理技术
预处理的目的是用一个「近似 $A$ 且易于求解」的矩阵 $M$ 改造谱分布。左预处理求解 $M^{-1}Ax = M^{-1}b$，右预处理求解 $AM^{-1}u = b$；对 SPD 情形常用分裂预处理 $M = LL^T$ 保持对称性。常用预处理子包括：Jacobi（$M = \mathrm{diag}(A)$）、Gauss-Seidel/SSOR、不完全 LU/Cholesky 分解（ILU/IC，通过控制填充水平平衡质量与成本）、稀疏近似逆（SPAI），以及针对 PDE 问题的多重网格与区域分解预处理——后者可使迭代步数不随网格加密而增长，达到最优复杂度。
数值方法与实现
下面用纯 NumPy 实现 CG 与 Jacobi 预处理 CG（PCG），作用于变系数扩散方程 $-(a(x)u’)’ = 1$（$a(x) = 1 + 4x$，Dirichlet 边界）的有限差分离散系统。变系数使主对角不再恒定，Jacobi 预处理才有实际意义：
import numpy as np
def var_coeff_1d(n):
"""-(a u')' = 1 的中心差分离散, a(x) = 1 + 4x"""
h = 1.0 / (n + 1)
ah = 1.0 + 4.0 * np.linspace(0.5 * h, 1.0 - 0.5 * h, n + 1)
d = (ah[:-1] + ah[1:]) / h**2
e = -ah[1:-1] / h**2
A = np.diag(d) + np.diag(e, 1) + np.diag(e, -1)
return A, np.ones(n)
def cg(A, b, M=None, tol=1e-10, maxit=5000):
"""(预条件)共轭梯度法, M 为预处理子(此处取对角向量)"""
x = np.zeros_like(b)
r = b - A @ x
z = r if M is None else r / M
p = z.copy()
rz = r @ z
for k in range(1, maxit + 1):
Ap = A @ p
alpha = rz / (p @ Ap)
x = x + alpha * p
r = r - alpha * Ap
if np.linalg.norm(r) < tol:
return x, k
z = r if M is None else r / M
rz_new = r @ z
p = z + (rz_new / rz) * p
rz = rz_new
return x, maxit
n = 2000
A, b = var_coeff_1d(n)
_, it_cg = cg(A, b)
_, it_pcg = cg(A, b, M=np.diag(A))
kappa = np.linalg.cond(A)
print(f"条件数 kappa(A) = {kappa:.3e}")
print(f"CG 迭代步数: {it_cg}, Jacobi-PCG 迭代步数: {it_pcg}")运行可见 PCG 步数明显少于 CG，量级差异符合 $\sqrt{\kappa}$ 规律：迭代步数大致正比于 $\sqrt{\kappa(M^{-1}A)}$。实际代码中 $A$ 绝不应显式存储稠密矩阵，而应以稀疏格式或「无矩阵」（matrix-free）方式只提供 $v \mapsto Av$ 接口。
工程案例分析
二维 Poisson 方程（静电场、稳态温度场、不可压流的压力投影）经五点差分离散后，$n = N^2$ 阶矩阵的条件数按 $\kappa \approx 4N^2/\pi^2$ 增长：$N = 1000$ 时 $\kappa \sim 4\times10^5$，裸 CG 需要上万步。工程中通常的选型是：Poisson 类问题用多重网格预处理 CG（步数与 $N$ 基本无关）；对流扩散等非对称问题用 GMRES 或 BiCGStab 配合 ILU；结构力学中的块状低秩或刚体模态问题则需更专门的 deflation 技术。下表对比主流 Krylov 方法：
| 方法 | 适用矩阵 | 递推长度 | 每步代价 | 极小化性质 | ||
|---|---|---|---|---|---|---|
| CG | 对称正定 | 短（3 项） | 1 次 $Av$ + 数次内积 | $ | e | _A$ 最小 |
| MINRES | 对称不定 | 短 | 1 次 $Av$ + 内积 | $ | r | _2$ 最小 |
| GMRES | 一般矩阵 | 全递推 | 1 次 $Av$ + $\mathcal{O}(km)$ 正交化 | $ | r | _2$ 最小 |
| BiCGStab | 一般矩阵 | 短 | 2 次 $Av$ | 无严格最优性 |
结论与展望
Krylov 子空间方法把 $n$ 维代数问题转化为 $m \ll n$ 维投影问题，其效率由特征值分布而非矩阵阶数决定；预处理则是改造谱分布、把理论优势兑现为实际加速的关键。当前研究热点包括：面向 GPU 的通信避免与流水线化 Krylov 变体（pipelined CG）、随机化预处理（randomized sketching）以及用神经网络学习问题相关预处理子的混合方法。理解「谱分布决定收敛」这一核心思想，是在工程中为具体问题正确选型迭代法与预处理子的前提。
发表回复