1. Krylov子空间方法到底在解决什么问题先说个最直接的场景。你在做有限元结构分析、流体求解或者图像处理里的某些反问题最后都会落到一个事上解一个大规模的稀疏线性方程组 Axb。规模有多大几十万阶、几百万阶甚至更大。这时候你用高斯消元直接分解先不说内存能不能装下稠密的 L 和 U光是计算量就能让人绝望。而这类矩阵往往是稀疏的大部分元素是零真正有用的信息集中在少数位置上。Krylov子空间方法就是在这种背景下被大量使用的一类迭代解法。它的核心思路其实很好理解不直接改矩阵 A而是从一个初始猜测出发不断构造一个维度逐次递增的子空间叫 Krylov 子空间然后把近似解投影到这个子空间里去找。你可能会问为什么偏偏是 Krylov 子空间因为它是只依赖矩阵乘法和向量内积就能生成的空间。这个特性决定了我们自始至终不需要接触矩阵 A 的显式分解只需要提供计算 A 乘以某个向量的能力这对于大规模稀疏问题来说是决定性的优势。如果你第一次接触这套东西可以把它类比成一个不断用有限信息逼近真实答案的过程。每一轮迭代你可以把当前已有的信息浓缩在一个低维空间里在这个低维空间里求一个局部最优解然后把这个解映射回高维原空间。迭代次数越多低维空间越丰富解就越逼近真实值。Krylov 子空间方法不是某一个具体算法而是一整族算法CG、MINRES、GMRES、BiCGSTAB、QMR 这些在数值线性代数里如雷贯耳的名字都属于这个家族只是针对不同的矩阵性质做了不同的优化。这篇文章面向的是正在学高等数值分析的研究生或者工作中需要自己写迭代求解器、又不想直接调库的工程师。我会从数学原理、方法族谱、收敛性分析、工程实现这些角度展开尽量把背后的为什么讲清楚而不是贴一堆公式就完事。2. Krylov子空间的构造逻辑与投影思想2.1 从 Cayley-Hamilton 定理到 Krylov 子空间先问个问题为什么自然就会想到用 Krylov 子空间来逼近解这背后有个非常漂亮的数学支撑——Cayley-Hamilton 定理。任何一个 n 阶方阵 A 都满足它自己的特征多项式 p(A)0这意味着 A⁻¹ 总可以表示为 A 的 0 到 n-1 次幂的一个多项式组合。于是 Axb 的精确解就可以写成x A⁻¹b q(A)b其中 q 是某个次数不超过 n-1 的多项式。这就告诉我们方程组的解一定落在由 b, Ab, A²b, ..., Aⁿ⁻¹b 张成的空间里这个空间就是 n 维的 Krylov 子空间。这个结论的意义非常深远。它把解一个大型线性系统这个看起来需要矩阵分解的问题转化成了寻找一个合适的多项式的问题。而多项式逼近是可以逐次迭代逼近的不需要一步到位。所以 Krylov 子空间方法本质上是在做一件事用 m 次多项式去近似那个理想的 n-1 次多项式m 远小于 n 时就用一个低维的、代价便宜的近似来换取可接受的误差。2.2 投影法在低维空间里找最优有了 Krylov 子空间之后接下来的问题是怎么在这个空间里选出最好的近似解。这里采用的框架是投影法也就是 Petrov-Galerkin 条件。假设当前 Krylov 子空间是 K_m我想要的近似解 x_m 满足两个条件第一x_m 属于 x₀ 加上 K_mx₀ 是初始猜测通常是零向量第二残差 r_m b - Ax_m 与某个约束子空间 L_m 正交。这里的关键选择在于 L_m 怎么取。取 L_m K_m这叫 Galerkin 投影适合对称正定矩阵CG 方法就是这条路。取 L_m A K_m这叫最小残差投影适合一般非对称矩阵GMRES 走的是这条路。两种路线的本质区别在于前者是想让残差被约束住后者是想让残差的范数直接最小化。你不需要把这两种选择背下来你只需要理解一件事——投影方向的选择直接决定了算法的稳定性和适用性这也解释了这个领域为什么会有那么多不同的方法流派。2.3 为什么要坚持矩阵-向量乘作为基本操作Krylov 子空间方法的实现只依赖两个基本操作矩阵乘向量MatVec和向量内积。这一点在工程上极其重要。你想象一个来自三维有限元分析的刚度矩阵它的非零元素可能只占全部元素的万分之一但每个非零元素都有明确的物理意义。如果用直接法fill-in 效应会把大量零元素变成非零矩阵的存储模式和稀疏结构瞬间被破坏。而用 Krylov 方法A 从头到尾只是一个可以被向量乘的黑盒可以是显式的稀疏矩阵也可以是某种算子比如快速傅里叶变换、多重网格甚至不需要在内存里完整存放。这一点保证了 Krylov 方法可以用于超大规模计算也让它成为多物理场耦合、全波形反演、深度学习里某些大型线性求解场景的首选。很多新接触这套方法的人会忽略这个优势他们觉得迭代法只是不用分解矩阵而已但实际上不用显式接触矩阵元素才是更本质的优势。3. 方法家族全景常见Krylov算法的选型逻辑3.1 CG对称正定矩阵的经典选择共轭梯度法CG是 Krylov 族里最经典、最被广泛使用的算法。它有两个前提矩阵 A 对称正定。在这个前提下Axb 的求解可以等价转化为二次泛函 φ(x) ½xᵀAx - bᵀx 的最小化问题CG 就是在 Krylov 子空间上对这个泛函做方向交替的最小化。CG 的一个显著优点是不需要长递归每步只需要保留上一步的方向向量存储开销极小几个向量而已。这意味着它能很好处理非常大的问题。但它的收敛性对矩阵条件数非常敏感收敛速度大致与 κ(A)⁻¹/² 相关的误差衰减因子有关条件数太大时收敛会很慢。我在实际项目中遇到过不少把 CG 用在不合适场合的情况。最常见的错误是矩阵根本不是对称正定的或者没有正确做预处理就硬跑。也有人说 CG 怎么发散——严格来说 CG 在数学上对对称正定矩阵是保证收敛的但如果浮点舍入误差累积导致方向向量失去共轭性收敛就会停滞甚至出现残差反弹这时候往往需要做重启或者更换预处理子。3.2 GMRES非对称问题的通用兜底方案如果你的矩阵不对称CG 的正交性基础瞬间崩塌。最直接的替代思路是用 Arnoldi 过程把矩阵投影到 Hessenberg 矩阵上然后用最小二乘在每一步求一个极小残差解这就是 GMRES广义最小残差法。GMRES 最大的优点是理论上保证残差每一步都不增是最稳妥的 Krylov 方法之一。但它的代价是存储和计算量随时间增长因为 Arnoldi 过程生成的正交基需要全部保存而且每步重正交化的代价也在增加。这就是为什么实际使用时几乎都会配一个重启策略设定一个重启步数 k迭代到 k 步之后把当前解当作新的初始猜测重新开始。重启蒙在工程上是必需品但数学上有个尴尬重启后的 GMRES 不保证收敛因为丢失了之前所有基的信息。实际测试中我见过重启太频繁导致收敛远慢于不重启版本的情况所以重启步长的选择需要针对问题调参。3.3 BiCGSTAB 和 QMR存储受限时的折中方案GMRES 的存储瓶颈在某些超大规模场景下是不可接受的。双共轭梯度法BiCG和它的稳定版本 BiCGSTAB 走的是另一条路线构造两个双正交的 Krylov 子空间用短递归逼近解存储开销和 CG 一样低同时适用于非对称矩阵。BiCGSTAB 是工程上使用频率最高的非对称 Krylov 方法因为它在存储和收敛性之间取得了一个很好的平衡点特别是对来自偏微分方程离散化的矩阵效果通常不错。它的缺点是容易出现 breakdown比如某个内积恰好为零而且残差曲线经常会有不规则震荡不像 GMRES 那样单调下降。实际调参时要注意BiCGSTAB 的收敛诊断不能光看残差值有时候残差骤降又骤升需要配合迭代次数上限来保护。表里我总结一下几个常用方法的选择逻辑方法适用矩阵存储开销收敛行为主要风险CG对称正定低固定平滑下降条件数过大时收敛慢MINRES对称不定低固定平滑下降对零特征值敏感GMRES一般非对称随迭代增长残差不增重启导致停滞BiCGSTAB一般非对称低固定可能震荡breakdown 风险QMR一般非对称低固定较平滑实现复杂工程用少3.4 怎么选几个实践判断标准面对实际问题时我不会一开始就急着选算法而是先做两件事第一确认矩阵的对称性和正定性第二估算一下能承受的存储量。之后按这样的优先级来做选择如果矩阵对称正定无脑 CG配合一个好的预处理子。如果对称但不定用 MINRES如果你没听过这个记住它是 CG 在对称不定情况的推广。如果矩阵非对称但矩阵规模很大、存储受限先试 BiCGSTAB跑不通再考虑其他。如果矩阵非对称且你有足够的存储GMRES 是最安全的选择。实际工程里我见过很多人图省事直接用 GMRES 解决一切问题这在小规模测试里可行但在超大规模场景下内存会爆掉所以从一开始就养成分场景选算法的习惯很重要。4. 收敛性分析与预处理真正拉开差距的地方4.1 特征值分布如何决定收敛速度Krylov 方法的收敛性不是玄学它跟矩阵的特征值分布有非常明确的数学联系。对 CG 来说收敛速度上界由条件数 κ 给出经过 m 步迭代后误差能量的衰减因子大约在 (√κ-1)/(√κ1) 的 m 次方左右。条件数为 1000 时粗略估算一下就知道收敛会非常吃力。但如果特征值聚集成少量几个簇收敛会快很多这就是预处理发挥作用的原因。GMRES 的收敛性分析更复杂它不单看条件数还看特征值的分布形状、矩阵是否可对角化等。实际经验中如果非对称矩阵的特征值分布靠近原点或者有接近零的特征值GMRES 的收敛曲线往往会出现一段平台期就是前若干步残差几乎不动然后突然开始下降。遇到这种情况不要急着怀疑程序写错了这其实是 Krylov 方法对特征值响应的正常现象。预处理本质上就是找一个矩阵 M把原方程改写成更容易让 Krylov 方法收敛的形式。比如用 M⁻¹ 左乘方程两边如果 M 近似于 A那么新矩阵 M⁻¹A 的特征值就聚集在 1 附近收敛自然会变快。但注意不能真的去求 M⁻¹只能是求解形如 Mzr 的辅助方程所以 M 的选择必须兼顾近似程度和求解 Mzr 的代价两个方面。4.2 预处理子的层次结构与工程选型预处理子的选择是 Krylov 方法实用化中最关键的一环没有之一。粗略可以分为三个层次第一层是不完全 LU 分解ILU。ILU 的做法是像做 LU 分解一样走流程但只保留原矩阵稀疏结构上的元素舍弃 fill-in 元素。它本质上是在假装分解得到一个近似的分解产物用起来效果通常不错。ILU 有个关键参数叫 fill-in level调大一点填元更多、分解更准确但内存和分解时间也上去了。实际工程里 ILU(0)不允许任何 fill-in是最常见的起点。第二层是稀疏近似逆预处理。它会直接计算一个稀疏矩阵 M⁻¹显式逼近 A⁻¹优势在于 M⁻¹ 可以直接做矩阵乘不需要回代求解。这类预处理子在某些并行场景下更友好因为求解 Mzr 这个步骤天然适合并行化。第三层是领域分解和多重网格类预处理子这类方法更接近求解器效果通常最强。比如流体计算里常用的代数多重网格AMG对椭圆型偏微分方程离散出来的矩阵效果极好但实现复杂不适合直接用数值分析课程里的小例子来评估。从我的经验看选择预处理子的逻辑是先用 ILU(0) BiCGSTAB 做一个基线观察迭代次数和单步代价然后根据矩阵来源决定是加强预处理比如 ILU(1)、ILUT还是换更高级的方法。当矩阵结构非常规整比如来自规则网格上的有限差分时AMG 往往能带来数量级级别的收敛加速。4.3 一个具体的收敛性调试案例有次我在处理一个来自流固耦合界面的非对称系统矩阵大约是 20 万阶稀疏度尚可但特征值分布很糟糕最大的特征值模到了 10⁸ 量级而最小的非零特征值只有 10⁻² 左右导致 GMRES(50) 重启了二十多次残差还在 10⁻³ 徘徊。当时我没有急着改算法而是把矩阵做了一次行均衡行范数归一化然后再配合 ILU(1) 预处理效果立刻改善同样的重启设置十几次迭代就降到了 10⁻⁸。这个案例告诉我们一个朴素却常被忽略的道理预处理之前的清洗步骤可能比预处理本身更关键。矩阵的尺度不一致会导致特征值分布变得极差有时候单纯的对角缩放diagonal scaling就能让收敛速度翻倍。很多教科书不会细讲这些 scale 的细节但在实际代码里它们往往决定算法的成败。5. 实操过程与关键环节从公式到能跑的代码5.1 Arnoldi 过程Krylov 方法的引擎不管哪种 Krylov 方法最底层的引擎都是 Arnoldi 过程它的作用是把矩阵 A 在 Krylov 子空间上的作用投影成一个上 Hessenberg 矩阵 H。你可以把它理解为是对矩阵 A 在特定子空间上的一个压缩表示。每一步Arnoldi 过程做的事情就是取当前基向量乘以 A然后把新得到的向量关于已经存在的所有基向量做正交化得到新的基向量。这个过程和 Gram-Schmidt 正交化很像但有一个需要留意的实际坑当迭代步数较多时简单 Gram-Schmidt 会因为舍入误差导致正交性严重丧失所以工程实现里几乎都要用 modified Gram-SchmidtMGS或者必要时做两次正交化DGKS 修正。我把这个问题的严重性再强调一遍正交性一旦丧失Arnoldi 关系 A V_m V_m H_m h_{m1,m} v_{m1}e_mᵀ 就不成立了后面所有基于它的算法推导都会失真。下面给一个用 Python/NumPy 实现的简化版 Arnoldi 过程它可以直接被 CG 以外的各种方法复用import numpy as np def arnoldi(A, v0, m): 简化版 Arnoldi 过程。 输入: 矩阵 A可以是线性算子初始向量 v0迭代步数 m 输出: 正交基矩阵 V (n x m1)Hessenberg 矩阵 H (m1 x m) n v0.shape[0] V np.zeros((n, m1), dtypenp.float64) H np.zeros((m1, m), dtypenp.float64) V[:, 0] v0 / np.linalg.norm(v0) for j in range(m): w A V[:, j] for i in range(j1): H[i, j] np.dot(V[:, i], w) w w - H[i, j] * V[:, i] H[j1, j] np.linalg.norm(w) if H[j1, j] 1e-14: # 此时 w 为零向量Krylov 子空间已经达到不变子空间 V[:, j1] 0.0 else: V[:, j1] w / H[j1, j] return V[:, :m], H # 简单的随机对称正定矩阵测试 n 100 A np.random.randn(n, n) A A A.T n * np.eye(n) v0 np.random.randn(n) V, H arnoldi(A, v0, 10) # 检查 Arnoldi 关系是否成立 residual np.linalg.norm(A V[:, :10] - V[:, :10] H[:10, :]) print(Arnoldi 关系残差:, residual)注意这段代码里对 H[j1,j] 接近零的处理。这个幸运 breakdown在理论上说明已经找到了一个不变子空间迭代可以终止实际代码里如果忽略它后面的除法就会出现除零警告或者产生无穷大。5.2 GMRES 的核心迭代循环有了 Arnoldi 过程GMRES 的框架就简单多了。每一步迭代实际上是在解一个最小二乘问题在 K_m 里找 x_m 使得 ||b - Ax_m||₂ 最小。利用 Arnoldi 关系这个最小二乘问题可以等价转化为一个 (m1)×m 的小规模最小二乘问题因为 H 的上 Hessenberg 结构可以用 Givens 旋转逐步做 QR 分解来高效求解。我自己在实际写 GMRES 时一开始犯的错误是直接用满秩最小二乘来解这个小问题结果每步都要做一次 O(m³) 的分解迭代步数一多就慢得离谱。换用 Givens 旋转后每次迭代只需要 O(m) 的旋转操作整体开销小得多这才是实际代码中应该采用的方式。以下是 GMRES 核心循环的关键伪代码def gmres(A, b, x0, restart, max_iter, tol): x x0.copy() r b - A x beta np.linalg.norm(r) for outer in range(max_iter): V, H arnoldi(A, r, restart) # 把最小二乘问题转化为标准形式 e1 np.zeros(restart 1) e1[0] beta # 用 Givens 旋转求解 min || beta e1 - H y ||_2 y, givens_resid solve_least_squares_givens(H, e1) # 更新解 x x V y # 计算新的残差 r b - A x if np.linalg.norm(r) tol: break return x这里的 solve_least_squares_givens 就是利用 H 的 Hessenberg 结构做连续 Givens 旋转把 H 变成上三角矩阵同时更新右端项。这个技巧本身不复杂但如果没有实现过你会觉得 GMRES 的代码怎么这么绕。5.3 停止准则与残差诊断的正确姿势实际工程中一个特别容易踩坑的地方是停止准则。很多人直接检查残差范数是否小于 tol但没注意残差范数小不代表真实误差小。当矩阵病态时一个很小的残差对应的解误差可能非常大。还有两个细节残差的定义要和算法内部一致比如很多实现检查的是预处理之后的残差如果你用预处理后的残差做停判但输出的是未预处理的残差你就会发现明明说已经收敛残差却很大的诡异现象。我在项目中习惯同时输出三类诊断信息每一次迭代的残差范数、相对残差、还有解的变化量。对于 CG 类方法残差曲线出现长期平台往往是预处理不够强对于 BiCGSTAB残差突然上升后又快速下降是正常现象不要一看到上升就中断。关于相对残差和绝对残差的选择我建议统一用相对残差||b - Ax_k||₂ / ||b||₂。这样做的好处是可以比较不同规模、不同右端项的问题之间的收敛行为。另外建议加上一个最大迭代次数作为硬保护避免预处理失效时程序无限循环。很多软件包默认还加了一条停滞检测如果连续若干步残差几乎不变干脆给出警告退出这在调试期能节省大量时间。5.4 MATLAB 和 Python 中的快速验证方式如果你只是想快速验证一个想法没必要从头写 Krylov 实现直接用现成的接口配合适当的监控手段就够了。MATLAB 里gmres 和 pcg 的函数签名都支持回调函数你可以传入一个 function handle 来记录每步残差也可以修改重启步长、容差和最大迭代次数。首次拿到一个新的稀疏矩阵时我的习惯是先用默认设置跑一次 GMRES确认矩阵本身没有符号错误再微调预处理和参数。Python 里最常用的工具是 SciPy 的 scipy.sparse.linalg里面提供了 gmres、bicgstab、cg 等迭代器接口返回的 info 字段很重要info0 表示成功收敛不为 0 则需要查具体含义。另外你可以把算子封装成 LinearOperator 对象传给这些函数矩阵不一定要显式给出这在某些只需要 MatVec 的应用场景里非常实用。6. 常见问题与排查技巧我踩过的那些坑6.1 残差长期不下降第一反应查什么残差长期不下降是最常遇到的问题我的排查顺序基本是固定的。先看有没有做预处理。如果完全没做预处理而矩阵条件数又很大不下降是正常现象先补一个对角缩放或者 ILU(0)。再看矩阵是否真的对称正定。有个项目里我碰到过一个矩阵从物理模型推导应该是对称的但因为网格编号和边界条件的处理数值上对称性被破坏了一个极小的量导致 CG 行为怪异后来用对称化处理解决。接着检查右端项是否在 A 的值域之外这个问题不一定代表代码有 bug可能只是问题本身无解或者有无穷多解需要检查相容性条件。最后如果用的是 GMRES看看是不是重启步长太小导致 K 空间无法捕捉到关键特征向量信息适当调大 restart比如从 20 调整到 50、100通常有帮助。6.2 内存和时间的平衡大规模场景下的隐形成本Krylov 方法虽然避免了矩阵分解的内存爆炸但迭代本身也有自己的内存成本。GMRES 的基向量随迭代步数线性增长restart 设为 500 时如果矩阵阶数百万只需要几十个基向量内存消耗就明显起来了。这里有个被我反复验证的经验处理超大规模问题优先考虑 BiCGSTAB 或带预处理的 CG而不是 GMRES除非你有充足的内存并且确实需要 GMRES 的稳定性。时间成本方面除了 MatVec 本身最大的开销往往来自向量内积和全局通信。在分布式的 MPI 环境里每一轮 Krylov 迭代通常需要至少一到两次全局归约all-reduce来计算内积这会导致一个分摊开销随进程数增加而恶化。很多大规模的并行迭代求解器性能瓶颈就在这里。如果追求极致扩展性可以考虑 s-step Krylov 方法或者通信避免算法不过这些方法数值稳定性还需要仔细评估目前还不是主流工程选项。6.3 矩阵特征值的快速估计不需要求特征值也能判断收敛很多人会问我想知道这个矩阵适不适合 Krylov 方法但求特征值本身太贵了怎么办其实有更廉价的替代方案用 Lanczos/Arnoldi 过程本身的中间结果做估计。GMRES 的 Hessenberg 矩阵 H_m 的特征值会随着迭代逐渐逼近 A 的部分特征值称之为 Ritz 值。你可以直接在求解过程中把它们打印出来如果看到某些 Ritz 值非常接近 0说明矩阵有接近奇异的模态这几乎一定会拖慢收敛。这个方法不需要额外代价因为 H_m 本来就计算好了。另外矩阵的迹和 Frobenius 范数可以给出特征值分布范围的粗略估计这些只需要 O(nnz) 的扫描即可获得。再配合幂迭代算一个最大特征值的近似你就能在开始正式迭代之前大致判断这个问题的难度级别。6.4 复数矩阵与多右端项问题Krylov 方法在复数域同样适用但要注意两点内积定义需要取共轭另外 GMRES 里的最小二乘问题要在复数意义下理解。实际代码中如果矩阵是 Hermitian 正定的CG 依然可以用只是内积和系数要对应修改。多右端项问题在电磁散射、地震成像等场景里非常常见多个右端项共享同一个矩阵。最朴素的做法是逐个用 Krylov 方法求解但这会浪费大量重复的预处理代价。可以尝试块 Krylov 方法block GMRES、block CG一次迭代同时处理多个右端项矩阵-向量乘也变成矩阵-块乘通常块宽度不大时比如 8-16效率提升明显。不过块方法对块间线性相关性敏感右端项如果有近似线性相关的情况需要做块降秩处理否则数值稳定性会出问题。6.5 舍入误差、除以零与幸运 breakdown浮点环境下的 Krylov 方法最让人头疼的问题之一是 breakdown。它表现为算法某一步突然出现除以零或者接近除以零的情况。在 Arnoldi 过程中这对应 h_{j1,j} 接近零在 BiCGSTAB 中这是 r₀ᵀv₀ 之类的内积接近零。处理原则很简单检测到接近零的情况时不要直接除以它要么提前终止如果这正好说明已经收敛要么做重启。另一种常见问题是 CG 的残差曲线出现微小的上升然后继续下降。这多半不是迭代发散而是舍入误差破坏了共轭性。解决办法是定期做一次重置操作重新计算真正的残差 b-Ax然后用它重新开始方向向量。很多生产级代码里会每隔 50 或 100 步做一次这样操作实际效果非常明显。7. 最后说点实在的我在很长一段时间里犯过一个认知错误觉得 Krylov 方法只要把迭代公式背熟、把库调对就够了。后来被真实问题教育过几次才明白真正的难点从来不在算法本身而在于根据具体的矩阵去选择合适的算法和预处理组合。接触一个新的稀疏矩阵时一定要先花时间弄清楚它的来源、结构特征、条件数和特征值分布再做选型决策这十几分钟的投入通常能让你后面少调几天的参。另外我建议每个学数值分析的人至少手写一遍 GMRES 和 CG。手写不是为了在生产里用而是因为只有亲手实现过 Arnoldi 过程、理解过 Givens 旋转为什么这样构造你才能对文档里那些晦涩的参数有体感。遇到问题的时候你的直觉才不会一片空白。这个领域后续还有不少可以深挖的方向比如针对大规模并行环境的通信避免算法、自适应预处理、机器学习辅助的预处理子选择都是目前很活跃的研究方向。但不管技术怎么演进理解清楚 Krylov 子空间的基本思想和工程约束始终是所有高阶技巧的基础。 Krylov子空间方法全解析：从数学原理到工程实践
1. Krylov子空间方法到底在解决什么问题先说个最直接的场景。你在做有限元结构分析、流体求解或者图像处理里的某些反问题最后都会落到一个事上解一个大规模的稀疏线性方程组 Axb。规模有多大几十万阶、几百万阶甚至更大。这时候你用高斯消元直接分解先不说内存能不能装下稠密的 L 和 U光是计算量就能让人绝望。而这类矩阵往往是稀疏的大部分元素是零真正有用的信息集中在少数位置上。Krylov子空间方法就是在这种背景下被大量使用的一类迭代解法。它的核心思路其实很好理解不直接改矩阵 A而是从一个初始猜测出发不断构造一个维度逐次递增的子空间叫 Krylov 子空间然后把近似解投影到这个子空间里去找。你可能会问为什么偏偏是 Krylov 子空间因为它是只依赖矩阵乘法和向量内积就能生成的空间。这个特性决定了我们自始至终不需要接触矩阵 A 的显式分解只需要提供计算 A 乘以某个向量的能力这对于大规模稀疏问题来说是决定性的优势。如果你第一次接触这套东西可以把它类比成一个不断用有限信息逼近真实答案的过程。每一轮迭代你可以把当前已有的信息浓缩在一个低维空间里在这个低维空间里求一个局部最优解然后把这个解映射回高维原空间。迭代次数越多低维空间越丰富解就越逼近真实值。Krylov 子空间方法不是某一个具体算法而是一整族算法CG、MINRES、GMRES、BiCGSTAB、QMR 这些在数值线性代数里如雷贯耳的名字都属于这个家族只是针对不同的矩阵性质做了不同的优化。这篇文章面向的是正在学高等数值分析的研究生或者工作中需要自己写迭代求解器、又不想直接调库的工程师。我会从数学原理、方法族谱、收敛性分析、工程实现这些角度展开尽量把背后的为什么讲清楚而不是贴一堆公式就完事。2. Krylov子空间的构造逻辑与投影思想2.1 从 Cayley-Hamilton 定理到 Krylov 子空间先问个问题为什么自然就会想到用 Krylov 子空间来逼近解这背后有个非常漂亮的数学支撑——Cayley-Hamilton 定理。任何一个 n 阶方阵 A 都满足它自己的特征多项式 p(A)0这意味着 A⁻¹ 总可以表示为 A 的 0 到 n-1 次幂的一个多项式组合。于是 Axb 的精确解就可以写成x A⁻¹b q(A)b其中 q 是某个次数不超过 n-1 的多项式。这就告诉我们方程组的解一定落在由 b, Ab, A²b, ..., Aⁿ⁻¹b 张成的空间里这个空间就是 n 维的 Krylov 子空间。这个结论的意义非常深远。它把解一个大型线性系统这个看起来需要矩阵分解的问题转化成了寻找一个合适的多项式的问题。而多项式逼近是可以逐次迭代逼近的不需要一步到位。所以 Krylov 子空间方法本质上是在做一件事用 m 次多项式去近似那个理想的 n-1 次多项式m 远小于 n 时就用一个低维的、代价便宜的近似来换取可接受的误差。2.2 投影法在低维空间里找最优有了 Krylov 子空间之后接下来的问题是怎么在这个空间里选出最好的近似解。这里采用的框架是投影法也就是 Petrov-Galerkin 条件。假设当前 Krylov 子空间是 K_m我想要的近似解 x_m 满足两个条件第一x_m 属于 x₀ 加上 K_mx₀ 是初始猜测通常是零向量第二残差 r_m b - Ax_m 与某个约束子空间 L_m 正交。这里的关键选择在于 L_m 怎么取。取 L_m K_m这叫 Galerkin 投影适合对称正定矩阵CG 方法就是这条路。取 L_m A K_m这叫最小残差投影适合一般非对称矩阵GMRES 走的是这条路。两种路线的本质区别在于前者是想让残差被约束住后者是想让残差的范数直接最小化。你不需要把这两种选择背下来你只需要理解一件事——投影方向的选择直接决定了算法的稳定性和适用性这也解释了这个领域为什么会有那么多不同的方法流派。2.3 为什么要坚持矩阵-向量乘作为基本操作Krylov 子空间方法的实现只依赖两个基本操作矩阵乘向量MatVec和向量内积。这一点在工程上极其重要。你想象一个来自三维有限元分析的刚度矩阵它的非零元素可能只占全部元素的万分之一但每个非零元素都有明确的物理意义。如果用直接法fill-in 效应会把大量零元素变成非零矩阵的存储模式和稀疏结构瞬间被破坏。而用 Krylov 方法A 从头到尾只是一个可以被向量乘的黑盒可以是显式的稀疏矩阵也可以是某种算子比如快速傅里叶变换、多重网格甚至不需要在内存里完整存放。这一点保证了 Krylov 方法可以用于超大规模计算也让它成为多物理场耦合、全波形反演、深度学习里某些大型线性求解场景的首选。很多新接触这套方法的人会忽略这个优势他们觉得迭代法只是不用分解矩阵而已但实际上不用显式接触矩阵元素才是更本质的优势。3. 方法家族全景常见Krylov算法的选型逻辑3.1 CG对称正定矩阵的经典选择共轭梯度法CG是 Krylov 族里最经典、最被广泛使用的算法。它有两个前提矩阵 A 对称正定。在这个前提下Axb 的求解可以等价转化为二次泛函 φ(x) ½xᵀAx - bᵀx 的最小化问题CG 就是在 Krylov 子空间上对这个泛函做方向交替的最小化。CG 的一个显著优点是不需要长递归每步只需要保留上一步的方向向量存储开销极小几个向量而已。这意味着它能很好处理非常大的问题。但它的收敛性对矩阵条件数非常敏感收敛速度大致与 κ(A)⁻¹/² 相关的误差衰减因子有关条件数太大时收敛会很慢。我在实际项目中遇到过不少把 CG 用在不合适场合的情况。最常见的错误是矩阵根本不是对称正定的或者没有正确做预处理就硬跑。也有人说 CG 怎么发散——严格来说 CG 在数学上对对称正定矩阵是保证收敛的但如果浮点舍入误差累积导致方向向量失去共轭性收敛就会停滞甚至出现残差反弹这时候往往需要做重启或者更换预处理子。3.2 GMRES非对称问题的通用兜底方案如果你的矩阵不对称CG 的正交性基础瞬间崩塌。最直接的替代思路是用 Arnoldi 过程把矩阵投影到 Hessenberg 矩阵上然后用最小二乘在每一步求一个极小残差解这就是 GMRES广义最小残差法。GMRES 最大的优点是理论上保证残差每一步都不增是最稳妥的 Krylov 方法之一。但它的代价是存储和计算量随时间增长因为 Arnoldi 过程生成的正交基需要全部保存而且每步重正交化的代价也在增加。这就是为什么实际使用时几乎都会配一个重启策略设定一个重启步数 k迭代到 k 步之后把当前解当作新的初始猜测重新开始。重启蒙在工程上是必需品但数学上有个尴尬重启后的 GMRES 不保证收敛因为丢失了之前所有基的信息。实际测试中我见过重启太频繁导致收敛远慢于不重启版本的情况所以重启步长的选择需要针对问题调参。3.3 BiCGSTAB 和 QMR存储受限时的折中方案GMRES 的存储瓶颈在某些超大规模场景下是不可接受的。双共轭梯度法BiCG和它的稳定版本 BiCGSTAB 走的是另一条路线构造两个双正交的 Krylov 子空间用短递归逼近解存储开销和 CG 一样低同时适用于非对称矩阵。BiCGSTAB 是工程上使用频率最高的非对称 Krylov 方法因为它在存储和收敛性之间取得了一个很好的平衡点特别是对来自偏微分方程离散化的矩阵效果通常不错。它的缺点是容易出现 breakdown比如某个内积恰好为零而且残差曲线经常会有不规则震荡不像 GMRES 那样单调下降。实际调参时要注意BiCGSTAB 的收敛诊断不能光看残差值有时候残差骤降又骤升需要配合迭代次数上限来保护。表里我总结一下几个常用方法的选择逻辑方法适用矩阵存储开销收敛行为主要风险CG对称正定低固定平滑下降条件数过大时收敛慢MINRES对称不定低固定平滑下降对零特征值敏感GMRES一般非对称随迭代增长残差不增重启导致停滞BiCGSTAB一般非对称低固定可能震荡breakdown 风险QMR一般非对称低固定较平滑实现复杂工程用少3.4 怎么选几个实践判断标准面对实际问题时我不会一开始就急着选算法而是先做两件事第一确认矩阵的对称性和正定性第二估算一下能承受的存储量。之后按这样的优先级来做选择如果矩阵对称正定无脑 CG配合一个好的预处理子。如果对称但不定用 MINRES如果你没听过这个记住它是 CG 在对称不定情况的推广。如果矩阵非对称但矩阵规模很大、存储受限先试 BiCGSTAB跑不通再考虑其他。如果矩阵非对称且你有足够的存储GMRES 是最安全的选择。实际工程里我见过很多人图省事直接用 GMRES 解决一切问题这在小规模测试里可行但在超大规模场景下内存会爆掉所以从一开始就养成分场景选算法的习惯很重要。4. 收敛性分析与预处理真正拉开差距的地方4.1 特征值分布如何决定收敛速度Krylov 方法的收敛性不是玄学它跟矩阵的特征值分布有非常明确的数学联系。对 CG 来说收敛速度上界由条件数 κ 给出经过 m 步迭代后误差能量的衰减因子大约在 (√κ-1)/(√κ1) 的 m 次方左右。条件数为 1000 时粗略估算一下就知道收敛会非常吃力。但如果特征值聚集成少量几个簇收敛会快很多这就是预处理发挥作用的原因。GMRES 的收敛性分析更复杂它不单看条件数还看特征值的分布形状、矩阵是否可对角化等。实际经验中如果非对称矩阵的特征值分布靠近原点或者有接近零的特征值GMRES 的收敛曲线往往会出现一段平台期就是前若干步残差几乎不动然后突然开始下降。遇到这种情况不要急着怀疑程序写错了这其实是 Krylov 方法对特征值响应的正常现象。预处理本质上就是找一个矩阵 M把原方程改写成更容易让 Krylov 方法收敛的形式。比如用 M⁻¹ 左乘方程两边如果 M 近似于 A那么新矩阵 M⁻¹A 的特征值就聚集在 1 附近收敛自然会变快。但注意不能真的去求 M⁻¹只能是求解形如 Mzr 的辅助方程所以 M 的选择必须兼顾近似程度和求解 Mzr 的代价两个方面。4.2 预处理子的层次结构与工程选型预处理子的选择是 Krylov 方法实用化中最关键的一环没有之一。粗略可以分为三个层次第一层是不完全 LU 分解ILU。ILU 的做法是像做 LU 分解一样走流程但只保留原矩阵稀疏结构上的元素舍弃 fill-in 元素。它本质上是在假装分解得到一个近似的分解产物用起来效果通常不错。ILU 有个关键参数叫 fill-in level调大一点填元更多、分解更准确但内存和分解时间也上去了。实际工程里 ILU(0)不允许任何 fill-in是最常见的起点。第二层是稀疏近似逆预处理。它会直接计算一个稀疏矩阵 M⁻¹显式逼近 A⁻¹优势在于 M⁻¹ 可以直接做矩阵乘不需要回代求解。这类预处理子在某些并行场景下更友好因为求解 Mzr 这个步骤天然适合并行化。第三层是领域分解和多重网格类预处理子这类方法更接近求解器效果通常最强。比如流体计算里常用的代数多重网格AMG对椭圆型偏微分方程离散出来的矩阵效果极好但实现复杂不适合直接用数值分析课程里的小例子来评估。从我的经验看选择预处理子的逻辑是先用 ILU(0) BiCGSTAB 做一个基线观察迭代次数和单步代价然后根据矩阵来源决定是加强预处理比如 ILU(1)、ILUT还是换更高级的方法。当矩阵结构非常规整比如来自规则网格上的有限差分时AMG 往往能带来数量级级别的收敛加速。4.3 一个具体的收敛性调试案例有次我在处理一个来自流固耦合界面的非对称系统矩阵大约是 20 万阶稀疏度尚可但特征值分布很糟糕最大的特征值模到了 10⁸ 量级而最小的非零特征值只有 10⁻² 左右导致 GMRES(50) 重启了二十多次残差还在 10⁻³ 徘徊。当时我没有急着改算法而是把矩阵做了一次行均衡行范数归一化然后再配合 ILU(1) 预处理效果立刻改善同样的重启设置十几次迭代就降到了 10⁻⁸。这个案例告诉我们一个朴素却常被忽略的道理预处理之前的清洗步骤可能比预处理本身更关键。矩阵的尺度不一致会导致特征值分布变得极差有时候单纯的对角缩放diagonal scaling就能让收敛速度翻倍。很多教科书不会细讲这些 scale 的细节但在实际代码里它们往往决定算法的成败。5. 实操过程与关键环节从公式到能跑的代码5.1 Arnoldi 过程Krylov 方法的引擎不管哪种 Krylov 方法最底层的引擎都是 Arnoldi 过程它的作用是把矩阵 A 在 Krylov 子空间上的作用投影成一个上 Hessenberg 矩阵 H。你可以把它理解为是对矩阵 A 在特定子空间上的一个压缩表示。每一步Arnoldi 过程做的事情就是取当前基向量乘以 A然后把新得到的向量关于已经存在的所有基向量做正交化得到新的基向量。这个过程和 Gram-Schmidt 正交化很像但有一个需要留意的实际坑当迭代步数较多时简单 Gram-Schmidt 会因为舍入误差导致正交性严重丧失所以工程实现里几乎都要用 modified Gram-SchmidtMGS或者必要时做两次正交化DGKS 修正。我把这个问题的严重性再强调一遍正交性一旦丧失Arnoldi 关系 A V_m V_m H_m h_{m1,m} v_{m1}e_mᵀ 就不成立了后面所有基于它的算法推导都会失真。下面给一个用 Python/NumPy 实现的简化版 Arnoldi 过程它可以直接被 CG 以外的各种方法复用import numpy as np def arnoldi(A, v0, m): 简化版 Arnoldi 过程。 输入: 矩阵 A可以是线性算子初始向量 v0迭代步数 m 输出: 正交基矩阵 V (n x m1)Hessenberg 矩阵 H (m1 x m) n v0.shape[0] V np.zeros((n, m1), dtypenp.float64) H np.zeros((m1, m), dtypenp.float64) V[:, 0] v0 / np.linalg.norm(v0) for j in range(m): w A V[:, j] for i in range(j1): H[i, j] np.dot(V[:, i], w) w w - H[i, j] * V[:, i] H[j1, j] np.linalg.norm(w) if H[j1, j] 1e-14: # 此时 w 为零向量Krylov 子空间已经达到不变子空间 V[:, j1] 0.0 else: V[:, j1] w / H[j1, j] return V[:, :m], H # 简单的随机对称正定矩阵测试 n 100 A np.random.randn(n, n) A A A.T n * np.eye(n) v0 np.random.randn(n) V, H arnoldi(A, v0, 10) # 检查 Arnoldi 关系是否成立 residual np.linalg.norm(A V[:, :10] - V[:, :10] H[:10, :]) print(Arnoldi 关系残差:, residual)注意这段代码里对 H[j1,j] 接近零的处理。这个幸运 breakdown在理论上说明已经找到了一个不变子空间迭代可以终止实际代码里如果忽略它后面的除法就会出现除零警告或者产生无穷大。5.2 GMRES 的核心迭代循环有了 Arnoldi 过程GMRES 的框架就简单多了。每一步迭代实际上是在解一个最小二乘问题在 K_m 里找 x_m 使得 ||b - Ax_m||₂ 最小。利用 Arnoldi 关系这个最小二乘问题可以等价转化为一个 (m1)×m 的小规模最小二乘问题因为 H 的上 Hessenberg 结构可以用 Givens 旋转逐步做 QR 分解来高效求解。我自己在实际写 GMRES 时一开始犯的错误是直接用满秩最小二乘来解这个小问题结果每步都要做一次 O(m³) 的分解迭代步数一多就慢得离谱。换用 Givens 旋转后每次迭代只需要 O(m) 的旋转操作整体开销小得多这才是实际代码中应该采用的方式。以下是 GMRES 核心循环的关键伪代码def gmres(A, b, x0, restart, max_iter, tol): x x0.copy() r b - A x beta np.linalg.norm(r) for outer in range(max_iter): V, H arnoldi(A, r, restart) # 把最小二乘问题转化为标准形式 e1 np.zeros(restart 1) e1[0] beta # 用 Givens 旋转求解 min || beta e1 - H y ||_2 y, givens_resid solve_least_squares_givens(H, e1) # 更新解 x x V y # 计算新的残差 r b - A x if np.linalg.norm(r) tol: break return x这里的 solve_least_squares_givens 就是利用 H 的 Hessenberg 结构做连续 Givens 旋转把 H 变成上三角矩阵同时更新右端项。这个技巧本身不复杂但如果没有实现过你会觉得 GMRES 的代码怎么这么绕。5.3 停止准则与残差诊断的正确姿势实际工程中一个特别容易踩坑的地方是停止准则。很多人直接检查残差范数是否小于 tol但没注意残差范数小不代表真实误差小。当矩阵病态时一个很小的残差对应的解误差可能非常大。还有两个细节残差的定义要和算法内部一致比如很多实现检查的是预处理之后的残差如果你用预处理后的残差做停判但输出的是未预处理的残差你就会发现明明说已经收敛残差却很大的诡异现象。我在项目中习惯同时输出三类诊断信息每一次迭代的残差范数、相对残差、还有解的变化量。对于 CG 类方法残差曲线出现长期平台往往是预处理不够强对于 BiCGSTAB残差突然上升后又快速下降是正常现象不要一看到上升就中断。关于相对残差和绝对残差的选择我建议统一用相对残差||b - Ax_k||₂ / ||b||₂。这样做的好处是可以比较不同规模、不同右端项的问题之间的收敛行为。另外建议加上一个最大迭代次数作为硬保护避免预处理失效时程序无限循环。很多软件包默认还加了一条停滞检测如果连续若干步残差几乎不变干脆给出警告退出这在调试期能节省大量时间。5.4 MATLAB 和 Python 中的快速验证方式如果你只是想快速验证一个想法没必要从头写 Krylov 实现直接用现成的接口配合适当的监控手段就够了。MATLAB 里gmres 和 pcg 的函数签名都支持回调函数你可以传入一个 function handle 来记录每步残差也可以修改重启步长、容差和最大迭代次数。首次拿到一个新的稀疏矩阵时我的习惯是先用默认设置跑一次 GMRES确认矩阵本身没有符号错误再微调预处理和参数。Python 里最常用的工具是 SciPy 的 scipy.sparse.linalg里面提供了 gmres、bicgstab、cg 等迭代器接口返回的 info 字段很重要info0 表示成功收敛不为 0 则需要查具体含义。另外你可以把算子封装成 LinearOperator 对象传给这些函数矩阵不一定要显式给出这在某些只需要 MatVec 的应用场景里非常实用。6. 常见问题与排查技巧我踩过的那些坑6.1 残差长期不下降第一反应查什么残差长期不下降是最常遇到的问题我的排查顺序基本是固定的。先看有没有做预处理。如果完全没做预处理而矩阵条件数又很大不下降是正常现象先补一个对角缩放或者 ILU(0)。再看矩阵是否真的对称正定。有个项目里我碰到过一个矩阵从物理模型推导应该是对称的但因为网格编号和边界条件的处理数值上对称性被破坏了一个极小的量导致 CG 行为怪异后来用对称化处理解决。接着检查右端项是否在 A 的值域之外这个问题不一定代表代码有 bug可能只是问题本身无解或者有无穷多解需要检查相容性条件。最后如果用的是 GMRES看看是不是重启步长太小导致 K 空间无法捕捉到关键特征向量信息适当调大 restart比如从 20 调整到 50、100通常有帮助。6.2 内存和时间的平衡大规模场景下的隐形成本Krylov 方法虽然避免了矩阵分解的内存爆炸但迭代本身也有自己的内存成本。GMRES 的基向量随迭代步数线性增长restart 设为 500 时如果矩阵阶数百万只需要几十个基向量内存消耗就明显起来了。这里有个被我反复验证的经验处理超大规模问题优先考虑 BiCGSTAB 或带预处理的 CG而不是 GMRES除非你有充足的内存并且确实需要 GMRES 的稳定性。时间成本方面除了 MatVec 本身最大的开销往往来自向量内积和全局通信。在分布式的 MPI 环境里每一轮 Krylov 迭代通常需要至少一到两次全局归约all-reduce来计算内积这会导致一个分摊开销随进程数增加而恶化。很多大规模的并行迭代求解器性能瓶颈就在这里。如果追求极致扩展性可以考虑 s-step Krylov 方法或者通信避免算法不过这些方法数值稳定性还需要仔细评估目前还不是主流工程选项。6.3 矩阵特征值的快速估计不需要求特征值也能判断收敛很多人会问我想知道这个矩阵适不适合 Krylov 方法但求特征值本身太贵了怎么办其实有更廉价的替代方案用 Lanczos/Arnoldi 过程本身的中间结果做估计。GMRES 的 Hessenberg 矩阵 H_m 的特征值会随着迭代逐渐逼近 A 的部分特征值称之为 Ritz 值。你可以直接在求解过程中把它们打印出来如果看到某些 Ritz 值非常接近 0说明矩阵有接近奇异的模态这几乎一定会拖慢收敛。这个方法不需要额外代价因为 H_m 本来就计算好了。另外矩阵的迹和 Frobenius 范数可以给出特征值分布范围的粗略估计这些只需要 O(nnz) 的扫描即可获得。再配合幂迭代算一个最大特征值的近似你就能在开始正式迭代之前大致判断这个问题的难度级别。6.4 复数矩阵与多右端项问题Krylov 方法在复数域同样适用但要注意两点内积定义需要取共轭另外 GMRES 里的最小二乘问题要在复数意义下理解。实际代码中如果矩阵是 Hermitian 正定的CG 依然可以用只是内积和系数要对应修改。多右端项问题在电磁散射、地震成像等场景里非常常见多个右端项共享同一个矩阵。最朴素的做法是逐个用 Krylov 方法求解但这会浪费大量重复的预处理代价。可以尝试块 Krylov 方法block GMRES、block CG一次迭代同时处理多个右端项矩阵-向量乘也变成矩阵-块乘通常块宽度不大时比如 8-16效率提升明显。不过块方法对块间线性相关性敏感右端项如果有近似线性相关的情况需要做块降秩处理否则数值稳定性会出问题。6.5 舍入误差、除以零与幸运 breakdown浮点环境下的 Krylov 方法最让人头疼的问题之一是 breakdown。它表现为算法某一步突然出现除以零或者接近除以零的情况。在 Arnoldi 过程中这对应 h_{j1,j} 接近零在 BiCGSTAB 中这是 r₀ᵀv₀ 之类的内积接近零。处理原则很简单检测到接近零的情况时不要直接除以它要么提前终止如果这正好说明已经收敛要么做重启。另一种常见问题是 CG 的残差曲线出现微小的上升然后继续下降。这多半不是迭代发散而是舍入误差破坏了共轭性。解决办法是定期做一次重置操作重新计算真正的残差 b-Ax然后用它重新开始方向向量。很多生产级代码里会每隔 50 或 100 步做一次这样操作实际效果非常明显。7. 最后说点实在的我在很长一段时间里犯过一个认知错误觉得 Krylov 方法只要把迭代公式背熟、把库调对就够了。后来被真实问题教育过几次才明白真正的难点从来不在算法本身而在于根据具体的矩阵去选择合适的算法和预处理组合。接触一个新的稀疏矩阵时一定要先花时间弄清楚它的来源、结构特征、条件数和特征值分布再做选型决策这十几分钟的投入通常能让你后面少调几天的参。另外我建议每个学数值分析的人至少手写一遍 GMRES 和 CG。手写不是为了在生产里用而是因为只有亲手实现过 Arnoldi 过程、理解过 Givens 旋转为什么这样构造你才能对文档里那些晦涩的参数有体感。遇到问题的时候你的直觉才不会一片空白。这个领域后续还有不少可以深挖的方向比如针对大规模并行环境的通信避免算法、自适应预处理、机器学习辅助的预处理子选择都是目前很活跃的研究方向。但不管技术怎么演进理解清楚 Krylov 子空间的基本思想和工程约束始终是所有高阶技巧的基础。