放牧代码和思想
专注自然语言处理、机器学习算法
    博主不用扣扣,公事请博客留言,私事请微博私信。开源项目一律GitHub见,发错地方恕不回复,谢谢。

数值优化:理解L-BFGS算法

目录

bfgs.png译自《Numerical Optimization: Understanding L-BFGS》,本来只想作为学习CRF的补充材料,读完后发现收获很多,把许多以前零散的知识点都串起来了。对我而言,的确比零散地看论文要轻松得多。原文并没有太多关注实现,对实现感兴趣的话推荐原作者的golang实现梯度.png

数值优化是许多机器学习算法的核心。一旦你确定用什么模型,并且准备好了数据集,剩下的工作就是训练了。估计模型的参数(训练模型)通常归结为最小化一个多元函数f(x).png,其中输入CodeCogsEqn (3).png是一个高维向量,也就是模型参数。换句话说,如果你求解出:

CodeCogsEqn (1).png

那么CodeCogsEqn (3).png*就是最佳的模型参数(当然跟你选择了什么目标函数有关系)。 

在这篇文章中,我将重点放在讲解L-BFGS算法的无约束最小化上,该算法在一些能用上批处理优化的ML问题中特别受欢迎。对于更大的数据集,则常用SGD方法,因为SGD只需要很少的迭代次数就能达到收敛。在以后的文章中,我可能会涉及这些技术,包括我个人最喜欢的AdaDelta 。

注 : 在整个文章中,我会假设你记得多元微积分。所以,如果你不记得什么是梯度或海森矩阵,你得先复习一下。

牛顿法

大多数数值优化算法都是迭代式的,它们产生一个序列,该序列最终收敛于CodeCogsEqn (4).png,使得f(x).png达到全局最小化。假设,我们有一个估计屏幕快照 2016-08-11 下午8.10.12.png,我们希望我们的下一个估计屏幕快照 2016-08-11 下午8.11.00.png有这种属性:屏幕快照 2016-08-11 下午8.11.27.png

牛顿的方法是在点屏幕快照 2016-08-11 下午8.10.12.png附近使用二次函数近似屏幕快照 2016-08-11 下午8.13.06.png。假设屏幕快照 2016-08-11 下午8.13.06.png是二次可微的,我们可以用屏幕快照 2016-08-11 下午8.13.06.png在点屏幕快照 2016-08-11 下午8.14.16.png的泰勒展开来近似屏幕快照 2016-08-11 下午8.13.06.png

屏幕快照 2016-08-11 下午8.15.12.png

其中,屏幕快照 2016-08-11 下午8.17.05.png屏幕快照 2016-08-11 下午8.17.24.png分别为目标函数屏幕快照 2016-08-11 下午8.13.06.png在点屏幕快照 2016-08-11 下午8.10.12.png处的梯度和Hessian矩阵。当屏幕快照 2016-08-11 下午8.18.28.png时,上面的近似展开式是成立的。你可能记得微积分中一维泰勒多项式展开,这是其推广。

为了简化符号,将上述二次近似记为屏幕快照 2016-08-11 下午8.21.18.png,我们把生成屏幕快照 2016-08-11 下午8.21.18.png这样的二次近似的迭代算法中的一些概念简记如下:

不失一般性,我们可以记屏幕快照 2016-08-11 下午8.19.33.png,那么上式可以写作:

屏幕快照 2016-08-11 下午8.22.34.png

其中屏幕快照 2016-08-11 下午8.23.12.png屏幕快照 2016-08-11 下午8.24.27.png分别表示目标函数在点屏幕快照 2016-08-11 下午8.10.12.png处的梯度和Hessian矩阵。

我们想找一个屏幕快照 2016-08-11 下午8.32.42.png,使得屏幕快照 2016-08-11 下午8.13.06.png屏幕快照 2016-08-11 下午8.10.12.png的二次近似最小。上式对屏幕快照 2016-08-11 下午8.32.42.png求导:

屏幕快照 2016-08-11 下午8.35.10.png

任何使得屏幕快照 2016-08-11 下午8.38.01.png屏幕快照 2016-08-11 下午8.38.31.png都是屏幕快照 2016-08-11 下午8.38.53.png的局部极值点,如果我们假设屏幕快照 2016-08-11 下午8.13.06.png是凸函数,则屏幕快照 2016-08-11 下午8.24.27.png正定的,那么局部极值点就是全局极值点(凸二次规划)。

解出屏幕快照 2016-08-11 下午8.38.31.png

屏幕快照 2016-08-11 下午8.54.55.png

这就得到了一个很好的搜索方向,在实际应用中,我们一般选择一个步长α,即按照下式更新屏幕快照 2016-08-11 下午8.11.00.png

屏幕快照 2016-08-11 下午8.58.16.png

使得屏幕快照 2016-08-11 下午8.59.06.png相比f(x).png的减小量最大化。

迭代算法伪码:

屏幕快照 2016-08-11 下午9.03.34.png

步长α的确定可以采用任何line search算法,其中最简单的一种是backtracking line search。该算法简单地选取越来越小的步长α,直到屏幕快照 2016-08-11 下午8.13.06.png的值小到满意为止。关于line search算法的详情请参考Line Search Methods.pdfLecture 5- Gradient Descent.pdf

在软件工程上,我们可以将牛顿法视作实现了下列Java接口的一个黑盒子:

public interface TwiceDifferentiableFunction
{
    // compute f(x)
    double valueAt(double[] x);

    // compute grad f(x)
    double[] gradientAt(double[] x);

    // compute inverse hessian H^-1
    double[][] inverseHessian(double[] x);
}

如果你有兴趣,你还可以通过一些枯燥无味的数学公式,证明对任意一个凸函数,上述算法一定可以收敛到一个唯一的最小值屏幕快照 2016-08-11 下午9.54.52.png,且不受初值屏幕快照 2016-08-11 下午9.55.27.png的影响。对于非凸函数,上述算法仍然有效,但只能保证收敛到一个局部极小值。在上述算法于非凸函数的实际应用中,用户需要注意初值的选取以及其他算法细节。

巨大的海森矩阵

牛顿法最大的问题在于我们必须计算海森矩阵的逆。注意在机器学习应用中,屏幕快照 2016-08-11 下午8.13.06.png的输入的维度常常与模型参数对应。十万维度的参数并不少见(SVM中文文本分类取词做特征的话,就在十万这个量级),在一些图像识别的场景中,参数可能上十亿。所以,计算海森矩阵或其逆并不现实。对许多函数而言,海森矩阵可能根本无法计算,更不用说表示出来求逆了。

所以,在实际应用中牛顿法很少用于大型的优化问题。但幸运的是,即便我们不求出屏幕快照 2016-08-11 下午8.13.06.png屏幕快照 2016-08-11 下午8.10.12.png的精确屏幕快照 2016-08-12 下午8.38.44.png,而使用一个近似的替代值,上述算法依然有效。

拟牛顿法

如果不求解屏幕快照 2016-08-11 下午8.13.06.png屏幕快照 2016-08-11 下午8.10.12.png的精确屏幕快照 2016-08-12 下午8.38.44.png,我们要使用什么样的近似呢?我们使用一种叫QuasiUpdate的策略来生成屏幕快照 2016-08-12 下午8.38.44.png的近似,先不管QuasiUpdate具体是怎么做的,有了这个策略,牛顿法进化为如下的拟牛顿法:

屏幕快照 2016-08-12 下午8.50.21.png

跟牛顿法相比,只是把屏幕快照 2016-08-12 下午8.38.44.png的计算交给了QuasiUpdate。为了辅助QuasiUpdate,计算了几个中间变量。QuasiUpdate只需要上个迭代的屏幕快照 2016-08-12 下午8.38.44.png、输入和梯度的变化量(屏幕快照 2016-08-12 下午8.54.12.png屏幕快照 2016-08-12 下午8.54.18.png)。如果QuasiUpdate能够返回精确的屏幕快照 2016-08-12 下午8.57.16.png的逆,则拟牛顿法等价于牛顿法。

在软件工程上,我们又可以写一个黑盒子接口,该接口不再需要计算海森矩阵的逆,只需要在内部更新它,再提供一个矩阵乘法的接口即可。事实上,内部如何处理,外部根本无需关心。用Java表示如下:

public interface DifferentiableFunction
{
    // compute f(x)
    double valueAt(double[] x);

    // compute grad f(x)
    double[] gradientAt(double[] x);
}

public interface QuasiNewtonApproximation
{
    // update the H^{-1} estimate (using x_{n+1}-x_n and grad_{n+1}-grad_n)
    void update(double[] deltaX, double[] deltaGrad);
    
    // H^{-1} (direction) using the current H^{-1} estimate
    double[] inverseHessianMultiply(double[] direction);
}

注意我们唯一用到海森矩阵的逆的地方就是求它与梯度的乘积,所以我们根本不需要在内存中将其显式地、完全地表示出来。这对接下来要阐述的L-BFGS特别有用。如果你对实现细节感兴趣,可以看看作者的golang实现。

近似海森矩阵

QuasiUpdate到底要如何近似海森矩阵呢?如果我们让QuasiUpdate忽略输入参数,直接返回单位矩阵,那么拟牛顿法就退化成了梯度下降法了,因为函数减小的方向永远是梯度屏幕快照 2016-08-11 下午8.17.05.png。梯度下降法依然能保证凸函数屏幕快照 2016-08-11 下午8.13.06.png收敛到全局最优对应的CodeCogsEqn (4).png,但直觉告诉我们,梯度下降法没有利用到屏幕快照 2016-08-11 下午8.13.06.png的二阶导数信息,收敛速度肯定更慢了。

我们先把屏幕快照 2016-08-11 下午8.13.06.png的二次近似写在下面,从这里找些灵感。

屏幕快照 2016-08-11 下午8.22.34.png

Secant Condition

屏幕快照 2016-08-12 下午9.28.30.png的一个性质是,它的梯度与屏幕快照 2016-08-11 下午8.13.06.png屏幕快照 2016-08-11 下午8.10.12.png处的梯度一致(近似函数的梯度与原函数的梯度一致,这才叫近似嘛)。也就是说我们希望保证:

屏幕快照 2016-08-12 下午9.31.14.png

我们做个减法:

屏幕快照 2016-08-12 下午9.33.53.png

由中值定理,我们有:

屏幕快照 2016-08-12 下午9.37.12.png

这个式子就是所谓的Secant Condition,该条件保证屏幕快照 2016-08-12 下午9.49.21.png至少对屏幕快照 2016-08-12 下午9.50.06.png而言是近似海森矩阵的。

等式两边同时乘以屏幕快照 2016-08-12 下午8.38.44.png,并且由于我们定义过:

屏幕快照 2016-08-12 下午9.39.06.png

于是我们得到:

屏幕快照 2016-08-12 下午9.40.50.png

对称性

由定义知海森矩阵是函数的二阶偏导数矩阵,即屏幕快照 2016-08-12 下午9.44.57.png,所以海森矩阵一定是对称的。

BFGS更新

形式化地讲,我们希望屏幕快照 2016-08-11 下午8.24.27.png至少满足上述两个条件:

  1. 屏幕快照 2016-08-12 下午8.54.12.png屏幕快照 2016-08-12 下午8.54.18.png而言满足Secant Condition

  2. 满足对称性

给定上述两个条件,我们还希望屏幕快照 2016-08-11 下午8.24.27.png相较于屏幕快照 2016-08-12 下午9.53.27.png的变化量最小。这类似“ MIRA 更新”,我们有许多满足条件的选项,但我们只选那个变化最小的。这种约束形式化地表述如下:

屏幕快照 2016-08-12 下午9.55.33.png

上面的范数屏幕快照 2016-08-12 下午9.57.15.png表示weighted frobenius norm。这个约束最小化问题的解是:

屏幕快照 2016-08-12 下午9.59.02.png

式中屏幕快照 2016-08-12 下午10.00.37.png。我不知道如何推导它,推导的话需要用很多符号,并且费时费力。

这种更新算法就是著名的Broyden–Fletcher–Goldfarb–Shanno (BFGS)算法,该算法是取发明者名字的首字母命名的。

bfgs.png

关于BFGS,有一些需要注意的点:

  • 只要屏幕快照 2016-08-12 下午8.38.44.png是正定的,屏幕快照 2016-08-12 下午10.11.17.png就一定是正定的。所以我们只需要选择一个正定的屏幕快照 2016-08-12 下午10.13.05.png即可,甚至可以选择单位矩阵。

  • 屏幕快照 2016-08-12 下午8.38.44.png屏幕快照 2016-08-12 下午10.11.17.png还存在简单的算术关系,给定屏幕快照 2016-08-12 下午10.11.17.png屏幕快照 2016-08-12 下午8.54.12.png屏幕快照 2016-08-12 下午8.54.18.png,我们就能倒推出屏幕快照 2016-08-12 下午8.38.44.png

把这些知识放到一起,我们就得出了BFGS更新的算法,给定方向d,该算法可以计算出屏幕快照 2016-08-12 下午10.17.59.png,却不需要求屏幕快照 2016-08-12 下午8.38.44.png矩阵,只需要按照上述表达式不断地递推即可:

屏幕快照 2016-08-12 下午10.19.22.png

由于屏幕快照 2016-08-12 下午8.38.44.png的唯一作用就是计算屏幕快照 2016-08-12 下午10.21.08.png,我们只需用该更新算法就能实现拟牛顿法。

L-BFGS:省内存的BFGS

BFGS拟牛顿近似算法虽然免去了计算海森矩阵的烦恼,但是我们仍然需要保存每次迭代的屏幕快照 2016-08-12 下午8.54.12.png屏幕快照 2016-08-12 下午8.54.18.png的历史值。这依然没有减轻内存负担,要知道我们选用拟牛顿法的初衷就是减小内存占用。

L-BFGS是limited BFGS的缩写,简单地只使用最近的m个屏幕快照 2016-08-12 下午8.54.12.png屏幕快照 2016-08-12 下午8.54.18.png记录值。也就是只储存屏幕快照 2016-08-12 下午10.29.13.png屏幕快照 2016-08-12 下午10.29.32.png,用它们去近似计算屏幕快照 2016-08-12 下午10.21.08.png。初值屏幕快照 2016-08-12 下午10.30.56.png依然可以选取任意对称的正定矩阵。

L-BFGS改进算法

在实际应用中有许多L-BFGS的改进算法。对不可微分的函数,可以用 othant-wise 的L-BFGS改进算法来训练屏幕快照 2016-08-12 下午10.34.38.png正则损失函数。

不在大型数据集上使用L-BFGS的原因之一是,在线算法可能收敛得更快。这里甚至有一个L-BFGS的在线学习算法,但据我所知,在大型数据集上它们都不如一些SGD的改进算法(包括 AdaGrad 或 AdaDelta)的表现好。

Reference

http://aria42.com/blog/2014/12/understanding-lbfgs


知识共享许可协议 知识共享署名-非商业性使用-相同方式共享码农场 » 数值优化:理解L-BFGS算法

分享到:更多 ()

评论 欢迎留言

  • 昵称 (必填)
  • 邮箱 (必填)
  • 网址

我的开源项目

HanLP自然语言处理包基于DoubleArrayTrie的Aho Corasick自动机