考虑无约束最优化问题

$ \min\limits_{x\in R^{n}} f(x), $

(1.1)

其中$f(x): R^n\rightarrow R$为连续可微函数.

拟牛顿法是解(1.1) 的最常用的方法之一, 它以拟牛顿方程为基础, 利用目标函数一阶导数的信息, 构造出目标函数的曲率近似, 从而避免了求Hesse阵或其逆矩阵.记$g_k=\nabla f(x_k)$, 其基本迭代格式为

$ x_{k+1}=x_k+\alpha_kd_k, $

(1.2)

其中$\alpha_k$为通过某种搜索策略确定的步长因子, $d_k$为搜索方向$-H_kg_k$, $H_k$为目标函数$f(x)$的Hesse阵的逆的近似, 要求正定且满足拟牛顿方程

$ H_{k+1}y_k=s_k, $

(1.3)

其中$y_k=g_{k+1}-g_k, s_k=x_{k+1}-x_k$.

传统的确定步长因子$\alpha_k$的Armijo线搜索策略要求$\alpha_k$满足

$ f(x_k+\alpha_kd_k)\leq f(x_k)+\gamma\alpha_kg_k^Td_k, $

(1.4)

其中$\gamma\in(0, 1).$显然这样搜索要求每次迭代目标函数值单调下降, 大大降低了算法的效率^[1].非单调线搜索技术由于其有利于求解全局最优解和算法的快速收敛而受到许多最优化爱好者的关注. Grippo等在文献[1]中提出的非单调线搜索技术放宽了对$\alpha_k$的选取范围, 每次迭代仅要求步长因子$\alpha_k$满足下式

$ f(x_k+\alpha_kd_k)\leq f(x_{l(k)})+\gamma\alpha_kg_k^Td_k, $

(1.5)

其中$f(x_{l(k)})=\max\limits_{0\leq j\leq m_k}\{f(x_{k-j})\}, m_0=0, 0\leq m_k\leq \min\{m_{k-1}+1, M\}(k\geq1)$, $M\geq0$为正整数.这种非单调技术在文献[2]中被推广到信赖域算法.但是Grippo等人的非单调技术有一些不足之处, 比如从式(1.5)可见, 由于$f(x_{l(k)})$是前面$M$个函数值里的最大者, 所以可能导致某步迭代得到的较好的函数值信息没有被用上, 另外算法的数值表现受$M$的取值的影响.于是, Zhang等在文献[3]中提出了一种新的非单调线搜索方法, 将式(1.5)中的$f(x_{l(k)})$换成某些函数值的加权平均$C_k$, 即步长因子$\alpha_k$要满足

$ f(x_k+\alpha_kd_k)\leq C_k+\gamma\alpha_kg_k^Td_k, $

(1.6)

其中

$ C_k=\left\{\begin{array}{ll} f(x_k), &k=0, \nonumber\\ \frac{\eta_{k-1}Q_{k-1}C_{k-1}+f(x_k)}{Q_k}, &k\geq1, \nonumber \end{array}\right. \, Q_k=\left\{\begin{array}{ll} 1, &k=0, \nonumber\\ \eta_{k-1}Q_{k-1}+1, &k\geq1, \nonumber \end{array}\right. $

$\eta_{k-1}\in [\eta_{\rm min}, \eta_{\rm max}], \eta_{\rm min}\in[0, 1)$及$\eta_{\rm max}\in [\eta_{\rm min}, 1)$为两个选定的参数.文献[4]中提出了采用该非单调策略的信赖域算法.随后, 文献[5]中又指出上述非单调方法每次迭代要计算$Q_k$以及$C_k$, 这在计算上是一种浪费, 于是提出了一种新的非单调信赖域方法, 将(1.6)式中的$C_k$替换为$D_k$, 这里$D_k$是$D_{k-1}$和$f_k$的一个简单凸组合, 即

$ D_k=\left\{\begin{array}{ll} f_k, &k=1,\\ \eta_kD_{k-1}+(1-\eta_k)f_k, &k\geq2, \end{array}\right. $

(1.7)

其中$\eta_k\in(0, 1)$可以是变量, 也可以取常数.数值试验表明文献[5]中的非单调策略更高效.

通常拟牛顿法每次迭代需要存储一个矩阵$H_k$, 而该矩阵一般是稠密的, 即算法的存储量至少为$O(n^2)$, 对中小规模的优化问题有较好的数值表现.一般来说, 大规模优化问题的目标函数Hesse阵是对称的, 且具有某种稀疏结构.经典拟牛顿法(如DFP方法、BFGS方法), 虽然校正矩阵能继承正定性和对称性, 却不能继承其稀疏性, 存储空间问题成为阻扰拟牛顿方法求解大规模问题的一大障碍.因此, 希望得到存储量小且收敛速度快、计算量小的算法, 需对传统的拟牛顿法进行改进.占据主要内存单元的Hesse阵是重点考虑对象, 若限制拟牛顿法中Hesse阵或其逆的近似为对角矩阵则可以大大降低存储量, 在一定的条件下也能保持相对较快的收敛速度.

1988年, Barzilai和Borwein^[6]提出的两点步长梯度法(BB方法), 利用当前迭代点以及前一点的信息来确定梯度方向的步长, 也可看成用具有一定拟牛顿性质的数量矩阵$H_{k+1}=\alpha_{k+1}I$作为Hesse阵的逆的近似, 其中$\alpha_{k+1}$为$\alpha_{k+1}=\displaystyle{\mbox{argmin}_{\alpha}}\|\alpha y_k-s_k\|^2$, 即

$ \alpha_{k+1}=\frac{s_k^Ty_k}{y_k^Ty_k}. $

BB方法对二维的二次严格凸的目标函数是超线性收敛的, 对于目标函数是一般的二次凸函数则是R -线性收敛的, 但不能保证全局收敛性.为克服这个缺点, 同时保留BB步长的性质, Raydan^[7]将BB方法与Grippo等人的非单调线搜索方法^[1]结合, 构造了一个具有全局收敛性的非单调线搜索算法(GBB方法).数值试验结果表明GBB方法对于求解大规模无约束优化问题的数值表现较好. 2006年时贞军^[8]受BB方法的启发, 在拟牛顿法中, 限制Hesse阵逆阵的近似$H_{k+1}$为正定对角矩阵且近似满足经典的拟牛顿方程(1.3), 通过求解

$ {\min}_{\underline{h}\leq h_{k+1}^{(i)}\leq \overline{h}}\|H_{k+1}y_k-s_k\|^2 $

而得到$H_{k+1}$(其中$\underline{h}$和$\overline{h}$为两个正数), 即$H_{k+1}$的第$i$个对角元

$ h_{k+1}^{(i)}=\left\{\begin{array}{llll} \frac{s_k^{(i)}}{y_k^{(i)}}, \quad\mbox{若}\quad \underline{h}\leq\frac{s_k^{(i)}}{y_k^{(i)}}\leq \overline{h};\\ \underline{h}, \quad\mbox{若}\quad\frac{s_k^{(i)}}{y_k^{(i)}}<\underline{h};\\ \overline{h}, \quad\mbox{若} \quad \frac{s_k^{(i)}}{y_k^{(i)}}>\overline{h};\\ h_k^i, \quad\mbox{若}\quad y_k^{(i)}=0.\end{array}\right. $

(1.8)

从而构造了对角稀疏拟牛顿法, 该方法为求解大规模问题提供了新的思路.文献[9,10]在此基础上设计了求解无约束最优化问题的非单调对角稀疏拟牛顿算法.

以上的对角拟牛顿法都是建立在经典拟牛顿方程(1.3) 的基础之上.近年来, 基于修正拟牛顿方程的修正拟牛顿法的研究吸引了不少国内外学者. 1991年, Yuan^[11]提出了一种修正BFGS方法, 随后, Zhang和Xu^[12], Wei、Li和Qi^[13]等相继又提出了几类修正拟牛顿算法.理论分析和数值试验表明, 这些改进后的拟牛顿法由于利用了迭代点的更多信息, 所以比经典拟牛顿法的数值表现更好.为获取Hesse阵的逆矩阵的更良好的对角逼近, 本文以逼近程度更好的修正拟牛顿方程为基础, 推导出一类新的对角校正公式, 并将文献[5]中的非单调策略应用于线搜索算法框架, 设计了一种求解无约束优化问题的广义对角拟牛顿算法.

考虑方程

$ B_{k+1}s_k=t_ky_k\, (t_k>0), $

(2.1)

其中$B_{k+1}$为Hesse阵的近似矩阵.当$t_k=1$时, 式(2.1) 为经典拟牛顿方程.当$t_k$选取适当的值时, 可以得到相关的修正拟牛顿方程^[11], 比如可取

$ t_k=t_k^*=\frac{2}{s_k^Ty_k}(f_k-f_{k+1}+s_k^Tg_{k+1}) $

或

$ t_k=t_k^{**}=\frac{1}{s_k^Ty_k}[(g_{k+1}-g_k)^Ts_k+6(f_k-f_{k+1})+3(g_k+g_{k+1})^Ts_k]. $

易证得当$t_k=t_k^*$或$t_k=t_k^{**}$时, 以式(2.1) 为基础得到的$B_{k+1}$较$t_k=1$时得到的$B_{k+1}$在某种意义上更逼近精确的Hesse阵(见文献[12, 13]).

由式(2.1) 得, Hesse阵的近似逆阵$H_{k+1}$需满足如下修正拟牛顿方程

$ H_{k+1}y_k=\frac{1}{t_k}s_k. $

(2.2)

本文中要求正定对角阵近似满足相应的弱拟牛顿方程^[14]

$ y_k^TH_{k+1}y_k=\frac{1}{t_k}y_k^Ts_k\equiv\rho_k. $

(2.3)

事实上, 满足式(2.3)的对角矩阵一般有无数多个, 在其中选取一个最近似满足经典拟牛顿方程(1.3)的一个, 即有如下子问题

$ \mbox{min}\, \frac{1}{2}\|H_{k+1}y_k-s_k\|^2, \nonumber\\ \mbox{s.t.}\, \, y_k^TH_{k+1}y_k=\rho_k, \\ H_{k+1} \mbox{为对角矩阵, }\nonumber $

(2.4)

其中$\|\cdot\|$为欧氏范数.然而, 问题(2.4)的解未必正定, 为保证其正定性, 需进一步限制$h_{k+1}^{(i)}$在一定的范围内, 即要求$0 < \underline{h}\leq h_{k+1}^{(i)}\leq\overline{h}.$下面给出新的正定对角阵的具体构造方法.

问题(2.4) 的目标函数是凸函数, 且约束集为凸集, 其Lagrangian函数为

$ L(H_{k+1}, \mu)=\frac{1}{2}\|H_{k+1}y_k-s_k\|^2+\mu(\rho_k-y_k^TH_{k+1}y_k), $

(2.5)

其中$\mu$为对应于约束条件的Lagrangian乘子.将$L$关于$H_{k+1}$的每个对角元求偏导数, 并置其为0, 得

$ h_{k+1}^{(i)}(y_k^{(i)})^2=\mu (y_k^{(i)})^2+ s_k^{(i)}y_k^{(i)} \, \, (i=1, 2, \cdots, n). $

(2.6)

将式(2.6) 中各式相加, 结合约束条件$y_k^TH_{k+1}y_k=\rho_k$, 得

$ \mu=\frac{\rho_k-s_k^Ty_k}{y_k^Ty_k}. $

(2.7)

将式(2.7) 代入式(2.6)得, 当$y_k^{(i)}\neq 0$时,

$ h_{k+1}^{(i)}=\frac{\rho_k-s_k^Ty_k}{y_k^Ty_k}+\frac{s_k^{(i)}}{y_k^{(i)}}. $

(2.8)

若记$\overline{h}_{k+1}^{(i)}=\frac{\rho_k-s_k^Ty_k}{y_k^Ty_k}+\frac{s_k^{(i)}}{y_k^{(i)}}$, 则

$ h_{k+1}^{(i)}=\left\{\begin{array}{llll} \overline{h}_{k+1}^{(i)}, &\mbox{若}\quad \underline{h}\leq\overline{h}_{k+1}^{(i)}\leq \overline{h};\\ \underline{h}, &\mbox{若}\quad\overline{h}_{k+1}^{(i)}<\underline{h};\\ \overline{h}, &\mbox{若} \quad\overline{h}_{k+1}^{(i)}>\overline{h};\\ h_k^i, &\mbox{若}\quad y_k^{(i)}=0.\end{array}\right. $

(2.9)

这样便得到了近似满足弱拟牛顿方程(2.3) 的正定对角矩阵$H_{k+1}$.显然, 若$t_k=1$, 即$\rho_k=s_k^Ty_k$时, 校正公式(2.9) 退化为文献[8]中的公式(1.8).

接下来给出求解无约束优化问题的非单调广义对角拟牛顿算法.

非单调广义对角拟牛顿算法(GDQN)

步 0  给出初始点$x_0\in R^n, 0 < \gamma < 1, 0 < \beta < 1, 0 < \underline{h} < \overline{h}, 0\leq\eta_{\rm min} < \eta_{\rm max} < 1.$置$k=0, H_0=I, D_0=f(x_0)$.

步 1  若满足算法终止准则, 则停止, 否则, 转步2.

步 2  计算$d_{k}=-H_kg_k$.

步 3  选取$\alpha_k$为$\{1, \beta, \beta^2, \cdots\}$中满足下式的最大者

$ f(x_k+\alpha d_{k})\leq D_k+\gamma\alpha g_k^Td_k, $

(2.10)

令$x_{k+1}=x_k+\alpha_kd_k$.

步 4  选取$\eta_{k+1}\in [\eta_{\rm min}, \eta_{\rm max}]$, 按(1.7) 式计算$D_{k+1}$.

步 5  计算$s_k$, $y_k$及$\rho_k$, 按(2.9) 式校正$H_{k+1}$, 令$k=k+1, $转步1.

注  (1) 若$\eta_k$恒取0且$\rho_k=s_k^Ty_k$, 算法GDQN即为文献[8]中的单调对角稀疏拟牛顿算法.

(2) 按照文献[11]中的修正拟牛顿方程的构造方法, $\rho_k$的选取可按如下两种方式选取:

(a) $t_k=\frac{2}{s_k^Ty_k}(f_k-f_{k+1}+s_k^Tg_{k+1})$, 即

$ \rho_k=\frac{(s_k^Ty_k)^2}{2(f_k-f_{k+1}+s_k^Tg_{k+1})}. $

(2.11)

(b) $t_k=\frac{1}{s_k^Ty_k}[(g_{k+1}-g_k)^Ts_k+6(f_k-f_{k+1})+3(g_k+g_{k+1})^Ts_k]$, 即

$ \rho_k=\frac{(s_k^Ty_k)^2}{(g_{k+1}-g_k)^Ts_k+6(f_k-f_{k+1})+3(g_k+g_{k+1})^Ts_k}. $

(2.12)

对一般的函数, 由式(2.11) 或(2.12) 定义的$\rho_k$可能是负的, 甚至可能溢出, 为解决这一问题, 按如下安全方式截取$\rho_k$,

$ \rho_k=\left\{\begin{array}{llll} \underline{h}\|y_k\|^2, &\mbox{如果}\quad \rho_k<\underline{h}\|y_k\|^2;\\ \overline{h}\|y_k\|^2, &\mbox{如果}\quad \rho_k>\overline{h}\|y_k\|^2.\end{array}\right. $

本节分析算法GDQN的收敛性, 首先作出如下假设:

假设 3.1  (1) 水平集$L(x_0)=\{x\, |\, f(x)\leq f(x_0)\}$有界; (2) $f(x)$在$L(x_0)$上二次连续可微.

引理 3.2  若$x_k$不是函数$f(x)$的稳定点, 则(1) $\|d_k\|\leq\overline{h}\|g_k\|$; (2) $g_k^Td_k\leq-\underline{h}\|g_k\|^2$.

证  由$d_k$的定义和$H_k$的取法易知.

引理 3.3  假设$\{x_k\}$是算法GDQN产生的序列, 则对任何$k$, 有

$ D_{k+1}\geq f_{k+1}. $

(3.1)

证  由$D_{k+1}$的定义得

$ D_{k+1}-f_{k+1}=\eta_{k+1}(D_k-f_{k+1}). $

(3.2)

由引理3.2及线搜索准则(2.10) 易得$D_k-f_{k+1}\geq -\gamma\alpha_kg_k^Td_k\geq0.$因此由式(3.2) 即得式(3.1) 成立.

引理 3.4  算法GDQN中步3中的线搜索必在有限步内终止, 即算法GDQN是适定的.

证  用反证法.假设结论不成立, 即存在某个$k$, 使得对于任何正整数$i$, 有

$ f(x_k+\beta^id_k)>D_k+\gamma\beta^ig_k^Td_k. $

(3.3)

由引理3.3得$f_k\leq D_k$.于是结合式(3.3) 得

$ \frac{f(x_k+\beta^id_k)-f_k}{\beta^i}> \gamma g_k^Td_k. $

(3.4)

由于$f(x)$二次连续可微, 令$i\rightarrow\infty$, 将式(3.4) 两边取极限得$g_k^Td_k\geq\gamma g_k^Td_k.$该式表明$g_k^Td_k\geq 0$, 与引理3.2矛盾.

引理 3.5  假设$\{x_k\}$是算法GDQN产生的序列, 则对任何$k$, 存在$M>0$, 使得步长因子$\alpha_k$满足

$ \alpha_k>\min\{\frac{\beta}{2}, \frac{2\beta(\gamma-1)g_k^Td_k}{M\|d_k\|^2}\}. $

(3.5)

证  若$\alpha_k/\beta>\frac{1}{2}$, 则$\alpha_k>\frac{\beta}{2}$.以下考虑$\alpha_k/\beta\leq\frac{1}{2}$的情形.由算法GDQN的步3得

$ f(x_k+\frac{\alpha_k}{\beta}d_k)>D_k+\gamma\frac{\alpha_k}{\beta}g_k^Td_k. $

(3.6)

由Taylor展开式得

$ f(x_k+\frac{\alpha_k}{\beta}d_k)=f_k+\frac{\alpha_k}{\beta}g_k^Td_k+\frac{1}{2}(\frac{\alpha_k}{\beta})^2d_k^T\nabla^2f(\xi_k)d_k, $

(3.7)

其中$\xi_k\in(x_k, x_k+\frac{\alpha_k}{\beta}d_k).$

由于$f(x)$在有界闭集$L(x_0)$上二次连续可微, 所以存在常数$M>0$, 使得$\|\nabla^2f(x)\|\leq M$.结合式(3.1), (3.6), (3.7) 以及假设3.1, 可得

$ f_k+\gamma\frac{\alpha_k}{\beta}g_k^Td_k<f_k+\frac{\alpha_k}{\beta}g_k^Td_k+\frac{1}{2}(\frac{\alpha_k}{\beta})^2d_k^T\nabla^2f(\xi_k)d_k\leq f_k+\frac{\alpha_k}{\beta}g_k^Td_k+\frac{1}{2}(\frac{\alpha_k}{\beta})^2M\|d_k\|^2. $

因此

$ \alpha_k>\frac{2\beta(\gamma-1)g_k^Td_k}{M\|d_k\|^2}. $

于是式(3.5) 得证.

引理 3.6  假设$\{x_k\}$是算法GDQN产生的序列, 则序列$\{D_k\}$单调不增.

证  由线搜索准则(2.10)、引理3.2和引理3.5得

$ f_{k+1}\quad\leq\quad D_k+\gamma\alpha_kg_k^Td_k\leq D_k-\gamma \underline{h}\min\{\frac{\beta}{2}, \frac{2\beta(\gamma-1)g_k^Td_k}{M\|d_k\|^2}\}\|g_k\|^2 \nonumber\\ \quad\quad\quad\leq\quad D_k-\gamma \min\{\frac{\beta\underline{h}}{2}, \frac{2\beta(1-\gamma)\underline{h}^2}{M \overline{h}^2}\}\|g_k\|^2 \equiv D_k-\delta \|g_k\|^2, $

(3.8)

其中常数$\delta=\gamma \min\{\frac{\beta\underline{h}}{2}, \frac{2\beta(1-\gamma)\underline{h}^2}{M \overline{h}^2}\}$.根据$D_{k+1}$的定义及式(3.8) 有

$ D_{k+1}\quad=\quad\eta_{k+1}D_k+(1-\eta_{k+1})f_{k+1}\nonumber\\ \quad\quad\quad\leq\quad\eta_{k+1}D_k+(1-\eta_{k+1})(D_k-\delta \|g_k\|^2)\nonumber\\ \quad\quad\quad=\quad D_k-(1-\eta_{k+1})\delta \|g_k\|^2. $

(3.9)

(3.9) 式表明序列$\{D_k\}$单调不增.

引理 3.7  若假设3.1成立, $\{x_k\}$是算法GDQN产生的序列, 则序列$\{D_k\}$收敛.

证  由假设3.1, 引理3.3, 引理3.6及$D_0=f(x_0)$易得算法GDQN产生的序列$\{x_k\}$包含在有界水平集$L(x_0)$中, 于是进一步可得$\{D_k\}$收敛.

定理 3.8  若假设3.1成立, 算法GDQN产生的无穷点列为$\{x_k\}$, 则

$ \mathop {\liminf}\limits_{k\rightarrow\infty}\|g_k\|=0. $

(3.10)

证  用反证法.假设式(3.10)不成立, 则存在$\varepsilon>0$使得对任何$k$, 有

$ \|g_k\|\geq\varepsilon. $

(3.11)

结合式(3.9) 和(3.11) 得

$ D_k-D_{k+1}\geq (1-\eta_{k+1})\delta \|g_k\|^2\geq (1-\eta_{k+1})\delta \varepsilon^2. $

(3.12)

于是

$ \sum\limits_{k=1}^\infty (D_k-D_{k+1})\geq \sum\limits_{k=1}^\infty (1-\eta_{k+1})\delta \varepsilon^2. $

(3.13)

由于$\{D_k\}$收敛, 因此

$ +\infty>\sum\limits_{k=1}^\infty (1-\eta_{k+1})\delta \varepsilon^2 \geq \sum\limits_{k=1}^\infty (1-\eta_{\rm max})\delta \varepsilon^2. $

(3.14)

而事实上$\displaystyle{\sum\limits_{k=1}^\infty} (1-\eta_{\rm max})\delta \varepsilon^2$发散, 这与(3.14)式矛盾, 因此该定理成立.

为验证算法GDQN的有效性, 本节选用文献[15]中的一些测试函数, 利用MATLAB编制程序进行数值实验, 并与采用文献[8]中的校正方法对应的算法进行比较.

算法GDQN中的有关参数选取如下: $\gamma=0.0001, \beta=0.5, \eta_k$取常数0.5, 选取

$ \underline{h}\quad=\quad\max\{0.5\cdot\frac{|(x_k-x_{k-1})^T(g_k-g_{k-1}))|}{\|g_k-g_{k-1}\|^2}, 0.0001 \}, \\ \overline{h}\quad=\quad\min\{5\cdot\frac{|(x_k-x_{k-1})^T(g_k-g_{k-1}))|}{\|g_k-g_{k-1}\|^2}, 10000 \}. $

根据算法GDQN中的$\rho_k$的不同取法, 得到如下的三个不同的算法:

(1) DQN: $\rho_k=s_k^Ty_k$, 即采用文献[8]中的校正方法的非单调对角拟牛顿法.

(2) GDQN1: $\rho_k$由式(2.11)确定.

(3) GDQN2: $\rho_k$由式(2.12)确定.

算法的终止准则为$\|g_k\|_\infty\leq 10^{-5}(1+|f(x_k)|).$

当迭代次数超过5000时迭代也终止.计算结果见表 1.表中NF表示函数估值次数, Iter表示迭代次数, CPU表示计算机执行算法所需的时间.

应用Dolan和Moré^[16]的剖面分析法, 进一步借助MATLAB软件分别比较了三个算法的函数估值次数、迭代次数(见图 1)以及CPU时间(见图 2).由图显见, 算法GDQN1和GDQN2的数值表现总体优于算法DQN, 其中算法GDQN2的表现最佳.初步的数值实验表明新的对角校正方法较文献[8]的校正方法有了较大的改进, 提高了对角拟牛顿算法的效率.

本文借助弱拟牛顿方程, 对文献[8]的对角稀疏拟牛顿法进行了推广与改进, 提出了广义对角拟牛顿法, 主要有以下几方面特点:

(1) 为获取更良好的Hesse阵逆矩阵的近似, 合理利用弱拟牛顿方程, 推导出了一类新的Hesse阵逆矩阵的对角校正公式;

(2) 用对角矩阵逼近Hesse阵的逆矩阵, 在计算过程中不用存储和计算矩阵, 有利于大型问题的求解;

(3) 引进了新的非单调技术提高算法的效率.