OWLQN算法介绍,及go语言实现logistic regression优化

qm1004 · · 5293 次点击 · · 开始浏览    
这是一个创建于 的文章,其中的信息可能已经有所发展或是发生改变。

用logistic regression做ctr预估或则分类的同学,大家应该都知道梯度下降能够求得目标问题的最优解。但是梯度下降迭代太慢,而且不一定能求得最优解。

owlqn不仅收敛速度比梯度下降算法快很多的,而且能对L1正则化的优化问题求解,L1正则化使目标函数在变量为0时不可微。OWLQN是LBFGS在L1问题上的解决方案(因为LBFGS只能用来优化目标函数为凸的情况)。L1相对L2正则化,具有特征选择的优势。下面github的链接是我用go实现的单机owlqn算法。

https://github.com/qm1004/OWLQN

理论方面的介绍网上已经有非常好的简介了,算法的介绍主要转自下面这篇文章。

转自

http://www.cnblogs.com/downtjs/p/3222643.html 

OWLQN算法

一、BFGS算法

      算法思想如下:

           Step1   取初始点x^{(0)} ,初始正定矩阵H_0,允许误差\epsilon>0,令k=0

           Step2   计算p^{(k)}=-H_k \nabla f(x^{(k)})

           Step3   计算\alpha_k>0,使得

                                               f(x^{(k)}+\alpha_kp^{(k)})={min}\limit_{\alpha \geq 0} f(x^{(k)}+\alpha p^{(k)})

          Step4    令x^{(k+1)}=x^{(k)}+\alpha_k p^{(k)}

          Step5    如果||\nabla f(x^{(k+1)})|| \leq \epsilon,则取x^{(k+1)}为近似最优解;否则转下一步;

          Step6    计算

                                s_k=x^{(k+1)}-x^{(k)}y_k=\nabla f(x^{(k+1)})-\nabla f(x^{(k)})

                         H_{k+1}=H_k+\frac{1}{s_k^Ty_k}(1+\frac{y_k^TH_ky_k}{s_k^Ty_k})s_ks_k^T-\frac{1}{s_k^Ty_k}(s_ky_k^TH_k+H_ky_ks_k)

          令k=k+1,转Step2.

优点:

1、不用直接计算Hessian矩阵;

2、通过迭代的方式用一个近似矩阵代替Hessian矩阵的逆矩阵。

缺点:

1、矩阵存储量为n^2,因此维度很大时内存不可接受;

2、矩阵非稀疏会导致训练速度慢。

 

二、L-BFGS算法

      针对BFGS的缺点,主要在于如何合理的估计出一个Hessian矩阵的逆矩阵,L-BFGS的基本思想是只保存最近的m次迭代信息,从而大大降低数据存储空间。对照BFGS,我重新整理一下用到的公式:

                                 \rho_k=\frac{1}{y_{k}^T s_k}     

                                 s_k=x_k-x_{k-1}            

                                 y_k=\nabla{f(x_k)}-\nabla{f(x_{k-1})}

                                 V_k=I-\rho_{k}y_{k}s_{k}^T

于是估计的Hessian矩阵逆矩阵如下:                         

                                H_k=(I-\rho_{k-1}s_{k-1}y_{k-1}^T)H_{k-1}(I-\rho_{k-1}y_{k-1}s_{k-1}^T)+s_{k-1}\rho_{k-1}s_{k-1}^T

                                       =V_{k-1}^TH_{k-1}V_{k-1}+ s_{k-1}\rho_{k-1}s_{k-1}^T

                                H_{k-1}=V_{k-2}^TH_{k-2}V_{k-2}+ s_{k-2}\rho_{k-2}s_{k-2}^T

带入上式,得:

                               H_k=V_{k-1}^TV_{k-2}^TH_{k-2}V_{k-2}V_{k-1}+ V_{k-1}^Ts_{k-2}\rho_{k-2}s_{k-2}^T V_{k-1}+s_{k-1}\rho_{k-1}s_{k-1}^T 

假设当前迭代为k,只保存最近的m次迭代信息,(即:从k-m~k-1),依次带入H,得到:

公式1:

                               H_k=(V_{k-1}^TV_{k-2}^T\ldots V_{k-m}^T) H_k^{0}(V_{k-m}\ldots V_{k-2}V_{k-1})

                                      + (V_{k-1}^TV_{k-2}^T\ldots V_{k-m+1}^T) S_{k-m}\rho_{k-m}S_{k-m}^T (V_{k-m}\ldots V_{k-2}V_{k-1})

                                      + (V_{k-1}^TV_{k-2}^T\ldots V_{k-m+2}^T) S_{k-m+1}\rho_{k-m+1}S_{k-m+1}^T (V_{k-m+1}\ldots V_{k-2}V_{k-1})

                                      +\ldots

                                      +V_{k-1}^T s_{k-2}\rho_{k-2} s_{k-2}^TV_{k-1}  

                                      +s_{k-1}\rho_{k-1}s_{k-1}^T

算法第二步表明了上面推导的最终目的:找到第k次迭代的可行方向,满足:

                                p_k=-H_k\nabla f(x_k)

为了求可行方向p,有下面的:

  two-loop recursion算法


                                  q=\nabla f(x_k)

                                  for (i=1 \ldots m) \quad \quad \quad do

                                          \alpha_i=\rho_{k-i}s_{k-i}^Tq

                                          q=q-\alpha_iy_{k-i}

                                   end \quad \quad \quad for

                                   r=H_k^{0}q

                                   for (i=m \ldots 1) \quad \quad \quad do

                                         \beta=\rho_{k-i}y_{k-i}^Tr

                                         r=r+s_{k-i}(\alpha_i-\beta)

                                    end \quad \quad \quad for

                                   return \quad \quad \quad r


该算法的正确性推导:

1、令:     q_0=\nabla f(x_k),递归带入q:

                                q_i=q_{i-1}-\rho_{k-i}y_{k-i}s_{k-i}^Tq_{i-1}

                                     =(I-\rho_{k-i}y_{k-i}s_{k-i}^T)q_{i-1}

                                     =V_{k-i}q_{i-1}

                                     =V_{k-i}V_{k-i+1}q_{i-2}

                                     =\ldots

                                     =V_{k-i}V_{k-i+1} \ldots V_{k-1} q_0

                                     =V_{k-i}V_{k-i+1} \ldots V_{k-1} \nabla f(x_k)

相应的:

                                \alpha_i=\rho_{k-i}s_{k-i}^Tq_{i-1}

                                     =\rho_{k-i}s_{k-i}^T V_{k-i+1}V_{k-i+2} \ldots V_{k-1} \nabla f(x_k)

2、令:r_{m+1}=H_{k-m}q=H_{k-m}V_{k-i}V_{k-i+1} \ldots V_{k-1} \nabla f(x_k)

                                r_i=r_{i+1}+s_{k-i}(\alpha_i-\beta) =r_{i+1}+s_{k-i}(\alpha_i-\rho_{k-i}y_{k-i}^Tr_{i+1})

                                     =s_{k-i}\alpha_i+(I-s_{k-i}\rho_{k-i}y_{k-i}^T)r_{i+1}

                                     =s_{k-i}\alpha_{i}+V_{k-i}^Tr_{i+1}

于是:

                                r_1=s_{k-1}\alpha_1+V_{k-1}^Tr_2 =s_{k-1}\rho_{k-1}s_{k-1}^T \nabla f(x_k)+V_{k-1}^Tr_2

                                     =s_{k-1}\rho_{k-1}s_{k-1}^T \nabla f(x_k)+V_{k-1}^T(s_{k-2}\alpha_2+V_{k-2}^Tr_3)

                                     =s_{k-1}\rho_{k-1}s_{k-1}^T \nabla f(x_k)+V_{k-1}^Ts_{k-2}\rho_{k-2}s_{k-2}^TV_{k-1}\nabla f(x_k)+V_{k-1}^T V_{k-2}^T r_3

                                     =\ldots

                                     =s_{k-1}\rho_{k-1}s_{k-1}^T \nabla f(x_k)           

                                      +V_{k-1}^T s_{k-2}\rho_{k-2} s_{k-2}^TV_{k-1}  \nabla f(x_k)

                                      +\ldots     

                                      + (V_{k-1}^TV_{k-2}^T\ldots V_{k-m+2}^T) S_{k-m+1}\rho_{k-m+1}S_{k-m+1}^T (V_{k-m+1}\ldots V_{k-2}V_{k-1}) \nabla f(x_k)

                                      + (V_{k-1}^TV_{k-2}^T\ldots V_{k-m+1}^T) S_{k-m}\rho_{k-m}S_{k-m}^T (V_{k-m}\ldots V_{k-2}V_{k-1}) \nabla f(x_k)

                                      +(V_{k-1}^TV_{k-2}^T\ldots V_{k-m}^T) H_{k-m}(V_{k-m}\ldots V_{k-2}V_{k-1}) \nabla f(x_k)

这个two-loop recursion算法的结果和公式1*初始梯度的形式完全一样,这么做的好处是:

1、只需要存储s_{k-i}y_{k-i} (i=1~m);

2、计算可行方向的时间复杂度从O(n*n)降低到了O(n*m),当m远小于n时为线性复杂度。

总结L-BFGS算法的步骤如下:

      Step1:       选初始点x_0,允许误差\epsilon >0,存储最近迭代次数m(一般取6);

      Step2:       k=0, \quad \quad \quad H_0=I ,\quad \quad \quad r=\nabla f(x_{0})

      Step3:       如果 ||\nabla f(x_{k+1})||\leq \epsilon 则返回最优解x_{k+1},否则转Step4;

      Step4:       计算本次迭代的可行方向:p_k=-r _k

      Step5:       计算步长\alpha_k>0,对下面式子进行一维搜索:

                         f(x_k+\alpha_kp_k)={min}\limits_{\alpha \geq 0} \quad \quad \quad f(x_k+\alpha p_k)

      Step6:       更新权重x:

                         x_{k+1}=x_k+\alpha_kp_k ;     

      Step7:      if k > m

                             只保留最近m次的向量对,需要删除(s_{k-m},y_{k-m});

      Step8:       计算并保存:

                        s_k=x_{k+1}=x_k

                        y_k=\nabla f(x_{k+1})-\nabla f(x_k)

      Step9:       用two-loop recursion算法求得:

                         r_k=H_k\nabla f(x_k)

      k=k+1,转Step3。

需要注意的地方,每次迭代都需要一个H_{k-m} ,实践当中被证明比较有效的取法为:

                          H_k^0=\gamma_k I

                          \gamma_k=\frac{s_{k-1}^Ty_{k-1}}{{y_{k-1}^Ty_{k-1}}

 

三、OWL-QN算法

1、问题描述

对于类似于Logistic Regression这样的Log-Linear模型,一般可以归结为最小化下面这个问题:

                          J(x)=l(x)+r(x)

其中,第一项为loss function,用来衡量当训练出现偏差时的损失,可以是任意可微凸函数(如果是非凸函数该算法只保证找到局部最优解),后者为regularization term,用来对模型空间进行限制,从而得到一个更“简单”的模型。

        根据对模型参数所服从的概率分布的假设的不同,regularization term一般有:L1-norm(模型参数服从Gaussian分布);L2-norm(模型参数服从Laplace分布);以及其他分布或组合形式。

L2-norm的形式类似于:

                         J(x)=l(x)+C\sum\limit_i{x_i^2}

L1-norm的形式类似于:

                         J(x)=l(x)+C\sum\limit_i{|x_i|}

L1-norm和L2-norm之间的一个最大区别在于前者可以产生稀疏解,这使它同时具有了特征选择的能力,此外,稀疏的特征权重更具有解释意义。

对于损失函数的选取就不在赘述,看两幅图:

image

图1 - 红色为Laplace Prior,黑色为Gaussian Prior 

 

image

图2 直观解释稀疏性的产生

 

        对LR模型来说损失函数选取凸函数,那么L2-norm的形式也是的凸函数,根据最优化理论,最优解满足KKT条件,即有:\nabla J(x^*)=0,但是L1-norm的regularization term显然不可微,怎么办呢?

2、Orthant-Wise Limited-memory Quasi-Newton

        OWL-QN主要是针对L1-norm不可微提出的,它是基于这样一个事实:任意给定一个维度象限,L1-norm 都是可微的,因为此时它是一个线性函数:

image

图3 任意给定一个象限后的L1-norm

OWL-QN中使用了次梯度决定搜索方向,凸函数不一定是光滑而处处可导的,但是它又符合类似梯度下降的性质,在多元函数中把这种梯度叫做次梯度,见维基百科http://en.wikipedia.org/wiki/Subderivative

举个例子:

File:Subderivative illustration.png

图4 次导数

对于定义域中的任何x0,我们总可以作出一条直线,它通过点(x0, f(x0)),并且要么接触f的图像,要么在它的下方。这条直线的斜率称为函数的次导数,推广到多元函数就叫做次梯度。

次导数及次微分:

凸函数f:I→R在点x0的次导数,是实数c使得:

                        f(x)-f(x_0)\ge c(x-x_0)

对于所有I内的x。可以证明,在点x0的次导数的集合是一个非空闭区间[a, b],其中a和b是单侧极限

              a=\lim_{x\to x_0^-}\frac{f(x)-f(x_0)}{x-x_0}
              b=\lim_{x\to x_0^+}\frac{f(x)-f(x_0)}{x-x_0}

它们一定存在,且满足a ≤ b。所有次导数的集合[a, b]称为函数f在x0的次微分。

OWL-QN和传统L-BFGS的不同之处在于:

1)、利用次梯度的概念推广了梯度

       定义了一个符合上述原则的虚梯度,求一维搜索的可行方向时用虚梯度来代替L-BFGS中的梯度:

                      image

                        image

                       DPB1$`%~T9U0]2%YSCVQ83K

怎么理解这个虚梯度呢?见下图:

对于非光滑凸函数,那么有这么几种情况:

image

图5  \partial_i^-f(x)>0

 

image

图6  \partial_i^+f(x)<0

 

image

图7  otherwise

2)、一维搜索要求不跨越象限

       要求更新前权重与更新后权重同方向:

           image

图8  OWL-QN的一次迭代

总结OWL-QN的一次迭代过程:

–Find vector of steepest descent

–Choose sectant

–Find L-BFGS quadratic approximation

–Jump to minimum

–Project back onto sectant

–Update Hessian approximation using gradient of loss alone

最后OWL-QN算法框架如下:

                                 8]O)3Y(ZOLX5C1AI9YXP72M

 

      与L-BFGS相比,第一步用虚梯度代替梯度,第二、三步要求一维搜索不跨象限,也就是迭代前的点与迭代后的点处于同一象限,第四步要求估计Hessian矩阵时依然使用loss function的梯度(因为L1-norm的存在与否不影响Hessian矩阵的估计)。

四、参考资料

1、Galen Andrew and Jianfeng Gao. 2007. 《Scalable training of L1-regularized log-linear models》. In Proceedings of ICML, pages 33–40.

2、http://freemind.pluskid.org/machine-learning/sparsity-and-some-basics-of-l1-regularization/#d20da8b6b2900b1772cb16581253a77032cec97e

3、http://research.microsoft.com/en-us/downloads/b1eb1016-1738-4bd5-83a9-370c9d498a03/default.aspx


有疑问加站长微信联系(非本文作者)

本文来自:CSDN博客

感谢作者:qm1004

查看原文:OWLQN算法介绍,及go语言实现logistic regression优化

入群交流(和以上内容无关):加入Go大咖交流群,或添加微信:liuxiaoyan-s 备注:入群;或加QQ群:692541889

5293 次点击  
加入收藏 微博
被以下专栏收入,发现更多相似内容
添加一条新回复 (您需要 登录 后才能回复 没有账号 ?)
  • 请尽量让自己的回复能够对别人有帮助
  • 支持 Markdown 格式, **粗体**、~~删除线~~、`单行代码`
  • 支持 @ 本站用户;支持表情(输入 : 提示),见 Emoji cheat sheet
  • 图片支持拖拽、截图粘贴等方式上传