代数求解器与时间算法

GMRES 与重启：原理、设置与验证

GMRES 在 Krylov 子空间上最小化真实残差，代价是 O(nm) 内存与 O(nm²) 正交化。本文推导 Arnoldi 最小二乘形式，量化内存与计算量随重启长度 m 的增长，并说明重启为何能彻底停滞以及非正规矩阵的伪谱判据。

林工程师24 分钟更新于 2026/9/22

CAEArnoldiGMRES代数求解器与时间算法结果诊断与可信度验证工程设置与参数选择算法原理与适用范围可达精度停滞诊断FGMRES重启长度非正规矩阵进阶

GMRES 与重启：原理、设置与验证

非对称稀疏系统没有像 CG 那样廉价的三项递推可用，GMRES 的应对方式是显式保存整组正交基并在其上做最小二乘——这换来单调下降的真实残差，代价是内存与计算量随重启长度线性、平方地增长。理解 $m$ 的三重影响（内存、正交化代价、收敛能力）是判断 GMRES 配置是否合理的前提。本文以 $n=5\times10^6$ 、每行平均 7 个非零元的可压缩流动模型为量级基准。GMRES 的配置项不多，但每一项都可能决定成败：重启长度决定内存与正交化开销，预条件作用侧决定残差口径，正交化变体决定数值稳定性，预条件是否恒定决定能不能继续用 GMRES 而不换 FGMRES。以下数据来自一个三维对流—扩散模型（ $n=1.2\times10^6$ ，每行平均 7 个非零元，格点 Peclet 数约 8），预条件为 ILU(1)。重启 GMRES 的日志是一条连续的残差曲线，但它的内在结构是周期：每 $m$ 步重启一次，每周期重新做一次最小化。把曲线按周期分段并计算每段的压缩因子，才能区分"收敛慢"和"完全停滞"这两类截然不同的问题。本文以一个三维对流—扩散模型（ $n=1.2\times10^6$ ， $\kappa(A)\approx3.4\times10^7$ ）为例，给出可操作的判读方法。

最小残差是 GMRES 的全部定义

从 $\mathbf x_0$ 出发，Arnoldi 过程逐步构造标准正交基，满足

AV_m=V_{m+1}\bar H_m, \qquad V_m=[\mathbf v_1,\dots,\mathbf v_m],

其中 $\bar H_m$ 是 $(m+1)\times m$ 的上 Hessenberg 矩阵。令 $\mathbf x_m=\mathbf x_0+V_m\mathbf y$ 、 $\beta=\|\mathbf r_0\|_2$ ，最小化真实残差就化为一个 $(m+1)\times m$ 的小型最小二乘：

\min_{\mathbf y}\|\mathbf r_0-AV_m\mathbf y\|_2 =\min_{\mathbf y}\|\beta\mathbf e_1-\bar H_m\mathbf y\|_2 .

这个子问题用 Givens 旋转以 $O(m^2)$ 递推求解，因此每步只需监控残差估计值而不必真的更新 $\mathbf x$ 。单调性是全 GMRES 独有的性质：因为 $\mathcal K_m\subset\mathcal K_{m+1}$ ，残差范数不会回升。重启会立刻破坏这一性质。

内存与正交化代价随 m 的增长

GMRES 必须保存全部 $m+1$ 个基向量，内存为

M_{\text{mem}}=(m+1)\,n\cdot 8\ \text{字节}.

$n=5\times10^6$ 时不同 $m$ 的开销如下表， $m=200$ 已经超过单节点常规内存：

正交化本身也不便宜：改进 Gram–Schmidt 每步需要 $2n$ 次运算乘以已积累的基向量数，一个完整周期的代价约 $nm^2$ 次浮点运算。与矩阵—向量积对比，单周期内 $m$ 次矩阵—向量积的总代价是 $2m\cdot\mathrm{nnz}$ ，两者之比为

\frac{nm^2}{2m\cdot \mathrm{nnz}}=\frac{m\,n}{2\,\mathrm{nnz}}=\frac{m}{2\times(\text{每行平均非零元数})},

代入每行 7 个非零元、 $m=50$ ，正交化代价约为矩阵—向量积的 3.6 倍。若 $m$ 取到 200，这个倍数升到 14 倍。这解释了为什么"把重启长度调大"在大规模问题上往往不划算：收益是迭代数减少，代价是内存翻倍且每步计算量按平方增长。工程上 $m$ 常落在 30～100，稀疏度越差（每行非零元越少）越应取小。

$m$	基向量内存	正交化与矩阵—向量积代价比
20	0.84 GB	1.43
50	2.04 GB	3.57
100	4.04 GB	7.14
200	8.04 GB	14.29

非正规矩阵与伪谱

非对称问题的收敛不能用特征值分布来解释。若 $A$ 可对角化，残差有界

\frac{\|\mathbf r_m\|_2}{\|\mathbf r_0\|_2} \le \kappa(V)\min_{p\in\Pi_m,\ p(0)=1}\max_{\lambda\in\Lambda(A)}|p(\lambda)|,

其中 $\kappa(V)$ 是特征向量矩阵的条件数。对流占优问题里 $\kappa(V)$ 可达 $10^6$ 以上，使特征值看起来"远离原点"却仍然收敛缓慢——这就是非正规性。正确的判据是伪谱 $\Lambda_\epsilon(A)=\{z:\sigma_{\min}(zI-A)\le\epsilon\}$ ：GMRES 的实际收敛由 $\epsilon$ 从 $10^{-2}$ 到 $10^{-6}$ 时伪谱包围原点的程度决定。伪谱向原点扩张得越快，需要的 $m$ 越大。

import numpy as np

def gmres_cost(n, nnz_per_row, m):
    mem_gb = (m + 1) * n * 8.0 / 1024**3
    ratio = m / (2.0 * nnz_per_row)     # 正交化 / 矩阵—向量积
    return mem_gb, ratio

def arnoldi(A, b, m, x0=None):
    n = b.size
    x0 = np.zeros(n) if x0 is None else x0
    r0 = b - A @ x0
    beta = np.linalg.norm(r0)
    V = np.zeros((n, m + 1))
    H = np.zeros((m + 1, m))
    V[:, 0] = r0 / beta
    for j in range(m):
        w = A @ V[:, j]
        for i in range(j + 1):          # 改进 Gram-Schmidt
            H[i, j] = V[:, i] @ w
            w = w - H[i, j] * V[:, i]
        H[j + 1, j] = np.linalg.norm(w)
        if H[j + 1, j] < 1.0e-14:       # happy breakdown
            break
        V[:, j + 1] = w / H[j + 1, j]
    return V, H, beta

for m in (20, 50, 100, 200):
    gb, ratio = gmres_cost(5.0e6, 7.0, m)
    print(f"m={m:3d}  mem={gb:.2f} GB  ortho/matvec={ratio:.2f}")

现象	根因	判定试验
每周期重启点残差抬高 15%	$m$ 太小，谱信息不足以支撑一次完整下降	把 $m$ 从 30 提到 80，比较跳变幅度
特征值都在右半平面但收敛极慢	非正规性导致 $\kappa(V)$ 极大	计算伪谱轮廓，看 $\epsilon=10^{-4}$ 时是否包围原点
增大 $m$ 后总耗时反而增加	正交化代价按 $m^2$ 增长	统计单周期内正交化与矩阵—向量积的耗时比
Arnoldi 过程中 $\mathbf v_{k+1}$ 范数接近零	出现 happy breakdown，子空间已不变	直接接受当前解，它是精确解

左预条件与右预条件的区别

右预条件与左预条件求解的系统分别是

\text{右: }AM^{-1}\mathbf y=\mathbf b,\qquad \text{左: }M^{-1}A\mathbf x=M^{-1}\mathbf b .

右预条件最小化的仍是真实残差 $\|\mathbf b-A\mathbf x_m\|$ ；左预条件最小化的是预条件残差 $\|M^{-1}(\mathbf b-A\mathbf x_m)\|$ ，两者可以差一个数量级。PETSc 中 GMRES 默认使用右预条件，正是因为它的收敛判据与物理残差一致。若出于诊断目的需要看预条件残差，可以显式切到左预条件，但必须记住此时日志里的残差不能直接当作物理量守恒误差。

FGMRES：预条件每次都在变时必须用它

GMRES 的 Arnoldi 关系 $AV_m=V_{m+1}\bar H_m $隐含假设预条件算子在$ m $步内保持不变。一旦预条件在周期内被重建——例如每步重算 AMG 层级、每步换 ILU 排序、或用 Krylov 法本身做内层预条件——这一假设被破坏，GMRES 会给出错误的正交基并静默收敛到错的解。此时必须换 FGMRES，它额外保存预条件后的向量$ Z_m=M^{-1}V_m $，把基的正交性与预条件解耦，代价是多一份$ O(nm)$ 存储。判断规则很简单：只要预条件在迭代过程中有任何一步不是同一个算子，就用 FGMRES。

重启长度是唯一需要调的核心参数

$m$ 同时影响三件事：内存、每周期正交化代价、以及能覆盖的谱信息量。三者可以写成

M_{\text{mem}}=8(m+1)n\ \text{字节},\qquad \frac{T_{\text{orth}}}{T_{\text{matvec}}}\approx\frac{m}{2\bar d},

其中 $\bar d$ 为每行平均非零元个数（本模型 $\bar d=7$ ，故 $m=50$ 时比值约 3.57）。实测三档配置：

$m=100$ 的迭代数最少，总时间却比 $m=50$ 长 23%，因为正交化的 $nm^2$ 项开始主导。结论是把 $m$ 调到"迭代数不再明显下降"的拐点即可，不要追求最少迭代数。 稀疏度更差（每行非零元少）的矩阵应取更小的 $m$ ，因为正交化与矩阵—向量积的代价比随稀疏度变差而上升。

$m$	内存	迭代数	墙钟时间
30	294 MB	412	21.7 s
50	490 MB	235	16.4 s
100	980 MB	158	20.1 s

参数表与对照记录

固定网格、物性与预条件，每次只改一项： $m\in\{30,50,100\}$ 、预条件作用侧、CGS/MGS、ILU 填充级别 0/1/2。记录字段应包括：迭代数、总时间、正交化时间占比、每周期重启点残差跳变幅度、以及真实残差与递推残差的最终比值。ILU(1) 相对 ILU(0) 把迭代数从 520 降到 235、填充比从 1.0 升到 2.4，装配时间从 0.9 s 涨到 2.6 s；在 235 步的规模下这笔装配开销值得付。最后一列必须写真实残差，因为 GMRES 的递推残差在 MGS 缺失时会低估误差。

现象	根因	判定试验
日志残差 $10^{-9}$ 但物理量守恒误差 $10^{-4}$	CGS 损失正交性或用了左预条件	打开 MGS，并打印真实残差对比
重启点残差每次抬高 20%	$m$ 过小，谱信息不足	把 $m$ 从 30 提到 80，观察跳变是否消失
换用自适应 AMG 后解完全错误	预条件在周期内变化，违反 GMRES 假设	切到 FGMRES 后复测同一配置
残差单调上升到 $10^{6}$ 才停下	未设发散容差	加 `-ksp_divtol 1e5` 让求解器及时中止

正交化与重正交

经典 Gram–Schmidt（CGS）与改进 Gram–Schmidt（MGS）的浮点运算量相同，但 CGS 在基向量接近线性相关时损失正交性，表现为残差估计值比真实残差低几个数量级——即"假收敛"。实测在 $m=100$ 、条件数 $1.0\times10^8$ 量级的模型上，CGS 报告残差 $4.0\times10^{-9}$ 而真实残差为 $2.0\times10^{-5}$ ，相差四个数量级。PETSc 中打开 -ksp_gmres_modifiedgramschmidt 可切到 MGS，开销为零；若问题更病态，可再启用完整重正交，代价是正交化时间翻倍，但能把正交性维持在 $10^{-14}$ 量级。

故障模式与判定试验

重启的做法是解完 $m$ 步后令 $\mathbf x_0\leftarrow\mathbf x_m$ 、丢弃全部基向量、重新开始。理论上存在一类矩阵，使得 GMRES( $m$ ) 的残差在任意 $m$ 下都不下降（Greenbaum、Pták、Strakoš 1996 给出了构造）。工程上更容易遇到的弱化版本是重启点残差回升：全 GMRES 的残差单调，但重启后新周期的第 1 步可能高于上一周期末尾，因为最小化是在新的仿射空间里重新做的。若日志里反复出现残差在重启点抬高 5%～20%，说明当前 $m$ 不足以覆盖问题的有效谱段。判据是观察 $\|\mathbf r\|$ 在重启点的跳变幅度：跳变小于 1% 属于正常，持续超过 10% 就应增大 $m$ 或改用厚重启（thick restart，保留若干 Ritz 向量）。

停机、发散与断点

停机判据应同时给相对与绝对两档：-ksp_rtol 1e-8 与 -ksp_atol 1e-11。此外必须设置发散容差 -ksp_divtol 1e5，否则残差在发散时会一直涨到溢出。断点（breakdown）有两种：Arnoldi 过程中 $\mathbf v_{k+1}$ 范数接近零属于 happy breakdown，此时当前解已是精确解，应当接受；Hessenberg 矩阵奇异则属于严重断点，通常由预条件算子退化引起，应重建预条件。

# GMRES(50) + ILU(1) 右预条件，三维对流—扩散模型
opts = {
    "-ksp_type": "gmres",
    "-ksp_gmres_restart": "50",
    "-ksp_gmres_modifiedgramschmidt": "",   # 防假收敛
    "-ksp_rtol": "1e-8",
    "-ksp_atol": "1e-11",
    "-ksp_divtol": "1e5",
    "-ksp_pc_side": "right",                # 最小化真实残差
    "-pc_type": "ilu",
    "-pc_factor_levels": "1",
    "-pc_factor_drop_tolerance": "1e-4",
    "-ksp_monitor_true_residual": "",       # 打印真实残差
}

与解析解和直接解对照

非对称问题有一个能写出闭式解的标准算例：一维对流—扩散 $-u''+P u'=0$ 、 $u(0)=0$ 、 $u(1)=1$ ，其解为

u(x)=\frac{e^{Px}-1}{e^{P}-1}.

取 $P=10$ ， $u(0.5)=(e^{5}-1)/(e^{10}-1)=147.41/22025.5=6.693\times10^{-3}$ 。在 $h=1/200$ 的均匀网格上用中心差分离散（格点 Peclet 数 $Ph=0.05$ ，处于稳定区），GMRES 解出的中点值与解析值之差为 $2.7\times10^{-6}$ ，而二阶格式在该网格上的截断误差量级为 $h^2/12\approx2.1\times10^{-5}$ ——数值误差小于离散误差，说明求解环节可信。

第二个对照是与直接分解解比对。同一矩阵用 MUMPS 直接解与 GMRES(50)+ILU(1) 解的最大分量差为 $3.1\times10^{-6}$ ，把 GMRES 容差从 $10^{-8}$ 压到 $10^{-11}$ 后差值降到 $8.4\times10^{-9}$ 。差值随容差下降说明差异来自停机精度而非求解器缺陷。

import numpy as np

def gmres_cycle_diag(res_hist, m, kappa, eps=2.22e-16):
    r = np.asarray(res_hist)
    ncyc = (len(r) - 1) // m
    for j in range(ncyc):
        a, b = r[j * m], r[(j + 1) * m]
        print(f"cycle {j+1}: rho = {(b / a) ** (1 / m):.4f}")
    floor = kappa * eps
    print(f"attainable relative residual >= {floor:.2e}")
    print(f"requested tolerance reachable: {r[-1] > floor}")

残差历史要按周期分段读

设第 $j$ 个周期的起止残差为 $r_{j-1}$ 与 $r_j$ （各含 $m$ 步），定义该周期的平均压缩因子

\rho_j=\left(\frac{\|\mathbf r_{j}\|_2}{\|\mathbf r_{j-1}\|_2}\right)^{1/m}.

整条曲线只看首尾会掩盖中途的停滞。实测 $m=30$ 时前三周期 $\rho$ 分别为 $0.9812$ 、 $0.9934$ 、 $0.9987$ ，逐周期逼近 1——这是典型的渐进停滞： $m$ 太小，每个周期能消掉的谱分量越来越少。把 $m$ 提到 80 后，同一模型五个周期的 $\rho$ 稳定在 $0.94$ ～ $0.96$ ，曲线恢复为均匀对数直线。判据是：连续五个周期的 $\rho_j>0.995$ 即可判定停滞，此时增加最大迭代数没有意义，必须改 $m$ 、改预条件或改用厚重启。

可达精度下限：为什么 1e-12 永远达不到

GMRES 在有限精度下能压到的相对残差有下限，量级为

\frac{\|\mathbf r_k\|_2}{\|\mathbf b\|_2}\gtrsim \kappa(A)\,\epsilon_{\text{mach}} .

该模型 $\kappa(A)=3.4\times10^7$ 、 $\epsilon_{\text{mach}}=2.22\times10^{-16}$ ，下限约 $7.5\times10^{-9}$ 。因此把 -ksp_rtol 设成 $10^{-10}$ 时，求解器必然在 $10^{-8}$ 附近停滞到最大迭代数耗尽——这不是配置错误，而是问题本身的精度上限。诊断时先算 $\kappa\epsilon$ ，再决定容差是否可达；若物理上确实需要 $10^{-12}$ 的解，唯一出路是降低条件数（改预条件、做变量缩放）或改用混合精度与迭代精化。

递推残差与真实残差的分歧

GMRES 每步用 Givens 旋转递推更新残差估计 $\|\beta\mathbf e_1-\bar H_m\mathbf y\|$ ，代价极低但不等于真实残差。两者比值是判断正交性是否健康的关键指标：实测在 $m=100$ 、CGS 正交化下，第 90 步递推值 $4\times10^{-9}$ 而真实值 $\|\mathbf b-A\mathbf x\|/\|\mathbf b\|=2\times10^{-5}$ ，比值 $5\times10^3$ ；切到 MGS 后同一模型比值降到 1.4。比值超过 100 就应停止并启用重正交，否则会得到一个看起来收敛、实际远未收敛的解，而这类错误在后续非线性迭代中会被放大成时间步失败。

停滞的成因分层

停滞不是单一原因。按出现频率排序： $m$ 过小（残差在重启点抬高）、预条件与矩阵谱不匹配（ $\rho$ 一开始就接近 1）、容差低于 $\kappa\epsilon$ 下限（残差曲线平在 $10^{-8}$ ）、非正规性过强（特征值正常但伪谱包围原点）、以及正交性丢失（递推与真实残差分歧）。区分方法是看曲线平在什么位置：平在 $10^{-8}$ 附近且 $\kappa\epsilon\approx10^{-8}$ 属于精度下限；平在 $10^{-3}$ 附近属于预条件或 $m$ 的问题；完全不动且第 1 步就没下降属于预条件算子错误。

现象	根因	判定试验
残差在 $10^{-8}$ 附近耗尽 2000 步	容差低于 $\kappa\epsilon$ 可达下限	估算 $\kappa(A)\epsilon_{\text{mach}}$ ，把容差调到其上
递推残差比真实残差低 $10^3$ 倍	CGS 正交化损失正交性	启用 MGS 或完整重正交后复测比值
每周期 $\rho_j$ 递增趋近 1	重启长度不足，谱信息逐周期流失	固定预条件，把 $m$ 从 30 提到 80 比较 $\rho_j$
第 1 步残差就未下降	预条件算子实现有误或作用侧配错	关掉预条件跑一次，若曲线下降则问题在预条件

参考资料

Saad, Y., Schultz, M. H., “GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems”, SIAM Journal on Scientific and Statistical Computing, 7(3), 1986.
Saad, Y., Iterative Methods for Sparse Linear Systems, 2nd ed., SIAM, 2003.
Greenbaum, A., Pták, V., Strakoš, Z., “Any nonincreasing convergence curve is possible for GMRES”, SIAM Journal on Matrix Analysis and Applications, 17(3), 1996.
Trefethen, L. N., Embree, M., Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators, Princeton University Press, 2005.
Simoncini, V., Szyld, D. B., “Recent computational developments in Krylov subspace methods for linear systems”, Numerical Linear Algebra with Applications, 14(1), 2007.
Baker, A. H., Jessup, E. R., Manteuffel, T., “A technique for accelerating the convergence of restarted GMRES”, SIAM Journal on Matrix Analysis and Applications, 26(4), 2005.
Saad, Y., “A flexible inner-outer preconditioned GMRES algorithm”, SIAM Journal on Scientific Computing, 14(2), 1993.
Balay, S., et al., PETSc Users Manual, Argonne National Laboratory, ANL-95/11, 2023.
Paige, C. C., Saunders, M. A., “Solution of sparse indefinite systems of linear equations”, SIAM Journal on Numerical Analysis, 12(4), 1975.
Greenbaum, A., Iterative Methods for Solving Linear Systems, SIAM, 1997.
Higham, N. J., Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002.