L-BFGS 算法

Reference

• 牛顿法与拟牛顿法学习笔记（五）L-BFGS 算法
• Limited-memory BFGS
• 一文读懂L-BFGS算法
• L-BFGS算法介绍
• 【技术分享】L-BFGS算法

【概述】

在 Sherman-Morrison 公式下的 BFGS 算法中，需要用到一个 $N\times N$ 的矩阵 $G_k$，当 $N$ 很大时，存储这个矩阵将十分消耗计算机的资源

以 $N=100,000$ 为例，其所需要消耗的内存空间如下：

$$\frac{N\:阶矩阵字节数\quad\quad}{1GB\:的字节数\quad\quad}=\frac{10^5 \cdot 10^5 \cdot 8}{2^{10} \cdot 2^{10} \cdot 2^{10}}=74.5\:GB$$

这对一般的服务器是难以承受的，虽然考虑到矩阵 $G_k$ 是对称阵，使用对称存储的方法内存可以降低一半，但是 $10W$ 规模的数据在机器学习问题中只能算是中小规模

为减少 BFGS 算法迭代过程中的内存开销，有了 L-BFGS 算法（Limited-memory BFGS），其对矩阵 $G_k$ 进行了近似，不再存储完整的矩阵 $G_k$，而是存储计算过程中的向量序列 $\{\boldsymbol{\delta_i}\}$ 和 $\{\mathbf{y_i}\}$，在需要矩阵 $G_k$ 时，利用这两个向量序列计算来代替

同时，向量序列 $\{\boldsymbol{\delta_i}\}$ 和 $\{\mathbf{y_i}\}$，也并非每次都存储，而是根据用户机器的内存每次各自存储最新的 $m$ 个，每次计算 $G_k$ 时，都利用这最新的 $m$ 个向量序列进行计算，这样存储就由原来的 $O(N^2)$ 降低到了 $O(mN)$

【算法原理】

L-BFGS 算法的出发点是 BFGS 算法中的迭代式：

$$G_{k+1} = (I-\frac{\boldsymbol{\delta_k}\mathbf{y_k}^T}{\mathbf{y_k}^T\boldsymbol{\delta_k}})G_k(I-\frac{\mathbf{y_k}\boldsymbol{\delta_k}^T}{\mathbf{y_k}^T\boldsymbol{\delta_k}})+\frac{\boldsymbol{\delta_k}\boldsymbol{\delta_k}^T}{\mathbf{y_k}^T\boldsymbol{\delta_k}}$$

记：$\rho_k=\frac{1}{\mathbf{y_k}^T\boldsymbol{\delta_k}}$，$V_k=I-\rho_k\mathbf{y_k}\boldsymbol{\delta_k}^T$，则上式可写为：

$$G_{k+1}=V_k^TG_kV_k+\rho_k\boldsymbol{\delta_k}\boldsymbol{\delta_k}^T$$

若给定初始矩阵 $G_0=I$，则依次可得：

$$\begin{align} G_1 &= V_0^TG_0V_0+\rho_0\boldsymbol{\delta_0}\boldsymbol{\delta_0}^T \notag \\ G_2 &= V_1^TG_1V_1+\rho_1\boldsymbol{\delta_1}\boldsymbol{\delta_1}^T \notag \\ &= V_1^T(V_0^TG_0V_0+\rho_0\boldsymbol{\delta_0}\boldsymbol{\delta_0}^T)V_1+\rho_1\boldsymbol{\delta_1}\boldsymbol{\delta_1}^T \notag \\ &= V_1^TV_0^TG_0V_0V_1+V_1^T\rho_0\boldsymbol{\delta_0}\boldsymbol{\delta_0}^TV_1+\rho_1\boldsymbol{\delta_1}\boldsymbol{\delta_1}^T \notag \\ G_3 &= V_2^TG_2V_2+\rho_2\boldsymbol{\delta_2}\boldsymbol{\delta_2}^T \notag \\ &= V_2^T(V_1^TV_0^TG_0V_0V_1+V_1^T\rho_0\boldsymbol{\delta_0}\boldsymbol{\delta_0}^TV_1+\rho_1\boldsymbol{\delta_1}\boldsymbol{\delta_1}^T)V_2+\rho_2\boldsymbol{\delta_2}\boldsymbol{\delta_2}^T \notag \\ &= V_2^TV_1^TV_0^TG_0V_0V_1V_2+V_2^TV_1^T\rho_0\boldsymbol{\delta_0}\boldsymbol{\delta_0}^TV_1V_2 + V_2^T\rho_1\boldsymbol{\delta_1}\boldsymbol{\delta_1}^TV_2 \notag \\ &\quad + \rho_2\boldsymbol{\delta_2}\boldsymbol{\delta_2}^T \notag \\ ... & \notag \end{align}$$

由此，可进行递推，有：

$$\begin{align} G_{k+1} &= \quad (V_k^TV_{k-1}^T...V_1^TV_0^T) G_0 (V_0V_1...V_{k-1}V_k) \notag \\ &\quad + (V_k^TV_{k-1}^T...V_2^TV_1^T)(\rho_0\boldsymbol{\delta_0}\boldsymbol{\delta_0}^T)(V_1V_2...V_{k-1}V_k) \notag \\ &\quad + (V_k^TV_{k-1}^T...V_3^TV_2^T)(\rho_1\boldsymbol{\delta_1}\boldsymbol{\delta_1}^T)(V_2V_3...V_{k-1}V_k) \notag \\ &\quad + \dots \notag \\ &\quad + (V_k^TV_{k-1}^T)(\rho_{k-2}\boldsymbol{\delta_{k-2}}\boldsymbol{\delta_{k-2}}^T)(V_{k-1}V_k) \notag \\ &\quad + V_k^T(\rho_{k-1}\boldsymbol{\delta_{k-1}}\boldsymbol{\delta_{k-1}}^T)V_k \notag \\ &\quad + \rho_k\boldsymbol{\delta_k}\boldsymbol{\delta_k}^T \notag \end{align}$$

可见，计算 $G_{k+1}$ 需要用到向量序列 $\{(\boldsymbol{\delta_i},\mathbf{y_i})\}_{i=0}^k$，若从 $\boldsymbol{\delta_0},\mathbf{y_0}$ 开始连续地存储 $m$ 组的话，只能依次计算到 $G_m$

也就是说，如果想要求 $G_{m+1},G_{m+2},…$ 的话，就要考虑丢弃一些最早生成的向量

举例来说，如果要计算 $G_{m+1}$，就保存 $\{(\boldsymbol{\delta_i},\mathbf{y_i})\}_{i=1}^{m}$，丢弃 $\{(\boldsymbol{\delta_0},\mathbf{y_0})\}$，如果要计算 $G_{m+2}$，就保存 $\{(\boldsymbol{\delta_i},\mathbf{y_i})\}_{i=2}^{m+1}$，丢弃 $\{(\boldsymbol{\delta_i},\mathbf{y_i})\}_{i=0}^{1}$

在舍弃一些向量后，就只能近似计算了，当 $k+1>m$ 时，按照上述的 $G_{k+1}$ 可以构造近似计算公式，即：

$$\begin{align} G_{k+1} &\approx \quad (V_k^TV_{k-1}^T...V_{k-m+2}^TV_{k-m+1}^T) G_0 (V_{k-m+1}V_{k-m+2}...V_{k-1}V_k) \notag \\ &\quad + (V_k^TV_{k-1}^T...V_{k-m+3}^TV_{k-m+2}^T)(\rho_0\boldsymbol{\delta_0}\boldsymbol{\delta_0}^T)(V_{k-m+2}V_{k-m+3}...V_{k-1}V_k) \notag \\ &\quad + (V_k^TV_{k-1}^T...V_{k-m+4}^TV_{k-m+3}^T)(\rho_1\boldsymbol{\delta_1}\boldsymbol{\delta_1}^T)(V_{k-m+3}V_{k-m+4}...V_{k-1}V_k) \notag \\ &\quad + \dots \notag \\ &\quad + (V_k^TV_{k-1}^T)(\rho_{k-2}\boldsymbol{\delta_{k-2}}\boldsymbol{\delta_{k-2}}^T)(V_{k-1}V_k) \notag \\ &\quad + V_k^T(\rho_{k-1}\boldsymbol{\delta_{k-1}}\boldsymbol{\delta_{k-1}}^T)V_k \notag \\ &\quad + \rho_k\boldsymbol{\delta_k}\boldsymbol{\delta_k}^T \notag \end{align}$$

若引入 $\hat{m}=\min\{k,m-1\}$ 则可以将上述的递推式与近似式进行合并，即：

$$\begin{align} G_{k+1} &= \quad (V_k^TV_{k-1}^T...V_{k-\hat{m}+1}^TV_{k-\hat{m}}^T) G_0 (V_{k-\hat{m}}V_{k-\hat{m}+1}...V_{k-1}V_k) \notag \\ &\quad + (V_k^TV_{k-1}^T...V_{k-\hat{m}+2}^TV_{k-\hat{m}+1}^T)(\rho_0\boldsymbol{\delta_0}\boldsymbol{\delta_0}^T)(V_{k-\hat{m}+1}V_{k-m+2}...V_{k-1}V_k) \notag \\ &\quad + (V_k^TV_{k-1}^T...V_{k-\hat{m}+3}^TV_{k-\hat{m}+2}^T)(\rho_1\boldsymbol{\delta_1}\boldsymbol{\delta_1}^T)(V_{k-\hat{m}+2}V_{k-\hat{m}+3}...V_{k-1}V_k) \notag \\ &\quad + \dots \notag \\ &\quad + (V_k^TV_{k-1}^T)(\rho_{k-2}\boldsymbol{\delta_{k-2}}\boldsymbol{\delta_{k-2}}^T)(V_{k-1}V_k) \notag \\ &\quad + V_k^T(\rho_{k-1}\boldsymbol{\delta_{k-1}}\boldsymbol{\delta_{k-1}}^T)V_k \notag \\ &\quad + \rho_k\boldsymbol{\delta_k}\boldsymbol{\delta_k}^T \notag \end{align}$$

事实上，根据 BFGS 算法流程可知，$G_k$ 的作用仅用于计算 $G_k\mathbf{g_k}$ 来获取搜索方向，因此，若能根据上式设计出一种能够快速计算 $G_k\mathbf{g_k}$ 的算法即可

【算法流程】

快速计算 $G_k\mathbf{g_k}$ 的算法流程如下：

Step 1：初始化

$$\mathbf{s}=\left\{\begin{array}{rl} 0 ,& k\leq m \\ k-m,& k>m \end{array} \right. ,\quad L=\left\{\begin{array}{rl} k ,& k\leq m \\ m,& k>m \end{array} \right. ,\quad \mathbf{q_L}=\mathbf{g_k}$$

Step 2：后向循环

$$\begin{align} FOR \quad i&=L-1,L-2,...,1,0 \quad \{ & \notag \\ j &=i+\mathbf{s}; &&& \notag \\ \alpha_i &=\rho_j\boldsymbol{\delta_j}^T\mathbf{q_{i+1}}; && \notag \\ \mathbf{q_i} &=\mathbf{q_{i+1}}-\alpha_i\mathbf{y_j}; \notag \\ \}\quad\quad\quad \notag \end{align}$$

Step 3：前向循环

$$\begin{align} \mathbf{r_0}=G_0\mathbf{q_0} \notag \\ FOR \quad i&=0,1,...,L-2,L-1 \quad \{ & \notag \\ j &=i+\mathbf{s}; &&& \notag \\ \beta_j &=\rho_j\mathbf{y_j}^T\mathbf{r_{i+1}}; && \notag \\ \mathbf{r_{i+1}} &=\mathbf{r_i}+(\alpha_i - \beta_i) \boldsymbol{\delta_j}; \notag \\ \}\quad\quad\quad\quad \notag \end{align}$$

最后求出的 $\mathbf{r_L}$ 即为 $H_k\mathbf{g_k}$ 的值