拟牛顿法（DFP、BFGS、L-BFGS）

拟牛顿法

一、牛顿法

1.1 基本介绍

牛顿法属于利用一阶和二阶导数的无约束目标最优化方法。基本思想是，在每一次迭代中，以牛顿方向为搜索方向进行更新。牛顿法对目标的可导性更严格，要求二阶可导，有Hesse矩阵求逆的计算复杂的缺点。XGBoost本质上就是利用牛顿法进行优化的。

1.2 基本原理

现在推导牛顿法。
假设无约束最优化问题是

minxf(x)$$\underset{x}{min}f(x)$$

对于一维 x$x$ 的情况，可以将 f(x(t+1))$f(x^{(t+1)})$ 在 x(t)$x^{(t)}$ 附近用二阶泰勒展开近似：

f(x(t+1))=f(x(t))+f′(x(t))Δx+12f″(x(t))Δx2$$f(x^{(t+1)})=f(x^{(t)})+f^{\prime }(x^{(t)})\mathrm{\Delta }x+\frac{1}{2}{f}^{″}(x^{(t)})\mathrm{\Delta }{x}^2$$

然后用泰勒展开的极值点近似 f(x)$f(x)$ 的极值点：

∂f(x(t+1))∂x(t+1)=f′(x(t))+f″(x(t))Δx=0$$\frac{\mathrm{\partial }f(x^{(t+1)})}{\mathrm{\partial }{x}^{(t+1)}}=f^{\prime }(x^{(t)})+f^{″}(x^{(t)})\mathrm{\Delta }x=0$$

因此

Δx=x(t+1)−x(t)=−f′(x(t))f″(x(t))=−gtht$$\mathrm{\Delta }x=x^{(t+1)}-x^{(t)}=-\frac{f^{\prime }(x^{(t)})}{f^{″}(x^{(t)})}=-\frac{g_t}{h_t}$$

于是得到迭代公式，g$g$和h$h$分别是目标在当前x$x$上的一阶和二阶导

x(t+1)=x(t)−gtht$$x^{(t+1)}=x^{(t)}-\frac{g_t}{h_t}$$

推广到x$x$是多维向量的情况，gt$g_t$ 仍然是向量，而 Ht$H_t$ 是Hesse矩阵

H=[∂2f∂xi∂xj]$$H=\left[\frac{{\mathrm{\partial }}^{2}f}{\mathrm{\partial }{x}_i\mathrm{\partial }{x}_j}\right]$$

以二维x=(x1,x2)$x=(x_1,x_2)$为例：

H=⎡⎣⎢⎢⎢∂2f∂x21∂2f∂x2x1∂2f∂x1x2∂2f∂x22⎤⎦⎥⎥⎥$$H=\left[\begin{array}{cc}\frac{{\mathrm{\partial }}^{2}f}{\mathrm{\partial }{x}_1^2}& \frac{{\mathrm{\partial }}^{2}f}{\mathrm{\partial }{x}_1{x}_2}\\ \frac{{\mathrm{\partial }}^{2}f}{\mathrm{\partial }{x}_2{x}_1}& \frac{{\mathrm{\partial }}^{2}f}{\mathrm{\partial }{x}_2^2}\end{array}\right]$$

参数更新方程推广为：

x(t+1)=x(t)−H−1tgt$$x^{(t+1)}=x^{(t)}-H_t^{-1}{g}_t$$

可见，每一次迭代的更新方向都是当前点的牛顿方向，步长固定为1。每一次都需要计算一阶导数g$g$以及Hesse矩阵的逆矩阵，对于高维特征而言，求逆矩阵的计算量巨大且耗时。

1.3 阻尼牛顿法

从上面的推导中看出，牛顿方向 −H−1g$-H^{-1}g$ 能使得更新后函数处于极值点，但是它不一定是极小点，也就是说牛顿方向可能是下降方向，也可能是上升方向，以至于当初始点远离极小点时，牛顿法有可能不收敛。因此提出阻尼牛顿法，在牛顿法的基础上，每次迭代除了计算更新方向（牛顿方向），还要对最优步长做一维搜索。

算法步骤

（1）给定给初始点 x(0)$x^{(0)}$，允许误差 ϵ$ϵ$
（2）计算点 x(t)$x^{(t)}$ 处梯度gt$g_t$和Hesse矩阵H$H$，若|gt|<ϵ$|g_t|<ϵ$则停止迭代
（3）计算点 x(t)$x^{(t)}$ 处的牛顿方向作为搜索方向：

d(t)=−H−1tgt$$d^{(t)}=-H_t^{-1}{g}_t$$

（4）从点 x(t)$x^{(t)}$ 出发，沿着牛顿方向 d(t)$d^{(t)}$ 做一维搜索，获得最优步长：

λt=argminλf(x(t)+λ⋅d(t))$${\lambda }_t=\arg \underset{\lambda }{min}f(x^{(t)}+\lambda \cdot d^{(t)})$$

（5）更新参数

x(t+1)=x(t)+λt⋅d(t)$$x^{(t+1)}=x^{(t)}+{\lambda }_t\cdot d^{(t)}$$

---

二、拟牛顿法

2.1 提出的初衷

牛顿法中的Hesse矩阵H$H$在稠密时求逆计算量大，也有可能没有逆（Hesse矩阵非正定）。拟牛顿法提出，用不含二阶导数的矩阵 Ut$U_t$ 替代牛顿法中的 H−1t$H_t^{-1}$，然后沿搜索方向 −Utgt$-U_t{g}_t$ 做一维搜索。根据不同的 Ut$U_t$ 构造方法有不同的拟牛顿法。
注意拟牛顿法的 关键词：

• 不用算二阶导数
• 不用求逆

2.2 拟牛顿条件

牛顿法的搜索方向是

d(t)=−H−1tgt$$d^{(t)}=-H_t^{-1}{g}_t$$

为了不算二阶导及其逆矩阵，设法构造一个矩阵 U$U$，用它来逼近 H−1$H^{-1}$
现在为了方便推导，假设 f(x)$f(x)$ 是二次函数，于是 Hesse 矩阵 H$H$ 是常数阵，任意两点 x(t)$x^{(t)}$ 和 x(t+1)$x^{(t+1)}$ 处的梯度之差是：

▽f(x(t+1))−▽f(x(t))=H⋅(x(t+1)−x(t))$$▽f(x^{(t+1)})-▽f(x^{(t)})=H\cdot (x^{(t+1)}-x^{(t)})$$

等价于

x(t+1)−x(t)=H−1⋅[▽f(x(t+1))−▽f(x(t))]$$x^{(t+1)}-x^{(t)}=H^{-1}\cdot [▽f(x^{(t+1)})-▽f(x^{(t)})]$$

那么对非二次型的情况，也仿照这种形式，要求近似矩阵 U$U$ 满足类似的关系：

x(t+1)−x(t)=Ut+1⋅[▽f(x(t+1))−▽f(x(t))]$$x^{(t+1)}-x^{(t)}=U_{t+1}\cdot [▽f(x^{(t+1)})-▽f(x^{(t)})]$$

或者写成

Δxt=Ut+1⋅Δgt$$\mathrm{\Delta }{x}_t=U_{t+1}\cdot \mathrm{\Delta }{g}_t$$

以上就是拟牛顿条件，不同的拟牛顿法，区别就在于如何确定 U$U$。

2.3 DFP法

为了方便区分，下面把U$U$称作D$D$（表示DFP）。

DFP推导

现在已知拟牛顿条件

Δxt=Dt+1⋅Δgt$$\mathrm{\Delta }{x}_t=D_{t+1}\cdot \mathrm{\Delta }{g}_t$$

假设已知Dt$D_t$，希望用叠加的方式求Dt+1$D_{t+1}$，即 Dt+1=Dt+ΔDt$D_{t+1}=D_t+\mathrm{\Delta }{D}_t$，代入得到

ΔDtΔgt=Δxt−DtΔgt$$\mathrm{\Delta }{D}_t\mathrm{\Delta }{g}_t=\mathrm{\Delta }{x}_t-D_t\mathrm{\Delta }{g}_t$$

假设满足这个等式的 ΔDt$\mathrm{\Delta }{D}_t$ 是这样的形式：

ΔDt=Δxt⋅qTt−DtΔgt⋅wTt$$\mathrm{\Delta }{D}_t=\mathrm{\Delta }{x}_t\cdot q_t^T-D_t\mathrm{\Delta }{g}_t\cdot w_t^T$$

首先，对照一下就能发现：

qTt⋅Δgt=wTt⋅Δgt=In$$q_t^T\cdot \mathrm{\Delta }{g}_t=w_t^T\cdot \mathrm{\Delta }{g}_t=I_n$$

其次，要保证 ΔDt$\mathrm{\Delta }{D}_t$ 是对称的，参照 ΔDt$\mathrm{\Delta }{D}_t$ 的表达式，最简单就是令

qt=αtΔxtwt=βtDtΔgt$$\begin{gathered} q_t={\alpha }_t\mathrm{\Delta }{x}_t \\ {w}_t={\beta }_t{D}_t\mathrm{\Delta }{g}_t \end{gathered}$$

第二个条件代入第一个得到：

αt=1ΔgTtΔxtβt=1ΔgTtDtΔgt$$\begin{gathered} {\alpha }_t=\frac{1}{\mathrm{\Delta }{g}_t^T\mathrm{\Delta }{x}_t} \\ {\beta }_t=\frac{1}{\mathrm{\Delta }{g}_t^T{D}_t\mathrm{\Delta }{g}_t} \end{gathered}$$

然后代入回ΔDt$\mathrm{\Delta }{D}_t$ 的表达式：

ΔDt=ΔxtΔxTtΔgTtΔxt−DtΔgtΔgTtDtΔgTtDtΔgt$$\mathrm{\Delta }{D}_t=\frac{\mathrm{\Delta }{x}_t\mathrm{\Delta }{x}_t^T}{\mathrm{\Delta }{g}_t^T\mathrm{\Delta }{x}_t}-\frac{D_t\mathrm{\Delta }{g}_t\mathrm{\Delta }{g}_t^T{D}_t}{\mathrm{\Delta }{g}_t^T{D}_t\mathrm{\Delta }{g}_t}$$

观察一下两项分式，第一项仅涉及向量乘法，时间复杂度是O(n)$O(n)$，第二项涉及矩阵乘法，时间复杂度是O(n2)$O(n^2)$，综合起来是O(n2)$O(n^2)$。

DFP算法步骤

（1）给定初始点 x(0)$x^{(0)}$，允许误差 ϵ$ϵ$，令 D0=In$D_0=I_n$（n$n$是x$x$的维数），t=0$t=0$
（2）计算搜索方向 d(t)=−D−1t⋅gt$d^{(t)}=-D_t^{-1}\cdot g_t$
（3）从点 x(t)$x^{(t)}$ 出发，沿着 d(t)$d^{(t)}$ 做一维搜索，获得最优步长并更新参数：

λt=argminλf(x(t)+λ⋅d(t))x(t+1)=x(t)+λt⋅d(t)$$\begin{gathered} {\lambda }_t=\arg \underset{\lambda }{min}f(x^{(t)}+\lambda \cdot d^{(t)}) \\ {x}^{(t+1)}=x^{(t)}+{\lambda }_t\cdot d^{(t)} \end{gathered}$$

（4）判断精度，若|gt+1|<ϵ$|g_{t+1}|<ϵ$则停止迭代，否则转（5）
（5）计算Δg=gt+1−gt$\mathrm{\Delta }g=g_{t+1}-g_t$，Δx=x(t+1)−x(t)$\mathrm{\Delta }x=x^{(t+1)}-x^{(t)}$，更新 H$H$

Dt+1=Dt+ΔxΔxTΔgTΔx−DtΔgΔgTDtΔgTDtΔg$$D_{t+1}=D_t+\frac{\mathrm{\Delta }x\mathrm{\Delta }{x}^T}{\mathrm{\Delta }{g}^T\mathrm{\Delta }x}-\frac{D_t\mathrm{\Delta }g\mathrm{\Delta }{g}^T{D}_t}{\mathrm{\Delta }{g}^T{D}_t\mathrm{\Delta }g}$$

（6）t=t+1$t=t+1$，转（2）

2.4 BFGS法

为了方便区分，下面把U$U$称作B−1$B^{-1}$（表示BFGS）。

BFGS推导

拟牛顿条件

Δxt=B−1t+1⋅ΔgtΔgt=Bt+1⋅Δxt$$\begin{gathered} \mathrm{\Delta }{x}_t=B_{t+1}^{-1}\cdot \mathrm{\Delta }{g}_t \\ \mathrm{\Delta }{g}_t=B_{t+1}\cdot \mathrm{\Delta }{x}_t \end{gathered}$$

推导与DFP相似，但是，可以看到BFGS这种拟牛顿条件的形式与BFP的是对偶的，所以迭代公式只要把 Δxt$\mathrm{\Delta }{x}_t$ 和 Δgt$\mathrm{\Delta }{g}_t$ 调换一下就好。

ΔBt=ΔgtΔgTtΔxTtΔgt−BtΔxtΔxTtBtΔxTtBtΔxt$$\mathrm{\Delta }{B}_t=\frac{\mathrm{\Delta }{g}_t\mathrm{\Delta }{g}_t^T}{\mathrm{\Delta }{x}_t^T\mathrm{\Delta }{g}_t}-\frac{B_t\mathrm{\Delta }{x}_t\mathrm{\Delta }{x}_t^T{B}_t}{\mathrm{\Delta }{x}_t^T{B}_t\mathrm{\Delta }{x}_t}$$

只不过有个问题，按照下面这个迭代公式，不也一样要求逆吗？这就要引入谢尔曼莫里森公式了。

Δxt=B−1t+1⋅Δgt$$\mathrm{\Delta }{x}_t=B_{t+1}^{-1}\cdot \mathrm{\Delta }{g}_t$$

Sherman-Morrison 公式

对于任意非奇异方阵A$A$，u,v∈Rn$u,v\in R^n$是n$n$维向量，若1+vTA−1u≠0$1+v^T{A}^{-1}u\ne 0$，则

(A+uvT)−1=A−1−(A−1u)(vTA−1)1+vTA−1u$$(A+u{v}^T{)}^{-1}=A^{-1}-\frac{(A^{-1}u)(v^T{A}^{-1})}{1+v^T{A}^{-1}u}$$

该公式描述了在矩阵A$A$发生某种变化时，如何利用之前求好的逆，求新的逆。
对迭代公式引入两次 Sherman-Morrison 公式就能得到

B−1t+1=(In−ΔxtΔgTtΔxTtΔgt)B−1t(In−ΔgtΔxTtΔxTtΔgt)+ΔxtΔxTtΔxTtΔgt$$B_{t+1}^{-1}=\left(I_n-\frac{\mathrm{\Delta }{x}_t\mathrm{\Delta }{g}_t^T}{\mathrm{\Delta }{x}_t^T\mathrm{\Delta }{g}_t}\right)B_t^{-1}\left(I_n-\frac{\mathrm{\Delta }{g}_t\mathrm{\Delta }{x}_t^T}{\mathrm{\Delta }{x}_t^T\mathrm{\Delta }{g}_t}\right)+\frac{\mathrm{\Delta }{x}_t\mathrm{\Delta }{x}_t^T}{\mathrm{\Delta }{x}_t^T\mathrm{\Delta }{g}_t}$$

就得到了逆矩阵之间的推导。可能有人会问，第一个矩阵不也要求逆吗？其实这是一个迭代算法，初始矩阵设为单位矩阵（对角阵也可以）就不用求逆了。
这个公式的详细推导可以参考这里或者这里。

BFGS算法步骤

虽然下面的矩阵写成B−1$B^{-1}$，但要明确，BFGS从头到尾都不需要算逆，把下面的B−1$B^{-1}$换成H$H$这个符号，也是一样的。
（1）给定初始点 x(0)$x^{(0)}$，允许误差 ϵ$ϵ$，设置 B−10$B_0^{-1}$，t=0$t=0$
（2）计算搜索 d(t)=−B−1t⋅gt$d^{(t)}=-B_t^{-1}\cdot g_t$
（3）从点 x(t)$x^{(t)}$ 出发，沿着 d(t)$d^{(t)}$ 做一维搜索，获得最优步长并更新参数：

λt=argminλf(x(t)+λ⋅d(t))x(t+1)=x(t)+λt⋅d(t)$$\begin{gathered} {\lambda }_t=\arg \underset{\lambda }{min}f(x^{(t)}+\lambda \cdot d^{(t)}) \\ {x}^{(t+1)}=x^{(t)}+{\lambda }_t\cdot d^{(t)} \end{gathered}$$

（4）判断精度，若|gt+1|<ϵ$|g_{t+1}|<ϵ$则停止迭代，否则转（5）
（5）计算Δg=gt+1−gt$\mathrm{\Delta }g=g_{t+1}-g_t$，Δx=x(t+1)−x(t)$\mathrm{\Delta }x=x^{(t+1)}-x^{(t)}$，更新 B−1$B^{-1}$，然后

B−1t+1=(In−ΔxtΔgTtΔxTtΔgt)B−1t(In−ΔgtΔxTtΔxTtΔgt)+ΔxtΔxTtΔxTtΔgt$$B_{t+1}^{-1}=\left(I_n-\frac{\mathrm{\Delta }{x}_t\mathrm{\Delta }{g}_t^T}{\mathrm{\Delta }{x}_t^T\mathrm{\Delta }{g}_t}\right)B_t^{-1}\left(I_n-\frac{\mathrm{\Delta }{g}_t\mathrm{\Delta }{x}_t^T}{\mathrm{\Delta }{x}_t^T\mathrm{\Delta }{g}_t}\right)+\frac{\mathrm{\Delta }{x}_t\mathrm{\Delta }{x}_t^T}{\mathrm{\Delta }{x}_t^T\mathrm{\Delta }{g}_t}$$

（6）t=t+1$t=t+1$，转（2）

2.5 L-BFGS法（Limited-memory BFGS）

对于d$d$维参数，BFGS算法需要保存一个O(d2)$O(d^2)$大小的B−1$B^{-1}$矩阵，实际上只需要每一轮的Δx$\mathrm{\Delta }x$和Δg$\mathrm{\Delta }g$，也可以递归计算出当前迭代的B−1$B^{-1}$矩阵，L-BFGS就是基于这种思想，实现了节省内存的BFGS。

L-BFGS推导

BFGS的递推公式：

B−1t+1=(In−ΔxtΔgTtΔxTtΔgt)B−1t(In−ΔgtΔxTtΔxTtΔgt)+ΔxtΔxTtΔxTtΔgt$$B_{t+1}^{-1}=\left(I_n-\frac{\mathrm{\Delta }{x}_t\mathrm{\Delta }{g}_t^T}{\mathrm{\Delta }{x}_t^T\mathrm{\Delta }{g}_t}\right)B_t^{-1}\left(I_n-\frac{\mathrm{\Delta }{g}_t\mathrm{\Delta }{x}_t^T}{\mathrm{\Delta }{x}_t^T\mathrm{\Delta }{g}_t}\right)+\frac{\mathrm{\Delta }{x}_t\mathrm{\Delta }{x}_t^T}{\mathrm{\Delta }{x}_t^T\mathrm{\Delta }{g}_t}$$

现在假设ρt=1ΔxTtΔgt${\rho }_t=\frac{1}{\mathrm{\Delta }{x}_t^T\mathrm{\Delta }{g}_t}$，Vt=In−ρtΔgtΔxTt$V_t=I_n-{\rho }_t\mathrm{\Delta }{g}_t\mathrm{\Delta }{x}_t^T$，则递推公式可以写成

B−1t+1=VTtB−1tVt+ρtΔxtΔxTt$$B_{t+1}^{-1}=V_t^T{B}_t^{-1}{V}_t+{\rho }_t\mathrm{\Delta }{x}_t\mathrm{\Delta }{x}_t^T$$

给定的初始矩阵B−10$B_0^{-1}$后，之后的每一轮都可以递推计算

B−11=VT0B−10V0+ρ0Δx0ΔxT0B−12=VT1B−10V1+ρ1Δx1ΔxT1=(VT1VT0)B−10(V0V1)+VT1ρ0Δx0ΔxT0V1+ρ1Δx1ΔxT1$$\begin{gathered} B_1^{-1}=V_0^T{B}_0^{-1}{V}_0+{\rho }_0\mathrm{\Delta }{x}_0\mathrm{\Delta }{x}_0^T \\ {B}_2^{-1}=V_1^T{B}_0^{-1}{V}_1+{\rho }_1\mathrm{\Delta }{x}_1\mathrm{\Delta }{x}_1^T \\ =(V_1^T{V}_0^T)B_0^{-1}(V_0{V}_1)+V_1^T{\rho }_0\mathrm{\Delta }{x}_0\mathrm{\Delta }{x}_0^T{V}_1+{\rho }_1\mathrm{\Delta }{x}_1\mathrm{\Delta }{x}_1^T \end{gathered}$$

一直到最后B−1k+1$B_{k+1}^{-1}$可以由t=0$t=0$到t=k$t=k$的Δxt$\mathrm{\Delta }{x}_t$和Δgt$\mathrm{\Delta }{g}_t$表示：

B−1t+1=++++(VTtVTt−1⋯VT1VT0)B−10(V0⋯Vt−1Vt)(VTtVTt−1⋯VT2VT1)(ρ0Δx0ΔxT0)(V1⋯Vt−1Vt)⋯VTt(ρt−1Δxt−1ΔxTt−1)VtρtΔxtΔxTt$$\begin{array}{cccc}{B}_{t+1}^{-1}& =& & (V_t^T{V}_{t-1}^T\cdots V_1^T{V}_0^T)B_0^{-1}(V_0\cdots V_{t-1}{V}_t)\\ & & +& (V_t^T{V}_{t-1}^T\cdots V_2^T{V}_1^T)({\rho }_0\mathrm{\Delta }{x}_0\mathrm{\Delta }{x}_0^T)(V_1\cdots V_{t-1}{V}_t)\\ & & +& \cdots \\ & & +& V_t^T({\rho }_{t-1}\mathrm{\Delta }{x}_{t-1}\mathrm{\Delta }{x}_{t-1}^T)V_t\\ & & +& {\rho }_t\mathrm{\Delta }{x}_t\mathrm{\Delta }{x}_t^T\end{array}$$

看起来很长，其实可以写成一个求和项

B−1t+1=(∏i=t0VTi)B−10(∏i=0tVi)+∑j=0t(∏i=tj+1VTi)(ρjΔxjΔxTj)(∏i=j+1tVi)$$B_{t+1}^{-1}=\left(\prod _{i=t}^0{V}_i^T\right)B_0^{-1}\left(\prod _{i=0}^t{V}_i\right)+\sum _{j=0}^t\left(\prod _{i=t}^{j+1}{V}_i^T\right)\left({\rho }_j\mathrm{\Delta }{x}_j\mathrm{\Delta }{x}_j^T\right)\left(\prod _{i=j+1}^t{V}_i\right)$$

这个求和项包含了从0$0$到t$t$的所有Δx$\mathrm{\Delta }x$和Δg$\mathrm{\Delta }g$，而根据实际需要，可以只取最近的m$m$个，也就是：

B−1t=(∏i=t−1t−mVTi)B−10(∏i=t−mt−1Vi)+∑j=t−1t−m(∏i=tj+1VTi)(ρjΔxjΔxTj)(∏i=j+1tVi)$$B_t^{-1}=\left(\prod _{i=t-1}^{t-m}{V}_i^T\right)B_0^{-1}\left(\prod _{i=t-m}^{t-1}{V}_i\right)+\sum _{j=t-1}^{t-m}\left(\prod _{i=t}^{j+1}{V}_i^T\right)\left({\rho }_j\mathrm{\Delta }{x}_j\mathrm{\Delta }{x}_j^T\right)\left(\prod _{i=j+1}^t{V}_i\right)$$

工程上的L-BFGS

我们关心的其实不是B−1t$B_t^{-1}$本身如何，算B−1t$B_t^{-1}$的根本目的是要算本轮搜索方向 B−1tgt$B_t^{-1}{g}_t$
以下算法摘自《Numerical Optimization》，它可以高效地计算出拟牛顿法每一轮的搜索方向。仔细观察一下，你会发现它实际上就是复现上面推导的那一堆很长的递推公式，你所需要的是最近m$m$轮的Δx$\mathrm{\Delta }x$和Δg$\mathrm{\Delta }g$，后向和前向算完得到最终的 r$r$ 就是搜索方向 B−1tgt$B_t^{-1}{g}_t$，之后要做一维搜索或者什么的都可以。
解释一下算法的符号和本文符号之间的对应关系，si=Δxi$s_i=\mathrm{\Delta }{x}_i$，yi=Δgi$y_i=\mathrm{\Delta }{g}_i$，Hk=B−1k$H_k=B_k^{-1}$
代码实现可以参考这里。

L-BFGS算法步骤

（1）给定初始点 x(0)$x^{(0)}$，允许误差 ϵ$ϵ$，预定保留最近m$m$个向量，设置 B−10$B_0^{-1}$，t=0$t=0$
（2）用Algorithm 9.1计算搜索方向 d(t)=−B−1t⋅gt$d^{(t)}=-B_t^{-1}\cdot g_t$
（3）从点 x(t)$x^{(t)}$ 出发，沿着 d(t)$d^{(t)}$ 做一维搜索，获得最优步长并更新参数：

λt=argminλf(x(t)+λ⋅d(t))x(t+1)=x(t)+λt⋅d(t)$$\begin{gathered} {\lambda }_t=\arg \underset{\lambda }{min}f(x^{(t)}+\lambda \cdot d^{(t)}) \\ {x}^{(t+1)}=x^{(t)}+{\lambda }_t\cdot d^{(t)} \end{gathered}$$

（4）判断精度，若|gt+1|<ϵ$|g_{t+1}|<ϵ$则停止迭代，否则转（5）
（5）判断 t>m$t>m$，删掉存储的 Δxt−m$\mathrm{\Delta }{x}_{t-m}$ 和 Δgt−m$\mathrm{\Delta }{g}_{t-m}$
（5）计算Δg=gt+1−gt$\mathrm{\Delta }g=g_{t+1}-g_t$，Δx=x(t+1)−x(t)$\mathrm{\Delta }x=x^{(t+1)}-x^{(t)}$，令t=t+1$t=t+1$，转（2）

---

最后，有时候你看不懂BFGS到底意味着什么，并不是你英文差，而是因为这个简称真的没有意义。。。。。

---

参考资料

1. 【博客】LBFGS方法推导-慢慢的回味
2. 【博客】数值优化：理解L-BFGS算法
3. 【博客】无约束优化算法——牛顿法与拟牛顿法（DFP，BFGS，LBFGS）
4. 【博客】无约束最优化方法——牛顿法、拟牛顿法、BFGS、LBFGS
5. 【博客】Numerical Optimization: Understanding L-BFGS
6. 【论文】A Stochastic Quasi-Newton Method for Online Convex Optimization
7. 【书籍】Numeric Optimization