GMRES 算法

科学和工程领域中许多最重大的挑战，最终都归结为求解一个巨大的线性方程组，其形式表示为 $A\boldsymbol{x}=\boldsymbol{b}$。当矩阵 $A$ 涉及数百万甚至数十亿个变量时，直接求解 $\boldsymbol{x}$ 在计算上是不可能的。我们无法用蛮力来攻克这些问题；相反，我们需要一种能够逐步找到解的智能而高效的方法。这正是广义最小残差 (GMRES) 算法所扮演的角色，它是一种强大的迭代方法，已成为现代科学计算的基石。

本文旨在解决直接方法失效时，求解大规模线性系统这一根本性挑战。文章全面概述了 GMRES 方法，引导读者从其基本原理走向其在现实世界中的影响。在“原理与机制”一章中，您将学习 GMRES 如何巧妙地构建一个特殊的搜索空间——Krylov 子空间——并利用优雅的 Arnoldi 迭代在每一步找到最佳可能解。随后，“应用与跨学科联系”一章将展示该方法的多功能性，探讨其在流体力学中求解偏微分方程的应用、其与其他计算技术的关系，以及其作为先进非线性求解器中关键组件的角色。

想象一下，你正在试图解决一个巨大的难题，比如一个包含数百万个线性方程的方程组，其紧凑形式为 $A\boldsymbol{x}=\boldsymbol{b}$。矩阵 $A$ 是一个庞然大物，一个 $n \times n$ 的巨型矩阵，其中 $n$ 可以是数百万甚至数十亿。试图通过直接对 $A$ 求逆来找到解 $\boldsymbol{x}$，不仅计算成本高昂；对于人类面临的最大规模问题，这在物理上是不可能的。这会比宇宙的年龄还要长，需要比地球上所有计算机加起来还要多的内存。我们无法用蛮力征服这个巨人。我们必须运用智慧。

这就是广义最小残差方法 (GMRES) 登场的时刻，它不像一根攻城槌，而是一位敏捷而聪明的探险家。它不会试图一次性解决整个难题。相反，它从一个猜测 $\boldsymbol{x}_0$ 开始，踏上一段旅程，采取一系列步骤，每一步都以智能的方式更接近真实解。

任何优秀的探险者都会问的第一个问题是：“我错了多少？”在线性代数的世界里，这由​​残差向量​​ $\boldsymbol{r} = \boldsymbol{b} - A\boldsymbol{x}$ 来衡量。如果我们的猜测 $\boldsymbol{x}$ 是完美的，$A\boldsymbol{x}$ 将等于 $\boldsymbol{b}$，残差 $\boldsymbol{r}$ 将是一个全零向量。这个向量的大小，即​​范数​​ $\|\boldsymbol{r}\|_2$，精确地告诉我们距离解有多远。GMRES 的目标简单而直观：在每一步，都将这个残差变得尽可能小。

让我们从初始猜测 $\boldsymbol{x}_0$ 开始。它可能不太好。它给了我们一个初始残差 $\boldsymbol{r}_0 = \boldsymbol{b} - A\boldsymbol{x}_0$。要改进我们的猜测，最明显的行进方向是什么？也许我们应该沿着残差本身的方向移动。我们可以将下一个猜测定义为 $\boldsymbol{x}_1 = \boldsymbol{x}_0 + \alpha \boldsymbol{r}_0$，其中 $\alpha$ 是我们选择的步长，以使新的残差 $\|\boldsymbol{b} - A(\boldsymbol{x}_0 + \alpha \boldsymbol{r}_0)\|_2$ 尽可能小。这是一个合理的开始，但这就像在浓雾中徒步，只沿着感觉最陡峭的下坡方向迈步。你可能会取得进展，但你忽略了大量关于地形的信息。

GMRES 的目标远不止于此。它决定不仅探索一个方向，而是探索一整个充满希望的方向的“子空间”。哪些方向是有希望的呢？嗯，$\boldsymbol{r}_0$ 是一个好的开始。但矩阵 $A$ 本身包含了我们问题的所有信息。将它应用于 $\boldsymbol{r}_0$ 会得到一个新向量 $A\boldsymbol{r}_0$，它告诉我们问题的几何结构如何“扭曲”残差。为什么不同时也在那个方向上搜索呢？那 $A(A\boldsymbol{r}_0)$，或者说 $A^2\boldsymbol{r}_0$ 呢？

这就引出了该方法的核心：​​Krylov 子空间​​。在第 $k$ 步，GMRES 不仅仅是沿着一条线搜索。它在整个仿射子空间 $\boldsymbol{x}_0 + \mathcal{K}_k(A, \boldsymbol{r}_0)$ 内搜索，其中 $\mathcal{K}_k(A, \boldsymbol{r}_0)$ 是由前 $k$ 个 Krylov 向量张成的空间： $\mathcal{K}_k(A, \boldsymbol{r}_0) = \text{span}\{\boldsymbol{r}_0, A\boldsymbol{r}_0, A^2\boldsymbol{r}_0, \dots, A^{k-1}\boldsymbol{r}_0\}$ 在每一步，GMRES 都会提出一个强有力的问题：在我能构建的所有向量中——即从 $\boldsymbol{x}_0$ 出发，加上这些前 $k$ 个 Krylov 向量的某种组合——哪一个能让我最接近解？“最接近”的定义是唯一重要的那个：哪一个能产生​​最小残差​​？这种最优性是 GMRES 的灵魂所在。

现在，我们有了一个绝妙的策略，但实践起来却是一场噩梦。原始的 Krylov 向量 $\{\boldsymbol{r}_0, A\boldsymbol{r}_0, \dots\}$ 是一套糟糕的工具。随着我们生成更多的向量，它们往往指向非常相似的方向，变得几乎线性相关。使用它们就像试图用一桶生锈弯曲的钉子进行精密工程。我们需要一套原始的工具：一组相互垂直（​​正交​​）且长度为一（​​单位化​​）的向量基。我们需要一个​​标准正交基​​。

这正是 ​​Arnoldi 迭代​​ 的任务，它是一个优美而高效的程序，扮演着 GMRES 的工匠大师的角色。它接收那些笨拙的 Krylov 向量，并逐一将它们锻造成一个完美的标准正交基 $\{\boldsymbol{v}_1, \boldsymbol{v}_2, \dots, \boldsymbol{v}_k\}$，用于同一个 Krylov 子空间 $\mathcal{K}_k$。

这个过程是一种优雅的提纯形式。

1. 它首先简单地将初始残差单位化：$\boldsymbol{v}_1 = \boldsymbol{r}_0 / \|\boldsymbol{r}_0\|_2$。
2. 为了得到第二个向量，它取下一个 Krylov 向量 $A\boldsymbol{v}_1$，并减去其中已经存在于 $\boldsymbol{v}_1$ 方向上的所有分量。剩下的部分，根据定义，与 $\boldsymbol{v}_1$ 正交。然后它将这个新向量单位化以得到 $\boldsymbol{v}_2$。
3. 为了得到 $\boldsymbol{v}_3$，它取 $A\boldsymbol{v}_2$，并减去所有存在于 $\boldsymbol{v}_1$ 和 $\boldsymbol{v}_2$ 方向上的分量。剩下的部分与两者都正交。将其单位化，你就得到了 $\boldsymbol{v}_3$。

这个过程持续进行，每个新向量 $\boldsymbol{v}_{j+1}$ 都是由 $A\boldsymbol{v}_j$ 通过仔细移除其在所有先前向量 $\{\boldsymbol{v}_1, \dots, \boldsymbol{v}_j\}$ 上的投影而创建的。值得注意的是，这是一个“即用即付”的系统。在第 $k$ 步，我们只需要计算一个新的矩阵-向量乘积 $A\boldsymbol{v}_k$，这通常是计算中最昂贵的部分。

但 Arnoldi 的真正魔力在于它能一箭双雕。在构建标准正交基向量 $V_k = [\boldsymbol{v}_1 | \dots | \boldsymbol{v}_k]$ 的同时，它还记录了减法过程中使用的系数。这些系数构成一个小的 $(k+1) \times k$ 矩阵 $\bar{H}_k$，称为​​上海森堡矩阵 (upper Hessenberg matrix)​​。这个矩阵具有特殊的结构，其第一副对角线下方均为零。它是关于巨大矩阵 $A$ 在我们子空间上作用的一个紧凑的“缩略图”。这种关系被一个单一、优雅的方程所捕捉：$AV_k = V_{k+1}\bar{H}_k$。这告诉我们，巨大且未知的 $A$ 在我们基上的作用，可以被微小且已知的 $\bar{H}_k$ 在相同基上的作用完美描述。举一个具体的例子，在一个示例问题上执行两步 Arnoldi 迭代后，我们可能会得到一个微小的 $3 \times 2$ 矩阵，就像在 中找到的那样。

有趣的是，如果原始矩阵 $A$ 恰好是对称的，问题的结构会得到极大的简化。Hessenberg 矩阵 $\bar{H}_k$ 会变成三对角矩阵，而 Arnoldi 过程会简化为著名的 ​​Lanczos 迭代​​，每一步所需的工作量要少得多。这是物理学和数学中一个反复出现的主题：对称性简化一切。

随着 Arnoldi 过程的完成，GMRES 现在可以施展其神来之笔。我们对最佳解 $\boldsymbol{x}_k$ 的搜索，是在表达式 $\boldsymbol{x}_k = \boldsymbol{x}_0 + V_k \boldsymbol{y}_k$ 中寻找最佳系数集 $\boldsymbol{y}_k$。问题是最小化残差范数： $\|\boldsymbol{r}_k\|_2 = \|\boldsymbol{b} - A \boldsymbol{x}_k\|_2 = \|\boldsymbol{b} - A(\boldsymbol{x}_0 + V_k \boldsymbol{y}_k)\|_2$ 利用我们的初始残差 $\boldsymbol{r}_0 = \boldsymbol{b} - A\boldsymbol{x}_0$ 和神奇的 Arnoldi 关系 $AV_k = V_{k+1}\bar{H}_k$，这个表达式经历了一次惊人的转变： $\|\boldsymbol{r}_k\|_2 = \|\boldsymbol{r}_0 - A V_k \boldsymbol{y}_k\|_2 = \|\beta \boldsymbol{v}_1 - V_{k+1}\bar{H}_k \boldsymbol{y}_k\|_2$ 其中 $\beta = \|\boldsymbol{r}_0\|_2$。由于 $\boldsymbol{v}_1$ 只是矩阵 $V_{k+1}$ 的第一列，我们可以将其写为 $V_{k+1} \boldsymbol{e}_1$，其中 $\boldsymbol{e}_1$ 是一个在第一个位置为 1、其余位置为零的向量。 $\|\boldsymbol{r}_k\|_2 = \|V_{k+1}(\beta \boldsymbol{e}_1 - \bar{H}_k \boldsymbol{y}_k)\|_2$ 因为 $V_{k+1}$ 是一个标准正交矩阵，乘以它不会改变向量的长度——它只是旋转了向量。所以，最小化上述表达式等同于最小化括号内向量的长度！

最初那个不可能解决的、在 $n$ 维空间中寻找最佳 $\boldsymbol{x}_k$ 的问题，已经转化为一个简单的问题：找到短向量 $\boldsymbol{y}_k$ 以最小化 $\|\beta \boldsymbol{e}_1 - \bar{H}_k \boldsymbol{y}_k\|_2$。这是一个可以高效求解的小型 $(k+1) \times k$ 最小二乘问题。这就是使 GMRES 如此强大的“技巧”。我们用一个微小的问题替代了一个极其庞大的问题，而这个小问题的解恰好能给我们所需的东西。

Krylov 子空间 $\mathcal{K}_k$ 随每次迭代而增长。在一个 $n$ 维空间中，这种增长不可能永远持续下去。Krylov 子空间的维度最多为 $n$。理论上，当 $k=n$ 时，子空间 $\mathcal{K}_n(A, \boldsymbol{r}_0)$ 必须张成整个空间 $\mathbb{R}^n$（或者如果在一个更小的子空间中找到解，则会更早终止）。如果你的搜索空间是整个宇宙，你保证能找到你正在寻找的物体。因此，在精确算术中，完整的 GMRES 方法保证在至多 $n$ 次迭代内找到精确解。

这个理论上的保证非常美妙，但代价高昂。在每一步 $k$，Arnoldi 过程都需要存储所有 $k$ 个基向量并执行 $k$ 次正交化步骤。如果 $n$ 是百万级别，我们根本无法承担运行算法数千步的成本。内存和计算成本将变得天文数字。

这导致了一个务实但哲学上代价高昂的妥协：​​重启动 GMRES(m)​​。在这里，我们让优雅的 GMRES 过程运行一个固定的、可管理的步数，比如 $m=50$。我们在这个有限的 50 维子空间中找到最佳解 $\boldsymbol{x}_m$。然后，我们做一个残酷的决定：我们将 $\boldsymbol{x}_m$ 声明为新的起始猜测，并丢弃整个基，即所有关于问题几何形状的累积知识。我们从新的残差 $\boldsymbol{b}-A\boldsymbol{x}_m$ 开始，从头构建一个全新的 Krylov 子空间。

这种“失忆”打破了理论上的收敛保证。通过丢弃旧的搜索方向，算法可能会变得短视。有一些著名的“病态”案例，其中重启动的 GMRES 可能会完全停滞，毫无进展。对于某些问题，取得进展所需的方向只有在探索超过 $m$ 步后才“可见”。如果我们在第 $m$ 步重启，我们就在找到那个关键方向之前将其丢弃，并且可能一次又一次地循环往复，永无止境 [@problemid:2183305]。这种理论完美与实际需求之间的张力是计算科学中的一个核心主题。

如果我们被迫使用短视的重启动 GMRES，我们至少能帮助它看得更清楚吗？答案是肯定的，通过​​预处理​​的艺术。其思想是找到一个“容易”的矩阵 $M$，它是我们“困难”矩阵 $A$ 的一个良好近似。所谓“容易”，是指用 $M$ 求解系统，即计算 $M^{-1}\boldsymbol{z}$，是快速的。然后我们用 $M$ 将原始问题转化为一个对 GMRES 来说更容易解决的问题。

主要有两种方法可以做到这一点：

1. ​​左预处理：​​ 我们求解系统 $M^{-1}A\boldsymbol{x} = M^{-1}\boldsymbol{b}$。GMRES 现在应用于矩阵 $M^{-1}A$，该矩阵有望比 $A$“更好”（其特征值更聚集在 1 附近）。问题在于，GMRES 现在最小化的是*预处理后*的残差范数 $\|M^{-1}(\boldsymbol{b}-A\boldsymbol{x}_k)\|_2$，而不是真实的残差。这可能会产生误导，因为一个小的预处理后残差并不总能保证一个小的真实残差。

2. ​​右预处理：​​ 我们引入一个新变量 $\boldsymbol{y}$ 并求解系统 $AM^{-1}\boldsymbol{y} = \boldsymbol{b}$ 以得到 $\boldsymbol{y}$。一旦我们有了 $\boldsymbol{y}_k$，我们通过 $\boldsymbol{x}_k = M^{-1}\boldsymbol{y}_k$ 恢复我们的解。这种方法的精妙之处既微妙又深刻。GMRES 算法最小化它正在求解的系统的残差，即 $\|\boldsymbol{b} - AM^{-1}\boldsymbol{y}_k\|_2$。但由于 $\boldsymbol{x}_k = M^{-1}\boldsymbol{y}_k$，这恰好等于 $\|\boldsymbol{b} - A\boldsymbol{x}_k\|_2$——即真实的残差！这意味着右预处理允许我们用一个性态更好的矩阵来引导算法，同时仍然监控我们朝向解的实际、真实进展。

归根结底，GMRES 不仅仅是一个算法；它是一种哲学。它告诉我们，通过结合一个最优但局部的策略（最小残差）、一个优美构造性的工具（Arnoldi 迭代）、一个巧妙的视角转换（求解小系统），以及实践智慧（重启动和预处理），我们能够成功地驾驭并解决真正天文数字规模的问题。

我们花了一些时间欣赏广义最小残差方法那精巧的机制。我们已经看到它如何巧妙地在一个高维空间的浩瀚中航行，每一步都构建一个特殊的“Krylov”子空间，以找到对我们解的最佳可能改进。这是一件优雅的数学艺术品。但在科学中，艺术不应被局限于画廊。它应该被使用，被应用，以帮助我们理解世界。现在，让我们把这个美丽的工具从盒子里拿出来，看看它能做什么。我们会发现，其逻辑的回响在令人惊讶的广泛的科学和工程学科中产生共鸣，其方式常常是出人意料且深刻的。

许多物理学的基本定律都是用偏微分方程 (PDE) 的语言写成的。这些方程描述了像热量、动量和电势这样的量如何在空间和时间中变化。要在计算机上求解它们，我们必须将它们从微积分的连续语言翻译成代数的离散语言。这种翻译，或称“离散化”，几乎总是导致一个线性方程组，其规模通常是巨大的。这个系统的矩阵，我们称之为 $A$，是其背后物理学的数字指纹。正是在这里，GMRES 找到了它最自然的归宿。

考虑流体的流动，这是计算流体力学 (CFD) 的核心问题。流体可能携带热量或污染物。这种“物质”由于扩散（一种随机行走过程）而散开，并且也被流体的整体运动所携带，这个过程称为对流。这两种效应之间的平衡由一个单一的无量纲数——Péclet 数——来捕捉。当扩散占主导地位时（低 Péclet 数），问题性态良好，产生的矩阵 $A$ 近似对称。但是当对流占主导地位时——如在一条有污染源的快速流动的河流中——矩阵会变得高度非对称。为什么？因为任何一点的信息现在都受到“上游”发生的事情的强烈影响。这种物理上的不对称性被矩阵直接继承。对于 GMRES 来说，一个高度非对称的矩阵是一个巨大的挑战。收敛速度可能会慢如蜗牛，这种现象被称为停滞。

这种停滞不仅仅是一个数值上的怪癖；它是矩阵特性的深刻反映。对于性态良好的对称矩阵，它们对向量的影响完全由其特征值描述。但对于在对流主导问题中常见的“非正规”矩阵，特征值并不能说明全部情况。这些矩阵可能导致奇怪的瞬态增长，在最终衰减某些向量之前先将其放大。GMRES，在其重启动变体 GMRES($m$) 的短时记忆下，可能只看到最初的放大而迷失方向，无法找到通往解的长期路径。问题的真实景观不是由特征值的精确岛屿描述的，而是由“伪谱”的幽灵大陆描述的。要克服这些问题，我们需要的不仅仅是基础算法。

PDE 的世界不仅限于流动。想象一下敲击鼓面或模拟雷达波。这些是波传播问题，由像 Helmholtz 方程这样的方程描述。当被离散化时，这些问题产生的矩阵不仅是非对称的（或者更普遍地，对于复数是非厄米 (non-Hermitian) 的），而且是“不定的”。一个不定矩阵同时拥有正负特征值，这意味着它不像一张拉伸的橡胶薄膜那样对应一个简单的能量最小化问题。解是相互竞争影响的微妙平衡。在这里，像共轭梯度法 (CG) 这样要求正定矩阵的方法是无用的。即使是能处理不定对称系统的 MINRES 也无能为力。GMRES，以其最低的要求——它只要求一种与矩阵相乘的方法——是少数足够稳健的通用工具之一，能够胜任这项工作。

一些物理系统涉及多个耦合现象。像蜂蜜这样的粘性流体的流动由 Stokes 方程控制，它将流体的速度与压力联系起来。离散化后会产生“鞍点”系统。由此产生的矩阵具有特定的分块结构，并且本质上是不定的。再一次，GMRES 是解决这些耦合系统的自然选择，通常与尊重问题潜在物理结构的巧妙的“分块感知”预处理器配对使用。

我们已经看到，GMRES 在处理困难矩阵时可能会遇到困难。这是否意味着我们必须放弃？完全不是！与其直接解决困难的问题 $A \boldsymbol{x} = \boldsymbol{b}$，我们可以解决一个更容易的相关问题。这就是预处理的艺术。一个好的预处理器 $M$ 是一个在某种程度上近似于 $A$（或其逆）但更容易求逆的矩阵。然后我们让 GMRES 解决 $A M^{-1} \boldsymbol{y} = \boldsymbol{b}$，并以 $\boldsymbol{x} = M^{-1} \boldsymbol{y}$ 的形式恢复我们的解。目标是使预处理后的矩阵 $A M^{-1}$“更好”——更接近单位矩阵，特征值聚集在一起，且非正规性不那么严重。

我们在哪里能找到这样一个神奇的矩阵 $M$？有时，我们通过代数方法构建它。对于一个稀疏矩阵 $A$，我们可以进行“不完全”分解，创建具有与 $A$ 相同稀疏模式的下三角和上三角因子 $L$ 和 $U$。这种不完全 LU (ILU) 分解给我们一个预处理器 $M=LU$，它是 $A$ 的一个廉价但有效的近似。

一个更优美的想法是使用一种简单的迭代方法作为一种更强大方法的预处理器。逐次超松弛 (SOR) 方法是一种经典的迭代求解器。单独使用时，它可能收敛缓慢。但是，一次 SOR 扫描可以用作 GMRES 内部的预处理步骤。这就像在使用重型电动工具之前，先用一个简单的手动工具进行快速调整。SOR 方法的松弛因子 $\omega$ 成为预处理器的调节旋钮，使我们能够优化外部 GMRES 算法的性能。这种分层使用方法是科学计算中一个强大且反复出现的主题。克服非正规性和停滞挑战的最佳策略，几乎总是涉及像 GMRES 这样的稳健 Krylov 求解器和一个精心选择的预处理器的组合。

当我们看到 GMRES 超越其最初的背景时，它的真正力量和美感变得显而易见。它不仅仅是用来求解 $\boldsymbol{x}$ 是一个数字列向量的 $A \boldsymbol{x} = \boldsymbol{b}$ 问题。

首先，Krylov 子空间中的“向量”可以是任何属于向量空间的东西。例如，在控制理论中，人们会遇到 Sylvester 方程 $AXB + CXD = E$，其中未知数 $X$ 本身就是一个矩阵。我们可以将 GMRES 应用于这个问题，此时它操作的“向量”是矩阵。在一个子空间上最小化残差的抽象逻辑保持不变，这展示了该方法令人难以置信的通用性。

也许最深刻的应用是当 GMRES 成为一个更宏大方案内部的一个组件时。科学中的许多（如果不是大多数）问题都是非线性的。解决非线性系统 $F(\boldsymbol{x}) = \boldsymbol{0}$ 的一种常用方法是牛顿法。在每一步，牛顿法都用一个线性问题来近似原问题：$J(\boldsymbol{x}_k) \delta \boldsymbol{x}_k = -F(\boldsymbol{x}_k)$，其中 $J$ 是一阶导数的雅可比矩阵。对于大规模问题，构建和存储雅可比矩阵 $J$ 的成本高得令人望而却步，甚至是不可能的。

这就是“无雅可比的牛顿-Krylov” (JFNK) 方法的魔力所在。我们使用 GMRES 在每一步求解线性牛顿系统。但 GMRES 不需要看到矩阵 $J$；它只需要知道如何计算任何向量 $\boldsymbol{v}$ 的乘积 $J\boldsymbol{v}$。而这个矩阵-向量乘积可以用有限差分来近似：

这太惊人了。我们可以利用雅可比矩阵的力量而无需计算它！我们只需要评估我们知道如何计算的原始非线性函数 $F$。GMRES 在牛顿法的非线性底盘内充当强大的线性引擎，探测系统的响应以找到通往解的方向。JFNK 方法是当今科学界一些最大、最复杂模拟背后的主力军。

最后，我们在意想不到的地方发现了 GMRES 的逻辑，揭示了计算科学中深刻的统一性。在量子化学中，自洽场 (SCF) 方法是一个寻找分子电子结构的迭代过程。一种加速其收敛的流行方法叫做 DIIS (迭代子空间直接求逆)。乍一看，DIIS 与 GMRES 完全不同。然而，如果将 DIIS 应用于一个线性问题，结果表明它在数学上等价于 GMRES。两种方法都从先验信息的子空间中构建解，并且都通过最小化残差来实现。它们是同一种残差最小化基本语言的两种不同方言。

从河流的流动和鼓的振动，到控制系统的设计和分子轨道的计算，GMRES 的印记无处不在。它不仅仅是一个算法；它是解决问题的基本原则——证明了在一个有限的、巧妙选择的可能性空间内寻找最佳答案的力量。