共轭梯度法是求解大规模光滑无约束优化问题$\min\{f(x)|\ x\in R^n\}$的有效方法之一, 其具有算法结构简单, 无需矩阵存储和二次终止性等优点, 特别适用于求解大规模优化问题, 其迭代公式一般形如

$ \begin{eqnarray} \label{1.1} x_{k+1}=x_k+\alpha_k d_k, \end{eqnarray} $

(1.1)

$ \begin{eqnarray}\label{1.2}d_k=\left\{ \begin{array} {lll}-g_k, &k=1;\\ -g_k+\beta_kd_{k-1}, &k\geq2, \end{array}\right. \end{eqnarray} $

(1.2)

其中$g_k=\nabla f(x_k)$, $\alpha_k$为搜索步长, $d_k$为搜索方向, $\beta_k$为搜索方向$d_k$的调控参数. 不同的$\beta_k$对应不同的共轭梯度法, 著名的$\beta_k$计算公式有:

$\begin{eqnarray*} \beta_k^{HS}&=&\frac{{g_k^T}(g_k-g_{k-1})}{{d_{k-1}^T}(g_k-g_{k-1})}, \ \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}, \, \, \, \, \ \beta_k^{DY}=\frac{\|g_k\|^2}{{d_{k-1}^T}(g_k-g_{k-1})}. \end{eqnarray*} $

与上述四个公式相对应的共轭梯度法分别称为HS方法[1]、FR方法[2]、PRP方法[3, 4]和DY方法[5]. 其中, FR方法和DY方法具有良好的收敛性质, 而HS方法和PRP方法则数值性能优良. 因此, 为寻求收敛性和数值效果都理想的方法, 基于上述调控参数所建立的共轭梯度法已被广泛研究，且获得很多理论性质和数值效果均优良的改进共轭梯度法[6-14]. 如, 文[9]中的算法改进了DY公式

$ \begin{eqnarray}\label{s1.3}\displaystyle \beta_k^{JMJ}=\frac{{\|g_k\|^2}-\frac{\|g_k\|}{{\|d_{k-1}\|}}|g_k^Td_{k-1}|}{{d_{k-1}^T}(g_k-g_{k-1})}, \end{eqnarray} $

(1.3)

在标准Wolfe线搜索准则下获得了算法的全局收敛性, 且具有较为理想的数值效果. 基于文[9]的思想, 文[10]修正了FR共轭梯度参数公式

$ \begin{eqnarray}\label{s1.4}\displaystyle \beta_k^{N}=\frac{{g_k^T}{(g_k-\frac{\|g_k\|}{\|d_{k-1}\|}d_{k-1})}}{{\|g_{k-1}\|}^2}, \end{eqnarray} $

(1.4)

由式(1.4)产生的搜索方向在强Wolfe线搜索准则下满足充分下降条件, 然而由此共轭梯度参数公式所建立的算法未能获得全局收敛性. 鉴于此, 受文[9, 10]思想的启发, 本文将公式(1.4)进一步修正为：

$ \begin{eqnarray}\label{s1.5} \beta_k^{MJJ}=\frac{{\|g_k\|^2}-\frac{({g_k^Td_{k-1}})^2}{\|d_{k-1}\|^2}}{{\|g_{k-1}\|}^2+u\max\{|g_k^Td_{k-1}|, |g_k^Tg_{k-1}|\}}, \end{eqnarray} $

(1.5)

其中参数$u>1$. 显然, 在精确线搜索条件下, 且目标函数为严格凸二次函数, 则式(1.5)退化为FR共轭梯度参数公式.

本文余下部分组织如下: 基于公式(1.5), 第二部分建立了算法框架, 不依赖于任何线搜索条件, 以新公式为方向调控参数所产生的搜索方向总是充分下降的. 第三部分在标准Wolfe非精确线搜索条件下, 获得了算法的全局收敛性; 最后对所提出的算法进行数值测试, 并与同类算法进行比对.

基于本文所改进的共轭参数公式(1.5), 建立算法框架A:

初始化  任取初始点${x_1}\in{R^n}$, 给定精度$\varepsilon>0$, 参数$u>1$, $d_1=-g_1, k=1$.

步骤1  若${\|g_k\|} < \varepsilon$, 则停止. 否则, 转步骤2.

步骤2   由某个非精确线搜索准则求步长$\alpha_k$.

步骤3   令${x_{k+1}={x_k}+{\alpha_k}{d_k}}$, 由式(1.5)计算${\beta_{k+1}^{MJJ}}$. 令$d_{k+1}=-g_{k+1}+\beta_{k+1}^{MJJ}d_k, \ k: =k+1$, 返回步骤1.

首先, 给出共轭参数公式(1.5)和算法框架A所产生的搜索方向具有如下的性质.

引理2.1   对于任意$k\geq1$, 总有$0\leq\beta_k^{MJJ}\leq\beta_k^{FR}$恒成立.

证  设$\theta_k$为$g_k^T$和$d_{k-1}$的夹角, 有

$ {{|{cos\theta_k}|}^2}=\frac{|{g_k^T}{d_{k-1}}|^2}{{\|g_k\|}^2{\|d_{k-1}\|}^2}. $

于是

$ \beta_k^{MJJ} %{\displaystyle =\frac{{\|g_k\|^2}-\frac{({g_k^Td_{k-1}})^2}{\|d_{k-1}\|^2}}{{\|g_{k-1}\|}^2+u\max\{|g_k^Td_{k-1}|, |g_k^Tg_{k-1}|\}}} {\displaystyle=\frac{{\|g_k\|^2}{[1-\frac{({g_k^Td_{k-1}})^2}{{\|g_k\|}^2\|d_{k-1}\|^2}]}}{{\|g_{k-1}\|}^2+u\max\{|g_k^Td_{k-1}|, |g_k^Tg_{k-1}|\}}} {\displaystyle =\frac{{\|g_k\|^2}{(1-{|cos\theta_k|^2})}}{{\|g_{k-1}\|}^2+u\max\{|g_k^Td_{k-1}|, |g_k^Tg_{k-1}|\}}, } $

进而, 可得

$ 0\leq\beta_k^{MJJ}\leq\frac{\|g_k\|^2}{{\|g_{k-1}\|}^2+u\max\{|g_k^Td_{k-1}|, |g_k^Tg_{k-1}|\}}\leq\frac{\|g_k\|^2}{\|g_{k-1}\|^2}=\beta_k^{FR}. $

故引理获证.

引理2.2  对所有的$k\geq1$和$u>1$, 由算法框架A所产生的搜索方向$d_k$满足

$g_k^Td_k\leq-(1-\frac{1}{u}){\|g_k\|}^2, $

即方向$d_k$是充分下降的.

证  由$\beta_k^{MJJ}$的定义, 易知

$ \begin{eqnarray} \label{2.1}0\leq\beta_k^{MJJ}=\frac{{\|g_k\|^2}-\frac{({g_k^Td_{k-1}})^2}{\|d_{k-1}\|^2}}{{\|g_{k-1}\|}^2+u\max\{|g_k^Td_{k-1}|, |g_k^Tg_{k-1}|\}}\leq\frac{\|g_k\|^2}{u\max\{|g_k^Td_{k-1}|, |g_k^Tg_{k-1}|\}}. \end{eqnarray} $

(2.1)

将(1.2)两端与$g_k$作内积, 利用(2.1), 得

$ \begin{eqnarray*} g_k^Td_k&&{={g_k^T}({-g_k}+\beta_k^{MJJ}d_{k-1})} {\leq{{-\|g_k\|}^2}+\beta_k^{MJJ}|g_k^Td_{k-1}|} \\&&{\displaystyle\leq{{-\|g_k\|}^2}+\frac{\|g_k\|^2}{u\max\{|g_k^Td_{k-1}|, |g_k^Tg_{k-1}|\}}|g_k^Td_{k-1}|}\\&&{\displaystyle\leq-(1-\frac{1}{u}){\|g_k\|}^2.} \end{eqnarray*} $

因此, 算法框架A所产生的搜索方向是充分下降的. 故引理获证.

由引理2.2可知, 不依赖于任何线搜索条件, 以新公式为方向调控参数所产生的搜索方向总是充分下降的.

为了便于后续的描述, 算法框架A采用如下标准$Wolfe$非精确线搜索准则求步长所产生的算法称为算法MJJ, 即步长$\alpha_k$满足:

$ f(x_k+\alpha_kd_k)\leqslant{f(x_k)+\delta\alpha_kg_k^Td_k}, $

(2.2)

$g(x_k+\alpha_kd_k)^Td_k\geqslant\sigma g_k^Td_k, $

(2.3)

其中参数$0 < \delta < \sigma < 1$.

为获得算法MJJ的全局收敛性, 首先给出算法所需的两个常规假设条件.

(A1)  目标函数$f(x)$在水平集$\Lambda=\{x\in R^n|f(x)\leqslant f(x_1)\}$上有下界, 其中$x_1$为算法的初始点.

(A2)  目标函数$f(x)$在水平集$\Lambda$的某邻域$U$内可微, 其梯度函数$g(x)=\nabla f(x)$满足Lipschitz条件, 即存在常数$L>0$, 使$\|g(x)-g(y)\|\leq L\|x-y\|, \forall\ x, y\in U$.

下面引理给出著名的Zoutendijk条件, 其在共轭梯度法全局收敛性分析中起着非常重要的作用.

引理3.1  设假设A1和A2成立, 考虑一般迭代方法(1.1)-(1.2), 若搜索方向$d_k$满足$g_k^Td_k < 0$, 步长因子$\alpha_k$满足标准Wolfe线搜索准则, 则$\sum\limits_{k=1}^\infty\frac{(g_k^Td_k)^2}{{\|d_k\|}^2} < \infty$成立. 特别地, 若方向$d_k$是充分下降的, 有$\sum\limits_{k=1}^\infty\frac{{\|g_k\|}^4}{{\|d_k\|}^2} < \infty$成立.

由引理2.1, 引理2.2和引理3.1可建立算法MJJ的全局收敛性定理.

定理3.2  设假设A1和A2成立, 点列$\{x_k\}$由(1.1)和(1.2)产生, 步长$\alpha_k$满足标准$Wolfe$非精确线搜索准则(2.2)-(2.3), 则$\lim\limits_{k\rightarrow\infty}\inf\|g_k\|=0$, 即算法MJJ是全局收敛的.

证  由反证法. 若定理3.2不成立, 不妨假设存在常数$\varepsilon>0$, 使得$\|g_k\|^2\geq\varepsilon, $$\ \forall\ k\geq0$. 将式(1.2)两端取平方, 得

$\begin{eqnarray}\label{3.1}\|d_k\|^2=\|g_k\|^2-2\beta_k^{MJJ}g_k^Td_{k-1}+(\beta_k^{MJJ})^2\|d_{k-1}\|^2.\end{eqnarray} $

(3.1)

由式(2.1)可得

$ -2\beta_k^{MJJ}g_k^Td_{k-1}\leq2\mid\beta_k^{MJJ}\mid\mid g_k^Td_{k-1}\mid\leq\frac{2\|g_k\|^2}{u\max\{|g_k^Td_{k-1}|, |g_k^Tg_{k-1}|\}}\mid g_k^Td_{k-1}\mid\leq\frac{2\|g_k\|^2}{u}. $

于是, 将上式代入(3.1), 并利用引理2.1, 有

$ \begin{eqnarray}\label{3.2}{\displaystyle\|d_k\|^2\leq\|g_k\|^2+\frac{2\|g_k\|^2}{u}+(\beta_k^{FR})^2\|d_{k-1}\|^2=\|g_k\|^2(1+\frac{2}{u})+\frac{\|g_k\|^4}{\|g_{k-1}\|^4}\|d_{k-1}\|^2.}\end{eqnarray} $

(3.2)

将式(3.2)两端同除以$\|g_k\|^4$, 得

$ \begin{eqnarray}\label{3.3}\frac{\|d_k\|^2}{\|g_k\|^4}\leq\frac{1}{\|g_k\|^2}(1+\frac{2}{u})+\frac{\|d_{k-1}\|^2}{\|g_{k-1}\|^4}.\end{eqnarray} $

(3.3)

注意到式(1.2)及$\|d_1\|^2=-g_1^Td_1=\|g_1\|^2$, 由(3.3)立得

$ \frac{\|d_k\|^2}{\|g_k\|^4}\leq(1+\frac{2}{u})\sum\limits_{i=1}^k\frac{1}{\|g_i\|^2}. $

又由$\|g_k\|^2\geq\varepsilon$, 有$\frac{\|g_k\|^4}{\|d_k\|^2}\geq\frac{\varepsilon}{k}\frac{u}{u+2}.$ 进而, 立得$\sum\limits_{k=1}^\infty\frac{\|g_k\|^4}{\|d_k\|^2}=+\infty$, 这与引理3.1矛盾. 故定理获证.

为检验本文所提出的算法实际数值效果, 数值试验共测试了43个问题, 所有算例均取自于无约束优化测试问题集[15], 并与文[9]中的JMJ方法, 文[10]中的NJJ方法, 文[2]中的FR方法及文[6]中的PRP+方法进行比较. 所有测试都采用标准Wolfe线搜索准则(2.2)和(2.3)获得步长. 测试的环境为MATLAB R2017b, Windows 10操作系统, 计算机硬件为Inter(R) Core(TM) i5-8250U CPU 1.80 GHz和8 GB RAM. 相关参数选取如下: $\delta=0.01, \sigma=0.1, u=2.5$. 算法的终止准则为以下两种情形之一: (1) $\|g_k\| < 10^{-5};$ (2) 迭代次数Itr$>2000$. 终止准则(2)出现时认为该方法对这个数值例子失效, 并记为"F".

在试验中, 对所测试算法的迭代次数(Itr), 目标函数函数值计算次数(NF), 梯度计算次数(NG), 计算时间(Tcpu)(单位为秒)及算法终止时目标函数梯度的2-范数($\|g_*\|$)等5个重要指标进行观测和比较, 数值结果详见表 1. 为进一步直观刻划试验效果, 我们还采用Dolan和Mor$\acute{e}$[16]性能图对试验效果进行比对. 图 1-4分别就算法MJJ、算法JMJ、算法NJJ、算法FR及算法PRP+在迭代次数、函数值计算次数、梯度计算次数及计算时间进行比较. 从图 1-4可以更直观地看出, 算法MJJ在所考查的4个指标中均明显优于其他4个算法. 因此, 算法MJJ不仅在解决问题的个数上占据优势, 并且鲁棒性最优. 综上所述, 本文所提出的算法是有效的.