拟牛顿法之L-BFGS算法

最新推荐文章于 2023-01-09 08:19:12 发布

大鹏的NLP博客

最新推荐文章于 2023-01-09 08:19:12 发布

阅读量1.8k

点赞数

分类专栏：数值计算文章标签：优化 l-bfgs

本文链接：https://blog.csdn.net/philosophyatmath/article/details/70343129

版权

数值计算专栏收录该内容

6 篇文章 0 订阅

订阅专栏

关于优化算法的求解，比如有梯度下降法，坐标下降法，牛顿法和拟牛顿法。
梯度下降法是基于目标函数梯度的，算法的收敛速度是线性的，并且当问题是病态时或者问题规模较大时，收敛速度尤其慢（几乎不适用）；
坐标下降法虽然不用计算目标函数的梯度，但是其收敛速度依然很慢，因此它的适用范围也有局限；
牛顿法是基于目标函数的二阶导数（海森矩阵）的，其收敛速度较快，迭代次数较少，尤其是在最优值附近时，收敛速度是二次的。但牛顿法的问题在于当海森矩阵稠密时，每次迭代的计算量比较大，因为每次都会计算目标函数的海森矩阵的逆，这样一来，当问题规模较大时，不仅计算量大（有时大到不可计算），而且需要的存储空间也多，因此牛顿法在面对海量数据时由于每一步迭代的开销巨大而变得不适用；
拟牛顿法是在牛顿法的基础上引入了海森矩阵的近似矩阵，避免每次迭代都要计算海森矩阵的逆，拟牛顿法的收敛速度介于梯度下降法和牛顿法之间，是超线性的。拟牛顿法的问题也是当问题规模很大时，近似矩阵变得很稠密，在计算和存储上也有很大的开销，因此变得不实用。
另外需要注意的是，牛顿法在每次迭代时不能总是保证海森矩阵是正定的，一旦海森矩阵不是正定的，优化方向就会“跑偏”，从而使得牛顿法失效，也说明了牛顿法的鲁棒性较差。拟牛顿法用海森矩阵的逆矩阵来替代海森矩阵，虽然每次迭代不能保证是最优的优化方向，但是近似矩阵始终是正定的，因此算法总是朝着最优值的方向在搜索。
本文将介绍一种在实际工程中解决大规模优化问题时必然会用到的优化算法：L-BFGS算法。
L-BFGS(Limited-Memory BFGS)是BFGS算法在受限内存时的一种近似算法，而BFGS是数学优化中一种无约束最优化算法。L-BFGS算法就是对拟牛顿算法的一个改进。它的名字已经告诉我们它是基于拟牛顿法BFGS算法的改进。L-BFGS算法的基本思想是：算法只保存并利用最近m次迭代的曲率信息来构造海森矩阵的近似矩阵。
在介绍L-BFGS算法之前，我们先来简单回顾下BFGS算法。

BFGS算法

BFGS校正公式为：

B k + 1 = B k - B k s k s T k B k s T k B k s k + y k y T k y T k s k

$B_{k+1}=B_k-\frac{B_ks_ks_k^TB_k}{s_k^TB_ks_k}+\frac{y_ky_k^T}{y_k^Ts_k}$
利用Sherman-Morrison公式可对上式进行变换，得到:

B - 1 k + 1 = (I - s k y T k y T k s k) T B - 1 k (I - y k s T k y T k s k) + s k s T k y T k s k

$B_{k+1}^{-1}=\left( I-\frac{s_ky_k^T}{y_k^Ts_k} \right)^TB_k^{-1}\left( I-\frac{y_ks_k^T}{y_k^Ts_k} \right)+\frac{s_ks_k^T}{y_k^Ts_k}$
令

Hk+1=B−1k+1 $H_{k+1}=B_{k+1}^{-1}$ ，则得到：

H - 1 k + 1 = (I - s k y T k y T k s k) T H - 1 k (I - y k s T k y T k s k) + s k s T k y T k s k

$H_{k+1}^{-1}=\left( I-\frac{s_ky_k^T}{y_k^Ts_k} \right)^TH_k^{-1}\left( I-\frac{y_ks_k^T}{y_k^Ts_k} \right)+\frac{s_ks_k^T}{y_k^Ts_k}$

BGFS算法存在的问题

在BFGS算法中，需要保存近似Hessian矩阵，特别在高维数据时，会占用大量的存储空间，而在实际运算中，我们只需要方向 $d_k$ ，因此L-BFGS产生了,L-BFGS是对BFGS算法的一种改进算法。在L-BFGS算法中，只保存最近的次迭代信息，以降低数据的存储空间。

L-BFGS算法思路

令 $\rho_k=\frac{1}{y_k^Ts_k},V_k=I-\frac{y_ks_k^T}{y_k^Ts_k}$ ,则BFGS算法中的 $H_{k+1}$ 可以表示为：

H k + 1 = V T k H k V k + ρ k s k s T k

$H_{k+1}=V_k^TH_kV_k+\rho_ks_ks_k^T$
若在初始时，假定初始的矩阵

H0=I $H_0=I$ ，则我们可以得到：

H 1 H 2 \dots H k + 1 = V T 0 H 0 V 0 + ρ 0 s 0 s T 0 = V T 1 H 1 V 1 + ρ 1 s 1 s T 1 = V T 1 (V T 0 H 0 V 0 + ρ 0 s 0 s T 0) V 1 + ρ 1 s 1 s T 1 = V T 1 V T 0 H 0 V 0 V 1 + V T 1 ρ 0 s 0 s T 0 V 1 + ρ 1 s 1 s T 1 = (V T k V T k - 1 \dots V T 1 V T 0) H 0 (V 0 V 1 \dots V k - 1 V k) + (V T k V T k - 1 \dots V T 1) ρ 1 s 1 s T 1 (V 1 \dots V k - 1 V k) + \dots + V T k ρ k - 1 s k - 1 s T k - 1 V k + ρ k s k s T k

$\begin{align} H_1&=V_0^TH_0V_0+\rho_0s_0s_0^T\\ H_2&=V_1^TH_1V_1+\rho_1s_1s_1^T\\ &=V_1^T(V_0^TH_0V_0+\rho_0s_0s_0^T)V_1+\rho_1s_1s_1^T\\ &=V_1^TV_0^TH_0V_0V_1+V_1^T\rho_0s_0s_0^TV_1+\rho_1s_1s_1^T\\ \cdots\\ H_{k+1}&=(V_k^TV_{k-1}^T\cdots V_1^TV_0^T)H_0(V_0V_1\cdots V_{k-1}V_k)\\ &+(V_k^TV_{k-1}^T\cdots V_1^T)\rho_1s_1s_1^T(V_1\cdots V_{k-1}V_k)\\ &+\cdots\\ &+V_k^T\rho_{k-1}s_{k-1}s_{k-1}^TV_k\\ &+\rho_ks_ks_k^T \end{align}$
若此时，只保留最近的m步：

H k + 1 = (V T k V T k - 1 \dots V T k - m) H 0 (V k - m \dots V k - 1 V k) + (V T k V T k - 1 \dots V T k - m) ρ 1 s 1 s T 1 (V k - m \dots V k - 1 V k) + \dots + V T k ρ k - 1 s k - 1 s T k - 1 V k + ρ k s k s T k

$\begin{align} H_{k+1}&=(V_k^TV_{k-1}^T\cdots V_{k-m}^T)H_0(V_{k-m}\cdots V_{k-1}V_k)\\ &+(V_k^TV_{k-1}^T\cdots V_{k-m}^T)\rho_1s_1s_1^T(V_{k-m}\cdots V_{k-1}V_k)\\ &+\cdots\\ &+V_k^T\rho_{k-1}s_{k-1}s_{k-1}^TV_k\\ &+\rho_ks_ks_k^T \end{align}$
这样在L-BFGS算法中，不再保存完整的Hessian矩阵，而是存储向量序列

{sk} $\{s_k\}$ 、

{yk} $\{y_k\}$ ，需要矩阵时，使用向量序列

{sk} $\{s_k\}$ 、

{yk} $\{y_k\}$ 计算就可以得到，而向量序列

{sk} $\{s_k\}$ 、

{yk} $\{y_k\}$ 也不是所有都要保存，只要保存最新的m步向量即可。

$H_k$ 一个在实践中经常用到的有效方法为：

$H_k=r_kI$
$r_k=\frac{s_{k-1}^T y_{k-1}}{y_{k-1}^Ty_{k-1}}$

利用最近一次的曲率信息来估计真实Hessian矩阵的大小，这样使得当前搜索方向较为理想，不至于跑得太偏。

见：
https://wenku.baidu.com/view/cd610728fe4733687e21aae3.html
http://www.tuicool.com/articles/EviQ32m
http://blog.csdn.net/acdreamers/article/details/44728041

L-BFGS算法中的方向的计算方法

q k \leftarrow f o r i α i q i e n d r k - m - 1 f o r β i r i e n d E n d, d k = - r \nabla f k = k - 1 t o k - m d o = ρ i s T i q i + 1 = q i + 1 - α i y i f o r = H 0 q k - m i = k - m, k - m + 1 t o k - 1 d o = ρ i y T i r i - 1 = r i - 1 + s i (α i - β i) f o r T h e r e s u l t i s H k + 1 \nabla f = r

$\begin{align} q_k\leftarrow &\nabla f_k\\ for \quad i &=k-1 \ to \ k-m \ do\\ \alpha_i&=\rho_is_i^Tq_{i+1}\\ q_i&=q_{i+1}-\alpha_iy_i\\ end \quad& for\\ r_{k-m-1}&=H_0q_{k-m}\\ \\ for \quad&i=k-m,k-m+1\ to \ k-1 \ do\\ \beta_i&=\rho_iy_i^Tr_{i-1}\\ r_i &=r_{i-1}+s_i(\alpha_i-\beta_i)\\ end \quad& for\\ End,\quad &The \quad result \quad is \quad H_{k+1}\nabla f=r\\ d_k=-r \end{align}$

算法（参照BFGS算法）

lbfgs.py

#coding:UTF-8
'''
Created on 2017年4月24日

@author: zhangdapeng
'''
from numpy import *
from function import *

def lbfgs(fun, gfun, x0):
    result = []#保留最终的结果
    maxk = 500#最大的迭代次数
    delta = 0.55
    sigma = 0.4

    H0 = eye(shape(x0)[0])

    #s和y用于保存最近m个，这里m取6
    s = []
    y = []
    m = 6
    epsilon=1e-10
    k = 1
    gk = mat(gfun(x0))#计算梯度
    dk = -H0 * gk
    while (k < maxk):             
        n = 0
        mk = 0
        gk = mat(gfun(x0))#计算梯度
        if linalg.norm(gk,1)<epsilon:
            break
        while (n < 20):
            newf = fun(x0 + delta ** n * dk)
            oldf = fun(x0)
            if (newf < oldf + sigma * (delta ** n) * (gk.T * dk)[0, 0]):
                mk = n
                break
            n = n + 1

        #LBFGS校正
        x = x0 + delta ** mk * dk
        #print x

        #保留m个
        if k > m:
            s.pop(0)
            y.pop(0)

        #计算最新的
        sk = x - x0
        yk = gfun(x) - gk

        s.append(sk)
        y.append(yk)

        #two-loop的过程
        t = len(s)
        qk = gfun(x)
        a = []
        for i in range(t):#i值从大到小
            alpha = (s[t - i - 1].T * qk) / (y[t - i - 1].T * s[t - i - 1])
            qk = qk - alpha[0, 0] * y[t - i - 1]
            a.append(alpha[0, 0])
        r = H0 * qk

        for i in range(t):
            beta = (y[i].T * r) / (y[i].T * s[i])
            r = r + s[i] * (a[t - i - 1] - beta[0, 0])


        if (yk.T * sk > 0):
            dk = -r            

        k = k + 1
        x0 = x
        result.append(fun(x0))

    return result

function.py

#coding:UTF-8
'''
Created on 2017年4月24日

@author: zhangdapeng
'''

from numpy import *

#fun  原始函数
def fun(x):  
    return 100 * (x[0,0] ** 2 - x[1,0]) ** 2 + (x[0,0] - 1) ** 2  

#对x1，x2求导后的函数  
def gfun(x):  
    result = zeros((2, 1))  
#     对x1求导
    result[0, 0] = 400 * x[0,0] * (x[0,0] ** 2 - x[1,0]) + 2 * (x[0,0] - 1)  
    result[1, 0] = -200 * (x[0,0] ** 2 - x[1,0])  #对x2求导
    return result

testLBFGS.py

#coding:UTF-8 
'''
Created on 2017年4月24日

@author: zhangdapeng
'''
from lbfgs import *

import matplotlib.pyplot as plt  

x0 = mat([[-1.2], [1]])
result = lbfgs(fun, gfun, x0)
print(result[-1])

n = len(result)
ax = plt.figure().add_subplot(111)
x = arange(0, n, 1)
y = result
ax.plot(x,y)

plt.show()

输出：

2.50088628478e-27

输出图片

这里写图片描述

https://en.wikipedia.org/wiki/Limited-memory_BFGS
http://blog.csdn.net/google19890102/article/details/46389869

大鹏的NLP博客

关注

0
点赞
踩
2

收藏

觉得还不错? 一键收藏
0
评论
拟牛顿法之L-BFGS算法

关于优化算法的求解，比如有梯度下降法，坐标下降法，牛顿法和拟牛顿法。梯度下降法是基于目标函数梯度的，算法的收敛速度是线性的，并且当问题是病态时或者问题规模较大时，收敛速度尤其慢（几乎不适用）；坐标下降法虽然不用计算目标函数的梯度，但是其收敛速度依然很慢，因此它的适用范围也有局限；牛顿法是基于目标函数的二阶导数（海森矩阵）的，其收敛速度较快，迭代次数较少，尤其是在最优值附近时，收敛速度是二
复制链接

扫一扫