数值优化（Numerical Optimization）学习系列-共轭梯度方法（Conjugate Gradient）

概述

共轭梯度算法在最优化问题中备受关注，有两层用途，一是可以求解线性方程Ax=bAx = b；二是可以求解最优化问题。相对于最速下降法，它没有额外的矩阵存储并且比更快，一般N步内收敛。实际收敛效率依赖于系数矩阵特征值的分布。

主要介绍一下内容：

1. 线性共轭梯度算法

2. 共轭方向算法

3. 共轭梯度算法

4. 收敛性

5. 非线性共轭梯度算法

线性共轭梯度算法

共轭梯度算法是一个求解线性方程的迭代方法

问题形式

CG算法求解问题的两种形式：

1. 线性方程Ax=bAx=b并且要求A是对称正定矩阵。

2. 最优化问题：minϕ(x)=12xTAx−bTxmin \phi(x) = \frac12 x^TAx-b^Tx，要求A对称正定，这样该问题是一个凸问题并且有最优解，根据最优解满足∇ϕ(x)=Ax−b=0，一般即为r(x)\nabla \phi(x)=Ax-b=0，一般即为r(x)。在迭代过程中第K步的残差为rk=Axk−br_k = Ax_k-b

共轭性

给定一个非零向量集合 {p0,p1,p2...pn−1{p_0, p_1,p_2...p_{n-1}}}和一个对称正定矩阵A；如果向量集合相对于A是共轭的当且仅当 pTiApj=0, i≠jp_i^TAp_j=0, \ i \ne j

如果向量集合是共轭的，则他们之间是相互线性独立。

证明：反证法。如果不是线性独立，则有 pi=λpjp_i=\lambda p_j ，则pTiApj=λpTjApj≠0p_i^TAp_j=\lambda p_j^TAp_j \ne 0，不满足共轭条件。

共轭方向算法

共轭方向算法（Conjugate Direction Method）不同于共轭梯度方法，共轭向量提前给出。

算法描述

给定共轭方向集合{p0,p1,p2...pn−1{p_0, p_1,p_2...p_{n-1}}}和任意初始点x0x_0

计算xk+1=xk+αkpkx_{k+1}=x_k + \alpha_kp_k

计算最优步长α\alpha，即优化ϕ(xk+αpk)\phi(x_k+\alpha p_k)，通过求解单变量最优化问题，可以容易得到最优步长αk=−rTkpkpTkApk\alpha_k=-\frac {r_k^Tp_k}{p_k^TAp_k}

关于αk\alpha_k的计算

αk=arg min(ϕ(xk+αpk))=argmin12(xk+αpk)TA(xk+αpk)−bT(xk+αpk)∇ϕ(α)=0→αk=−rTkpkpTkApk\begin{align}
\alpha_k &=arg \ min(\phi(x_k+\alpha p_k)) \\
&= arg min \frac12 (x_k+\alpha p_k)^TA(x_k+\alpha p_k)-b^T(x_k+\alpha p_k) \\
& \nabla \phi(\alpha)=0 \to \alpha_k=-\frac {r_k^Tp_k}{p_k^TAp_k}
\end{align}

对于任何初始点x0x_0共轭方向算法最多N步可以收敛到最优解x∗x^*

证明：思路通过上述算法可以表示出最优解。

1. 由于pi{p_i}线性独立，则在N维空间中可以生成整个空间，因此x∗−x0=δ0p0+δ1p1+...+δn−1pn−1x^*-x_0=\delta_0p_0+\delta_1p_1+...+\delta_{n-1}p_{n-1}左边同时乘上pTkAp_k^TA得到pTkA(x∗−x0)=pTkA(δ0p0+δ1p1+...+δn−1pn−1)=δkpTkApkp_k^TA(x^*-x_0)=p_k^TA(\delta_0p_0+\delta_1p_1+...+\delta_{n-1}p_{n-1})=\delta_k p_k^TAp_k，可以推出δk=pTkA(x∗−x0)pTkApk\delta_k=\frac{p_k^TA(x^*-x_0)}{p_k^TAp_k}

2. 根据共轭方向算法可以得到xk=x0+α0p0+α1p1+...+αk−1pk−1x_k=x_0+\alpha_0p_0+\alpha_1p_1+...+\alpha_{k-1}p_{k-1}，同样两边同时乘上pTkAp_k^TA可以得到pTkA(xk−x0)=0p_k^TA(x_k-x_0)=0即pTkAxk=pTkAx0p_k^TAx_k=p_k^TAx_0

3. 因此有pTkA(x∗−x0)=pTkA(x∗−xk)=pTk(Ax∗−Axk)=pTk(b−Axk)=−pTkrk=δkpTkApk\begin{align}p_k^TA(x^*-x_0)&=p_k^TA(x^*-x_k) \\ &=p_k^T(Ax^*-Ax_k) \\ &=p_k^T(b-Ax_k) \\ &= -p_k^Tr_k \\ &=\delta_k p_k^TAp_k \end{align}

4. 可见αk=δk\alpha_k=\delta_k，从而任意起点都可以再N步内得到最优解。

算法理解

如果假设矩阵A是对角阵，则我们可以沿着坐标轴方向依次优化各个方向的值。坐标轴方向为e0,e1...en−1e_0,e_1...e_{n-1}。二维情况如下：


但是如果A是非对角阵，如果还沿着坐标轴方向进行优化，有可能不会收敛。

但是可以通过转换的方式将A变成对角阵。

定义x^=S−1x\hat x=S^{-1}x，其中S的定义为S=[p0,p1...pn−1]S=[p_0,p_1...p_{n-1}]

则原问题可以表示为:ϕ^(x^)=ϕ(Sx^)=12x^T(STAS)x−(STb)Tx^\hat \phi(\hat x)=\phi(S\hat x)=\frac12 \hat x^T(S^TAS)x-(S^Tb)^T \hat x

由于共轭性，则有(STAS)(S^TAS)是对角阵，可以按照坐标轴方向进行优化。其中的关键点是

在x^空间内的ei方向等价于在x空间内的pi方向\hat x 空间内的e_i方向等价于在x空间内的p_i方向，从而进一步验证该算法的正确性。

算法性质

从上述算法理解中可以看出，在计算到eke_k方向时，此时的xkx_k是扩展子空间e1...eke_1...e_k中的最优解。类比于共轭方向算法有

对于任意初始值x0x_0，解序列x0,x1...xkx_0,x_1...x_k由共轭方向算法生成，则有以下性质

1. rTkpi=0i=0,1,...,k−1r_k^Tp_i=0 \quad i=0,1,...,k-1

2. xk在以下集合中是ϕ(x)的最优解，{x|x=x0+span(p0,p1...pk−1)}x_k在以下集合中是\phi(x)的最优解，\{x|x=x_0+span(p_0,p_1...p_{k-1})\}

在证明之间简单理解一下，第K步的残差rkr_k和前面K-1个共轭方向是正交的。

证明：用归纳法可以得到。

1. 当k=1时，需要证明rT1p0=0r_1^Tp_0=0，由于x1=x0+α0p0，α0=−rT0p0pT0Ap0x_1=x_0+\alpha_0p_0，\alpha_0=-\frac {r_0^Tp_0}{p_0^TAp_0}，又rT1p0=(Ax1−b)Tp0=(Ax0+α0Ap0−b)Tp0=(r0+α0Ap0)Tp0=rT0p0+α0pT0Ap0=0\begin{align}r_1^Tp_0 &=(Ax_1-b)^Tp_0 \\ &=(Ax_0+\alpha_0Ap_0-b)^Tp_0 \\ &=(r_0+\alpha_0Ap_0)^Tp_0 \\ &=r_0^Tp_0 + \alpha_0p_0^TAp_0 \\ &=0\end{align}

2. 假设k-1时也成立，即rTk−1pi=0; i=0,1,2...k−2r_{k-1}^Tp_i=0 ; \ i=0,1,2...k-2

3. 证明k时也成立，由于rk=Axk−b=A(xk−1+αk−1pk−1)−b=rk−1+αk−1Apk−1r_k=Ax_k-b=A(x_{k-1}+\alpha_{k-1}p_{k-1})-b=r_{k-1}+\alpha_{k-1}Ap_{k-1}

4. pTk−1rk=pk−1rk−1+αk−1pTk−1Apk−1=0p_{k-1}^Tr_k=p_{k-1}r_{k-1}+\alpha_{k-1}p_{k-1}^TAp_{k-1}=0，根据αk−1\alpha_{k-1}的定义可以得到

5. pTirk=pirk−1+αk−1pTiApk−1=0i=0,1...k−2p_i^Tr_k=p_{i}r_{k-1}+\alpha_{k-1}p_{i}^TAp_{k-1}=0 \quad i=0,1...k-2根据归纳假设可以得到。

共轭方向计算

可以直接采用特征向量，但是特征向量计算复杂度比较高

可以采用Gram_Schmidt求解正交分解的方式得到共轭方向，对于大规模数据计算复杂度也高。

共轭梯度算法

共轭梯度算法是共轭方向算法的一种，它提供了一种计算共轭方向集合的方法，并且当前方向的计算仅和上一步方向相关，从而减少矩阵存储。

当前共轭方向pkp_k选择为当前残差方向和上一步方向pk−1p_{k-1}的线性表示，即pk=−rk+βkpk−1p_k=-r_k+\beta_k p_{k-1}。由于要满足pkApk−1=0p_kAp_{k-1}=0，从而可以推出：βk=rTkApk−1pk−1TApk−1\beta_k=\frac{r_k^TAp_{k-1}}{p_{k-1^TAp_{k-1}}}

初始点x0x_0可以随机选择，初始搜索方向p0p_0选择为最速下降方向−r0-r_0

CG-Preliminary算法

初始化x0x_0

r0=Ax0−b,p0=r0,k=0r_0=Ax_0-b,p_0=r_0,k=0

while rk≠0r_k \ne 0

αk=−rTkpkpTkApk\alpha_k=-\frac {r_k^Tp_k}{p_k^TAp_k}

xk+1=xk+αkpkx_{k+1}=x_k + \alpha_kp_k

rk+1=Axk+1−br_{k+1}=Ax_{k+1}-b

βk+1=rTk+1ApkpTkApk\beta_{k+1}=\frac{r_{k+1}^TAp_{k}}{p_{k}^TAp_{k}}

pk+1=−rk+1+βk+1pkp_{k+1}=-r_{k+1}+\beta_{k+1}p_k

k=k+1k=k+1

end while

算法证明

现在需要解决的问题是上述算法的正确性？如何更好的优化？

证明正确性需要解决的问题是根据上述算法求解得到的pkp_k是否满足相对于A是共轭的。

定理：共轭算法产生的解序列{xkx_k}满足一下性质：

1. rTkri=0i=0,1,2...k−1r_k^Tr_i=0 \quad i=0,1,2...k-1

2. pTkApi=0i=0,1,2...k−1p_k^TAp_i=0 \quad i=0,1,2...k-1

即当前步骤产出的残差和之前所有残差正交；当前搜索方向和之前所有方向共轭。

证明：首先回顾一下在共轭方向算法中的一个性质rTkpi=0i=0,1,2...k−1r_k^Tp_i=0 \quad i=0,1,2...k-1即当前步骤的残差方向和之前所有的搜索方向正交。

性质1证明：由于pk=−rk+βkpk−1p_k=-r_k+\beta_kp_{k-1}则有rk=−pk+βkpk−1r_k=-p_k+\beta_kp_{k-1}，两边同时乘上rir_i可以推出rTirk=rTi(−pk+βkpk−1)r_i^Tr_k=r_i^T(-p_k+\beta_kp_{k-1})根据上述性质，显而易见结果为0.

性质2证明：可以采用归纳方法证明。

1. k=0时，p0=−r0p_0=-r_0

2. k=1时，p1=−r1+β1p0p_1=-r_1+\beta_1p_0，根据定义肯定满足pT1Ap0=0p_1^TAp_0=0

3. k=2时，p2=−r2+β2p1p_2=-r_2+\beta_2 p_1，根据定义pT2Ap1=0p_2^TAp_1=0，下面需要证明pT2Ap0=0p_2^TAp_0=0。pT2Ap0=(−r2+β2p1)TAp0=−rT2Ap0+β2p1Ap0=−rT2Ap0+0=−(r1+α1Ap1)TAp0=−(r0+α0Ap0+α1Ap1)TAp0=代入α0 and α1=0\begin{align}p_2^TAp_0 &=(-r_2+\beta_2 p_1)^TAp_0 \\ &=-r_2^TAp_0+\beta_2 p_1Ap_0 \\ &=-r_2^TAp_0+0 \\ &=-(r_1+\alpha_1Ap_1)^TAp_0 \\ &=-(r_0+\alpha_0Ap_0+\alpha_1Ap_1)^TAp_0 \\ &= 代入\alpha_0 \ and \ \alpha_1 \\ &=0\end{align}

CG算法

下面正式介绍CG算法，根据上述性质可以对基础算法进行优化改进。

改进思路1：αk=−rTkpkpTkApk=−rTk(−rk+βkpk−1)pTkApk=rTkrkpTkApk\begin{align}\alpha_k&=-\frac {r_k^Tp_k}{p_k^TAp_k} \\ &=-\frac {r_k^T(-r_k+\beta_kp_{k-1})}{p_k^TAp_k} \\ &=\frac {r_k^Tr_k}{p_k^TAp_k}\end{align}

改进思路2：由于rk+1=rk+αkApkr_{k+1}=r_k+\alpha_kAp_k可以得到Apk=rk+1−rkαkAp_k=\frac{r_{k+1}-r_k}{\alpha_k} βk+1=rTk+1ApkpTkApk=rTk+1ApkpTkApk=rTk+1(rk+1−rk)αkpTkApk=rTk+1rk+1rTkrk\begin{align}\beta_{k+1}&=\frac{r_{k+1}^TAp_{k}}{p_{k}^TAp_{k}} \\ &=\frac{r_{k+1}^TAp_{k}}{p_{k}^TAp_{k}}\\ &=\frac{r_{k+1}^T(r_{k+1}-r_k)}{\alpha_kp_{k}^TAp_{k}} \\ &=\frac{r_{k+1}^Tr_{k+1}}{r_{k}^Tr_{k}}\end{align}最后一步变换根据改进思路1和上述性质1。

因此最后的算法为：

1. 初始化x0x_0

2. r0=Ax0−b,p0=r0,k=0r_0=Ax_0-b,p_0=r_0,k=0

3. while rk≠0r_k \ne 0

αk=rTkrkpTkApk\alpha_k=\frac {r_k^Tr_k}{p_k^TAp_k}

xk+1=xk+αkpkx_{k+1}=x_k + \alpha_kp_k

rk+1=Axk+1−br_{k+1}=Ax_{k+1}-b

βk+1=rTk+1rk+1rTkrk\beta_{k+1}=\frac{r_{k+1}^Tr_{k+1}}{r_{k}^Tr_{k}}

pk+1=−rk+1+βk+1pkp_{k+1}=-r_{k+1}+\beta_{k+1}p_k

k=k+1k=k+1

4. end while

CG算法收敛性质

如果矩阵A有r个不同的特征值，则最多循环r步。

算法的收敛速度和A的特征值分布相关，如果A的条件数越大收敛越慢。

根据收敛性质，在实际应用中可以先对A进行预处理，从而使得A特征值分布更均匀一些。

非线性共轭梯度算法

将共轭的思想应用到一般化的最优化问题中，甚至非线性问题中去。

Fletcher-Reeves方法

算法如下：


主要改变如上图红框所示

1. 计算步长αk\alpha_k不一定能找到最优步长，可以采用wolf 条件进行计算。

2. 残差的计算可以直接用梯度代替。

3. 重点是需要保证每一步的搜索方向都满足下降，否则不保证收敛。选择强wolf条件可以保证上述约束

FR其他变种

主要思路优化βk\beta_k的计算。

总结

共轭梯度算法可以认为解决了一类特殊最优化问题:minϕ(x)=12xTAx−bTxmin\phi(x)=\frac12x^TAx-b^Tx，并且矩阵A对称且正定。

解决思路和lineSearch比较相似，首先确定搜索方向然后计算步长。CG特别的地方在于搜索方向相对于A共轭，并且步长均选择为最优步长。其收敛速度线性。通过本节内容需要了解

1. CG算法解决问题的形式

2. 共轭性

3. 共轭方向算法

4. 共轭梯度算法以及正确性证明

5. 收敛性