基于拟牛顿方程一个改进的非线性共轭梯度算法

A Modified Nonlinear Conjugate Gradient Algorithm Using Secant Conditions

摘要: 基于拟牛顿方程，借鉴文[1]的思想，本文给出一个改进的具有全局收敛非线性共轭梯度算法。该文的算法对文[1]中算法进行了改进，所产生的搜索方向具有充分下降性。在温和的假设下，新算法具有全局收敛性。最后，数值结果检验其有效性。

Abstract: In this paper, based on the idea of Ref. [1], we propose a modified conjugate gradient method using secant conditions for unconstrained optimization problems. The proposed algorithm improves the method in Ref. [1], which possesses the following properties: the search direction has the sufficient descent property; the global convergence of the given algorithm will be established under suitable assumptions; numerical results are reported to test its efficiency.

1. 引言

本文考虑求解如下无约束优化问题

$\underset{x\in R^n}{\min }f\left(x\right),$ (1.1)

其中 $f\left(x\right):R^n\to R$ 连续可微，此问题在经济管理、生产运作、工程设计等领域具有广泛的应用背景。共轭梯度法是求解问题(1.1)的有效方法之一，该方法具有算法结构简单、计算存储量少、数值效果明显的优点，其迭代公式的一般形式为：

$x_{k+1}=x_k+{\alpha }_k{d}_k,k=0,1,2,\cdots$

其中 $x_k$ 为第 $k$ 次迭代点， ${\alpha }_k>0$ 为搜索步长， $d_k$ 为搜索方向且其更新方式为

$d_k=\{\begin{array}{l}-g_k+{\beta }_k{d}_{k-1},k\ge 1,\\ -g_k,k=0,\end{array}$ (1.2)

其中 ${\beta }_k\in \Re$ 是方向调控参数。数值优化方法是根据搜索方向的选择来定义其方法类型，从上式可以看出，搜索方向 $d_k$ 是由参数 ${\beta }_k$ 决定的，所以根据 ${\beta }_k$ 的选取方式可以定义不同的共轭梯度公式，下面给出几个著名的公式：

${\beta }_k^{PRP}=\frac{g_{k+1}^T\left(g_{k+1}-g_k\right)}{{\Vert g_k\Vert }^2}$ [2] [3] ， ${\beta }_k^{HS}=\frac{g_{k+1}^T\left(g_{k+1}-g_k\right)}{{\left(g_{k+1}-g_k\right)}^T{d}_k}$ [4] ， ${\beta }_k^{FR}=\frac{g_{k+1}^T{g}_{k+1}}{{\Vert g_k\Vert }^2}$ [5] ，

${\beta }_k^{DY}=\frac{g_{k+1}^T{g}_{k+1}}{{\left(g_{k+1}-g_k\right)}^T{d}_k}$ [6] ， ${\beta }_k^{CD}=\frac{g_{k+1}^T{g}_{k+1}}{-d_k^T{g}_k}$ [7] ， ${\beta }_k^{LS}=\frac{g_{k+1}^T\left(g_{k+1}-g_k\right)}{-g_k{}^T{d}_k}$ [8] ，

其中 $g_k=\nabla f\left(x_k\right)$ 和 $g_{k+1}=\nabla f\left(x_{k+1}\right)$ 表示函数 $f\left(x\right)$ 在 $x_k$ 和 $x_{k+1}$ 的梯度值， $\Vert .\Vert$ 是欧氏向量范数。众多优化专家都希望能找到理论性质良好且数值表现又好的算法，已取得众多成果(见 [9] - [25] )，其中Yuan [9] 给出如下PRP修改公式：

${\beta }_k^{MPRP}={\beta }_k^{PRP}-\min \left\{{\beta }_k^{PRP},\frac{\mu {\Vert y_k\Vert }^2}{{\Vert g_k\Vert }^4}{g}_{k+1}^T{d}_k\right\},$

其中 $\mu >\frac{1}{4}$ 是常数， $y_k=g_{k+1}-g_k$ ，同时又将此公式推广到了其它5个公式中，得到了修改的FR公式、修改的HS公式、修改的DY公式、修改的CD公式和修改的LS公式，其中文 [1] 修改的LS为：

${\beta }_k^{MMLS+}={\beta }_k^{MMLS}-\min \left\{{\beta }_k^{MMLS},\frac{\mu {\Vert y_{{}_k}^m\Vert }^2}{{\left(-d_k^T{g}_k\right)}^2}{g}_{k+1}^T{d}_k\right\}$ ， (1.3)

其中 $y_k^m=g_{k+1}-\frac{\Vert g_{k+1}\Vert }{\Vert g_k\Vert }{g}_k$ ， ${\beta }_k^{MMLS}=\frac{g_{k+1}^T{y}_k^m}{-d_k^T{g}_k}$ 。

本文将在此公式的基础上进一步研究，从上式可以看出公式中不含函数值信息，这势必影响此方法的有效性。众所周知，拟牛顿方法不仅含有梯度值信息也含有函数值信息，且具有更好的理论和实际表现 [10] [11] 。也有许多学者将拟牛顿法的技术思想应用到共轭梯度公式中，从而获得更好的理论性质和实际数值效果 [12] [13] 。Zhang等 [11] 给出了一个如下形式的拟牛顿方程

$B_{k+1}{S}_k=y_k^m$ ， (1.4)

其中 $y_k^m=y_k+{\gamma }_k^1{s}_k$ ， $s_k=x_{k+1}-x_k$ 及

${\gamma }_k^1=\frac{3{\left(g_{k+1}+g_k\right)}^T{s}_k+6\left(f\left(x_k\right)-f\left(x_{k+1}\right)\right)}{{\Vert s_k\Vert }^2},$

此公式包含有梯度值信息和函数值信息，理论上能更高阶地近似目标函数的Hesse矩阵。基于公式(1.3)，(1.4)和文 [12] [13] 的思想，本文给出如下修改的LS共轭梯度公式

${\beta }_k^{MMLS*}={\beta }_k^{MMLS}-\min \left\{{\beta }_k^{MMLS},\frac{\mu {\Vert y_k^m\Vert }^2}{{\left(-d_k^T{g}_k\right)}^2}{g}_{k+1}^T{d}_k\right\}$ ，(1.5)

其中 ${\beta }_k^{MMLS}=\frac{g_{k+1}^T{y}_k^m}{-d_k^T{g}_k}$ 。相对于原LS方法，本文的创新点主要有：1) 新公式能保证方向调控参数 ${\beta }_k^{MMLS*}\ge 0$ ；2) 新方法具有充分下降性；3) 新方法不仅使用了梯度值信息还有函数值信息。

2. 算法的描述

基于搜索方向，给出一个基于拟牛顿方程改进的非线性共轭梯度算法。

算法1：

步骤0：给定 $x_0\in {\Re }^n,\delta \in \left(0,1/2\right),\sigma \in \left(\delta ,1\right)$ ， $\epsilon >0$ ，令 $d_0=-g_0=-\nabla f\left(x_0\right)$ ，置 $k:=0$ 。

步骤1：若 $\Vert g_k\Vert \le \epsilon$ ，停止。否则，进入步骤2。

步骤2：利用弱Wolfe-Powell(WWP)线搜索技术产生步长 ${\alpha }_k$

$f\left(x_k+{\alpha }_k{d}_k\right)\le f_k+\delta {\alpha }_k{g}_k^T{d}_k,$ (2.1)

$g{\left(x_k+{\alpha }_k{d}_k\right)}^T{d}_k\ge \sigma g_k^T{d}_k.$ (2.2)

步骤3：令 $x_{k+1}=x_k+{\alpha }_k{d}_k$ 。如果 $\Vert g_{k+1}\Vert \le \epsilon$ ，停止。

步骤4：利用公式(1.5)计算 ${\beta }_k^{MMLS*}$ ，并按如下公式计算搜索方向

$d_{k+1}=-g_{k+1}+{\beta }_k^{MMLS*}{d}_k,$ (2.3)

步骤5：置 $k:=k+1$ ，转步骤2。

3. 算法的充分下降性和全局收敛性

引理3.1：对所有 $k\ge 0$ ，搜索方向 $d_k$ 满足下面两式

$d_k^T{g}_k\le -c{\Vert g_k\Vert }^2$ (3.1)

和

$d_k^T{y}_k\ge c\left(1-\sigma \right){\Vert g_k\Vert }^2$ (3.2)

$c>0$ 是常数。

证明：首先证明(3.1)成立。根据算法1，如果 $k=0$ ，则 $g_0^T{d}_0=-{\Vert g_0\Vert }^2$ ，式(3.1)显然成立。下面利用归纳法，当 $k\ge 1$ ，假如公式(2.3)满足(3.1)，对 $k+1$ ，有

$\begin{array}{c}{g}_{k+1}^T{d}_{k+1}=-{\Vert g_{k+1}\Vert }^2+{\beta }_k^{MMLS*}{d}_k^T{g}_{k+1}\\ =-{\Vert g_{k+1}\Vert }^2+\left[\frac{g_{k+1}^T{y}_k^m}{-d_k^T{g}_k}-\min \left\{\frac{g_{k+1}^T{y}_k^m}{-d_k^T{g}_k},\frac{\mu {\Vert y_k^m\Vert }^2}{{\left(-d_k^T{g}_k\right)}^2}{g}_{k+1}^T{d}_k\right\}\right]d_k^T{g}_{k+1}\end{array}$ (3.3)

由 $-d_k^T{g}_k>0.$ 取 $u=\frac{\sqrt{-d_k^T{g}_k}}{\sqrt{2\mu }}{g}_{k+1},v=\frac{\sqrt{2\mu }{g}_{k+1}^T{d}_k}{\sqrt{-d_k^T{g}_k}}{y}_k^m$ 。下面我们分两种情形分别讨论(3.3)式：

情形I：如果 $\frac{g_{k+1}^T{y}_k^m}{-d_k^T{g}_k}<\frac{\mu {\Vert y_k^m\Vert }^2}{{\left(-d_k^T{g}_k\right)}^2}{g}_{k+1}^T{d}_k$ ，显然 $g_{k+1}^T{d}_{k+1}=-{\Vert g_{k+1}\Vert }^2$ 成立。

情形II：如果 $\frac{g_{k+1}^T{y}_k^m}{-d_k^T{g}_k}\ge \frac{\mu {\Vert y_k^m\Vert }^2}{{\left(-d_k^T{g}_k\right)}^2}{g}_{k+1}^T{d}_k$ ，(3.3)式可化为

$\begin{array}{c}{g}_{k+1}^T{d}_{k+1}=-{\Vert g_{k+1}\Vert }^2+\left(\frac{g_{k+1}^T{y}_k^m}{-d_k^T{g}_k}-\frac{\mu {\Vert y_k^m\Vert }^2}{{\left(-d_k^T{g}_k\right)}^2}{g}_{k+1}^T{d}_k\right)d_k^T{g}_{k+1}\\ =\frac{d_k^T{g}_{k+1}{g}_{k+1}^T{y}_k^m-{\Vert g_{k+1}\Vert }^2\left(-d_k^T{g}_k\right)-\frac{\mu {\Vert y_k^m\Vert }^2}{-d_k^T{g}_k}{\left(g_{k+1}^T{d}_k\right)}^2}{-d_k^T{g}_k}\\ =\frac{u^{T}v-\frac{1}{2}\left({\Vert u\Vert }^2+{\Vert v\Vert }^2\right)}{-d_k^T{g}_k}+\frac{-\left(1-\frac{1}{4\mu }\right){\Vert g_{k+1}\Vert }^2{d}_k^T{y}_k^m}{-d_k^T{g}_k}\\ \le -\left(1-\frac{1}{4\mu }\right){\Vert g_{k+1}\Vert }^2,\end{array}$

上述不等式利用了 $u^{T}v\le \frac{1}{2}\left({\Vert u\Vert }^2+{\Vert v\Vert }^2\right)$ 的关系。取 $c=1-\frac{1}{4\mu }$ ，则(3.1)成立。根据线搜索技术(2.2)可推出(3.2)也成立。证毕。

我们将采用反证法来证明算法1的全局收敛性，假设存在常数 ${\epsilon }_0>0$ 满足

$\Vert g_k\Vert \ge {\epsilon }_0,\forall k\ge 0.$ (3.4)

根据(3.4)导出矛盾，从而证明 $\Vert g_k\Vert \to 0,k\to \infty$ 成立。

为证明算法1的收敛性，还需要下面的常规假设条件。

假设A：1) 定义水平集 $\Omega =\left\{x\in {\Re }^n:f\left(x\right)\le f\left(x_0\right)\right\}$ 且有界，存在 $L_0>0$ 满足

$\Vert {\nabla }^{2}f\left(x\right)\Vert \le L_0,x\in \Omega .$

2) 目标函数 $f$ 在水平集 $\Omega$ 上连续可微并有下界，目标函数梯度满足Lipschitz条件，所谓Lipschitz条件就是存在常数 $L>0$ 满足下式

$\Vert g\left(x\right)-g(y)\Vert \le L\Vert x-y\Vert ,\forall x,y\in \Omega$ 。 (3.5)

下面引理3.2和引理3.3在共轭梯度算法文献中有类似的结论，本文只给出引理3.3证明过程。

引理3.2：设假设A成立，序列 $\left\{g_k\right\}$ 和 $\left\{d_k\right\}$ 由算法1产生，则 $d_k\ne 0$ 及 ${\displaystyle \underset{k=0}{\overset{\infty }{\sum }}{\Vert u_{k+1}-u_k\Vert }^2}<\infty$ ，其中 $u_k=\frac{d_k}{\Vert d_k\Vert }$ 。

Gilbert和Nocedal [14] 给出下面性质(P)，具体内容是：

性质(P)对于两参数 ${\gamma }_1,{\gamma }_2>0$ ，且满足

$\text{0}<{\gamma }_1\le \Vert g_k\Vert \le {\gamma }_2$ 。 (3.6)

若对所有 $k$ ，存在常数 $b>1$ 和 $\lambda >0$ 满足 $\left|{\beta }_k\right|\le b$ 和 $\Vert s_k\Vert \le \lambda$ ，有 $\left|{\beta }_k\right|\le \frac{1}{2b}$ 成立。

引理3.3：如果假设A满足且序列 $\left\{g_k\right\}$ 和 $\left\{d_k\right\}$ 由算法1产生，且存在常数 $M>0$ 使得 $\Vert d_k\Vert \le M$ 成立，则 ${\beta }_k^{MMLS*}$ 满足性质(P)。

证明：对于公式(2.3)，如果 $\frac{g_{k+1}^T{y}_k^m}{-d_k^T{g}_k}\le \frac{\mu {\Vert y_k^m\Vert }^2}{{\left(-d_k^T{g}_k\right)}^2}{g}_{k+1}^T{d}_k$ 成立，结论显然成立。否则，利用假设A(1)，可推出存在常数 $M_1>0$ 使下式

$\Vert s_k\Vert \le M_1$ (3.7)

成立。利用 $\Vert d_k\Vert \le M$ ，(3.1)，(3.2)，(3.5)~(3.7)，有

$\begin{array}{c}\left|{\beta }_k^{MMLS*}\right|\le \left|\frac{g_{k+1}^T{y}_k^m}{d_k^T{g}_k}\right|+\frac{\mu {\Vert y_k^m\Vert }^2}{{\left(-d_k^T{g}_k\right)}^2}\left|g_{k+1}^T{d}_k\right|\le \frac{\Vert g_{k+1}\Vert \Vert \Vert g_{k+1}-g_k\Vert +3\left[\Vert g_{k+1}-g_k\Vert +\Vert {\nabla }^{2}f\left({\varsigma }_k\right)s_k\Vert \right]\Vert }{c{\Vert g_k\Vert }^2}\\ +\frac{\mu \Vert g_{k+1}\Vert \Vert d_k\Vert {\Vert \Vert g_{k+1}-g_k\Vert +3\left[\Vert g_{k+1}-g_k\Vert +\Vert {\nabla }^{2}f\left({\varsigma }_k\right)s_k\Vert \right]\Vert }^2}{c^2{\Vert g_k\Vert }^4}\\ \le \frac{\left[L{\gamma }_2+3\left(L+L_0\right)\right]\Vert s_k\Vert }{c{\gamma }_1^2}+\frac{\left[\mu {\gamma }_{2}M{\left(L+3\left(L+L_0\right)\right)}^2{M}_1\right]\Vert s_k\Vert }{c^2{\gamma }_1^4}\\ =\left(\frac{\left[L{\gamma }_2+3\left(L+L_0\right)\right]c{\gamma }_1^2+[\mu {\gamma }_{2}M{\left(L+3\left(L+L_0\right)\right)}^2{M}_1]}{c^2{\gamma }_1^4}\right)\Vert s_k\Vert ,\end{array}$

取 $b=\max \left\{2,\frac{\left[L{\gamma }_2+3\left(L+L_0\right)\right]c{\gamma }_1^2+\left[\mu {\gamma }_{2}M{\left(L+3\left(L+L_0\right)\right)}^2{M}_1\right]}{c^2{\gamma }_1^4}\right\}$ 和

$\lambda =\frac{c^2{\gamma }_1^4}{b\left\{\left[L{\gamma }_2+3\left(L+L_0\right)\right]c{\gamma }_1^2+\left[\mu {\gamma }_{2}M{\left(L+3\left(L+L_0\right)\right)}^2{M}_1\right]\right\}}.$

利用(3.8)， $\lambda$ 和 $b>1$ ，则 $\left|{\beta }_k^{MMLS*}\right|\le b$ 和

$\begin{array}{c}\left|{\beta }_k^{MMLS*}\right|\le \frac{\left[L{\gamma }_2+3\left(L+L_0\right)\right]c{\gamma }_1^2+\left[\mu {\gamma }_{2}M{\left(L+3\left(L+L_0\right)\right)}^2{M}_1\right]}{c^2{\gamma }_1^4}\Vert s_k\Vert \\ \le \frac{\left[L{\gamma }_2+3\left(L+L_0\right)\right]c{\gamma }_1^2+\left[\mu {\gamma }_{2}M{\left(L+3\left(L+L_0\right)\right)}^2{M}_1\right]}{c^2{\gamma }_1^4}\lambda =\frac{1}{2b}.\end{array}$

因此公式(2.3)满足性质(P)，引理得证。

利用假设A，引理3.1~3.3，与文献 [15] 中的定理3.2的证明类似，可得到算法1的全局收敛性定理，定理的证明请参照文献 [15] 中的定理3.2完成。

定理3.1：序列 $\left\{d_k,g_k\right\}$ 由算法产生，若假设A成立，则 ${\lim }_{k\to \infty }\inf \Vert g_k\Vert =0$ 成立。

4. 数值结果

这部分将给出数值检验，Benchmar问题 [16] 在工程领域具有广泛应用背景。

1) Sphere函数： $f_{Sph}\left(x\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{x}_i^2},x_i\in \left[-5.12,5.12\right],f_{Sph}\left(x^{*}\right)=0.$

2) Schwefel’s函数： $f_{SchDS}\left(x\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left({\displaystyle \underset{j=1}{\overset{i}{\sum }}{x}_j}\right)}^2},x_i\in \left[-65.536,65.536\right],f_{SchDS}\left(x^{*}\right)=0.$

3) Rastrigin函数： $f_{Ras}\left(x\right)=10n+{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(x_i^2-10\cos \left(2\text{π}{x}_i\right)\right)},x_i\in \left[-5.12,5.12\right],f_{Ras}\left(x^{*}\right)=0.$

4) Griewank函数： $f_{Gri}\left(x\right)=1+{\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{x_i^2}{4000}}-{\displaystyle \underset{i=1}{\overset{n}{\prod }}\cos \frac{x_i}{i}},x_i\in \left[-600,600\right],f_{Gri}\left(x^{*}\right)=0.$

参数选取如下： $\epsilon ={10}^{-5},\delta =0.1,\sigma =0.9$ ，停止准则采用Himmeblau准则：如果 $\left|f\left(x_k\right)\right|>e_1$ ，令 $\text{stop}1=\frac{\left|f\left(x_k\right)-f\left(x_{k+1}\right)\right|}{\left|f\left(x_k\right)\right|}$ ；否则，令 $\text{stop}1=\left|f\left(x_k\right)-f\left(x_{k+1}\right)\right|$ 。如果 $\Vert g\left(x\right)\Vert <\epsilon$ 或者 $\text{stop}1<e_2$ 满足，算法终止， $e_1=e_2={10}^{-5}$ 。若迭代次数超过1000次，程序也将终止。为比较算法1的有效性，也给出通常LS方法的数值检验，并进行比较。以下是各代码含义：

$x_0$ ：初始点；Dim：问题的维数；NI：迭代次数；NFG：函数值与梯度值迭代次数和；Time：以秒为单位的计算机CPU时间； $f\left(x\right)$ ：算法终止时的函数值。

为比较算法1与原LS算法的效率，我们采用下述技术即两算法的所有NFG的乘积的比值再开二十四分之一(结果的总个数)次方，以通常LS算法作为底，定义为1，比值见表2，此技术见文献 [17] 。

Table 1. Numerical results

表1. 数值结果

Table 2. The Efficiency of NFG

表2. NFG的效率

从表1的结果可看，算法1和LS算法对上述4个问题都能有效地求解，迭代次数可以接受，最终函数值很接近最优函数值，新算法的CPU时间相对少一些，更具竞争力。表2的效率表明，算法1相对于原LS算法，效率提高2%。

基金项目

2017年国家级大学生创新创业训练计划项目(201710606041)，广西自然科学基金(2018JJA110039)和玉林师范学院校级科研项目(2019YJKY16)。

NOTES

^*第一作者。

^#通讯作者。