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

拟牛顿法

一、牛顿法

1.1 基本介绍

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

1.2 基本原理

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

$$\min_x f(x)$$
对于一维 $x$ 的情况，可以将 $f(x^{(t+1)})$ 在 $x^{(t)}$ 附近用二阶泰勒展开近似：
$$f(x^{(t+1)})=f(x^{(t)})+f'(x^{(t)})\Delta x+\frac{1}{2}f''(x^{(t)})\Delta x^2$$
然后用泰勒展开的极值点近似 $f(x)$ 的极值点：
$$\frac{\partial f(x^{(t+1)})}{\partial x^{(t+1)}}=f'(x^{(t)})+f''(x^{(t)})\Delta x=0$$
因此
$$\Delta x = x^{(t+1)}-x^{(t)}=-\frac{f'(x^{(t)})}{f''(x^{(t)})}=-\frac{g_t}{h_t}$$
于是得到迭代公式， $g$和$h$分别是目标在当前 $x$上的一阶和二阶导
$$x^{(t+1)}=x^{(t)}-\frac{g_t}{h_t}$$
推广到 $x$是多维向量的情况，$g_t$ 仍然是向量，而 $H_t$ 是Hesse矩阵
$$H=\left [ \frac{\partial^2f}{\partial x_i\partial x_j} \right ]$$
以二维 $x=(x_1,x_2)$为例：
$$H=\left [\begin{matrix} \frac{\partial^2 f}{\partial x_1^2} & \frac{\partial^2 f}{\partial x_1x_2} \\ \frac{\partial^2 f}{\partial x_2x_1} & \frac{\partial^2 f}{\partial x_2^2} \end{matrix} \right ]$$
参数更新方程推广为：
$$x^{(t+1)}=x^{(t)}-H_t^{-1}g_t$$
可见，每一次迭代的更新方向都是当前点的牛顿方向，步长固定为1。每一次都需要计算一阶导数 $g$以及Hesse矩阵的逆矩阵，对于高维特征而言，求逆矩阵的计算量巨大且耗时。

1.3 阻尼牛顿法

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

算法步骤

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

$$d^{(t)}=-H_t^{-1}g_t$$
（4）从点 $x^{(t)}$ 出发，沿着牛顿方向 $d^{(t)}$ 做一维搜索，获得最优步长： $$\lambda_t = \arg \min_{\lambda} f(x^{(t)}+\lambda\cdot d^{(t)})$$
（5）更新参数 $$x^{(t+1)}=x^{(t)}+\lambda_t\cdot d^{(t)}$$

---

二、拟牛顿法

2.1 提出的初衷

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

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

2.2 拟牛顿条件

牛顿法的搜索方向是

$$d^{(t)}=-H_t^{-1}g_t$$
为了不算二阶导及其逆矩阵，设法构造一个矩阵 $U$，用它来逼近 $H^{-1}$
现在为了方便推导，假设 $f(x)$ 是二次函数，于是 Hesse 矩阵 $H$ 是常数阵，任意两点 $x^{(t)}$ 和 $x^{(t+1)}$ 处的梯度之差是：
$$\bigtriangledown f(x^{(t+1)}) - \bigtriangledown f(x^{(t)}) = H\cdot (x^{(t+1)}-x^{(t)})$$
等价于
$$x^{(t+1)}-x^{(t)} = H^{-1}\cdot [\bigtriangledown f(x^{(t+1)}) - \bigtriangledown f(x^{(t)})]$$
那么对非二次型的情况，也仿照这种形式，要求近似矩阵 $U$ 满足类似的关系：
$$x^{(t+1)}-x^{(t)}=U_{t+1}\cdot [\bigtriangledown f(x^{(t+1)})-\bigtriangledown f(x^{(t)})]$$
或者写成
$$\Delta x_t=U_{t+1}\cdot \Delta g_t$$
以上就是 拟牛顿条件，不同的拟牛顿法，区别就在于如何确定 $U$。

2.3 DFP法

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

DFP推导

现在已知拟牛顿条件

$$\Delta x_t=D_{t+1}\cdot \Delta g_t$$
假设已知 $D_t$，希望用叠加的方式求 $D_{t+1}$，即 $D_{t+1}=D_{t}+\Delta D_t$，代入得到
$$\Delta D_t \Delta g_t=\Delta x_t - D_t \Delta g_t$$
假设满足这个等式的 $\Delta D_t$ 是这样的形式：
$$\Delta D_t=\Delta x_t \cdot q_t^T-D_t\Delta g_t\cdot w_t^T$$
首先，对照一下就能发现：
$$q_t^T\cdot \Delta g_t=w_t^T \cdot \Delta g_t = I_n$$
其次，要保证 $\Delta D_t$ 是对称的，参照 $\Delta D_t$ 的表达式，最简单就是令
$$q_t=\alpha_t \Delta x_t\\ w_t=\beta_t D_t\Delta g_t$$
第二个条件代入第一个得到： $$\alpha_t=\frac{1}{\Delta g_t^T\Delta x_t} \\\beta_t=\frac{1}{\Delta g_t^TD_t\Delta g_t}$$
然后代入回 $\Delta D_t$ 的表达式：
$$\Delta D_t = \frac{\Delta x_t\Delta x_t^T}{\Delta g_t^T\Delta x_t}-\frac{D_t\Delta g_t\Delta g_t^TD_t}{\Delta g_t^TD_t\Delta g_t}$$
观察一下两项分式，第一项仅涉及向量乘法，时间复杂度是 $O(n)$，第二项涉及矩阵乘法，时间复杂度是 $O(n^2)$，综合起来是 $O(n^2)$。

DFP算法步骤

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

$$\lambda_t = \arg \min_{\lambda} f(x^{(t)}+\lambda\cdot d^{(t)})\\x^{(t+1)}=x^{(t)}+\lambda_t\cdot d^{(t)}$$
（4）判断精度，若 $|g_{t+1}|<\epsilon$则停止迭代，否则转（5）
（5）计算 $\Delta g=g_{t+1}-g_t$， $\Delta x=x^{(t+1)}-x^{(t)}$，更新 $H$
$$D_{t+1}=D_{t}+\frac{\Delta x\Delta x^T}{\Delta g^T\Delta x}-\frac{D_t\Delta g\Delta g^TD_t}{\Delta g^TD_t\Delta g}$$
（6） $t=t+1$，转（2）

2.4 BFGS法

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

BFGS推导

拟牛顿条件

$$\Delta x_t=B_{t+1}^{-1}\cdot \Delta g_t\\ \Delta g_t=B_{t+1} \cdot \Delta x_t$$
推导与DFP相似，但是，可以看到BFGS这种拟牛顿条件的形式与BFP的是对偶的，所以迭代公式只要把 $\Delta x_t$ 和 $\Delta g_t$ 调换一下就好。
$$\Delta B_t = \frac{\Delta g_t\Delta g_t^T}{\Delta x_t^T\Delta g_t}-\frac{B_t\Delta x_t\Delta x_t^TB_t}{\Delta x_t^TB_t\Delta x_t}$$
只不过有个问题，按照下面这个迭代公式，不也一样要求逆吗？这就要引入谢尔曼莫里森公式了。
$$\Delta x_t=B_{t+1}^{-1}\cdot \Delta g_t$$

Sherman-Morrison 公式

对于任意非奇异方阵$A$，$u,v\in R^n$是$n$维向量，若$1+v^TA^{-1}u\neq 0$，则

$$(A+uv^T)^{-1} = A^{-1}-\frac{(A^{-1}u)(v^TA^{-1})}{1+v^TA^{-1}u}$$
该公式描述了在矩阵 $A$发生某种变化时，如何利用之前求好的逆，求新的逆。
对迭代公式引入两次 Sherman-Morrison 公式就能得到
$$B^{-1}_{t+1}=\left (I_n-\frac{\Delta x_t \Delta g_t^T}{\Delta x_t^T \Delta g_t}\right )B_{t}^{-1}\left (I_n-\frac{\Delta g_t \Delta x_t^T}{\Delta x_t^T \Delta g_t}\right )+\frac{\Delta x_t \Delta x_t^T}{\Delta x_t^T \Delta g_t}$$
就得到了逆矩阵之间的推导。可能有人会问，第一个矩阵不也要求逆吗？其实这是一个迭代算法，初始矩阵设为单位矩阵（对角阵也可以）就不用求逆了。
这个公式的详细推导可以参考 这里或者 这里。

BFGS算法步骤

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

$$\lambda_t = \arg \min_{\lambda} f(x^{(t)}+\lambda\cdot d^{(t)})\\x^{(t+1)}=x^{(t)}+\lambda_t\cdot d^{(t)}$$
（4）判断精度，若 $|g_{t+1}|<\epsilon$则停止迭代，否则转（5）
（5）计算 $\Delta g=g_{t+1}-g_t$， $\Delta x=x^{(t+1)}-x^{(t)}$，更新 $B^{-1}$，然后
$$B^{-1}_{t+1}=\left (I_n-\frac{\Delta x_t \Delta g_t^T}{\Delta x_t^T \Delta g_t}\right )B_{t}^{-1}\left (I_n-\frac{\Delta g_t \Delta x_t^T}{\Delta x_t^T \Delta g_t}\right )+\frac{\Delta x_t \Delta x_t^T}{\Delta x_t^T \Delta g_t}$$
（6） $t=t+1$，转（2）

2.5 L-BFGS法（Limited-memory BFGS）

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

L-BFGS推导

BFGS的递推公式：

$$B^{-1}_{t+1}=\left (I_n-\frac{\Delta x_t \Delta g_t^T}{\Delta x_t^T \Delta g_t}\right )B_{t}^{-1}\left (I_n-\frac{\Delta g_t \Delta x_t^T}{\Delta x_t^T \Delta g_t}\right )+\frac{\Delta x_t\Delta x_t^T}{\Delta x_t^T \Delta g_t}$$
现在假设 $\rho_t = \frac{1}{\Delta x_t^T \Delta g_t}$， $V_t = I_n-\rho_t \Delta g_t \Delta x_t^T$，则递推公式可以写成
$$B^{-1}_{t+1}=V_t^TB^{-1}_{t}V_t+\rho_t \Delta x_t \Delta x_t^T$$
给定的初始矩阵 $B^{-1}_{0}$后，之后的每一轮都可以递推计算
$$B^{-1}_{1}=V_0^TB^{-1}_{0}V_0+\rho_0 \Delta x_0 \Delta x_0^T\\B^{-1}_{2}=V_1^TB^{-1}_{0}V_1+\rho_1 \Delta x_1 \Delta x_1^T\\ =(V_1^TV_0^T)B^{-1}_{0}(V_0V_1)+V_1^T\rho_0\Delta x_0 \Delta x_0^TV_1+\rho_1 \Delta x_1 \Delta x_1^T$$
一直到最后 $B^{-1}_{k+1}$可以由 $t=0$到 $t=k$的 $\Delta x_t$和 $\Delta g_t$表示：
$$\begin{matrix} B^{-1}_{t+1} & = & & (V_t^TV_{t-1}^T\cdots V_1^TV_0^T)B^{-1}_0 (V_0\cdots V_{t-1}V_t) \\ & & + & (V_t^TV_{t-1}^T\cdots V_2^TV_1^T)(\rho_0\Delta x_0\Delta x_0^T)(V_1\cdots V_{t-1}V_t)\\ & & + & \cdots\\ & & + & V_t^T(\rho_{t-1} \Delta x_{t-1} \Delta x_{t-1}^T)V_t\\ & & + & \rho_t \Delta x_t \Delta x_t^T \end{matrix}$$
看起来很长，其实可以写成一个求和项
$$B^{-1}_{t+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\Delta x_j \Delta x_j^T\right ) \left (\prod_{i=j+1}^t V_i \right )$$
这个求和项包含了从 $0$到$t$的所有 $\Delta x$和 $\Delta g$，而根据实际需要，可以只取最近的 $m$个，也就是：
$$B^{-1}_{t} = \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\Delta x_j \Delta x_j^T\right ) \left (\prod_{i=j+1}^t V_i \right )$$

工程上的L-BFGS

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

L-BFGS算法步骤

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

$$\lambda_t = \arg \min_{\lambda} f(x^{(t)}+\lambda\cdot d^{(t)})\\x^{(t+1)}=x^{(t)}+\lambda_t\cdot d^{(t)}$$
（4）判断精度，若 $|g_{t+1}|<\epsilon$则停止迭代，否则转（5）
（5）判断 $t>m$，删掉存储的 $\Delta x_{t-m}$ 和 $\Delta g_{t-m}$
（5）计算 $\Delta g=g_{t+1}-g_t$， $\Delta x=x^{(t+1)}-x^{(t)}$，令 $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