Krylov 子空间方法：从多项式逼近到 GMRES

前面几章里，我们一直在讨论如何求解或者近似求解

$$\mathbf{A}\mathbf{x}=\mathbf{b}$$

如果 $\mathbf{A}$ 是一个规模不大的矩阵，那么我们可以考虑直接分解它，比如使用 LU decomposition、QR decomposition，或者在最小二乘问题中使用 SVD 和伪逆

但是在很多数值计算问题里，$\mathbf{A}$ 可能非常大，而且通常是 sparse matrix

这时直接求

$$\mathbf{x}=\mathbf{A}^{-1}\mathbf{b}$$

往往不是一个好主意

第一，显式求逆本身很贵；第二，即使 $\mathbf{A}$ 是稀疏矩阵，$\mathbf{A}^{-1}$ 通常也会变成稠密矩阵；第三，我们很多时候并不真的需要整个 $\mathbf{A}^{-1}$，而只需要它作用在当前右端项 $\mathbf{b}$ 上的结果

也就是说，我们真正想要的是

$$\mathbf{A}^{-1}\mathbf{b}$$

而不是完整的矩阵 $\mathbf{A}^{-1}$

Krylov subspace method 的核心思想就是：不直接构造 $\mathbf{A}^{-1}$，而是用矩阵-向量乘法和一个低维优化问题，隐式构造出 $\mathbf{A}^{-1}\mathbf{b}$ 的近似

$$\mathbf{A}\mathbf{v}$$

这里矩阵-向量乘法 $\mathbf{A}\mathbf{v}$ 很关键，因为对于稀疏矩阵来说，它通常比矩阵分解或者显式求逆便宜得多

为了看清楚这种方法为什么可能成立，最好先从一个理想化的谱分解视角开始

从多项式逼近逆矩阵开始

先假设 $\mathbf{A}$ 是可对角化的，也就是

$$\mathbf{A} = \mathbf{V}\mathbf{\Lambda}\mathbf{V}^{-1}$$

其中

$$\mathbf{\Lambda} = \operatorname{diag}(\lambda_1,\lambda_2,\dots,\lambda_n)$$

先看为什么逆矩阵会对应特征值的倒数。若

$$\mathbf{A}\mathbf{v}_i = \lambda_i\mathbf{v}_i$$

并且 $\lambda_i\neq 0$，那么两边同时左乘 $\mathbf{A}^{-1}$，得到

$$\begin{align*} \mathbf{A}^{-1}\mathbf{A}\mathbf{v}_i &= \mathbf{A}^{-1}(\lambda_i\mathbf{v}_i) \\ \mathbf{v}_i &= \lambda_i\mathbf{A}^{-1}\mathbf{v}_i \end{align*}$$

所以

$$\mathbf{A}^{-1}\mathbf{v}_i = \frac{1}{\lambda_i}\mathbf{v}_i$$

也就是说，$\mathbf{A}$ 在方向 $\mathbf{v}_i$ 上放大 $\lambda_i$ 倍，而 $\mathbf{A}^{-1}$ 在同一个方向上就要做相反的事情，也就是缩放 $1/\lambda_i$ 倍

因此如果 $\mathbf{A}$ 可对角化，且所有 $\lambda_i\neq 0$，就有

$$\mathbf{A}^{-1} = \mathbf{V}\mathbf{\Lambda}^{-1}\mathbf{V}^{-1}$$

其中

$$\mathbf{\Lambda}^{-1} = \operatorname{diag} \left( \frac{1}{\lambda_1}, \frac{1}{\lambda_2}, \dots, \frac{1}{\lambda_n} \right)$$

现在我们想做的事情是：不用真的构造 $\mathbf{A}^{-1}$，而是找一个更容易计算的对象来扮演它

一种自然的想法是构造某个 矩阵函数 matrix function

$$f(\mathbf{A})$$

希望它能够近似

$$\mathbf{A}^{-1}$$

这里说的矩阵函数，不是指输入一个矩阵然后输出一个 子空间，而是指输入一个矩阵，输出另一个矩阵

例如对于标量函数

$$f(t)=t^2+3t+1$$

输入一个数 $t$ 会输出一个数；对应地，把矩阵 $\mathbf{A}$ 代进去，就得到矩阵函数

$$f(\mathbf{A}) = \mathbf{A}^2+3\mathbf{A}+\mathbf{I}$$

它输入的是矩阵 $\mathbf{A}$，输出的仍然是一个矩阵

如果我们想让某个矩阵函数 $f(\mathbf{A})$ 像 $\mathbf{A}^{-1}$，最自然的一类候选就是 多项式矩阵函数

$$p(t) = c_0+c_1t+c_2t^2+\cdots+c_kt^k$$

把矩阵 $\mathbf{A}$ 代入这个多项式，得到

$$p(\mathbf{A}) = c_0\mathbf{I} +c_1\mathbf{A} +c_2\mathbf{A}^2 +\cdots +c_k\mathbf{A}^k$$

先看它对单个特征方向的作用。因为

$$\mathbf{A}\mathbf{v}_i = \lambda_i\mathbf{v}_i$$

而对于 多项式 的每一项，实际上就是反复作用 $\mathbf{A}$ ，我们知道这在 特征向量 方向上会得到

$$\mathbf{A}^k\mathbf{v}_i = \lambda_i^k\mathbf{v}_i$$

于是

$$\begin{align*} p(\mathbf{A})\mathbf{v}_i &= \left( c_0\mathbf{I} +c_1\mathbf{A} +c_2\mathbf{A}^2 +\cdots +c_k\mathbf{A}^k \right)\mathbf{v}_i \\ &= \left( c_0+c_1\lambda_i+c_2\lambda_i^2+\cdots+c_k\lambda_i^k \right)\mathbf{v}_i \\ &= p(\lambda_i)\mathbf{v}_i \end{align*}$$

也就是说，$p(\mathbf{A})$ 在特征方向 $\mathbf{v}_i$ 上的作用，就是把这个方向缩放 $p(\lambda_i)$ 倍

如果我们让

$$p(\lambda_i) \approx \lambda_i$$

那么 $p(\mathbf{A})$ 在每个特征方向上的缩放就接近 $\mathbf{A}$，于是可以认为 $p(\mathbf{A})$ 也逼近 $\mathbf{A}$

$$p(\mathbf{A}) \approx \mathbf{A}$$

因为如果一个多项式在所有相关特征值上都给出和 $\mathbf{A}$ 一样的缩放，那么它就在这些特征方向上表达了和 $\mathbf{A}$ 一样的变换

当然，这件事本身没有太大意义，因为 $\mathbf{A}$ 已经在手里了

真正有用的是：如果我们想逼近的不是 $\mathbf{A}$，而是 $\mathbf{A}^{-1}$，那么就要看 $\mathbf{A}^{-1}$ 在同一个特征方向上做了什么

前面已经知道

$$\mathbf{A}^{-1}\mathbf{v}_i = \frac{1}{\lambda_i}\mathbf{v}_i$$

所以如果希望 $p(\mathbf{A})$ 像 $\mathbf{A}^{-1}$，就应该让

$$p(\lambda_i) \approx \frac{1}{\lambda_i}$$

把这件事写成完整的矩阵形式时，要先说明为什么可以把多项式直接作用到 $\mathbf{\Lambda}$ 上

由于

$$\mathbf{A} = \mathbf{V}\mathbf{\Lambda}\mathbf{V}^{-1}$$

所以 $\mathbf{A}$ 的平方是

$$\begin{align*} \mathbf{A}^2 &= (\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{-1}) (\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{-1}) \\ &= \mathbf{V}\mathbf{\Lambda} (\mathbf{V}^{-1}\mathbf{V}) \mathbf{\Lambda}\mathbf{V}^{-1} \\ &= \mathbf{V}\mathbf{\Lambda}^2\mathbf{V}^{-1} \end{align*}$$

同理，反复相乘会得到

$$\mathbf{A}^m = \mathbf{V}\mathbf{\Lambda}^m\mathbf{V}^{-1}$$

于是把 $\mathbf{A}$ 代入多项式

$$p(\mathbf{A}) = c_0\mathbf{I} +c_1\mathbf{A} +c_2\mathbf{A}^2 +\cdots +c_k\mathbf{A}^k$$

就有

$$\begin{align*} p(\mathbf{A}) &= c_0\mathbf{I} +c_1\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{-1} +c_2\mathbf{V}\mathbf{\Lambda}^2\mathbf{V}^{-1} +\cdots +c_k\mathbf{V}\mathbf{\Lambda}^k\mathbf{V}^{-1} \\ &= \mathbf{V} \left( c_0\mathbf{I} +c_1\mathbf{\Lambda} +c_2\mathbf{\Lambda}^2 +\cdots +c_k\mathbf{\Lambda}^k \right) \mathbf{V}^{-1} \\ &= \mathbf{V}p(\mathbf{\Lambda})\mathbf{V}^{-1} \end{align*}$$

所以

$$p(\mathbf{A}) = \mathbf{V}p(\mathbf{\Lambda})\mathbf{V}^{-1}$$

其中 $p(\mathbf{\Lambda})$ 只是把多项式逐个作用到特征值上

$$p(\mathbf{\Lambda}) = \begin{bmatrix} p(\lambda_1) & 0 & \cdots & 0 \\ 0 & p(\lambda_2) & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & p(\lambda_n) \end{bmatrix}$$

如果把每个 $p(\lambda_i)$ 展开，就是

$$p(\mathbf{\Lambda}) = \begin{bmatrix} c_0+c_1\lambda_1+\cdots+c_k\lambda_1^k & 0 & \cdots & 0 \\ 0 & c_0+c_1\lambda_2+\cdots+c_k\lambda_2^k & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & c_0+c_1\lambda_n+\cdots+c_k\lambda_n^k \end{bmatrix}$$

而逆矩阵对应的是

$$\mathbf{\Lambda}^{-1} = \begin{bmatrix} \frac{1}{\lambda_1} & 0 & \cdots & 0 \\ 0 & \frac{1}{\lambda_2} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \frac{1}{\lambda_n} \end{bmatrix}$$

所以如果想让

$$p(\mathbf{A}) = \mathbf{V}p(\mathbf{\Lambda})\mathbf{V}^{-1}$$

逼近

$$\mathbf{A}^{-1} = \mathbf{V}\mathbf{\Lambda}^{-1}\mathbf{V}^{-1}$$

就需要让每一个对角位置上的多项式值都逼近对应特征值的倒数

$$\begin{align*} c_0+c_1\lambda_1+\cdots+c_k\lambda_1^k &\approx \frac{1}{\lambda_1} \\ c_0+c_1\lambda_2+\cdots+c_k\lambda_2^k &\approx \frac{1}{\lambda_2} \\ &\vdots \\ c_0+c_1\lambda_n+\cdots+c_k\lambda_n^k &\approx \frac{1}{\lambda_n} \end{align*}$$

注意这里不是对每个 $\lambda_i$ 分别选择一组新的系数，而是同一组

$$c_0,c_1,\dots,c_k$$

必须同时让这个多项式在所有相关特征值上逼近倒数函数。也就是说，我们是在用一个多项式 $p(t)$ 同时拟合这些点

$$(\lambda_i,\frac{1}{\lambda_i})$$

如果真的能做到这一点，那么由同一组系数 $c_0,c_1,\dots,c_k$ 构造出来的多项式矩阵函数

$$p(\mathbf{A}) = c_0\mathbf{I} +c_1\mathbf{A} +\cdots +c_k\mathbf{A}^k$$

就会在所有这些特征方向上模拟 $\mathbf{A}^{-1}$ 的作用，因此可以认为

$$p(\mathbf{A}) \approx \mathbf{A}^{-1}$$

这就是 Krylov 方法背后的第一个直觉：

如果一个多项式在 $\mathbf{A}$ 的特征值上像倒数函数 $1/t$，那么把 $\mathbf{A}$ 代入这个多项式，就会像 $\mathbf{A}^{-1}$

严格来说，只有当

$$p(\lambda_i) = \frac{1}{\lambda_i}$$

对所有相关特征值都成立时，才有

$$p(\mathbf{A}) = \mathbf{A}^{-1}$$

但实际计算里，我们通常不需要完整恢复 $\mathbf{A}^{-1}$，而只是想求当前这个线性系统的解

$$\mathbf{x} = \mathbf{A}^{-1}\mathbf{b}$$

所以更现实的目标是

$$p(\mathbf{A})\mathbf{b} \approx \mathbf{A}^{-1}\mathbf{b}$$

这一点很重要，因为它说明 Krylov 方法通常是在近似 $\mathbf{A}^{-1}$ 对当前 $\mathbf{b}$ 的作用，而不是显式构造整个逆矩阵

换句话说，如果想让 $p(\mathbf{A})$ 对所有可能的右端项都像 $\mathbf{A}^{-1}$，那就必须在所有相关特征方向上都逼近得很好；但如果只要求当前这个 $\mathbf{b}$ 成立，那么主要需要照顾的是 $\mathbf{b}$ 分量较大的那些特征方向

多项式逼近倒数函数

接下来还有一个自然问题：为什么可以指望一个多项式去逼近倒数函数？

原因是 $1/t$ 虽然不是多项式，但只要讨论的区域远离 $0$，它就是一个光滑函数。光滑函数在有限区间上通常可以被多项式很好地近似

例如在 $t=1$ 附近，可以写成

$$\frac{1}{t} = \frac{1}{1+(t-1)}$$

如果令

$$s=t-1$$

接下来使用几何级数公式

$$\frac{1}{1-r} = 1+r+r^2+r^3+\cdots, \qquad |r|<1$$

由于

$$\frac{1}{1+s} = \frac{1}{1-(-s)}$$

相当于令 $r=-s$，所以

$$\frac{1}{t} = \frac{1}{1+s} = 1-s+s^2-s^3+\cdots$$

把 $s=t-1$ 代回去，就得到一个关于 $t$ 的多项式级数

$$\frac{1}{t} = 1-(t-1)+(t-1)^2-(t-1)^3+\cdots$$

取前有限项，就是一个多项式近似

当然，这个直觉有一个关键限制：特征值不能太靠近 $0$。如果某些 $\lambda_i$ 很接近 $0$，那么 $1/\lambda_i$ 会非常大，低阶多项式就很难稳定逼近它。这也正是病态矩阵难求解的原因之一

当前右端项决定哪些方向重要

继续使用可对角化的直觉。若 $\mathbf{A}\mathbf{v}_i=\lambda_i\mathbf{v}_i$，并且把右端项展开为

$$\mathbf{b} = \sum_i \alpha_i\mathbf{v}_i$$

即将 $\mathbf{b}$ 表示在 $\mathbf{v}_i$ 这组基向量上，$\alpha_i$ 是在这组基上的 坐标

那么真实解是

$$\mathbf{A}^{-1}\mathbf{b} = \sum_i \frac{\alpha_i}{\lambda_i} \mathbf{v}_i$$

而多项式近似给出的是

$$p(\mathbf{A})\mathbf{b} = \sum_i \alpha_i p(\lambda_i)\mathbf{v}_i$$

于是误差为

$$\begin{align*} \mathbf{A}^{-1}\mathbf{b} - p(\mathbf{A})\mathbf{b} &= \sum_i \alpha_i \left( \frac{1}{\lambda_i} - p(\lambda_i) \right) \mathbf{v}_i \end{align*}$$

这说明：并不是所有特征值都同等重要

如果 $\mathbf{b}$ 在某个特征方向上的系数 $\alpha_i$ 很小，那么即使 $p(\lambda_i)$ 没有很好逼近 $1/\lambda_i$，它对当前解的影响也可能很小

所以实际 Krylov 方法更准确地说，是为当前输入 $\mathbf{b}$ 找一个合适的低阶多项式，而不是为所有可能的右端项一次性构造完整的 $\mathbf{A}^{-1}$

为什么不显式做特征值插值

从上面的讨论看，如果能知道所有特征值 $\lambda_i$，似乎可以直接构造插值多项式，让

$$p(\lambda_i) = \frac{1}{\lambda_i}$$

然后得到 $p(\mathbf{A})\approx \mathbf{A}^{-1}$

但这在大型问题里通常不可行

第一，显式求出全部特征值和特征向量本身就很贵，往往已经接近直接求解问题的代价

第二，即使原矩阵 $\mathbf{A}$ 是稀疏的，完整的 $\mathbf{A}^{-1}$ 通常也是稠密的，显式构造它会带来很大的计算和存储代价

第三，我们多数时候只需要当前这个 $\mathbf{b}$ 对应的解 $\mathbf{A}^{-1}\mathbf{b}$，并不需要得到能处理所有右端项的完整逆矩阵

所以 Krylov 方法换了一种做法：不显式求特征值，也不显式构造 $\mathbf{A}^{-1}$，而是只通过反复计算

$$\mathbf{A}\mathbf{b}, \quad \mathbf{A}^2\mathbf{b}, \quad \mathbf{A}^3\mathbf{b}, \quad \dots$$

隐式地寻找一个好的多项式 $p(\mathbf{A})$

这就是 Krylov 子空间出现的地方

Krylov 子空间

前面说过，我们希望用一个低阶多项式矩阵函数

$$p_{j-1}(\mathbf{A}) = c_0\mathbf{I} +c_1\mathbf{A} +\cdots +c_{j-1}\mathbf{A}^{j-1}$$

去近似 $\mathbf{A}^{-1}$ 对 $\mathbf{b}$ 的作用。把它作用到 $\mathbf{b}$ 上，就得到

$$\begin{align*} p_{j-1}(\mathbf{A})\mathbf{b} &= \left( c_0\mathbf{I} +c_1\mathbf{A} +\cdots +c_{j-1}\mathbf{A}^{j-1} \right)\mathbf{b} \\ &= c_0\mathbf{b} +c_1\mathbf{A}\mathbf{b} +\cdots +c_{j-1}\mathbf{A}^{j-1}\mathbf{b} \end{align*}$$

因此任何形如 $p_{j-1}(\mathbf{A})\mathbf{b}$ 的向量，一定都落在下面这些向量张成的空间里

第 $j$ 阶 Krylov 子空间定义为

$$K_j(\mathbf{A},\mathbf{b}) = \text{span} \{ \mathbf{b}, \mathbf{A}\mathbf{b}, \mathbf{A}^2\mathbf{b}, \dots, \mathbf{A}^{j-1}\mathbf{b} \}$$

这个空间里的任意向量都可以写成

$$\begin{align*} \mathbf{x}_j &= c_0\mathbf{b} +c_1\mathbf{A}\mathbf{b} +c_2\mathbf{A}^2\mathbf{b} +\cdots +c_{j-1}\mathbf{A}^{j-1}\mathbf{b} \\ &= p_{j-1}(\mathbf{A})\mathbf{b} \end{align*}$$

其中

$$p_{j-1}(t) = c_0+c_1t+c_2t^2+\cdots+c_{j-1}t^{j-1}$$

所以“在 Krylov 子空间里找近似解”，本质上就是“在所有次数不超过 $j-1$ 的多项式里，找一个最适合当前 $\mathbf{b}$ 的多项式”

换句话说，Krylov 方法不是显式求特征值再插值，而是通过

$$\mathbf{b},\quad \mathbf{A}\mathbf{b},\quad \mathbf{A}^2\mathbf{b},\quad \dots$$

这些矩阵-向量乘法，隐式生成一族候选多项式近似

上面的讨论默认了初始解是 $\mathbf{x}_0=\mathbf{0}$，所以我们直接用 $p_{j-1}(\mathbf{A})\mathbf{b}$ 去近似 $\mathbf{A}^{-1}\mathbf{b}$

但在实际迭代方法里，我们经常已经有一个当前近似解 $\mathbf{x}_0$。这时新的问题不是“从零开始找整个解”，而是“在 $\mathbf{x}_0$ 的基础上再补一个修正量”

所以接下来要把

$$\mathbf{b}$$

换成当前还没有被解释掉的那部分，也就是初始 residual

从 b 推广到一般初始解

如果初始解不是 $\mathbf{0}$，而是某个已有近似解 $\mathbf{x}_0$，那么我们不再直接从 $\mathbf{b}$ 出发，而是先看当前还差多少

定义初始 residual

$$\mathbf{r}_0 = \mathbf{b}-\mathbf{A}\mathbf{x}_0$$

这个式子的含义很直接：当前近似解 $\mathbf{x}_0$ 经过矩阵 $\mathbf{A}$ 作用后，得到的是 $\mathbf{A}\mathbf{x}_0$，但我们真正希望得到的是 $\mathbf{b}$，两者之间还差的那一部分就是 $\mathbf{r}_0$

如果要把当前近似解 $\mathbf{x}_0$ 修正成真解，可以设还需要补上的量为

$$\mathbf{e}_0$$

并写成

$$\mathbf{x}_\ast = \mathbf{x}_0+\mathbf{e}_0$$

把它代回原方程

$$\mathbf{A}\mathbf{x}_\ast=\mathbf{b}$$

得到

$$\begin{align*} \mathbf{A}(\mathbf{x}_0+\mathbf{e}_0) &= \mathbf{b} \\ \mathbf{A}\mathbf{x}_0+\mathbf{A}\mathbf{e}_0 &= \mathbf{b} \\ \mathbf{A}\mathbf{e}_0 &= \mathbf{b}-\mathbf{A}\mathbf{x}_0 \\ &= \mathbf{r}_0 \end{align*}$$

所以真实修正量满足

$$\mathbf{e}_0 = \mathbf{A}^{-1}\mathbf{r}_0$$

此时 Krylov 方法要近似的不是 $\mathbf{A}^{-1}\mathbf{b}$，而是这个修正量

$$\mathbf{A}^{-1}\mathbf{r}_0$$

因此搜索空间变成

$$K_j(\mathbf{A},\mathbf{r}_0) = \text{span} \{ \mathbf{r}_0, \mathbf{A}\mathbf{r}_0, \mathbf{A}^2\mathbf{r}_0, \dots, \mathbf{A}^{j-1}\mathbf{r}_0 \}$$

当我们选择用 $K_j(\mathbf{A},\mathbf{r}_0)$ 作为 $j$ 维搜索空间时，候选近似解也就被限制在

$$\mathbf{x}_j \in \mathbf{x}_0+K_j(\mathbf{A},\mathbf{r}_0)$$

也就是说，先在 $\mathbf{x}_0$ 的基础上加一个 Krylov 子空间里的修正量

$$\mathbf{x}_j = \mathbf{x}_0+\mathbf{s}_j, \qquad \mathbf{s}_j\in K_j(\mathbf{A},\mathbf{r}_0)$$

由于 $\mathbf{s}_j$ 在 Krylov 子空间里，所以它一定可以写成

$$\mathbf{s}_j = p_{j-1}(\mathbf{A})\mathbf{r}_0$$

这就是 Krylov 方法和多项式之间的关系：

Krylov 子空间里的每一个候选修正量，都是某个低阶多项式 $p_{j-1}(\mathbf{A})$ 作用在初始 residual 上得到的向量

为什么需要 Gram-Schmidt

到这里为止，我们已经把 Krylov 方法理解成了一个多项式搜索问题：候选修正量都来自

$$p_{j-1}(\mathbf{A})\mathbf{r}_0$$

也就是来自

$$\mathbf{r}_0,\quad \mathbf{A}\mathbf{r}_0,\quad \mathbf{A}^2\mathbf{r}_0,\quad \dots$$

这些向量的线性组合

因此，一个自然想法是直接使用

$$\mathbf{r}_0,\quad \mathbf{A}\mathbf{r}_0,\quad \mathbf{A}^2\mathbf{r}_0,\quad \dots$$

作为基底

但这样做在数值上很不稳定

这些向量很容易变得接近线性相关。尤其当 $\mathbf{A}$ 反复作用以后，某些方向会被不断放大，另一些方向会被压小，最后得到的一组向量可能几乎都指向相近的方向

这会导致两个问题

• 子空间的表示不稳定
• 在这个空间里求系数时会变得病态

所以我们不直接使用原始 Krylov 序列，而是用 Gram-Schmidt 把它们正交化，得到一组正交归一基

$$\mathbf{q}_1,\mathbf{q}_2,\dots,\mathbf{q}_j$$

并记

$$\mathbf{Q}_j = [ \mathbf{q}_1,\dots,\mathbf{q}_j ]$$

它们满足

$$\mathbf{Q}_j^\top\mathbf{Q}_j = \mathbf{I}$$

同时张成同一个 Krylov 子空间

$$K_j(\mathbf{A},\mathbf{r}_0) = \text{span}\{\mathbf{q}_1,\dots,\mathbf{q}_j\}$$

因此 Gram-Schmidt 的作用并非直接求出解，而是稳定地组织搜索空间

GS 造空间的正交基，优化问题决定在这个空间里的具体位置

Arnoldi 过程

在 Krylov 方法里，对 Krylov 序列做 Gram-Schmidt 正交化的标准过程叫做 Arnoldi iteration

前面已经说明，一般初始解 $\mathbf{x}_0$ 下的搜索空间从初始 residual 开始：

$$K_j(\mathbf{A},\mathbf{r}_0) = \text{span} \{ \mathbf{r}_0, \mathbf{A}\mathbf{r}_0, \dots, \mathbf{A}^{j-1}\mathbf{r}_0 \}$$

所以 Arnoldi 要做的第一件事，就是把这个空间里的第一个方向 $\mathbf{r}_0$ 变成单位向量。这个长度记作

先令

$$\beta=\|\mathbf{r}_0\|_2, \qquad \mathbf{q}_1=\frac{\mathbf{r}_0}{\beta}$$

也就是说，$\mathbf{q}_1$ 和 $\mathbf{r}_0$ 指向同一个方向，只是长度被归一化成 $1$

如果 $\beta=0$，说明

$$\mathbf{r}_0=\mathbf{b}-\mathbf{A}\mathbf{x}_0=\mathbf{0}$$

此时 $\mathbf{x}_0$ 已经是精确解，不需要继续构造 Krylov 子空间。下面默认 $\beta\neq 0$

假设现在已经有了

$$\mathbf{q}_1,\mathbf{q}_2,\dots,\mathbf{q}_j$$

下一步不是把当前近似解 $\mathbf{x}_j$ 乘以 $\mathbf{A}$，而是把最后一个 Krylov 正交基向量 $\mathbf{q}_j$ 乘以 $\mathbf{A}$

$$\mathbf{v} = \mathbf{A}\mathbf{q}_j$$

然后把 $\mathbf{v}$ 对所有已有基向量都做正交化

$$h_{ij} = \mathbf{q}_i^\top \mathbf{v}, \qquad i=1,\dots,j$$

并更新

$$\mathbf{v} \leftarrow \mathbf{v} - \sum_{i=1}^j h_{ij}\mathbf{q}_i$$

注意这里不是只和上一个向量 $\mathbf{q}_j$ 正交化，而是要和所有已有向量

$$\mathbf{q}_1,\dots,\mathbf{q}_j$$

正交化

这样才能保证新的方向同时垂直于整个旧空间

最后令

$$h_{j+1,j} = \|\mathbf{v}\|_2$$

如果 $h_{j+1,j}\neq 0$，就归一化得到

$$\mathbf{q}_{j+1} = \frac{\mathbf{v}}{h_{j+1,j}}$$

所以一轮 Arnoldi 的结构可以概括为

$$\mathbf{q}_j \xrightarrow{\times\mathbf{A}} \mathbf{A}\mathbf{q}_j \xrightarrow{\text{GS against } \mathbf{q}_1,\dots,\mathbf{q}_j} \mathbf{q}_{j+1}$$

Arnoldi 过程还会得到一个重要关系

$$\mathbf{A}\mathbf{Q}_j = \mathbf{Q}_{j+1}\mathbf{H}_j$$

这个式子其实只是把前面每一步 Gram-Schmidt 正交化的结果，用矩阵形式压缩写在一起

先看单独一列。为了生成下一个 Krylov 基向量，我们从已有的最后一个基向量 $\mathbf{q}_j$ 出发，先计算

$$\mathbf{A}\mathbf{q}_j$$

然后把它在已有基底

$$\mathbf{q}_1,\dots,\mathbf{q}_j$$

上的分量减掉，最后留下一个新的正交方向 $\mathbf{q}_{j+1}$

所以这一整步等价于把 $\mathbf{A}\mathbf{q}_j$ 写成

$$\mathbf{A}\mathbf{q}_j = h_{1j}\mathbf{q}_1 +h_{2j}\mathbf{q}_2 +\cdots +h_{jj}\mathbf{q}_j +h_{j+1,j}\mathbf{q}_{j+1}$$

这里前 $j$ 项是它在旧空间里的投影，而最后一项是新产生的正交方向

同样地，对更早的列也成立。例如

$$\mathbf{A}\mathbf{q}_1 = h_{11}\mathbf{q}_1+h_{21}\mathbf{q}_2$$

以及

$$\mathbf{A}\mathbf{q}_2 = h_{12}\mathbf{q}_1+h_{22}\mathbf{q}_2+h_{32}\mathbf{q}_3$$

可以看出，第 $k$ 列一般满足

$$\mathbf{A}\mathbf{q}_k = \sum_{i=1}^{k+1} h_{ik}\mathbf{q}_i$$

也就是说，$\mathbf{A}$ 作用在第 $k$ 个基向量上之后，结果仍然落在

$$\text{span}\{\mathbf{q}_1,\dots,\mathbf{q}_{k+1}\}$$

里

现在把这些列并排放在一起。左边变成

$$\mathbf{A} [ \mathbf{q}_1,\dots,\mathbf{q}_j ] = \mathbf{A}\mathbf{Q}_j$$

右边则变成

$$[ \mathbf{q}_1,\dots,\mathbf{q}_{j+1} ] \begin{bmatrix} h_{11} & h_{12} & \cdots & h_{1j} \\ h_{21} & h_{22} & \cdots & h_{2j} \\ 0 & h_{32} & \cdots & h_{3j} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & h_{j+1,j} \end{bmatrix} = \mathbf{Q}_{j+1}\mathbf{H}_j$$

于是就得到

$$\mathbf{A}\mathbf{Q}_j = \mathbf{Q}_{j+1}\mathbf{H}_j$$

这里 $\mathbf{H}_j$ 之所以是 upper Hessenberg matrix，原因也很直接：第 $k$ 列只会用到

$$\mathbf{q}_1,\dots,\mathbf{q}_{k+1}$$

所以在第 $k+1$ 行以下全都是 $0$

这个关系的意义是：虽然原矩阵 $\mathbf{A}$ 可能很大，但它作用在当前 Krylov 正交基上之后，完全可以用一个很小的 upper Hessenberg matrix $\mathbf{H}_j$ 来记录这些作用结果在基底中的坐标

这可以理解为：巨大的矩阵 $\mathbf{A}$ 被投影到 Krylov 子空间后，变成了小矩阵 $\mathbf{H}_j$

在子空间里求最优系数

有了正交基 $\mathbf{Q}_j$ 以后，近似解写成

$$\mathbf{x}_j = \mathbf{x}_0+\mathbf{Q}_j\mathbf{c}$$

这是因为前面已经规定，当我们使用 $K_j(\mathbf{A},\mathbf{r}_0)$ 作为 $j$ 维搜索空间时，近似解只允许从下面这个仿射空间中选取

$$\mathbf{x}_0+K_j(\mathbf{A},\mathbf{r}_0)$$

也就是说，任何候选解都必须能写成

$$\mathbf{x}_j=\mathbf{x}_0+\mathbf{s}_j$$

其中

$$\mathbf{s}_j\in K_j(\mathbf{A},\mathbf{r}_0)$$

而现在 $\mathbf{Q}_j=[\mathbf{q}_1,\dots,\mathbf{q}_j]$ 的列向量正好就是这个 Krylov 子空间的一组正交基，所以子空间里的任意向量都可以写成

$$\mathbf{s}_j = c_1\mathbf{q}_1+\cdots+c_j\mathbf{q}_j = \mathbf{Q}_j\mathbf{c}$$

因此

$$\mathbf{x}_j = \mathbf{x}_0+\mathbf{Q}_j\mathbf{c}$$

这里的 $\mathbf{c}$ 并不是原问题中的解，而是“在 Krylov 基底坐标系下，这个近似解位于什么位置”的系数向量

这里 $\mathbf{c}\in\mathbb{R}^j$ 是一个小向量

原始问题里 $\mathbf{x}$ 可能有几百万个分量，但如果我们只使用 $K_j$ 这个 $j$ 维搜索空间，那么 $\mathbf{c}$ 只有 $j$ 个分量

所以问题从

$$\text{在 } \mathbb{R}^n \text{ 中找 } \mathbf{x}$$

变成了

$$\text{在 } \mathbb{R}^j \text{ 中找 } \mathbf{c}$$

接下来必须选择最好的 $\mathbf{c}$，而不是构造完 $\mathbf{Q}_j$ 就自动得到解

从多项式视角看，选择 $\mathbf{c}$ 就是在选择这些 Krylov 基向量的线性组合，也就是在当前低维空间里选出最合适的多项式近似

对 GMRES 来说，“最好” 的意思是让 residual 的 2-范数最小

$$\mathbf{c}_j = \arg\min_{\mathbf{c}} \| \mathbf{b} - \mathbf{A} (\mathbf{x}_0+\mathbf{Q}_j\mathbf{c}) \|_2$$

然后

$$\mathbf{x}_j = \mathbf{x}_0+\mathbf{Q}_j\mathbf{c}_j$$

如果 $\mathbf{x}_0=\mathbf{0}$，这就变成

$$\mathbf{c}_j = \arg\min_{\mathbf{c}} \| \mathbf{b} - \mathbf{A}\mathbf{Q}_j\mathbf{c} \|_2$$

也就是在 Krylov 子空间中寻找残差最小的那个向量

利用 Arnoldi 关系

$$\mathbf{A}\mathbf{Q}_j = \mathbf{Q}_{j+1}\mathbf{H}_j$$

这个大问题还可以进一步化成一个小的 least squares 问题

关键点在于：我们虽然仍然在最小化一个很大的 residual 向量，但这个 residual 已经被限制在 Arnoldi 构造出的低维结构里，因此可以把它改写成一个小坐标向量的长度

先从 residual 本身开始。把

$$\mathbf{x}_j = \mathbf{x}_0+\mathbf{Q}_j\mathbf{c}$$

代回去，就有

$$\begin{align*} \mathbf{b}-\mathbf{A}\mathbf{x}_j &= \mathbf{b}-\mathbf{A}(\mathbf{x}_0+\mathbf{Q}_j\mathbf{c}) \\ &= (\mathbf{b}-\mathbf{A}\mathbf{x}_0)-\mathbf{A}\mathbf{Q}_j\mathbf{c} \\ &= \mathbf{r}_0-\mathbf{A}\mathbf{Q}_j\mathbf{c} \end{align*}$$

所以我们真正要最小化的对象，其实就是

$$\mathbf{r}_0-\mathbf{A}\mathbf{Q}_j\mathbf{c}$$

接下来要做的事情，是把这一项完全改写到 Arnoldi 基底里

因为

$$\mathbf{r}_0 = \beta\mathbf{q}_1 = \mathbf{Q}_{j+1} \left( \beta\mathbf{e}_1 \right)$$

这里的 $\beta=\|\mathbf{r}_0\|_2$，而

$$\mathbf{q}_1=\frac{\mathbf{r}_0}{\beta}$$

所以自然有

$$\mathbf{r}_0=\beta\mathbf{q}_1$$

另一方面，$\mathbf{Q}_{j+1}$ 的第一列就是 $\mathbf{q}_1$

$$\mathbf{Q}_{j+1} = [ \mathbf{q}_1,\mathbf{q}_2,\dots,\mathbf{q}_{j+1} ]$$

而标准基向量

$$\mathbf{e}_1 = \begin{bmatrix} 1 \\ 0 \\ \vdots \\ 0 \end{bmatrix}$$

右乘之后只会取出第一列，因此

$$\mathbf{Q}_{j+1}\mathbf{e}_1=\mathbf{q}_1$$

于是

$$\mathbf{r}_0 = \beta\mathbf{q}_1 = \mathbf{Q}_{j+1}(\beta\mathbf{e}_1)$$

这一步的本质只是：把初始 residual 写成 Arnoldi 基底坐标的形式

所以 residual 可以写成

$$\begin{align*} \mathbf{b} - \mathbf{A} (\mathbf{x}_0+\mathbf{Q}_j\mathbf{c}) &= \mathbf{r}_0-\mathbf{A}\mathbf{Q}_j\mathbf{c} \\ &= \mathbf{Q}_{j+1} \left( \beta\mathbf{e}_1-\mathbf{H}_j\mathbf{c} \right) \end{align*}$$

其中第二步只是同时用了两件事：

$$\mathbf{r}_0=\mathbf{Q}_{j+1}(\beta\mathbf{e}_1)$$

以及

$$\mathbf{A}\mathbf{Q}_j=\mathbf{Q}_{j+1}\mathbf{H}_j$$

于是整个 residual 被写成

$$\mathbf{Q}_{j+1} \left( \beta\mathbf{e}_1-\mathbf{H}_j\mathbf{c} \right)$$

也就是说，高维空间里的 residual 向量，现在等价地表示成了一个低维坐标向量

$$\beta\mathbf{e}_1-\mathbf{H}_j\mathbf{c}$$

由于 $\mathbf{Q}_{j+1}$ 的列正交归一，最小化原 residual 等价于最小化

$$\min_{\mathbf{c}} \| \beta\mathbf{e}_1-\mathbf{H}_j\mathbf{c} \|_2$$

原因是：对任意向量 $\mathbf{y}$，只要 $\mathbf{Q}_{j+1}$ 的列正交归一，就有

$$\|\mathbf{Q}_{j+1}\mathbf{y}\|_2=\|\mathbf{y}\|_2$$

所以左边那个大 residual 的长度

$$\left\| \mathbf{Q}_{j+1} \left( \beta\mathbf{e}_1-\mathbf{H}_j\mathbf{c} \right) \right\|_2$$

和右边这个小向量的长度完全相同

$$\| \beta\mathbf{e}_1-\mathbf{H}_j\mathbf{c} \|_2$$

因此，最小化高维 residual，等价于最小化一个维度只有 $j+1$ 的小 least squares 问题

这就是 GMRES 的核心形式

大的线性系统没有被直接求逆，而是被投影成了一个低维 least squares 问题

子空间维度越大一定越好吗

Krylov 子空间是嵌套的

$$K_1(\mathbf{A},\mathbf{r}_0) \subset K_2(\mathbf{A},\mathbf{r}_0) \subset K_3(\mathbf{A},\mathbf{r}_0) \subset \cdots$$

因为从 $K_j$ 扩展到 $K_{j+1}$，只是向原来的搜索空间里多加入一个方向

$$\mathbf{A}^{j}\mathbf{r}_0$$

所以如果每次扩展子空间后都在当前子空间里做同一种意义下的最优近似，那么理论上的最佳误差不会变差

例如 GMRES 中

$$\mathbf{x}_j = \arg\min_{\mathbf{x}\in \mathbf{x}_0+K_j} \|\mathbf{b}-\mathbf{A}\mathbf{x}\|_2$$

而

$$K_j\subset K_{j+1}$$

所以在更大的空间 $K_{j+1}$ 中优化时，至少还可以选择原来 $K_j$ 中的解

因此理论上有

$$\|\mathbf{b}-\mathbf{A}\mathbf{x}_{j+1}\|_2 \leq \|\mathbf{b}-\mathbf{A}\mathbf{x}_{j}\|_2$$

但是这不意味着每次增加子空间维度都会明显变好

收敛速度取决于很多因素，包括

• $\mathbf{A}$ 的特征值分布
• 特征值是否远离 $0$
• 条件数是否很大
• residual 在不同特征方向上的分量
• 数值正交化误差
• 是否使用了 preconditioner

从多项式逼近角度看，如果 $\mathbf{A}$ 的重要特征值集中在远离 $0$ 的区域，那么 $1/t$ 比较容易被低阶多项式近似，Krylov 方法通常收敛较快

如果特征值非常分散，或者有很多特征值接近 $0$，那么低阶多项式很难同时照顾所有重要方向，收敛就会变慢

和直接分解法的区别

直接分解法的思路是先对矩阵做一次比较昂贵的分解

$$\mathbf{A}=\mathbf{L}\mathbf{U}$$

或者

$$\mathbf{A}=\mathbf{Q}\mathbf{R}$$

之后如果有很多个右端项

$$\mathbf{A}\mathbf{x}_1=\mathbf{b}_1, \qquad \mathbf{A}\mathbf{x}_2=\mathbf{b}_2, \qquad \mathbf{A}\mathbf{x}_3=\mathbf{b}_3$$

就可以重复使用同一套分解

因此直接分解法适合：

• 矩阵规模没有大到无法分解
• 需要求很多个右端项 RHS
• 需要较高精度的确定性求解

而 Krylov 方法更像是为当前的 $\mathbf{b}$ 构造一个专门的低维搜索空间

$$K_j(\mathbf{A},\mathbf{r}_0)$$

如果换一个右端项 $\mathbf{b}$，初始 residual 通常也会变

$$\mathbf{r}_0 = \mathbf{b}-\mathbf{A}\mathbf{x}_0$$

于是 Krylov 子空间一般也要重新构造

所以 Krylov 方法适合：

• $\mathbf{A}$ 非常大
• $\mathbf{A}$ 很稀疏
• 矩阵分解太贵
• 只需要少数右端项 RHS
• 可以接受迭代近似解
• 能够快速计算 $\mathbf{A}\mathbf{v}$

这也是它常被称为 matrix-free 或 matvec-based 方法的原因

总结

Krylov 方法可以压缩成三层直觉

第一，从谱分解角度看，如果一个多项式在 $\mathbf{A}$ 的特征值上满足

$$p(\lambda_i) \approx \frac{1}{\lambda_i}$$

那么就可以期待

$$p(\mathbf{A})\mathbf{b} \approx \mathbf{A}^{-1}\mathbf{b}$$

这说明 Krylov 方法背后是在用低阶多项式扮演 $\mathbf{A}^{-1}$ 对当前 $\mathbf{b}$ 的作用

第二，Krylov 子空间正是这些低阶多项式自然生成的搜索空间

$$\mathbf{b},\mathbf{A}\mathbf{b},\mathbf{A}^2\mathbf{b},\dots$$

或者更一般地构造

$$\mathbf{r}_0,\mathbf{A}\mathbf{r}_0,\mathbf{A}^2\mathbf{r}_0,\dots$$

这些向量张成的 Krylov 子空间等价于一族多项式近似

$$p_{j-1}(\mathbf{A})\mathbf{r}_0$$

第三，Gram-Schmidt / Arnoldi 负责给这个空间构造稳定的正交基，GMRES 这类方法再在这个低维空间里解一个小优化问题，选出最好的系数

所以最重要的分工是：

$$\boxed{ \text{乘 } \mathbf{A} \Rightarrow \text{扩展 Krylov 空间} }$$ $$\boxed{ \text{Gram-Schmidt} \Rightarrow \text{得到稳定的正交基} }$$ $$\boxed{ \text{优化 } \Rightarrow \text{选择最好的系数并得到 } \mathbf{x}_j }$$

它不是显式求特征值、显式插值、再显式构造 $\mathbf{A}^{-1}$，而是通过矩阵-向量乘法、Krylov 投影和低维优化，隐式完成对 $\mathbf{A}^{-1}\mathbf{b}$ 或 $\mathbf{A}^{-1}\mathbf{r}_0$ 的近似

这就是 Krylov 子空间方法为什么能用相对少的矩阵-向量乘法，去近似一个本来需要求逆才能得到的解