从大二数值分析课那个学期开始，我就对数线性方程组的解法很着迷。当时课本里最核心的就是高斯消元、LU分解，解法完完整整，矩阵一扔进去，一顿操作出一个“精确解”。后来真正用有限差分、有限元去解偏微分方程，才意识到现实世界里的矩阵有多大：三维网格一加密，未知量轻松到百万级，直接法那套“保存在内存里慢慢消元”的思路直接崩。于是Krylov子空间方法成了绕不开的主角。这篇文章我尽量把高等数值分析里关于Krylov子空间方法的核心逻辑讲清楚，不只是摆一遍CG和GMRES的公式，而是把“为什么需要子空间”“每个算法在优化什么”“收敛性为什么波动”“预处理为什么不可或缺”这些底层问题一次说明白。适合正在学数值分析、要写有限元求解器，或者被大型稀疏矩阵折磨的工程师参考。

1. 直接法在大规模问题上的天花板在哪里

1.1 高斯消元的存储噩梦与填充现象

先从一个具体例子看起。求解二维泊松方程 (-u_{xx}-u_{yy}=f) 时，用五点差分格式离散，(N\times N) 的内部网格点对应 (n=N^2) 个未知量，离散矩阵是三对角块加边结构，稀疏性非常好，每行只有5个非零元。按n=10^6估算，存储稀疏矩阵本体只要约4000万个浮点数，也就是几百MB级别，这还在现代机器承受范围内。

但直接法一做LU分解，情况就完全变了。消元过程中原本是零的位置会被填充（fill-in），二维问题填充后矩阵带宽大概是O(N)，算下来需要存储的非零元数量级从O(n)一下跳到O(n^{3/2})；三维问题更夸张，会接近O(n^2)。用1亿个未知量的三维网格做一次完整LU分解，内存需求会到几十TB乃至更高，这显然不是普通工作站能扛住的。

所以在“稀疏大规模”这个场景下，直接法的主要瓶颈从来不是浮点运算次数，而是对内存的贪婪消耗。这也是我后来在课堂上反复强调的一点：评估一个线性求解器是否适用，第一个问题应该是“矩阵结构和规模如何”，而不是“精度有多高”。

1.2 经典迭代法为何慢到无法接受

既然直接法在大规模下吃亏，那就想到迭代法。最基础的Jacobi、Gauss-Seidel迭代，每次迭代只需做几次矩阵向量乘，单步成本低，内存消耗也小。但问题在于收敛速度。

以模型泊松问题为例，Jacobi迭代的收敛因子约为 (\cos(\pi h))，网格步长 (h) 越小，收敛因子越接近1，需要的迭代次数按 (O(N^2)) 增长。换句话说，你为了精确解加密网格，没想到迭代次数反而暴涨，最终整体耗时根本不比直接法省心。这类方法适合作为光滑器或预处理子，却扛不起求解器的大梁。

这个瓶颈背后的原因很有意思：经典迭代法每一轮只是在“局部修正误差”，它通过矩阵分裂 (A=D-L-U) 将问题转化成不动点迭代，步与步之间没有任何全局信息。我们要的下一步迭代修正量，是整个空间里让残差下降最多的方向，这需要从矩阵全局提取信息。

1.3 Krylov线索的出现：幂向量里藏着解的信息

现在我们换个视角。给定一个初始近似解 (x_0)，可以得到初始残差 (r_0=b-Ax_0)。如果我们不断用矩阵去作用这个残差，就生成一串向量：

[ r_0,\quad Ar_0,\quad A^2r_0,\quad A^3r_0,\ldots ]

这些向量构成所谓Krylov子空间的自然生成元：

[ \mathrm{span}{r_0, Ar_0, A^2r_0,\ldots,A^{m-1}r_0} ]

为什么这个奇怪的向量序列里藏着解的线索？因为如果解 (x) 可以写成误差形式 (x=x_0+z)，则残差满足 (r_0=Az)。当矩阵可逆时 (z=A^{-1}r_0)，而根据Hamilton-Cayley定理，(A^{-1}) 可以表示为 (A) 的次数不超过 (n-1) 的多项式。也就是说，理论上存在一组系数，使得精确解被 (A) 作用在 (r_0) 上的多项式组合精确表示。

既然如此，我们完全可以放弃在高维空间里硬碰硬，转而到Krylov子空间里找近似解。矩阵越大，这种“降维打击”的价值越明显：你不需要知道全空间的几何，只需要顺着矩阵自己的“流向”一步步走。这就是Krylov子空间方法的出发点。

2. 搭好子空间的骨架：Arnoldi过程与投影原理

2.1 如何生成子空间里的正交基底

Krylov子空间的自然生成元 ({r_0, Ar_0, A^2r_0,\ldots}) 在数值上并不适合直接用：随着幂次增加，这些向量会快速趋向矩阵最大特征值对应的特征向量方向，导致所有向量挤成一个方向，失去了表征子空间的能力。

解决办法是逐次正交化。Arnoldi过程就是标准Gram-Schmidt思想在Krylov序列上的应用：先归一化 (v_1=r_0/|r_0|)，然后对每一个 (k=1,2,\ldots,m)，计算 (w=A v_k)，把 (w) 在之前所有基底上的投影减掉，再归一化得到 (v_{k+1})。整个过程写成矩阵关系非常漂亮：

[ A V_m = V_m H_m + h_{m+1,m}v_{m+1}e_m^T ]

其中 (V_m=[v_1,\ldots,v_m]) 是正交基底，(H_m) 是 (m\times m) 的上Hessenberg矩阵，也就是下三角部分除次对角线外全为零的那种结构。这个公式的意义在于：高维矩阵A作用到子空间上的效果，被完全压缩到这个小很多的H_m上，我们后面所有的优化都是在H_m上做的。

2.2 Arnoldi递推的完整实现

Arnoldi过程的算法骨架并不复杂，核心循环可以写成这样：

import numpy as np def arnoldi(A, v0, m): n = A.shape[0] V = np.zeros((n, m+1)) H = np.zeros((m+1, m)) v0_norm = np.linalg.norm(v0) if v0_norm == 0: raise ValueError("初始向量为零") V[:, 0] = v0 / v0_norm for j in range(m): w = A @ V[:, j] for i in range(j+1): H[i, j] = np.dot(V[:, i], w) w -= H[i, j] * V[:, i] H[j+1, j] = np.linalg.norm(w) if H[j+1, j] < 1e-15: # 这是happy breakdown，子空间已经不变，可以提前停机 break V[:, j+1] = w / H[j+1, j] return V, H

这个代码里最容易被忽略的一个数值细节是：教材里写经典Gram-Schmidt是三行伪代码，但直接这样写，在大规模计算中会因为舍入误差丢失正交性，导致整个Arnoldi过程失真。经验做法是用修正Gram-Schmidt（MGS），也就是每减掉一个投影分量，立即对剩余向量做归一化前的修正。

2.3 理解Hessenberg矩阵的“扁平投影”

很多初学者困惑：Arnoldi过程跑了一堆正交化，最后得到的H_m为什么是Hessenberg结构？道理其实在他的正交化策略里：每步只拿新向量 (w=A v_k) 与已经生成的前 (k) 个基底做正交化，减去这些投影之后，剩下分量与 (v_1,\ldots,v_k) 都正交，它唯一可能与下一个新基底 (v_{k+1}) 对齐。因此在H_m的列上，第 (k) 列只有前 (k+1) 行非零，下标恰好呈阶梯型。

这个结构价值很大。它意味着我们可以用递推方式一点一点扩展子空间，而不必每次重新生成整组基底。后面用GMRES求解残差最小化问题时，H_m的Hessenberg结构也能让最小二乘子问题保持在很低的计算复杂度。

3. GMRES与CG：两者真正差异在于优化目标不同

3.1 GMRES：在Krylov子空间里极小化残差范数

有了Arnoldi给出的正交基底 (V_m)，我们自然想问：在这个m维子空间里，哪个向量让当前残差范数最小？既然解限定为 (x=x_0+V_m y)，那么新残差为：

[ r_m = r_0 - A V_m y = \beta v_1 - V_{m+1}\tilde{H}_m y ]

其中 (\beta=|r_0|)，(\tilde{H}m) 是包含最后一行的 ((m+1)\times m) 上Hessenberg矩阵。由于 (V{m+1}) 列正交，极小化 (|r_m|) 等价于极小化下面的小规模最小二乘问题：

[ \min_y |\beta e_1 - \tilde{H}_m y| ]

这个问题用QR分解或Givens旋转就能稳定求解，代价远小于原问题的规模。GMRES名字里的“Generalized Minimal Residual”说的就是这个过程——在每一步都取当前子空间里的最优残差。

GMRES最直观的优点是对任意非奇异矩阵都能走。但它有一个工程痛点：随着m增大，需要保存全部m个正交基底向量，每次Arnoldi还要和所有已存向量做正交化，存储和计算成本按 (O(m^2)) 增长速度上升。实际工程里几乎不会跑完整GMRES，而是用重启动版本GMRES(m)，每迭代m步，把得到的近似解作为新初值重新开始一组Krylov子空间。但重启动会让GMRES失去全局最优性——你只优化了当前窗口，无法记住更早的子空间信息。这个特性在病态矩阵上体现得很明显，后面会专门讲。

3.2 对称正定矩阵：Lanczos退化和CG的出现

如果A是对称矩阵，Arnoldi过程中的H_m会变成对称的三对角矩阵，正交基的递推关系也从需要跟所有历史向量正交，退化成只跟前两个基底正交的Lanczos三递推。这也是对称问题里往往用Lanczos而不是Arnoldi的原因，存储成本从O(mn)降到O(n)，好用太多。

在此基础上进一步假设A对称正定（SPD），就得到经典共轭梯度法（CG）。很多教材直接扔出CG递推公式，学生背会了却不知道它在优化什么。CG实际上是在Krylov子空间里极小化如下能量范数误差：

[ |x-x_k|_A^2 = (x-x_k)^T A (x-x_k) ]

而不是残差二范数。这意味着CG每次迭代得到的是在当前子空间里“按A范数理解”的最优解。这个区别直接导致CG与GMRES在同样条件下的收敛曲线形态很不一样。CG的实现非常简单：

def cg(A, b, x0=None, tol=1e-10, max_iter=1000): x = np.zeros_like(b) if x0 is None else x0.copy() r = b - A @ x p = r.copy() rs_old = np.dot(r, r) for k in range(max_iter): Ap = A @ p alpha = rs_old / np.dot(p, Ap) x += alpha * p r -= alpha * Ap rs_new = np.dot(r, r) if np.sqrt(rs_new) < tol: return x, k + 1 p = r + (rs_new / rs_old) * p rs_old = rs_new return x, max_iter

这段代码里最反直觉的一行是方向更新 (p=r+\beta p)。为什么每次新的搜索方向要加上旧方向的修正？因为单纯沿负梯度方向走会重复做无用功，加上旧方向修正后，能保证每个新搜索方向是关于A共轭的，也就是满足 (p_i^T A p_j=0)，这样理论上有限步就能在精确算术下收敛到n步内。

3.3 其他Krylov家族成员与适用边界

一旦理解GMRES是“无脑对残差求最小”，CG是“对称正定下的能量最优”，就能理解整个Krylov方法家族为什么长成这副模样。

MINRES可以从“对称版本的GMRES”理解，它不做A共轭方向，而是用Lanczos基底直接极小化残差。BiCGSTAB设计思路更取巧，它不用保存全部基向量，但用双正交条件构造了两个方向的迭代，代价是放弃了每一步严格的单调残差下降。实际用起来BiCGSTAB偶发震荡，如果残差不降反升，就要警惕数值稳定性问题。

4. 收敛性分析：Krylov方法为什么时快时慢

4.1 收敛不单纯由条件数决定

很多资料会给出CG收敛上界：

[ |x-x_k|_A \le 2\left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^k |x-x_0|_A ]

其中 (\kappa=\kappa(A)) 是谱条件数。这个公式容易给人错觉：只要条件数大，迭代就一定慢。但实际上这只是“最坏情况”的估计，真实收敛速度通常远好于它，原因是决定Krylov方法收敛速度的核心是矩阵特征值的分布形状，而不只是最大与最小特征值的比值。

举个例子，下面两个矩阵条件数都是1000，但特征值分布截然不同：一个是特征值在区间[0.001, 1]内均匀散布；另一个是100个特征值聚集在0.001附近，剩下9900个全部挤在1周围。后者对于Krylov方法来说收敛会非常快，因为A的谱被“聚簇”了，Krylov多项式只需要把簇内特征值一起处理。这也是为什么预处理的目标总是“让特征值聚拢”而不是单纯压条件数。

4.2 非对称矩阵的困难与“最坏情况”

GMRES的收敛理论更复杂，对正规矩阵（满足 (AA^T=A^TA)）可以用谱分布给出类似上界。但对非正规矩阵，仅凭特征值分布无法解释实际收敛行为，伪谱（pseudospectrum）比谱更能说明问题。我见过工程里有些非对称矩阵，特征值全在右半平面，看起来人畜无害，实际GMRES迭代几十步残差纹丝不动，这正是伪谱效应在作怪。

所以对非对称问题，我通常建议先跑几十步GMRES看残差曲线，再判断是继续硬跑、换预处理，还是改用其他算法族。纸上算收敛阶，在非对称领域参考价值远低于对称正定领域。

4.3 超线性收敛和停滞期：残差曲线的“心跳”

实际观察CG或GMRES的残差曲线，经常发现它不是一条平滑下降的直线，而是先慢吞吞下降、甚至中途出现一个平台期，然后突然加速。这个现象被称作超线性收敛。原因是Krylov多项式在前若干步还没有“认全”特征值的信息，对某些特征值方向的衰减很弱；随着迭代继续，多项式在这些方向上的取值被逐渐压低，整体收敛速度就会变快。

这一点在日常调试中有直接指导意义：看到残差曲线进入平台，不一定要立刻判定方法失效，可以观察平台长度。如果平台出现在迭代早期，多迭代几步可能会突降；如果平台出现在重启动后的每一轮，说明重启动m太小，Krylov子空间根本来不及积累足够信息，这时候应当增大m。

4.4 理想中断与数值边界

还有个现象叫happy breakdown，在Arnoldi或Lanczos过程中表现为 (h_{m+1,m}) 恰好为零。这意味着当前Krylov子空间在A的作用下不再向外扩展，子空间已经是A的不变子空间，此时如果 (r_0\in K_m) 且矩阵在这个子空间上非奇异，那么GMRES得到的解将直接是精确解。这个“意外之喜”判断代码里很简单：一旦当前残差小于阈值，立即退出循环，避免无意义的后续迭代。

不过实际计算中更常遇到的是 (h_{m+1,m}) 非常小但不为零，比如降到1e-14量级。此时再除以它归一化v_{m+1}会让数值噪声被放大，后续基底质量急剧恶化。稳妥的做法是引入软阈值：如果 (h_{m+1,m}) 小于某个经验值（我常用的下限是1e-14乘以当前矩阵范数的量级），就把它视为零并终止Arnoldi递推，把当前解当作已收敛处理。

5. 预处理：同等代码量下收益最大的操作

5.1 为什么必须预处理

一个简单的椭圆型问题，五点差分离散后矩阵条件数约为 (O(h^{-2}))。网格从100×100细化到1000×1000，条件数增长一万倍，直接跑CG，按最坏上界看迭代次数会成比例上升。实际虽然没那么夸张，但网格加密几倍后迭代次数增加一个数量级是很正常的。

物理解释更直观：网格越密，相邻未知量之间的耦合越强，低频误差分量的衰减就越慢。这些低频分量数学上对应A的小特征值方向，正是它拉高了条件数。预处理的作用就是把A变成另一个矩阵，让新矩阵的特征值分布更集中，把低频分量也变成快速收敛的方向。

5.2 三个层次的预处理：对角、ILU、矩阵分裂

最简单的预处理是Jacobi预处理（对角缩放）：取 (M=\mathrm{diag}(A))，求解 (M^{-1}Ax=M^{-1}b)。对对角占优矩阵效果立竿见影，代码只需要一行，几乎零风险。但它的能力上限很低，对强耦合的PDE离散问题帮助有限。

再往上走是ILU（不完全LU分解）族。它模拟直接法的消元过程，但严格控制填充：ILU(0)只保留原矩阵非零结构上的因子；ILUT基于阈值丢弃小于阈值的填充元。这类预处理在实际工程里非常有价值，特别是非对称问题。一个关键经验是ILU预处理的有效性强烈依赖矩阵排序，用反向Cuthill-McKee（RCM）重排序后，ILU(0)的表现通常会有肉眼可见的提升。

还有一类矩阵分裂预处理可以看作Gauss-Seidel思想的近亲，比如SSOR预处理选取 (M=\frac{1}{2-\omega}(D+\omega L)D^{-1}(D+\omega U))。这类预处理胜在完全基于分裂，不需要额外分解，适合并行实现，但效果一般不如ILU锐利。

5.3 左预处理与右预处理的取舍

预处理不是随便把 (M^{-1}) 乘到方程左边就叫完事。左预处理：

[ M^{-1}Ax=M^{-1}b ]

虽然理论上和右预处理

[ AM^{-1}y=b,\qquad x=M^{-1}y ]

有相近的收敛性质，但注意：左预处理改变了残差定义。GMRES在最小化残差范数时，实际上最小化的是 (|M^{-1}(b-Ax)|)，这个“加权残差”不同于原始残差。好处是如果 (M^{-1}A) 谱分布越好，收敛越快；坏处是如果M构造得很差，加权后的收敛标准会和真实残差脱节。

右预处理没有这个问题，因为GMRES可以保持对原始残差 (|b-Ax|) 的监控，实现上更加透明。所以我的默认选择是右预处理，除非有特殊需要才会用左预处理版本。

5.4 一个直观的小实验

为了说明预处理的效果，我经常在课堂上演示这样一个二维Poisson例子：100×100网格，用五点差分生成SPD矩阵A，分别统计未预处理CG与Jacobi预处理CG（PCG）的迭代次数。

未预处理CG大约需要两百多次迭代收敛到 (10^{-10})，而Jacobi预处理后降到五十多次，ILU(0)预处理后一般二十次以内就能收敛。这组数字的价值不在于“Jacobi预处理多么神奇”，而在于它说明了一个机制：哪怕只是把对角线拉平，Krylov方法也能节省一大半工作量。现实中面对三维复杂几何模型，一个好的预处理器带来的加速可以达到几十倍甚至上百倍，远超优化CPU指令的收益。

6. 课堂之外的排坑经验与实用建议

6.1 停止准则不要只盯相对残差

工程代码里最常见的错误是只设置相对残差阈值 (|r_k|/|b| < \text{tol})。听起来合理，但如果初始残差本身就远大于解的范数，这个准则可能在迭代早期就误判收敛。我见过一个流体求解器因为这样的停止准则，把“看起来收敛”的近似解拿去做后处理，结果压力场全是震荡。

更稳健的做法是同时监控两个指标：绝对残差 (|r_k|) 和相对残差 (|r_k|/|r_0|)，以二者都满足为收敛条件。如果二者无法同时达标，优先相信绝对残差，结合解的物理量量纲判断可行性。

6.2 重启动m的选择策略

GMRES(m)里的m不是越大越好，也不是越小越省。太小，比如m=10，对很多非对称问题基本会陷进周期性的停滞，每轮重启后残差只下降一点点，甚至原地踏步。太大，每步正交化要跟所有历史向量算内积，成本呈二次方增长。

我的实验习惯是先用m=30跑一个算例，输出残差曲线观察。如果曲线在接近收敛时突然“断崖”说明m刚好够用；如果残差在某个高位反复震荡，就依次尝试50、80、100。大部分非对称工程问题，m=50到80这个区间是个甜点区。

6.3 检查对称性再选求解器

“对称”在工程矩阵里未必像教科书那样干净。网格生成、边界条件处理、装配顺序稍有差池，就可能出现类似 (10^{-14}) 量级的非对称扰动。这时如果直接用CG，结果可能完全跑飞，因为CG的整个推导建立在A对称正定之上。

所以在决定用CG之前，建议在代码里显式检查：

sym_error = np.linalg.norm(A - A.T) / np.linalg.norm(A)

这个值大于1e-12就要警惕，考虑换MINRES或GMRES。这个检查代码只要两行，却能避免一大类莫名其妙的调试问题，值得长期保留。

6.4 矩阵向量的性能陷阱

Krylov方法九十多步都在做矩阵向量乘。同一个矩阵，用CSR格式和坐标格式做乘法，性能可能差一个数量级。尤其在大规模稀疏矩阵上，真正耗时的是内存带宽而不是计算量，所以但凡是自己实现求解器，尽量用成熟的稀疏线性代数库，比如Eigen、PETSc、SuiteSparse，而不是自己写循环。

如果不得不自己实现，一个重要原则是矩阵向量乘要保证对稀疏非零元素顺序的连续访问。比如CSR行压缩格式下，内层循环应该按列索引扫描非零元，避免随机跳跃落到内存里。

6.5 奇异与半正定矩阵的处理

Krylov方法大多默认矩阵非奇异。半正定或奇异矩阵会导致CG在某个搜索方向上出现零除，GMRES的最小二乘问题也可能变得病态。工程上如果遇到这类问题，常见手段是先加一个小正则化项 (\epsilon I)，把问题转化为良态问题；或者在物理建模阶段找出奇异性的源头（比如浮体刚体模态），显式约束掉零空间。

从课程角度说，这部分会让初学同学觉得超纲，但实际工程里几乎一定会遇到，提前在代码里留好检查逻辑，能省很多善后时间。

最后再讲一点个人体会。高等数值分析这门课教Krylov子空间方法，我觉得最重要不是背下每个算法的迭代公式，而是建立起“把大型问题投影到低维子空间来观察”的思维方式。GMRES和CG是这套思维方式最经典的两个落点，后面遇到特征值问题、矩阵方程、模型降阶，你都会看到同一个思想在不同场景下反复出现。如果这篇文章能让你在读完以后，拿到一个陌生的大规模稀疏矩阵时愿意先想想“它的谱分布可能是什么样、应该选哪类Krylov方法、预处理怎么搭”，那写它的目的就达到了。