Numerical Optimization Ch5. Conjugate Gradient Methods

第五章: 共轭梯度法

共轭梯度法的作用主要体现在两个方面:

• 它们可高效求解大规模的线性方程组;
• 经调整, 它们可用于求解非线性优化问题.

线性共轭梯度法最早由Hestenes和Steifel在上世纪五十年代提出, 用于迭代求解带正定系数矩阵的线性系统. 由于它在大规模问题上的优越性, 它很快受到了广泛的欢迎. 而第一个非线性共轭梯度法则由Fletcher和Reeves于上市就六十年代提出. 而这也是最早的用于求解大规模非线性优化问题的算法之一. 迄今, 已出现了共轭梯度法的许多变体. 而它们所共同的关键特征, 就是无需矩阵存储且要快于最速下降法. 这为处理大规模问题提供了方便.

1. 线性共轭梯度法

线性共轭梯度法主要用于求解下面的线性系统: $$Ax=b,$$其中$A$为$n\times n$对称正定矩阵. 该问题可等价地描述成泛函极值问题, 即$$\min \varphi (x)\stackrel{def}{=}\frac{1}{2}{x}^{T}Ax-b^{T}x.$$ 因此线性共轭梯度法即可用于求解线性系统, 也可用于求解一类凸二次规划问题. 注意$\varphi$的梯度等于线性系统的残差, 即$$\nabla \varphi (x)=Ax-b\stackrel{def}{=}r(x),$$对于某一特定$x=x_k$, 简记$$r_k=A{x}_k-b.$$

1.1 共轭方向法

为了获取一些广而泛的特征与直观, 我们先来讨论共轭方向法. 实际上线性共轭梯度法是共轭方向法的一种特殊情形. 所谓"共轭", 此处指的是对于对称正定矩阵$A$的共轭, 即对于一组非零向量 { p 0 , p 1 , … , p l } \{p_0,p_1,\ldots,p_l\} {p0​,p1​,…,pl​}, 有$$p_i^{T}A{p}_j=0,\forall i̸=j.$$ 共轭方向法的好处在于, 它可在$n$步内达到$\varphi$的最小. 给定初始$x_0\in {\mathbb{R}}^n$, 并先不考虑生成方案地给定一组共轭方向 { p 0 , p 1 , … , p n − 1 } \{p_0,p_1,\ldots,p_{n-1}\} {p0​,p1​,…,pn−1​}, 则有生成序列$$x_{k+1}=x_k+{\alpha }_k{p}_k,$$其中${\alpha }_k$为二次函数$\varphi$沿着$x_k+\alpha p_k$的一维最小点, 我们可以解析地算出其表达式: $$\frac{\mathrm{d}}{\mathrm{d}\alpha }\varphi (x_k+\alpha p_k)=\nabla \varphi (x_k+\alpha p_k{)}^T{p}_k=(r_k+\alpha A{p}_k{)}^T{p}_k\stackrel{set}{=}0$$$$\Rightarrow {\alpha }_k=-\frac{r_k^T{p}_k}{p_k^{T}A{p}_k}.$$ (进一步验算$\frac{{\mathrm{d}}^2}{{\mathrm{d}}^2\alpha }\varphi (x_k+\alpha p_k)=p_k^{T}A{p}_k>0$, 因此为最小点.)

定理1 对于任意$x_0\in {\mathbb{R}}^n$, 由共轭方向产生的序列 { x k } \{x_k\} {xk​}经至多$n$步收敛到线性系统的解$x^{\ast }$.
证明: 注意到共轭方向 { p i } \{p_i\} {pi​}必定是线性无关的(这与$A$是正定的有直接的关系), 因此它们张成${\mathbb{R}}^n$. 从而存在坐标$({\sigma }_0,{\sigma }_1,\dots ,{\sigma }_{n-1})$, $$x^{\ast }-x_0={\sigma }_0{p}_0+{\sigma }_1{p}_1+\cdots +{\sigma }_{n-1}{p}_{n-1}.$$$\forall k$, 在上式两边左乘$p_k^{T}A$, 并利用共轭性质, 得到$${\sigma }_k=\frac{p_k^{T}A(x^{\ast }-x_0)}{p_k^{T}A{p}_k}.$$下面证明这些${\sigma }_k$恰好就是步长${\alpha }_k$. 由迭代公式, $$x_k=x_0+{\alpha }_0{p}_0+{\alpha }_1{p}_1+\cdots +{\alpha }_{k-1}{p}_{k-1}.$$同样两边左乘$p_k^{T}A$得到$$p_k^{T}A(x_k-x_0)=0.$$因此$$p_k^{T}A(x^{\ast }-x_0)=p_k^{T}A(x^{\ast }-x_k)=p_k^T(b-A{x}_k)=-p_k^T{r}_k.$$从而$${\alpha }_k=-\frac{r_k^T{p}_K}{p_k^{T}A{p}_k}=\frac{p_k^{T}A(x^{\ast }-x_0)}{p_k^{T}A{p}_k}={\sigma }_k.$$证毕.

对于共轭方向法有一种简单的直观解释: 设$A$是对角阵, 从而$\varphi$的等高线为椭球, 其轴向与坐标轴平行. 令 { p i } \{p_i\} {pi​}为标准正交基, 从而有图示如下 (以$n=2$为例):

图中$x^{\ast }$经过两步沿着坐标轴的搜索即可达到. 不过可以想象, 当$A$不再是对角阵时, 情形应当会有所不同. 图示如下:

图中$\varphi$的等高线依然是椭圆, 但轴向不再与坐标轴平行. 此时再沿着标准正交基构成的方向集搜索, 将不再能保证在$n$步以内 (甚至有限步内) 达到$x^{\ast }$. 这是因为 { e i } \{e_i\} {ei​}不再是其共轭方向. 不过我们可以将$A$转化成对角阵后再做搜索. 具体说, 由于$A$是对称正定矩阵, 所以可以正交对角化, 即存在正交矩阵$S$, 使得$S^{T}AS$为对角阵$\Lambda$. 作变量代换$y=S^{-1}x$, 则对应的泛函极值问题化为$$\min \hat{\varphi }(y)\stackrel{def}{=}\frac{1}{2}{y}^T\Lambda y-(S^{T}b{)}^{T}y.$$对应于变换后的坐标系, { p i = S e i } \{p_i=Se_i\} {pi​=Sei​}就是原坐标系下的对于$A$的共轭方向 (其实就是$S$的$n$列). 我们也可以在新的坐标系下得到$y^{\ast }$, 再令$x^{\ast }=S{y}^{\ast }$得到原问题的解.
在上面提到的对角阵的情形中, 我们可以发现: 每一次迭代都将确定$x^{\ast }$的一个分量. 换句话说, 再做完$k$次一维最小化后, $\varphi$已经在 s p a n { e 1 , e 2 , … , e k } \mathrm{span}\{e_1,e_2,\ldots,e_k\} span{e1​,e2​,…,ek​}上达到最小了. 下面我们将对任一对称正定矩阵证明这一性质. 为方便说明, 我们指出$r_k$也有类似于$x_k$的更新公式:$$r_{k+1}=A{x}_{k+1}-b=A{x}_k-b+A(x_{k+1}-x_k)=r_k+{\alpha }_{k}A{p}_k.$$

定理2 (扩张子空间最小化) 设$x_0\in {\mathbb{R}}^n$为任意选定的初始点, { x k } \{x_k\} {xk​}为生成的迭代序列. 则$$r_k^T{p}_i=0,i=0,1,\dots ,k-1,$$且$x_k$为$\varphi (x)$在 S k = d e f { x ∣ x = x 0 + s p a n { p 0 , p 1 , … , p k − 1 } } S_k\xlongequal{def}\{x|x=x_0+\mathrm{span}\{p_0,p_1,\ldots,p_{k-1}\}\} Sk​def {x∣x=x0​+span{p0​,p1​,…,pk−1​}}上的极小点.

证明: 我们先证明$\tilde{x}$在$S_k$上极小化$\varphi$当且仅当$r(\tilde{x}{)}^T{p}_i=0,i=0,1,\dots ,k-1.$ 定义$h(\sigma )=\varphi (x_0+{\sigma }_0{p}_0+\cdots +{\sigma }_{k-1}{p}_{k-1})$, 其中$\sigma =({\sigma }_0,{\sigma }_1,\dots ,{\sigma }_{k-1}{)}^T$. 由于$h(\sigma )$是严格凸的二次函数, 所有它有唯一极小点${\sigma }^{\ast }$, 且满足一阶稳定条件$${\frac{\partial h(\sigma )}{\partial {\sigma }_i}|}_{\sigma ={\sigma }^{\ast }}=0,i=0,1,\dots ,k-1$$$$\Rightarrow \nabla \varphi (x_0+{\sigma }_0^{\ast }{p}_0+\cdots +{\sigma }_{k-1}^{\ast }{p}_{k-1}{)}^T{p}_i=0,i=1,2,\dots ,k-1$$$$\Rightarrow r(\tilde{x}{)}^T{p}_i=0,i=1,2,\dots ,k-1.$$下面归纳证明$r_k^T{p}_i=0$, 从而有$x_k$为$\varphi (x)$在$S_k$上的极小点.

• 对于$k=1$, 我们由$x_1=x_0+{\alpha }_0{p}_0$为$\varphi$沿 { x 0 + α p 0 } \{x_0+\alpha p_0\} {x0​+αp0​}的极小得到$r_1^T{p}_0=0$.
• 假设$r_{k-1}^T{p}_i=0,i=0,1,\dots ,k-2$. 从而$$p_{k-1}^T{r}_k=p_{k-1}^T{r}_{k-1}+{\alpha }_{k-1}{p}_{k-1}^{T}A{p}_{k-1}\stackrel{{\alpha }_{k-1}的定义}{=}0.$$而对于$p_i,i=0,1,\dots ,k-2$, 由归纳假设以及共轭性质有$$p_i^T{r}_k=p_i^T{r}_{k-1}+{\alpha }_{k-1}{p}_i^{T}A{p}_{k-1}=0.$$从而$r_k^T{p}_i=0,i=0,1,\dots ,k-1$, 证毕.

这个定理指出: 当前的残差$r_k$与之前所有的搜索方向都正交. 这一性质在本章将被大量应用.
共轭方向法的优良性质不言而喻, 但随之而来的问题是: 这样好的共轭方向如何产生 (以及储存). 一种方法是之前提到的正交相似, 其中得到的正交矩阵的列向量组即是共轭方向. 然而对于大规模问题, 要得到所有的特征向量耗费巨大, 而且储存也成问题 (当然我们可以算出一个扔掉一个). 一种替代的方法是修正Gram-Schmidt正交化从而产生相互共轭的方向集, 而这样储存就成了主要的问题, 因为我们势必要存下全部的方向集. 基于以上, 所谓的线性共轭梯度法应运而生.

1.2 线性共轭梯度法的基本性质

线性共轭梯度法是特殊的共轭方向法, 只不过在产生其共轭方向时采用了特殊的方法: 在产生$p_k$时只需要用到$p_{k-1}$. 这一点已经说明了其在计算与存储上的巨大优势 (实际计算上还具有其他优势, 下面会提到).
具体说来, $p_k$被假定为负残差$-r_k$ (也即$\varphi$在$x_k$处的最速下降方向)和前一个方向$p_{k-1}$的线性组合:$$p_k=-r_k+{\beta }_k{p}_{k-1},$$其中${\beta }_k$被选定为用于满足共轭性质. 上式两端左乘$p_{k-1}^{T}A$并利用共轭性质得到$${\beta }_k=\frac{r_k^{T}A{p}_{k-1}}{p_{k-1}^{T}A{p}_{k-1}}.$$当然这里仅保证了相邻两方向的共轭性, 而完整的共轭性将在后面给出证明. 从这里我们会发现"共轭梯度法"实际上是一种误称: 不是梯度相互共轭, 而是搜索方向相互共轭. 注意$p_0$选定为$x_0$处的最速下降方向. ${\alpha }_k$仍采用一维搜索的方式解析确定. 从而得到如下算法:

算法1 (CG-Preliminary Version)
Given $x_0$;
Set $r_0\leftarrow A{x}_0-b,p_0\leftarrow -r_0,k\leftarrow 0$;
while $r_k̸=0$
${\alpha }_k\leftarrow -\frac{r_k^T{p}_k}{p_k^{T}A{p}_k}$;
$x_{k+1}\leftarrow x_k+{\alpha }_k{p}_k$;
$r_{k+1}\leftarrow A{x}_{k+1}-b$;
${\beta }_{k+1}\leftarrow \frac{r_{k+1}^{T}A{p}_k}{p_k^{T}A{p}_k}$;
$p_{k+1}\leftarrow -r_{k+1}+{\beta }_{k+1}{p}_k$;
$k\leftarrow k+1$;
end (while)

既然是Preliminary Version, 就说明我们后面会给出更加简洁且易于实施的版本. 下面的定理将说明以下问题:

• $p_0,p_1,\dots ,p_{n-1}$相互共轭, 从而CG能在至多$n$步内终止;
• 残差$r_i$与部分先前的$p_k$共轭;
• 残差$r_i$相互正交;
• 搜索方向$p_k$和残差$r_k$都包含在$r_0$的$k$阶Krylov子空间, 其定义为 K ( r 0 ; k ) = d e f s p a n { r 0 , A r 0 , … , A k r 0 } . \mathscr{K}(r_0;k)\xlongequal{def}\mathrm{span}\{r_0,Ar_0,\ldots,A^kr_0\}. K(r0​;k)def span{r0​,Ar0​,…,Akr0​}.

定理3 设共轭梯度法生成的第$k$个迭代项不是$x^{\ast }$. 则以下性质成立:

1. $r_k^T{r}_i=0,i=0,1,\dots ,k-1$;
2. s p a n { r 0 , r 1 , … , r k } = K ( r 0 ; k ) \mathrm{span}\{r_0,r_1,\ldots,r_k\}=\mathscr{K}(r_0;k) span{r0​,r1​,…,rk​}=K(r0​;k);
3. s p a n { p 0 , p 1 , … , p k } = K ( r 0 ; k ) \mathrm{span}\{p_0,p_1,\ldots,p_k\}=\mathscr{K}(r_0;k) span{p0​,p1​,…,pk​}=K(r0​;k);
4. $p_k^{T}A{p}_i=0,i=0,1,\dots ,k-1$;
5. $r_k^{T}A{p}_i=0,i=0,1,\dots ,k-2.$

从而序列 { x k } \{x_k\} {xk​}经至多$n$步收敛到$x^{\ast }$.

证明: 用数学归纳法证明第2,3,4条, 第5条为证明过程中的副产品. 其中第2,3条对于$k=0$显然成立, 第4条对于$k=1$显然成立. 假设这三条对于某个$k$仍然成立, 下面证明$k+1$的情形.

• 第2条的证明. 先证明左包含于右. 由归纳假设, $$r_k\in \mathcal{K}(r_0;k),p_k\in \mathcal{K}(r_0;k).$$在第二个式子两边左乘$A$得到 A p k ∈ s p a n { A r 0 , … , A k + 1 r 0 } . Ap_k\in\mathrm{span}\{Ar_0,\ldots,A^{k+1}r_0\}. Apk​∈span{Ar0​,…,Ak+1r0​}.因为$r_{k+1}=r_k+{\alpha }_{k}A{p}_k$, 所以 r k + 1 ∈ K ( r 0 ; k + 1 ) ⇒ s p a n { r 0 , r 1 , … , r k , r k + 1 } ⊂ K ( r 0 ; k + 1 ) . r_{k+1}\in\mathscr{K}(r_0;k+1)\Rightarrow \mathrm{span}\{r_0,r_1,\ldots,r_k,r_{k+1}\}\subset\mathscr{K}(r_0;k+1). rk+1​∈K(r0​;k+1)⇒span{r0​,r1​,…,rk​,rk+1​}⊂K(r0​;k+1).为证明相反方向, 利用第3条的归纳假设, 得到 A k + 1 r 0 = A ( A k r 0 ) ∈ s p a n { A p 0 , A p 1 , … , A p k } . A^{k+1}r_0=A(A^kr_0)\in\mathrm{span}\{Ap_0,Ap_1,\ldots,Ap_k\}. Ak+1r0​=A(Akr0​)∈span{Ap0​,Ap1​,…,Apk​}.由残差的更新公式我们可以反过来得到$A{p}_i=(r_{i+1}-r_i)/{\alpha }_i,i=0,1,\dots ,k$, 因此 A k + 1 r 0 ∈ s p a n { r 0 , r 1 , … , r k + 1 } ⇒ K ( r 0 ; k + 1 ) ⊂ s p a n { r 0 , r 1 , … , r k + 1 } . A^{k+1}r_0\in\mathrm{span}\{r_0,r_1,\ldots,r_{k+1}\}\Rightarrow\mathscr{K}(r_0;k+1)\subset\mathrm{span}\{r_0,r_1,\ldots,r_{k+1}\}. Ak+1r0​∈span{r0​,r1​,…,rk+1​}⇒K(r0​;k+1)⊂span{r0​,r1​,…,rk+1​}.因此第二条对于$k+1$仍然成立.
• 第3条的证明. s p a n { p 0 , p 1 , … , p k , p k + 1 } = s p a n { p 0 , p 1 , … , p k , r k + 1 } ( ∵ p k + 1 = − r k + 1 + β k + 1 p k ) = s p a n { r 0 , A r 0 , … , A k r 0 , r k + 1 ) ( 由 归 纳 假 设 ) = s p a n { r 0 , r 1 , … , r k , r k + 1 ) ( 由 第 2 条 ) = K ( r 0 ; k + 1 ) . \begin{aligned} &\mathrm{span}\{p_0,p_1,\ldots,p_k,p_{k+1}\}\\ &=\mathrm{span}\{p_0,p_1,\ldots,p_k,r_{k+1}\}(\because p_{k+1}=-r_{k+1}+\beta_{k+1}p_k)\\ &=\mathrm{span}\{r_0,Ar_0,\ldots,A^kr_0,r_{k+1})(由归纳假设)\\ &=\mathrm{span}\{r_0,r_1,\ldots,r_k,r_{k+1})(由第2条)\\ &=\mathscr{K}(r_0;k+1).\end{aligned} ​span{p0​,p1​,…,pk​,pk+1​}=span{p0​,p1​,…,pk​,rk+1​}(∵pk+1​=−rk+1​+βk+1​pk​)=span{r0​,Ar0​,…,Akr0​,rk+1​)(由归纳假设)=span{r0​,r1​,…,rk​,rk+1​)(由第2条)=K(r0​;k+1).​
• 第4条的证明. 在$$p_{k+1}=-r_{k+1}+{\beta }_{k+1}{p}_k$$两端右乘$A{p}_i,i=0,1,\dots ,k$, 得到$$p_{k+1}^{T}A{p}_i=-r_{k+1}^{T}A{p}_i+{\beta }_{k+1}{p}_k^{T}A{p}_i.$$由${\beta }_k$的定义, 上式右端当$i=k$时为0. 而对于$i\le k-1$则需要另做讨论. 由归纳假设, $p_0,p_1,\dots ,p_k$相互共轭, 因此上式右端第二项为0, 只需证明第一项为0 (此即第5条). 由定理2推出$$r_{k+1}^T{p}_i=0,i=0,1,\dots ,k.$$反复利用第3条, 我们得到对于$i=0,1,\dots ,k-1$, 以下成立: A p i ∈ A K ( r 0 ; i ) = s p a n { A r 0 , A 2 r 0 , … , A i + 1 r 0 } ⊂ s p a n { p 0 , p 1 , … , p i + 1 } . Ap_i\in A\mathscr{K}(r_0;i)=\mathrm{span}\{Ar_0,A^2r_0,\ldots,A^{i+1}r_0\}\subset\mathrm{span}\{p_0,p_1,\ldots,p_{i+1}\}. Api​∈AK(r0​;i)=span{Ar0​,A2r0​,…,Ai+1r0​}⊂span{p0​,p1​,…,pi+1​}.因此$$r_{k+1}^{T}A{p}_i=0,i=0,1,\dots ,k-1.$$从而第4条 (顺便第5条) 得证. 这也得出共轭梯度法求解线性系统至多$n$步终止.
• 第1条的证明. 这里我们不再使用数学归纳法. 从方向集的共轭性质以及定理2我们得到$r_k^T{p}_i=0,i=0,1,\dots ,k-1,k=1,2,\dots ,n-1$. 由残差的更新公式我们得到 r i ∈ s p a n { p i , p i − 1 } , i = 1 , … , k − 1 r_i\in\mathrm{span}\{p_i,p_{i-1}\},i=1,\ldots,k-1 ri​∈span{pi​,pi−1​},i=1,…,k−1. 因此$r_k^T{r}_i=0,i=1,\dots ,k-1$. 而$i=0$的情形是显然的: $r_k^T{r}_0=-r_k^T{p}_0=0$. 这就完成了全部的证明.

值得注意的是, 定理的证明要依赖于$p_0$的选取, 即至少要与初始负梯度$r_0$同方向, 否则定理的第3条对$k=0$就不成立了, 从而会影响第4条的证明, 进而"共轭"不再共轭.

1.3 实用共轭梯度法

我们可以利用定理2和定理3的结论推出更加经济高效的共轭梯度法. 这主要通过改变一些量的计算方法得到.

• ${\alpha }_k$的计算: 利用$p_{k+1}=-r_{k+1}+{\beta }_{k+1}{p}_k$以及残差与搜索方向的正交性, 得到$${\alpha }_k=\frac{p_k^T(r_{k+1}-r_k)}{p_k^{T}A{p}_k}=\frac{-p_k^T{r}_k}{p_k^{T}A{p}_k}=\frac{r_k^T{r}_k}{p_k^{T}A{p}_k}.$$
• ${\beta }_{k+1}$的计算: 由残差更新公式$r_{k+1}=r_k+{\alpha }_{k}A{p}_k$以及上面用到的搜索方向更新公式、正交性, 得到$${\beta }_{k+1}=\frac{r_{k+1}^{T}A{p}_k}{p_k^{T}A{p}_k}=\frac{r_{k+1}^T(r_{k+1}-r_k)/{\alpha }_k}{p_k^T(r_{k+1}-r_k)/{\alpha }_k}=\frac{r_{k+1}^T{r}_{k+1}}{r_k^T{r}_k}.$$

利用以上就可以得到实用版本的共轭梯度法.

算法2 (CG)
Given $x_0$;
Set $r_0\leftarrow A{x}_0-b,p_0\leftarrow -r_0,k\leftarrow 0$;
while $r_k̸=0$
${\alpha }_k\leftarrow \frac{r_k^T{r}_k}{p_k^{T}A{p}_k}$;
$x_{k+1}\leftarrow x_k+{\alpha }_k{p}_k$;
$r_{k+1}\leftarrow r_k+{\alpha }_{k}A{p}_k$;
${\beta }_{k+1}\leftarrow \frac{r_{k+1}^T{r}_{k+1}}{r_k^T{r}_k}$;
$p_{k+1}\leftarrow -r_{k+1}+{\beta }_{k+1}{p}_k$;
$k\leftarrow k+1$;
end (while)

对于算法2, 我们在每一步都只需要$x,r,p$的至多两项. 因此 (正如前述), 共轭梯度法对于存储量的需求是很小的. 而主要的计算量就落在矩阵-向量乘积$A{p}_k,p_k^{T}A{p}_k,r_{k+1}^T{r}_{k+1}$上, 后二者都可在$O(n)$计算量内解决, 而前者则需$O(n^2)$.
注意, 共轭梯度法只在大规模问题中具有优越性, 否则诸如Gauss消去或者其他的如奇异值分解的矩阵分解方法将更加有效, 这主要是因为它们对于舍入误差更加强健. 而对于大型问题, 共轭梯度法还有一个优点是不会改变系数矩阵以及一定程度上保证矩阵的稀疏性. 另一个优点就得从收敛速度上说了.

1.4 收敛速度

我们知道, 在无舍入误差的计算下共轭梯度法将在至多$n$步终止. 不过更加有意思的是, 当系数矩阵$A$的特征值分布具有一定的特点时, 算法往往仅需远少于$n$步终止. 我们先来以一种稍微不同的角度审视定理2, 从而说明算法2在某种程度上是最优的.
由迭代公式和定理3中的第3条可知, 存在${\gamma }_i,i=0,1,\dots ,k$使得$$\begin{array}{rl}{x}_{k+1}& =x_0+{\alpha }_0{p}_0+\cdots +{\alpha }_k{p}_k\\ & =x_0+{\gamma }_0{r}_0+{\gamma }_{1}A{r}_0+\cdots +{\gamma }_k{A}^k{r}_0.\end{array}$$定义$P_k^{\ast }(\cdot )$为带系数${\gamma }_0,{\gamma }_1,\dots ,{\gamma }_k$的$k$次多项式, 从而$$x_{k+1}=x_0+P_k^{\ast }(A)r_0.$$下面我们证明, 在Krylov空间$\mathcal{K}(r_0;k)$内寻得$k$步的所有算法中, 算法2表现最好, 即: $k$步之后算法2达到的迭代项离$x^{\ast }$是最近的. 这里的距离为加权范数诱导的距离, $$\Vert z{\Vert }_A^2=z^{T}Az.$$(回忆我们在最速下降法的分析中也用到了这种带权的距离, 我们还会在后面比较共轭梯度法与最速下降法.) 使用这一范数, 我们可得到$$\frac{1}{2}\Vert x-x^{\ast }{\Vert }_A^2=\frac{1}{2}(x-x^{\ast }{)}^{T}A(x-x^{\ast })=\varphi (x)-\varphi (x^{\ast }).$$也就是说自变量的带权距离反映了函数值的差.
定理2是说, $x_{k+1}$在$S_k$上极小化$\varphi$, 也就是$\Vert x-x^{\ast }{\Vert }_A^2$. 换句话说, 就是$P_k^{\ast }$是如下问题的解:$$\underset{P_k}{\min }\Vert x_0+P_k(A)r_0-x^{\ast }{\Vert }_A.$$我们对上述问题做一些变换. 由于$$r_0=A{x}_0-b=A{x}_0-A{x}^{\ast }=A(x_0-x^{\ast }),$$所以$$x_{k+1}-x^{\ast }=x_0+P_k^{\ast }(A)r_0-x^{\ast }=[I+P_k^{\ast }(A)A](x_0-x^{\ast }).$$设$0<{\lambda }_1\le {\lambda }_2\le \cdots \le {\lambda }_n$为$A$的特征值, $v_1,v_2,\dots ,v_n$为对应的相互标准正交的特征向量, 因此有$A$的谱分解$$A=\sum _{i=1}^n{\lambda }_i{v}_i{v}_i^T.$$因为特征向量全体张成整个${\mathbb{R}}^n$, 所以存在${\xi }_i,i=1,\dots ,n$使得$$x_0-x^{\ast }=\sum _{i=1}^n{\xi }_i{v}_i.$$从而$$x_{k+1}-x^{\ast }=\sum _{i=1}^n[1+{\lambda }_i{P}_k^{\ast }({\lambda }_i)]{\xi }_i{v}_i.$$由于$\Vert z{\Vert }_A^2=z^{T}Az={\sum }_{i=1}^n{\lambda }_i(v_i^{T}z{)}^2$, 所以$$\Vert x_{k+1}-x^{\ast }{\Vert }_A^2=\sum _{i=1}^n{\lambda }_i[1+{\lambda }_i{P}_k^{\ast }({\lambda }_i){]}^2{\xi }_i^2.$$由于由共轭梯度法诱导的多项式$P_k^{\ast }$在此范数下是最优的, 因此有$$\Vert x_{k+1}-x^{\ast }{\Vert }_A^2=\underset{P_k}{\min }\sum _{i=1}^n{\lambda }_i[1+{\lambda }_i{P}_k({\lambda }_i){]}^2{\xi }_i^2.$$提出最大的$[1+{\lambda }_i{P}_k({\lambda }_i){]}^2$, 我们有$$\begin{array}{rl}\Vert x_{k+1}-x^{\ast }{\Vert }_A^2& \le \underset{P_k}{\min }\underset{1\le i\le n}{\max }[1+{\lambda }_i{P}_k({\lambda }_i){]}^2\left(\sum _{j=1}^n{\lambda }_j{\xi }_j^2\right)\\ & =\underset{P_k}{\min }\underset{1\le i\le n}{\max }[1+{\lambda }_i{P}_k({\lambda }_i){]}^2\Vert x_0-x^{\ast }{\Vert }_A^2.\end{array}$$上式使我们可通过估计$$\underset{P_k}{\min }\underset{1\le i\le n}{\max }[1+{\lambda }_i{P}_k({\lambda }_i){]}^2$$来量化共轭梯度法的收敛速度. 换句话说, 我们可以找到一个$P_k$ (不见得是$P_k^{\ast }$) 使得这一表达式尽可能地小. 在某些特殊情形下, 我们可以解析地找到这个多项式并得到关于共轭梯度法的一些有趣的性质. 下面就是一个例子.

定理4 若$A$仅有$r$个不同的特征值, 则共轭梯度法至多$r$步终止.
证明: 假设${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n$中有$r$个不同的特征值${\tau }_1<{\tau }_2<\cdots <{\tau }_r$. 定义多项式$Q_r(\lambda )$为$$Q_r(\lambda )=\frac{(-1{)}^r}{{\tau }_1{\tau }_2\cdots {\tau }_r}(\lambda -{\tau }_1)(\lambda -{\tau }_2)\cdots (\lambda -{\tau }_r).$$注意$Q_r({\lambda }_i)=0,i=1,2,\dots ,n,Q_r(0)=1$. 从而$Q_r(\lambda )-1$为$r$次多项式且由一根$\lambda =0$. 由此定义$r-1$次多项式${\bar{P}}_{r-1}$为$${\bar{P}}_{r-1}=\frac{Q_r(\lambda )-1}{\lambda }.$$在$$\underset{P_k}{\min }\underset{1\le i\le n}{\max }[1+{\lambda }_i{P}_k({\lambda }_i){]}^2$$中令$k=r-1$我们有$$0\le \underset{P_{r-1}}{\min }\underset{1\le i\le n}{\max }[1+{\lambda }_i{P}_{r-1}({\lambda }_i){]}^2\le \underset{1\le i\le n}{\max }[1+{\lambda }_i{\bar{P}}_{r-1}({\lambda }_i){]}^2=\underset{1\le i\le n}{\max }{Q}_r^2({\lambda }_i)=0.$$因此对于$k=r-1$, $$\underset{P_k}{\min }\underset{1\le i\le n}{\max }[1+{\lambda }_i{P}_k({\lambda }_i){]}^2=0.$$所以$\Vert x_r-x^{\ast }{\Vert }_A^2=0\Rightarrow x_r=x^{\ast }$. 证毕.

更有甚者, 我们有以下估计:
定理5 若$A$有特征值${\lambda }_1\le {\lambda }_2\le \cdots \le {\lambda }_n$, 则我们有$$\Vert x_{k+1}-x^{\ast }{\Vert }_A^2\le {\left(\frac{{\lambda }_{n-k}-{\lambda }_1}{{\lambda }_{n-k}+{\lambda }_1}\right)}^2\Vert x_0-x^{\ast }{\Vert }_A^2.$$
证明 (Sketch): 选取多项式${\bar{P}}_k$使得多项式$Q_{k+1}(\lambda )=1+\lambda {\bar{P}}_k(\lambda )$以$k$个最大的特征值${\lambda }_n,{\lambda }_{n-1},\dots ,{\lambda }_{n-k+1}$和${\lambda }_1$与${\lambda }_{n-k}$的中点为根. 可以证明$Q_{k+1}$在余下的特征值${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_{n-k}$上的最大值就是$({\lambda }_{n-k}-{\lambda }_1)/({\lambda }_{n-k}+{\lambda }_1)$.

比起定理5的详细证明, 其对于共轭梯度法在一些问题上表现的解释更加吸引人. 假设特征值如下图一般分布:

其中$A$有$m$个较大的特征值和$n-m$个较小的在1附近的特征值. 定义$ϵ={\lambda }_{n-m}-{\lambda }_1$, 则定理5告诉我们在$m+1$步之后$$\Vert x_{m+1}-x^{\ast }{\Vert }_A\approx ϵ\Vert x_0-x^{\ast }{\Vert }_A.$$对于较小的$ϵ$, 我们有信心说共轭梯度法可在$m+1$步之后就提供解的较好估计.
我们也可以反过来用共轭梯度法的表现得到一些关于矩阵$A$特征值分布的信息. 比如下图:

其中实线表示的是$A$具有5个最大的特征值, 而剩下的较小特征值在0.95~1.05之间的情形, 虚线表示的是$A$特征值随机分布的情形. 对于实线情形 (这时我们也说$A$的谱呈簇状), 定理5会说误差对数会在第6步迭代得到一次显著的下降. 但实际上 (图中), 在第5步就已经下降了许多. 这也说明定理5只是给出了估计上界, 收敛速度可能会更快些. 进一步, 我们观察到下一步迭代又带来了误差对数的显著下降. 从定理4我们可以说, $A$几乎只有6个不同的特征值: 5个较大的特征值和1. 此时定理5就会说, 我们差不多可以在6步迭代后得到收敛. 不过实际上却是7步, 这是特征值在1附近的分散带来的效果. 相比之下, 不成簇的谱带来的收敛曲线是平缓均匀的.
另一方面, 我们可以得到基于$A$的条件数的共轭梯度法的收敛估计. 由2-范数诱导的2-条件数的性质, 我们有$$\kappa (A)=\Vert A{\Vert }_2\Vert A^{-1}{\Vert }_2={\lambda }_n/{\lambda }_1.$$可以证明$$\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.$$这由定理5可以轻松证出, 这里我们给出不一样的证明方法.

定理6 共轭梯度法的误差有估计$$\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.$$证明: 事实上, 从前面的推导中我们得到$$\begin{array}{rl}\Vert x_k-x^{\ast }{\Vert }_A^2& \le \underset{P_{k-1}}{\min }\underset{1\le i\le n}{\max }[1+{\lambda }_i{P}_{k-1}({\lambda }_i){]}^2\Vert x_0-x^{\ast }{\Vert }_A^2\\ & \le \underset{P_{k-1}}{\min }\underset{a\le \lambda \le b}{\max }[1+\lambda P_{k-1}(\lambda ){]}^2\Vert x_0-x^{\ast }{\Vert }_A^2,\end{array}$$其中$a={\lambda }_1,b={\lambda }_n$. 由Chebyshev多项式逼近定理知, 最优化问题$$\underset{P_{k-1}}{\min }\underset{a\le \lambda \le b}{\max }[1+\lambda P_{k-1}(\lambda ){]}^2$$具有唯一解$$1+\lambda {\bar{P}}_{k-1}(\lambda )=\frac{T_k\left(\frac{b+a-2\lambda }{b-a}\right)}{T_k\left(\frac{b+a}{b-a}\right)},$$其中$T_k(x)$是$k$次Chebyshev多项式. 由Chebyshev多项式的性质知$$\underset{a\le \lambda \le b}{\max }|1+\lambda {\bar{P}}_{k-1}(\lambda )|=\frac{1}{T_k\left(\frac{b+a}{b-a}\right)}\le 2{\left(\frac{\sqrt{\kappa (A)}-1}{\sqrt{\kappa (A)}+1}\right)}^k,$$于是就有$$\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$的2-条件数, 或具体说依赖于$A$最大特征值与最小特征值的比值. 若$A$本身的2-条件数足够接近1或者等价地说$A$的特征值较集聚, 则此定理充分地推出共轭梯度法的收敛速度将非常快. 否则, 我们应当在共轭梯度法实施之前对$A$做一些"手脚", 这就是下面要谈及的预处理问题.

这里再提一下共轭梯度法与最速下降法收敛速度之比较. 回忆第三章中最速下降法的收敛速度$$\Vert x_{k+1}-x^{\ast }{\Vert }_Q\le \left(\frac{{\lambda }_n-{\lambda }_1}{{\lambda }_n+{\lambda }_1}\right)\Vert x_k-x^{\ast }{\Vert }_Q=\left(\frac{\kappa (Q)-1}{\kappa (Q)+1}\right)\Vert x_k-x^{\ast }{\Vert }_Q.$$递推得到$$\Vert x_k-x^{\ast }{\Vert }_Q\le {\left(\frac{\kappa (Q)-1}{\kappa (Q)+1}\right)}^k\Vert x_0-x^{\ast }{\Vert }_Q.$$因为$\frac{\lambda -1}{\lambda +1}$在$\lambda >1$时是单调递增的, 因此共轭梯度法的 (渐进)收敛速度是优于最速下降法的. 注意这里共轭梯度法的误差上界还是保守估计, 而且共轭梯度法还具有计算和存储的双重优势. 可以说在求解线性系统上, 共轭梯度法要优于最速下降法.

1.5 预处理

1.4中我们谈到共轭梯度法的误差上界与系数矩阵$A$的2-条件数紧密相关, 具体说, 若$A$的2-条件数越接近于1, 则理论上共轭梯度法的收敛速度会更快; 极端情况下, $A$的2-条件数就是1, 此时共轭梯度法将一步收敛.
而对于一般的$A$, 我们不能先验地苛求其特征值具有何种特殊的分布. 相比之下, 后验地对其处理才更加易于实施. 受$A$的原有正交相似的启发, 我们可寻求可逆变换矩阵$C$使得$C^{-T}A{C}^{-1}$的特征值的分布要优于$A$, 从而我们只需求解以$C^{-T}A{C}^{-1}$为系数矩阵的问题, 再最后变换为原解即可. 具体说来, 设有可逆矩阵$C$, 做变量代换$y=Cx$, 则$\varphi$此时变成$$\hat{\varphi }(y)=\frac{1}{2}{y}^T(C^{-T}A{C}^{-1})y-(C^{-T}b{)}^{T}y.$$使用算法2极小化$\hat{\varphi }$等价于求解线性系统$$(C^{-T}A{C}^{-1})y=C^{-T}b.$$此时收敛速度就取决于$C^{-T}A{C}^{-1}\triangleq \tilde{A}$的特征值分布.
$C$的作用应当有:

• 使得$\tilde{A}$的条件数显著小于$A$的条件数; 或者
• 使得$\tilde{A}$的特征值更加集聚成簇, 从而共轭梯度法能在更少的步数下收敛.

1.5.1 算法2的预处理版本

我们不需要显式地作变换$y=Cx$. 相比之下, 将变换内嵌入算法2显得更加经济实惠 (否则需要矩阵-矩阵乘积). 我们也不显式地使用$C$, 而是设$M=C^{T}C$, 此矩阵也称作预处理子. 实际上若$C^{-T}A{C}^{-1}$的特征值都接近于1, 即有$$C^{-T}A{C}^{-1}\approx I\Rightarrow M=C^{T}C\approx A.$$此时$M$相当于是$A$的逆.

算法3 (Preconditioned CG)
Given $x_0$, preconditioner $M$;
Set $r_0\leftarrow A{x}_0-b$;
Solve $M{y}_0=r_0$ for $y_0$;
Set $p_0=-y_0,k\leftarrow 0$;
while $r_k̸=0$
${\alpha }_k\leftarrow \frac{r_k^T{y}_k}{p_k^{T}A{p}_k}$;
$x_{k+1}\leftarrow x_k+{\alpha }_k{p}_k$;
$r_{k+1}\leftarrow r_k+{\alpha }_{k}A{p}_k$;
$\mathrm{S}\mathrm{o}\mathrm{l}\mathrm{v}\mathrm{e}M{y}_{k+1}=r_{k+1}$;
${\beta }_{k+1}\leftarrow \frac{r_{k+1}^T{y}_{k+1}}{r_k^T{y}_k}$;
$p_{k+1}\leftarrow -y_{k+1}+{\beta }_{k+1}{p}_k$;
$k\leftarrow k+1$;
end (while)

若在算法3中令$M=I$, 则回到算法2. 此时残差的正交性变成了$$r_i^T{M}^{-1}{r}_j=0,i̸=j.$$就计算量而言, 算法3与算法2的主要区别在于是否需要求解$My=r$ (当然若有$C$的显式表达, 我们就可以设计更好的算法求解之).

1.5.2 实用预处理子

没有任何一种预处理的方式是普适最佳的. 我们在设计预处理的时候需要关注许多方面的问题, 例如

• $M$的有效性, 即是否是$A$的较好近似, 亦或能够使$\tilde{A}$的特征值聚集成簇;
• $M$的计算与存储;
• 线性系统$My=r$的求解.

往往这些目标不能同时达到, 此时我们需要权衡利弊 (tradeoff) 以获得更佳方案. 预处理子的设计没有统一的方案, 随着问题的不同而不同. 例如在数值求解偏微分方程时, 我们选取的$M$往往使得$My=r$为系统$Ax=b$较为粗糙的近似 (例如ILU预处理). 而在其他的领域, 问题的结构和起源则会成为设计预处理子的关键.
尽管如此, 还是有一些一般性的预处理子可供选择, 它们在各个问题的表现有所不同. 其中包括: 对称逐次超松弛(SSOR)预处理子、不完全Cholesky分解预处理子与带状预处理子. 一般不完全Cholesky分解预处理子是最高效的. 其想法很简单, 我们不精确地计算$A$的Cholesky分解$A=L{L}^T$, 而是计算近似的因子$\tilde{L}$ (比$L$更加稀疏), 使得$A\approx \tilde{L}{\tilde{L}}^T$. 令$C={\tilde{L}}^T$, 就有$M=\tilde{L}{\tilde{L}}^T$以及$$C^{-T}A{C}^{-1}={\tilde{L}}^{-1}A{\tilde{L}}^{-T}\approx I.$$此时具有显式$C$的表达, 且$C$还是上三角矩阵, 从而我们不需要在算法3中显式计算$M$; $My=r$的求解也可用回代-前代的方式经济地求解.
不完全Cholesky分解也有一定的缺陷. 例如, 其导出的$M$可能不是(充分)正定的, 因此需要我们增加一些防护措施 (safeguard). 再比如, 由于我们在$\tilde{L}$上强加了稀疏性的限制, 这可能会带来数值上的不稳定. 当然我们可以使用更稠密的$\tilde{L}$增加稳定性, 但这样一来计算量也就有所提升.

2. 非线性共轭梯度法

之前指出, 线性共轭梯度法可用于求解凸二次函数的规划问题. 那么我们自然要问: 该算法能否推广至求解一般的凸函数问题或者非线性问题. 本节我们就来讨论非线性情形下线性共轭梯度法的变体——非线性共轭梯度法——及其优良的性质.

2.1 Fletcher-Reeves(FR)方法

FR方法想法很直接, 其只在算法2上改动两处:

1. ${\alpha }_k$的计算: 先前我们能在凸二次函数的基础上解析地计算${\alpha }_k$的表达式. 但在一般情形下, 解析计算难以实施. 我们应当使用某种手段近似地获取$\varphi$在搜索方向上的最小值点;
2. 残差与梯度: 先前在凸二次函数上残差就等于梯度. 一般情形下, 我们需要额外估计梯度.

算法4 (FR)
Given $x_0$;
Evaluate $f_0=f(x_0),\nabla f_0=\nabla f(x_0)$;
Set $p_o\leftarrow -\nabla f_0,k\leftarrow 0$;
while $\nabla f_k̸=0$
\quad\quad Compute ${\alpha }_k$ and set $x_{k+1}=x_k+{\alpha }_k{p}_k$;
\quad\quad Evaluate $\nabla f_{k+1}$;
${\beta }_{k+1}^{\mathrm{F}\mathrm{R}}\leftarrow \frac{\nabla f_{k+1}^T\nabla f_{k+1}}{\nabla f_k^T\nabla f_k}$;
$p_{k+1}\leftarrow -\nabla f_{k+1}+{\beta }_{k+1}^{\mathrm{F}\mathrm{R}}{p}_k$;
$k\leftarrow k+1$;
end (while)

注意当$f$为严格凸二次函数时, 算法4就回到了算法2. 算法4在处理大型非线性优化时具有独到的优势:

• 每步迭代只需要函数值和梯度;
• 每步计算均不涉及矩阵运算;
• 只需存储少量的向量.

下面我们来完善算法4中的细节——${\alpha }_k$的选取. 若${\alpha }_k$选取不善, 则$p_k$就可能不是个下降方向. 对$$p_k=-\nabla f_k+{\beta }_k^{\mathrm{F}\mathrm{R}}{p}_{k-1}$$两边内积上向量$\nabla f_k$, 得到$$\nabla f_k^T{p}_k=-\Vert \nabla f_k{\Vert }^2+{\beta }_k^{\mathrm{F}\mathrm{R}}\nabla f_k^T{p}_{k-1}.$$

• 若对于${\alpha }_{k-1}$的线搜索精确, 则${\alpha }_{k-1}$就是$f$沿着方向$p_{k-1}$的局部极小点, 从而$$\nabla f_k^T{p}_{k-1}=0.$$此时有$\nabla f_k^T{p}_k<0$, 从而$p_k$确为下降方向.
• 若对于${\alpha }_{k-1}$的线搜索非精确, 则${\beta }_k^{\mathrm{F}\mathrm{R}}\nabla f_k^T{p}_{k-1}$就可能决定右端的正负, 进而可能会导致$\nabla f_k^T{p}_k>0$. 此时若在线搜索上加强Wolfe条件, 即$$\begin{array}{rl}f(x_k+{\alpha }_k{p}_k)& \le f(x_k)+c_1{\alpha }_k\nabla f_k^T{p}_k,\\ |\nabla f(x_k+{\alpha }_k{p}_k{)}^T{p}_k|& \le -c_2\nabla f_k^T{p}_k,\end{array}$$其中$0<c_1<c_2<\frac{1}{2}$ (这要比之前$c_2<1$要强), 利用我们接下来要证明的定理7, 就可以说$p_k$是下降方向.

2.2 Polak-Ribiere(PR)方法及其变体

FR方法具有许多变体, 它们的主要区别都在于${\beta }_k$的选取上.

• PR方法: Polak与Ribiere提出的参数为$${\beta }_{k+1}^{\mathrm{P}\mathrm{R}}=\frac{\nabla f_{k+1}^T(\nabla f_{k+1}-\nabla f_k)}{\Vert \nabla f_k{\Vert }^2}.$$注意到, 当$f$时强凸二次函数且线搜索精确时, 由于不同迭代步上的梯度是相互正交的 (因为残差相互正交), 所以${\beta }_{k+1}^{\mathrm{P}\mathrm{R}}={\beta }_{k+1}^{\mathrm{F}\mathrm{R}}$. 然而当用于一般的非线性函数且线搜索非精确时, FR与PR的表现差异就相当大了. 数值实验表示, PR方法一般要更加强健和高效.
  不过PR方法也有缺陷: 尽管加强Wolfe条件, 我们仍不能保证$p_k$是下降方向. 不过如果我们修正$\beta$为 β k + 1 + = max ⁡ { β k + 1 P R , 0 } , \beta_{k+1}^+=\max\{\beta_{k+1}^{\mathrm{PR}},0\}, βk+1+​=max{βk+1PR​,0},从而得到PR+方法, 此时只需稍稍改变强Wolfe条件即可产生下降方向 (后面我们会提到这相当于一种重启策略).
• HS方法: 还有许多其他的${\beta }_{k+1}$可供选择, 在目标函数强凸二次且线搜索精确时, 它们都与FR方法中的${\beta }_{k+1}^{\mathrm{F}\mathrm{R}}$相契合. 例如Hestenes-Stiefel(HS)方法, 其中$${\beta }_{k+1}^{\mathrm{H}\mathrm{S}}=\frac{\nabla f_{k+1}^T(\nabla f_{k+1}-\nabla f_k)}{(\nabla f_{k+1}-\nabla f_k{)}^T{p}_k}.$$HS方法的理论收敛性质和实际表现与PR方法都很相近. 实际上, HS方法的公式可通过要求相邻搜索方向对$[x_k,x_{k+1}]$上的平均Hessian相互共轭得到, 其中平均Hessian定义为$${\bar{G}}_k={\int }_0^1[{\nabla }^{2}f(x_k+\tau {\alpha }_k{p}_k)]\mathrm{d}\tau .$$由Taylor定理我们知道, $\nabla f_{k+1}=\nabla f_k+{\alpha }_k{\bar{G}}_k{p}_k$, 从而若对$p_{k+1}=-\nabla f_{k+1}+{\beta }_{k+1}{p}_k$要求$p_{k+1}^T{\bar{G}}_k{p}_k=0$就可以得到${\beta }_{k+1}^{\mathrm{H}\mathrm{S}}$.
• FR-PR方法: 后面我们会看到, 所有满足$$|{\beta }_k|\le {\beta }_k^{\mathrm{F}\mathrm{R}},\forall k\ge 2$$的${\beta }_k$都能保证全局收敛性. 这一命题启发我们在PR方法上做一些修正, 即在$k\ge 2$上令$${\beta }_k=\{\begin{array}{rl}-{\beta }_k^{\mathrm{F}\mathrm{R}}& \mathrm{i}\mathrm{f}{\beta }_k^{\mathrm{P}\mathrm{R}}<-{\beta }_k^{\mathrm{F}\mathrm{R}},\\ {\beta }_k^{\mathrm{P}\mathrm{R}}& \mathrm{i}\mathrm{f}|{\beta }_k^{\mathrm{P}\mathrm{R}}|\le {\beta }_k^{\mathrm{F}\mathrm{R}},\\ {\beta }_k^{\mathrm{F}\mathrm{R}}& \mathrm{i}\mathrm{f}{\beta }_k^{\mathrm{P}\mathrm{R}}>{\beta }_k^{\mathrm{F}\mathrm{R}}.\end{array}$$这样的${\beta }_k$在某些应用上表现很好, 我们称对应的方法为FR-PR方法.
• 其他: 以下二者具有良好的理论性质和数值表现. 而且仅需在线搜索上加Wolfe条件即可保证下降.
  1. ${\beta }_{k+1}=\frac{\Vert f_{k+1}{\Vert }^2}{(\nabla f_{k+1}-\nabla f_k{)}^T{p}_k}$;
  2. ${\beta }_{k+1}=\left({\hat{y}}_k-2p_k\frac{\Vert {\hat{y}}_k{\Vert }^2}{{\hat{y}}_k^T{p}_k}\right)\frac{\nabla f_{k+1}}{{\hat{y}}_k^T{p}_k},{\hat{y}}_k=\nabla f_{k+1}-\nabla f_k$.

2.3 二次终止性与重启

二次终止性是指, 当非线性共轭梯度法用于处理严格凸二次函数时, 算法退化为线性情形下的算法2, 从而同样能达到在至多$n$步内终止. 值得注意的是, 为了保证这一点, 我们需要在步长${\alpha }_k$的选取上加一些条件. 例如, 使用第三章中的插值方法获取${\alpha }_k$.
重启是非线性共轭梯度法的一种修正手段, 其操作即为每$n$步在$p_{k+1}=-\nabla f_{k+1}+{\beta }_{k+1}{p}_k$中令${\beta }_{k+1}=0$, 从而取负梯度为当前搜索方向 (从而一定是下降方向, 这一点在算法的进行会使得$p_k$逐渐与负梯度正交时显得尤为重要). 重启的作用在于, 定期地清理算法中的无益的过期信息. 我们甚至可以证明关于重启的一个较强的结论: 重启的共轭梯度法是$n$-步二次收敛的, 即$$\Vert x_{k+n}-x^{\ast }\Vert =O(\Vert x_k-x^{\ast }{\Vert }^2).$$这一点并不奇怪, 我们从以下两种情形说明:

• 当解附近$f$为强凸二次函数, 但在其他地方并不是二次函数. 假设算法产生的迭代点收敛于解, 从而最终会进入二次区域. 从某一点开始, 算法重启, 而往后就直接是在实施线性共轭梯度法. 特别地, 算法在$n$步内终止. 此时重启是重要的, 因为线性共轭梯度法的$n$步终止的前提, 是初始搜索方向为负梯度.
• 当解附近并不是二次函数. 由Taylor定理, 目标函数仍然可以用二次函数近似. 不过因此, 我们也不能再苛求$n$步收敛.

尽管$n$-步二次收敛在理论上很好, 但并不能用于实际操作. 这是因为, 共轭梯度法仅在大规模问题上应用良好 (中小规模的问题上, 共轭梯度法对舍入误差的敏感性会有影响. 这种问题上, 可能远没到$n$步算法就已经收敛, 来不及重启了). 因此, 非线性共轭梯度法在实际应用中要么不重启, 要么以不同于迭代计数的方式重启. 例如, 我们知道在$f$为二次函数时, 有梯度之间的正交关系. 因此规定, 当两相邻梯度"相当"不正交时, 即$$\frac{|\nabla f_k^T\nabla f_{k-1}|}{\Vert \nabla f_k{\Vert }^2}\ge v,$$我们就重启算法. 这里$v$的常用取值为0.1.
之前提到, PR+方法中的${\beta }_{k+1}$实际上就是一种重启策略. 回忆 β k + 1 + = max ⁡ { β k + 1 P R , 0 } , \beta_{k+1}^+=\max\{\beta_{k+1}^{\mathrm{PR}},0\}, βk+1+​=max{βk+1PR​,0},以及$${\beta }_{k+1}^{\mathrm{P}\mathrm{R}}=\frac{\nabla f_{k+1}^T(\nabla f_{k+1}-\nabla f_k)}{\Vert \nabla f_k{\Vert }^2}.$$PR+中, 只要${\beta }_{k+1}^{\mathrm{P}\mathrm{R}}<0$即取负梯度为下一个搜索方向, 实现重启. 而${\beta }_{k+1}^{\mathrm{P}\mathrm{R}}<0$当且仅当$\nabla f_{k+1}^T\nabla f_k$充分大 (至少大于$\Vert f_{k+1}{\Vert }^2$), 而这也就说明梯度之间的正交性被破坏了. 等价地就是上面$v$取1的情形. 因此相比于上面陈述的重启方案, PR+方法的重启更加不频繁.

2.4 FR方法的收敛性

2.4.1 FR方法的性质

下面这个定理说的是, 在一定条件下, FR方法产生的搜索方向都是下降方向.

定理7 假设算法4中的步长选取满足强Wolfe条件 (其中$0<c_2<\frac{1}{2}$) (只需假设水平集 L = { x : f ( x ) ≤ f ( x 0 ) } \mathcal{L}=\{x:f(x)\le f(x_0)\} L={x:f(x)≤f(x0​)}有界以及$f$二阶连续可微, 由第三章的定理即可得存在性), 则算法产生的下降方向$p_k$满足以下不等式:$$-\frac{1}{1-c_2}\le \frac{\nabla f_k^T{p}_k}{\Vert \nabla f_k{\Vert }^2}\le \frac{2c_2-1}{1-c_2},k=0,1,\dots .$$
证明: 首先注意函数$t(\xi )\stackrel{def}{=}(2\xi -1)/(1-\xi )$在区间$[0,1/2]$上单调递增, 且$t(0)=-1,t(1/2)=0$. 因此对于$c_2\in (0,1/2)$, 有$$-1<\frac{2c_2-1}{1-c_2}<0.$$因此若不等式成立, 则立得$p_k$是下降方向.
下面使用数学归纳法证明不等式.

1. $k=0$时, 不等式中间为-1, 从而成立;
2. 假设对于某个$k\ge 1$成立不等式. 由搜索方向公式和${\beta }_{k+1}$的计算公式, 我们有$$\frac{\nabla f_{k+1}^T{p}_{k+1}}{\Vert \nabla f_{k+1}{\Vert }^2}=-1+{\beta }_{k+1}\frac{\nabla f_{k+1}^T{p}_k}{\Vert \nabla f_{k+1}{\Vert }^2}=-1+\frac{\nabla f_{k+1}^T{p}_k}{\Vert \nabla f_k{\Vert }^2}.$$由强Wolfe条件的曲率条件, 我们有$$|\nabla f_{k+1}^T{p}_k|\le -c_2\nabla f_k^T{p}_k,$$
   与之前的等式结合, 就推出$$-1+c_2\frac{\nabla f_k^T{p}_k}{\Vert \nabla f_k{\Vert }^2}\le \frac{\nabla f_{k+1}^T{p}_{k+1}}{\Vert \nabla f_{k+1}{\Vert }^2}\le -1-c_2\frac{\nabla f_k^T{p}_k}{\Vert \nabla f_k{\Vert }^2}.$$由归纳假设, $$-1-\frac{c_2}{1-c_2}\le \frac{\nabla f_{k+1}^T{p}_{k+1}}{\Vert \nabla f_{k+1}{\Vert }^2}\le -1+\frac{c_2}{1-c_2}.$$因此不等式对$k+1$也成立. 证毕.

注意上述结论只用到了强Wolfe条件的曲率条件, 而Armijo条件则用于证明后面的全局收敛性. $\nabla f_k^T{p}_k$的界表明了$\Vert p_k\Vert$增长速度的限制, 进而在下面的分析中能起到关键性的作用.
定理7同样可用来解释FR方法的一个缺陷: 我们说, 如果FR产生了一个"较坏"的搜索方向以及一个较小的步长, 则会产生连锁效应. 同第三章一样, 令${\theta }_k$表示$p_k$和负梯度方向$-\nabla f_k$之间的夹角, 即$$\cos {\theta }_k=\frac{-\nabla f_k^T{p}_k}{\Vert \nabla f_k\Vert \Vert p_k\Vert }.$$假设$p_k$"不好", 即${\theta }_k$将近$90^{\circ }$. 在之前不等式中同乘$\Vert \nabla f_k\Vert /\Vert p_k\Vert$, 得到$$\frac{1-2c_2}{1-c_2}\frac{\Vert \nabla f_k\Vert }{\Vert p_k\Vert }\le \cos {\theta }_k\le \frac{1}{1-c_2}\frac{\Vert \nabla f_k\Vert }{\Vert p_k\Vert },k=0,1,\dots .$$从这, 我们得出$\cos {\theta }_k\approx 0$当且仅当$\Vert \nabla f_k\Vert \ll \Vert p_k\Vert$. 由于$p_k$与梯度几近正交, 所以$x_{k+1}\approx x_k$, 从而$\nabla f_{k+1}\approx \nabla f_k$, 推出${\beta }_{k+1}^{\mathrm{F}\mathrm{R}}\approx 1$. 再由$\Vert \nabla f_{k+1}\Vert \approx \Vert \nabla f_k\Vert \ll \Vert p_k\Vert$, 所以$p_{k+1}\approx p_k$. 于是新产生的搜索方向相较于原来并没有多大改善. 因此若${\theta }_k\approx 90^{\circ }$在某个$k$成立且接下来的步长很小, 则后续会产生一系列几近无效的迭代项. 此时, 我们需要在FR方法中嵌入一些重启的策略以保证收敛.

相比之下, PR方法在这样的情形下会有十分不同的表现. 由PR方法中的公式, 我们反而由${\beta }_{k+1}^{\mathrm{P}\mathrm{R}}\approx 0$. 从而下一个搜索方向$p_{k+1}$将与负梯度相近, $\cos {\theta }_{k+1}\approx 1$. 也就是说, PR方法在遇到了一个"坏"方向后会选择重启. PR+和HS方法也有同样的效果. 对于FR-PR方法, 我们已经有${\beta }_{k+1}^{\mathrm{F}\mathrm{R}}\approx 1,{\beta }_{k+1}^{\mathrm{P}\mathrm{R}}\approx 0$, 从而${\beta }_{k+1}={\beta }_{k+1}^{\mathrm{P}\mathrm{R}}$.

2.4.2 全局收敛性

不想线性共轭梯度法, 非线性共轭梯度法往往具有较为奇特的收敛性质. 下面我们对于FR和PR方法做一些说明. 为方便起见, 我们不过分地假设目标函数具有以下性质:

假设

1. 水平集 L : = { x ∣ f ( x ) ≤ f ( x 0 ) } \mathcal{L}:=\{x|f(x)\le f(x_0)\} L:={x∣f(x)≤f(x0​)}有界;
2. 在$\mathcal{L}$的某个开邻域$\mathcal{N}$中, $f$Lipschitz连续可微, 即存在常数$\bar{\gamma }$使得$$\Vert \nabla f(x)\Vert \le \bar{\gamma },\forall x\in \mathcal{L}.$$

我们主要的工具是第三章的Zoutendijk定理. 它表示, 在以上假设下, 有$$\sum _{k=0}^{\infty }{\cos }^2{\theta }_k\Vert \nabla f_k{\Vert }^2<\infty ,$$其中$p_k$为下降方向, ${\alpha }_k$为满足Wolfe条件的搜索步长. 我们可以利用这一结论证明周期重启的算法的弱全局收敛性: 若$k_1,k_2,\dots$表示重启发生的位置, 则从Zoutendijk条件得$$\sum _{k=k_1,k_2,\dots }\Vert \nabla f_k{\Vert }^2<\infty .$$从而$$\underset{j\to \infty }{\lim }\Vert \nabla f_{k_j}\Vert =0\Rightarrow \underset{k\to \infty }{\mathrm{liminf}}\Vert \nabla f_k\Vert =0.$$比起这个, 我们更在意不重启的算法的收敛性. 下面证明FR方法的弱全局收敛性.

定理8 在假设条件满足, 且算法4中线搜索满足强Wolfe条件 (其中$0<c_1<c_2<\frac{1}{2}$)时, 有$$\underset{k\to \infty }{\mathrm{liminf}}\Vert \nabla f_k\Vert =0.$$
证明: 用反证法证明. 假设不然, 存在常数$\gamma >0$使得对于充分大的$k$, $$\Vert \nabla f_k\Vert \ge \gamma .$$将$$\cos {\theta }_k\ge \frac{1-2c_2}{1-c_2}\frac{\Vert \nabla f_k\Vert }{\Vert p_k\Vert }$$代入Zoutendijk条件, 我们得到$$\sum _{k=0}^{\infty }\frac{\Vert \nabla f_k{\Vert }^4}{\Vert p_k{\Vert }^2}<\infty .$$由曲率条件和定理7, 我们有$$|\nabla f_k^T{p}_{k-1}|\le -c_2\nabla f_{k-1}^T{p}_{k-1}\le \frac{c_2}{1-c_2}\Vert \nabla f_{k-1}{\Vert }^2.$$由搜索方向公式, 有$$\begin{array}{rl}\Vert p_k{\Vert }^2& \le \Vert \nabla f_k{\Vert }^2+2{\beta }_k^{\mathrm{F}\mathrm{R}}|\nabla f_k^T{p}_{k-1}|+({\beta }_k^{\mathrm{F}\mathrm{R}}{)}^2\Vert p_{k-1}{\Vert }^2\\ & \le \Vert \nabla f_k{\Vert }^2+\frac{2c_2}{1-c_2}{\beta }_k^{\mathrm{F}\mathrm{R}}\Vert \nabla f_{k-1}{\Vert }^2+({\beta }_k^{\mathrm{F}\mathrm{R}}{)}^2\Vert p_{k-1}{\Vert }^2\\ & =\left(\frac{1+c_2}{1-c_2}\right)\Vert \nabla f_k{\Vert }^2+({\beta }_k^{\mathrm{F}\mathrm{R}}{)}^2\Vert p_{k-1}{\Vert }^2.\end{array}$$定义$c_3\stackrel{def}{=}(1+c_2)/(1-c_2)\ge 1$, 重复应用上述, 可得$$\begin{array}{rl}\Vert p_k{\Vert }^2& \le c_3\Vert \nabla f_k{\Vert }^2+({\beta }_k^{\mathrm{F}\mathrm{R}}{)}^2(c_3\Vert \nabla f_{k-1}{\Vert }^2+({\beta }_{k-1}^{\mathrm{F}\mathrm{R}}{)}^2(c_3\Vert f_{k-2}{\Vert }^2+\cdots +({\beta }_1^{\mathrm{F}\mathrm{R}}{)}^2\Vert p_0{\Vert }^2))\cdots )\\ & =c_3\Vert \nabla f_k{\Vert }^4\sum _{j=0}^k\Vert \nabla f_j{\Vert }^{-2}.(\because ({\beta }_k^{\mathrm{F}\mathrm{R}}{)}^2({\beta }_{k-1}^{\mathrm{F}\mathrm{R}}{)}^2\cdots ({\beta }_{k-i}^{\mathrm{F}\mathrm{R}}{)}^2=\frac{\Vert \nabla f_k{\Vert }^4}{\Vert \nabla f_{k-i-1}{\Vert }^4},p_0=-\nabla f_0)\end{array}$$由Lipshcitz连续可微条件以及一开始的反证条件, 可得$$\Vert p_k{\Vert }^2\le \frac{c_3{\bar{\gamma }}^4}{{\gamma }^2}(k+1),$$于是$$\sum _{k=1}^{\infty }\frac{1}{\Vert p_k{\Vert }^2}\ge {\gamma }_4\sum _{k=1}^{\infty }\frac{1}{k+1}=\infty ,$$其中${\gamma }_4$为一正常数.
另一方面, 从$$\Vert \nabla f_k\Vert \ge \gamma$$和$$\sum _{k=0}^{\infty }\frac{\Vert \nabla f_k{\Vert }^4}{\Vert p_k{\Vert }^2}<\infty$$可得$$\sum _{k=1}^{\infty }\frac{1}{\Vert p_k{\Vert }^2}<\infty .$$因而矛盾! 所以FR方法的弱全局收敛性成立.

从定理的证明过程中可知, 此全局收敛性对于所有$|{\beta }_k|\le {\beta }_k^{\mathrm{F}\mathrm{R}}$均成立, 例如FR-PR方法.
特别地, 我们可以增强条件以获得更强的全局收敛性质. 例如若存在常数$c_4,c_5>0$使得$$\cos {\theta }_k\ge c_4\frac{\Vert \nabla f_k\Vert }{\Vert p_k\Vert },\frac{\Vert \nabla f_k\Vert }{\Vert p_k\Vert }\ge c_5>0,k=1,2,\dots ,$$则从Zoutendijk条件即知$$\underset{k\to \infty }{\lim }\Vert \nabla f_k\Vert =0.$$事实上, 在$f$强凸时使用带精确线搜索的PR方法即可得到全局收敛性. 而对一般的(非凸)函数, 即使使用了"理想"的步长${\alpha }_k$ (这里"理想"指的是${\alpha }_k$是一阶稳定点), PR方法仍然可能会陷入循环无法到达解.

定理9 考虑带理想线搜索的PR方法. 存在二阶连续可微函数$f:{\mathbb{R}}^3\to \mathbb{R}$以及初始点$x_0\in {\mathbb{R}}^3$使得梯度范数序列 { ∥ ∇ f k ∥ } \{\Vert\nabla f_k\Vert\} {∥∇fk​∥}始终位于0的某个邻域以外.

定理9的证明中需要两相邻搜索方向几近相反, 即$p_{k+1}\approx -\alpha p_k,\alpha >0$. 而在理想线搜索的前提下, 有$\nabla f_{k+1}^T{p}_k=0$. 由搜索方向公式, $p_{k+1}=-\nabla f_{k+1}+{\beta }_{k+1}{p}_k\Rightarrow p_k=\frac{1}{\alpha +{\beta }_{k+1}}\nabla f_{k+1}$. 代入前一个等式得$\Vert \nabla f_{k+1}{\Vert }^2\cdot \frac{1}{\alpha +{\beta }_{k+1}}=0$. 从而必有${\beta }_{k+1}<0$. 这里就需要PR+方法重启. 我们之前谈到, 在PR+方法的线搜索上只需对Wolfe条件加微小修正即可保证PR+方法的所有搜索方向均为下降方向. 因此, 我们可以对PR+方法证明类似定理8的全局收敛性. 而对于那两个比较奇怪的${\beta }_k$公式, 令人惊讶的是, 它们的全局收敛性甚至不需要对Wolfe条件作任何修正即可获取.