数值优化 | 详解拟牛顿法之SR1、DFP、BFGS、L-BFGS(附案例分析与Python实现)

0 专栏介绍

🔥课设、毕设、创新竞赛必备！🔥本专栏涉及更高阶的运动规划算法轨迹优化实战，包括：曲线生成、碰撞检测、安全走廊、优化建模(QP、SQP、NMPC、iLQR等)、轨迹优化(梯度法、曲线法等)，每个算法都包含代码实现加深理解

🚀详情：运动规划实战进阶：轨迹优化篇

---

1 拟牛顿法框架

在数值优化 | 图解牛顿法、阻尼牛顿法与高斯牛顿法(附案例分析与Python实现)中我们介绍了经典牛顿法的迭代公式

$${\mathbf{x}}_{k+1}={\mathbf{x}}_k-{\mathbf{H}}^{-1}\left({\mathbf{x}}_k\right)\nabla F\left({\mathbf{x}}_k\right)$$

可以看到，经典牛顿法需要对目标函数求Hessian矩阵的逆矩阵，矩阵求逆的时间复杂度是$O(n^3)$，其中$n$为矩阵维度，因此牛顿法计算复杂度高、数值稳定性较差。同时，无法保证对于任意$f\left(\mathbf{x}\right)$的Hessian矩阵可逆且正定。

为了克服牛顿法的缺陷，引入拟牛顿法(quasi-Newton method)：根据迭代点处的梯度信息求$\mathbf{H}\left(\mathbf{x}\right)$的近似值${\mathbf{B}}_k\left(\mathbf{x}\right)$或${\mathbf{H}}^{-1}\left(\mathbf{x}\right)$的近似值${\mathbf{H}}_k\left(\mathbf{x}\right)$来降低计算量，算法流程如图所示。根据近似矩阵更新算法的不同，拟牛顿法有若干子算法

2 SR1算法推导

秩一公式(Symmetric Rank-One, SR1)是结构最简单的拟牛顿矩阵更新算法，考虑

$${\mathbf{H}}_{k+1}={\mathbf{H}}_k+\Delta {\mathbf{H}}_k$$

期望${\mathbf{H}}_{k+1}\approx \mathbf{H}\left({\mathbf{x}}_{k+1}\right)$、${\mathbf{H}}_k\approx \mathbf{H}\left({\mathbf{x}}_k\right)$，所以$\Delta {\mathbf{H}}_k\approx \mathbf{H}\left({\mathbf{x}}_{k+1}\right)-\mathbf{H}\left({\mathbf{x}}_k\right)$。由于Hessian矩阵对称，所以$\Delta {\mathbf{H}}_k$也应对称，不妨令$\Delta {\mathbf{H}}_k=\beta {\mathbf{uu}}^T$；接着，对$\nabla F\left(\mathbf{x}\right)$在${\mathbf{x}}_{k+1}$处一阶泰勒展开可得

$$\nabla F\left(\mathbf{x}\right)=\nabla F\left({\mathbf{x}}_{k+1}\right)+\mathbf{H}\left({\mathbf{x}}_{k+1}\right)\left(\mathbf{x}-{\mathbf{x}}_{k+1}\right)$$

记${\mathbf{s}}_k={\mathbf{x}}_{k+1}-{\mathbf{x}}_k$，${\mathbf{y}}_k=\nabla F\left({\mathbf{x}}_{k+1}\right)-\nabla F\left({\mathbf{x}}_k\right)$，则

$${\mathbf{y}}_k={\mathbf{B}}_{k+1}{\mathbf{s}}_k\Rightarrow {\mathbf{H}}_{k+1}{\mathbf{y}}_k={\mathbf{s}}_k$$

注意到${\mathbf{s}}_k^T{\mathbf{B}}_{k+1}{\mathbf{s}}_k={\mathbf{s}}_k^T{\mathbf{y}}_k$，为了使${\mathbf{H}}_k$或${\mathbf{B}}_k$正定，需要在迭代过程中保证${\mathbf{s}}_k^T{\mathbf{y}}_k>0$。基于上式，在$\Delta {\mathbf{H}}_k=\beta {\mathbf{uu}}^T$两边同时右乘${\mathbf{y}}_k$

$$\mathbf{u}=\gamma \left({\mathbf{s}}_k-{\mathbf{H}}_k{\mathbf{y}}_k\right)$$

其中

$$\gamma =\frac{1}{\beta {\mathbf{u}}^T{\mathbf{y}}_k}$$

将$\mathbf{u}=\gamma \left({\mathbf{s}}_k-{\mathbf{H}}_k{\mathbf{y}}_k\right)$反代回${\mathbf{H}}_{k+1}{\mathbf{y}}_k={\mathbf{H}}_k{\mathbf{y}}_k+\beta {\mathbf{uu}}^T{\mathbf{y}}_k$可得

$${\mathbf{s}}_k-{\mathbf{H}}_k{\mathbf{y}}_k=\beta {\gamma }^2{\left({\mathbf{s}}_k-{\mathbf{H}}_k{\mathbf{y}}_k\right)}^T{\mathbf{y}}_k\left({\mathbf{s}}_k-{\mathbf{H}}_k{\mathbf{y}}_k\right)$$

显见$\beta {\gamma }^2=1/{\left({\mathbf{s}}_k-{\mathbf{H}}_k{\mathbf{y}}_k\right)}^T{\mathbf{y}}_k$，从而得到SR1的更新公式

$${\mathbf{H}}_{k+1}={\mathbf{H}}_k+\frac{\left({\mathbf{s}}_k-{\mathbf{H}}_k{\mathbf{y}}_k\right){\left({\mathbf{s}}_k-{\mathbf{H}}_k{\mathbf{y}}_k\right)}^T}{{\left({\mathbf{s}}_k-{\mathbf{H}}_k{\mathbf{y}}_k\right)}^T{\mathbf{y}}_k}$$

代码实现为

gx = g_func(x)
d = (H @ gx).squeeze()
alpha = self.line_searcher_.search(func, g_func, x, -d)

x_next = x - alpha * d
s = x_next - x
y = g_func(x_next) - gx

item = s - H @ y
H = H + item @ item.T / (item.T @ y)

3 DFP算法推导

DFP算法是SR1算法的拓展版本，引入了更大的自由度

$${\mathbf{H}}_{k+1}={\mathbf{H}}_k+\beta {\mathbf{uu}}^T+\gamma {\mathbf{vv}}^T$$

与SR1算法推导类似，方程两边同时右乘${\mathbf{y}}_k$得到

$${\mathbf{s}}_k={\mathbf{H}}_{k+1}{\mathbf{y}}_k={\mathbf{H}}_k{\mathbf{y}}_k+\beta {\mathbf{uu}}^T{\mathbf{y}}_k+\gamma {\mathbf{vv}}^T{\mathbf{y}}_k$$

令$\mathbf{u}={\mathbf{s}}_k$、$\mathbf{v}={\mathbf{H}}_k{\mathbf{y}}_k$，则${\mathbf{s}}_k-{\mathbf{H}}_k{\mathbf{y}}_k=\beta {\mathbf{s}}_k{\mathbf{s}}_k^T{\mathbf{y}}_k+\gamma {\mathbf{H}}_k{\mathbf{y}}_k{\mathbf{y}}_k^T{\mathbf{H}}_k^T{\mathbf{y}}_k$，从而

$$\{\begin{array}{l}\beta =\frac{1}{{\mathbf{s}}_k^T{\mathbf{y}}_k}\\ \gamma =-\frac{1}{{\mathbf{y}}_k^T{\mathbf{H}}_k^T{\mathbf{y}}_k}\end{array}$$

得到DFP的更新公式

$${\mathbf{H}}_{k+1}={\mathbf{H}}_k+\frac{{\mathbf{s}}_k{\mathbf{s}}_k^T}{{\mathbf{s}}_k^T{\mathbf{y}}_k}-\frac{{\mathbf{H}}_k{\mathbf{y}}_k{\mathbf{y}}_k^T{\mathbf{H}}_k^T}{{\mathbf{y}}_k^T{\mathbf{H}}_k^T{\mathbf{y}}_k}$$

代码实现为

gx = g_func(x)
d = (H @ gx).squeeze()
alpha = self.line_searcher_.search(func, g_func, x, -d)

x_next = x - alpha * d
s = x_next - x
y = g_func(x_next) - gx

sTy = float(s.T @ y)
H = H + s @ s.T / sTy - H @ y @ y.T @ H.T / float(y.T @ H.T @ y)

4 BFGS算法推导

BFGS算法是DFP的对偶版本，假设

$${\mathbf{B}}_{k+1}={\mathbf{B}}_k+\beta {\mathbf{uu}}^T+\gamma {\mathbf{vv}}^T$$

类比DFP推导过程可得到

$${\mathbf{B}}_{k+1}={\mathbf{B}}_k+\frac{{\mathbf{y}}_k{\mathbf{y}}_k^T}{{\mathbf{y}}_k^T{\mathbf{s}}_k}-\frac{{\mathbf{B}}_k{\mathbf{s}}_k{\mathbf{s}}_k^T{\mathbf{B}}_k^T}{{\mathbf{s}}_k^T{\mathbf{B}}_k^T{\mathbf{s}}_k}$$

为了进一步计算${\mathbf{H}}_k={\mathbf{B}}_k^{-1}$，应用下面的引理

Shermann-Morrison-Woodbury公式 设$\mathbf{A}$是$n\times n$可逆矩阵，向量$\mathbf{u}$、$\mathbf{v}\in {\mathbb{R}}^n$ ，则$\mathbf{A}+{\mathbf{uv}}^T$可逆当且仅当$1+{\mathbf{v}}^T{\mathbf{A}}^{-1}\mathbf{u}\ne 0$，且

$${\left(\mathbf{A}+{\mathbf{uv}}^T\right)}^{-1}={\mathbf{A}}^{-1}-\frac{{\mathbf{A}}^{-1}{\mathbf{uv}}^T{\mathbf{A}}^{-1}}{1+{\mathbf{v}}^T{\mathbf{A}}^{-1}\mathbf{u}}$$

从而得到BFGS更新公式

$${\mathbf{H}}_{k+1}=\left(\mathbf{I}-\frac{{\mathbf{s}}_k{\mathbf{y}}_k^T}{{\mathbf{s}}_k^T{\mathbf{y}}_k}\right){\mathbf{H}}_k\left(\mathbf{I}-\frac{{\mathbf{y}}_k{\mathbf{s}}_k^T}{{\mathbf{s}}_k^T{\mathbf{y}}_k}\right)+\frac{{\mathbf{s}}_k{\mathbf{s}}_k^T}{{\mathbf{s}}_k^T{\mathbf{y}}_k}$$

代码实现为

gx = g_func(x)
d = (H @ gx).squeeze()
alpha = self.line_searcher_.search(func, g_func, x, -d)

x_next = x - alpha * d
s = x_next - x
y = g_func(x_next) - gx

sTy = float(s.T @ y)
H = (I - s @ y.T / sTy) @ H @ (I - y @ s.T / sTy) + s @ s.T / sTy

5 L-BFGS算法推导

在BFGS算法中，每次迭代需要上一轮的${\mathbf{H}}_k$，而保存${\mathbf{H}}_k$至少需要$O(n^2)$的空间复杂度，这对于大规模数值优化问题而言无法接受。L-BFGS对BFGS的内存占用率进行了限制，记${\rho }_k=1/{\mathbf{s}}_k^T{\mathbf{y}}_k$、${\mathbf{V}}_k=\mathbf{I}-{\rho }_k{\mathbf{y}}_k{\mathbf{s}}_k^T$，则递推地

$$\begin{array}{rl}{\mathbf{H}}_{k+1}& ={\mathbf{V}}_k^T{\mathbf{H}}_k{\mathbf{V}}_k+{\rho }_k{\mathbf{s}}_k{\mathbf{s}}_k^T\\ & =\left({\mathbf{V}}_k^T\cdots {\mathbf{V}}_0^T\right){\mathbf{H}}_0\left({\mathbf{V}}_0\cdots {\mathbf{V}}_k\right)+\left({\mathbf{V}}_k^T\cdots {\mathbf{V}}_1^T\right){\rho }_1{\mathbf{s}}_1{\mathbf{s}}_1^T\left({\mathbf{V}}_1\cdots {\mathbf{V}}_k\right)+\cdots +{\rho }_k{\mathbf{s}}_k{\mathbf{s}}_k^T\end{array}$$

仅保留最近$m$步的结果

$$\begin{array}{rl}{\mathbf{H}}_{k+1}& =\left({\mathbf{V}}_k^T\cdots {\mathbf{V}}_{k-m+1}^T\right){\mathbf{H}}_0\left({\mathbf{V}}_{k-m+1}\cdots {\mathbf{V}}_k\right)\\ & +\left({\mathbf{V}}_k^T\cdots {\mathbf{V}}_{k-m+2}^T\right){\rho }_{k-m+1}{\mathbf{s}}_{k-m+1}{\mathbf{s}}_{k-m+1}^T\left({\mathbf{V}}_{k-m+2}\cdots {\mathbf{V}}_k\right)\\ & +\cdots \\ & +{\rho }_k{\mathbf{s}}_k{\mathbf{s}}_k^T\end{array}$$

因此在L-BFGS算法中无需保留完整的${\mathbf{H}}_k$，只需存储向量序列 { s i } i = k − m k \left\{ \boldsymbol{s}_i \right\} _{i=k-m}^{k} {si​}i=k−mk​、 { y i } i = k − m k \left\{ \boldsymbol{y}_i \right\} _{i=k-m}^{k} {yi​}i=k−mk​，占用存储空间$2mn$，空间复杂度为$O(n)$，下降了一个量级

代码实现为

 gx = g_func(x)

# two-iteration
a = []
d = g_func(x).reshape(-1, 1)
for i in range(len(S) - 1, -1, -1):
    s_i, y_i, p_i = S[i], Y[i], P[i]
    a_i = p_i * s_i.T @ d
    d = d - a_i * y_i
    a.append(a_i)
for i in range(len(S)):
    s_i, y_i, p_i = S[i], Y[i], P[i]
    a_i = a[len(S) - i - 1]
    beta = p_i * y_i.T @ d
    d = d + s_i * (a_i - beta)
d = d.squeeze()
alpha = self.line_searcher_.search(func, g_func, x, -d)

x_next = x - alpha * d

s = (x_next - x).reshape(-1, 1)
y = (g_func(x_next) - gx).reshape(-1, 1)

if len(S) == self.m:
    S.pop(0)
S.append(s)

if len(Y) == self.m:
    Y.pop(0)
Y.append(y)

if len(P) == self.m:
    P.pop(0)
P.append(1.0 / (s.T @ y))

6 案例分析

以曲线拟合为例展示三种牛顿类优化算法的效果

考虑如下方程：
$$\mathbf{y}=e^{a{\mathbf{x}}^2+b\mathbf{x}+c}+\mathbf{w}$$
其中$a$、$b$、$c$是待估计的曲线参数，$\mathbf{w}$是高斯噪声，假设有$n$个观测样本$(x_i,y_i)$，希望根据这组样本得到最佳的曲线参数，构造如下的优化问题
$$F(\mathbf{\theta })=\min \frac{1}{2n}\sum _{i=1}^n(y_i-e^{a{x}_i^2+b{x}_i+c}{)}^2$$

下面来看结果

• SR1算法

  ===========================================================
  Method: SR1
  Iterations: 10, Using time: 0.27718114852905273
  Start Point: [ 2 -1  5]
  Start Value: 16427.297702943684
  End Point: [1.58470727 1.19606329 1.24535202]
  End Value: 0.492582013221927
  ===========================================================
  

• DFP算法

  ===========================================================
  Method: DFP
  Iterations: 68, Using time: 1.8977277278900146
  Start Point: [ 2 -1  5]
  Start Value: 16427.297702943684
  End Point: [1.32228464 1.58602539 1.11512083]
  End Value: 0.4622062689277505
  ===========================================================
  

• BFGS算法

  ===========================================================
  Method: BFGS
  Iterations: 19, Using time: 0.5240802764892578
  Start Point: [ 2 -1  5]
  Start Value: 16427.297702943684
  End Point: [1.32275944 1.58517781 1.11542359]
  End Value: 0.46220542143350907
  ===========================================================
  

  

• L-BFGS

  ===========================================================
  Method: LBFGS
  Iterations: 19, Using time: 0.5343437194824219
  Start Point: [ 2 -1  5]
  Start Value: 16427.297702943684
  End Point: [1.32204779 1.58592409 1.11533423]
  End Value: 0.46220620676321367
  ===========================================================
  

从结果来看，由于LBFGS在内存和计算上都有显著优势，因此它在许多实际问题中（如大规模数据优化、深度学习训练等）得到了广泛应用。它是许多优化库（如SciPy、TensorFlow、PyTorch等）中常用的优化算法之一

完整工程代码请联系下方博主名片获取

---

🔥 更多精彩专栏：

• 《ROS从入门到精通》
• 《Pytorch深度学习实战》
• 《机器学习强基计划》
• 《运动规划实战精讲》
• …

👇源码获取 · 技术交流 · 抱团学习 · 咨询分享 请联系👇

原文地址：https://blog.csdn.net/FRIGIDWINTER/article/details/144013803

免责声明：本站文章内容转载自网络资源，如本站内容侵犯了原著者的合法权益，可联系本站删除。更多内容请关注自学内容网（zxcms.com）！