GMRES算法数学原理及推导过程广义最小残差法（Generalized Minimal Residual, GMRES）

背景介绍

广义最小残差法（Generalized Minimal Residual, GMRES）广泛应用于计算流体力学和电磁仿真等领域，是科学计算中求解大型稀疏线性方程组 $A\mathbf{x}=b$ 的主流迭代算法之一，适用的问题比较广，其中A为稀疏矩阵，b通常称为右端向量。与传统的直接解法相比，GMRES通过Krylov子空间投影逐步逼近方程的解，从而减少计算量并提高数值稳定性。

传统的解法一般就是求 $A$ 的逆矩阵 $A^{-1}$ ， $\mathbf{x} = A^{-1}b$

基本数学知识

Cayley-Hamilton Theorem

对于 $n \times n$ 矩阵 $A$ ，特征值 $\lambda$ 和特征向量 $\mathbf{x}$ 满足：

A\mathbf{x} = \lambda \mathbf{x} \;\; (\mathbf{x} \neq \mathbf{0})

移项可得：

(\lambda I_n - A)\mathbf{x}=\mathbf{0}

其中 $I_n$ 是 $n \times n$ 单位矩阵，可以理解成向量中的1，也是标量 $\lambda$ 变成向量参与运算的媒介。由于 $\mathbf{x}$ 是非零向量，齐次线性方程组 $(1)$ 有非零解的充要条件是系数矩阵的行列式为0：

det(\lambda I_n - A) = 0

即特征方程，展开左边：

\lambda I_n - A = \begin{pmatrix} \lambda - a_{11} & -a_{12} & \cdots & -a_{1n} \\ -a_{21} & \lambda - a_{22} & \cdots & -a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ -a_{n1} & -a_{n2} & \cdots & \lambda - a_{nn} \end{pmatrix}

根据Leibniz公式， $\lambda$ 的最高次项为 $n$ 即主对角线相乘 $(\lambda - a_{11})(\lambda - a_{22})\cdots(\lambda - a_{nn})$ ，且系数为1。这个式子中每包含一个非对角线元素，则必然会替换掉一个主对角线的 $\lambda$ 项，从而得到 $\lambda^{n-1},\lambda^{n-2},\cdots,\lambda^{0}$ 。当 $\lambda=0$ 时， $det(0 \cdot I_n - A) = (-1)^ndet(A)$ ，因此常数项 $c_0 = (-1)^ndet(A)$ 。

最终得到特征多项式方程，即：

p(\lambda) = \lambda^n + c_{n-1}\lambda^{n-1} + c_{n-2}\lambda^{n-2} + \cdots + c_1\lambda + c_0 = 0

而Cayley-Hamilton Theorem的核心结论就是：矩阵 A 满足其自身的特征多项式方程：

p(A) = A^n + c_{n-1}A^{n-1} + c_{n-2}A^{n-2} + \cdots + c_1A + c_0 = 0 \;\;\;\; (1)

Krylov Subspace

对式子(1)两边同时乘以一个 $A^{-1}$ :

A^{n-1} + c_{n-1}A^{n-2} + c_{n-2}A^{n-3} + \cdots + c_1I + c_0A^{-1} = 0

再同除一个 $c_0$ 并移项：

A^{-1} = \frac {-1}{c_0} A^{n-1} + \frac {-c_{n-1}}{c_0}A^{n-2} + \cdots + \frac {-c_1}{c_0}I

于是我们可以得到：

A^{-1} \in span\{A^{n-1},A^{n-2},\ldots,I\}

即 $A^{-1}$ 在这些向量张成的子空间里，但这并不是最终的Krylov Subspace，我们还要新设一个向量 $\mathbf{r}$ ，然后对上面的式子两边同乘 $\mathbf{r}$ ，得到：

A^{-1}\mathbf{r} \in span\{\mathbf{r},A^{1}\mathbf{r},A^{2}\mathbf{r},\ldots,A^{n-1}\mathbf{r}\} \;\;\;\; (2)

这才是Krylov Subspace，至于为什么要乘以一个 $mathsub{r}$ ，后面将会解释。

Gram-Schmidt and Arnoldi Iteration

一般解大型稀疏矩阵 $A$ ，维度 $n$ 是非常大的，虽然我们已经得到了Krylov Subspace，但是我们无法保证张成这个子空间的基向量都是线性独立的，以至于如果直接在这个子空间里搜寻近似解 $\mathbf{x}$ ，效率其实是很低的。为此，我们要想办法用更少的基向量来张成这个子空间。于是我们很容易想到，如果这个子空间里的基向量都相互正交并且范数都为1的话，所需要的基向量数肯定是最少的，也就大大提高了我们搜寻正确答案的效率。

我们从第一个元素 $\mathbf{r}$ 开始，易得：

q_0 = \frac{\mathbf{r}}{\mathbf{||r||}}

再来看第二个元素 $A^1\mathbf{r}$ ，很简单，我们只需要找到它垂直于 $q_0$ 的分量，具体如图：

我们设 $k_1=A\mathbf{r}$ ， $u_1$ 为 $k_1$ 垂直于 $q_0$ 的分量， $proj\_q_0(k_1)$ 为 $k_1$ 在 $q_0$ 上的投影，易得：

u_1 = k_1 - proj\_q_0(k_1) \\ q_1 = \frac{u_1}{||u_1||}

以此类推：

u_2 = k_2 - proj\_q_0(k2) - proj\_q_1(k2), \; q_2 = \frac{u2}{||u2||} \\ u_3 = k_3 - proj\_q_0(k3) - proj\_q_1(k3) - proj\_q_2(k3), \; q_3 = \frac{u3}{||u3||} \\ \cdots \\ u_m = k_m - \sum_{i=0}^{m-1} proj\_q_i(k_{m}), \; q_m = \frac{u_m}{||u_m||} \;\;\;\; (3)

计算 $u$ 的过程被称为Gram-Schmidt Orthogonalization(正交化)，计算 $q$ 的过程被称为Gram-Schmidt Orthonomalization(正交归一化)，每一个被计算出来的 $q$ 又被称为Arnoldi Vector，整个过程我们称为Arnoldi Iteration

需要注意的是：这里的 $m$ 一定是小于 $n$ 的，除非最开始的Krylov Subspace里的向量就全部互相正交，但这显然不太可能，当我们计算出的 $q_i = 0$ 的时候，我们就可以认为 $Q$ 里已经包含了所有的维度，不再需要迭代

GMRES算法过程

回到 $A\mathbf{x} = \mathbf{b}$ 问题本身，大多数时候，我们并不关注正确答案本身，而是假设一个答案 $\mathbf{x_0}$ ，当这个答案与正确答案 $\mathbf{x}$ 之间的差值 $\mathbf{r}$ 足够小的时候，我们就认为我们找到了正确答案。

设初始答案为 $\mathbf{x_0}$ ，则有：

\begin{align*} &\mathbf{r_0} = A\mathbf{x_0} - \mathbf{b} \\ &A^{-1}\mathbf{r_0} = \mathbf{x_0} - A^{-1}\mathbf{b} \\ &A^{-1}\mathbf{r_0} = \mathbf{x_0} - \mathbf{x} \\ &\mathbf{x} = \mathbf{x_0} - A^{-1}\mathbf{r_0} \;\;\;\; \end{align*}

看到这里，我相信你已经知道为什么前面的Krylov Subspace，即 $A^{-1}\mathbf{r} \in span\{\mathbf{r},A^{1}\mathbf{r},A^{2}\mathbf{r},\ldots,A^{n-1}\mathbf{r}\}$ 需要这样构造了吧。再经过Arnoldi Iteration，我们就能得到维度更小，但是空间一致的子空间 $span\{q_0,q_1,q_2,\cdots,q_m\}$ 。可以表示为 $Q_m = [q_0 \; q_1 \; q_2 \; \cdots \; q_m]$ 。

于是我们可以得到：

\mathbf{x} = \mathbf{x_0} - Q_my \;\;\;\;

为了使等式成立，我们引进了新的变量 $y$ ，这个 $y$ 可以理解成 $x$ 在各个维度上的分量。

所以:

\begin{align*} &\mathbf{r} = A\mathbf{x} - \mathbf{b} \\ &\mathbf{r} = A(\mathbf{x_0} - Q_iy) - \mathbf{b} \\ &\mathbf{r} = \mathbf{r_0} - AQ_iy \;\;\;\; (4) \end{align*}

这里的 $\mathbf{x}$ 指的不是正确答案，而是在误差范围内我们要寻找的答案， $i$ 指的是迭代次数，即：

\begin{aligned} &for \;\; i \;\; in \;\; range(m): \\ &\mathbf{r} = \mathbf{r_0} - AQiy \end{aligned}

这就是GMRES算法的迭代过程，当 $r < resid\_limit$ ，这里的 $resid\_limit$ 是我们规定的误差范围，一般为 $1e^{-6}$ ，迭代终止，也就认为我们找到了最终答案。为了方便我们的迭代，我们引入了一个新的矩阵，Hessenberg Matrix(海森堡矩阵)，它的格式固定为：

H_{i+1,i} = \begin{pmatrix} h_{0,0} & h_{0,1} & \cdots & h_{0,i} \\ h_{1,0} & h_{1,1} & \cdots & h_{1,i} \\ 0 & h_{2,1} & \cdots & h_{2,i} \\ 0 & 0 & \cdots & h_{3,i} \\ \cdots & \cdots & \cdots & \cdots \\ 0 & 0 & \cdots & h_{i+1,i} \\ \end{pmatrix}

我们将这个矩阵记为 $\widetilde{H}_i$ ，他是用来搭建 $AQ_i$ 和 $Q_{i+1}$ 之间关系的桥梁，即： $AQ_i = Q_{i+1} \widetilde{H}$

回忆Arnoldi Iteration**(注意：以下的m指的不是Krylov Subspace维度数，而是迭代次数)**：

u_m = k_m - \sum_{i=0}^{m-1} proj\_q_i(k_{m}), \; q_m = \frac{u_m}{||u_m||}

因为 $proj\_q_i(k_m)$ 是 $k_m$ 在 $q_i$ 上的投影，也就可以理解成一个常数乘以 $q_i$ 。很明显， $q_m$ 和 $u_m$ 也是这样的关系，所以：

u_m = a_mq_m = k_m - \sum_{i=0}^{m-1}b_iq_i \\ q_m = c_0k_m + c_1q_0 + c_2q_1 + \cdots + c_mq_{m-1}

也就是说 $q_m \in span\{k_m,\;q_0,\;q_1,\;\cdots\;,\;q_{m-1}\}$ ，而我们又很容易知道 $span\{q_0,\;q_1,\;\cdots\;,\;q_{m-1}\} = span\{k_0,\;k_1,\;\cdots\;,\;k_{m-1}\}$ ，因为前者只是用互相正交且范数为1的基重新去表示了Krylov Subspace，张成的空间肯定是一致的。所以：

q_m \in span\{k_0 \cdots k_m\} \\ Aq_m \in span\{Ak_0 \cdots Ak_m\}

又因为 $k_m = A^mr$ ，所以：

Aq_m \in span\{k_1 \cdots k_{m+1}\}

这里的 $Aq_m$ 是新加入的向量，于是我们就会发现，与其按照式子(2)的方式来得到新的维度 $k_{m+1}$ ，不如直接计算 $Aq_m$ 来得到 $k_{m+1}$ ，然后就能得到新的 $q_{m+1}$ 维度，加入 $Q_m$ 变成 $Q_{m+1}$ 进入下一个迭代，我们记通过这种方式得到的子空间 $span\{k_1 \cdots k_{m+1}\}$ 为 $k_{m+1}^{’}$ 。

所以我们可以展开 $AQ_m$ :

AQ_m = [Aq_0,\cdots,Aq_m] = [k_1^{’},\cdots,k_{m+1}^{’}] = K_{m+1}^{’}

又因为： $k_{m+1}^{’} \in span\{q_0 \cdots q_{m+1}\}$ ，也就是 $Q_{m+1}$ 的一个线性组合，所以：

\begin{align*} K_{m+1}^{’} &= \begin{bmatrix} h_0q_0,\;h_1q_1,\;\cdots ,\;h_{m+1}q_{m+1} \end{bmatrix} \\ &= \begin{bmatrix} q_0,\;q_1,\;\cdots,\;q_{m+1} \end{bmatrix} \begin{bmatrix} h_{0,0} & h_{0,1} & \cdots & h_{0,m} \\ h_{1,0} & h_{1,1} & \cdots & h_{1,m} \\ 0 & h_{2,1} & \cdots & h_{2,m} \\ 0 & 0 & \cdots & h_{3,m} \\ \cdots & \cdots & \cdots & \cdots \\ 0 & 0 & \cdots & h_{m+1,m} \\ \end{bmatrix} \\ &=Q_{m+1}\widetilde{H} \end{align*}

相当于我们的 $\widetilde{H}$ 矩阵里的每一列 $h_i$ 实际上记录着的就是如何将 $q_0 \cdots q_{i+1}$ 转变为 $k_{i+1}^{’}$ 的信息。所以，式子(3)可以重写为：

\begin{align*} k_{m+1}^{’}&=q_{m+1} * ||u_{m+1}|| + \sum_{i=0}^{m}proj\_q_i(k_{m+1}^{’})\\ &=q_{m+1} * ||u_{m+1}|| + \sum_{i=0}^{m}(k_{m+1}^{’}·q_i)q_i\\ &=u_{m+1} + \sum_{i=0}^{m}h_iq_i \\ &=Aq_m \end{align*}

然后我们就只需要解出 $u_{m+1}$ ，也就得到了对应的 $q_{m+1}$ 和 $||u_{m+1}||$ 。关系为：

q_{m+1} = \frac{Aq_m - \sum_{i=0}^{m}h_iq_i}{||Aq_m - \sum_{i=0}^{m}h_iq_i||} = \frac{u_{m+1}}{||u_{m+1}||}

而 $||u_{m+1}||$ 恰好就是我们想要的 $h_{m+1,m}$ ，因为我们看 $Q_{m+1}\widetilde{H}$ 的最后一列：

\begin{align*} Aq_m &= h_{m+1,m}q_{m+1} + \sum_{i=0}^{m}h_{i,m}q_i \\ &=||u_{m+1}||q_{m+1} + \sum_{i=0}^{m}h_iq_i \end{align*}

这也就是完整的一次迭代过程，每次迭代实际上都是在计算 $\widetilde{H}$ 的第 $m$ 列。

式子(4)也就能改写为：

\begin{align*} &r = r_0 - Q_{m+1}\widetilde{H}_{m+1,m}y\\ &r = \beta q_0 - Q_{m+1}\widetilde{H}_{m+1,m}y \\ &r = Q_{m+1}(\beta e - \widetilde{H}_{m+1,m}y) \\ &||r||_2 = ||Q_{m+1}(\beta e - \widetilde{H}_{m+1,m}y)||_2 \end{align*}

又因为 $Q_{m+1}$ 的所有矢量都是单位矢量且互相正交，所以：

||r||_2 = ||(\beta e - \widetilde{H}_{m+1,m}y)||_2

到这里，问题也就简化为了求最小二乘解的问题了，一般用QR分解求解。在这里，由于GMRES算法求解的大多数是稀疏矩阵乘，采用Givens变换比较好。