線性收斂的隨機L-BFGS算法
以下皆爲翻譯A Linearly-Convergent Stochastic L-BFGS Algorithm 原文鏈接
六級沒過,莫怪
概要
我們提出了一個新的隨機L-BFGS算法,並證明其對強凸連續函數可以達到線性收斂.我們的算法很大程度上基於最近在Byrd(2014)提出的L-BFGS算法的隨機變種,以及Johnson和Zhang最近所提出的對於隨機梯度下降的方差約減.我們用實驗證明了我們的算法在大規模凸優化以及非凸問題中表現良好,展現了線性收斂率以及能夠快速求出優化問題中的高精度解.此外,我們的算法能夠適用於廣泛的不同數量級的學習速率.
1 介紹
機器學習的一個趨勢是使用更多的參數來模擬更大的數據集。因此,爲這些大規模問題設計優化算法非常重要。這種情況下出現的典型優化問題是經驗風險最小化。形式如下:
minw1N∑i=1Nfi(w),(28)
其中ω∈
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 算法
我們考慮最小化下列函數的問題
f(w)=1N∑i=1Nfi(w)(29)
其中ω∈
Rd
.對於子集S⊂{1,2,…,N},我們定義下列子採樣目標函數
fS
fS(w)=1|S|∑i∈Sfi(w).(30)
我們的更新,將採用梯度
∇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估計.
在算法1中,我們用I來指代單位矩陣.我們用
Fk,t
表示當計數器k和t滿足一定條件時隨機變量所產生的sigma代數.形式如下
Fk,t=σ({Sk′,t′:k′<kork′=kandt′<t}∪{Tr:rL≤mk+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/sTjyTj
,然後當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}
初始化
H(r−M)r=(s⊤ryr/∥yr∥2)I
以及
Hr=H(r)r
注意上述等式(4)中的更新保留了正定性(注意
ρj
>0),這表明了
Hr
以及每個
Hjr
都是正定的,並且他們的逆也是正定的.
3 前言
我們的分析使用了下列假設
假設1
當1≤i≤N,函數
fi:Rn
→R 是凸函數,並且二階可微.
假設2
存在常熟λ 和Λ, 使下式對於所有ω∈
Rd
和非空子集T⊆{1,…,N}.需要注意的是,在正則化的情況下,下界是平凡下界.
λI⪯∇2fT(w)⪯ΛI(40)
我們通過對我們的目標函數添加一個強凸的正則化,來使之保持強凸性.這些假設保證函數f有唯一最小值,我們用
ω∗
表示.
引理 3
假如假設1和假設2成立.記
Br=H−1r.
則有
tr(Br)≤(d+M)Λdet(Br)≥λ(M+d)/((d+M)Λ)M
我們在7.1節證明了引理3
引理 4
假如假設1和假設2成立.則存在常數0<γ≤Γ使對所有r≥1,
Hr
滿足下列條件
γI⪯Hr⪯ΓI(41)
在7.2節,我們用如下特定值證明了引理4
γ=1(d+M)ΛandΓ=((d+M)Λ)d+M−1λd+M.
我們將使用引理5, 這是一個很簡單的強凸函數的結論.我們包含有完整的證明.
引理 5 假設函數f是連續可微並且λ強凸.記
ω∗
是函數f的唯一的最小值.那麼對於x∈
Rd
,我們有
∥∇f(x)∥2≥2λ(f(x)−f(w∗)).
證明.對於強凸函數f有
f(w∗)≥f(x)+∇f(x)⊤(w∗−x)+λ2∥w∗−x∥2≥f(x)+minv(∇f(x)⊤v+λ2∥v∥2)=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
的期望,我們有
Ek,t[∥vt∥2]≤4Λ(f(xt)−f(w∗)+f(wk)−f(w∗)).(42)
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連續性,我們有
≤=f(xt+1)f(xt)+∇f(xt)⊤(xt+1−xt)+Λ2∥xt+1−xt∥2f(xt)−η∇f(xt)⊤Hrvt+η2Λ2∥Hkvt∥2.(33)
基於
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,t∥vt∥2.
現在我們用引理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=0m−1E[f(xt)]−mf(w∗))E[f(wk)]+2mη2Γ2Λ2E[f(wk)−f(w∗)]−2mη(γλ−ηΓ2Λ2)E[f(wk+1)−f(w∗)].
重排上面的順序,則有
0≤≤=E[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,然後展示最優的一個.矩陣填充問題爲了打破對稱性,我們用縮放了
10−5
的標準正態分佈的隨機變量來初始化實驗,除此之外所有實驗都初始化爲0向量.
首先,我們展示在包含了
4.6×105
樣本的百萬歌曲數據集上的嶺迴歸.我們初始化參數
λ=10−3
.在這個實驗中,SLBFGS和SVRG都快速解決問題,達到高水平精度.接着,我們在RCV1(Lewis et al.,2004)上訓練了支持向量機,該樣本大概有
7.8×106
個樣本.我們令正則化係數λ=0.在這個實驗裏,SGD和SQN如預期在開始獲得了很好的進展,但是SLBFGS找到了更優解.然後,我們在Netfix Prize
數據集上求解了非凸的矩陣補全問題,如Recht and Re (2013)所闡述,大概有
108
個樣本.我們讓正則化係數
λ=10−4
.SVRG和SGD在這個問題上表現出很差的性能,可能是因爲算法初始化爲全0向量,而這正好是一個馬鞍點(儘管不是最優).可能使用曲率信息幫助SLBFGS和SQN對於SVRG和SGD更快的逃離全0向量的領域.
圖1繪製了三種方法在三種問題上的比較.對於凸問題,我們繪製出相對於預先計算好的解的優化誤差的對數.對於非凸問題,我們由於不知道全局最優解,所以簡單的繪製出目標函數的值.
6.1 對步長選擇的魯棒性
在這一節,我們將解釋SLBFGS在凸問題大範圍的步長下都能夠表現優秀.而SVRG,SQN,SGD表現優秀下的步長範圍則窄了很多.在圖2,我們繪製了SLBFGS, SVRG,SQN,SGD在百萬歌曲數據集上嶺迴歸的性能,步長在幾個數量級上變化.在圖3,我們繪製了支持向量機在RCV1數據集與上面相似的圖.在這兩個例子中, SLBFGS表現良好,在大範圍的步長下都能解決問題,並達到高水平精度,然而SVRG,SQN,SGD在最壞步長下性能下降很快.

7 前言的證明
7.1 引理3的證明
以下分析緊密遵循許多用在L-BFGS(Nocedal and Wright, 2006; Byrd et al., 2014;Mokhtari and Ribeiro, 2014a,b)中Hessian逆的近似的分析,我們包含了他證明的完整形式.
令
sTjyj=sj∇2f(Tj)(μj)sj
,根據假設2有
λ∥sj∥2≤s⊤jyj≤Λ∥sj∥2.(34)
類似的令
zj=(∇2fTj(uj))1/2sj
,那麼就有
∥yj∥2s⊤jyj=z⊤j∇2fTj(uj)zjz⊤jzj,
再一次應用假設2,推導可得
λ≤∥yj∥2s⊤jyj≤Λ.(35)
注意這裏使用Sherman-Morrison-Woodbury定律,那麼對於Hessian的近似B_r=H_r^(-1),我們就可以等效寫出式子(4)如下
B(j)r=B(j−1)r−B(j−1)rsjs⊤jB(j−1)rs⊤jB(j−1)rsj+yjy⊤jy⊤jsj.(36)
我們將開始限定B_r的特徵值範圍.我們將通過他的跡和行列式來間接限定特徵值.我們可以得到
tr(B(j)r)=tr(B(j−1)r)−tr(B(j−1)rsjs⊤jB(j−1)r)s⊤jB(j−1)rsj+tr(yjy⊤j)y⊤jsj=tr(B(j−1)r)−∥B(j−1)rsj∥2s⊤jB(j−1)rsj+∥yj∥2y⊤jsj≤tr(B(j−1)r)+∥yj∥2y⊤jsj≤tr(B(j−1)r)+Λ.
第一個不等式可以從矩陣跡的線性得出.第二個不等式可以從tr(AB)=tr(AB)得出.第四個關係可以從式子(12)得出.既然
tr(B(0)r)=d∥yr∥2s⊤ryr≤dΛ,
由上可以歸納出下列結論
tr(Bk)≤(d+M)Λ.
現在開始限定行列式的範圍,我們得到
det(B(j)r)=det(B(j−1)r)det⎛⎝I−sjs⊤jB(j−1)rs⊤jB(j−1)rsj+(B(j−1)r)−1yjy⊤jy⊤jsj⎞⎠=det(B(j−1)r)y⊤jsjs⊤jB(j−1)rsj=det(B(j−1)r)y⊤jsj∥sj∥2∥sj∥2s⊤jB(j−1)rsj≥det(B(j−1)r)λλmax(B(j−1)r)≥det(B(j−1)r)λtr(λλmax(B(j−1)r)≥det(B(j−1)r)λtr(B(j−1)r)≥det(B≥det(B(j−1)r)λtr(B(j−1)r)≥det(B(j−1)r)λ(d+M)Λ−1)r)