## 线性方程组的预处理技术 ¶

首先看一个例子（线性方程组）：

$$ \underbrace{\begin{bmatrix} 1 & 10^4 \\ 1 & 1 \end{bmatrix}}_{A} \quad \underbrace{\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}}_{x} = \underbrace{\begin{bmatrix} 10^4 \\ 2 \end{bmatrix}}_{b} $$

计算 $A$ 的 条件数（以行和范数） 不知道矩阵的范数和条件数是啥？请看 向量和矩阵的范数 和 矩阵的条件数 ：

$$ A^{-1} = \frac{1}{10^4 - 1} \begin{bmatrix} -1 & 10^4 \\ 1 & -1 \end{bmatrix} $$

$$ \kappa_{\infty}(A) = \lVert A \rVert _{\infty} \lVert A^{-1} \rVert _{\infty} = (10^4 + 1) \times \frac{1}{10^4 - 1} (10^4 + 1) = \frac{(10^4 + 1)^2}{10^4 - 1} \approx 10^4 $$

$A$ 是病态矩阵，主要原因是 $A$ 元素的绝对值相差很大，如果引入因子 $10^4$，除 $A$ 的第一行（同时也要除 $b$ 的第一行），那么得到：

$$ \underbrace{\begin{bmatrix} 10^{-4} & 1 \\ 1 & 1 \end{bmatrix}}_{A^{\prime}} \quad \underbrace{\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}}_{x} = \underbrace{\begin{bmatrix} 1 \\ 2 \end{bmatrix}}_{b^{\prime}} $$

用矩阵形式表示，就是两边同时左乘了一个矩阵 $M$：

$$ \begin{aligned} \underbrace{\begin{bmatrix} 1 & 10^4 \\ 1 & 1 \end{bmatrix}}_{A} \quad \underbrace{\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}}_{x} &= \underbrace{\begin{bmatrix} 10^4 \\ 2 \end{bmatrix}}_{b} \\ \underbrace{\begin{bmatrix} 10^{-4} & 0 \\ 0 & 1 \end{bmatrix}}_{M} \quad \underbrace{\begin{bmatrix} 1 & 10^4 \\ 1 & 1 \end{bmatrix}}_{A} \quad \underbrace{\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}}_{x} &= \underbrace{\begin{bmatrix} 10^{-4} & 0 \\ 0 & 1 \end{bmatrix}}_{M} \quad \underbrace{\begin{bmatrix} 10^4 \\ 2 \end{bmatrix}}_{b} \\ \underbrace{\begin{bmatrix} 10^{-4} & 1 \\ 1 & 1 \end{bmatrix}}_{A^{\prime}} \quad \underbrace{\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}}_{x} &= \underbrace{\begin{bmatrix} 1 \\ 2 \end{bmatrix}}_{b^{\prime}} \end{aligned} $$

引入 $M$ 的意义何在呢？

此时，尽管在代数上，$Ax = b$ 和 $A^\prime x = b^\prime$ 的解是一一对应的， 解这两个线性方程组是等价的。 但是，在数值上，$A^\prime = MA$ 却比 $A$ 性质更好，因为它的条件数远低于 $A$：

$$ (A^\prime)^{-1} = \frac{1}{10^{-4} - 1} \begin{bmatrix} 1 & -1 \\ -1 & 10^{-4} \end{bmatrix} $$

$$ \kappa_{\infty}(A^\prime) = (1 + 1) \times \frac{1}{1 - 10^{-4}} (1 + 1) = \frac{4}{1 - 10^{-4}} \approx 4 $$

可以看出：

$$ \begin{aligned} \kappa_{\infty}(A^\prime) = \kappa_{\infty}(MA) &\approx 4 \\ \kappa_{\infty}(A) &\approx 10^4 \\ \kappa_{\infty}(MA) &\ll \kappa_{\infty}(A) \end{aligned} $$

矩阵 $M$ 的引入，可以将方程组 $Ax = b$ 转换成等价的 $MAx = Mb$， 在保证解的正确性的前提下，改善方程系数矩阵的条件数，让寻找数值解的过程更稳定。

倒霉的是，大规模线性方程组大多是病态的，会给数值方法带来困难：

1. 使迭代收敛缓慢。
2. 对数值解的精度有很大影响。

对于第一点，当前最有效的处理手段是预处理， 预处理用来改善系数矩阵的条件数，上面引入 $M$ 的例子就是在做预处理。

令 $P = M^{-1}$，矩阵 $P$ 一般被称为预处理子。从另一个角度说，预处理子 $P$ 应当尽可能地 “逼近”/“近似”/“像” $A$。 如果 $P \approx A$，那么 $P^{-1} A \approx I$，单位矩阵的条件数是 1，预处理后的迭代只需要一步就能解出来！这是最完美的情况。

线性方程组的预处理

对于线性方程组 $$ Ax = b $$ 其中 $A \in \mathbb{R}^{n \times n}$ 非奇异，$b \in \mathbb{R}^n$， 我们可以对其进行 预处理 (Preconditioning) 操作：

• 左预处理 (Left Preconditioning)：$A$ 左乘矩阵 $P^{-1}$，同时右边 $b$ 也左乘 $P^{-1}$，这个方式最常用。 $$ \color{blue} \boxed{\begin{aligned} \text{左预处理} &\quad P^{-1} A x = P^{-1} b \\[8pt] \text{预处理子} &\quad P \end{aligned}} $$ $P$ 称为 预处理子 (Preconditioner)，使得条件数 $\kappa(P^{-1} A) \ll \kappa(A)$

• 右预处理 (Right Preconditioning)：$A$ 右乘矩阵 $P^{-1}$，为了保证同解，$x$ 需要左乘 $P$ $$ \begin{array}{c} A I x = b \\ A (P^{-1} P) x = b \\ (A P^{-1}) (P x) = b \end{array} $$

  那么先解 $A P^{-1} y = b$ 得到 $y$，再解 $P x = y$ 得到 $x$

  $$ \color{blue} \boxed{\begin{aligned} &\text{右预处理} \quad \begin{aligned} A P^{-1} y &= b \\ x &= P^{-1} y \end{aligned} \\[8pt] &\text{预处理子} \quad P \end{aligned}} $$

• 两边预处理 (Two-sided / Split Preconditioning)：$A$ 左乘 $L^{-1}$ 并右乘 $R^{-1}$，是上面两种方式的结合。 $$ \begin{array}{c} L^{-1} A x = L^{-1} b \\ L^{-1} A (R^{-1} R) x = L^{-1} b \\ (L^{-1} A R^{-1}) (R x) = L^{-1} b \end{array} $$

  先解 $L^{-1} A R^{-1} y = L^{-1} b$ 得到 $y$，再解 $R x = y$ 得到 $x$

  $$ \color{blue} \boxed{\begin{aligned} &\text{两边预处理} \quad \begin{aligned} L^{-1} A R^{-1} y &= L^{-1} b \\ x &= R^{-1} y \end{aligned} \\[8pt] &\text{预处理子} \qquad P = LR \end{aligned}} $$

  这里预处理子 $P$ 是通过某种矩阵分解方式，分解成 $LR$ 的形式。

经过预处理的方程组称为 带预处理的线性方程组 (preconditioned linear systems)， 相应的有 带预处理的求解器 (preconditioned iterative solvers)。

关于预处理子

预处理子不是随便搞一个都可以的，它的构造和选取是有讲究的：

1. $P \approx A$，即 $P^{-1} \approx A^{-1}$，这点和 2 在本质上是一致的
2. $P^{-1}A$ 或 $AP^{-1}$ 的条件数一定（远）低于 $A$（或者特征值分布更好），否则失去意义
3. $Pu = v$ 容易求解，或者 $P$ 易于求逆，否则 $P$ 的使用成本高

预处理本质上是构造一个 $A$ 的近似 $P \approx A$，同时易于求解（$Pu = v$）

预处理技术特别对 Krylov 子空间迭代的收敛性质有明显改善。 通常用 Krylov 子空间法求解预处理后的方程组时，每步迭代都需要额外求解一个以 $P$ 为系数矩阵的方程组。

## 带预处理的 GMRES ¶

### 左预处理的 GMRES ¶

如果 GMRES 方法应用左预处理，预处理子为 $P$ $$ P^{-1} A x = P^{-1} b $$

第 $j$ 步迭代的残量为 $$ \begin{aligned} \tilde{r}_j &= P^{-1} (b - A x_j) \\ &= P^{-1} r_j \end{aligned}$$

Arnoldi 迭代构造的是左预处理的 Krylov 子空间 $\mathcal{K}_j(P^{-1} A, r_0)$ 的一组正交基 $$ \text{span} \lbrace \tilde{r}_0, P^{-1} A \tilde{r}_0, \dots, (P^{-1} A)^{j-1} \tilde{r}_0 \rbrace $$

下面是左预处理的 GMRES，与原始 GMRES 算法区别只在第 7 行和 第 12-13 行 预处理子 $\boldsymbol{P}$ 的逆通常不是显式计算得到，而是要构造 $\boldsymbol{P}$，使得线性系统 $\boldsymbol{P} \boldsymbol{u} = \boldsymbol{v}$ 易于求解、代价小。可以看到，算法第 13 行需要解一个以预处理子 $\boldsymbol{P}$ 为系数矩阵的线性方程组。 ：初始残量的计算，和每次迭代构造子空间 $\mathcal{K}_j(P^{-1} A, r_0)$ 新的基向量 $\boldsymbol{v}_j$ 时候。

注意 $\rho_{j+1}$ 不是原始（真实）残量的范数，如果非要用到真实残量，比如判定是否终止迭代，没什么好办法，只能计算 $\boldsymbol{r}_j = \boldsymbol{P} \boldsymbol{z}_j$。

    \begin{algorithm}
\caption{Left Preconditioned GMRES}
\begin{algorithmic}
\PROCEDURE{givens}{$x, y$}
    \STATE $r := \sqrt{x^2 + y^2}$
    \STATE $\displaystyle c, s := x / r, y / r$
    \RETURN $c, s$
\ENDPROCEDURE
\PROCEDURE{left-preconidtioend-gmres}{$\boldsymbol{P}$, $\boldsymbol{A}$, $\boldsymbol{b}$, $\boldsymbol{x}_0$, $\epsilon$}
    \STATE Solve the linear system $\boldsymbol{P} \boldsymbol{r}_0 := \boldsymbol{b} - \boldsymbol{A} \boldsymbol{x}_0$ for $\boldsymbol{r}_0$
    \STATE $\beta := \lVert \boldsymbol{r}_0 \rVert _2$
    \STATE $\boldsymbol{v}_1 := \boldsymbol{r}_0 \ / \ \beta$
    \STATE $\rho_1 := \beta$
    \FOR{$j := 1$ \TO $m$}
        \STATE $\boldsymbol{u}_j := \boldsymbol{A} \boldsymbol{v}_j$
        \STATE Solve the linear system $\boldsymbol{P} \boldsymbol{w}_j = \boldsymbol{u}_j$ for $\boldsymbol{w}_j$
        \FOR{$i := 1$ \TO $j$}
            \STATE $h_{ij} := \langle \boldsymbol{w}_j, \boldsymbol{v}_i \rangle$
            \STATE $\boldsymbol{w}_j := \boldsymbol{w}_j - h_{ij} \boldsymbol{v}_i$
        \ENDFOR
        \STATE $h_{j+1,j} := \lVert \boldsymbol{w}_j \rVert_2$
        \FOR{$i := 1$ \TO $j-1$}
            \STATE $\begin{bmatrix}h_{i,j} \\ h_{i+1,j}\end{bmatrix} = \begin{bmatrix}c_i & s_i \\ -s_i & c_i\end{bmatrix} \begin{bmatrix}h_{i,j} \\ h_{i+1,j}\end{bmatrix}$
        \ENDFOR
        \IF{$h_{j+1,j} = 0$}
            \STATE $m = j$
            \BREAK
        \ENDIF
        \STATE $\boldsymbol{v}_{j+1} := \boldsymbol{w}_j \ / \ h_{j+1,j}$
        \STATE $c_j, s_j :=$ \CALL{givens}{$h_{j,j}, h_{j+1,j}$}
        \STATE $h_{j,j} = c_j h_{j,j} + s_j h_{j+1,j} \enspace , \enspace h_{j+1,j} = 0$
        \STATE $\rho_{j+1} = -s_j \rho_j \enspace , \enspace \rho_j = c_j \rho_j$
        \IF{$\rho_{j+1} \ / \ \beta < \epsilon$}
            \STATE $m = j$
            \BREAK
        \ENDIF
    \ENDFOR
    \STATE Let $\boldsymbol{H}_m := \begin{bmatrix} h_{ij} \end{bmatrix}_{\ i, j \in [1, m]}$, $\boldsymbol{\rho}_m := \begin{bmatrix}\rho_j\end{bmatrix}_{j \in [1, m]}$ and $\boldsymbol{V}_m := \begin{bmatrix}\boldsymbol{v}_j\end{bmatrix}_{j \in [1, m]}$
    \STATE Solve the upper triangular linear system $\boldsymbol{H}_m \boldsymbol{y}_m = \boldsymbol{\rho}_m$ for $\boldsymbol{y}_m$
    \STATE Compute the solution $\boldsymbol{x}_m := \boldsymbol{x}_0 + \boldsymbol{V}_m \boldsymbol{y}_m$
\ENDPROCEDURE
\end{algorithmic}
\end{algorithm}

### 右预处理的 GMRES ¶

右预处理线性方程组 $Ax = b$，预处理子为 $P$ $$ \begin{aligned} A P^{-1} z &= b \\ P x &= z \end{aligned} $$

我们主要看第一个方程，因为我们求解 $x$ 必须要先得到 $z$，它在第 $j$ 步迭代的残量为

$$ \begin{aligned} \tilde{r}_k &= b - A P^{-1} z_j \\ &= b - A x_j \qquad (x = P^{-1} z)\\ &= r_j \end{aligned} $$

什么情况？它的残量和原始方程 $Ax = b$ 是一样的！所以我们不需要对它应用右预处理！也不用担心终止条件的判定问题，初始残量直接计算 $b - A x_0$ 拿来用就行，完全不用操心了喵，好幸福！

但由于系数矩阵变成了 $A P^{-1}$，那么 Arnoldi 迭代构造的是右处理的 Krylov 子空间 $\mathcal{K}_j(A P^{-1}, r_0)$ 的一组正交基

$$ \text{span} \lbrace r_0, A P^{-1} r_0, \dots, (A P^{-1})^{j-1} r_0 \rbrace $$

根据原始 GMRES，我们在最后得到的是 $z_m$

$$ z_m = z_0 + V_m y_m $$

其中 $z_0 = P x_0$，$z_m$ 是方程组 $A P^{-1} z_m = b$ 的解，这时再解方程组

$$ P x_m = z_m $$

就可以得到原始方程组 $Ax = b$ 的解 $x_m$

下面是右预处理的 GMRES，和原始 GMRES 的区别在 11-12 行和 31-32 行 同上，这里算法第 11 行和第 32 行需要解一个以预处理子 $\boldsymbol{P}$为系数矩阵的线性方程组。 ：每次迭代 Arnoldi 构造 Krylov 子空间新的基向量 $\boldsymbol{v}_j$ 的时候需要应用右预处理， 以及在最后需要多解一个方程组 $\boldsymbol{P} \boldsymbol{x}_m = \boldsymbol{z}_m$

    \begin{algorithm}
\caption{Right Preconditioned GMRES}
\begin{algorithmic}
\PROCEDURE{givens}{$x, y$}
    \STATE $r := \sqrt{x^2 + y^2}$
    \STATE $\displaystyle c, s := x / r, y / r$
    \RETURN $c, s$
\ENDPROCEDURE
\PROCEDURE{right-preconidtioend-gmres}{$\boldsymbol{P}$, $\boldsymbol{A}$, $\boldsymbol{b}$, $\boldsymbol{x}_0$, $\epsilon$}
    \STATE $\boldsymbol{r}_0 := \boldsymbol{b} - \boldsymbol{A} \boldsymbol{x}_0$
    \STATE $\beta := \lVert \boldsymbol{r}_0 \rVert _2$
    \STATE $\boldsymbol{v}_1 := \boldsymbol{r}_0 \ / \ \beta$
    \STATE $\rho_1 := \beta$
    \FOR{$j := 1$ \TO $m$}
        \STATE Solve the linear system $\boldsymbol{P} \boldsymbol{u}_j = \boldsymbol{v}_j$ for $\boldsymbol{u}_j$
        \STATE $\boldsymbol{w}_j := \boldsymbol{A} \boldsymbol{u}_j$
        \FOR{$i := 1$ \TO $j$}
            \STATE $h_{ij} := \langle \boldsymbol{w}_j, \boldsymbol{v}_i \rangle$
            \STATE $\boldsymbol{w}_j := \boldsymbol{w}_j - h_{ij} \boldsymbol{v}_i$
        \ENDFOR
        \STATE $h_{j+1,j} := \lVert \boldsymbol{w}_j \rVert _2$
        \FOR{$i := 1$ \TO $j-1$}
            \STATE $\begin{bmatrix}h_{i,j} \\ h_{i+1,j}\end{bmatrix} = \begin{bmatrix}c_i & s_i \\ -s_i & c_i\end{bmatrix} \begin{bmatrix}h_{i,j} \\ h_{i+1,j}\end{bmatrix}$
        \ENDFOR
        \IF{$h_{j+1,j} = 0$}
            \STATE $m = j$
            \BREAK
        \ENDIF
        \STATE $\boldsymbol{v}_{j+1} := \boldsymbol{w}_j \ / \ h_{j+1,j}$
        \STATE $c_j, s_j :=$ \CALL{givens}{$h_{j,j}, h_{j+1,j}$}
        \STATE $h_{j,j} = c_j h_{j,j} + s_j h_{j+1,j} \enspace , \enspace h_{j+1,j} = 0$
        \STATE $\rho_{j+1} = -s_j \rho_j \enspace , \enspace \rho_j = c_j \rho_j$
        \IF{$\rho_{j+1} \ / \ \beta < \epsilon$}
            \STATE $m = j$
            \BREAK
        \ENDIF
    \ENDFOR
    \STATE Let $\boldsymbol{H}_m := \begin{bmatrix} h_{ij} \end{bmatrix}_{\ i, j \in [1, m]}$, $\boldsymbol{\rho}_m := \begin{bmatrix}\rho_j\end{bmatrix}_{j \in [1, m]}$ and $\boldsymbol{V}_m := \begin{bmatrix}\boldsymbol{v}_j\end{bmatrix}_{j \in [1, m]}$
    \STATE Solve the upper triangular linear system $\boldsymbol{H}_m \boldsymbol{y}_m = \boldsymbol{\rho}_m$ for $\boldsymbol{y}_m$
    \STATE Compute the solution $\boldsymbol{z}_m := \boldsymbol{z}_0 + \boldsymbol{V}_m \boldsymbol{y}_m$ for $\boldsymbol{z}_m$,  where $\boldsymbol{z}_0 := \boldsymbol{P} \boldsymbol{x}_0$
    \STATE Solve the linear system $\boldsymbol{P} \boldsymbol{x}_m = \boldsymbol{z}_m$ for the final solution $\boldsymbol{x}_m$
\ENDPROCEDURE
\end{algorithmic}
\end{algorithm}

## Flexible GMRES (FGMRES) ¶

前面讲的带预处理的 GMRES 默认预处理子 $P$ 是固定的，迭代过程中不会变。 但是，有些情况下，我们找不到合适的 $P$。

在右处理的 GMRES 中，每步迭代中的 $P^{-1}v_j$ 可能来自于未指定的计算 ── 甚至是其他的迭代过程！ 一般我们将这种迭代过程称为 "Flexible" 迭代。 假设预处理子每步迭代时都可能改变 $P_j$，那么 $u_j$ 变为

$$ \begin{aligned} &\text{Solve the linear system} \ \boldsymbol{P}_j \boldsymbol{u}_j = \boldsymbol{v}_j \ \text{for} \ \boldsymbol{u}_j \\ \rightarrow\quad &\text{Solve the linear system} \ \boldsymbol{P}_j \boldsymbol{u}_j = \boldsymbol{v}_j \ \text{for} \ \boldsymbol{u}_j \end{aligned} $$

最后计算近似解 $x_m$ 变为

$$ \begin{aligned} &\boldsymbol{x}_m := \boldsymbol{x}_0 + \boldsymbol{V}_m \boldsymbol{y}_m \\ \rightarrow\quad &\boldsymbol{x}_m := \boldsymbol{x}_0 + \boldsymbol{U}_m \boldsymbol{y}_m \end{aligned} $$

其中 $\boldsymbol{U}_m := [ \boldsymbol{u}_j ] _{j \in [1, m]}$，$\boldsymbol{y}_m$ 的计算不变。

FGMRES 就是右处理 GMRES 的变体，除了 $\boldsymbol{v}_j$，FGMRES 还需要额外记录迭代中的 $\boldsymbol{u}_j$，即预处理后的 $\boldsymbol{u}_j = \boldsymbol{P}^{-1}_j \boldsymbol{v}_j$，所以 FGMRES 相比 GMRES 内存占用将近 double

    \begin{algorithm}
\caption{Flexible GMRES (FGMRES)}
\begin{algorithmic}
\PROCEDURE{givens}{$x, y$}
    \STATE $r := \sqrt{x^2 + y^2}$
    \STATE $\displaystyle c, s := x / r, y / r$
    \RETURN $c, s$
\ENDPROCEDURE
\PROCEDURE{fgmres}{$\boldsymbol{P}$, $\boldsymbol{A}$, $\boldsymbol{b}$, $\boldsymbol{x}_0$, $\epsilon$}
    \STATE $\boldsymbol{r}_0 := \boldsymbol{b} - \boldsymbol{A} \boldsymbol{x}_0$
    \STATE $\beta := \lVert \boldsymbol{r}_0 \rVert _2$
    \STATE $\boldsymbol{v}_1 := \boldsymbol{r}_0 \ / \ \beta$
    \STATE $\rho_1 := \beta$
    \FOR{$j := 1$ \TO $m$}
        \STATE Solve the linear system $\boldsymbol{P}_j \boldsymbol{u}_j = \boldsymbol{v}_j$ for $\boldsymbol{u}_j$
        \STATE $\boldsymbol{w}_j := \boldsymbol{A} \boldsymbol{u}_j$
        \FOR{$i := 1$ \TO $j$}
            \STATE $h_{ij} := \langle \boldsymbol{w}_j, \boldsymbol{v}_i \rangle$
            \STATE $\boldsymbol{w}_j := \boldsymbol{w}_j - h_{ij} \boldsymbol{v}_i$
        \ENDFOR
        \STATE $h_{j+1,j} := \lVert \boldsymbol{w}_j \rVert _2$
        \FOR{$i := 1$ \TO $j-1$}
            \STATE $\begin{bmatrix}h_{i,j} \\ h_{i+1,j}\end{bmatrix} = \begin{bmatrix}c_i & s_i \\ -s_i & c_i\end{bmatrix} \begin{bmatrix}h_{i,j} \\ h_{i+1,j}\end{bmatrix}$
        \ENDFOR
        \IF{$h_{j+1,j} = 0$}
            \STATE $m = j$
            \BREAK
        \ENDIF
        \STATE $\boldsymbol{v}_{j+1} := \boldsymbol{w}_j \ / \ h_{j+1,j}$
        \STATE $c_j, s_j :=$ \CALL{givens}{$h_{j,j}, h_{j+1,j}$}
        \STATE $h_{j,j} = c_j h_{j,j} + s_j h_{j+1,j} \enspace , \enspace h_{j+1,j} = 0$
        \STATE $\rho_{j+1} = -s_j \rho_j \enspace , \enspace \rho_j = c_j \rho_j$
        \IF{$\rho_{j+1} \ / \ \beta < \epsilon$}
            \STATE $m = j$
            \BREAK
        \ENDIF
    \ENDFOR
    \STATE Let $\boldsymbol{H}_m := \begin{bmatrix} h_{ij} \end{bmatrix}_{\ i, j \in [1, m]}$, $\boldsymbol{\rho}_m := \begin{bmatrix}\rho_j\end{bmatrix}_{j \in [1, m]}$ and $\boldsymbol{U}_m := \begin{bmatrix}\boldsymbol{u}_j\end{bmatrix}_{j \in [1, m]}$
    \STATE Solve the upper triangular linear system $\boldsymbol{H}_m \boldsymbol{y}_m = \boldsymbol{\rho}_m$ for $\boldsymbol{y}_m$
    \STATE Compute the solution by $\boldsymbol{x}_m := \boldsymbol{x}_0 + \boldsymbol{U}_m \boldsymbol{y}_m$
\ENDPROCEDURE
\end{algorithmic}
\end{algorithm}

非常重要的是，第 11 行虽允许预处理子变化，但是每次需要存储预处理后的向量 $\boldsymbol{u}_j$，在最后计算近似解 $\boldsymbol{x}_m$ 时，用到的是 $\boldsymbol{U}_m$（而不是 $\boldsymbol{V}_m$）。

🐮🍺的是，可以把这一步 $\boldsymbol{u}_j = \boldsymbol{P}^{-1}_j \boldsymbol{v}_j$ 当成黑箱，只要弄一个 $\boldsymbol{u}_j$ 出来就行，为啥？因为这个东东是用来扩充 Krylov 子空间的，理论上可以自由选取，只要保证每步的 $\boldsymbol{u}_j$ 和 $\boldsymbol{u}_i (i < j)$ 是线性无关的即可（厚礼谢）！

这使得 FGMRES 有很大的灵活性，甚至可以通过任意的内迭代得到 $\boldsymbol{u}_j$（比如嵌套 GMRES），在这个意义上，迭代器也可以当成预处理子。当然，如果你选取的不合适，算法可能提前终止。

关于 Inner-Outer Iterations

FGMRES 是一种典型的 inner-outer 迭代结构：对于外迭代，当然是用来解原始线性方程组的；对于内迭代，这看起来有点套娃、令人迷惑，用来代替预处理，是什么意思呢？该怎么做呢？是不是 $\boldsymbol{u}_j$ 随便构造一个都可以呢？

FGMRES 第 11-12 行： $$ \begin{array}{l} 11. \quad \text{Solve the linear system} \ \boldsymbol{P}_j \boldsymbol{u}_j = \boldsymbol{v}_j \ \text{for} \ \boldsymbol{u}_j \\ 12. \quad \boldsymbol{w}_j := \boldsymbol{A} \boldsymbol{u}_j \end{array} $$

还记得吗，预处理子 $\boldsymbol{P}$ 是用来 “近似” $\boldsymbol{A}$ 的，第 11 行本质上是在找一个向量 $\boldsymbol{u}_j$，使得 $\boldsymbol{A} \boldsymbol{u}_j \approx \boldsymbol{v}_j$

什么？你是说 $\boldsymbol{u}_j$ 其实是方程组 $\boldsymbol{A} \boldsymbol{u}_j = \boldsymbol{v}_j$ 一个不太 “近似” 的近似解吗？那么，我们为什么不可以用一个迭代器去解 $\boldsymbol{A} \boldsymbol{u}_j = \boldsymbol{v}_j$ 这个线性系统呢？

所以，为了维持收敛性，$\boldsymbol{u}_j$ 应当使 $\boldsymbol{A} \boldsymbol{u}_j \approx \boldsymbol{v}_j$，或者说使残差 $\lVert \boldsymbol{v}_j - \boldsymbol{A} \boldsymbol{u}_j \rVert$ 足够小。 这一步不管用什么方法，无论是先构造一个近似 $\boldsymbol{A}$ 的矩阵 $\boldsymbol{P}$，还是跑几步迭代方法比如 Jacobi、GMRES，甚至是用神经网络猜一个，都是 OK 的。

如果用原始 GMRES 方法作为预处理器，这里第 11 行，以一个内迭代的方式该如何做呢？

GMRES 为内迭代的 FGMRES

GMRES 为内迭代的 FGMRES 第 11-12 行： $$ \begin{array}{l} 11. \quad \mathbf{GMRES}(\boldsymbol{A}=\boldsymbol{A}, \ \boldsymbol{b}=\boldsymbol{v}_j, \ \boldsymbol{x}=\boldsymbol{u}_j, \ m = m_{\text{inner}}) \\ 12. \quad \boldsymbol{w}_j := \boldsymbol{A} \boldsymbol{u}_j \end{array} $$

其中 $\boldsymbol{u}_j$ 可以全 0 初始化，$m_{\text{inner}}$ 是内迭代的次数，一般比较小不会超过 10。

Until next time!

## 参考资料 ¶

• https://en.wikipedia.org/wiki/Preconditioner
• Iterative Methods for Sparse Linear Systems, 2nd Edition, Yousef Saad.
  • Chapter 9.3 Preconditioned GMRES
  • Chapter 9.4.1 Flexible GMRES
• 第五讲：预处理方法，迭代方法与预处理，华东师范大学，潘建瑜
• Thanks to these chatbots:
  • Google Gemini 3 Pro
  • Alibaba Qwen3-Max
  • Kimi K2