求解无约束问题的一种新的扰动BFGS方法

A New Perturbed BFGS Method for Unconstrained Optimization Problems

摘要: 基于文献[1]的扰动思想和文献[2]中的BFGS型方法，本文提出了一种新的扰动BFGS方法并证明了其在Wolfe搜索下求解非凸优化问题具有全局收敛性。数值结果表明该方法比较有效。

Abstract: Based on the idea of [1] and the BFGS method in [2], this paper presents a new perturbed BFGS method for the unconstrained optimization. We prove that the proposed method has global con-vergence for nonconvex optimization problems. Numerical results show that this method is effi-cient.

1. 引言

本文考虑如下的无约束问题：

$\min f\left(x\right),x\in R^n,$ (1)

其中 $$ 是连续可微函数。

BFGS方法是求解无约束优化问题的一种非常有效的拟牛顿方法，BFGS方法在Wolfe搜索下对求解凸优化问题具有全局收敛性，但对非凸问题不一定收敛。在文献 [1] 中刘陶文提出了一种扰动的BFGS方法，其中对标准的BFGS方法的迭代矩阵进行了扰动，该方法证明了在Wolfe搜索下对求解非凸问题具有全局收敛性。Zhang，Deng和Chen在文献 [2] 中基于新的拟牛顿矩阵提出了一种新的拟牛顿算法，该方法满足更好的割线方程。本文的目的就是基于文献 [1] 的扰动思想，同时结合文献 [2] 中的BFGS迭代公式，提出了一种新的扰动的BFGS方法。

本文安排如下，在第二节中我们提出新的扰动BFGS方法的算法，在第三节中我们分析了该算法的全局收敛性，第四节中我们进行了数值实验，验证了该算法的有效性，在第五节中我们对论文进行了总结。

2. 算法

求解问题(1)新的扰动BFGS方法基本的迭代格式为

$x_{k+1}=x_k+{\alpha }_k{d}_k,k\ge 1,$ (2)

其中 $$ 是当前的近似解， $d_k$ 为搜索方向， ${\alpha }_k$ 是步长。

基于文献 [1] 的扰动思想，通过求解下列线性方程组得到搜索方向 $d_k$

$\left(B_k+{\mu }_{k}Q\right)d_k+\nabla f\left(x_k\right)=0,$ (3)

其中 $\nabla f\left(x\right)$ 表示其梯度， $B_k$ 为拟牛顿矩阵，也就是 $f\left(x\right)$ 在点x处的Hessian矩阵的近似，并且它要求是正定的。Q是某个给定的对称正定矩阵， $u_k>0$ 且它的值依赖于k，依照某种规则取值为 ${\epsilon }_k$ 或者 ${\epsilon }_k{\Vert B_k\Vert }_F$ ，其中 ${\epsilon }_k>0$ 是某一个参数且它的值与k有关，当 $k\to +\infty$ 时， ${\epsilon }_k\to 0$ 。

步长 ${\alpha }_k$ 是采用Wolfe型非精确搜索确定的，即通过求解

$\{\begin{array}{l}f\left(x_k+{\alpha }_k{d}_k\right)-f\left(x_k\right)\le {\sigma }_1{\alpha }_k\nabla f{\left(x_k\right)}^{\text{T}}{d}_k\\ \nabla f{\left(x_{k+1}\right)}^{\text{T}}{d}_k\ge {\sigma }_2\nabla f{\left(x_k\right)}^{\text{T}}{d}_k\end{array}$ (4)

来计算步长 ${\alpha }_k$ ，则条件 $y_k^{\text{T}}{s}_k>0$ 被满足，其中 ${\sigma }_1$ 和 ${\sigma }_2$ 是正的常数，且满足 ${\sigma }_1<{\sigma }_2<1$ 。

新的拟牛顿矩阵由Zhang，Deng和Chen文献 [2] 提出，由下列BFGS公式得到

$B_{k+1}=B_k-\frac{B_k{s}_k{s}_k^{\text{T}}{B}_k}{s_k^{\text{T}}{B}_k{s}_k}+\frac{{\bar{y}}_k{\bar{y}}_k^{\text{T}}}{{\bar{y}}_k^{\text{T}}{s}_k},$ (5)

其中

${\bar{y}}_k=\nabla f\left(x_{k+1}\right)-\nabla f\left(x_k\right)+r_k{s}_k,$ (6)

$s_k=x_{k+1}-x_k,$ (7)

$r_k=\frac{3{\left(\nabla f\left(x_{k+1}\right)+\nabla f\left(x_k\right)\right)}^{\text{T}}{s}_k-6\left(f\left(x_{k+1}\right)-f\left(x_k\right)\right)}{{\Vert s_k\Vert }^2}.$ (8)

再由Wolfe搜索下则条件被满足，而由Zhang，Deng和Chen文献 [2] 中的引理4可知当 $y_k^{\text{T}}{s}_k>0$ ， ${\bar{y}}^{\text{T}}{}_k{s}_k>0$ 。

下面我们给出求解问题(1)的新的扰动的BFGS方法的具体步骤：

算法1.1：(新的扰动的BFGS方法)

步0：选择初始点 $x_1\in R^n$ 和初始的对称正定矩阵以及。选用常数 $0<{\sigma }_1<{\sigma }_2<1$ ， $\rho$ ， $\tau$ ， $\eta \in \left(0,1\right)$ ， ${\epsilon }_1>0$ 以及 $M_B>1$ 。令 $u_1={\epsilon }_1$ 以及 $\delta =\Vert \nabla f\left(x_1\right)\Vert$ 。令 $k=1$ 。

步1：解线性方程组(3)得到 $d_k$ ，转步2。

步2：由条件(4)来计算 ${\alpha }_k$ ，令 $x_{k+1}=x_k+{\alpha }_k{d}_k$ 。

步3：由修正公式(5)，(6)，(7)，(8)确定 $B_{k+1}$ 。

步4：如果 $\Vert \nabla f\left(x_{k+1}\right)\Vert /\delta \le \eta$ ，令 ${\epsilon }_{k+1}=\tau {\epsilon }_k$ ， $u_{k+1}={\epsilon }_{k+1}$ ， $\delta =\Vert \nabla f\left(x_{k+1}\right)\Vert$ 。

否则令 ${\epsilon }_{k+1}={\epsilon }_k$ 以及

$u_{k+1}=\{\begin{array}{l}{\epsilon }_{k+1}{\Vert B_{k+1}\Vert }_F,如果{\Vert B_{k+1}\Vert }_F\ge M_{k+1}\\ {\epsilon }_{k+1},其他\end{array}$

其中 $M_k=\max \left\{M_B,1/\Vert \nabla f\left(x_k\right)\Vert \right\}$ 。

步5：令 $k=k+1$ ，然后转步1。

在算法1.1中 ${\epsilon }_k$ 按下面的规则逐步被减少：

其中 $\tau$ ， $\eta \in \left(0,1\right)$ 是正的常数，并且 $\delta =\Vert \nabla f\left(x_i\right)\Vert$ ，i取中的某个指标，它的值随着迭代过程被不断的修改，这样得到一个序列 $\left\{\Vert \nabla f\left(x_k\right)\Vert \right\}$ 中的一个公比不超过 $\eta \left(0<\eta <1\right)$ 的几何子列。如果

$\underset{k\to \infty }{\lim }\inf \Vert \nabla f\left(x_k\right)\Vert =0$

成立，则 $\left\{{\epsilon }_k\right\}$ 单调非增并且当 $k\to \infty$ 时 ${\epsilon }_k\to 0$ 。由此可以得到当 ${\Vert B_k\Vert }_F$ 有界时， ${\epsilon }_k{\Vert B_k\Vert }_F\to 0$ 。

3. 算法全局收敛性证明

为了证明新的扰动BFGS方法的全局收敛性，我们做出如下的假设

假设3.1：水平集

$\Omega =\left\{x\in R^n:f\left(x\right)\le f\left(x_1\right)\right\}$

有界，函数 $f\left(x\right)$ 在 $\Omega$ 上Lipschitz连续的梯度，即存在常数 $L>0$ 满足

$\Vert \nabla f\left(x\right)-\nabla f\left(y\right)\Vert \le L\Vert x-y\Vert ,\forall x,y\in \Omega .$ (3.1)

需要注意的是采用线性搜索(4)，由(5)产生的矩阵总是正定的，从而序列 $\left\{f\left(x_k\right)\right\}$ 是单调递减的，而且由线性搜索(4)的第一个不等式以及假设3.1可以得到

$-{\displaystyle \underset{k=1}{\overset{\infty }{\sum }}{\alpha }_k}\nabla f{\left(x_k\right)}^{\text{T}}{d}_k<+\infty ,$

由此可推出

$-\underset{k\to \infty }{\lim }{\alpha }_k\nabla f{\left(x_k\right)}^{\text{T}}{d}_k=-\underset{k\to \infty }{\lim }\nabla f{\left(x_k\right)}^{\text{T}}{s}_k=0.$ (9)

为了表示的方便，根据算法1.1中的步4我们定义指标集合

$K_{\tau }=\left\{k:{\epsilon }_{k+1}=\tau {\epsilon }_k\right\},$ (10)

以及

$K_B=\left\{k:{\Vert B_k\Vert }_F\ge M_k\right\}.$ (11)

为了证明算法1.1的全局收敛性我们需要先证明如下的引理：

引理3.1：假设 $K_{\tau }$ 是有限集，则算法1.1产生的序列 $\left\{{\alpha }_k\right\}$ 有非零的下界，并且当 $k\to \infty$ 时， $d_k\to 0$ 。

证明：如果 $K_{\tau }$ 是有限集，因此 ${\epsilon }_k$ 被减少有限次，也就是说存在一个整数 $k_1$ 使得对所有的 $k\ge k_1$ 有 ${\epsilon }_k={\epsilon }_{k_1}$ ，由于 $B_k$ 是正定的，对于所有的k有

$-\nabla f{\left(x_k\right)}^{\text{T}}{d}_k=d_k^{\text{T}}\left(B_k+u_{k}I\right)d_k\ge {\epsilon }_{k_1}\min \left\{1,M_B\right\}{\Vert d_k\Vert }^2.$ (12)

由(4)中的第二个不等式以及假设3.1，我们得到对所有的 $k\ge 1$ ，不等式

$L{\alpha }_k{\Vert d_k\Vert }^2\ge {\left(\nabla f\left(x_k+{\alpha }_k{d}_k\right)-\nabla f\left(x_k\right)\right)}^{\text{T}}{d}_k\ge -\left(1-{\sigma }_2\right)\nabla f{\left(x_k\right)}^{\text{T}}{d}_k$

成立。因此，由(12)，我们可以推出

$\begin{array}{c}{\alpha }_k\ge \frac{-\left(1-{\sigma }_2\right)\nabla f{\left(x_k\right)}^{\text{T}}{d}_k}{L{\Vert d_k\Vert }^2}=\frac{\left(1-{\sigma }_2\right)d_k^{\text{T}}\left(B_k+u_{k}I\right)d_k}{L{\Vert d_k\Vert }^2}\\ \ge \left(1-{\sigma }_2\right)L^{-1}{\epsilon }_{k_1}\min \left\{1,M_B\right\}>0\end{array}$ (13)

而且由(13)，(9)和(12)易知当 $k\to \infty$ 时， $d_k\to 0$ 。

引理3.2：如果 $K_B$ 是无限集，而 $K_{\tau }$ 是有限集，则 $K_{\tau }$ 必有

$\underset{k\to \infty }{\lim }\inf \Vert \nabla f\left(x_k\right)\Vert =0$ ． (14)

证明：令 $\widehat{d_k}={\Vert B_k\Vert }_F{d}_k$ ， $\widehat{{\alpha }_k}={\alpha }_k/{\Vert B_k\Vert }_F$ 。则。当 $$ 是有限集时，我们从集合 $K_B$ 的定义(11)以及算法1.1的步4知，对于任何的 $k\ge k_1$ 且有， $u_k={\epsilon }_{k_1}{\Vert B_k\Vert }_F$ 。即满足：

$B_k{d}_k+{\epsilon }_{k_1}{\Vert B_k\Vert }_F{d}_k+\nabla f\left(x_k\right)=0.$

所以

$-\nabla f{\left(x_k\right)}^{\text{T}}\widehat{d_k}=d_k^{\text{T}}{B}_k\widehat{d_k}+{\epsilon }_{k_1}{\Vert \widehat{d_k}\Vert }^2.$ (15)

再由(4)的第二个不等式和 $\nabla f\left(x\right)$ 的Lipschitz连续性得

$L\widehat{{\alpha }_k}{\Vert \widehat{d_k}\Vert }^2=L{\alpha }_k{d}_k^{\text{T}}\widehat{d_k}\ge {\left(\nabla f\left(x_k+{\alpha }_k{d}_k\right)-\nabla f\left(x_k\right)\right)}^{\text{T}}\widehat{d_k}\ge -\left(1-{\sigma }_2\right)\nabla f{\left(x_k\right)}^{\text{T}}\widehat{d_k}.$

由上面不等式以及(15)可推出

$\widehat{{\alpha }_k}\ge \left(1-{\sigma }_2\right){\epsilon }_{k_1}{L}^{-1}>0.$

进而由上面的关系以及(15)和(9)，我们得

$\underset{k\to \infty ,k\in K_B}{\lim }\widehat{d_k}=0.$ (16)

而且对所有的 $k\in K_B$ ，有

$\Vert \nabla f\left(x_k\right)\Vert \le \Vert B_k\Vert \Vert \widehat{d_k}\Vert /{\Vert B_k\Vert }_F+{\epsilon }_{k_1}\Vert \widehat{d_k}\Vert \le \left(1+{\epsilon }_{k_1}\right)\Vert \widehat{d_k}\Vert .$

因而由(16)可得(14)。

引理3.3：如果 $K_B$ 是有限集且 $K_{\tau }$ 也是有限的，那么(14)成立。

证明：我们用反证法假设存在一个某个正数 $\gamma$ 使得对所有k有 $\Vert \nabla f\left(x_k\right)\Vert \ge \gamma$ 成立。由算法的1.1的步4可得，当 $k\notin K_B$ 时， $u_k={\epsilon }_k$ ，并且

${\Vert B_k\Vert }_F<\max \left\{M_B，1/\Vert \nabla f\left(x_k\right)\Vert \right\}\le \max \left\{M_B，1/\gamma \right\}.$ (17)

当 $K_{\tau }$ 是有限集时，则(12)对所有的 $k\ge k_1$ 成立。由引理3.1以及(9)可得 $d_k\to 0\left(k\to \infty \right)$ 。由(3)和(17)可推出，对所有 $k\notin K_B$ 有下面的不等式成立

$\Vert \nabla f\left(x_k\right)\Vert \le \left(\Vert B_k\Vert +{\epsilon }_k\right)\Vert d_k\Vert \le \left({\Vert B_k\Vert }_F+{\epsilon }_1\right)\Vert d_k\Vert \le \left(\max \left\{M_B,1/\gamma \right\}+{\epsilon }_1\right)\Vert d_k\Vert .$

这就意味着

$\underset{k\to \infty ,k\notin K_B}{\lim }\Vert \nabla f\left(x_k\right)\Vert =0.$

这与假设相矛盾，所以必有(14)成立。

在引理3.2和3.3的基础上我们证明了算法1.1的全局收敛性。

定理3.1：设序列 $\left\{x_k\right\}$ 由算法1.1产生并且满足假设3.1则(14)成立。

证明：我们考虑两种情形：(i)是无限集和(ii) $K_{\tau }$ 是有限集。

情形(i)：由 $K_{\tau }$ 的定义(10)，必存在 $\left\{\Vert \nabla f\left(x_k\right)\Vert \right\}$ 中的一个公比不超过 $\eta \left(0<\eta <1\right)$ 几何子序列，此时必有(14)成立。

情形(ii)：我们需要考虑 $K_B$ 是无限还是有限两种情况。由引理3.2和3.3可知，上面两种情况的任意一种情况中，(14)总是成立的，也就是说 ${\epsilon }_k$ 被减少无穷多次，由 $K_{\tau }$ 的定义(10)中可知情形(ii)是不会出现的。

由定理3.1的证明可以看出， ${\epsilon }_k$ 一定被减少无穷多次，所以当 $k\to \infty$ 时， ${\epsilon }_k\to 0$ 。如果 $B_k$ 是有界的，那么当k充分大时， $u_k={\epsilon }_k$ ，并且当 $k\to \infty$ 时， $u_k\to 0$ 。

4. 数值实验与结果分析

在CPU为1.6 HZ、内存为512 MB的电脑上分别对新的扰动的BFGS算法(算法1.1)和文献 [3] 中的CBFGS等算法进行了数值实验，其中测试问题的函数都来源于文献 [4] 。程序采用Matlab语言编写，在算法中我们使用 $$ 作为算法的终止条件，算法中有关常数与参数的选择如下： ${\sigma }_1=0.001$ ， ${\sigma }_2=0.9$ ， $\eta =0.5$ ， $\tau =0.7$ ， ${\epsilon }_1=1.0$ ， $M_B={10}^{10}$ ， $B_1=I$ (单位矩阵)。

首先选择了3个二维问题，1个三维问题，2个四维问题，对算法进行了分析，问题分别是：

rose

$\min f\left(x\right)=100{\left(x_1^2-x_2\right)}^2+{\left(x_1-1\right)}^2,x={\left(x_1,x_2\right)}^{\text{T}}\in R^2,$

初始点： $x_0={\left(-1.2,1\right)}^T$ ，精确解： $x^{*}={\left(1,1\right)}^{\text{T}}$ ， $f\left(x^{\ast }\right)=0$ 。

badscp

$\min f\left(x\right)={\left({10}^4{x}_1{x}_2-1\right)}^2+{\left({\text{e}}^{-x_1}+{\text{e}}^{-x_2}-1.0001\right)}^2,x={\left(x_1,x_2\right)}^{\text{T}}\in R^2$

初始点： $x_0={\left(0,1\right)}^{\text{T}}$ ，精确解： $x^{\ast }={\left(1.098,9.106\right)}^{\text{T}}$ ， $f\left(x^{\ast }\right)=0$ 。

badscb

$\min f\left(x\right)={\left(x_1-{10}^6\right)}^2+{\left(x_2-2\times {10}^{-6}\right)}^2+{\left(x_1{x}_2-2\right)}^2,x={\left(x_1,x_2\right)}^{\text{T}}\in R^2$

初始点： $x_0={\left(1,1\right)}^{\text{T}}$ ，精确解： $x^{\ast }={\left({10}^6,2\times {10}^{-6}\right)}^{\text{T}}$ ， $f\left(x^{\ast }\right)=0$ 。

helix

$\min f\left(x\right)=100{\left[x_3-10\theta \left(x_1,x_2\right)\right]}^2+100{\left[{\left(x_1^2+x_2^2\right)}^{\frac{1}{2}}-1\right]}^2+x_3^2,$

其中 $\theta \left(x_1,x_2\right)=\{\begin{array}{l}\frac{1}{\text{2π}}\arctan \left(\frac{x_2}{x_1}\right),x_1>0\\ \frac{1}{\text{2π}}\arctan \left(\frac{x_2}{x_1}\right)+0.5,x_1<0\end{array}$ ， $x={\left(x_1,x_2,x_3\right)}^{\text{T}}\in R^3$ ，

初始点： ${\left(-1,0,0\right)}^{\text{T}}$ ，精确解： $x^{*}={\left(1,0,0\right)}^{\text{T}}$ ， $f\left(x^{\ast }\right)=0$ 。

sing

$\min f\left(x\right)={\left(x_1+10x_2\right)}^2+5{\left(x_3-x_4\right)}^2+{\left(x_2-2x_3\right)}^4+10{\left(x_1-x_4\right)}^4,x={\left(x_1,x_2,x_3,x_4\right)}^{\text{T}}\in R^4,$

初始点： $x_0={\left(3,-1,0,1\right)}^{\text{T}}$ ，精确解： $x^{*}={\left(0,0,0,0\right)}^{\text{T}}$ ， $f\left(x^{\ast }\right)=0$ 。

wood

$\begin{array}{c}\begin{array}{c}\min f\left(x\right)=100{\left(x_2-x_1^2\right)}^2+{\left(1-x_1\right)}^2+90{\left(x_4-x_3^2\right)}^2+{\left(1-x_3\right)}^2\\ +10{\left(x_2+x_4-2\right)}^2+\frac{1}{10}{\left(x_2-x_4\right)}^2\end{array},x={\left(x_1,x_2,x_3,x_4\right)}^{\text{T}}\in R^4,\\ \end{array}$

初始点： $x_0={\left(-3,-1,-3,-1\right)}^{\text{T}}$ ，精确解 $x^{*}={\left(1,1,1,1\right)}^{\text{T}}$ ， $f\left(x^{\ast }\right)=0$ 。

通过编程计算得到表1。

Table 1. Result analysis of solving several kinds of problems

表1. 求解几类问题的结果分析

通过表1的结果，得出我们的算法能有效地解决二维和三维问题，以及四维问题，说明我们的算法是有效的。

Table 2. Comparisons with other methods

表2. 与其他几类方法的比较分析

通过表2对不同的BFGS方法的迭代步进行分析，我们可以得出我们扰动的算法在求解badscp问题时需要188步，比其他三种方法的迭代步要少，同时在求解其他问题时的迭代步骤相差不大。

5. 小结

通过数值实验结果和全局收敛性的证明，说明我们提出的新扰动的BFGS方法对求解非凸优化问题是有效的。

参考文献