广义最小残量法

广义最小残量法（generalized minimal residual method，GMRES）在 Krylov 子空间中选择使线性方程残差二范数最小的近似解，适用于一般非对称非奇异系统。它用 Arnoldi 正交化建立基，再解一个较小的最小二乘问题。未重启时搜索空间不断扩大，存储和正交化成本也随迭代增长。

一个非对称系统的首步

取 $A=\left(\begin{array}{cc}3& 1\\ 0& 2\end{array}\right)$、$b=(1,1{)}^T$、$x_0=0$。第一搜索空间只有 $b$ 的倍数，令 $x=\alpha b$。残差为 $(1-4\alpha ,1-2\alpha {)}^T$，其平方范数为 $$(1-4\alpha {)}^2+(1-2\alpha {)}^2=20{\alpha }^2-12\alpha +2.$$ 求极小得 $\alpha =3/10$，所以 $x_1=(0.3,0.3{)}^T$、$r_1=(-0.2,0.4{)}^T$、$\Vert r_1{\Vert }_2=\sqrt{0.2}$。第二搜索空间由 $b,Ab$ 张成整个平面，故可得到 $x_{\ast }=(1/6,1/2{)}^T$。与 CG 不同，这里最小化的是残差二范数，矩阵也不要求对称正定。

Arnoldi 关系

记初始残差 $r_0=b-A{x}_0$，$\beta =\Vert r_0{\Vert }_2$。若 $\beta =0$ 就已解出，否则令 $v_1=r_0/\beta$。依次对 $A{v}_j$ 减去它在已有正交基上的投影： $$h_{ij}=v_i^{T}A{v}_j,w=A{v}_j-\sum _{i=1}^j{h}_{ij}{v}_i,h_{j+1,j}=\Vert w{\Vert }_2.$$ 若最后范数非零，令 $v_{j+1}=w/h_{j+1,j}$。收集为 $V_m=[v_1,\dots ,v_m]$，便有 $$A{V}_m=V_{m+1}{\overline{H}}_m,$$ 其中 ${\overline{H}}_m$ 为 $(m+1)\times m$ 上 Hessenberg 矩阵。这个关系是每列投影恒等式的集合，不是对 $A$ 作完整稠密分解。SciPy 的 GMRES 文档说明了一般系统、重启和真实残差检验的接口约定。

小最小二乘问题与最优性证明

在 $x_0+{𝒦}_m(A,r_0)$ 中写 $x=x_0+V_{m}y$。利用 $r_0=V_{m+1}\beta e_1$ 和 Arnoldi 关系，得 $$r=b-Ax=V_{m+1}(\beta e_1-{\overline{H}}_{m}y).$$ 列正交性给出 $\Vert V_{m+1}z{\Vert }_2=\Vert z{\Vert }_2$，因此大系统的最小残差问题等价于 $$\underset{y\in {ℝ}^m}{min}\Vert \beta e_1-{\overline{H}}_{m}y{\Vert }_2.$$ 这是 GMRES 最优性的完整来源。搜索空间嵌套，所以下一步可行集合包含上一步，精确算术下未重启残差范数不增加；它可能暂时不变，并不保证每轮严格下降。

首步在 A乘K₁ 所在直线上投影 b。点 Ax₁=(1.2,0.6)，残差 b−Ax₁=(−0.2,0.4) 与方向 (4,2) 垂直。

主例中 $v_1=(1,1{)}^T/\sqrt{2}$，$A{v}_1=(4,2{)}^T/\sqrt{2}$，故 $h_{11}=3,h_{21}=1,v_2=(1,-1{)}^T/\sqrt{2}$。首个小问题为最小化 $\Vert (\sqrt{2},0{)}^T-(3,1{)}^{T}y\Vert$，其解 $y=3\sqrt{2}/10$，再次得到 $x_1=(0.3,0.3{)}^T$。第二列满足 $A{v}_2=2v_2$，新方向范数为零：Krylov 空间已封闭。对非奇异 $A$，这种精确封闭意味着当前空间中已有精确解，称为 happy breakdown；浮点中仍须核对真实残差。

重启、预条件与收敛边界

完整 GMRES 保存越来越多的基向量，通常每 $m$ 步重启，用当前近似重新开始。重启限制存储，却丢弃此前搜索方向；即使原矩阵非奇异，也可能出现停滞或很慢的收敛。非正规矩阵的行为也不能只看特征值分布。

左预条件最小化的是预条件后的残差；右预条件可保持原始残差最小化的解释，但解空间变量随之改变。无论采用哪一种，验收最终都应查看 $b-Ax$ 的约定范数。若预条件子随步改变，需使用 flexible GMRES 等相应变体。

透明实现与复算

以下 NumPy 演示不重启，对实矩阵使用两遍修正 Gram–Schmidt，并每步直接解小最小二乘。它返回近似解及真实残差历史，不修改输入。若达到封闭空间却未满足容限，会报告失败，避免将数值 breakdown 当作成功。

矩阵为稠密时，$m$ 步矩阵向量积需 $O(m{n}^2)$，正交化需 $O(n{m}^2)$，存储基需 $O(nm)$。本教学实现反复调用通用最小二乘，在 $m$ 增大时另有约 $O(m^4)$ 累积小问题成本；生产实现通常用增量 Givens 旋转避免这一重复。

Python 3

import numpy as np


def gmres_demo(A, b, tol=1e-12):
    A = np.array(A, dtype=float, copy=True)
    b = np.array(b, dtype=float, copy=True)
    if b.ndim != 1 or not b.size or A.shape != (b.size, b.size):
        raise ValueError("square matrix and matching vector required")
    if not np.isfinite(A).all() or not np.isfinite(b).all() or not np.isfinite(tol) or tol <= 0:
        raise ValueError("finite data and positive tolerance required")
    n = b.size; beta = np.linalg.norm(b)
    if not np.isfinite(beta):
        raise ArithmeticError("initial residual norm outside floating-point range")
    target = tol*max(1.0, beta)
    if not np.isfinite(target) or not np.isfinite(np.linalg.norm(A)):
        raise ArithmeticError("residual threshold or matrix norm outside floating-point range")
    if beta <= target:
        return np.zeros(n), [float(beta)]
    V = np.zeros((n, n+1)); H = np.zeros((n+1, n))
    V[:, 0] = b/beta; history = [float(beta)]
    for j in range(n):
        w = A @ V[:, j]
        for _ in range(2):
            for i in range(j+1):
                hij = np.dot(V[:, i], w)
                H[i, j] += hij
                w -= hij*V[:, i]
        H[j+1, j] = np.linalg.norm(w)
        if not np.isfinite(H[:j+2, :j+1]).all() or not np.isfinite(w).all():
            raise ArithmeticError("nonfinite Arnoldi coefficients")
        rhs = np.zeros(j+2); rhs[0] = beta
        y = np.linalg.lstsq(H[:j+2, :j+1], rhs, rcond=None)[0]
        x = V[:, :j+1] @ y
        residual = np.linalg.norm(b-A@x)
        if not np.isfinite(residual):
            raise ArithmeticError("nonfinite residual")
        history.append(float(residual))
        if residual <= target:
            return x, history
        if H[j+1, j] <= np.finfo(float).eps*max(1.0, np.linalg.norm(A)):
            raise RuntimeError("Krylov breakdown before verified convergence")
        V[:, j+1] = w/H[j+1, j]
    raise RuntimeError("true residual did not reach tolerance")


x, history = gmres_demo([[3, 1], [0, 2]], [1, 1])
print(x)        # [0.16666667 0.5]
print(history)  # sqrt(2), sqrt(0.2), roundoff-sized final residual

有解变式取 $A=2I$、任意非零 $b$。第一 Krylov 方向就含精确解 $x=b/2$，Arnoldi 新向量为零且残差也为零；这是一种正确的有限终止。若改成奇异不相容系统 $A=0,b\ne 0$，同样出现新方向为零，残差却无法下降，代码应报告失败。这区分了封闭空间与已解出原问题。