无约束优化问题的两步Broyden秩一校正方法

A Two-Step Broyden Symmetric Rank-One Correction Method for Unconstrained Optimization Problems

摘要: 本文提出了一种求解无约束优化问题的两步Broyden对称秩一校正方法，并证明了其全局收敛性。数值结果表明，两步方法相较于单步方法具有明显的数值效果，同时也说明两步方法具有加速效果。

Abstract: This paper proposes a two-step Broyden symmetric rank-one correction method for solving unconstrained optimization problems and establishes its global convergence. Numerical results demonstrate that the two-step method exhibits significant improvements over the single-step method in terms of computational efficiency, confirming the acceleration effect of the proposed approach.

1. 引言

牛顿法与拟牛顿法是求解无约束优化问题的两类经典方法[1]。传统的牛顿法因其局部二次收敛性而备受关注，但需要计算目标函数的Hessian矩阵及其逆，在高维问题中计算成本较高[2]。拟牛顿法通过构造Hessian矩阵的近似，直接计算二阶导数，从而减少了计算量[3]。

对称秩一(SR1, Symmetric Rank-1)更新是通过秩一校正更新Hessian近似矩阵，具有存储需求低和收敛速度快的优势[4]。Conn等人[5]证明了SR1更新生成的矩阵序列在一定条件下能够全局收敛到真实的Hessian矩阵，尤其适用于非凸优化问题。Khalfan等人[6]的理论与实验研究表明，SR1更新在捕捉Hessian矩阵特征方面具有高效性，然而，SR1方法存在正定性不足的问题，可能影响算法的稳定性。为解决这一问题，Leong和Hassan [7]提出了一种自适应正定缩放的重启SR1方法，有效提升了非凸优化中的性能。此外，Modarres等人[8] [9]通过修正割线方程改进了SR1方法的Hessian近似精度，为本文算法的理论基础提供了支持。

与此同时，两步Broyden方法作为一种加速技术受到广泛关注。Magrenan和Argyros [10]研究了两步Broyden方法的收敛行为，指出额外的修正步骤可以显著提高局部收敛效率。Zhou和Zhang [11]将两步修正思想应用于拟牛顿法，提出了一种改进的Broyden-like方法，并证明了其全局收敛性和超线性收敛性。此外，Ali等人[12]提出了一种多步SR1更新方法，通过插值多项式改进Hessian近似，为两步策略与SR1的结合提供了理论支持。

尽管SR1和两步Broyden方法各有优势，但二者的结合尚未得到充分探索。传统SR1方法因正定性不足可能导致收敛不稳定，而两步Broyden方法通过额外修正可提高局部收敛速度。基于此，本文提出了一种结合SR1更新和两步Broyden方法的算法，利用SR1生成高效Hessian近似，并通过两步修正克服正定性缺陷，从而兼顾计算效率和稳定性。

与MSSR1方法[12]相比，TSSR1-B在以下方面有所不同：首先是MSSR1采用插值多项式构造多步割线方程，而TSSR1-B通过两步Broyden修正构造更简洁的割线方程，避免了高次插值带来的复杂性；其次是TSSR1-B引入自适应缩放因子和重启机制，更灵活地保持矩阵正定性和数值稳定性。

本文的贡献在于：在SR1理论框架中引入两步修正步骤，提升了算法的收敛速度；通过理论分析证明了算法的全局收敛性；数值实验验证了TSSR1-B在迭代次数和计算时间上的优越性。

本文结构如下：第二节介绍无约束优化问题的数学模型及算法框架，第三节详细描述修正割线方程的构造，第四节提出结合SR1更新和两步修正的TSSR1-B算法，第五节分析算法的全局收敛性，第六节通过数值实验对比TSSR1-B与传统方法的性能，第七节总结全文并展望未来研究方向。

2. 无约束优化问题的数学模型及算法框架

本节提出一种无约束优化问题的两步Broyden秩一校正方法(TSSR1-B)，用于求解如下形式的问题：

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

其中$f:R^n\to R$ 是一个光滑函数。假设$f$ 是连续的，并且至少是二次可微的。TSSR1-B方法的基本框架可以概括如下：

给定$x_k,\nabla f\left(x_k\right)$ 或其有限差分近似，以及$B_k\in R^{n\times n}$ (${\nabla }^{2}f\left(x_k\right)$ 的割线近似)，我们通过下式确定TSSR1-B方向$d_k$

$B_k{d}_k+\nabla f\left(x_k\right)=0,$ (1.2)

通过基于以下等式的线搜索方法选择新的迭代$x_{k+1}$

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

其中

$d_k={\widehat{d}}_k+{\overline{d}}_k,{\widehat{d}}_k=-H_k\nabla f\left(x_k\right)，{\overline{d}}_k=-H_k\nabla f\left(x_k+{\widehat{d}}_k\right),$ (1.4)

将$B_k$ 更新为$B_{k+1}$ ，使$B_{k+1}$ 对称并满足正割方程

$B_{k+1}{s}_k=y_k,$ (1.5)

其中$s_k=x_{k+1}-x_k,y_k=f\left(x_{k+1}\right)-f\left(x_k\right)$ 。

2.1. ${\tilde{y}}_k$ 的构造

将目标函数$f$ 的良好曲率信息合并到更新矩阵的方法之一是强制矩阵满足一般正割方程(1.5)。为了生成新的正割更新，我们可以从TSSR1-B方向(1.2)生成方向$d_k$ ，其中$B_k$ 被矩阵

${\tilde{B}}_k=B_k+{\mu }_{k}I,$ (2.1)

其中$I$ 是单位矩阵，为更好地近似Hessian矩阵的曲率信息，引进修正项$u_k$ ，$u_k$ 以适当的方式选择($u_k$ 的实际选择将在下一小节中讨论)。由于$B_{k+1}$ 满足

$B_{k+1}\left(x_{k+1}-x_k\right)\approx \nabla f\left(x_{k+1}\right)-\nabla f\left(x_k\right),$

矩阵$B_{k+1}$ 满足以下关系：

${\tilde{B}}_{k+1}\left(x_{k+1}-x_k\right)=\left(B_{k+1}+{\mu }_{k}I\right)\left(x_{k+1}-x_k\right)\approx \nabla f\left(x_{k+1}\right)-\nabla f\left(x_k\right),$

或等效地

$B_{k+1}{s}_k\approx y_k+{\mu }_k{s}_k,$

所以我们可以把修改后的正割方程写成如下

$B_{k+1}{s}_k=y_k^{\ast }\text{;}{y}_k^{\ast }=y_k+{\mu }_k{s}_k.$ (2.2)

很明显，如果$k$ 趋近于无穷大，$u_k$ 趋近于零，那么修正的正割方程就接近于一般的正割方程。

现在假设目标函数$f$ 足够光滑。设$x_{k+1}=x_k+s_k,$ 通过沿着$s_k$ 的泰勒级数，我们有

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

或者

$f\left(x_k\right)\cong f\left(x_{k+1}\right)-s_k^T\nabla f\left(x_{k+1}\right)+\frac{1}{2}{s}_k^T{\nabla }^{2}f\left(x_{k+1}\right)s_k,$

因此，

$s_k^T{\nabla }^{2}f\left(x_{k+1}\right)s_k=2\left(f\left(x_k\right)-f\left(x_{k+1}\right)\right)+2{\left(\nabla f\left(x_{k+1}\right)\right)}^T{s}_k$

$=2\left(f\left(x_k\right)-f\left(x_{k+1}\right)\right)+{\left(\nabla f\left(x_{k+1}\right)+\nabla f\left(x_k\right)\right)}^T{s}_k+s_k^T{y}_k,$ (2.3)

通过(2.2)，我们有

$s_k^T{B}_{k+1}{s}_k=s_k^T{y}_k^{\ast }=s_k^T{y}_k+{\mu }_k{s}_k^T{s}_k,$ (2.4)

通过比较(2.3)与(2.4)，我们可以以合理的方式选择$u_k$ ，使其满足

${\mu }_k{s}_k^T{s}_k={\psi }_k,$ (2.5)

其中${\psi }_k=2\left(f\left(x_k\right)-f\left(x_{k+1}\right)\right)+{\left(\nabla f\left(x_{k+1}\right)+\nabla f\left(x_k\right)\right)}^T{s}_k.$

根据公式(2.5)，我们可以很好的选择$u_k$

${\mu }_k{s}_k=w_k;w_k=\frac{{\psi }_k}{s_k^T{u}_k}{u}_k,$ (2.6)

且$u_k\in R^n$ 是使得$s_k^{{\rm T}}{u}_k\ne 0$ 。

Modarres等人基于上述关系提出了以下等式

$B_{k+1}{s}_k=y_k^{\ast };y_k^{\ast }=y_k+\frac{{\psi }_k}{s_k^T{u}_k}{u}_k.$ (2.7)

本文所考虑的是通常的$l_2$ -范数向量和相应的诱导范数矩阵。

下面的定理显示了割线方程(2.7)的好的性质，可以在[8] [9]中找到。

定理1. [8] 假设函数$f$ 足够光滑，且$u_k$ 满足式(2.5)。如果$\Vert s_k\Vert$ 足够小，则以下估计保持

$s_k^T\left({\nabla }^{2}f\left(x_{k+1}\right)s_k-y_k^{\ast }\right)=\frac{1}{3}{s}_k^T\left(T_{k+1}{s}_k\right)s_k+O\left({\Vert s_k\Vert }^4\right).$ (2.8)

$s_k^T\left({\nabla }^{2}f\left(x_{k+1}\right)s_k-y_k\right)=\frac{1}{2}{s}_k^T\left(T_{k+1}{s}_k\right)s_k+O\left({\Vert s_k\Vert }^4\right).$ (2.9)

其中$y_k^{\ast }=y_k+u_k{s}_k.$ $T_{k+1}$ 是$f$ 在$x_{k+1}$ 处的张量，并且

$s_k^T\left(T_{k+1}{s}_k\right)s_k={\displaystyle \sum _{i,j,l=1}^n\frac{{\partial }^{3}f\left(x_{k+1}\right)}{\partial x^i\partial x^j\partial x^l}}{s}_k^i{s}_k^j{s}_k^l.$

上面的定理表明，$y_k^{\ast }$ 比$y_k$ 更好地近似于${\nabla }^{2}f\left(x_{k+1}\right)s_k$ 。

2.2. 对正割方程(2.7)的一个简单修改

由于$s_k^{{\rm T}}{y}_k^{\ast }=s_k^{{\rm T}}{y}_k+u_k，$ 且$u_k$ 可能为负，如果$u_k<0$ 大于$s_k^{{\rm T}}{y}_k$ ，并且$s_k^{{\rm T}}{y}_k^{\ast }<0.$ 。为了解决这个缺陷，我们进一步修改向量$y_k^{\ast }$ ，如下所示

${\tilde{y}}_k=y_k+\text{sgn}\left({\psi }_k\right)\frac{{\psi }_k}{s_k^T{u}_k}{u}_k,$ (2.10)

因此很容易看出，我们总是有$s_k^{{\rm T}}{y}_k^{\ast }>0，$ 当且仅当$s_k^{{\rm T}}{y}_k>0$ ，

公式(2.10)中代入$u_k=s_k$ ，我们得到

${\tilde{y}}_k=y_k+\text{sgn}\left({\psi }_k\right)\frac{{\psi }_k}{s_k^T{u}_k}{s}_k.$ (2.11)

3. 正割方程的SR1修正

在这一节中，我们考虑用修正的割线方程(2.10)对SR1进行修正。此外，我们将介绍一种策略来保持修正的SR1公式中的正定性。

3.1. 修改的SR1更新

现在我们提出基于修改的正割方程(2.10)的修改的SR1更新：

$B_{k+1}=B_k+\frac{\left({\tilde{y}}_k-B_k{s}_k\right){\left({\tilde{y}}_k-B_k{s}_k\right)}^T}{{\left({\tilde{y}}_k-B_k{s}_k\right)}^T{s}_k},$ (3.1)

其中${\tilde{y}}_k=y_k+\text{sgn}\left({\psi }_k\right)\frac{{\psi }_k}{s_k^T{u}_k}{u}_k$。如果用$H_k$ 近似Hessian的逆，则修改的割线方程被重写为

$s_k=H_{k+1}{\tilde{y}}_k,$ (3.2)

基于该等式，我们具有逆SR1更新

$H_{k+1}=H_k+\frac{\left(s_k-H_k{\tilde{y}}_k\right){\left(s_k-H_k{\tilde{y}}_k\right)}^T}{{\left(s_k-H_k{\tilde{y}}_k\right)}^T{\tilde{y}}_k},$ (3.3)

我们将更新(3.1)称为修改的SR1更新，而(3.3)是其逆版本。

3.2. 修改的SR1更新的正定性

在本小节中，我们将对对称秩一更新应用重新启动过程，以避免正定性的丢失。重新开始方法和常用的调整方法之间的重要区别在于，调整被合并到每个SR1更新中以保持正定性，而重新开始方法在SR1更新不是正定时用单位矩阵的倍数重新开始SR1算法。我们不喜欢这种方法的主要原因是，它需要计算$B_k$ (Hessian的近似)。如果该迭代中的SR1矩阵是正定的，则大小调整将是多余的。

为了克服这个缺陷，我们采用Leong和Hassan [6]提出的最佳缩放因子，并考虑由${\tilde{\lambda }}_{k-1}I$ 定义的缩放单位矩阵，其中

${\tilde{\lambda }}_{k-1}=\frac{s_{k-1}^2{s}_{k-1}}{{\tilde{y}}_{k-1}^2{s}_{k-1}}-{\left(\frac{{\left(s_{k-1}^2{s}_{k-1}\right)}^2}{{\left({\tilde{y}}_{k-1}^2{s}_{k-1}\right)}^2}-\frac{s_{k-1}^2{s}_{k-1}}{{\tilde{y}}_{k-1}^2{\tilde{y}}_{k-1}}\right)}^{1/2}.$ (3.4)

我们的方法是结合重新启动程序，当我们失去了正定性的修改SR1或修改SR1的分母接近零或$\left\{H_k\right\}$ 是无界的。现在，通过这种修改，我们提出了我们的新算法。

4. 算法描述

算法4.1：无约束优化问题的两步Broyden秩一校正方法(TSSR1-B)。

Step 0：给定初始点$x_0\in R^n,$ 初始正定矩阵$H_0=I$ ，令$k=0$ 。

Step 1：若$\Vert \nabla f\left(x_k\right)\Vert \le \epsilon \cdot \max \left(1,\Vert x_k\Vert \right)$ ，则停止迭代。

Step 2：通过如下方程式来求解${\widehat{d}}_k$

${\widehat{d}}_k=-H_k\nabla f\left(x_k\right),$ (4.1)

其中令：

$z_k=x_k+{\widehat{d}}_k,$ (4.2)

通过下列方程式求解${\overline{d}}_k$

${\overline{d}}_k=-H_k\nabla f\left(z_k\right),$ (4.3)

令：

$d_k={\widehat{d}}_k+{\overline{d}}_k,$ (4.4)

该方向结合了当前梯度信息和历史曲率信息，能更有效地引导搜索方向。

Step 3：若以下任一条件成立：

1) $s_k^T{\tilde{y}}_k-{\tilde{y}}_k^T{H}_k{\tilde{y}}_k<0$ (搜索方向非下降)， (4.5)

2) $\left|{\tilde{y}}_k^T\left(s_k-H_k{\tilde{y}}_k\right)\right|<r\Vert {\tilde{y}}_k\Vert \Vert s_k-H_k{\tilde{y}}_k\Vert$ (SR1更新分母接近零)， (4.6)

3) ${\Vert H_k\Vert }_{\infty }>L$ (近似矩阵范数过大)， (4.7)

则重启：$H_k={\tilde{\lambda }}_{k-1}I,$ 其中${\tilde{\lambda }}_{k-1}$ 由式(3.4)计算。

Step 4：尝试${\alpha }_k=1,$ 若满足Wolfe条件：

$f\left(x_k+{\alpha }_k{d}_k\right)\le f(x_k)+{\delta }_1{\alpha }_k\nabla f{\left(x_k\right)}^T{d}_k,\nabla f{\left(x_k+{\alpha }_k{d}_k\right)}^T{d}_k\ge {\delta }_2\nabla f{\left(x_k\right)}^T{d}_k,$ (4.8)

则接受；否则通过回溯法(如${\alpha }_k\leftarrow r{\alpha }_k$ )找到满足条件的步长。

Step 5：更新迭代点：

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

Step 6：通过修正割线方程计算${\tilde{y}}_k$

${\tilde{y}}_k=y_k+\text{sgn}\left({\psi }_k\right)\frac{{\psi }_k}{s_k^T{u}_k}{u}_k,$ (4.9)

其中

$w_k=2\left(f\left(x_k\right)-f\left(x_{k+1}\right)\right)+{\left(\nabla f\left(x_{k+1}\right)+\nabla f\left(x_k\right)\right)}^T{s}_k,$ (4.10)

Step 7：通过(3.3)使用对称秩一更新逆Hessian近似。

Step 8：令$k=k+1$ ，返回Step 1。

5. 全局收敛性

为了证明算法的全局收敛性，我们需要如下假设：

(1) $f$ 在水平集$D=\left\{x\in R^n:f\left(x\right)\le f\left(x_0\right)\right\}$ 的领域内连续可微，且有下界。

(2) 梯度$\nabla f$ 在$D$ 上Lipschitz连续，即存在常数$L>0,$ 使得：

$$ $\Vert \nabla f\left(x_2\right)-\nabla f\left(x_1\right)\Vert \le L\Vert x_2-x_1\Vert ,\forall x_1,x_2\in D.$ (5.1)

定义：

$\cos \left({\theta }_k\right)=-\frac{\nabla f{\left(x_k\right)}^T{d}_k}{\Vert \nabla f\left(x_k\right)\Vert \Vert d_k\Vert }.$ (5.2)

下面的定理表明算法4.1具有全局收敛性。

定理1：若$x_0$ 是满足上述假设的起点，算法生成的序列$\left\{x_k\right\}$ 满足强Wolfe条件，则：

${\displaystyle \sum _{k=0}^{\infty }{\cos }^2}\left({\theta }_k\right){\Vert \nabla f\left(x_k\right)\Vert }^2<\infty .$ (5.3)

证明：由强Wolfe条件：

${\left(\nabla f\left(x_{k+1}\right)-\nabla f\left(x_k\right)\right)}^T{d}_k\ge \left({\sigma }_2-1\right)\nabla f{\left(x_k\right)}^T{d}_k,$

结合Lipschitz条件：

${\left(\nabla f\left(x_{k+1}\right)-\nabla f\left(x_k\right)\right)}^T{d}_k\le {\alpha }_{k}L{\Vert d_k\Vert }^2,$

得：

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

由第一个Wolfe条件：

$f\left(x_{k+1}\right)\le f\left(x_k\right)+{\sigma }_1{\alpha }_k\nabla f{\left(x_k\right)}^{{\rm T}}{d}_k,$

代入${\alpha }_{k,}$ 得：

$f\left(x_{k+1}\right)\le f\left(x_k\right)+{\sigma }_1\frac{{\sigma }_2-1}{L}\frac{{\left(\nabla f{\left(x_k\right)}^{{\rm T}}{d}_k\right)}^2}{{\Vert d_k\Vert }^2}，$

使用$\cos ({\theta }_k),$ 可写为：

$f\left(x_{k+1}\right)\le f\left(x_k\right)+c{\cos }^2\left({\theta }_k\right){\Vert \nabla f\left(x_k\right)\Vert }^2,c={\sigma }_1\frac{{\sigma }_2-1}{L}，$

由于$f$ 有下界，对上式求和：

${\displaystyle \sum _{k=0}^{\infty }{\cos }^2}\left({\theta }_k\right){\Vert \nabla f\left(x_k\right)\Vert }^2<\infty ，$

由于$H_k$ 有界(由跳跃更新和正定性修正保证)，存在常数$\Delta >0$ ，使得：

$\cos \left({\theta }_k\right)\ge \frac{1}{\Delta }，$

因此：

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

即算法全局收敛。

6. 数值实验

本小节给出一些数值结果，以表明所提出的无约束优化问题的两步Broyden秩一校正方法(TSSR1-B)的有效性，并与[9]中的MSR1方法、传统的SR1方法进行比较，算法在MATLAB2016A编程实现，数值实验在windows系统中进行。当$\Vert \nabla f\left(x_k\right)\Vert <{10}^{-8}$ 或达到最大迭代次数500时，终止实验。在同样的参数设置下(具体的参数设置为：线搜索参数${\delta }_1\text{=}{10}^{-4},{\delta }_2=0.9$ ，重启参数$r={10}^{-12},L={10}^{10}$ )。下表1展示不同算法在四种测试函数上的平均迭代次数和CPU时间(单位：秒)：

Table 1. Performance comparison of different algorithms on test functions

表1. 不同算法在测试函数上的性能对比

函数名称	维度	TSSR1B (迭代/时间)	MSR1 (1) (迭代/时间)	SR1 (迭代/时间)
强病态二次函数	30	12/0.007	26/0.032	35/0.012
	50	8/0.006	26/0.041	53/0.022
	80	10/0.008	27/0.055	89/0.042
	100	18/0.010	29/0.073	106/0.075
标准Rastrigin函数	30	46/0.024	72/0.020	98/0.270
	50	83/0.033	112/0.198	77/0.016
	80	104/0.062	500/0.230	104/0. 029
	100	24/0.028	500/0.330	113/0.040
非光滑过渡函数	30	4/0.002	8/0.005	10/0.012
	50	4/0.002	10/0.002	11/0.008
	80	6/0.003	9/0.003	10/0.005
稀疏病态函数	30	10/0.010	40/0.047	178/0.088
	50	16/0.011	42/0.110	213/0.021
	80	23/0.016	45/0.154	500/0.317
	100	8/0.020	42/0.205	500/0.476

表格列举出了算法TSSR1-B、MSR1和SR1在四种测试函数上的不同初始点的数值结果。由表格可知，虽然三种算法都能求解测试问题，但是算法TSSR1-B的数值优于MSR1和SR1。在相同的条件下，算法TSSR1-B在大多数问题上进一步减少15%~40%的迭代次数。两步方向在病态问题上贡献更为显著，但在极低维简单问题上，两步优势不明显，单步方法可能更高效。

7. 结论

本文提出的无约束优化问题的两步Broyden秩一校正方法(TSSR1-B)通过结合SR1更新和两步修正策略，提升了无约束优化问题的求解效率。理论分析和数值实验验证了其优越性。未来研究可探索该方法在随机优化和大规模分布式场景中的应用。