简介这是数值线性代数领域关于Krylov子空间迭代算法的教学讲稿面向需要求解大型稀疏线性方程组的本科生、研究生及科研人员。内容从子空间迭代思想切入给出Krylov子空间定义详细推导Arnoldi过程与Lanczos过程的矩阵表达和三项递推公式并系统讲解GMRES算法与共轭梯度法的构造原理、最优性条件及收敛性分析。资源为1个PDF文档大小364KB结构紧凑、公式清晰适合作为课堂讲义或自学笔记使用。目前已有610人学习读者可借此深入理解Krylov子空间方法的数学基础与算法实现框架。1. 第七讲Krylov 子空间迭代算法到底在迭代什么如果你手头有一个 100 万阶的稀疏线性方程组直接做 LU 分解填充量会把内存撑爆用经典 Jacobi 迭代三千步下去残差还赖在 1e-2 不肯走。这时候 Krylov 子空间迭代算法是少数几个既不需要显式存储矩阵因子、又能在可接受步数内收敛的选择。它的核心逻辑和教科书里教的“迭代逼近解”不太一样不是逐分量修正而是在一个不断扩张的子空间里寻找最优近似解每步只做矩阵向量乘MatVec矩阵本身可以是隐式的甚至不需要能访问到元素。这一讲把 Krylov 子空间迭代算法拆成三层来讲先建立子空间逼近的数学直觉再比较四种主流算法在什么条件下该选谁最后落到底层实现里的预处理子、停机准则和参数调优。适合正在写有限元或 CFD 求解器、被收敛曲线折磨的工程师也适合想把 GMRES 的“重启”和“预处理”一次弄明白的研究生。文中所有代码基于 Python 生态但结论和参数设置在 C/Fortran 实现里同样成立。2. 从投影角度看 Krylov 子空间迭代算法为什么是多项式逼近2.1 用多项式逼近解释子空间扩张的本质Krylov 子空间迭代算法的起点很朴素。对Ax b取初始解x0定义初始残差r0 b - Ax0那么第 k 步的搜索空间定义为K_k(A, r0) span{ r0, A r0, A^2 r0, ..., A^(k-1) r0 }这里的关键理解是A^i r0不是一个需要显式构造的向量序列而是每一步把上一次 MatVec 的结果再乘一次 A。于是x_k一定可以写成x_k x0 q_{k-1}(A) r0其中q_{k-1}是一个次数不超过 k-1 的多项式。对应的残差是r_k p_k(A) r0p_k是另一个多项式且满足p_k(0) 1。所以 Krylov 子空间迭代算法的本质是把“解线性方程组”转化为“在所有满足p_k(0)1的多项式里找一个让||p_k(A) r0||最小的 p_k”。这个视角解释了为什么算法族叫“子空间迭代”——每一步只是让多项式次数加一。迭代步数上限由极小多项式次数决定。如果极小多项式次数是 m那么理论上最多 m 步必然停机。但实际计算中舍入误差导致极小多项式性质丢失所以步数经常超过理论值这也直接引出后面要讲的预处理和重启策略。2.2 三种投影方式决定算法分类Krylov 子空间迭代算法的具体实现取决于约束条件——也就是“在子空间里找近似解”的最优性准则。常见有三类第一类是 Galerkin 条件要求残差与搜索子空间正交r_k ⊥ K_k。CG 和 FOM 属于这一类。第二类是极小残差条件直接最小化||b - Ax||_2在 Krylov 子空间上做最小二乘。GMRES、MINRES 属于这一类。第三类是 Petrov-Galerkin 条件残差与另一个辅助子空间正交BiCG 和 QMR 走这条路线。算法适用矩阵投影方式每步存储每步代价CG对称正定GalerkinO(n)1 次 MatVecMINRES对称不定极小残差O(n)1 次 MatVecGMRES非对称极小残差O(n·k)1 次 MatVecBiCGSTAB非对称Petrov-GalerkinO(n)2 次 MatVec选型逻辑很简单对称正定用 CG对称不定用 MINRES因为此时 CG 可能除以零或收敛不稳定非对称且能承受存储增长用 GMRES非对称且矩阵规模大到存不下 Krylov 基用 BiCGSTAB 这类短递推方法。这里每步代价和存储必须是选型时的第一道筛子。2.3 快速检验收敛性的三个矩阵性质在动手实现前先用三个性质预估这个矩阵适不适合 Krylov 子空间迭代算法。第一个是条件数。cond(A)直接决定 CG 的收敛速度上界理论上步数与sqrt(cond)成正比矩阵越病态收敛越慢。第二个是特征值分布。特征值越聚集尤其是不含零的紧簇多项式逼近越容易构造收敛越快。第三个是正规性。非正规矩阵即使特征值分布很好看收敛曲线也会出现“先停滞、后骤降”的假象上界估计可能完全失效。我一般会先用eigvalsh或eigs快速算一下极端特征值再跑 20 步不带预处理的迭代看残差下降趋势。如果 20 步残差只降了一个数量级说明矩阵病态或特征值散直接跳到预处理章节。3. 四种 Krylov 子空间迭代算法的内存特征与选型边界3.1 CG 的局部性优势与破裂风险CG 是内存效率最高的 Krylov 子空间迭代算法只需要保存当前解、残差和搜索方向各一个向量每步一次 MatVec。它能在对称正定矩阵上保证残差范数单调不增。但 CG 的优势也来自对矩阵性质的强依赖一旦矩阵不对称或不定可能在迭代中直接除以零术语叫“破裂”或者收敛曲线出现剧烈震荡。工程上最常见的误用是把 CG 用在非对称问题上。对此我一般这样处理先检查A - A^T的范数是否明显大于零如果是直接考虑 GMRES 或 BiCGSTAB不要在 CG 上浪费调参时间。另一个容易忽略的问题是舍入误差下的“延迟收敛”——理论上最多 n 步收敛但浮点环境下极小多项式性质被破坏实际可能需要 2n 到 3n 步。CG 的稳定实现有一个关键细节避免直接计算x_k而是通过递推更新。这种方式能减少一次矩阵向量乘的误差累积同时为后面积累误差分析留出余地。3.2 GMRES 的存储墙和重启阈值GMRES 的收敛性在理论上是最稳的每一步都保证残差范数不增极小残差性质使得它几乎不会发散。但代价是必须保存全部 Krylov 基向量来做最小二乘k 步后存储量为 O(n·k)。当矩阵规模到百万阶、迭代到上千步时内存占用会失控。重启是解决存储问题的标准手段。重启后的 GMRES(m) 每 m 步清空基向量重新开始但代价是失去全局最优性收敛可能停滞。我一般把重启阈值设在 30 到 50 之间。设大则内存压力大设小则收敛变慢。更好的做法是先跑一次不重启的 GMRES观察“残差显著下降需要多少步”然后把这个值乘 1.5 作为重启阈值。GMRES 的另一个隐藏成本是每步的 Arnoldi 正交化。Modified Gram-Schmidt 是经典选择但对大规模并行计算要小心其同步开销如果有 GPU 或 MPI 环境可以考虑使用 TSQR 或 Cholesky QR 替代。3.3 BiCGSTAB 的低内存优势与“伪收敛”陷阱BiCGSTAB 用两组双正交基做短递推每步两次 MatVec内存占用 O(n)。它最大的问题是对舍入误差敏感容易出现“伪收敛”残差范数已经降到 1e-12但实际误差||x - x_true||还停留在 1e-3。原因是双正交过程在舍入误差下会丢失正交性导致残差和误差脱钩。规避方法有两个层面。一是从算法参数入手采用稳定变体 BiCGSTAB(l) 或 BiCGSTAB2这些变体通过周期性地对残量进行额外修正来恢复精度。二是从工程层面入手把停机准则从“只看残差”改为“残差与解变化量同时满足条件”。如果最终要的是高精度解我建议用 BiCGSTAB 快速迭代到 1e-6 附近再切换到 GMRES 做最后几步精化。3.4 选型决策表与测试脚本用下面这段 Python 脚本可以一次性对比四种算法对同一矩阵的收敛行为import numpy as np from scipy.sparse.linalg import cg, gmres, minres, bicgstab n 2000 A np.random.randn(n, n) * 0.01 A A A.T n * np.eye(n) # 对称正定条件数可控 b np.random.randn(n) # 关闭 scipy 内部的预处理和重启观察原始算法表现 x_cg, info_cg cg(A, b, rtol1e-8, maxiter500, atol0) x_gmres, info_gmres gmres(A, b, rtol1e-8, maxiter500, atol0) x_minres, info_minres minres(A, b, rtol1e-8, maxiter500, atol0) x_bicg, info_bicg bicgstab(A, b, rtol1e-8, maxiter500, atol0) for name, info in [(CG, info_cg), (GMRES, info_gmres), (MINRES, info_minres), (BiCGSTAB, info_bicg)]: print(f{name}: converged{info 0})info返回值是理解 scipy 迭代器行为的关键info0表示收敛info0表示达到最大迭代步数未收敛info0表示输入非法或算法破裂。测试对称正定矩阵时CG 和 MINRES 都应该在几十步内收敛GMRES 会收敛但每步存储增长BiCGSTAB 慢一些。换用非对称矩阵后CG 和 MINRES 会失败GMRES 和 BiCGSTAB 则继续工作。对比这些行为能帮你建立对算法边界的直觉。4. 让 Krylov 子空间迭代算法真正快的预处理实现4.1 预处理子的选择顺序从对角到不完全分解预处理是 Krylov 子空间迭代算法从“能跑”到“跑得快”的分水岭。核心思想是对M^{-1}Ax M^{-1}b做迭代其中 M 是对 A 的近似。M 越接近 A预处理后的矩阵特征值越聚集收敛越快但 M 的构造和每次应用代价也越高。这条权衡曲线是所有预处理选择的出发点。预处理子的选择顺序从廉价到昂贵依次为Jacobi仅对角、SSOR、ILU(0)、ILU(k)、多水平方法。我的经验法是条件数在 1e3 以内Jacobi 足够到 1e5 量级SSOR 或 ILU(0) 是首选超过 1e7必须上多水平预处理或领域分解配合。M^{-1}不需要显式构造只需要能在每步迭代中计算M^{-1}v。这意味着预处理子的实现复杂度与 Krylov 迭代本身解耦你可以把任意预处理逻辑封装成黑盒。4.2 ILU(0) 的填充阈值与实现要点ILU(0) 是最常用的通用预处理子——它做 LU 分解但强制保持与 A 相同的稀疏模式不产生任何填充元素。实现上关键点在于不填充导致近似精度有限所以 ILU(0) 适合中等问题对更难的矩阵需要 ILU(k) 允许有限层填充或 ILU(t) 按数值大小淘汰小元素。import scipy.sparse.linalg as spla # 构造 ILU(0) 预处理子对象 ilu spla.spilu(A, drop_tol1e-4, fill_factor10) # 把预处理子包装成 M^{-1} 的调用形式 def apply_precond(v): return ilu.solve(v) # 在 GMRES 中传入预处理函数 x, info spla.gmres(A, b, rtol1e-10, maxiter200, Mapply_precond)这里的drop_tol控制丢弃小元素的门槛fill_factor控制允许的填充量上限。两者共同决定预处理子的质量和构建成本。实际调参时先固定drop_tol1e-4跑一遍观察迭代步数如果步数降不下来把fill_factor从 5 提到 15如果构建时间过久则增大drop_tol减少填充。这个试参顺序比直接乱试更高效。4.3 块预处理与并行化当矩阵来自多物理场耦合或多组分工况时标量预处理子往往不充分。块预处理把 A 分块只对对角块做不完全分解块间耦合在迭代中处理。这类预处理在流体力学和电磁场模拟中表现突出实现上可以直接用scipy.linalg.lu_factor对每个块做密分解或用spla.spilu对每个块做稀疏不完全分解。并行环境下要特别注意ILU 的求解过程是串行依赖的大规模并行时可能成为瓶颈。此时考虑三类替代方案多项式预处理用 A 的多项式近似M^{-1}只用 MatVec天然并行、多色 SSOR按图着色重排让同一颜色的点互不依赖可并行、或用域分解法把问题切成子域各算各的。5. 收敛性诊断如何在 Krylov 子空间迭代算法中定位停滞和破裂5.1 看对的量残差、相对残差与 A 范数Krylov 子空间迭代算法的停机准则设置不当会导致两种情况过早停止得到错误解或过晚停止浪费算力。需要区分三个量真实残差||b - Ax_k||、递推残差迭代过程中通过递推隐式维持的量、以及误差||x_k - x_true||。在浮点环境下递推残差可能远小于真实残差这是所有 Krylov 迭代共有的问题。我建议在迭代中同时输出真实残差和相对残差并在接近停机阈值时切换到真实残差验证。对于 CG 还有专门的 A 范数——||x_k - x_true||_A——它直接度量能量误差比二范数更适合同类问题。5.2 停滞的三种模式与各自对策Krylov 收敛曲线出现停滞通常是三种模式之一。第一种是“平台期”残差在一个水平持久不变然后突然下降这常见于非正规矩阵GMRES 理论上不该出现但有限精度下仍可能发生。对策是检查是否因为重启导致丢失了关键方向或者换更高质量的预处理子。第二种是“渐进爬坡”残差持续但极其缓慢下降这是预处理不足的典型信号应该考虑加强预处理而不是增大迭代步数。第三种是“发散振荡”残差上下波动不收敛通常是矩阵非正规、预处理不稳定或迭代格式不适合问题类型。区分这三种模式最快的方法是同时画半对数坐标下的残差曲线和每步的||x_k - x_{k-1}||。如果解的变化量已经小到机器精度但残差还在高位问题出在预处理上如果解变化量很大但残差不降问题可能出在算法的正交性丢失上。5.3 参数调整的最短路径面对不收敛的迭代很多人的第一反应是调大maxiter——但这通常是无效的。正确的调整路径是先看矩阵条件数和特征值分布确认问题是否出在矩阵本身的病态性然后依次尝试从无预处理到 Jacobi、再到 ILU(0)、再到 ILU(k)观察每步收敛步数的变化率如果 ILU 没有带来数量级改善再考虑换算法族比如把 GMRES 换成 BiCGSTAB(l) 或重启 GMRES。一个实用的做法是把预处理质量作为第一优先级。ILU 的参数选择优先于迭代算法的参数选择。迭代参数中rtol的设定要考虑最终用途做特征值问题的内迭代rtol1e-3也许就够做高精度结构分析可能需要rtol1e-12。不要在需求未知时盲目追求高精度这会平白增加几十步迭代。6. 病态矩阵处理与 Krylov 子空间迭代算法的规模化实践6.1 矩阵重排序Cuthill-McKee 对 ILU 的影响预处理子的质量取决于矩阵结构而结构可以通过重排改变。Reverse Cuthill-McKeeRCM是一种带宽缩减算法它把矩阵重新排序成轮廓更窄的形式。对于有限元网格生成的矩阵RCM 重排后 ILU 的填充更集中在对角线附近不完全分解的丢弃误差更小预处理质量显著提升。from scipy.sparse.csgraph import reverse_cuthill_mckee from scipy.sparse import csr_matrix # 假设 A 是 CSR 格式的稀疏矩阵 A_csr csr_matrix(A) perm reverse_cuthill_mckee(A_csr, symmetric_modeTrue) # 重排矩阵和右端项 A_rcm A_csr[perm][:, perm] b_rcm b[perm] # 在重排后的系统上做迭代最后把解映射回原次序 x_rcm, info spla.gmres(A_rcm, b_rcm, rtol1e-10, Milu_precond) x np.empty_like(b) x[perm] x_rcmperm是重排后的下标映射数组A_csr[perm][:, perm]同时对行和列重排保持对称性。RCM 不改变特征值所以 Krylov 收敛的理论性质不变但它改变了预处理子可用的结构信息实际效果往往是迭代步数下降一半以上。值得注意一个容易出错的地方解必须映射回原顺序才能和其他模块对接忘记这一步是接错数据的高频原因。6.2 特征值平移对大规模问题收敛性的改善A的零特征值或接近零的特征值是 Krylov 子空间迭代算法收敛的最大障碍。对平移后的系统(A σI)x b σx做迭代可以改善特征值聚集性但代价是需要额外处理 x 项。更常见的是把它作为预处理手段把平移合并到预处理子中M A σI的近似。这记录了“越接近零的特征值对 Krylov 方法的危害越大”的工程直觉。对 Near-null 空间明显的矩阵比如结构分析中的刚体模态平移量 σ 取一个比最小非零特征值小一两个数量级的数值效果最好。6.3 与多层方法的分工配合当 Krylov 子空间迭代算法在千万自由度级别的结构分析中遭遇收敛瓶颈时可以考虑与多层方法分工协作。底层思路是Krylov 负责处理平滑的高频误差粗网格修正负责处理低频误差。在规模化实践中比较经典的方案是“多层预处理 Krylov 加速”的组合模式——用多重网格或领域分解做预处理子用 GMRES 或 CG 驱动整体收敛。参数设置的经验法则是多层方法的层数取决于网格细化层次和特征值分布Krylov 侧的重启阈值通常设为 30 或 50如果配合代数多重网格预处理CG 通常在 20 步以内收敛。这个组合也是当前主流商用有限元软件在大规模分析时的默认选择。最后提一个容易被忽略的验证细节无论预处理多复杂最终都要额外做一次真实残差检查。因为经过重排、平移和预处理之后递推残差与真实残差的差距可能被放大只有直接计算||b - Ax||才能确认你的 Krylov 子空间迭代算法真的收敛到了该有的精度。本文还有配套的精品资源点击获取 
长河君 · 2026/9/17 21:31:40 · 阅读 1.2k Krylov子空间迭代算法全解析：从投影原理到预处理实战
简介这是数值线性代数领域关于Krylov子空间迭代算法的教学讲稿面向需要求解大型稀疏线性方程组的本科生、研究生及科研人员。内容从子空间迭代思想切入给出Krylov子空间定义详细推导Arnoldi过程与Lanczos过程的矩阵表达和三项递推公式并系统讲解GMRES算法与共轭梯度法的构造原理、最优性条件及收敛性分析。资源为1个PDF文档大小364KB结构紧凑、公式清晰适合作为课堂讲义或自学笔记使用。目前已有610人学习读者可借此深入理解Krylov子空间方法的数学基础与算法实现框架。1. 第七讲Krylov 子空间迭代算法到底在迭代什么如果你手头有一个 100 万阶的稀疏线性方程组直接做 LU 分解填充量会把内存撑爆用经典 Jacobi 迭代三千步下去残差还赖在 1e-2 不肯走。这时候 Krylov 子空间迭代算法是少数几个既不需要显式存储矩阵因子、又能在可接受步数内收敛的选择。它的核心逻辑和教科书里教的“迭代逼近解”不太一样不是逐分量修正而是在一个不断扩张的子空间里寻找最优近似解每步只做矩阵向量乘MatVec矩阵本身可以是隐式的甚至不需要能访问到元素。这一讲把 Krylov 子空间迭代算法拆成三层来讲先建立子空间逼近的数学直觉再比较四种主流算法在什么条件下该选谁最后落到底层实现里的预处理子、停机准则和参数调优。适合正在写有限元或 CFD 求解器、被收敛曲线折磨的工程师也适合想把 GMRES 的“重启”和“预处理”一次弄明白的研究生。文中所有代码基于 Python 生态但结论和参数设置在 C/Fortran 实现里同样成立。2. 从投影角度看 Krylov 子空间迭代算法为什么是多项式逼近2.1 用多项式逼近解释子空间扩张的本质Krylov 子空间迭代算法的起点很朴素。对Ax b取初始解x0定义初始残差r0 b - Ax0那么第 k 步的搜索空间定义为K_k(A, r0) span{ r0, A r0, A^2 r0, ..., A^(k-1) r0 }这里的关键理解是A^i r0不是一个需要显式构造的向量序列而是每一步把上一次 MatVec 的结果再乘一次 A。于是x_k一定可以写成x_k x0 q_{k-1}(A) r0其中q_{k-1}是一个次数不超过 k-1 的多项式。对应的残差是r_k p_k(A) r0p_k是另一个多项式且满足p_k(0) 1。所以 Krylov 子空间迭代算法的本质是把“解线性方程组”转化为“在所有满足p_k(0)1的多项式里找一个让||p_k(A) r0||最小的 p_k”。这个视角解释了为什么算法族叫“子空间迭代”——每一步只是让多项式次数加一。迭代步数上限由极小多项式次数决定。如果极小多项式次数是 m那么理论上最多 m 步必然停机。但实际计算中舍入误差导致极小多项式性质丢失所以步数经常超过理论值这也直接引出后面要讲的预处理和重启策略。2.2 三种投影方式决定算法分类Krylov 子空间迭代算法的具体实现取决于约束条件——也就是“在子空间里找近似解”的最优性准则。常见有三类第一类是 Galerkin 条件要求残差与搜索子空间正交r_k ⊥ K_k。CG 和 FOM 属于这一类。第二类是极小残差条件直接最小化||b - Ax||_2在 Krylov 子空间上做最小二乘。GMRES、MINRES 属于这一类。第三类是 Petrov-Galerkin 条件残差与另一个辅助子空间正交BiCG 和 QMR 走这条路线。算法适用矩阵投影方式每步存储每步代价CG对称正定GalerkinO(n)1 次 MatVecMINRES对称不定极小残差O(n)1 次 MatVecGMRES非对称极小残差O(n·k)1 次 MatVecBiCGSTAB非对称Petrov-GalerkinO(n)2 次 MatVec选型逻辑很简单对称正定用 CG对称不定用 MINRES因为此时 CG 可能除以零或收敛不稳定非对称且能承受存储增长用 GMRES非对称且矩阵规模大到存不下 Krylov 基用 BiCGSTAB 这类短递推方法。这里每步代价和存储必须是选型时的第一道筛子。2.3 快速检验收敛性的三个矩阵性质在动手实现前先用三个性质预估这个矩阵适不适合 Krylov 子空间迭代算法。第一个是条件数。cond(A)直接决定 CG 的收敛速度上界理论上步数与sqrt(cond)成正比矩阵越病态收敛越慢。第二个是特征值分布。特征值越聚集尤其是不含零的紧簇多项式逼近越容易构造收敛越快。第三个是正规性。非正规矩阵即使特征值分布很好看收敛曲线也会出现“先停滞、后骤降”的假象上界估计可能完全失效。我一般会先用eigvalsh或eigs快速算一下极端特征值再跑 20 步不带预处理的迭代看残差下降趋势。如果 20 步残差只降了一个数量级说明矩阵病态或特征值散直接跳到预处理章节。3. 四种 Krylov 子空间迭代算法的内存特征与选型边界3.1 CG 的局部性优势与破裂风险CG 是内存效率最高的 Krylov 子空间迭代算法只需要保存当前解、残差和搜索方向各一个向量每步一次 MatVec。它能在对称正定矩阵上保证残差范数单调不增。但 CG 的优势也来自对矩阵性质的强依赖一旦矩阵不对称或不定可能在迭代中直接除以零术语叫“破裂”或者收敛曲线出现剧烈震荡。工程上最常见的误用是把 CG 用在非对称问题上。对此我一般这样处理先检查A - A^T的范数是否明显大于零如果是直接考虑 GMRES 或 BiCGSTAB不要在 CG 上浪费调参时间。另一个容易忽略的问题是舍入误差下的“延迟收敛”——理论上最多 n 步收敛但浮点环境下极小多项式性质被破坏实际可能需要 2n 到 3n 步。CG 的稳定实现有一个关键细节避免直接计算x_k而是通过递推更新。这种方式能减少一次矩阵向量乘的误差累积同时为后面积累误差分析留出余地。3.2 GMRES 的存储墙和重启阈值GMRES 的收敛性在理论上是最稳的每一步都保证残差范数不增极小残差性质使得它几乎不会发散。但代价是必须保存全部 Krylov 基向量来做最小二乘k 步后存储量为 O(n·k)。当矩阵规模到百万阶、迭代到上千步时内存占用会失控。重启是解决存储问题的标准手段。重启后的 GMRES(m) 每 m 步清空基向量重新开始但代价是失去全局最优性收敛可能停滞。我一般把重启阈值设在 30 到 50 之间。设大则内存压力大设小则收敛变慢。更好的做法是先跑一次不重启的 GMRES观察“残差显著下降需要多少步”然后把这个值乘 1.5 作为重启阈值。GMRES 的另一个隐藏成本是每步的 Arnoldi 正交化。Modified Gram-Schmidt 是经典选择但对大规模并行计算要小心其同步开销如果有 GPU 或 MPI 环境可以考虑使用 TSQR 或 Cholesky QR 替代。3.3 BiCGSTAB 的低内存优势与“伪收敛”陷阱BiCGSTAB 用两组双正交基做短递推每步两次 MatVec内存占用 O(n)。它最大的问题是对舍入误差敏感容易出现“伪收敛”残差范数已经降到 1e-12但实际误差||x - x_true||还停留在 1e-3。原因是双正交过程在舍入误差下会丢失正交性导致残差和误差脱钩。规避方法有两个层面。一是从算法参数入手采用稳定变体 BiCGSTAB(l) 或 BiCGSTAB2这些变体通过周期性地对残量进行额外修正来恢复精度。二是从工程层面入手把停机准则从“只看残差”改为“残差与解变化量同时满足条件”。如果最终要的是高精度解我建议用 BiCGSTAB 快速迭代到 1e-6 附近再切换到 GMRES 做最后几步精化。3.4 选型决策表与测试脚本用下面这段 Python 脚本可以一次性对比四种算法对同一矩阵的收敛行为import numpy as np from scipy.sparse.linalg import cg, gmres, minres, bicgstab n 2000 A np.random.randn(n, n) * 0.01 A A A.T n * np.eye(n) # 对称正定条件数可控 b np.random.randn(n) # 关闭 scipy 内部的预处理和重启观察原始算法表现 x_cg, info_cg cg(A, b, rtol1e-8, maxiter500, atol0) x_gmres, info_gmres gmres(A, b, rtol1e-8, maxiter500, atol0) x_minres, info_minres minres(A, b, rtol1e-8, maxiter500, atol0) x_bicg, info_bicg bicgstab(A, b, rtol1e-8, maxiter500, atol0) for name, info in [(CG, info_cg), (GMRES, info_gmres), (MINRES, info_minres), (BiCGSTAB, info_bicg)]: print(f{name}: converged{info 0})info返回值是理解 scipy 迭代器行为的关键info0表示收敛info0表示达到最大迭代步数未收敛info0表示输入非法或算法破裂。测试对称正定矩阵时CG 和 MINRES 都应该在几十步内收敛GMRES 会收敛但每步存储增长BiCGSTAB 慢一些。换用非对称矩阵后CG 和 MINRES 会失败GMRES 和 BiCGSTAB 则继续工作。对比这些行为能帮你建立对算法边界的直觉。4. 让 Krylov 子空间迭代算法真正快的预处理实现4.1 预处理子的选择顺序从对角到不完全分解预处理是 Krylov 子空间迭代算法从“能跑”到“跑得快”的分水岭。核心思想是对M^{-1}Ax M^{-1}b做迭代其中 M 是对 A 的近似。M 越接近 A预处理后的矩阵特征值越聚集收敛越快但 M 的构造和每次应用代价也越高。这条权衡曲线是所有预处理选择的出发点。预处理子的选择顺序从廉价到昂贵依次为Jacobi仅对角、SSOR、ILU(0)、ILU(k)、多水平方法。我的经验法是条件数在 1e3 以内Jacobi 足够到 1e5 量级SSOR 或 ILU(0) 是首选超过 1e7必须上多水平预处理或领域分解配合。M^{-1}不需要显式构造只需要能在每步迭代中计算M^{-1}v。这意味着预处理子的实现复杂度与 Krylov 迭代本身解耦你可以把任意预处理逻辑封装成黑盒。4.2 ILU(0) 的填充阈值与实现要点ILU(0) 是最常用的通用预处理子——它做 LU 分解但强制保持与 A 相同的稀疏模式不产生任何填充元素。实现上关键点在于不填充导致近似精度有限所以 ILU(0) 适合中等问题对更难的矩阵需要 ILU(k) 允许有限层填充或 ILU(t) 按数值大小淘汰小元素。import scipy.sparse.linalg as spla # 构造 ILU(0) 预处理子对象 ilu spla.spilu(A, drop_tol1e-4, fill_factor10) # 把预处理子包装成 M^{-1} 的调用形式 def apply_precond(v): return ilu.solve(v) # 在 GMRES 中传入预处理函数 x, info spla.gmres(A, b, rtol1e-10, maxiter200, Mapply_precond)这里的drop_tol控制丢弃小元素的门槛fill_factor控制允许的填充量上限。两者共同决定预处理子的质量和构建成本。实际调参时先固定drop_tol1e-4跑一遍观察迭代步数如果步数降不下来把fill_factor从 5 提到 15如果构建时间过久则增大drop_tol减少填充。这个试参顺序比直接乱试更高效。4.3 块预处理与并行化当矩阵来自多物理场耦合或多组分工况时标量预处理子往往不充分。块预处理把 A 分块只对对角块做不完全分解块间耦合在迭代中处理。这类预处理在流体力学和电磁场模拟中表现突出实现上可以直接用scipy.linalg.lu_factor对每个块做密分解或用spla.spilu对每个块做稀疏不完全分解。并行环境下要特别注意ILU 的求解过程是串行依赖的大规模并行时可能成为瓶颈。此时考虑三类替代方案多项式预处理用 A 的多项式近似M^{-1}只用 MatVec天然并行、多色 SSOR按图着色重排让同一颜色的点互不依赖可并行、或用域分解法把问题切成子域各算各的。5. 收敛性诊断如何在 Krylov 子空间迭代算法中定位停滞和破裂5.1 看对的量残差、相对残差与 A 范数Krylov 子空间迭代算法的停机准则设置不当会导致两种情况过早停止得到错误解或过晚停止浪费算力。需要区分三个量真实残差||b - Ax_k||、递推残差迭代过程中通过递推隐式维持的量、以及误差||x_k - x_true||。在浮点环境下递推残差可能远小于真实残差这是所有 Krylov 迭代共有的问题。我建议在迭代中同时输出真实残差和相对残差并在接近停机阈值时切换到真实残差验证。对于 CG 还有专门的 A 范数——||x_k - x_true||_A——它直接度量能量误差比二范数更适合同类问题。5.2 停滞的三种模式与各自对策Krylov 收敛曲线出现停滞通常是三种模式之一。第一种是“平台期”残差在一个水平持久不变然后突然下降这常见于非正规矩阵GMRES 理论上不该出现但有限精度下仍可能发生。对策是检查是否因为重启导致丢失了关键方向或者换更高质量的预处理子。第二种是“渐进爬坡”残差持续但极其缓慢下降这是预处理不足的典型信号应该考虑加强预处理而不是增大迭代步数。第三种是“发散振荡”残差上下波动不收敛通常是矩阵非正规、预处理不稳定或迭代格式不适合问题类型。区分这三种模式最快的方法是同时画半对数坐标下的残差曲线和每步的||x_k - x_{k-1}||。如果解的变化量已经小到机器精度但残差还在高位问题出在预处理上如果解变化量很大但残差不降问题可能出在算法的正交性丢失上。5.3 参数调整的最短路径面对不收敛的迭代很多人的第一反应是调大maxiter——但这通常是无效的。正确的调整路径是先看矩阵条件数和特征值分布确认问题是否出在矩阵本身的病态性然后依次尝试从无预处理到 Jacobi、再到 ILU(0)、再到 ILU(k)观察每步收敛步数的变化率如果 ILU 没有带来数量级改善再考虑换算法族比如把 GMRES 换成 BiCGSTAB(l) 或重启 GMRES。一个实用的做法是把预处理质量作为第一优先级。ILU 的参数选择优先于迭代算法的参数选择。迭代参数中rtol的设定要考虑最终用途做特征值问题的内迭代rtol1e-3也许就够做高精度结构分析可能需要rtol1e-12。不要在需求未知时盲目追求高精度这会平白增加几十步迭代。6. 病态矩阵处理与 Krylov 子空间迭代算法的规模化实践6.1 矩阵重排序Cuthill-McKee 对 ILU 的影响预处理子的质量取决于矩阵结构而结构可以通过重排改变。Reverse Cuthill-McKeeRCM是一种带宽缩减算法它把矩阵重新排序成轮廓更窄的形式。对于有限元网格生成的矩阵RCM 重排后 ILU 的填充更集中在对角线附近不完全分解的丢弃误差更小预处理质量显著提升。from scipy.sparse.csgraph import reverse_cuthill_mckee from scipy.sparse import csr_matrix # 假设 A 是 CSR 格式的稀疏矩阵 A_csr csr_matrix(A) perm reverse_cuthill_mckee(A_csr, symmetric_modeTrue) # 重排矩阵和右端项 A_rcm A_csr[perm][:, perm] b_rcm b[perm] # 在重排后的系统上做迭代最后把解映射回原次序 x_rcm, info spla.gmres(A_rcm, b_rcm, rtol1e-10, Milu_precond) x np.empty_like(b) x[perm] x_rcmperm是重排后的下标映射数组A_csr[perm][:, perm]同时对行和列重排保持对称性。RCM 不改变特征值所以 Krylov 收敛的理论性质不变但它改变了预处理子可用的结构信息实际效果往往是迭代步数下降一半以上。值得注意一个容易出错的地方解必须映射回原顺序才能和其他模块对接忘记这一步是接错数据的高频原因。6.2 特征值平移对大规模问题收敛性的改善A的零特征值或接近零的特征值是 Krylov 子空间迭代算法收敛的最大障碍。对平移后的系统(A σI)x b σx做迭代可以改善特征值聚集性但代价是需要额外处理 x 项。更常见的是把它作为预处理手段把平移合并到预处理子中M A σI的近似。这记录了“越接近零的特征值对 Krylov 方法的危害越大”的工程直觉。对 Near-null 空间明显的矩阵比如结构分析中的刚体模态平移量 σ 取一个比最小非零特征值小一两个数量级的数值效果最好。6.3 与多层方法的分工配合当 Krylov 子空间迭代算法在千万自由度级别的结构分析中遭遇收敛瓶颈时可以考虑与多层方法分工协作。底层思路是Krylov 负责处理平滑的高频误差粗网格修正负责处理低频误差。在规模化实践中比较经典的方案是“多层预处理 Krylov 加速”的组合模式——用多重网格或领域分解做预处理子用 GMRES 或 CG 驱动整体收敛。参数设置的经验法则是多层方法的层数取决于网格细化层次和特征值分布Krylov 侧的重启阈值通常设为 30 或 50如果配合代数多重网格预处理CG 通常在 20 步以内收敛。这个组合也是当前主流商用有限元软件在大规模分析时的默认选择。最后提一个容易被忽略的验证细节无论预处理多复杂最终都要额外做一次真实残差检查。因为经过重排、平移和预处理之后递推残差与真实残差的差距可能被放大只有直接计算||b - Ax||才能确认你的 Krylov 子空间迭代算法真的收敛到了该有的精度。本文还有配套的精品资源点击获取