CAPCG 复盘（一）

PCG 简介

预处理共轭梯度法 (Preconditioned Conjugate Gradient) 从共轭梯度法 (Conjugate Gradient) 衍生而来，是常用的用于求解线性方程组 $Ax = b$ 的数值解法，目前主要应用于数值计算领域。

Conjugate Gradient

共轭梯度法主要适用于系数矩阵 $A$ 较为稀疏的情况下。如果矩阵 $A$ 不是稀疏的，那么最好的求解方法是对系数矩阵进行分解，这样也可以快速对不同的 $b$ 求得答案，在系数矩阵很大而且稀疏的情况下，分解产生的矩阵可能含有比 $A$ 更多的非零元素，并且在时间和空间上均不具有优势，所以使用迭代法是一种较好的选择。

引理 1

我们称

$$f(x) = \frac{1}{2}x^{T}Ax - b^Tx + c$$

是一个二次型，这是一个标量，其中 $A$ 是系数矩阵，$x$ 和 $b$ 均为向量，$c$ 为常数。

下面将说明如果 $A$ 是一个对称正定矩阵，那么 $f(x)$ 将在 $Ax = b$ 的解处取到最小值。

如果 $A$ 是对称正定矩阵，那么 $f(x)$ 的图像是一个碗状抛物面（有最小值且最小值唯一），我们可以看出：

$$\begin{aligned} \nabla^{2} f(x) &=\frac{1}{2} A^{T}+\frac{1}{2} A \\ &=A \end{aligned}$$

$f(x)$ 的 hessian 矩阵是正定的，那么根据凸函数的二阶充分必要条件，$f(x)$ 是严格凸函数，有最小值且最小值唯一（当然也可以用朴素的方法证明 $x = A^{-1}b$ 处有唯一最小值）。

所以 $f(x)$ 可以在梯度为 0 处取到最小值：

$$\begin{aligned} \nabla f(x) &=\frac{1}{2} A^{T}x+\frac{1}{2} Ax - b \\ &=Ax - b = 0 \end{aligned}$$

这也就证明了引理，即如果 $A$ 是一个对称正定矩阵，那么 $f(x)$ 将在 $Ax = b$ 的解处取到最小值（注意如果 $A$ 不是对称的，那么将会在 $\frac{1}{2}(A^T+A)x=b$ 处取到最小值）。

最速梯度法 (Steepest Descent)

首先定义误差向量 $\mathbf{e_{(i)}} = \mathbf{x_{(i)}} - \mathbf{x}$，这个向量衡量了目前解到正确解的距离。然后定义残差向量 $\mathbf{r_{(i)}} = \mathbf{b} - A\mathbf{x_{(i)}}$，这个向量衡量 $A\mathbf{x_{(i)}}$ 与期望值 $\mathbf{b}$ 之间的误差。注意，我们可以看出 $\mathbf{r_{(i)}} = -A\mathbf{e_{(i)}}$，这可以看作是误差通过线性变换之后到了和 $\mathbf{b}$ 相同的空间中。并且根据上文中 $\nabla f(\mathbf{x}) |_{\mathbf{x} = \mathbf{x_{(i)}}} = A\mathbf{x_{(i)}} - \mathbf{b}$ ，可以发现 $\mathbf{r_{(i)}} = -\nabla f(\mathbf{x}) |_{\mathbf{x}= \mathbf{x_{(i)}}}$ ，而这对应了梯度的反方向，也就是函数 $f(\mathbf{x})$ 下降最快的方向，所以以后我们会将残差和函数梯度下降最快的方向对应起来。

使用最速梯度法，从任意位置出发，我们可以通过迭代到达 $f(x)$ 的最小值点，迭代每次选取的方向是与梯度相反的方向（梯度是函数增长最快的方向），使用公式：

$$\mathbf{x_{(i+1)}} = \mathbf{x_{(i)}} + \alpha \mathbf{r_{i}}$$

其中 $\alpha$ 是步长，那么如何决定步长呢？

我们希望每次在残差的方向上的走一步之后尽可能使 $f$ 接近最小值，所以:

$$\begin{aligned} \frac{d}{d \alpha} f\left(\mathbf{x_{(i+1)}}\right) &= (\frac{d f}{d \mathbf{x_{i+1}}})^T \frac{d \mathbf{x_{i+1}}}{d \alpha}\\ &= -\mathbf{r_{i+1}}^T r_{i} = 0 \end{aligned}$$

这样我们可以看出需要 ${r_{i+1}}, r_{i}$ 相互正交，我们也可以来看一下两次搜索方向相互正交的几何含义：

实线方向即为本次搜索的方向（本次的梯度反方向），我们希望本次搜索到达的这条线上离最小值最近的点，这个点的特征是它的梯度（下次搜索方向）与本次搜索的方向（本次梯度）垂直。如果不垂直（比如直线上的其他点），我们可以按照梯度在直线上投影（虚线）的反方向继续往前走，可以达到更小的值（因为梯度的反方向是函数值变小的方向），以上就是一点 intuition。

进一步推导：

$$\begin{aligned} r_{(i+1)}^{T} r_{(i)} &=0 \\ \left(b-A x_{(i+1)}\right)^{T} r_{(i)} &=0 \\ \left(b-A\left(x_{(i)}+\alpha r_{(i)}\right)\right)^{T} r_{(i)} &=0 \\ \left(b-A x_{(i)}\right)^{T} r_{(i)}-\alpha\left(A r_{(i)}\right)^{T} r_{(i)} &=0 \\ \left(b-A x_{(i)}\right)^{T} r_{(i)} &=\alpha\left(A r_{(i)}\right)^{T} r_{(i)} \\ r_{(i)}^{T} r_{(i)} &=\alpha r_{(i)}^{T}\left(A r_{(i)}\right) \\ \alpha &=\frac{r_{(i)}^{T} r_{(i)}}{r_{(i)}^{T} A r_{(i)}} \end{aligned}$$

总结上述过程可以得到：

$$\begin{aligned} r_{(i)} &=b-A x_{(i)} \\ \alpha_{(i)} &=\frac{r_{(i)}^{T} r_{(i)}}{r_{(i)}^{T} A r_{(i)}} \\ x_{(i+1)} &=x_{(i)}+\alpha_{(i)} r_{(i)} \end{aligned}$$

这就是最速下降法的迭代公式，我们也可以将最后一个公式两侧同时乘 $-A$ 再加上 $b$，得到：

$$r_{(i+1)}=r_{(i)}-\alpha_{(i)} A r_{(i)}$$

我们可以通过最后三个式子迭代直到收敛，这样做比较方便，但有由于 $r$ 每次的迭代都没有 $x$ 的反馈，有可能最后存在一定的误差（浮点数计算误差），我们可以通过周期性用 $r_{(i)} =b-A x_{(i)}$ 更新 $r_{(i)}$ 来校正误差。

共轭梯度法 (Conjugate Gradient Method)

Conjugate Direction

使用最速梯度法的过程中，可能会遇到重复走之前走过的方向这种情况：

我们希望在搜索的方向上一次走合适的距离，之后就不再重复搜索这个方向，提高搜索的效率（有的时候最速下降法会在靠近解的时候来回震荡，这不是我们想要的结果）。所以我们将这样进行迭代

$$x_{(i+1)}=x_{(i)}+\alpha_{(i)} d_{(i)}$$

其中，$\mathbf{d_{(i)}}$ 是每次搜索的方向，如果要求 $\mathbf{e_{(i+1)}}$ 和 $\mathbf{d_{(i)}}$ 正交（这样就将在 $\mathbf{d_{(i)}}$ 的误差降为 0，之后不会再搜索这个方向）

$$\begin{aligned} x_{(i+1)} &=x_{(i)}+\alpha_{(i)} d_{(i)} \\ x_{(i+1)} - x &=x_{(i)} - x +\alpha_{(i)} d_{(i)} \\ e_{(i+1)} &= e_{(i)} + \alpha_{(i)} d_{(i)} \\ \end{aligned}$$

所以有：

$$\begin{aligned} d_{(i)}^{T} e_{(i+1)} &=0 \\ d_{(i)}^{T}\left(e_{(i)}+\alpha_{(i)} d_{(i)}\right) &=0 \quad \\ \alpha_{(i)} &=-\frac{d_{(i)}^{T} e_{(i)}}{d_{(i)}^{T} d_{(i)}} \end{aligned}$$

但这个方法显然存在问题，如果我们知道了 $e_{(i)}$，那么也就相当于知道了问题的答案，也就不用迭代法求解了😂。

为了解决这个问题，我们选择让本次搜索方向和之后的误差关于矩阵 $A$ 正交，而不是直接正交：

$$d_{(i)}^{T} A e_{(i+1)}=0$$

本质上，这就是在搜索方向上找到函数值最小的点，通过推导我们也能看出这个约束：

$$\begin{aligned} \frac{d}{d \alpha} f\left(x_{(i+1)}\right) &=0 \\ f^{\prime}\left(x_{(i+1)}\right)^{T} \frac{d}{d \alpha} x_{(i+1)} &=0 \\ -r_{(i+1)}^{T} d_{(i)} &=0 \\ d_{(i)}^{T} A e_{(i+1)} &=0 \end{aligned}$$

同上面的推导，可以推出：

$$\begin{aligned} \alpha_{(i)} &=-\frac{d_{(i)}^{T} A e_{(i)}}{d_{(i)}^{T} A d_{(i)}} \\ &=\frac{d_{(i)}^{T} r_{(i)}}{d_{(i)}^{T} A d_{(i)}} \end{aligned}$$

这次，表达式变成了我们可以求解的值。本质上，我们通过关于矩阵 $A$ 正交的这个条件，将原先 $\alpha$ 中的 $e_{(i)}$ 转化为了 $Ae_{(i)} = r_{(i)}$ 也就是我们可以通过迭代求解的值。注意，如果某一个 $d_{(i)}$ 选择 $r_{(i)}$，那么迭代公式就和最速梯度法中的公式相同，但实际上参考后面的构建规则，$d_{(i)}$ 和 $r_{(i)}$ 不会一直完全相等。

下面将证明共轭梯度法可以在 n 步之内得到 $x$ 的值，我们记误差为各个搜索方向的线性组合：

$$e_{(0)}=\sum_{j=0}^{n-1} \delta_{j} d_{(j)}$$ $$\begin{aligned} d_{(k)}^{T} A e_{(0)}&=\sum_{j} \delta_{(j)} d_{(k)}^{T} A d_{(j)} \\ d_{(k)}^{T} A e_{(0)}&=\delta_{(k)} d_{(k)}^{T} A d_{(k)} \text{ by A-orthogonality of $d$ vectors}\\ \delta_{(k)} &=\frac{d_{(k)}^{T} A e_{(0)}}{d_{(k)}^{T} A d_{(k)}} \\ &=\frac{d_{(k)}^{T} A\left(e_{(0)}+\sum_{i=0}^{k-1} \alpha_{(i)} d_{(i)}\right)}{d_{(k)}^{T} A d_{(k)}} \quad \text { (by } A \text { -orthogonality of } d \text { vectors) }\\ &=\frac{d_{(k)}^{T} A e_{(k)}}{d_{(k)}^{T} A d_{(k)}}\\ &=-\alpha_{(k)} \end{aligned}$$

事实上，我们可以看出迭代过程就是一项一项地将每个方向上的误差移除掉，最终得到期望值。

构建共轭方向

可以使用 Gram-Schmidt 进行构建，大概的过程就是在一组 n 个线性无关的向量 $u_0, u_1, \dots, u_n$ 中，$d_i$ 由 $u_i$ 减去与所有之前的 $d_j$ 不是 A 正交的分量，即：

$$d_{(i)}=u_{(i)}+\sum_{k=0}^{i-1} \beta_{i k} d_{(k)}$$

为了求出系数 $\beta_{i k}$，我们可以根据定义令其与之前的方向 A-正交，即：

$$\begin{aligned} d_{(i)}^{T} A d_{(j)} &=u_{i}^{T} A d_{(j)}+\sum_{k=0}^{i-1} \beta_{i k} d_{(k)}^{T} A d_{(j)} \\ 0 &=u_{i}^{T} A d_{(j)}+\beta_{i j} d_{(j)}^{T} A d_{(j)}, \quad i>j \quad \text{by A-orthogonality}\\ \beta_{i j} &=-\frac{u_{i}^{T} A d_{(j)}}{d_{(j)}^{T} A d_{(j)}} \end{aligned}$$

一个缺点是这样需要保存之前的所有方向用来构建当前的方向，时间复杂度也较高。值得注意的是，参考文献中表示，当初始的搜索方向选定为坐标轴方向时，共轭方向法进行求解退化为高斯消元（未证明），但共轭梯度的出现解决了这些不足。但无论是共轭梯度还是共轭方向法，都是在 A 矩阵的空间中，进行每次迭代都全部移除该搜索方向上的误差的操作（A-正交）。

共轭方向的重要性质是每一步它都在 $e_0$ 和 $span\{d_{0}, d_{1}, d_{2}, \dots, d_{i} \} = \mathcal{D}_{i}$ 中寻找误差最小的点，即 $\left|e_{(i)}\right|_{A}$ 在 $e_{(0)} + \mathcal{D}_{i}$ 中误差最小：

$$\begin{aligned} e_{(i)} &= \sum_{j=i}^{n-1} \delta_{(j)} d_{(j)}^{T}\\ \left\|e_{(i)}\right\|_{A} &=\sum_{j=i}^{n-1} \sum_{k=i}^{n-1} \delta_{(j)} \delta_{(k)} d_{(j)}^{T} A d_{(k)} \\ &=\sum_{j=i}^{n-1} \delta_{(j)}^{2} d_{(j)}^{T} A d_{(j)} \end{aligned}$$

可以看出 $\left |e_{(i)}\right|_{A}$ 的误差项全部存在于 $\mathcal{R}^n - \mathcal{D}_{i}$，在 $e_{(0)} + \mathcal{D}_{i}$ 中误差最小。

另外一个重要性质是残差 $r_(i)$ 和之前所有的搜索方向正交，因为所有的搜索方向我们只走一次，这意味着误差 $e_{(i)}$ 与之前所有的搜索方向是 A-正交 的，$r_{(i)} = -Ae_{(i)}$，所以可以推出 $r_{(i)}$ 与之前搜索方向正交。也可以从数学上证明：

$$\begin{aligned} e_{(j)} &= \sum_{k=j}^{n-1} \delta_{(k)} d_{(k)}^{T}\\ -d_{(i)}^{T} A e_{(j)} &=-\sum_{k=j}^{n-1} \delta_{(k)} d_{(i)}^{T} A d_{(k)} \\ d_{(i)}^{T} r_{(j)} &=0, \quad i<j \quad \text {A-orthogonality} \end{aligned}$$

由于搜索方向是由一组线性无关的向量 $u_{(i)}$ 构建而来的，我们能够继续推出：

$$\begin{aligned} d_{(i)}&=u_{(i)}+\sum_{k=0}^{i-1} \beta_{i k} d_{(k)}\\ d_{(i)}^{T} r_{(j)} &=u_{i}^{T} r_{(j)}+\sum_{k=0}^{i-1} \beta_{i k} d_{(k)}^{T} r_{(j)} \\ 0 &=u_{i}^{T} r_{(j)}, \quad i<j \end{aligned}$$

特别地，当上式中 $j=i$ 时，有：

$$d_{(i)}^{T} r_{(i)}=u_{i}^{T} r_{(i)}$$

以及，我们的残差迭代公式可以进一步由矩阵和误差相乘改为由自身迭代产生：

$$\begin{aligned} e_{(i+1)} &=e_{(i)}+\alpha_{(i)} d_{(i)} \\ r_{(i+1)} &=-A e_{(i+1)} \\ &=-A\left(e_{(i)}+\alpha_{(i)} d_{(i)}\right) \\ &=r_{(i)}-\alpha_{(i)} A d_{(i)} \end{aligned}$$

参考文献

1. painless-conjugate-gradient.