共轭梯度法是求解大规模无约束优化问题$ \min\{f(x)|x\in {{\Bbb R}}^n\} $的有效方法之一, 其迭代公式为

$ x_{k+1} = x_k+\alpha_k d_k, $

$ d_k = \left\{\begin{array}{ll} -g_k, & k = 0, \\ -g_k+\beta _k d_{k-1}, & k\geq1, \end{array}\right. $

其中, $ g_k = \nabla f(x_k) $, $ d_k $为搜索方向, $ \alpha_k $为步长, $ \beta_k $为搜索方向$ d_k $的调控参数.不同的$ \beta_k $对应不同的共轭梯度法.比较经典的共轭梯度法有Hestenes-Stiefel(HS)方法^[1], Fletcher-Reeves(FR)方法^[2], Polak-Ribiére-Polyak(PRP)方法^{[3, 4]}和Dai-Yuan(DY)方法^[5], 它们的调控参数$ \beta_k $为

$ \beta_k^{HS} = \frac{g_k^T(g_k-g_{k-1})}{d_{k-1}^T(g_k-g_{k-1})}, {\quad} \beta_k^{FR} = \frac{\|g_{k}\|^2}{\|g_{k-1}\|^2}, $

$ \beta_k^{PRP} = \frac{g_k^T(g_k-g_{k-1})}{\|g_{k-1}\|^2}, {\quad} \beta_k^{DY} = \frac{\|g_{k}\|^2}{d_{k-1}^T(g_k-g_{k-1})}. $

当采用精确线性搜索并且目标函数为二次凸函数时, 以上四种共轭梯度法是等价的.然而不采用精确线搜索或者目标函数为一般函数时, 这四种方法的理论结果和数值效果各不相同, 其中分子为$ \|g_{k}\|^2 $的FR方法以及DY方法数值结果较差, 但具有良好的理论结果; 分子为$ g_k^T(g_k-g_{k-1}) $的PRP方法以及HS方法的数值效果比较好, 但是理论结果比较差.为得到数值实验和理论结果都很好的共轭梯度方法, 许多学者对这些经典方法做了不少修正.

PRP方法之所以数值实验效果比较好, 是因为对一般函数来说, PRP方法具有自动重启的优良性能.但是针对它的理论结果较差这一点, Powell在文献[6]中建议将$ \beta_k $限制为非负, 基于此建议Gilbert和Nocedal^[7]考虑如下的$ PRP^+ $方法

$ \begin{eqnarray*} \beta^{PRP^+} = \max\{\beta^{PRP}, 0\}, \end{eqnarray*} $

得出$ PRP^+ $方法在假设充分下降条件满足时, 对一般函数采用Wolfe线搜索或者精确线搜索具有全局收敛性.

为了保证PRP公式非负, Wei、Yao和Liu在文献[8]提出了PRP方法的一个变型, 称为WYL方法或者VPRP方法, 其中调整参数$ \beta_k $为

$ \begin{eqnarray*} \beta_k^{WYL} = \frac{g_k^T(g_k-\frac{\|g_k\|}{\|g_{k-1}\|}g_{k-1})}{\|g_{k-1}\|^2}. \end{eqnarray*} $

WYL方法不仅具有PRP方法的良好数值性能而且在强Wolfe线搜索条件下是全局收敛的.之后, Huang、Wei以及Yao在文献[9]中证明了当$ \sigma<\frac{1}{4} $时($ \sigma $为Wolfe线搜索的参数), 该方法在强Wolfe线搜索下具有充分下降性以及全局收敛性.同时, 也有不少学者在WYL方法基础上对调整参数$ \beta_k $加以修改, 如文献[10-14].

Zhang^[15]提出了一个基于PRP的三项共轭梯度法, 称为MPRP方法, 其公式如下

$ d_k = \left\{\begin{array}{ll} -g_k, & k = 0, \\ -g_k+\beta_k^{PRP} d_{k-1}-\theta_k y_{k-1}, & k\geq1, \end{array}\right. $

$ \theta_k = \frac{g_k^Td_{k-1}}{\|g_{k-1}\|^2}, $

其中, 步长$ \alpha_k $是在修改的Armijo线搜索下, 令$ \alpha_k = \max\{\rho^j, j = 0, 1, 2, \cdots \} $, 满足

(1.1) $ \begin{equation} f(x_k+\alpha_k d_k)\leq f(x_k)-\delta \alpha_k^2\|d_k\|^2, \; \delta \in (0, 1), \end{equation} $

是充分下降且全局收敛的.有不少学者在两项共轭梯度法的基础上加以修改, 提出三项共轭梯度法, 如文献[16-20].修改之后的三项共轭梯度法在数值效果和理论结果上比两项有很大改善.

文献[21]提出了一类新的WYL型三项共轭梯度法, 称为MWYL方法, 其公式如下

$ d_k^{MWYL} = -g_k+\beta_k^{MWYL} d_{k-1}+\theta_k^{MWYL} g_{k-1}, $

$ \beta_k^{MWYL} = \frac{g_k^T(g_k-\frac{\|g_k\|}{\|g_{k-1}\|}g_{k-1})}{\|g_{k-1}\|^2+\max\{\alpha_{k-1}, c\}\|d_{k-1}\|^2}, $

$ \theta_k^{MWYL} = \frac{g_k^T d_{k-1}}{\|g_{k-1}\|^2+\max\{\alpha_{k-1}, c\}\|d_{k-1}\|^2}\cdot \frac{\|g_k\|}{\|g_{k-1}\|}, $

参数$ c\geq0 $为常数, 且在强Wolfe线搜索方式下具有充分下降及全局收敛性.

受此启发, 下面给出一个新的WYL型三项共轭梯度算法并利用Zhang^[15]的线搜索方式(1.1), 算法公式为

(1.2) $ \begin{equation} d_k = \left\{\begin{array}{ll} -g_k, & k = 0, \\ { } -g_k+\beta_k^{TWYL} d_{k-1}-\theta_k (g_k-\frac{\|g_k\|}{\|g_{k-1}\|}g_{k-1}), & k\geq1, \end{array}\right. \end{equation} $

(1.3) $ \begin{equation} \beta_k^{TWYL} = \frac{g_k^T(g_k-\frac{\|g_k\|}{\|g_{k-1}\|}g_{k-1})}{(1+v\frac{|g_k^Tg_{k-1}|}{\|g_k\|\|g_{k-1}\|})\|g_{k-1}\|^2}, \end{equation} $

(1.4) $ \begin{equation} \theta_k = \frac{g_k^Td_{k-1}}{\|g_{k-1}\|^2}, \end{equation} $

其中步长$ \alpha_k $满足(1.1)式.当$ g_k^Td_{k-1}>0 $时, 令参数$ v = 0 $; 否则$ v>0 $.

本文组织如下:第一部分引言介绍了本文提出算法的思路; 接下来第二部分将介绍算法的流程以及证明算法是充分下降的; 第三部分证明算法在修改的Armijo线搜索条件下是全局收敛的; 第四部分对算法进行数值实验, 并与四个相关的经典算法进行比较.

本文所提的新方法描述如下:

步骤0   选定初始步$ x_0 $, 参数$ \rho, \delta\in(0, 1) $. $ v\geq0 $, $ \varepsilon>0 $.并令$ k = 0 $;

步骤1  计算$ g_k = \nabla f(x_k) $.若$ \|g_k\|\leq \varepsilon $, 则停止计算; 否则由(1.2)–(1.4)式计算$ d_k $;

步骤2  按照(1.1)式计算步长$ \alpha_k $;

步骤3  计算$ x_{k+1} = x_k+\alpha_k d_k $, 令$ k: = k+1 $, 转步骤1.

为表示方便, 记为算法1.此外, 关于方向$ d_k $, 可以得到如下引理.

引理2.1  若$ d_k $由公式(1.2)–(1.4)定义, 则$ d_k $满足下降条件$ g_k^Td_k\leq -\|g_k\|^2 $.

证  当$ k = 0 $时, 容易得$ g_0^Td_0 = -\|g_0\|^2 $; 当$ k>0 $时, 在(1.2)式两边同时与$ g_k $做内积得

$ g_k^Td_k = -\|g_k\|^2+ \frac{g_k^T(g_k-\frac{\|g_k\|}{\|g_{k-1}\|}g_{k-1})g_k^Td_{k-1}}{(1+v\frac{|g_k^Tg_{k-1}|}{\|g_k\|\|g_{k-1}\|})\|g_{k-1}\|^2} - \frac{g_k^Td_{k-1}g_k^T(g_k-\frac{\|g_k\|}{\|g_{k-1}\|}g_{k-1})}{\|g_{k-1}\|^2}. $

由算法1的定义知, 当$ g_k^Td_{k-1}>0 $时$ v = 0 $, 故

(2.1) $ \begin{equation} g_k^Td_k = -\|g_k\|^2; \end{equation} $

否则$ v>0 $, 有

(2.2) $ \begin{equation} g_k^Td_k\leq-\|g_k\|^2. \end{equation} $

引理2.1得证, 表明算法1在不依赖于任何搜索方式下满足充分下降条件.

本节分析算法的收敛性, 首先给出相关的引理.

引理3.1  如果$ f $是下有界的, 有$ \sum\limits_{k = 0}^\infty \alpha_k^2 \|d_k\|^2<\infty $, 从而$ \lim\limits_{k\rightarrow \infty}\alpha_kd_k = 0 $.

证  由引理2.1看出算法1的搜索方向是充分下降的.此外, 从(1.1)式可以看出函数值序列$ \{f(x_k)\} $是递减的.因为$ f $是下有界的, 从式(1.1)中可以得到

$ \begin{eqnarray*} \sum\limits_{k = 0}^\infty \alpha_k^2 \|d_k\|^2<\infty, \end{eqnarray*} $

因此有

(3.1) $ \begin{equation} \lim\limits_{k\rightarrow \infty}\alpha_k\|d_k\| = 0.\\ \end{equation} $

引理3.1证毕.

为证明算法1的全局收敛性, 首先做出如下必要假设$ \mathrm{H} $:

($ \mathrm{H} $1) $ f(x) $在水平集$ \Omega = \{x\in R^n| f(x)\leq f(x_0)\} $是有界的.

($ \mathrm{H} $2) $ f(x) $在水平集$ \Omega $的一个邻域$ N $内连续可微, 且梯度函数$ g(x) = \nabla f(x) $满足Lipschitz条件, 即存在常数$ L>0 $, 使得

(3.2) $ \begin{equation} \|g(x)-g(y)\|\leq L\|x-y\|, \; \; \forall \; x, y\in N. \end{equation} $

由于$ \{f(x_k)\} $是递减的, 由算法1产生的序列$ \{x_k\}\subset \Omega $.此外, 从假设H中可知, 存在一个常数$ \gamma_1>0 $, 使得

(3.3) $ \begin{equation} \|g(x)\|\leq \gamma_1, \; \; \; \; \; \forall x \in \Omega. \end{equation} $

引理3.2  若存在$ \varepsilon >0 $, 使得对任意的$ k $有

(3.4) $ \begin{equation} \|g_k\|>\varepsilon, \end{equation} $

从而存在一个常数$ M>0 $, 使得对任意的$ k $有

(3.5) $ \begin{equation} \|d_k\|\leq M. \end{equation} $

证  根据$ d_k $的公式、Cauchy-Schwarz不等式以及公式(3.2)–(3.4)有

$ \begin{eqnarray*} \|d_k\|&\leq &\|g_k\| + \frac{\|g_k^T(g_k-\frac{\|g_k\|}{\|g_{k-1}\|}g_{k-1})\| \|d_{k-1}\|}{(1+v\frac{|g_k^Tg_{k-1}|}{\|g_k\|\|g_{k-1}\|})\|g_{k-1}\|^2} + \frac{\|g_k\| \|d_{k-1}\| \|g_k-\frac{\|g_k\|}{\|g_{k-1}\|}g_{k-1}\|}{\|g_{k-1}\|^2} \\ &\leq& \|g_k\| + \frac{4\|g_k\|\|g_k-g_{k-1}\|\|d_{k-1}\|}{\|g_{k-1}\|^2} \\ &\leq &\gamma_1+\frac{4\gamma_1 L \alpha_{k-1}\|d_{k-1}\|}{\varepsilon^2}\|d_{k-1}\|. \end{eqnarray*} $

因为$ \lim\limits_{k\rightarrow \infty}\alpha_k \|d_k\| = 0 $, 所以存在一个常数$ \gamma \in (0, 1) $, 并且对于一个正整数$ K $, 当$ k>K $时有下面不等式成立

$ 4\gamma_1 L \alpha_{k-1}\|d_{k-1}\| \leq \gamma, $

因此, 对于任意的$ k>K $, 有

$ \|d_k\|\leq \gamma_1 + \gamma \|d_{k-1}\| \leq \gamma_1(1+\gamma+\gamma^2+\cdots + \gamma^{k-K-1})+\gamma^{k-K} \|d_K\| \leq \frac{\gamma_1}{1-\gamma} + \|d_K\|. $

令

$ M = \max \{\|d_1\|, \|d_2\|, \cdots , \|d_K\|, \frac{\gamma_1}{1-\gamma} + \|d_K\|\}, $

故(3.5)式成立.证毕.

现在我们给出算法1的全局收敛性.

定理3.1  若假设$ \mathrm{H} $成立, 则有

(3.6) $ \begin{equation} \lim\limits_{k\rightarrow \infty} \inf \|g_k\| = 0. \end{equation} $

证  反证法.假设结论不成立, 因此存在一个常数$ \varepsilon >0 $, 使得

(3.7) $ \begin{equation} \|g_k\|>\varepsilon, \; \; \; \forall k. \end{equation} $

如果$ \lim\limits_{k\rightarrow \infty} \inf \alpha_k>0 $, 由(2.1)、(2.2)以及(3.1)式中得到$ \lim\limits_{k\rightarrow \infty} \inf \|g_k\| = 0 $, 与(3.7)式矛盾.如果$ \lim\limits_{k\rightarrow \infty} \inf \alpha_k = 0 $, 有一个无限集合$ N $, 使得

$ \begin{eqnarray*} \lim\limits_{k\rightarrow \infty}\inf\limits_{k\in N}\alpha_k = 0. \end{eqnarray*} $

根据算法1中的步骤2, 当$ k\in N $充分大时, $ \rho^{-1}\alpha_k $并不满足(1.1)式, 即

(3.8) $ \begin{equation} f(x_k+\rho^{-1}\alpha_kd_k)-f(x_k)>-\delta \rho^{-2}\alpha_k^2\|d_k\|^2. \end{equation} $

再利用中值定理, 以及引理3.2, (2.1)、(2.2)以及(3.2)式, 存在$ h_k\in(0, 1) $, 使得

$ \begin{eqnarray*} f(x_k+\rho^{-1}\alpha_kd_k)-f(x_k)& = &\rho^{-1}\alpha_kg(x_k+h_k\rho^{-1}\alpha_kd_k)^Td_k \\ & = &\rho^{-1}\alpha_kg_k^Td_k+\rho^{-1}\alpha_k(g(x_k+h_k\rho^{-1}\alpha_kd_k)-g_k)^Td_k \\ &\leq &\rho^{-1}\alpha_kg_k^Td_k+L \rho^{-2}\alpha_k^2\|d_k\|^2. \end{eqnarray*} $

将上面的不等式代入(3.8)式, 故当$ \; \forall k\in N $充分大时, 有

$ \|g_k\|^2\leq \rho^{-1}(L+\delta)\alpha_k\|d_k\|^2. $

因为$ \{d_k\} $是有界的, 并且$ \lim\limits_{k\rightarrow \infty}\inf\alpha_k = 0 $, 得到$ \lim\limits_{k\in N, k\rightarrow \infty}\|g_k\| = 0 $.矛盾.故定理成立.

本节展示算法的数值结果.在文献[22]中选取了29个函数, 分别比较了该文所参考的MPRP方法、WYL方法、MWYL方法、同类经典的HZ方法^[23]和提出的TWYL方法的数值性能.其中MPRP、TWYL方法采用修改的Armijo线搜索(1.3), WYL、MWYL采用强Wolfe线搜索, HZ采用修改的Wolfe线搜索求步长.测试环境为Win10操作系统, Intel(R) Core(TM) i7-9750H CPU 2.60GHZ 16.0GB内存.强Wolfe线搜索相关参数: $ \rho = 0.45 $, $ \delta = 0.0002 $, $ \sigma = 0.4 $. Armijo线搜索中参数为: $ \rho = 0.5 $, $ \delta = 0.2 $.算法的终止条件准则为以下二者情形之一:

(1) $ \|g_k\|\leq\varepsilon(1+f(x_k)) $, ($ \varepsilon = 10^{-6} $);

(2) $ k>10000. $

表 1给出了29个测试函数的编号及名称.

接下来采用Dolan和More^[24]提出的方法, 画出数据的性能图 1–4.从四个图中可以直观看到MPRP方法、WYL方法、MWYL方法、TWYL方法、HZ方法分别在算法迭代次数、函数迭代次数、梯度迭代次数以及CPU运行时间四个方面的数值结果的比较.明显可以看出本文提出的TWYL方法在这四个方面都优于MPRP方法、WYL方法、MWYL方法、HZ方法.其中, 在迭代次数以及时间性能方面MPRP方法与MWYL方法比较接近.

表 2表示的是MPRP方法、WYL方法、MWYL方法、TWYL方法、HZ方法对算法迭代次数$ (k) $、目标函数迭代次数$ (q) $、梯度迭代次数$ (r) $以及算法的CPU运行时间$ (t) $四个重要指标的数值结果.其中表头中N表示函数的编号, D表示函数的维度, 加黑的数据表示五种方法中效果最好.从表中可以看出所提的新方法大部分情况下都能得到最好的结果:较少的迭代次数和CPU时间.这和性能图表现的一致.因此, 数值结果表明新算法是有效的.