线性收敛的随机L-BFGS算法

线性收敛的随机L-BFGS算法

以下皆为翻译A Linearly-Convergent Stochastic L-BFGS Algorithm 原文链接
六级没过,莫怪

概要

我们提出了一个新的随机L-BFGS算法,并证明其对强凸连续函数可以达到线性收敛.我们的算法很大程度上基于最近在Byrd(2014)提出的L-BFGS算法的随机变种,以及Johnson和Zhang最近所提出的对于随机梯度下降的方差约减.我们用实验证明了我们的算法在大规模凸优化以及非凸问题中表现良好,展现了线性收敛率以及能够快速求出优化问题中的高精度解.此外,我们的算法能够适用于广泛的不同数量级的学习速率.

1 介绍

机器学习的一个趋势是使用更多的参数来模拟更大的数据集。因此,为这些大规模问题设计优化算法非常重要。这种情况下出现的典型优化问题是经验风险最小化。形式如下:

(28)minw1Ni=1Nfi(w),

其中ω∈Rd可以指定为机器学习模型中的参数, fi(ω)量化模型ω拟合数据的程度.当想要求解式子(1)的时候会面临两个问题.第一点,维度d可能非常大.第二点,N可能非常大.

当d很小的时候, 由于牛顿法其(理论和实践的中)所展现出的快速的收敛性,经常选择他作为优化算法.然而, ,在高纬度数据中,牛顿法要求计算Hessian矩阵以及他的逆矩阵,计算可能太过昂贵.结果,从业者通常仅使用一阶优化方法,他只需要计算目标的梯度,每次迭代O(d)的复杂度.梯度下降是最简单的一阶优化算法,但是已经有大量工作用来构造拟牛顿法,拟牛顿法能够不计算二阶导数但却可以提取目标函数的曲率信息. L-BFGS(Liu and Nocedal,1989),是经典BFGS算法的限制内存版本, 是该领域最成功的算法之一.不精确的拟牛顿法是另外一种利用二阶信息进行大规模优化的算法.他们在O(d)步中近似的翻转Hessian.这可以通过共轭梯度法在固定步数中完成.

当N很大的时候,批处理算法,比如梯度算法,每轮迭代需要计算全部目标函数的梯度,更新模型之前必须处理一遍所有的数据,故而会变得很慢.随机优化算法, 通过每处理一部分数据的时候就更新模型克服了这个问题,这让他们能够在一次全梯度的时间里取得更多进展.

对于大多数维度d和N都很大的机器学习问题,随机梯度下降(SGD)以及他的相关变种算法是应用的最广的算法,大多因为他们是少数几种能够实际应用于这种场景下的算法.

鉴于这种情况,已经有很多优化领域的研究致力于设计更好的一阶随机优化算法.有关部分,可以参考(Kingma and Ba, 2015; Sutskever et al., 2013; Duchi et al., 2011;Shalev-Shwartz and Zhang, 2013; Johnson and Zhang,2013; Roux et al.,2012; Wang et al., 2013; Nesterov,2009; Frostig et al., 2015; Agarwal et al., 2014).

不像梯度下降算法,L-BFGS算法不能够立即转换成随机版本.随机梯度下降算法中的更新从期望上平均值能够产生一个下降的方向.然而,正如Byrd所指出,L-BFGS中用构造Hessian逆的近似的更新是覆盖了另一个量,而不是平均.我们的算法采用相同的方法,通过更大的minibatch来计算Hessian vector乘积解决了这个问题,从而解决了这个问题.

尽管随机方法在早期通常会取得迅速进展,但是梯度估计的方差减慢了他们在最优解附近的收敛.现在来解释这种现象,即使SGD初始解在最优解附近,他都会立刻移动到目标函数值更差的一个点.鉴于这种原因,他的收敛通常需要逐渐减少的步长来保证. 通过减少梯度估计的方差是加速一阶随机方法收敛性的标志性工作.

我们引入了随机L-BFGS并结合方差约减想法的变种算法,它具有两个理想特性.首先,他保证了强凸条件下的线性收敛率.特别是,他不需要一个递减步长来保证收敛(部分事实证明,如果我们的算法在最优解初始化,那么他会留在那里).第二点,他展现了高质量的线性收敛率,在大规模机器学中表现得很好.

2 算法

我们考虑最小化下列函数的问题

(29)f(w)=1Ni=1Nfi(w)

其中ω∈Rd.对于子集S⊂{1,2,…,N},我们定义下列子采样目标函数fS
(30)fS(w)=1|S|iSfi(w).

我们的更新,将采用梯度fS作为随机梯度估计,也将随机来近似Hessian的逆.和Byrd et al.(2014)保持一致,我们采用不同的子集S,T⊂{1,…,N}来将梯度的估计从Hessian的估计中解耦.记b=|S|, bH=|T|.采用Johnson and Zhang (2013)一样, 我们隔段时间来计算一次全梯度, 以此来减少梯度估计中的方差.

我们算法中的更新规则将采用如下形式

ωk+1=ωkηkHkvk

在梯度下降算法中,Hk是单位矩阵.在牛顿法中,他是Hessian矩阵的逆(2f(wk))1.在我们的算法中, 和L-BFGS一样,H_k是对Hessian矩阵逆的拟合.不同于一般随机梯度估计,νk是方差约减后的随机梯度估计.

我们算法的伪代码在Algorithm 1中给出.他由若干个参数确定.要求的参数有学习率η 内存大小M, 正整数m和L.每m次迭代,算法求一次全梯度,用来减少随机梯度估计的方差.每L次迭代,算法会更新Hessian逆的拟合.向量s_r记录算过过去2L次迭代中的平均方向.向量y_r是由s_r和Hessian矩阵逆的拟合矩阵相乘得到的.需要注意的是这里和传统的L-BFGS不同,传统的L-BFGS中的y_r是由两个相邻梯度的差获得.我们发现这种方法在随机的背景下效果更好.Hessian矩阵逆的近似矩阵H_r是由2.1中所描述的标准L-BFGS算法中(s_j,y_j )所确定,其中r-M+1≤j≤r.使用者还必须选择批量大小b和b_H,来构造随机梯度估计和随机Hessian估计.

SLBFGS
在算法1中,我们用I来指代单位矩阵.我们用Fk,t表示当计数器k和t满足一定条件时随机变量所产生的sigma代数.形式如下

Fk,t=σ({Sk,t:k<kork=kandt<t}{Tr:rLmk+t}).

我们将用Ek,t 来表示Fk,t的条件期望.
我们在2.1节定义Hessian矩阵逆的近似矩阵为H_r.需要注意的是我们不直接构造矩阵H_r,因为这样要求O(d2)的计算复杂度.实践中,我们直接用两层循环计算H_r v的乘积.

2.1 构造Hessian逆的近似H_r

我们根据传统的L-BFGS方法来用(sj,yj)键值对来构造Hessian.记pj=1/sjTyjT,然后当r-M+1≤j≤r时候,递归定义:

\begin{equation} \label{eq:inv_hess_update}  H_r^{(j)} = (I - \rho_j s_j y_j^{\top})^{\top}H_r^{(j-1)}(I - \rho_j s_j y_j^{\top}) + \rho_j s_j s_j^{\top} ,\end{equation}

初始化Hr(rM)=(sryr/yr2)I 以及Hr=Hr(r)
注意上述等式(4)中的更新保留了正定性(注意ρj>0),这表明了Hr以及每个Hrj都是正定的,并且他们的逆也是正定的.

3 前言

我们的分析使用了下列假设

假设1

当1≤i≤N,函数fi:Rn→R 是凸函数,并且二阶可微.

假设2

存在常熟λ 和Λ, 使下式对于所有ω∈Rd和非空子集T⊆{1,…,N}.需要注意的是,在正则化的情况下,下界是平凡下界.

(40)λI2fT(w)ΛI

我们通过对我们的目标函数添加一个强凸的正则化,来使之保持强凸性.这些假设保证函数f有唯一最小值,我们用ω表示.

引理 3

假如假设1和假设2成立.记Br=Hr1.则有

tr(Br)(d+M)Λdet(Br)λ(M+d)/((d+M)Λ)M

我们在7.1节证明了引理3

引理 4

假如假设1和假设2成立.则存在常数0<γ≤Γ使对所有r≥1,Hr满足下列条件

(41)γIHrΓI

在7.2节,我们用如下特定值证明了引理4
γ=1(d+M)ΛandΓ=((d+M)Λ)d+M1λd+M.

我们将使用引理5, 这是一个很简单的强凸函数的结论.我们包含有完整的证明.

引理 5 假设函数f是连续可微并且λ强凸.记ω是函数f的唯一的最小值.那么对于x∈Rd,我们有

f(x)22λ(f(x)f(w)).

证明.对于强凸函数f有
f(w)f(x)+f(x)(wx)+λ2wx2f(x)+minv(f(x)v+λ2v2)=f(x)12λf(x)2.

最后一个等式当v=f(x)/λ.的时候达到最小值.
在引理6中,我们限定了我们方差约减后的梯度估计的范围.引理6的证明在7.3节中给出,和Johnson and Zhang (2013, Theorem 1)很接近.
引理6 记ω是函数f 的唯一的最小元.并记 vt=fS(xt)fS(wk)+μk, νt表示的是方差约减的随机梯度.基于条件Fk,t 并基于对S的期望,我们有
(42)Ek,t[vt2]4Λ(f(xt)f(w)+f(wk)f(w)).

4 收敛分析

定理7阐述了我们的主要结论.

定理7

假如假设1和假设2成立. 并记ω_*是函数f的唯一最小元.那么对于所有k ≥0,我们有

E[f(wk)f(w)]αkE[f(w0)f(w)]

其中收敛率α由下式给出
α=1/(2mη)+ηΓ2Λ2γληΓ2Λ2<1

假设我们选择η<γλ/(2Γ2Λ2),并选择m足够大使之满足下式
\begin{align}    \gamma \lambda & > \frac{1}{2m\eta} + 2\eta\Gamma^2\Lambda^2 . \label{eq:eta_m_size}  \end{align}

证明 基于假设2,使用∇f 的Lipschitz连续性,我们有

(33)f(xt+1)f(xt)+f(xt)(xt+1xt)+Λ2xt+1xt2=f(xt)ηf(xt)Hrvt+η2Λ2Hkvt2.

基于Fk,t的条件,并对式(9)期望化,然后就有

\begin{align} \label{eq:expansion_with_expectations}  & \,\,\mathbb E_{k,t}[f(x_{t+1})]  \\  \le & \,\, f(x_{t}) - \eta \nabla f(x_{t})^{\top} H_r \nabla f(x_{t}) + \frac{\eta^2 \Lambda}{2} \mathbb E_{k,t}\|H_{k}v_{t} \|^2 , \nonumber\end{align}

这里我们用到了Ek,t[νt]=f(xt). 我们用引理4来限定式(10)下中的第二项和第三项,并得到
Ek,t[f(xt+1)]f(xt)ηγf(xt)2+η2Γ2Λ2Ek,tvt2.

现在我们用引理6限定了Ek,t||νt||2.并且我们用引理5限定了||f(xt)||2.这样我们就有
Ek,t[f(xt+1)]f(xt)2ηγλ(f(xt)f(w))+2η2Γ2Λ2(f(xt)f(w)+f(wk)f(w))=f(xt)2η(γληΓ2Λ2)(f(xt)f(w))+2η2Γ2Λ2(f(wk)f(w)).

考虑所有随机变量的期望,并对所有t=0,…,m-1求和,利用裂项求和则有如下
E[f(xm)]E[f(x0)]+2mη2Γ2Λ2E[f(wk)f(w)]2η(γληΓ2Λ2)(t=0m1E[f(xt)]mf(w))=E[f(wk)]+2mη2Γ2Λ2E[f(wk)f(w)]2mη(γληΓ2Λ2)E[f(wk+1)f(w)].

重排上面的顺序,则有
0E[f(wk)f(xm)]+2mη2Γ2Λ2E[f(wk)f(w)]2mη(γληΓ2Λ2)E[f(wk+1)f(w)]E[f(wk)f(w)]+2mη2Γ2Λ2E[f(wk)f(w)]2mη(γληΓ2Λ2)E[f(wk+1)f(w)]=(1+2mη2Γ2Λ2)E[f(wk)f(w)]2mη(γληΓ2Λ2)E[f(wk+1)f(w)].

第二个不等式实际上利用了f(ω)f(xm).由于η<γλ/(2Γ2Λ2),可以得到
E[f(wk+1)f(w)]1+2mη2Γ2Λ22mη(γληΓ2Λ2)E[f(wk)f(w)].

由于我们选择m和η使得式(8)成立, 那么就可推出收敛率α是小于1的.这就完成了证明

5 相关工作

现在已经有大量工作试着通过方差约减来提升随机梯度下降. Shalev-Shwartz and Zhang (2013)提出了随机对偶梯度下降. Roux et al. (2012)提出了随机平均梯度方法(SAG). Johnson and Zhang (2013)提出方差约减梯度下降(SVRG). Wang et al.(2013) 开发了一种基于构造控制变量的方法.最近, Frostig et al. (2015)提出了在线版本的SVRG版本,他使用流估计的梯度来获得方差的减少.
同样的,一些随机拟牛顿法也被提出. Bordes et al.(2009)提出了一种能利用二阶信息的随机梯度下降的变种算法. Mokhtari and Ribeiro (2014a)分析了L-BFGS在随机下的直接应用,并在强凸条件下给予了O(1/k)的收敛性证明. Byrd et al. (2014) 提出了修改版本的随机L-BFGS,并在强凸条件下证明了O(1/k)的收敛率. Sohl-Dickstein et al. (2014) 提出了基于最小化函数和形式的拟牛顿法,他通过分别维护函数和中每个函数的Hessian逆的近似. Schraudolph et al. (2007) 开发了在线强凸下的随机L-BFGS. Wang et al.(2014) 证明了多种随机的拟牛顿法在非凸问题下的收敛率.我们的研究工作和上述不同在于,我们保证了线性的收敛率.
Lucchi et al.(2015) 独立提出了方差约减来加速随机拟牛顿法,并能够达到线性收敛率.他们更新Hessian逆的近似矩阵和L-BFGS很相似,然而我们的方法利用了Hessian-vector乘积来稳定了近似.

6实验结果

为了探究我们的理论结果,我们比较了算法1(SLFBFGS)算法和随机方差约减算法(SVRG) (Johnson and Zhang, 2013),随机拟牛顿法(SQN) (Byrd et al., 2014),以及随机梯度下降算法(SGD).我们在一些主流的机器学习模型上评估这些算法,包括岭回归,支持向量机,矩阵补全.我们的实验展示了算法在现实问题上的有效性,而这些问题不需要是(强)凸的.
由于SLBFGS和SVRG要求计算全体度,所以每个epoch要求额外遍历一边数据.除此之外,SLBFGS和SQN要求Hessian-vector计算,每一次Hessian-vector都相当于一次梯度计算Pearlmutter (1994).每个epoch需要计算Hessian-vector的次数引用中是(b_H N)\/(bL),而我们实验中要么是N要么是2N.为了能够合并这些额外的计算开销,我们的策略是展示基于遍历所有数据次数下的错误率(也就是计算梯度或者Hessian-vector次数除以N).由于这个原因,SLBFGS,SVRG,SQN,SGD第一次迭代的时间不同,SGD最先迭代而SLBFGS最后才开始迭代.
对于所有实验,我们设置批处理大小b的为20或者100,我们设置Hessian的批处理大小b_H为10b或者20b,设置Hessian更新间隔L为10,设置存储大小M为10,设置随机更新次数m为N/b.我们通过网格搜索来优化学习率.SLFBGS和SVRG使用了固定大小的步长.对于SQN和SGD,我们用了三种不同步长方案:固定值,1\/√t和1/t,然后展示最优的一个.矩阵填充问题为了打破对称性,我们用缩放了105的标准正态分布的随机变量来初始化实验,除此之外所有实验都初始化为0向量.
首先,我们展示在包含了4.6×105样本的百万歌曲数据集上的岭回归.我们初始化参数λ=103.在这个实验中,SLBFGS和SVRG都快速解决问题,达到高水平精度.接着,我们在RCV1(Lewis et al.,2004)上训练了支持向量机,该样本大概有7.8×106个样本.我们令正则化系数λ=0.在这个实验里,SGD和SQN如预期在开始获得了很好的进展,但是SLBFGS找到了更优解.然后,我们在Netfix Prize
数据集上求解了非凸的矩阵补全问题,如Recht and Re (2013)所阐述,大概有108个样本.我们让正则化系数λ=104.SVRG和SGD在这个问题上表现出很差的性能,可能是因为算法初始化为全0向量,而这正好是一个马鞍点(尽管不是最优).可能使用曲率信息帮助SLBFGS和SQN对于SVRG和SGD更快的逃离全0向量的领域.
figure
图1绘制了三种方法在三种问题上的比较.对于凸问题,我们绘制出相对于预先计算好的解的优化误差的对数.对于非凸问题,我们由于不知道全局最优解,所以简单的绘制出目标函数的值.

6.1 对步长选择的鲁棒性

fifure2
在这一节,我们将解释SLBFGS在凸问题大范围的步长下都能够表现优秀.而SVRG,SQN,SGD表现优秀下的步长范围则窄了很多.在图2,我们绘制了SLBFGS, SVRG,SQN,SGD在百万歌曲数据集上岭回归的性能,步长在几个数量级上变化.在图3,我们绘制了支持向量机在RCV1数据集与上面相似的图.在这两个例子中, SLBFGS表现良好,在大范围的步长下都能解决问题,并达到高水平精度,然而SVRG,SQN,SGD在最坏步长下性能下降很快.
figure3

7 前言的证明

7.1 引理3的证明

以下分析紧密遵循许多用在L-BFGS(Nocedal and Wright, 2006; Byrd et al., 2014;Mokhtari and Ribeiro, 2014a,b)中Hessian逆的近似的分析,我们包含了他证明的完整形式.
sjTyj=sj2f(Tj)(μj)sj,根据假设2有

(34)λsj2sjyjΛsj2.

类似的令zj=(2fTj(uj))1/2sj,那么就有

yj2sjyj=zj2fTj(uj)zjzjzj,

再一次应用假设2,推导可得
(35)λyj2sjyjΛ.

注意这里使用Sherman-Morrison-Woodbury定律,那么对于Hessian的近似B_r=H_r^(-1),我们就可以等效写出式子(4)如下
(36)Br(j)=Br(j1)Br(j1)sjsjBr(j1)sjBr(j1)sj+yjyjyjsj.

我们将开始限定B_r的特征值范围.我们将通过他的迹和行列式来间接限定特征值.我们可以得到
tr(Br(j))=tr(Br(j1))tr(Br(j1)sjsjBr(j1))sjBr(j1)sj+tr(yjyj)yjsj=tr(Br(j1))Br(j1)sj2sjBr(j1)sj+yj2yjsjtr(Br(j1))+yj2yjsjtr(Br(j1))+Λ.

第一个不等式可以从矩阵迹的线性得出.第二个不等式可以从tr(AB)=tr(AB)得出.第四个关系可以从式子(12)得出.既然
tr(Br(0))=dyr2sryrdΛ,

由上可以归纳出下列结论
tr(Bk)(d+M)Λ.

现在开始限定行列式的范围,我们得到
det(Br(j))=det(Br(j1))det(IsjsjBr(j1)sjBr(j1)sj+(Br(j1))1yjyjyjsj)=det(Br(j1))yjsjsjBr(j1)sj=det(Br(j1))yjsjsj2sj2sjBr(j1)sjdet(Br(j1))λλmax(Br(j1))det(Br(j1))λtr(Br(j1))det(Br(j1))λ(d+M)Λ.

第一个等式用到了det⁡(AB)=det⁡(A)det⁡(B).
第二个等式从下列恒等式得出
(37)det(I+u1v1+u2v2)=(1+u1v1)(1+u2v2)(u1v2)(v1u2)

通过令u1=sj, v1=(Br(j1)sj)/(sjBr(j1)sj), u2=(Br(j1))1yj, 以及 v2=yj/(yjsj). 具体证明参见Dennis and More (1977, Lemma 7.6), 或者可以把式子(14)看成是当I+μ1ν1可逆的时候det(A+μνT)=(1+νTA(1)μ)det(A)的两次应用.第三个等式是分子分母同时乘以||sj||2得到的.第四个关系式基于式子(11)并结合 sjBr(j1)sjλmax(Br(j1))sj2 得到.
第五个关系式用到了正定矩阵的最大特征值小于他的迹这样的性质.
第六个关系式用到之前tr(Brj1的上界.
既然有
det(Br(0))=(yr2sryr)dλd,

归纳推导可以得到
det(Br)λd+M((d+M)Λ)M.

7.2 引理4的证明
利用引理3和Hr是正定矩阵,我们可以得到
λmax(Br)tr(Br)(d+M)Λ.

以及
λmin(Br)det(Br)λmax(Br)d1λd+M((d+M)Λ)d+M1.

由于我们定义了Br=Hr1,可以得出
1(d+M)ΛIHr((d+M)Λ)d+M1λd+MI.

7.3 引理6的证明
定义函数gS(w)=fS(w)fS(w)fS(w)(ww) 来得到fS在最小元ω_*附近的线性近似,需要注意的是在ω这点时,gS达到最小值.对于任意ω我们有
0=gS(w)gS(w1ΛgS(w))gS(w)12ΛgS2.

整理上列式子,可以得到

fS(w)fS(w)22Λ(fS(w)fS(w)fS(w)(ww)).

对上面所有minibatch S⊆{1,…,N}以b为大小取平均,再用∇f(ω_* )=0,我们可以看到下列形式
(38)(Nb)1|S|=bfS(w)fS(w)22Λ(f(w)f(w)).

现在,令μk=f(wk) 以及vt=fS(xt)fS(wk)+μk. 在Fk,t的条件下对S取期望,我们可以发现

(39)Ek,t[vt2]2Ek,t[fS(xt)fS(w)2]+2Ek,t[fS(wk)fS(w)μk2]2Ek,t[fS(xt)fS(w)2]+2Ek,t[fS(wk)fS(w)2]4Λ(f(xt)f(w)+f(wk)f(w)).

第一个不等式用到了a+b22a2+2b2.第二个不等式,注意用到了μk=Ek,t[fS(wk)fS(w)]以及对于任意随机变量ξ有E[ξE[ξ]2]E[ξ2]第三个不等式来自式子(15).

8 探讨

这篇论文介绍了随机版本的L-BFGS,并提供了在强凸条件下的线性收敛率.理论7保证了SLBFGS算法线性收敛的质量,这在我们的实验结果中有体现出来. 当处于的糟糕的条件,曲率信息很有价值的时候,以及在解决高精度优化问题时,我们希望SLBFGS算法能够优胜于其他随机一阶优化方法.
在未来的工作中,有很多有趣的地需要解决. 定理7的证明和用于分析拟牛顿法的许多证明类似,导致常数随着问题大小缩小.在更深层次上,研究拟牛顿法的重点在于设计一些介于梯度下降和牛顿法之间的算法,以此能够获得梯度下降的计算优势以及牛顿法的快速收敛性.文献中的许多证明(包括定理7的证明),通过限制Hessain逆的近似偏离单位矩阵的程度,来限制了拟牛顿法偏离梯度下降的程度.然后用这些限制来保证准牛顿法不会比梯度下降差得多.未来的研究途径是研究是否可以设计随机拟牛顿法,如同在非随机情况下那样,证明其可证明超线性收敛.

译文原文出处:
Moritz, P., Nishihara, R., & Jordan, M. (2016, May). A linearly-convergent stochastic L-BFGS algorithm. In Artificial Intelligence and Statistics (pp. 249-258).

阅读更多
想对作者说点什么?