CS205A Lecture 13 Conjugate gradients II; Formulation; preconditioning; and variants

课程主页：https://graphics.stanford.edu/courses/cs205a-13-fall/schedule.html

这次回顾第十三讲，这一讲完整介绍了共轭梯度法。

首先回顾梯度下降的残差：

$$\vec{r}_{k} = \vec{b}-A \vec{x}_{k}$$

下面讨论$\vec r_{k+1}$和$\vec r_k$的关系。我们知道更新公式为

$$\vec{x}_{k+1}=\vec{x}_{k}+\alpha_{k+1} \vec{v}_{k+1}$$

其中

$$\alpha_{k+1}=\frac{\vec{v}_{k+1}^{\top} \vec{r}_{k} }{\vec{v}_{k+1}^{\top} A \vec{v}_{k+1} }$$

那么$\vec r_{k+1}$计算方式如下：

$$\begin{aligned} \vec{r}_{k+1} &=\vec{b}-A \vec{x}_{k+1} \\ &=\vec{b}-A\left(\vec{x}_{k}+\alpha_{k+1} \vec{v}_{k+1}\right) \\ &=\left(\vec{b}-A \vec{x}_{k}\right)-\alpha_{k+1} A \vec{v}_{k+1} \\ &=\vec{r}_{k}-\alpha_{k+1} A \vec{v}_{k+1} \end{aligned}$$

特别的，在梯度下降法中，我们取

$$\vec{v}_{k+1} =\vec{r}_{k}$$

那么

$$\vec{r}_{k+1}=\vec{r}_{k}-\alpha_{k+1} A \vec{r}_{k}$$

递推上式得到如下命题：

命题10.3

$$如果对f使用梯度下降法，那么\operatorname{span}\left\{\vec{r}_{0}, \ldots, \vec{r}_{k}\right\}=\operatorname{span}\left\{\vec{r}_{0}, A \vec{r}_{0}, \ldots, A^{k} \vec{r}_{0}\right\}$$

其中

$$f(\vec{x}) = \frac{1}{2} \vec{x}^{\top} A \vec{x}-\vec{b}^{\top} \vec{x}+c$$

注意上述命题不能保证如下事实成立：

$$\vec{x}_{k}-\vec{x}_{0} \neq {\arg \min} _{\vec{v} \in \operatorname{span}\left\{\vec{r}_{0}, A \vec{r}_{0}, \ldots, A^{k-1} \vec{r}_{0}\right\} } f\left(\vec{x}_{0}+\vec{v}\right)$$

因为第$k$步的更新结果能会破坏前$k-1$次的结果。但是，如果找到$A$共轭向量组$\{\vec v_1,\ldots ,\vec v_n\}$，那么结合之前的讨论不难得到如下命题：

命题10.4

$$假设\vec x_k是沿着共轭向量组迭代更新k次的结果，那么\\ \vec{x}_{k}-\vec{x}_{0} = {\arg \min} _ {\vec{v} \in \operatorname{span}\left\{\vec{v}_{1}, \vec{v}_{2}, \ldots, \vec{v}_{k}\right\} } f\left(\vec{x}_{0}+\vec{v}\right)$$

所以接下来只要讨论如何找到$A$共轭向量组即可。

生成$A$共轭方向

利用Gram-Schmidt的思路，不难得到最原始的做法：

$$\begin{aligned} \text{Update search direction (bad Gram-Schmidt step): }&\vec{v}_{k} =A^{k-1} \vec{r}_{0}-\sum_{i<k} \frac{\left\langle A^{k-1} \vec{r}_{0}, \vec{v}_{i}\right\rangle_{A} }{\left\langle\vec{v}_{i}, \vec{v}_{i}\right\rangle_{A} } \vec{v}_{i} \\ \text{Line search: }&\alpha_{k} =\frac{\vec{v}_{k}^{\top} \vec{r}_{k-1} }{\vec{v}_{k}^{\top} A \vec{v}_{k} } \\ \text{Update estimate: }&\vec{x}_{k} =\vec{x}_{k-1}+\alpha_{k} \vec{v}_{k} \\ \text{Update residual: }&\vec{r}_k =\vec{r}_{k-1}-\alpha_{k} A \vec{v}_{k} \end{aligned}$$

不难看出该算法有如下性质：

$$\operatorname{span}\left\{\vec{v}_{1}, \ldots, \vec{v}_{k}\right\}=\operatorname{span}\left\{\vec{r}_{0}, A \vec{r}_{0}, \ldots, A^{k-1} \vec{r}_{0}\right\} \tag 1$$

这个算法有如下两个问题：

1. 当迭代多次以后，$A^{k-1}\vec r_0​$趋近于$A​$最大特征值对应的特征向量，这使得投影的效果越来越差。
2. 为了计算$\vec v_k$，我们需要存储$\vec v_1,\ldots ,\vec v_{k-1}$；因此每一轮比上一轮都需要更多的内存。

上述算法的第一步是将$A^{k-1}\vec r_0 $减去其在$\operatorname{span}\left\{\vec{v}_{1}, \ldots, \vec{v}_{k-1}\right\}$上的投影，结合(1)，不难发现我们只要找到$\vec w$，使其满足如下条件即可

$$\vec{w} \in \operatorname{span}\left\{\vec{r}_{0}, A \vec{r}_{0}, \ldots, A^{k-1} \vec{r}_{0}\right\} \backslash \operatorname{span}\left\{\vec{r}_{0}, A \vec{r}_{0}, \ldots, A^{k-2} \vec{r}_{0}\right\}$$

因为

$$\vec{r}_{k-1} =\vec{r}_{k-2}-\alpha_{k-1} A \vec{v}_{k-1}$$

所以结合(1)，我们有

$$\vec{r}_{k-1} \in \operatorname{span}\left\{\vec{r}_{0}, A \vec{r}_{0}, \ldots, A^{k-1} \vec{r}_{0}\right\} \backslash \operatorname{span}\left\{\vec{r}_{0}, A \vec{r}_{0}, \ldots, A^{k-2} \vec{r}_{0}\right\}$$

所以可以选择$\vec r_k$进行如上操作，于是得到第二个算法：

$$\begin{aligned} \text{Update search direction (bad Gram-Schmidt step): }&\vec{v}_{k} =\vec{r}_{k-1}-\sum_{i<k} \frac{\left\langle \vec{r}_{k-1}, \vec{v}_{i}\right\rangle_{A} }{\left\langle\vec{v}_{i}, \vec{v}_{i}\right\rangle_{A} } \vec{v}_{i} \\ \text{Line search: }&\alpha_{k} =\frac{\vec{v}_{k}^{\top} \vec{r}_{k-1} }{\vec{v}_{k}^{\top} A \vec{v}_{k} } \\ \text{Update estimate: }&\vec{x}_{k} =\vec{x}_{k-1}+\alpha_{k} \vec{v}_{k} \\ \text{Update residual: }&\vec{r}_k =\vec{r}_{k-1}-\alpha_{k} A \vec{v}_{k} \end{aligned}$$

该算法解决了问题1，即不用计算$A^{k-1}\vec r_0$，但是内存问题似乎还没解决。事实上，该算法已经解决内存的问题，因为我们有如下命题：

命题10.5

$$在上述算法的条件下，对\ell<k,我们有\\ \left\langle\vec{r}_{k}, \vec{v}_{\ell}\right\rangle_{A}=0$$

证明：利用数学归纳法证明。

当$k=1$时，无需证明。

当$k>1$时，假设结论对$k’<k$成立，那么因为

$$\vec r_k = \vec r_{k-1} -\alpha_k A\vec v_k$$

所以我们有

$$\begin{aligned} \left\langle\vec{r}_{k}, \vec{v}_{\ell}\right\rangle_{A} &=\left\langle \vec{r}_{k-1}-\alpha_{k} A \vec{v}_{k} , \vec{v}_{\ell}\right\rangle_{A} \\ &=\left\langle\vec{r}_{k-1}, \vec{v}_{\ell}\right\rangle_{A}-\alpha_{k}\left\langle A \vec{v}_{k}, \vec{v}_{\ell}\right\rangle_{A}\\ &=\left\langle\vec{r}_{k-1}, \vec{v}_{\ell}\right\rangle_{A}-\alpha_{k} (A \vec{v}_{k})^{\top} A \vec{v}_{\ell}\\ &=\left\langle\vec{r}_{k-1}, \vec{v}_{\ell}\right\rangle_{A}-\alpha_{k} \vec{v}_{k}^{\top} A^{\top}A \vec{v}_{\ell}\\ &=\left\langle\vec{r}_{k-1}, \vec{v}_{\ell}\right\rangle_{A}-\alpha_{k} (\vec{v}_{k}^{\top}) A(A \vec{v}_{\ell})\\ &=\left\langle\vec{r}_{k-1}, \vec{v}_{\ell}\right\rangle_{A}-\alpha_{k}\left\langle\vec{v}_{k}, A \vec{v}_{\ell}\right\rangle_{A} \end{aligned}$$

如果$\ell <k-1$，那么由归纳假设，第一项为$0$。注意由(1)可得

$$\begin{aligned} \vec v_{\ell} &\in \operatorname{span}\left\{\vec{r}_{0}, A \vec{r}_{0}, \ldots, A^{\ell-1} \vec{r}_{0}\right\}\\ A\vec v_{\ell} &\in \operatorname{span}\left\{A\vec{r}_{0}, A^2 \vec{r}_{0}, \ldots, A^{\ell} \vec{r}_{0}\right\} \subseteq \operatorname{span}\left\{\vec{r}_{0},A\vec{r}_{0}, A^2 \vec{r}_{0}, \ldots, A^{\ell} \vec{r}_{0}\right\}= \operatorname{span}\left\{\vec{v}_{1}, \ldots, \vec{v}_{\ell+1}\right\} \end{aligned}$$

由$\vec v_k$的构造（$A$共轭），我们可得第二项为$0​$。

如果$\ell =k-1$，直接计算可得

$$\left\langle\vec{r}_{k}, \vec{v}_{\ell}\right\rangle_{A} = \left\langle\vec{r}_{k}, \vec{v}_{k-1}\right\rangle_{A}$$

由残差公式可得

$$A \vec{v}_{k-1}=\frac{1}{\alpha_{k-1} }\left(\vec{r}_{k-2}-\vec{r}_{k-1}\right)$$

和$\vec r_k​$做内积得到

$$\left\langle\vec{r}_{k}, \vec{v}_{k-1}\right\rangle_{A}=\frac{1}{\alpha_{k-1} } \vec{r}_{k}^{\top}\left(\vec{r}_{k-2}-\vec{r}_{k-1}\right)$$

由(1)可得

$$\vec{r}_{k-2}-\vec{r}_{k-1}\in \operatorname{span}\left\{\vec{r}_{0}, A \vec{r}_{0}, \ldots, A^{k-1} \vec{r}_{0}\right\}$$

因为

$$\vec r_k = -\nabla f(\vec x_k)$$

结合命题10.4可得我们必然有

$$\vec{r}_{k} \perp \operatorname{span}\left\{\vec{r}_{0}, A \vec{r}_{0}, \ldots, A^{k-1} \vec{r}_{0}\right\}$$

因为如果上式不成立，那么子空间存在一个方向可以减少$f ​$，因此

$$\left\langle\vec{r}_{k}, \vec{v}_{k-1}\right\rangle_{A}=\frac{1}{\alpha_{k-1} } \vec{r}_{k}^{\top}\left(\vec{r}_{k-2}-\vec{r}_{k-1}\right)=0$$

有了上述命题，我们可以对第二个算法的第一步进行化简：

$$\begin{aligned} \vec{v}_{k} &=\vec{r}_{k-1}-\sum_{i<k} \frac{\left\langle\vec{r}_{k-1}, \vec{v}_{i}\right\rangle_{A} }{\left\langle\vec{v}_{i}, \vec{v}_{i}\right\rangle_{A} } \vec{v}_{i} \\ &=\vec{r}_{k-1}-\frac{\left\langle\vec{r}_{k-1}, \vec{v}_{k-1}\right\rangle_{A} }{\left\langle\vec{v}_{k-1}, \vec{v}_{k-1}\right\rangle_A} \vec{v}_{k-1} \end{aligned}$$

将上述内容总结即可得到最终算法：

$$\begin{aligned} \text{Update search direction: }&\vec{v}_{k} =\vec{r}_{k-1}-\frac{\left\langle\vec{r}_{k-1}, \vec{v}_{k-1}\right\rangle_{A} }{\left\langle\vec{v}_{k-1}, \vec{v}_{k-1}\right\rangle_A} \vec{v}_{k-1} \\ \text{Line search: }&\alpha_{k} =\frac{\vec{v}_{k}^{\top} \vec{r}_{k-1} }{\vec{v}_{k}^{\top} A \vec{v}_{k} } \\ \text{Update estimate: }&\vec{x}_{k} =\vec{x}_{k-1}+\alpha_{k} \vec{v}_{k} \\ \text{Update residual: }&\vec{r}_k =\vec{r}_{k-1}-\alpha_{k} A \vec{v}_{k} \end{aligned}$$

最后补充一个命题：

命题10.6

$$\forall k \leq i, \vec{r}_{i} \cdot \vec{v}_{k}=0$$

证明：由之前讨论可得

$$\vec r_k=-\nabla f\left(\vec{x}_{k}\right) \perp \operatorname{span}\left\{\vec{r}_{0}, A \vec{r}_{0}, \ldots, A^{k-1} \vec{r}_{0}\right\}= \operatorname{span}\left\{\vec{v}_{1}, \ldots, \vec{v}_{k}\right\}$$

所以结论得证。

共轭梯度法分析

现在对上述算法稍作分析，由命题10.6可得

$$\vec v_{k-1} ^{\top} \vec r_{k-1}=0$$

因此由$\vec v_k$的定义，可得

$$\alpha_{k} =\frac{\vec{v}_{k}^{\top} \vec{r}_{k-1} }{\vec{v}_{k}^{\top} A \vec{v}_{k} } =\frac{ \vec{r}_{k-1}^{\top} \vec{r}_{k-1} }{\vec{v}_{k}^{\top} A \vec{v}_{k} } \ge 0$$

记

$$\begin{aligned} \beta_{k} & =-\frac{\left\langle\vec{r}_{k-1}, \vec{v}_{k-1}\right\rangle_{A} }{\left\langle\vec{v}_{k-1}, \vec{v}_{k-1}\right\rangle_{A} } \\ \vec{v}_{k} &=\vec{r}_{k-1}+\beta_{k} \vec{v}_{k-1} \end{aligned}$$

化简$\beta_k$可得

$$\begin{aligned} \beta_{k} & =-\frac{\left\langle\vec{r}_{k-1}, \vec{v}_{k-1}\right\rangle_{A} }{\left\langle\vec{v}_{k-1}, \vec{v}_{k-1}\right\rangle_{A} } \\ &=-\frac{\vec{r}_{k-1}^{\top} A \vec{v}_{k-1} }{\vec{v}_{k-1}^{\top} A \vec{v}_{k-1} } \\ &=-\frac{\vec{r}_{k-1}^{\top}\left(\vec{r}_{k-2}-\vec{r}_{k-1}\right)}{\alpha_{k-1} \vec{v}_{k-1}^{\top} A \vec{v}_{k-1} }& \vec{r}_{k-1}=\vec{r}_{k-2}-\alpha_{k-1} A \vec{v}_{k-1}\\ &=\frac{\vec{r}_{k-1}^{\top} \vec{r}_{k-1} }{\alpha_{k-1} \vec{v}_{k-1}^{\top} A \vec{v}_{k-1} } & \vec{r}_{k-1}^{\top}\vec{r}_{k-2} =0,后续证明这点 \\ &=\frac{\vec{r}_{k-1}^{\top} \vec{r}_{k-1} }{\vec{r}_{k-2}^{\top} \vec{r}_{k-2} } & \alpha_{k} \vec{v}_{k}^{\top} A \vec{v}_{k} = \vec{r}_{k-1}^{\top} \vec{r}_{k-1},\alpha_k的定义 \end{aligned}$$

所以

$$\beta_k \ge 0$$

最后补充证明

$$\vec{r}_{k-1}^{\top}\vec{r}_{k-2} =0$$

原因如下

$$\begin{aligned} \vec{r}_{k-2}^{\top} \vec{r}_{k-1} &=\vec{r}_{k-2}^{\top}\left(\vec{r}_{k-2}-\alpha_{k-1} A \vec{v}_{k-1}\right)\\ &=\vec{r}_{k-2}^{\top} \vec{r}_{k-2}-\frac{\vec{r}_{k-2}^{\top} \vec{r}_{k-2} }{\vec{v}_{k-1}^{\top} A \vec{v}_{k-1} } \vec{r}_{k-2}^{\top} A \vec{v}_{k-1} &\alpha_{k-1}=\frac{ \vec{r}_{k-2}^{\top} \vec{r}_{k-2} }{\vec{v}_{k-1}^{\top} A \vec{v}_{k-1} }\\ & = \vec{r}_{k-2}^{\top} \vec{r}_{k-2}-\frac{\vec{r}_{k-2}^{\top} \vec{r}_{k-2} } {\vec{v}_{k-1}^{\top} A \vec{v}_{k-1} } (\vec{v}_{k-1}-\beta_{k-1} \vec{v}_{k-2})^{\top} A \vec{v}_{k-1} &\vec{v}_{k-1} =\vec{r}_{k-2}+\beta_{k-1} \vec{v}_{k-2}\\ &=\vec{r}_{k-2}^{\top} \vec{r}_{k-2}-\frac{\vec{r}_{k-2}^{\top} \vec{r}_{k-2} }{\vec{v}_{k-1}^{\top} A \vec{v}_{k-1} } \vec{v}_{k-1}^{\top} A \vec{v}_{k-1} & \vec v_{k-2},\vec v_{k-1}为A共轭\\ &=0 \end{aligned}$$

将上述内容总结成如下算法：

$$\begin{aligned} \text{Update search direction: }\beta_{k} &=\frac{\vec{r}_{k-1}^{\top} \vec{r}_{k-1} }{\vec{r}_{k-2}^{\top} \vec{r}_{k-2} } \\ \vec{v}_{k} &=\vec{r}_{k-1}+\beta_{k} \vec{v}_{k-1} \\ \text{Line search: }\alpha_{k} &=\frac{\vec{r}_{k-1}^{\top} \vec{r}_{k-1} }{\vec{v}_{k}^{\top} A \vec{v}_{k} } \\ \text{Update estimate: }\vec{x}_{k} &=\vec{x}_{k-1}+\alpha_{k} \vec{v}_{k} \\ \text{Update residual: }\vec{r}_{k} &=\vec{r}_{k-1}-\alpha_{k} A \vec{v}_{k} \end{aligned}$$

收敛性

假设$\vec x’$为极值点，关于收敛性有如下结论：

$$\frac{f\left(\vec{x}_{k}\right)-f\left(\vec{x}'\right)}{f\left(\vec{x}_{0}\right)-f\left(\vec{x}'\right)} \leq 2\left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^{k}$$

其中

$$\kappa = \text{cond }A$$

Preconditioning

我们知道，对于线性系统

$$A\vec x = \vec b$$

如果$P$可逆，那么如下线性系统等价于原始线性系统

$$P A\vec x = P \vec b$$

如果

$$\text{cond }PA\ll \text{cond }A$$

那么新的系统比原始系统更好求解，这种条件下，称$P$为preconditioner。Preconditioner的想法很有吸引力，但是有如下两个问题：

1. 如果$A$正定，但$PA$不一定正定。
2. 我们需要$P$比$A^{-1}$更好计算。

CG with Preconditioning

假设Preconditioner $P$为正定对称矩阵，那么使用Cholesky分解，我们有

$$P^{-1} =EE^{\top}$$

在这个假设下，我们有如下命题：

命题10.7

$$PA的条件数和E^{-1}AE^{-{\top} }相同$$

证明：首先考虑$E^{-1}AE^{-{\top} }$的奇异值，因为$E^{-1}AE^{-{\top} }$对称，所以奇异值等于特征值，所以讨论特征值即可。假设

$$E^{-1}AE^{-{\top} } \vec x =\lambda \vec x \tag 2$$

注意

$$P =E^{-{\top} } E^{-1}$$

所以对(2)左乘$E^{-{\top} }​$得到

$$E^{-\top }E^{-1}AE^{-\top } \vec x =PAE^{-\top } \vec x =\lambda E^{-\top } \vec x$$

定义

$$\vec y =E^{-\top} \vec x$$

那么有

$$PA\vec y =\lambda \vec y$$

所以$PA$特征值和$E^{-1}AE^{-{\top} }$的特征值相同，假设最大最小特征值为$\lambda_1,\lambda_2$，对应的特征向量为$\vec x _1,\vec x_2 $，那么

$$\begin{aligned} \lambda_1\| \vec x \| \le {\| PA \vec x \|} \le \lambda_2\| \vec x \|\\ \lambda_1 \le \frac{\| PA \vec x \|}{\| \vec x \|} \le \lambda_2 \end{aligned}$$

因此

$$\text{cond }PA = \frac{\lambda_2}{\lambda_1} =\text{cond }E^{-1}AE^{-{\top} }$$

最后一个等号是因为奇异值和条件数的关系。

上述事实告诉我们，如果我们求解

$$E^{-1}AE^{-{\top} } \vec x =E^{-1}\vec b$$

在条件数方面等价于求解

$$PA \vec x =P\vec b$$

结合上述讨论，我们可以用如下方式近似求解$A\vec x=\vec b​$：

1. 求解$E^{-1}AE^{-{\top} } \vec y =E^{-1}\vec b​$
2. 求解$\vec x = E^{-\top} \vec y$

将这个步骤和CG结合得到如下算法，将$A$用$E^{-1}AE^{-{\top} } $替换：

$$\begin{aligned} \text{Update search direction: }\beta_{k}&=\frac{\vec{r}_{k-1}^{\top} \vec{r}_{k-1} }{ {\vec{r}_{k-2}^{\top} \vec{r}_{k-2} }} \\ \vec{v}_{k}&=\vec{r}_{k-1}+\beta_{k} \vec{v}_{k-1} \\ \text{Line search: }\alpha_{k}&=\frac{\vec{r}_{k-1} \vec{r}_{k-1} }{\vec{v}_{k}^{\top} E^{-1} A E^{-\top} \vec{v}_{k} } \\ \text{Update estimate: }\vec{y}_{k}&=\vec{y}_{k-1}+\alpha_{k} E^{-1} A E^{-\top} \vec{v}_{k} \\ \text{Update residual: }\vec{r}_{k}&=\vec{r}_{k-1}-\alpha_{k} E^{-1} A E^{-\top} \vec{v}_{k} \end{aligned}$$

如果定义

$$\tilde{r}_{k} = E \vec{r}_{k}, \tilde{v}_{k} = E^{-\top} \vec{v}_{k},\vec{x}_{k}= E \vec{y}_{k}$$

结合

$$P =E^{-{\top} } E^{-1}$$

我们可以将上述算法改写为：

$$\begin{aligned} \text{Update search direction: }\beta_{k}&=\frac{\tilde{r}_{k-1}^{\top} P \tilde{r}_{k-1} }{\tilde{r}_{k-2}^{\top} P \tilde{r}_{k-2} } \\ \tilde{v}_{k}&=P \tilde{r}_{k-1}+\beta_{k} \tilde{v}_{k-1} \\ \text{Line search: }\alpha_{k}&=\frac{\vec{r}_{k-1}^{\top} P \vec{r}_{k-1} }{\tilde{v}_{k}^{\top} A \tilde{v}_{k} }\\ \text{Update estimate: }\vec{x}_{k}&=\vec{x}_{k-1}+\alpha_{k} \tilde{v}_{k} \\ \text{Update residual: }\tilde{r}_{k}&=\tilde{r}_{k-1}-\alpha_{k} A \tilde{v}_{k} \end{aligned}$$

常见的Preconditioners

比较常见Preconditioner为选择对角元或者用系数矩阵近似逆矩阵，具体的见讲义。