从大二数值分析课那个学期开始我就对数线性方程组的解法很着迷。当时课本里最核心的就是高斯消元、LU分解解法完完整整矩阵一扔进去一顿操作出一个“精确解”。后来真正用有限差分、有限元去解偏微分方程才意识到现实世界里的矩阵有多大三维网格一加密未知量轻松到百万级直接法那套“保存在内存里慢慢消元”的思路直接崩。于是Krylov子空间方法成了绕不开的主角。这篇文章我尽量把高等数值分析里关于Krylov子空间方法的核心逻辑讲清楚不只是摆一遍CG和GMRES的公式而是把“为什么需要子空间”“每个算法在优化什么”“收敛性为什么波动”“预处理为什么不可或缺”这些底层问题一次说明白。适合正在学数值分析、要写有限元求解器或者被大型稀疏矩阵折磨的工程师参考。1. 直接法在大规模问题上的天花板在哪里1.1 高斯消元的存储噩梦与填充现象先从一个具体例子看起。求解二维泊松方程 (-u_{xx}-u_{yy}f) 时用五点差分格式离散(N\times N) 的内部网格点对应 (nN^2) 个未知量离散矩阵是三对角块加边结构稀疏性非常好每行只有5个非零元。按n10^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)) 增长。换句话说你为了精确解加密网格没想到迭代次数反而暴涨最终整体耗时根本不比直接法省心。这类方法适合作为光滑器或预处理子却扛不起求解器的大梁。这个瓶颈背后的原因很有意思经典迭代法每一轮只是在“局部修正误差”它通过矩阵分裂 (AD-L-U) 将问题转化成不动点迭代步与步之间没有任何全局信息。我们要的下一步迭代修正量是整个空间里让残差下降最多的方向这需要从矩阵全局提取信息。1.3 Krylov线索的出现幂向量里藏着解的信息现在我们换个视角。给定一个初始近似解 (x_0)可以得到初始残差 (r_0b-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) 可以写成误差形式 (xx_0z)则残差满足 (r_0Az)。当矩阵可逆时 (zA^{-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_1r_0/|r_0|)然后对每一个 (k1,2,\ldots,m)计算 (wA v_k)把 (w) 在之前所有基底上的投影减掉再归一化得到 (v_{k1})。整个过程写成矩阵关系非常漂亮[ A V_m V_m H_m h_{m1,m}v_{m1}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, m1)) H np.zeros((m1, 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(j1): H[i, j] np.dot(V[:, i], w) w - H[i, j] * V[:, i] H[j1, j] np.linalg.norm(w) if H[j1, j] 1e-15: # 这是happy breakdown子空间已经不变可以提前停机 break V[:, j1] w / H[j1, j] return V, H这个代码里最容易被忽略的一个数值细节是教材里写经典Gram-Schmidt是三行伪代码但直接这样写在大规模计算中会因为舍入误差丢失正交性导致整个Arnoldi过程失真。经验做法是用修正Gram-SchmidtMGS也就是每减掉一个投影分量立即对剩余向量做归一化前的修正。2.3 理解Hessenberg矩阵的“扁平投影”很多初学者困惑Arnoldi过程跑了一堆正交化最后得到的H_m为什么是Hessenberg结构道理其实在他的正交化策略里每步只拿新向量 (wA v_k) 与已经生成的前 (k) 个基底做正交化减去这些投影之后剩下分量与 (v_1,\ldots,v_k) 都正交它唯一可能与下一个新基底 (v_{k1}) 对齐。因此在H_m的列上第 (k) 列只有前 (k1) 行非零下标恰好呈阶梯型。这个结构价值很大。它意味着我们可以用递推方式一点一点扩展子空间而不必每次重新生成整组基底。后面用GMRES求解残差最小化问题时H_m的Hessenberg结构也能让最小二乘子问题保持在很低的计算复杂度。3. GMRES与CG两者真正差异在于优化目标不同3.1 GMRES在Krylov子空间里极小化残差范数有了Arnoldi给出的正交基底 (V_m)我们自然想问在这个m维子空间里哪个向量让当前残差范数最小既然解限定为 (xx_0V_m y)那么新残差为[ r_m r_0 - A V_m y \beta v_1 - V_{m1}\tilde{H}_m y ]其中 (\beta|r_0|)(\tilde{H}m) 是包含最后一行的 ((m1)\times m) 上Hessenberg矩阵。由于 (V{m1}) 列正交极小化 (|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, x0None, tol1e-10, max_iter1000): 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这段代码里最反直觉的一行是方向更新 (pr\beta p)。为什么每次新的搜索方向要加上旧方向的修正因为单纯沿负梯度方向走会重复做无用功加上旧方向修正后能保证每个新搜索方向是关于A共轭的也就是满足 (p_i^T A p_j0)这样理论上有限步就能在精确算术下收敛到n步内。3.3 其他Krylov家族成员与适用边界一旦理解GMRES是“无脑对残差求最小”CG是“对称正定下的能量最优”就能理解整个Krylov方法家族为什么长成这副模样。方法适用矩阵每步存储优化的对象典型场景CG对称正定O(n)能量范数误差椭圆型PDE、结构力学MINRES对称不定O(n)残差范数鞍点问题、约束优化GMRES一般非奇异O(mn)残差范数非对称对流扩散BiCGSTAB一般非奇异O(n)残差双正交条件非对称问题兼顾内存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^TA^TA)可以用谱分布给出类似上界。但对非正规矩阵仅凭特征值分布无法解释实际收敛行为伪谱pseudospectrum比谱更能说明问题。我见过工程里有些非对称矩阵特征值全在右半平面看起来人畜无害实际GMRES迭代几十步残差纹丝不动这正是伪谱效应在作怪。所以对非对称问题我通常建议先跑几十步GMRES看残差曲线再判断是继续硬跑、换预处理还是改用其他算法族。纸上算收敛阶在非对称领域参考价值远低于对称正定领域。4.3 超线性收敛和停滞期残差曲线的“心跳”实际观察CG或GMRES的残差曲线经常发现它不是一条平滑下降的直线而是先慢吞吞下降、甚至中途出现一个平台期然后突然加速。这个现象被称作超线性收敛。原因是Krylov多项式在前若干步还没有“认全”特征值的信息对某些特征值方向的衰减很弱随着迭代继续多项式在这些方向上的取值被逐渐压低整体收敛速度就会变快。这一点在日常调试中有直接指导意义看到残差曲线进入平台不一定要立刻判定方法失效可以观察平台长度。如果平台出现在迭代早期多迭代几步可能会突降如果平台出现在重启动后的每一轮说明重启动m太小Krylov子空间根本来不及积累足够信息这时候应当增大m。4.4 理想中断与数值边界还有个现象叫happy breakdown在Arnoldi或Lanczos过程中表现为 (h_{m1,m}) 恰好为零。这意味着当前Krylov子空间在A的作用下不再向外扩展子空间已经是A的不变子空间此时如果 (r_0\in K_m) 且矩阵在这个子空间上非奇异那么GMRES得到的解将直接是精确解。这个“意外之喜”判断代码里很简单一旦当前残差小于阈值立即退出循环避免无意义的后续迭代。不过实际计算中更常遇到的是 (h_{m1,m}) 非常小但不为零比如降到1e-14量级。此时再除以它归一化v_{m1}会让数值噪声被放大后续基底质量急剧恶化。稳妥的做法是引入软阈值如果 (h_{m1,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}AxM^{-1}b)。对对角占优矩阵效果立竿见影代码只需要一行几乎零风险。但它的能力上限很低对强耦合的PDE离散问题帮助有限。再往上走是ILU不完全LU分解族。它模拟直接法的消元过程但严格控制填充ILU(0)只保留原矩阵非零结构上的因子ILUT基于阈值丢弃小于阈值的填充元。这类预处理在实际工程里非常有价值特别是非对称问题。一个关键经验是ILU预处理的有效性强烈依赖矩阵排序用反向Cuthill-McKeeRCM重排序后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}AxM^{-1}b ]虽然理论上和右预处理[ AM^{-1}yb,\qquad xM^{-1}y ]有相近的收敛性质但注意左预处理改变了残差定义。GMRES在最小化残差范数时实际上最小化的是 (|M^{-1}(b-Ax)|)这个“加权残差”不同于原始残差。好处是如果 (M^{-1}A) 谱分布越好收敛越快坏处是如果M构造得很差加权后的收敛标准会和真实残差脱节。右预处理没有这个问题因为GMRES可以保持对原始残差 (|b-Ax|) 的监控实现上更加透明。所以我的默认选择是右预处理除非有特殊需要才会用左预处理版本。5.4 一个直观的小实验为了说明预处理的效果我经常在课堂上演示这样一个二维Poisson例子100×100网格用五点差分生成SPD矩阵A分别统计未预处理CG与Jacobi预处理CGPCG的迭代次数。未预处理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不是越大越好也不是越小越省。太小比如m10对很多非对称问题基本会陷进周期性的停滞每轮重启后残差只下降一点点甚至原地踏步。太大每步正交化要跟所有历史向量算内积成本呈二次方增长。我的实验习惯是先用m30跑一个算例输出残差曲线观察。如果曲线在接近收敛时突然“断崖”说明m刚好够用如果残差在某个高位反复震荡就依次尝试50、80、100。大部分非对称工程问题m50到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方法、预处理怎么搭”那写它的目的就达到了。