共轭梯度法（CG）详解

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

文章目录

共轭梯度法（CG）详解
- 线性共轭梯度法
- - 共轭方向
  - 共轭方向法
  - CG 方法
  - 收敛率
  - 预条件
- 非线性共轭梯度法
- - FR 方法
  - 其他非线性 CG
  - PR+ 方法
  - 重启动

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

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

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

线性共轭梯度法

线性共轭梯度方法是 Hestenes 和 Stiefel 在 20 世纪 50 年代提出来的的一个迭代方法，用于求解正定系数矩阵的线性系统。
假定 AAA 是对称正定的矩阵，求解线性方程组
Ax=bA x=b Ax=b
等价于求解如下凸优化问题：
min⁡ϕ(x)≡12xTAx−bTx\min \phi(x) \equiv \frac{1}{2} x^{T} A x-b^{T} x minϕ(x)≡21xTAx−bTx
该问题的梯度便是原线性系统的残差，
∇ϕ(x)=Ax−b≡r(x)\nabla \phi(x)=A x-b \equiv r(x) ∇ϕ(x)=Ax−b≡r(x)
在 x=xkx=x_kx=xk 点， rk=Axk−br_{k}=A x_{k}-brk=Axk−b。

共轭方向

定义对于非零向量集合 {p0,p1,⋯,pt}\left\{p_{0}, p_{1}, \cdots, p_{t}\right\}{p0,p1,⋯,pt} 关于对称正定矩阵 AAA 是共轭的，若
piTApj=0,for all i≠j.p_{i}^{T} A p_{j}=0, \quad \text { for all } i \neq j . piTApj=0, for all i=j.

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

假设已经有了一组共轭向量，我们把未知量表示为它们的线性组合 x=∑i=1nαipix=\sum_{i=1}^{n} \alpha^{i} p_{i}x=∑i=1nαipi，我们希望能够寻找一组系数，去极小化
ϕ(x)=∑i=1n((αi)22piTApi−αipiTb)\phi(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) ϕ(x)=i=1∑n(2(αi)2piTApi−αipiTb)
求和中的每一项都是独立的，极小化之，那么我们就可以得到
αi=piTbpiTApi\alpha^{i}=\frac{p_{i}^{T} b}{p_{i}^{T} A p_{i}} αi=piTApipiTb

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

从矩阵的角度来看这个问题，我们把自变量做一个变换，
x^=S−1x\hat{x}=S^{-1} x x^=S−1x
其中，SSS 由共轭向量张成，
S=[p0,p1,⋯,pn−1]S=\left[p_{0}, p_{1}, \cdots, p_{n-1}\right] S=[p0,p1,⋯,pn−1]
那么二次问题变为，
ϕ^(x^)≡ϕ(Sx^)=12x^T(STAS)x^−(STb)Tx^\hat{\phi}(\hat{x}) \equiv \phi(S \hat{x})=\frac{1}{2} \hat{x}^{T}\left(S^{T} A S\right) \hat{x}-\left(S^{T} b\right)^{T} \hat{x} ϕ^(x^)≡ϕ(Sx^)=21x^T(STAS)x^−(STb)Tx^
由共轭性，我们知道矩阵 STASS^{T} A SSTAS 是个对角矩阵，那么久变成了一个对角矩阵系数的极简方程。

共轭方向法

所谓的共轭方向法，就是给定初值点 x0x_0x0 和一组共轭方向，我们通过如下方式迭代更新 xkx_kxk：
xk+1=xk+αkpkx_{k+1}=x_{k}+\alpha_{k} p_{k} xk+1=xk+αkpk
αk=−rkTpkpkTApk\alpha_{k}=-\frac{r_{k}^{T} p_{k}}{p_{k}^{T} A p_{k}} αk=−pkTApkrkTpk

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

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

定理： x0∈ℜnx_{0} \in \Re^{n}x0∈ℜn 是任意起点， {xk}\left\{x_{k}\right\}{xk} 通过共轭方向法生成，那么
rkTpi=0,for i=0,1,⋯,k−1,r_{k}^{T} p_{i}=0, \text { for } i=0,1, \cdots, k-1, rkTpi=0, for i=0,1,⋯,k−1,
且 xkx_{k}xk 在集合
{x∣x=x0+span⁡{p0,p1,⋯,pk−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}}.
上，关于 ϕ(x)=12xTAx−bTx\phi(x)=\frac{1}{2} x^{T} A x-b^{T} xϕ(x)=21xTAx−bTx 的极小化。

CG 方法

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

它的思想在于，选取当前共轭方向为负梯度方向和前一个共轭方向的线性组合，
pk=−rk+βkpk−1p_{k}=-r_{k}+\beta_{k} p_{k-1} pk=−rk+βkpk−1
将其左乘 pk−1TAp_{k-1}^{T} Apk−1TA，由 pkp_kpk 与 pk−1p_{k-1}pk−1 的共轭性，可以得到组合系数：
βk=rkTApk−1pk−1TApk−1\beta_{k}=\frac{r_{k}^{T} A p_{k-1}}{p_{k-1}^{T} A p_{k-1}} βk=pk−1TApk−1rkTApk−1
在这个过程中，选择 p0p_0p0 为 x0x_0x0 处负梯度方向，结合前面的介绍，就可以得到线性共轭梯度方法。

注意到梯度和共轭方向的一些关系：
rkTri=0,∀i=0,⋯,k−1span⁡{r0,r1,⋯,rk}=span⁡{r0,Ar0,⋯,Akr0}span⁡{p0,p1,⋯,pk}=span⁡{r0,Ar0,⋯,Akr0}pkTApi=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} rkTrispan{r0,r1,⋯,rk}span{p0,p1,⋯,pk}pkTApi=0,∀i=0,⋯,k−1=span{r0,Ar0,⋯,Akr0}=span{r0,Ar0,⋯,Akr0}=0,∀i=0,1,⋯,k−1.
通过一些简单的推导，替换掉 CG 算法中的一些表达，就得到了如下的 CG 方法的更加经济的实用形式，

收敛率

定义条件数：
κ(A)=∥A∥2∥A−1∥2=λnλ1\kappa(A)=\|A\|_{2}\left\|A^{-1}\right\|_{2}=\frac{\lambda_{n}}{\lambda_{1}} κ(A)=∥A∥2∥∥A−1∥∥2=λ1λn
那么，CG 的收敛率可以表达为：
∥xk−x∗∥A≤2(κ(A)−1κ(A)+1)k∥x0−x∗∥A\left\|x_{k}-x^{*}\right\|_{A} \leq 2\left(\frac{\sqrt{\kappa(A)}-1}{\sqrt{\kappa(A)}+1}\right)^{k}\left\|x_{0}-x^{*}\right\|_{A} ∥xk−x∗∥A≤2(κ(A)+1κ(A)−1)k∥x0−x∗∥A

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

预条件

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

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

注意到，这里没有显式用到 CCC，而是用到了
M=CTCM = C^TCM=CTC
性质中的残差的正交性表达也发生了改变，
riTM−1rj=0for all i≠jr_{i}^{T} M^{-1} r_{j}=0 \text { for all } i \neq j riTM−1rj=0 for all i=j

非线性共轭梯度法

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

fff 此时是非线性函数。

FR 方法

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

对于步长 αk\alpha_kαk，我们需要采取一种线搜索方法沿着 pkp_kpk 去逼近非线性目标函数 fff 的极小。（满足所谓的强 wolfe 条件的步长）
残差 rrr 原来是线性 CG 方法的梯度，现在需要用 fff 的梯度来替代它。

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

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

其他非线性 CG

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

PR 方法：
βk+1PR=∇fk+1T(∇fk+1−∇fk)∥∇fk∥2.\beta_{k+1}^{P R}=\frac{\nabla f_{k+1}^{T}\left(\nabla f_{k+1}-\nabla f_{k}\right)}{\left\|\nabla f_{k}\right\|^{2}} . βk+1PR=∥∇fk∥2∇fk+1T(∇fk+1−∇fk).

HS 方法：
βk+1HS=∇fk+1T(∇fk+1−∇fk)(∇fk+1−∇fk)Tpk\beta_{k+1}^{H S}=\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}} βk+1HS=(∇fk+1−∇fk)Tpk∇fk+1T(∇fk+1−∇fk)

DY 方法：
βk+1DY=∇fk+1T∇fk+1(∇fk+1−∇fk)Tpk\beta_{k+1}^{D Y}=\frac{\nabla f_{k+1}^{T} \nabla f_{k+1}}{\left(\nabla f_{k+1}-\nabla f_{k}\right)^{T} p_{k}} βk+1DY=(∇fk+1−∇fk)Tpk∇fk+1T∇fk+1

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

我们指出以下几点：

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

PR+ 方法

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

重启动

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

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