共轭梯度法（CG）详解

共轭梯度法（CG）详解

这篇文章写得不错，建议收藏。想要了解 CG，把它认认真真读一遍，就很清楚了。

之前写过几个关于共轭梯度法的注记，譬如：

• https://blog.csdn.net/lusongno1/article/details/78550803
• https://blog.csdn.net/lusongno1/article/details/68942821

但事实上很多人反应，看得一头雾水，基于此，本篇文章旨在对于共轭梯度方法从优化的角度给一个干净的描述。

线性共轭梯度法

线性共轭梯度方法是 Hestenes 和 Stiefel 在 20 世纪 50 年代提出来的的一个迭代方法，用于求解正定系数矩阵的线性系统。
假定 $A$ 是对称正定的矩阵，求解线性方程组
$$Ax=b$$
等价于求解如下凸优化问题：
$$\min \varphi (x)\equiv \frac{1}{2}{x}^{T}Ax-b^{T}x$$
该问题的梯度便是原线性系统的残差，
$$\nabla \varphi (x)=Ax-b\equiv r(x)$$
在 $x=x_k$ 点， $r_k=A{x}_k-b$。

共轭方向

定义 对于非零向量集合 { p 0 , p 1 , ⋯   , p t } \left\{p_{0}, p_{1}, \cdots, p_{t}\right\} {p0​,p1​,⋯,pt​} 关于对称正定矩阵 $A$ 是共轭的，若
$$p_i^{T}A{p}_j=0,\text{ for all }i\ne j.$$

容易证明，共轭向量之间是线性独立的。

假设已经有了一组共轭向量，我们把未知量表示为它们的线性组合 $x={\sum }_{i=1}^n{\alpha }^i{p}_i$，我们希望能够寻找一组系数，去极小化
$$\varphi (x)=\sum _{i=1}^n\left(\frac{{\left({\alpha }^i\right)}^2}{2}{p}_i^{T}A{p}_i-{\alpha }^i{p}_i^{T}b\right)$$
求和中的每一项都是独立的，极小化之，那么我们就可以得到
$${\alpha }^i=\frac{p_i^{T}b}{p_i^{T}A{p}_i}$$

通过共轭方向，把一个 n 维问题，拆解成了 n 个一维问题。

从矩阵的角度来看这个问题，我们把自变量做一个变换，
$$\hat{x}=S^{-1}x$$
其中，$S$ 由共轭向量张成，
$$S=\left[p_0,p_1,\cdots ,p_{n-1}\right]$$
那么二次问题变为，
$$\hat{\varphi }(\hat{x})\equiv \varphi (S\hat{x})=\frac{1}{2}{\hat{x}}^T\left(S^{T}AS\right)\hat{x}-{\left(S^{T}b\right)}^T\hat{x}$$
由共轭性，我们知道矩阵 $S^{T}AS$ 是个对角矩阵，那么久变成了一个对角矩阵系数的极简方程。

共轭方向法

所谓的共轭方向法，就是给定初值点 $x_0$ 和一组共轭方向，我们通过如下方式迭代更新 $x_k$：
$$x_{k+1}=x_k+{\alpha }_k{p}_k$$
$${\alpha }_k=-\frac{r_k^T{p}_k}{p_k^{T}A{p}_k}$$

1、这里的步长 ${\alpha }_k$ 是二次函数 $\varphi$ 沿着 $x_k+\alpha p_k$ 的一维的极小化，我们一般称之为精确线搜索步长。
2、理论上，精确线搜索方法至多 n 步收到到线性系统的解。忽略证明。

对于共轭方向法来说，有如下定理。

定理： $x_0\in {\Re }^n$ 是任意起点， { x k } \left\{x_{k}\right\} {xk​} 通过共轭方向法生成，那么
$$r_k^T{p}_i=0,\text{ for }i=0,1,\cdots ,k-1,$$
且 $x_k$ 在集合
{ x ∣ x = x 0 + span ⁡ { p 0 , p 1 , ⋯   , p k − 1 } } . \left\{x \mid x=x_{0}+\operatorname{span}\left\{p_{0}, p_{1}, \cdots, p_{k-1}\right\}\right\} . {x∣x=x0​+span{p0​,p1​,⋯,pk−1​}}.
上，关于 $\varphi (x)=\frac{1}{2}{x}^{T}Ax-b^{T}x$ 的极小化。

CG 方法

共轭方向法的共轭方向如何得到呢？共轭梯度方法（Conjugate Gradient，CG）方法是一个特别的共轭方向法：它的共轭方向是在 $x_k$ 的迭代中一个一个生成出来的，并且 $p_k$ 的计算只用到 $p_{k-1}$。

它的思想在于，选取当前共轭方向为负梯度方向和前一个共轭方向的线性组合，
$$p_k=-r_k+{\beta }_k{p}_{k-1}$$
将其左乘 $p_{k-1}^{T}A$，由 $p_k$ 与 $p_{k-1}$ 的共轭性，可以得到组合系数：
$${\beta }_k=\frac{r_k^{T}A{p}_{k-1}}{p_{k-1}^{T}A{p}_{k-1}}$$
在这个过程中，选择 $p_0$ 为 $x_0$ 处负梯度方向，结合前面的介绍，就可以得到线性共轭梯度方法。

注意到梯度和共轭方向的一些关系：
r k T r i = 0 , ∀ i = 0 , ⋯   , k − 1 span ⁡ { r 0 , r 1 , ⋯   , r k } = span ⁡ { r 0 , A r 0 , ⋯   , A k r 0 } span ⁡ { p 0 , p 1 , ⋯   , p k } = span ⁡ { r 0 , A r 0 , ⋯   , A k r 0 } p k T A p i = 0 , ∀ i = 0 , 1 , ⋯   , k − 1. \begin{aligned} r_{k}^{T} r_{i} &=0, \quad \forall i=0, \cdots, k-1 \\ \operatorname{span}\left\{r_{0}, r_{1}, \cdots, r_{k}\right\} &=\operatorname{span}\left\{r_{0}, A r_{0}, \cdots, A^{k} r_{0}\right\} \\ \operatorname{span}\left\{p_{0}, p_{1}, \cdots, p_{k}\right\} &=\operatorname{span}\left\{r_{0}, A r_{0}, \cdots, A^{k} r_{0}\right\} \\ p_{k}^{T} A p_{i} &=0, \quad \forall i=0,1, \cdots, k-1 . \end{aligned} rkT​ri​span{r0​,r1​,⋯,rk​}span{p0​,p1​,⋯,pk​}pkT​Api​​=0,∀i=0,⋯,k−1=span{r0​,Ar0​,⋯,Akr0​}=span{r0​,Ar0​,⋯,Akr0​}=0,∀i=0,1,⋯,k−1.​
通过一些简单的推导，替换掉 CG 算法中的一些表达，就得到了如下的 CG 方法的更加经济的实用形式，

收敛率

定义条件数：
$$\kappa (A)=\Vert A{\Vert }_2{\Vert A^{-1}\Vert }_2=\frac{{\lambda }_n}{{\lambda }_1}$$
那么，CG 的收敛率可以表达为：
$${\Vert x_k-x^{\ast }\Vert }_A\le 2{\left(\frac{\sqrt{\kappa (A)}-1}{\sqrt{\kappa (A)}+1}\right)}^k{\Vert x_0-x^{\ast }\Vert }_A$$

由表达式可以看出，当 $A$ 条件数很大的时候，前面的系数趋近于 1，收敛速度无法保证。

预条件

所谓的预条件，就是希望对矩阵 $A$ 做一个改造，改进特征值分布，让它的条件数小一些。

具体地，引入一个非奇异矩阵 $C$，做变量替换，
$$\hat{x}=Cx.$$
二次问题就变为了，
$$\hat{\varphi }(\hat{x})=\frac{1}{2}{\hat{x}}^T{\left(C^{-T}A{C}^{-1}\right)}^{-1}\hat{x}-{\left(C^{-T}b\right)}^T\hat{x}$$
其对应的线性系统是，
$$\left(C^{-T}A{C}^{-1}\right)\hat{x}=C^{-T}b$$
我们要做的，就是找一个逆比较好求的 $C$，使得 $C^{-T}A{C}^{-1}$ 特征值分布更集中。落实到实用算法上，得到：

注意到，这里没有显式用到 $C$，而是用到了
$$M=C^{T}C$$
性质中的残差的正交性表达也发生了改变，
$$r_i^T{M}^{-1}{r}_j=0\text{ for all }i\ne j$$

非线性共轭梯度法

求解非线性极小化问题：
$$\min f(x)$$

$f$ 此时是非线性函数。

FR 方法

相对于共轭梯度法，我们做两点改动：

• 对于步长 ${\alpha }_k$，我们需要采取一种线搜索方法沿着 $p_k$ 去逼近非线性目标函数 $f$ 的极小。（满足所谓的强 wolfe 条件的步长）

• 残差 $r$ 原来是线性 CG 方法的梯度，现在需要用 $f$ 的梯度来替代它。

那么我们就得到了第一个非线性共轭梯度法，它是 Fletcher 和 Reeves 在 20 世纪 60 年代搞的。

对于 FR 方法，如果某步的方向不太好或者步长太小，那么下一步的方向和步长也会很糟糕。

其他非线性 CG

除了 PR 方法，我们选取不同的组合系数 $\beta$，就能得到不同的非线性 CG 方法。

PR 方法：
$${\beta }_{k+1}^{PR}=\frac{\nabla f_{k+1}^T\left(\nabla f_{k+1}-\nabla f_k\right)}{{\Vert \nabla f_k\Vert }^2}.$$

HS 方法：
$${\beta }_{k+1}^{HS}=\frac{\nabla f_{k+1}^T\left(\nabla f_{k+1}-\nabla f_k\right)}{{\left(\nabla f_{k+1}-\nabla f_k\right)}^T{p}_k}$$

DY 方法：
$${\beta }_{k+1}^{DY}=\frac{\nabla f_{k+1}^T\nabla f_{k+1}}{{\left(\nabla f_{k+1}-\nabla f_k\right)}^T{p}_k}$$

容易观察到，这四种方法无非是两个分子和两个分母的四种组合。

我们指出以下几点：

• DY 方法是我们所的戴彧虹和袁亚湘老师提出的。
• 对于 $f$ 是强凸的二次问题，若采用精确想搜索，那么 PR-CG 和 HS-CG 是一个东西。
• 数值实验表明，PR 更鲁棒更有效。
• PR 方法其实就是在 FR 的基础上，当遇到前后两步梯度变化比较小的坏条件的时候，重新开始梯度下降的 “重启动” 方法。
• PR 方法可能不收敛。

PR+ 方法

若要保证 $p_k$ 是下降方向，我们只需要为 PR 的 $\beta$ 进行微调：
$${\beta }_{k+1}^{+}=\max \left\{{\beta }_{k+1}^{PR},0\right\}$$
称之为 PR+ 方法。

重启动

一个重启动的方式是，每迭代 $n$ 步，就设置 ${\beta }_k=0$，即取迭代方向为最速下降方向。重启动能抹掉一些旧的信息。但是这种重启动，没有什么实际的意义，只能说是一种理论的贡献。因为大部分情况下 $n$ 很大，可能不需要迭代多少个 $n$ 步，差不多就达到了比较好的逼近解。

另外一个重新启动是基于前后两步的梯度不正交的考虑，即当遇到
$$\frac{\left|\nabla f_k^T\nabla f_{k-1}\right|}{{\Vert \nabla f_k\Vert }^2}\ge 0.1$$
我们进行一个重启动。