一种带记忆的阻尼牛顿法

A Damping Newton Method with Memory

摘要: 本文在引进阻尼技术的基础上构造了新的拟牛顿方程，该方程通过利用相邻两次迭代点之间的梯度差和变化量的线性组合，实现了曲率对的更新，以此改进校正公式。为求解高维问题，引入有限内存记忆框架，提出了D-LBFGS算法。利用D-LBFGS算法对神经网络中的权值阈值参数进行优化的算法，进一步提出了D-LBFGS—CNN算法，在解决无创光谱血糖浓度的预测问题中迭代351次时精度达到99.999%。所提出的D-LBFGS算法在一定程度上优于标准的LBFGS算法(迭代474次时精度99.999%)、阻尼牛顿法(精度99.359%)及文献中的LBFGS算法(精度99.744%)。

Abstract: This paper constructs a new quasi-Newton equation based on the introduction of damping technology. This equation uses the linear combination of the gradient difference and change between two adjacent iteration points to update the curvature pair, thereby improving the correction formula. In order to solve high-dimensional problems, the limited memoryframe work is introduced and the D-LBFGS algorithm is proposed. The D-LBFGS algorithm is used to optimize the weight threshold parameters in the neural network, and the D-LBFGS—CNN algorithm is further proposed. The accuracy reaches 99.999% when iterating 351 times to solve the prediction problem of non-invasive spectral blood glucose concentration. The proposed D-LBFGS algorithm is superior to the standard LBFGS algorithm (accuracy 99.999% at 474 iterations), the damped Newton method (accuracy 99.359%), and the LBFGS algorithm in the literature (accuracy 99.744%) to a certain extent.

1. 引言

考虑无约束优化问题

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

其中，$f:R^n\to R$ 为连续可微函数。牛顿法的优点是收敛速度快，局部收敛性强，适用于求解方程组或函数最小值问题。但牛顿法需要计算二阶导数，计算复杂性高，而且目标函数的Hesse矩阵可能出现非正定的情况。

为了克服牛顿法中出现的问题，1959年由Davidon提出了第一个拟牛顿方法，是在牛顿法的基础上减少对目标函数的二阶导数信息的需求[1]。相对于传统的牛顿法需要存储目标函数的精确Hesse矩阵的逆，拟牛顿法使用一个近似的Hesse矩阵来替代真正的Hesse矩阵，降低了计算复杂性和存储成本。它通过观察目标函数的一阶导数的变化来构造这个近似矩阵，因此拟牛顿法中，只需要目标函数的一阶信息。这种存储成本的降低在处理大规模数据或高维问题时尤为重要，不同构造近似矩阵的方法使得出现不同的拟牛顿法，使得拟牛顿法成为一种高效的优化算法选择。BFGS算法是在1970年由Broyden，Fletcher，Goldfard和Shanno所提出的，通过数值实验表明，该算法表现效果更好，因此得到广泛应用。1980年，Nocedal提出了有限记忆BFGS算法(简记为LBFGS)，该方法没有直接储存n阶Hesse矩阵，而是储存m对曲率对[2]。1989年Liu等在一致凸问题上证明了L-BFGS算法的全局收敛性，同时实验结果表明，LBFGS算法能够求解大规模无约束问题，更进一步的降低了存储和计算成本，并且在计算上具有更高的效率[3]。

2017年Wang等人提出了使用经典LBFGS算法框架的非凸随机优化的SdLBGFS算法[4]；2018年，Ioannis E. Livieris等人提出了一种新的混合共轭梯度方法，总体上优于经典的共轭梯度方法。基于自缩放无记忆BFGS方法的强大的理论特性和良好的计算性能，通常被认为是可以解决大规模优化问题的最有效方法之一[5]；2019年，陈昱含等人提出了自调比正则化LBFGS (RLSBFGS)方法，该方法通过增加少量的存储来改善了LBFGS方法在病态问题上的计算效率。数值结果表明，RLSBFGS方法是略优于LBFGS [6]；2020年，Albert S. Berahas等人，提出了采样LBFGS方法，在当前迭代周围随机采样点，所构建的近似值利用了更可靠的信息，没有依赖于可能明显过时的过去信息[7]；2020年，Daniela di Serafino等人对最速下降(SD)方法重新产生了兴趣，将最速下降(SD)策略应用于牛顿法和BFGS方法[8]；2022年，Neculai Andrei等人提出了无内存BFGS方法，该方法在每次迭代时通过使用单位矩阵进行初始化[9]。

下面提出了一种改进的LBFGS算法，通过改进曲率对的更新方式，更多的利用相邻迭代点的数据信息，提出新的拟牛顿方程。在此基础上提出了一种利用D-LBFGS算法对神经网络中的权值阈值参数进行优化的算法，记为D-LBFGS—CNN算法，在反向更新过程中使用D-LBFGS算法更新权值阈值参数，优化神经网络的性能，解决无创光谱血糖浓度的预测问题，提升了无创血糖浓度预测问题的精度。

2. 拟牛顿法

拟牛顿法大多在定义拟牛顿方程、计算搜索方向以及参数更新方面有所不同。常见的拟牛顿法有SR1算法，DFP算法，BFGS算法等，其中BFGS算法是求解无约束优化问题的最有效的拟牛顿法之一。

2.1. BFGS算法

考虑无约束优化问题(1)，BFGS算法有如下的迭代格式

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

其中，$x_k$ 是第k次的迭代点，${\alpha }_k$ 是步长，$d_k$ 是搜索方向。

令$g_k$ 表示$f\left(x\right)$ 在$x_k$ 处的梯度，$d_k$ 是通过求解$B_k{d}_k=-g_k$ 得到的搜索方向，其中$B_k$ 是$f\left(x\right)$ 在$x_k$ 处的Hesse矩阵的近似，满足以下拟牛顿方程：

$B_{k+1}{s}_k=y_k.$ (3)

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

在BFGS算法中，$B_k$ 更新公式为

$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{y_k{y}_k^{\text{T}}}{y_k^{\text{T}}{s}_k}.$ (4)

运用两次Sherman-Morrison-Woodbury公式，可以推出Hesse矩阵的逆矩阵的近似迭代公式：

$H_{k+1}=\left(I-\frac{s_k{y}_k^{\text{T}}}{y_k^{\text{T}}{s}_k}\right)H_k\left(I-\frac{s_k{y}_k^{\text{T}}}{y_k^{\text{T}}{s}_k}\right)+\frac{s_k{s}_k^{\text{T}}}{y_k^{\text{T}}{s}_k}.$ (5)

2.2. 有限记忆的拟牛顿算法

为更多地减少矩阵的存储量，降低计算成本，保持更快的收敛速率，Liu和Nocedal等人，提出了有限记忆的拟牛顿法(LBFGS算法)，只需要存储最近的(已知给定的有限记忆数) m对曲率对$\left\{\left(s_i,y_i\right)\right\}$ ，而不需要显式地存储完整的Hesse矩阵或其逆[3]。这使得该算法非常适用于大规模优化问题。

BFGS算法中$H_{k+1}$ 的校正公式为：

$H_{k+1}=V_k^{\text{T}}{H}_k{V}_k+{\rho }_k{s}_k{s}_k^{\text{T}}.$

其中，$V_k=I-{\rho }_k{y}_k{s}_k^{\text{T}}$ ，${\rho }_k=\frac{1}{s_k^{\text{T}}{y}_k}$ 。当$k\ge m$ 时，为第k次迭代选择一个初始矩阵$H_k^{\left(0\right)}$ ，代入上式重复迭代m次，且只保留最近m步，得到LBFGS算法的$H_{k+1}$ 校正公式：

$H_k^{\left(j+1\right)}=V_{k-m+j}^{\text{T}}{H}_k^{\left(j\right)}{V}_{k-m+j}+{\rho }_{k-m+j}{s}_{k-m+j}{s}_{k-m+j}^{\text{T}},j=0,1,\cdots ,m-1.$

令$H_{k+1}=H_k^{\left(m\right)}$ 。

Jorge Nocedal等人提出了使用双循环递归算法实现LBFGS算法，当$k\ge m$ 时，利用双循环递归算法计算得到r即$H_k{g}_k$ ，由$d_k=-r$ 得到搜索方向，步骤如下：

算法1：双循环递归(two-loop recursion)算法

初始化：$q=g_k$ ，有限记忆数m，$H_k^{\left(0\right)}$ ，曲率对${\left\{\left(s_i,y_i\right)\right\}}_{i=k-m}^{k-1}$ ；

for $i=k-1,k-2,\cdots ,k-m$ ；

计算${\rho }_i=\frac{1}{y_i^{\text{T}}{s}_i}$ ；

计算并保存${\alpha }_i={\rho }_i{\left(s_i\right)}^{\text{T}}q$ ；

更新$q=q-{\alpha }_i{y}_i$ ；

end for；

初始化$r=H_k^{\left(0\right)}q$ ；

for $i=k-m,\cdots ,k-1$ ；

计算$\beta ={\rho }_i{\left(y_i\right)}^{\text{T}}r$ ；

更新$r=r+\left({\alpha }_i-\beta \right)s_i$ ；

end for；

输出：r即$\left(H_k{g}_k\right)$ 。

在实际应用中，可取$H_k^{\left(0\right)}={\xi }_{k}I$ ，其中，I表示单位矩阵，${\xi }_k=\frac{s_{k-1}^{\text{T}}{y}_{k-1}}{y_{k-1}^{\text{T}}{y}_{k-1}}$ 。

3. 阻尼牛顿法

1978年，Powell定义了第一个阻尼牛顿校正公式，对BFGS算法进行修正[10]。2006年，Nocedal在第18.3节中介绍了阻尼牛顿技术[11]，令$y_k$ 与$B_k{s}_k$ 的组合形式记作${\gamma }_k$ ：

${\gamma }_k={\theta }_k{y}_k+\left(1-{\theta }_k\right)B_k{s}_k.$ (6)

其中阻尼参数${\theta }_k$ 的取值如下：

${\theta }_k=\{\begin{array}{l}\frac{0.8s_k^{\text{T}}{B}_k{s}_k}{s_k^{\text{T}}{B}_k{s}_k-s_k^{\text{T}}{y}_k},若s_k^{\text{T}}{y}_k<0.2s_k^{\text{T}}{B}_k{s}_k,\\ 1,其他,\end{array}$ (7)

可知${\theta }_k\in \left(0,1\right]$ 。用${\gamma }_k$ 代替$y_k$ ，则带有阻尼技术的BFGS方法的Hesse矩阵近似的校正公式为：

$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{{\gamma }_k{\gamma }_k^{\text{T}}}{s_k^{\text{T}}{\gamma }_k}.$ (8)

2021年，赵睿颖提出一种改进的 LBFGS算法，利用相邻迭代点梯度之差，以及由此定义的矩阵与向量之间的乘积的线性组合来更新向量对，构造新的Hesse矩阵近似[12]。算法将迭代方向进行单位化，以保证所提算法的稳定性，给出了一个新的拟牛顿方程：

$H_{k+1}{s}_k={\widehat{y}}_k={\theta }_k{y}_k+\left(1-{\theta }_k\right)M_k{s}_k.$ (9)

其中，$M_k$ 为如下定义的对角矩阵

$M_k=\left(\begin{array}{ccc}\sqrt{{\left(y_{k-1}\right)}_1^2}& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & \sqrt{{\left(y_{k-1}\right)}_n^2}\end{array}\right)\text{,}$ (10)

其中，${\left(y_{k-1}\right)}_i$ 表示向量$y_{k-1}$ 的第i个分量，

${\theta }_k=\{\begin{array}{l}0.25s_k^{\text{T}}{M}_k{s}_k,若s_k^{\text{T}}{y}_k<0.25s_k^{\text{T}}{s}_k,\\ s_k^{\text{T}}{y}_k,其他.\end{array}$ (11)

当沿着迭代方向$d_k$ 进行更新时，迭代方向$d_k$ 进行单位化仍记为$d_k=\frac{d_k}{\Vert d_k\Vert }$ ，使改进的算法更加稳定。

4. D-LBFGS算法提出

综合考虑阻尼牛顿技术和有限记忆的思想，提出一种带记忆的阻尼牛顿法记为“D-LBFGS算法”，是一种基于L-BFGS算法框架的变种，其中引入了对角矩阵。通过利用相邻两次迭代点之间的梯度差和两次迭代点之间的变化量之间乘积的线性组合，实现曲率对的更新。这一改进旨在更有效地使用当前信息适应问题的特性，提高算法的收敛速度和性能。构建的拟牛顿方程如下：

$H_{k+1}{s}_k={\widehat{y}}_k={\theta }_k{y}_k+\left(1-{\theta }_k\right)M_k{s}_k.$ (12)

其中$M_k$ 定义为下面的对角矩阵，${\left(y_k\right)}_i$ 表示向量$y_k$ 的第i个分量，

$M_k=\left(\begin{array}{ccc}\sqrt{{\left(y_k\right)}_1^2}& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & \sqrt{{\left(y_k\right)}_n^2}\end{array}\right)$ (13)

${\theta }_k$ 的取值如下：

${\theta }_k=\{\begin{array}{l}\frac{0.8s_k^{\text{T}}{B}_k{s}_k}{s_k^{\text{T}}{B}_k{s}_k-s_k^{\text{T}}{y}_k},若s_k^{\text{T}}{y}_k<0.2s_k^{\text{T}}{B}_k{s}_k,\\ 1,其他,\end{array}$ (14)

Hesse矩阵的近似矩阵校正公式为：

$B_{k+1}=B_k+\frac{{\widehat{y}}_k{\widehat{y}}_k^{\text{T}}}{{\widehat{y}}_k^{\text{T}}{s}_k}-\frac{B_k{s}_k{s}_k^{\text{T}}{B}_k}{s_k^{\text{T}}{B}_k{s}_k}.$ (15)

D-LBFGS算法的Hesse矩阵的逆的近似矩阵校正公式为：

$H_{k+1}=\left(I-\frac{s_k{\widehat{y}}_k^{\text{T}}}{{\widehat{y}}_k^{\text{T}}{s}_k}\right)H_k\left(I-\frac{{\widehat{y}}_k{s}_k^{\text{T}}}{{\widehat{y}}_k^{\text{T}}{s}_k}\right)+\frac{s_k{s}_k^{\text{T}}}{{\widehat{y}}_k^{\text{T}}{s}_k}.$ (16)

令${\widehat{\rho }}_k=\frac{1}{s_k^{\text{T}}{\widehat{y}}_k}$，再令$V_k=I-{\widehat{\rho }}_k{\widehat{y}}_k{s}_k^{\text{T}}$，则(16)式可化为：

$H_{k+1}=V_k^{\text{T}}{H}_k{V}_k+{\widehat{\rho }}_k{s}_k{s}_k^{\text{T}}.$

考虑到LBFGS算法能减少矩阵的存储量，降低计算成本的思想，给定有限记忆次数非负整数m，给定初始$H_k^{\left(0\right)}$ ，使用上述公式m次，令$H_{k+1}=H_k^{\left(m\right)}$ 得到D-LBFGS算法的Hesse矩阵的逆的近似矩阵校正公式为：

$\begin{array}{c}{H}_{k+1}=H_k^{\left(m\right)}\\ =V_k^{\text{T}}{V}_{k-1}^{\text{T}}\cdots V_{k-m+1}^{\text{T}}{H}_k^{\left(0\right)}{V}_{k-m+1}\cdots V_{k-1}{V}_k\\ +V_k^{\text{T}}{V}_{k-1}^{\text{T}}\cdots V_{k-m+2}^{\text{T}}{\widehat{\rho }}_{k-m+1}{s}_{k-m+1}{s}_{k-m+1}^{\text{T}}{V}_{k-m+2}\cdots V_{k-1}{V}_k\\ ⋮\\ +V_k^{\text{T}}{\widehat{\rho }}_{k-1}{s}_{k-1}{s}_{k-1}^{\text{T}}{V}_k\\ +{\widehat{\rho }}_k{s}_k{s}_k^{\text{T}}\end{array}$ (17)

D-LBFGS算法通过weak Wolfe-Powell (WWP)非精确线搜索得到步长${\alpha }_k$ ，步长规则：

$f\left(x_k+{\alpha }_k{d}_k\right)\le f\left(x_k\right)+\delta {\alpha }_k{g}_k^{\text{T}}{d}_k.$ (18)

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

其中，$\delta \in \left(\text{0,}\frac{\text{1}}{\text{2}}\right)$ ，$\sigma \in \left(\delta ,\text{1}\right)$ 。

D-LBFGS算法的计算步骤如下：

算法2：D-LBFGS算法

Step 1：给定初始点$x_0\in R^n$ ，$H_0=B_0=I$ ，有限记忆数$m>0$ ，运行误差$\epsilon$ ，$\delta$ ，$\sigma$ ，$k=0$ ，令$g_k=\nabla f\left(x_k\right)$ ；

Step 2：如果$\Vert g_k\Vert \le \epsilon$ ，则停止迭代，记$x^{*}=x_k$ ；否则转入Step 3；

Step 3：当$k<m$ 时，计算更新方向$d_k=-H_k{g}_k$ ；当$k\ge m$ 时，计算${\widehat{\xi }}_k=\frac{s_{k-1}^{\text{T}}{\widehat{y}}_{k-1}}{{\widehat{y}}_{k-1}^{\text{T}}{\widehat{y}}_{k-1}}$，得到$H_k^{\left(0\right)}={\widehat{\xi }}_{k}I$ ，通过算法1输出的r即$H_k{g}_k$ ，计算$d_k=-r$ ；

Step 4：利用WWP线搜索(18) (19)得到步长${\alpha }_k$ ；

Step 5：计算$x_{k+1}=x_k+{\alpha }_k{d}_k$ ，计算$g_{k+1}=\nabla f\left(x_{k+1}\right)$ ；

Step 6：计算$s_k=x_{k+1}-x_k$ ，$y_k=g_{k+1}-g_k$ ，$M_k$ ，${\theta }_k$ ，由(12)式求得${\widehat{y}}_k$ ；若$k<m$ ，则储存曲率对${\left\{\left(s_i,{\widehat{y}}_i\right)\right\}}_{i=0}^k$ ，若$k\ge m$ ，保留最近的${\left\{\left(s_i,{\widehat{y}}_i\right)\right\}}_{i=k-m+1}^k$ ，删除已经储存的$\left(s_{k-m},{\widehat{y}}_{k-m}\right)$ ；

Step 7：由(15)式计算得到$B_{k+1}$ ；

Step 8：$k=k+1$ ，转入Step 2。

D-LBFGS算法与LBFGS算法相比，在构造曲率对时，不仅用到了相邻两次迭代点之间的梯度差，利用了迭代点之间的差值，还构造了关于相邻两次梯度差的对角矩阵，因此，在逆Hesse矩阵的构造中，该方法更充分地利用了目标函数在当前迭代点的曲率信息，这种优势使得逆Hesse矩阵的近似更为准确，从而有效提高了优化算法的性能。

与BFGS算法相比，在计算Hesse的逆的近似矩阵的过程中，不需要储存全部的曲率对和计算得到的近似矩阵信息，只需要储存2个$m\times n$ 的矩阵，以较小的存储成本获得相对高的计算效率。

5. D-LBFGS算法收敛性分析

下面研究D-LBFGS算法的收敛性质，在拟牛顿算法的收敛性的理论证明基础之上，给出本文的假设条件：

假设5.1 目标函数$f:R^n\to R$ 是二次连续可微的，$x_0$ 是初始迭代点。$f\left(x\right)$ 在$R^n$ 上有下界，$x^{*}$ 是$f\left(x\right)$ 的极小点。

假设5.2 $f\left(x\right)$ 是一致凸函数，即对所有$z\in R^n$ ，存在正常数m、M，使得

$m{\Vert z\Vert }^2\le z^{\text{T}}{G}_k\left(x\right)z\le M{\Vert z\Vert }^2.$ (20)

假设5.3 $\nabla f\left(x\right)$ 是全局Lipschitz连续的，即存在$L>0$ ，使得对于任意的$x,y\in R^n$ ，都有

$\Vert \nabla f\left(x\right)-\nabla f\left(y\right)\Vert \le L\Vert x-y\Vert .$ (21)

假设5.4 存在常数$p>0$ ，$q>0$ ，使得对$\forall k$ 都满足

$pI\underset{\_}{\prec }{M}_k\underset{\_}{\prec }qI.$ (22)

假设5.5 ${\displaystyle \sum _{k=1}^{+\infty }{\alpha }_k}=+\infty$ ，${\displaystyle \sum _{k=1}^{+\infty }{\alpha }_k^2}<+\infty$ 。

定理5.1 如果$H_k$ 和$B_k$ 正定，则由D-LBFGS算法构造出的$H_{k+1}$ 和$B_{k+1}$ 也是正定的。

证明：用归纳法证明。

在算法2中，$B_0=I$ 是正定矩阵。设$B_k$ 是正定的，下面证明$B_{k+1}$ 也是正定的。

对任意的非零向量$z\in R^n$ ，有

$\begin{array}{c}{z}^{\text{T}}{B}_{k+1}z=z^{\text{T}}{B}_{k}z+\frac{z^{\text{T}}{\widehat{y}}_k{\widehat{y}}_k^{\text{T}}z}{s_k^{\text{T}}{\widehat{y}}_k}-\frac{z^{\text{T}}{B}_k{s}_k{s}_k^{\text{T}}{B}_{k}z}{s_k^{\text{T}}{B}_k{s}_k}\\ =z^{\text{T}}{B}_{k}z+\frac{{\left(z^{\text{T}}{\widehat{y}}_k\right)}^2}{s_k^{\text{T}}{\widehat{y}}_k}-\frac{{\left(z^{\text{T}}{B}_k{s}_k\right)}^2}{s_k^{\text{T}}{B}_k{s}_k}.\end{array}$ (23)

因为$B_k$ 是正定矩阵，因此存在正定矩阵$B_k^{\frac{1}{2}}$ ，使得$B_k=B_k^{\frac{1}{2}}{B}_k^{\frac{1}{2}}$ ，令$u=B_k^{\frac{1}{2}}z$ ，$s=B_k^{\frac{1}{2}}{s}_k$ ，则(23)式可写成

$\begin{array}{c}{z}^{\text{T}}{B}_{k+1}z=u^{\text{T}}u-\frac{{\left(u^{\text{T}}s\right)}^2}{s^{\text{T}}s}+\frac{{\left(z^{\text{T}}{\widehat{y}}_k\right)}^2}{s_k^{\text{T}}{\widehat{y}}_k}\\ =\frac{\left(u^{\text{T}}u\right)\left(s^{\text{T}}s\right)-{\left(u^{\text{T}}s\right)}^2}{\left(s^{\text{T}}s\right)}+\frac{{\left(z^{\text{T}}{\widehat{y}}_k\right)}^2}{s_k^{\text{T}}{\widehat{y}}_k}.\end{array}$

根据Schwartz不等式，有$\left(u^{\text{T}}u\right)\left(s^{\text{T}}s\right)\ge {\left(u^{\text{T}}s\right)}^2$ ，因此得到

$\frac{\left(u^{\text{T}}u\right)\left(s^{\text{T}}s\right)-{\left(u^{\text{T}}s\right)}^2}{\left(s^{\text{T}}s\right)}\ge 0,$

故第1项非负。

当${\theta }_k=1$ 时，由(16)式得到，${\widehat{y}}_k=y_k$ ，由BFGS算法的理论分析可知$s_k^{\text{T}}{y}_k>0$ ，故此时$B_{k+1}$ 为正定矩阵[11]。

当${\theta }_k\ne 1$ 时，

$s_k^{\text{T}}{\widehat{y}}_k={\theta }_k{s}_k^{\text{T}}{y}_k+\left(1-{\theta }_k\right)s_k^{\text{T}}{M}_k{s}_k.$

${\theta }_k=\frac{0.8s_k^{\text{T}}{M}_k{s}_k}{s_k^{\text{T}}{M}_k{s}_k-s_k^{\text{T}}{y}_k}.$

可以得到

${\theta }_k\left(s_k^{\text{T}}{M}_k{s}_k-s_k^{\text{T}}{y}_k\right)=0.8s_k^{\text{T}}{M}_k{s}_k.$

$\begin{array}{c}{\theta }_k{s}_k^{\text{T}}{y}_k={\theta }_k{s}_k^{\text{T}}{M}_k{s}_k-0.8s_k^{\text{T}}{M}_k{s}_k+s_k^{\text{T}}{M}_k{s}_k-s_k^{\text{T}}{M}_k{s}_k\\ =\left({\theta }_k-1\right)s_k^{\text{T}}{M}_k{s}_k+0.2s_k^{\text{T}}{M}_k{s}_k.\end{array}$

$0.2s_k^{\text{T}}{M}_k{s}_k={\theta }_k{s}_k^{\text{T}}{y}_k+\left(1-{\theta }_k\right)s_k^{\text{T}}{M}_k{s}_k=s_k^{\text{T}}{\widehat{y}}_k.$

由假设5.4可知$s_k^{\text{T}}{\widehat{y}}_k>0$。所以$B_{k+1}$ 为正定矩阵。

在算法2中，$H_0=I$ 是正定矩阵。设$H_k$ 是正定的，下面证明$H_{k+1}$ 也是正定的。由(16)式，对任意的非零向量$z\in R^n$ ，有

$\begin{array}{c}{z}^{\text{T}}{H}_{k+1}z=z^{\text{T}}\left(I-\frac{s_k{\widehat{y}}_k^{\text{T}}}{{\widehat{y}}_k^{\text{T}}{s}_k}\right)H_k\left(I-\frac{s_k{\widehat{y}}_k^{\text{T}}}{{\widehat{y}}_k^{\text{T}}{s}_k}\right)z+z^{\text{T}}\frac{s_k{s}_k^{\text{T}}}{{\widehat{y}}_k^{\text{T}}{s}_k}z\\ =z^{\text{T}}\left(I-\frac{s_k{\widehat{y}}_k^{\text{T}}}{{\widehat{y}}_k^{\text{T}}{s}_k}\right)H_k\left(I-\frac{s_k{\widehat{y}}_k^{\text{T}}}{{\widehat{y}}_k^{\text{T}}{s}_k}\right)z+\frac{{\left(s_k^{\text{T}}z\right)}^2}{{\widehat{y}}_k^{\text{T}}{s}_k}\\ >0\end{array}$

因此$H_{k+1}$ 为正定矩阵。

定理5.2 若假设5.1~5.3成立，$\left\{x_k\right\}$ 是由D-LBFGS算法产生的序列，则$\underset{k\to \infty }{\lim }\Vert s_k\Vert =0$ 。

证明：由Taylor公式展开可得：

$f\left(x_{k+1}\right)=f\left(x_k\right)+g_k^{\text{T}}{s}_k+{\displaystyle {\int }_0^1\left(1-\phi \right)s_k^{\text{T}}G\left(x_k+\phi s_k\right)s_k\text{d}\phi }.$

由假设条件可知，

${\displaystyle {\int }_0^1\left(1-\phi \right)s_k^{\text{T}}G\left(x_k+\phi s_k\right)s_k\text{d}\phi }\ge {\displaystyle {\int }_0^1\left(1-\phi \right)m{\Vert s_k\Vert }^2\text{d}\phi }=\frac{m{\Vert s_k\Vert }^2}{2}.$

由WWP线搜索可知$f\left(x_k+{\alpha }_k{d}_k\right)\le f\left(x_k\right)+\delta {\alpha }_k{g}_k^{\text{T}}{d}_k$ ，$0<\delta <1$ ，由此可得

${\displaystyle {\int }_0^1\left(1-\phi \right)s_k^{\text{T}}G\left(x_k+\phi s_k\right)s_k\text{d}\phi }=f\left(x_{k+1}\right)-f\left(x_k\right)-g_k^{\text{T}}{s}_k\le \left(1-\delta \right)\left(-g_k^{\text{T}}{s}_k\right).$

又因为

$-g_k^{\text{T}}{s}_k={\alpha }_k{g}_k^{\text{T}}{H}_k{g}_k>0,$

${\displaystyle \sum _{k=0}^{+\infty }\left(-g_k^{\text{T}}{s}_k\right)}={\displaystyle \sum _{k=0}^{+\infty }\left(g_k^{\text{T}}{\alpha }_k{d}_k\right)}\le \frac{1}{\delta }{\displaystyle \sum _{k=0}^{+\infty }\left[f\left(x_k\right)-f\left(x_{k+1}\right)\right]}\le \frac{1}{\delta }\left[f\left(x_0\right)-f\left(x^{*}\right)\right].$

所以

${\displaystyle \sum _{k=0}^{+\infty }{\Vert s_k\Vert }^2}\le \frac{2}{m}{\displaystyle \sum _{k=1}^{+\infty }\left(1-\delta \right)}\left(-g_k^{\text{T}}{s}_k\right).$

故函数项级数${\displaystyle \sum _{k=1}^{+\infty }{\Vert s_k\Vert }^2}$ 收敛，故

$\underset{k\to \infty }{\lim }\Vert s_k\Vert =0.$

定理5.3 若假设5.1~5.5成立，$\left\{x_k\right\}$ 是由D-LBFGS算法产生的序列，则$\underset{k\to \infty }{\lim \inf }\Vert g_k\Vert =0$ 。

证明：由$d_k=-H_k{g}_k$ ，因为$H_k$ 正定，故有

$g_k^{\text{T}}{d}_k=-g_k^{\text{T}}{H}_k{g}_k\le 0.$

由WWP线搜索可知$f\left(x_k+{\alpha }_k{d}_k\right)\le f\left(x_k\right)+\delta {\alpha }_k{g}_k^{\text{T}}{d}_k$ ，即

$f\left(x_k\right)-f\left(x_k+{\alpha }_k{d}_k\right)\ge -\delta {\alpha }_k{g}_k^{\text{T}}{d}_k\ge 0.$

则$\left\{f\left(x_k\right)\right\}$ 是递减序列，又由假设5.1，则f是有界的。假设有收敛序列$\left\{x_{n_k}\right\}$ ，存在正整数K，当$n_i,n_m>K$ 时，有

$\underset{n_i,n_m\to \infty }{\lim }\Vert x_{n_i}-x_{n_m}\Vert =0.$

又由假设5.3可知，$\nabla f$ 是Lipschitz连续的，因此

$\Vert \nabla f\left(x_{n_i}\right)-\nabla f\left(x_{n_m}\right)\Vert \le L\Vert x_{n_i}-x_{n_m}\Vert .$

因此，$\underset{n_i,n_m\to \infty }{\lim \inf }\Vert \nabla f\left(x_{n_i}\right)-\nabla f\left(x_{n_m}\right)\Vert =0$ ，即$\underset{k\to \infty }{\lim \inf }\Vert g_k\Vert =0$ 。

6. 使用D-LBFGS—CNN算法解决无创血糖浓度的预测问题

BP神经网络是一种按照误差反向传播算法训练的多层前馈神经网络，有输入层、隐藏层和输出层，包括前向传播、反向传播两个部分，通过前向传播计算误差，反向传播更新权值和阈值。早在1983年，Dennis和Schnabel就将牛顿法与神经网络相结合，提出了基于牛顿法的前馈型神经网络算法[13]。

6.1. D-LBFGS—CNN算法

在BP神经网络的反向传播过程中，参数更新主要基于梯度下降类算法，记每一层的权重$\omega$ 和偏置b为参数向量θ，在迭代的每一轮中对参数θ进行更新。设E为损失函数，${\theta }_t$ 为第t次迭代参数。参数θ的更新公式一般记为：${\theta }_{t+1}={\theta }_t+\alpha g_t$ ，其中，$g_t$ 为损失函数E关于θ的负梯度，$\alpha$ 为学习率。通过大量的样本对BP神经网络进行训练，使得预测输出与实际标签间的误差越来越小。

Figure 1. D-LBFGS—CNN algorithm flow chart

图1. D-LBFGS—CNN算法流程图

下面提出了一种利用D-LBFGS算法对神经网络中的权值阈值参数进行优化的算法，记为D-LBFGS—CNN算法，在反向更新过程中使用D-LBFGS算法更新θ，D-LBFGS—CNN算法的计算步骤如下：

Step 1：初始化BP神经网络结构及超参数，包括网络结构参数、步长等，初始权值阈值参数${\theta }_0$ ；初始化D-LBFGS算法中的参数，包括有限记忆数m、$\delta$ 、$\sigma$ 等。令$t=0$ 。

Step 2：输入训练集，使用当前已知固定权值阈值参数${\theta }_t$ 的神经网络计算损失函数值，检查损失是否满足设定的终止条件。如果满足条件，则输出当前计算的权值阈值，否则，进行Step 3。

Step 3：把神经网络的损失作为拟牛顿法优化的目标函数，以初始参数${\theta }_t$ 开始，使用D-LBFGS算法进行迭代计算，得到最优参数解记为${\theta }_{t+1}$ 。令$t=t+1$ ，返回Step 2继续训练网络向前传播。

Step 4：保存在Step 2或Step 3中得到的最优权值和阈值参数，得到最优神经网络模型。

D-LBFGS—CNN算法流程图如图1。

6.2. 实验数据

糖尿病是一种由于胰岛素分泌不足或身体无法有效利用胰岛素引起的内分泌紊乱类疾病。目前医疗手段尚无法完全治愈糖尿病，但实时监测血糖浓度并及时调整药物剂量是减缓并发症发生最有效的方法之一。传统血糖监测分为有创和微创两种，但频繁的有创监测会给患者带来不便，例如指尖采血可能导致指尖硬化或感染。因此，无创检测成为一种理想选择。无创检测速度快、方便，并吸引了学者们的关注。手指因易于取样、毛细血管丰富成为无创检测部位的首选。

光谱分析是无创检测中备受关注的方法，特别是近红外光谱法[14]。它有操作简便、提供信息量大、速度快等优点。实验中的光谱数据是使用Hyperspec^TM NIP近红外光谱扫描成像仪对手指指腹进行采集的，波长范围为955 nm至1700 nm，积分时间为35 ms，光谱分辨率为5 nm。血糖浓度采用标准化酶法检测，符合CE标准。口服葡萄糖耐量实验(OGTT)用于测量被试者的血糖浓度，评估其血糖调节能力。血糖浓度呈现先升后降的变化趋势，每10分钟采集血糖样本，以光谱数据预测血糖浓度。实验包含多次OGTT，通过口服葡萄糖溶液和手指指腹位置的光谱扫描获取了光谱和血糖浓度数据。通过测量测试者血液中39个不同的血糖浓度，获得对应于每个浓度的原始光谱信号，用近红外光谱仪测定指尖的局部区域，获得100帧近红外光谱图像，每个浓度提取100帧图像，以及光谱对应的总共1500个波长吸光度。

利用近红外光谱法和血糖样本测得的真实血糖数据获取了光谱和血糖浓度数据，下面使用D-LBFGS—CNN算法通过光谱数据来预测血糖浓度。使用80%的血糖数据进行训练，剩下的20%用于测试。主要网络结构为三层神经网络，分别是输入层(150个节点)、隐藏层(1000个节点)、输出层(39个节点)，激活函数采用ReLU，损失函数采用交叉熵函数，迭代次数设定为500次。拟牛顿算法优化网络权值阈值参数时，参数设置为$\epsilon ={10}^{-2}$ ，$\delta ={10}^{-4}$ ，$\sigma =0.9$ ，使用WWP线搜索条件得到步长。

在编程中，选择固定的随机种子数1188，通过固定随机种子，可以确保在神经网络梯度下降算法中初始的权重值是一致的[15]。这意味着，多次运行，网络有相同的初始权重，使得训练过程的初始状态具有可重复性，有助于验证模型性能和调试代码，增加实验结果可信度。

该模型程序使用Python3.9.0编写，在搭载AMD Ryzen 5 6600H with Radeon Graphics处理器，16 GB内存，运行Windows 11操作系统的电脑上执行。

实验结果

在上述D-LBFGS—CNN算法步骤中，反向更新过程分别使用LBFGS算法、阻尼牛顿法、文献[13]的LBFGS算法和D-LBFGS算法，对数据进行预测，如图2，通过对应的训练损失，测试损失和测试b精度结果比较算法的性能。

图2(a)表示四种算法在血糖数据集上的训练损失，图2(b)表示四种算法在血糖数据集上的测试损失；由图2可知在训练500个迭代周期时，LBFGS算法、阻尼牛顿法、文献[13]的LBFGS算法和D-LBFGS算法的训练损失分别是0.0155、0.0240、0.0188和0.0086；在测试500个迭代周期时，LBFGS算法、阻尼牛顿法、文献[13]的LBFGS算法和D-LBFGS算法的测试损失分别是0.0104、0.0194、0.0147和0.0078。

(a) 训练损失 (b) 测试损失

Figure 2. The losses of the four algorithms on the data set

图2. 四种算法在数据集上的损失

由图3可以看出，D-LBFGS算法优化神经网络可以得到较好的训练测试效果，测试集精度较高。相较于LBFGS算法在第474次Epoch时达到99.99%的测试精度，D-LBFGS算法在第351次Epoch时就已经达到所需要的测试精度。

Figure 3. Test accuracy of four algorithms on data sets

图3. 四种算法在数据集上的测试精度

由表1可以看出，四种算法在固定50次迭代次数间隔时的预测精度，D-LBFGS算法优化神经网络得到较好的训练测试效果，测试集精度达到99.999%。D-LBFGS算法在第351次Epoch时就已经达到所需要的测试精度，优于LBFGS算法在474次迭代精度99.999%，阻尼牛顿法精度99.359%和文献[13]中LBFGS算法精度99.359%。由此可以看出，D-LBFGS算法在更少的迭代次数内就达到了预期的测试精度，有更快的收敛速度。

Table 1. Prediction accuracy of four algorithms when the number of iterations is fixed

表1. 四种算法固定迭代次数时的预测精度

迭代次数	LBFGS	阻尼牛顿法	文献[13]的LBFGS	D-LBFGS
50	90.000	85.769	89.231	89.487
100	92.821	89.878	94.231	93.205
150	95.385	92.564	96.154	95.385
200	97.693	95.128	98.205	98.590
250	98.334	97.821	97.949	98.590
300	99.103	97.821	98.718	99.615
350	99.744	98.462	99.359	99.872
400	98.872	99.103	99.487	99.615
450	99.615	98.718	98.615	99.872
500	99.999	99.359	99.359	99.999

Figure 4. Heat map of non-invasive blood glucose concentration prediction results

图4. 无创血糖浓度预测结果热力图

通过热力图可以分析数据集中各个特征之间的相似度，图4清晰地呈现了经D-LBFGS算法优化后的BP神经网络对测试集测试结果与实际血糖值之间的相似程度，显示了高达99.99%的测试准确率，凸显了该优化模型在血糖浓度测试中的出色表现。这样的准确率表明模型对于预测血糖浓度具有相当高的可靠性，对未知数据的预测也具备很好的实用性。

7. 结论

基于阻尼技术和有限记忆的思想，通过利用相邻两次迭代点之间的梯度差和变化量的线性组合，实现曲率对的更新，以改进的校正公式为基础提出了D-LBFGS算法，研究基于WWP线搜索条件证明了算法的全局收敛性。提出利用D-LBFGS算法优化神经网络权值阈值参数的算法，记D-LBFGS—CNN算法，用于解决无创血糖浓度预测问题，并对计算结果进行了分析。研究结果表明，算法在无创血糖浓度预测问题中具有较好的测试成果，训练速度和准确率都有相应的提升，是有效且可行的。本文提出的D-LBFGS算法是以目标函数都为凸函数的基础上展开的，对于目标函数非凸的情况尚未涉及。以及参数的设定多参考文献中其他学者的设置条件或多次实验来选择的，实验中可能有一定的盲目性，可以结合其他优化算法，可能有更好的实验效果。

基金项目

吉林省自然科学基金(NO. YDZJ202201ZYTS519)。

NOTES

^*通讯作者。