企业官网建设流程全解析

简介：这是数值线性代数领域关于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) r0，p_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；对称不定用 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, rtol=1e-8, maxiter=500, atol=0) x_gmres, info_gmres = gmres(A, b, rtol=1e-8, maxiter=500, atol=0) x_minres, info_minres = minres(A, b, rtol=1e-8, maxiter=500, atol=0) x_bicg, info_bicg = bicgstab(A, b, rtol=1e-8, maxiter=500, atol=0) for name, info in [("CG", info_cg), ("GMRES", info_gmres), ("MINRES", info_minres), ("BiCGSTAB", info_bicg)]: print(f"{name}: converged={info == 0}")

info返回值是理解 scipy 迭代器行为的关键：info=0表示收敛，info>0表示达到最大迭代步数未收敛，info<0表示输入非法或算法破裂。测试对称正定矩阵时，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_tol=1e-4, fill_factor=10) # 把预处理子包装成 M^{-1} 的调用形式 def apply_precond(v): return ilu.solve(v) # 在 GMRES 中传入预处理函数 x, info = spla.gmres(A, b, rtol=1e-10, maxiter=200, M=apply_precond)

这里的drop_tol控制丢弃小元素的门槛，fill_factor控制允许的填充量上限。两者共同决定预处理子的质量和构建成本。实际调参时，先固定drop_tol=1e-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的设定要考虑最终用途：做特征值问题的内迭代，rtol=1e-3也许就够；做高精度结构分析，可能需要rtol=1e-12。不要在需求未知时盲目追求高精度，这会平白增加几十步迭代。

6. 病态矩阵处理与 Krylov 子空间迭代算法的规模化实践

6.1 矩阵重排序：Cuthill-McKee 对 ILU 的影响

预处理子的质量取决于矩阵结构，而结构可以通过重排改变。Reverse Cuthill-McKee（RCM）是一种带宽缩减算法，它把矩阵重新排序成轮廓更窄的形式。对于有限元网格生成的矩阵，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_mode=True) # 重排矩阵和右端项 A_rcm = A_csr[perm][:, perm] b_rcm = b[perm] # 在重排后的系统上做迭代，最后把解映射回原次序 x_rcm, info = spla.gmres(A_rcm, b_rcm, rtol=1e-10, M=ilu_precond) x = np.empty_like(b) x[perm] = x_rcm

perm是重排后的下标映射数组，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 子空间迭代算法真的收敛到了该有的精度。

本文还有配套的精品资源，点击获取