Line Search and Quasi-Newton Methods 線性搜索與擬牛頓法

Gradient Descent

機器學習中很多模型的參數估計都要用到優化算法，梯度下降是其中最簡單也用得最多的優化算法之一。梯度下降(Gradient Descent)[3]也被稱之為最快梯度(Steepest Descent)，可用於尋找函數的局部最小值。梯度下降的思路為，函數值在梯度反方向下降是最快的，只要沿着函數的梯度反方向移動足夠小的距離到一個新的點，那么函數值必定是非遞增的，如圖1所示。

梯度下降思想的數學表述如下：

b=a−α∇F(a)⇒f(a)≥f(b)(1)(1)b=a−α∇F(a)⇒f(a)≥f(b)

其中f(x)f(x)為存在下界的可導函數。根據該思路，如果我們從x0x0為出發點，每次沿着當前函數梯度反方向移動一定距離αkαk，得到序列x0,x1,⋯,xnx0,x1,⋯,xn：

xk+1=xk−αk∇f(xk),0≤k≤n(2)(2)xk+1=xk−αk∇f(xk),0≤k≤n

對應的各點函數值序列之間的關系為：

f(x0)≥f(x1)≥f(x2)≥⋯≥f(xn)(3)(3)f(x0)≥f(x1)≥f(x2)≥⋯≥f(xn)

很顯然，當nn達到一定值時，函數f(x)f(x)是會收斂到局部最小值的。算法1簡單描述了一般化的梯度優化方法。在算法1中，我們需要選擇一個搜索方向dkdk滿足以下關系：

f(xk+αdk)<f(xk)∀α∈(0,ϵ](4)(4)f(xk+αdk)<f(xk)∀α∈(0,ϵ]

當dk=−∇f(x)dk=−∇f(x)時f(x)f(x)下降最快，但是只要滿足∇f(xk)Tdk<0∇f(xk)Tdk<0的dkdk都可以作為搜素方向。一般搜索方向表述為如下形式：

dk=−Bk∇f(xk)(5)(5)dk=−Bk∇f(xk)

其中BkBk為正定矩陣。當Bk=IBk=I時對應最快梯度下降算法；當Bk=H(xk)−1Bk=H(xk)−1時對應牛頓法，如果H(xk)=∇2f(xk)H(xk)=∇2f(xk)為正定矩陣。 在迭代過程中用於更新xkxk的步長αkαk可以是常數也可以是變化的。如果αkαk足夠小，收斂是可以得到保證的，但這意味這迭代次數nn要很大時函數才會收斂(圖2(a))；如果αkαk比較大，更新后的點很可能越過局部最優解(圖2(b))。有什么方法可以幫助我們自動確定最優步長呢？下面要說的線性搜索就包含一組解決方案。

 

Line Search

在給定搜索方向dkdk的前提下，線性搜索要解決的問題如下：

α=argminα≥0h(α)=argminα≥0f(xk+αdk)(6)(6)α=argminα≥0h(α)=argminα≥0f(xk+αdk)

如果h(α)h(α)是可微的凸函數，我們能通過解析解直接求得上式最優的步長；但非線性的優化問題需要通過迭代形式求得近似的最優步長。對於上式，局部或全局最優解對應的導數為h′(α)=∇f(xk+αdk)Tdk=0h′(α)=∇f(xk+αdk)Tdk=0。因為dkdk與f(xk)f(xk)在xkxk處的梯度方向夾角大於90度，因此h′(0)≤0h′(0)≤0，如果能找到α^α^使得h′(α^)>0h′(α^)>0，那么必定存在α⋆∈[0,α^)α⋆∈[0,α^)使得h′(α⋆)=0h′(α⋆)=0。有多種迭代算法可以求得α⋆α⋆的近似值，下面選擇幾種典型的介紹。

 

Bisection Search

二分線性搜索(Bisection Line Search)[2]可用於求解函數的根，其思想很簡單，就是不斷將現有區間划分為兩半，選擇必定含有使h′(α)=0h′(α)=0的半個區間作為下次迭代的區間，直到尋得h′(α⋆)≈0h′(α⋆)≈0為止，算法描述見2。二分線性搜素可以確保h(α)h(α)是收斂的，只要h(α)h(α)在區間(0,α^)(0,α^)上是連續的且h′(0)h′(0)和h(α^)h(α^)異號。經歷nn次迭代后，當前區間[αl,αh][αl,αh]的長度為：

L=(12)nα^(7)(7)L=(12)nα^

由迭代的終止條件之一αh−αl≥ϵαh−αl≥ϵ知迭代次數的上界為：

L≤ϵ⇒k≤[log2(α^ϵ)](8)(8)L≤ϵ⇒k≤[log2⁡(α^ϵ)]

下面給出二分搜索的Python代碼

 

 1 def bisection(dfun,theta,args,d,low,high,maxiter=1e4):
 2     """
 3     #Functionality:find the root of the function(fun) in the interval [low,high]
 4     #@Parameters
 5     #dfun:compute the graident of function f(x)
 6     #theta:Parameters of the model
 7     #args:other variables needed to compute the value of dfun
 8     #[low,high]:the interval which contains the root
 9     #maxiter:the max number of iterations
10     """
11     eps=1e-6
12     val_low=np.sum(dfun(theta+low*d,args)*d.T)
13     val_high=np.sum(dfun(theta+high*d,args)*d.T)
14     if val_low*val_high>0:
15         raise Exception('Invalid interval!')
16     iter_num=1
17     while iter_num<maxiter:
18         mid=(low+high)/2
19         val_mid=np.sum(dfun(theta+mid*d,args)*d.T)
20         if abs(val_mid)<eps or abs(high-low)<eps:
21             return mid
22         elif val_mid*val_low>0:
23             low=mid
24         else:
25             high=mid
26         iter_num+=1

Backtracking

回溯線性搜索(Backing Line Search)[1]基於Armijo准則計算搜素方向上的最大步長，其基本思想是沿着搜索方向移動一個較大的步長估計值，然后以迭代形式不斷縮減步長，直到該步長使得函數值f(xk+αdk)f(xk+αdk)相對與當前函數值f(xk)f(xk)的減小程度大於期望值(滿足Armijo准則)為止。Armijo准則(見圖3)的數學描述如下：

f(xk+αdk)≤f(xk)+c1αf′(xk)Tdk(9)(9)f(xk+αdk)≤f(xk)+c1αf′(xk)Tdk

其中f:Rn→Rf:Rn→R，c1∈(0,1)c1∈(0,1)，αα為步長，dk∈Rndk∈Rn為滿足f′(xk)Tdk<0f′(xk)Tdk<0的搜索方向。但是僅憑Armijo准則不足以求得較好的步長，根據前面的梯度下降的知識可知，只要αα足夠小就能滿足Armijo准則。因此常用的策略就是從較大的步長開始，然后以τ∈(0,1)τ∈(0,1)的速度縮短步長，直到滿足Armijo准則為止，這樣選出來的步長不至於太小，對應的算法描述見3。前面介紹的二分線性搜索的目標是求得滿足h′(α)≈0h′(α)≈0的最優步長近似值，而回溯線性搜索放松了對步長的約束，只要步長能使函數值有足夠大的變化即可。前者可以少計算幾次搜索方向，但在計算最優步長上花費了不少代價；后者退而求其次，找到一個差不多的步長即可，那么代價就是要多計算幾次搜索方向。接下來，我們要證明回溯線性搜索在Armijo准則下的收斂性問題[6]。因為h′(0)=f′(xk)Tdk<0h′(0)=f′(xk)Tdk<0，且0<c1<10<c1<1，則有

h′(0)<c1h′(0)<0(10)(10)h′(0)<c1h′(0)<0

根據導數的基本定義，結合上式，有如下關系：

h′(0)=limα→0h(α)−h(0)α=limα→0f(xk+αdk)−f(xk)α<ch′(0)(11)(11)h′(0)=limα→0h(α)−h(0)α=limα→0f(xk+αdk)−f(xk)α<ch′(0)

因此，存在一個步長α^>0α^>0，對任意的α∈(0,α^)α∈(0,α^),下式均成立

f(xk+αdk)−f(xk)α<cf′(xk)Tdk(12)(12)f(xk+αdk)−f(xk)α<cf′(xk)Tdk

即∀α∈(0,α^),f(xk+αdk)<f(xk)+cαf′(xk)Tdk∀α∈(0,α^),f(xk+αdk)<f(xk)+cαf′(xk)Tdk。 下面給出基於Armijo准則的線性搜索Python代碼：

 

 1 def ArmijoBacktrack(fun,dfun,theta,args,d,stepsize=1,tau=0.5,c1=1e-3):
 2     """
 3     #Functionality:find an acceptable stepsize via backtrack under Armijo rule
 4     #@Parameters
 5     #fun:compute the value of objective function
 6     #dfun:compute the gradient of objective function
 7     #theta:a vector of parameters of the model
 8     #stepsize:initial step size
 9     #c1:sufficient decrease Parameters
10     #tau:rate of shrink of stepsize
11     """
12     slope=np.sum(dfun(theta,args)*d.T)
13     obj_old=costFunction(theta,args)
14     theta_new=theta+stepsize*d
15     obj_new=costFunction(theta_new,args)
16     while obj_new>obj_old+c1*stepsize*slope:
17         stepsize*=tau
18         theta_new=theta+stepsize*d
19         obj_new=costFunction(theta_new,args)
20     return stepsize

Interpolation

基於Armijo准則的回溯線性搜索的收斂速度無法得到保證，特別是要回退很多次后才能落入滿足Armijo准則的區間。如果我們根據已有的函數值和導數信息，采用多項式插值法(Interpolation)[12,6,5,9]擬合函數，然后根據該多項式函數估計函數的極值點，這樣選擇合適步長的效率會高很多。 假設我們只有xkxk處的函數值f(xk)f(xk)及其倒數f′(xk)f′(xk)，且第一次嘗試的步長為α0α0。如果α0α0不滿足條件，那么我們根據這些信息可以構造一個二次近似函數hq(α)hq(α)

hq(α)=(h(α0)−h(0)−α0h′(0)α20)α2+h′(0)α+h(0)(13)(13)hq(α)=(h(α0)−h(0)−α0h′(0)α02)α2+h′(0)α+h(0)

注意，該二次函數滿足hq(0)=h(0)hq(0)=h(0)，h′q(0)=h′(0)hq′(0)=h′(0)和hq(α0)=h(α0)hq(α0)=h(α0)，如圖4(a)所示。接下來，根據hq(α)hq(α)的最小值估計下一個步長：

α1=h′(0)α202[h(0)+h′(0)α0−h(α0)](14)(14)α1=h′(0)α022[h(0)+h′(0)α0−h(α0)]

如果α1α1仍然不滿足條件，我們可以繼續重復上述過程，直到得到的步長滿足條件為止。假設我們在整個線性搜索過程中都用二次插值函數，那么最好有c1∈(0,0.5]c1∈(0,0.5]，為什么呢？簡單證明一下：如果α0α0不滿足Armijo准則，那么必定存在比α0α0小的步長滿足該准則，所以利用二次插值函數估算的步長α1<α0α1<α0才合理。結合α0α0不滿足Armijo准則和α1<α0α1<α0，可知c1≤0.5c1≤0.5。 如果我們已經嘗試了多個步長，卻每次只用上一次步長的相關信息構造二次函數，未免是對計算資源的浪費，其實我們可以利用多個步長的信息構造信息量更大更准確的插值函數的。在計算導數的代價大於計算函數值的代價時，應盡量避免計算h′(α)h′(α)，下面給出一個三次插值函數hc(α)hc(α)，如圖4(b)所示

hc(α)=aα3+bα2+h′(0)α+h(0)(15)(15)hc(α)=aα3+bα2+h′(0)α+h(0)

其中

[ab]=1α2i−1α2i(αi−αi−1)[α2i−1−α3i−1−α2iα3i][h(αi)−h(0)−h′(0)αih(αi−1)−h(0)−h′(0)αi−1](16)(16)[ab]=1αi−12αi2(αi−αi−1)[αi−12−αi2−αi−13αi3][h(αi)−h(0)−h′(0)αih(αi−1)−h(0)−h′(0)αi−1]

對hc(α)hc(α)求導，可得極值點αi+1∈[0,αi]αi+1∈[0,αi]的形式如下：

αi+1=−b+b2−3ah′(0)−−−−−−−−−−√3a(17)(17)αi+1=−b+b2−3ah′(0)3a

利用以上的三次插值函數求解下一個步長的過程不斷重復，直到步長滿足條件為止。如果出現a=0a=0的情況，三次插值函數退化為二次插值函數，在實現該算法時需要注意這點。在此過程中，如果αiαi太小或αi−1αi−1與αiαi太接近，需要重置αi=αi−1/2αi=αi−1/2，該保護措施(safeguards)保證下一次的步長不至於太小[6,5]。為什么會有這個作用呢？1）因為αi+1∈[0,αi]αi+1∈[0,αi]，所以當αiαi很小時αi+1αi+1也很小；2）當αi−1αi−1與αiαi太靠近時有a≈b≈∞a≈b≈∞，根據αi+1αi+1的表達式可知αi+1≈0αi+1≈0。 但是，在很多情況下，計算函數值后只需付出較小的代價就能順帶計算出導數值或其近似值，這使得我們可以用更精確的三次Hermite多項式[6]進行插值，如圖4(c)所示

H3(α)=[1+2αi−ααi−αi−1][α−αi−1αi−αi−1]2h(αi)+[1+2α−αi−1αi−αi−1][αi+1−ααi−αi−1]2h(αi−1)+(α−αi)[α−αi−1αi−αi−1]2h′(αi)+(α−αi−1)[αi−ααi−αi−1]2h′(αi−1)(18)(18)H3(α)=[1+2αi−ααi−αi−1][α−αi−1αi−αi−1]2h(αi)+[1+2α−αi−1αi−αi−1][αi+1−ααi−αi−1]2h(αi−1)+(α−αi)[α−αi−1αi−αi−1]2h′(αi)+(α−αi−1)[αi−ααi−αi−1]2h′(αi−1)

其中，三次Hermite多項式滿足H3(αi−1)=h(αi−1)H3(αi−1)=h(αi−1)，H3(αi)=h(αi)H3(αi)=h(αi)，H′3(αi−1)=h′(αi−1)H3′(αi−1)=h′(αi−1)和H′3(αi)=h′(αi)H3′(αi)=h′(αi)。H(α)H(α)的最小值只可能在兩側的端點或者極值點處，令H′3(α)=0H3′(α)=0可得下一個步長：

αi+1=αi−(αi−αi−1)[h′(αi)+d2−d1h′(αi)−h′(αi−1)+2d2](19)(19)αi+1=αi−(αi−αi−1)[h′(αi)+d2−d1h′(αi)−h′(αi−1)+2d2]

其中

d1=h′(αi)+h′(αi−1)−3[h(αi)−h(αi−1)αi−αi−1](20)(20)d1=h′(αi)+h′(αi−1)−3[h(αi)−h(αi−1)αi−αi−1]

d2=sign(αi−αi−1)d21−h′(αi−1)h′(αi)−−−−−−−−−−−−−−−√(21)(21)d2=sign(αi−αi−1)d12−h′(αi−1)h′(αi)

下面給出二次插值及三次插值的Python代碼：

 

 1 def quadraticInterpolation(a,h,h0,g0):
 2     """
 3     #Functionality:Approximate h(a) with a quadratic function and return its stationary point
 4     #@Parameters
 5     #a:current stepsize
 6     #h:a function value about stepsize,h(a)=f(x_k+a*d)
 7     #h:h(0)=f(x_k)
 8     #g0:h'(0)=f'(0)
 9     """
10     numerator=g0*a**2
11     denominator=2*(g0*a+h0-h)
12     if abs(denominator)<1e-12:#indicates that a is almost 0
13         return a
14     return numerator/denominator

def cubicInterpolation(a0,h0,a1,h1,h,g):
    """
    #Functionality:Approximate h(x) with a cubic function and return its stationary point
    #This version of cubic interpolation computes h'(x) as few as possible,suitable for the case in which computing derivative is more expensive than computing function values
    #@Parameters
    #a0 and a1 are stepsize it previous two iterations
    #h0:h(a0)
    #h1:h(a1)
    #h:h(0)=f(x)
    #g:h'(0)
    """
    mat=matlib.matrix([[a0**2,-a1**2],[-a0**3,a1**3]])
    vec=matlib.matrix([[h1-h-g*a1],[h0-h-g*a0]])
    ab=mat*vec/(a0**2*a1**2*(a1-a0))
    a=ab[0,0]
    b=ab[1,0]
    if abs(a)<1e-12:#a=0 and cubic function is a quadratic one
        return -g/(2*b)
    return (-b+np.sqrt(b**2-3*a*g))/(3*a)

def cubicInterpolationHermite(a0,h0,g0,a1,h1,g1):
    """
    #Functionality:Approximate h(a) with a cubic Hermite polynomial function and return its stationary point
    #This version of cubic interpolation computes h(a) as few as possible,suitable for the case in which computing derivative is easier than computing function values
    #@Parameters
    #a0 and a1 are stepsize it previous two iterations
    #h0:h(a0)
    #g0:h'(a0)
    #h1:h(a1)
    #g1:h'(a1)
    """
    d1=g0+g1-3*(h1-h0)/(a1-a0)
    d2=np.sign(a1-a0)*np.sqrt(d1**2-g0*g1)
    res=a1-(a1-a0)*(g1+d2-d2)/(g1-g0+2*d2)
    return res

基於Armijo准則的線性搜索的算法描述如下[4]對應的Armijo線性搜索的Python代碼如下：

 1 def ArmijoLineSearch(fun,dfun,theta,args,d,a0=1,c1=1e-3,a_min=1e-7,max_iter=1e5):
 2     """
 3     #Functionality:Line search under Armijo condition with quadratic and cubic interpolation
 4     #@Parameters
 5     #fun:objective Function
 6     #dfun:compute the gradient of fun
 7     #theta:a vector of parameters of the model
 8     #args:other variables needed for fun and func
 9     #d:search direction
10     #a0:initial stepsize
11     #c1:constant used in Armijo condition
12     #a_min:minimun of stepsize
13     #max_iter:maximum of the number of iterations
14     """
15     eps=1e-6
16     c1=min(c1,0.5)#c1 should<=0.5
17     a_pre=h_pre=g_pre=0
18     a_cur=a0
19     f_val=fun(theta,args) #h(0)=f(0)
20     g_val=np.sum(dfun(theta,args)*d.T) #h'(0)=f'(x)^Td
21     h_cur=g_cur=0
22     k=0
23     while a_cur>a_min and k<max_iter:
24         h_cur=fun(theta+a_cur*d,args)
25         g_cur=np.sum(dfun(theta+a_cur*d,args)*d.T)
26         if h_cur<=f_val+c1*a_cur*g_val: #meet Armijo condition
27             return a_cur
28         if not k: #k=0,use quadratic interpolation
29             a_new=quadraticInterpolation(a_cur,h_cur,f_val,g_val)
30         else: #k>0,use cubic Hermite interpolation
31             a_new=cubicInterpolationHermite(a_pre,h_pre,g_pre,a_cur,h_cur,g_cur)
32         if abs(a_new-a_cur)<eps or abs(a_new)<eps: #safeguard procedure
33             a_new=a_cur/2
34         a_pre=a_cur
35         a_cur=a_new
36         h_pre=h_cur
37         g_pre=g_cur
38         k+=1
39     return a_min #failed search

Wolfe Search

前面說到單憑Armijo准則(不考慮回溯策略)選出的步長可能太小，為了排除這些微小的步長，我們加上曲率的約束條件(如圖5所示)

h′(α)=f′(xk+αdk)Tdk≥c2f′(xk)Tdk(22)(22)h′(α)=f′(xk+αdk)Tdk≥c2f′(xk)Tdk

其中c2∈(c1,1)c2∈(c1,1)，c1c1為Armijo准則中的常量。當h′(α)h′(α)為很小的負數甚至為正數時，說明從xkxk沿着dkdk移動αα后的函數梯度方向與搜索方向的夾角接近90度，繼續向前移動已經不能很明顯減小函數值了，此時可以停止沿着dkdk繼續搜索；反之，說明繼續減小函數值的空間還是很大的，可以繼續向前搜索。Armijo准則與曲率約束兩者合起來稱為Wolfe准則[5]：

{f(xk+αdk)f′(xk+αdk)Tdk≤f(xk)+c1αf′(xk)Tdk≥c2f′(xk)Tdk(23)(23){f(xk+αdk)≤f(xk)+c1αf′(xk)Tdkf′(xk+αdk)Tdk≥c2f′(xk)Tdk

其中0<c1<c2<10<c1<c2<1。如圖6所示，滿足Wolfe准則的步長也許離h(α)h(α)的極值點較遠。我們可以修改曲率約束條件使得步長落到h(α)h(α)的極值點的一個較寬的領域中，強Wolfe准則對步長αα的約束如下：

{f(xk+αdk)|f′(xk+αdk)Tdk|≤f(xk)+c1αf′(xk)Tdk≥c2|f′(xk)Tdk|(24)(24){f(xk+αdk)≤f(xk)+c1αf′(xk)Tdk|f′(xk+αdk)Tdk|≥c2|f′(xk)Tdk|

強Wolfe准則不允許h′(α)h′(α)為太大的正數，可以排除遠離極值點的區間。那么到底是否存在滿足強Wolfe准則的步長呢？假設h(α)=f(xk+αdk)h(α)=f(xk+αdk)連續可微，在整個α>0α>0的定義域上存在下界。因為0<c1<10<c1<1，所以l(α)=f(xk)+αc1f′(xk)Tdkl(α)=f(xk)+αc1f′(xk)Tdk必然與h(α)h(α)至少有一個交點。假設α′α′為最小的交點對應的步長，則有

f(xk+α′dk)=f(xk)+α′c1f′(xk)Tdk(25)(25)f(xk+α′dk)=f(xk)+α′c1f′(xk)Tdk

那么對於滿足α∈(0,α′)α∈(0,α′)的步長必然都滿足Armijo准則。根據零值定理，存在α′′∈(0,α′)α″∈(0,α′)滿足

f(xk+α′dk)−f(xk)=α′f′(xk+α′′dk)Tdk(26)(26)f(xk+α′dk)−f(xk)=α′f′(xk+α″dk)Tdk

結合上面兩個關系式，由0<c1<c20<c1<c2和f′(xk)Tdk<0f′(xk)Tdk<0，可得

f′(xk+α′′dk)Tdk=c1f′(xk)Tdk>c2f′(xk)Tdk(27)(27)f′(xk+α″dk)Tdk=c1f′(xk)Tdk>c2f′(xk)Tdk

由此可知，α′′α″滿足強Wolfe准則。如果h(α)h(α)是一個較為平滑的函數，那么包含α′′α″的較小領域都會滿足強Wolfe准則。 如果在線性搜索過程中利用強Wolfe准則，可以更精確得找到更靠近極值點的步長，在目前線性搜索中用得很多。基於強Wolfe准則的線性搜索包含兩個階段：第一個階段從初始步長開始，不斷增大步長，直到找到一個滿足強Wolfe准則的步長或包含該步長的區間為止；第二個階段是在已知包含滿足強Wolfe准則步長區間的基礎上，不斷縮減區間，直到找到滿足強Wolfe准則的步長為止。基於強Wolfe准則的線性搜索算法描述如下5

 

在算法5中，αnewαnew的更新是因為在區間(αpre,αcur)(αpre,αcur)內沒有滿足Wolfe准則的步長，所以要選取一個大於αcurαcur的步長αnewαnew。在算法中，我們是用二次插值函數計算αnewαnew的，所以要求0<c1<0.50<c1<0.5。當然，也可以用其他方法，比如讓αcurαcur乘以一個大於1的常數，只要能較快找到一個包含滿足Wolfe的區間即可。所以，該算法每次嘗試的步長αcurαcur的逐漸遞增的；一旦找到了包含滿足Wolfe准則的步長的區間，立即調用zoomzoom函數不斷縮短區間，並返回滿足Wolfe的步長。根據算法邏輯，我們可以推斷出αpreαpre滿足Armijo准則，但違背曲率約束，而且導數為負數。由上述三個條件，可知αpreαpre必定位於滿足Wolfe准則的區間的左側的呈下降趨勢的曲線上，只要αcurαcur位於該區間的右側即可。那么怎樣判斷區間(αpre,αcur)(αpre,αcur)是否包含滿足Wolfe准則的步長呢？下面給出三種αcurαcur位於該區間的右側的充分條件：

1. αcurαcur不滿足Armijo准則；
2. h(αcur)≥h(αpre)h(αcur)≥h(αpre);
3. h′(αcur)≥0h′(αcur)≥0

這一點結合圖7就很容易理解了，我在圖中分別用紅色和綠色點標注了αpreαpre和αcurαcur可能的位置，藍色帶數字的圓圈注明了αcurαcur滿足哪些條件。基於Wolfe准則的線性搜索Python代碼如下：

 1 def WolfeLineSearch(fun,dfun,theta,args,d,a0=1,c1=1e-4,c2=0.9,a_min=1e-7,max_iter=1e5):
 2     """
 3     #Functionality:find a stepsize meeting Wolfe condition
 4     #@Parameters
 5     #fun:objective Function
 6     #dfun:compute the gradient of fun
 7     #theta:a vector of parameters of the model
 8     #args:other variables needed for fun and func
 9     #d:search direction
10     #a0:intial stepsize
11     #c1:constant used in Armijo condition
12     #c2:constant used in curvature condition
13     #a_min:minimun of stepsize
14     #max_iter:maximum of the number of iterations
15     """
16     eps=1e-16
17     c1=min(c1,0.5)
18     a_pre=0
19     a_cur=a0
20     f_val=fun(theta,args) #h(0)=f(x)
21     g_val=np.sum(dfun(theta,args)*d.T)
22     h_pre=f_val #h'(0)=f'(x)^Td
23     k=0
24     while k<max_iter and abs(a_cur-a_pre)>=eps:
25         h_cur=fun(theta+a_cur*d,args) #f(x+ad)
26         if h_cur>f_val+c1*a_cur*g_val or h_cur>=h_pre and k>0:
27             return zoom(fun,dfun,theta,args,d,a_pre,a_cur,c1,c2)
28         g_cur=np.sum(dfun(theta+a_cur*d,args)*d.T)
29         if abs(g_cur)<=-c2*g_val:#satisfy Wolfe condition
30             return a_cur
31         if g_cur>=0:
32             return zoom(fun,dfun,theta,args,d,a_pre,a_cur,c1,c2)
33         a_new=quadraticInterpolation(a_cur,h_cur,f_val,g_val)
34         a_pre=a_cur
35         a_cur=a_new
36         h_pre=h_cur
37         k+=1
38     return a_min

zoom函數的算法描述見6。zoom函數中需要傳入搜尋區間[αlow,αhigh][αlow,αhigh],其中αlow<αhighαlow<αhigh。本文中的zoom函數與文獻[5]中的內容略有差異，但是本文的zoom函數思路更簡單和清晰。由算法5中分析得到的調用zoom函數的條件，知道αlowαlow必須滿足Armijo准則，且位於所有滿足Wolfe准則的步長的左側。我們先取[αlow,αhigh][αlow,αhigh]區間的中值作為下一個測試的步長αnewαnew，如果恰好滿足Wolfe准則，則直接返回。如果αnewαnew違反Armijo准則或大於h(αlow)h(αlow)，顯然區間[αlow,αnew][αlow,αnew]包含滿足Wolfe准則的步長，因此用αnewαnew替換αhighαhigh以縮短區間長度；否則，αnewαnew必然也滿足Armijo准則，如果h′(αnew)>0h′(αnew)>0，則αnewαnew與αhighαhigh都在滿足Wolfe准則的區間右側，用αnewαnew替代αhighαhigh，反之則用αnewαnew替代αlowαlow。上述的迭代過程不斷縮短步長，知道求得滿足Wolfe准則的步長為止。如果在有限迭代次數內搜索失敗，則返回必然滿足Armijo准則的步長αlowαlow。

zoom函數對應的Python代碼如下：

 1 def zoom(fun,dfun,theta,args,d,a_low,a_high,c1=1e-3,c2=0.9,max_iter=1e4):
 2     """
 3     #Functionality:enlarge the interval to find a stepsize meeting Wolfe condition
 4     #@Parameters
 5     #fun:objective Function
 6     #dfun:compute the gradient of fun
 7     #theta:a vector of parameters of the model
 8     #args:other variables needed for fun and func
 9     #d:search direction
10     #[a_low,a_high]:interval containing a stepsize satisfying Wolfe condition
11     #c1:constant used in Armijo condition
12     #c2:constant used in curvature condition
13     #max_iter:maximum of the number of iterations
14     """
15     if a_low>a_high:
16         print('low:%f,high:%f'%(a_low,a_high))
17         raise Exception('Invalid interval of stepsize in zoom procedure')
18     eps=1e-16
19     h=fun(theta,args) #h(0)=f(x)
20     g=np.sum(dfun(theta,args)*d.T) #h'(0)=f'(x)^Td
21     k=0
22     h_low=fun(theta+a_low*d,args)
23     h_high=fun(theta+a_high*d,args)
24     if h_low>h+c1*a_low*g:
25         raise Exception('Left endpoint violates Armijo condition in zoom procedure')
26     while k<max_iter and abs(a_high-a_low)>=eps:
27         a_new=(a_low+a_high)/2
28         h_new=fun(theta+a_new*d,args)
29         if h_new>h+c1*a_new*g or h_new>h_low:
30             a_high=a_new
31             h_high=h_new
32         else:
33             g_new=np.sum(dfun(theta+a_new*d,args)*d.T)
34             if abs(g_new)<=-c2*g: #satisfy Wolfe condition
35                 return a_new 
36             if g_new*(a_high-a_low)>=0:
37                 a_high=a_new
38                 h_high=h_new
39             else:
40                 a_low=a_new
41                 h_low=h_new
42         k+=1
43     return a_low #a_low definitely satisfy Armijo condition

Newton's Method

牛頓法(Newton's method)[8]以迭代方式求解函數的根，其基本思想是從一個初始點出發，不斷在當前點xkxk處用切線近似函數f(x)f(x)，並求得該切線與xx軸的交點作為下一次的迭代初始點xk+1xk+1，直到找到f(x)=0f(x)=0的近似解為止。Newton法可用於二次可微函數f(x)f(x)的最優化問題。 在xkxk處用二階泰勒展開來對f(xk)f(xk)其進行逼近。

f(xk+△x)≈f(xk)+f′(xk)△x+12△xTBk△x(28)(28)f(xk+△x)≈f(xk)+f′(xk)△x+12△xTBk△x

現在，我們的目標是在xkxk附近求得使f(x)f(x)取得極小值的△x△x。將上式對△x△x求導可得函數f′(x)f′(x)在xk+1=xk+△xxk+1=xk+△x處的線性近似如下：

f′(xk+1)=f′(xk)+Bk(xk+1−xk)(29)(29)f′(xk+1)=f′(xk)+Bk(xk+1−xk)

其中Bk=∇2f(xk)Bk=∇2f(xk)為f(x)f(x)在xkxk處對應的Hessian矩陣。由於函數的極值點一般都對應f′(x)=0f′(x)=0，令f′(xk+1)=0f′(xk+1)=0並化簡可得迭代公式為：

xk+1=xk−B−1kf′(xk)(30)(30)xk+1=xk−Bk−1f′(xk)

牛頓迭代法收斂速度很快，對於二次函數可以一次性找到最優解。但用於求解優化問題時，需要付出很大的代價求得函數的一階導數、二階導數及其逆矩陣。此外，有的函數還存在不可導、Hessian矩陣不可逆、迭代點之間存在循環(即xk+t=xkxk+t=xk)等情形，這些都成為了牛頓法的應用障礙。牛頓迭代法用於求解極值點f′(x)=0f′(x)=0的步驟見算法7。當然，也可用牛頓法求最優步長，只需將算法7中的函數f(x)f(x)替換為關於步長的函數h(α)h(α)即可。

 

Quasi-Newton Method

擬牛頓(Quasi-Newton)[11]算法可用於求解函數的局部最優解，也就是那些導數為0的駐點。牛頓法用於解決優化問題時，事先假設原函數可用二次函數近似，然后用一階和二階導數尋找局部最優解。而在擬牛頓算法中，不需要准確計算Hessian矩陣，取而代之的是運用下面的擬牛頓條件分析連續兩個梯度向量得到的近似值矩陣BkBk：

f′(xk+1)−f′(xk)≈Bk+1(xk+1−xk)(31)(31)f′(xk+1)−f′(xk)≈Bk+1(xk+1−xk)

擬牛頓算法的算法流程見8，其基本思想是利用矩陣BkBk計算牛頓方向的近似值dkdk。各種擬牛頓算法的主要差異在於近似Hessian矩陣的更新策略，表1列出了部分主流的擬牛頓算法的迭代更新規則，其中sk=xk+1−xk=−αkB−1kf′(xk)sk=xk+1−xk=−αkBk−1f′(xk),yk=f′(xk+1)−f′(xk)yk=f′(xk+1)−f′(xk)。擬牛頓算法中最常用的是BFGS，其針對有限內存的機器的算法變種L-BFGS[4]在機器學習領域又備受青睞。BFGS需要存儲n×nn×n的矩陣HkHk用於近似Hessian矩陣的逆矩陣；而L-BFGS僅需要存儲過去mm(mm一般小於10)對nn維的更新數據(x,f′(xi))(x,f′(xi))即可。L-BFGS的空間復雜度為線性，特別適用於變量非常多的優化問題。BFGS的算法描述很容易寫出來，如下：

 

 1 def BFGS(fun,dfun,theta,args,H=None,mode=0,eps=1e-12,max_iter=1e4):
 2     """
 3     #Functionality:find the minimum of objective function f(x)
 4     #@Parameters
 5     #fun:objective function f(x)
 6     #dfun:compute the gradient of f(x)
 7     #args:parameters needed by fun and dfun
 8     #theta:start vector of parameters of the model
 9     #H:initial inverse Hessian approximation
10     #mode:index of line search algorithm
11     """
12     x_pre=x_cur=theta
13     g=dfun(x_cur,args)
14     I=matlib.eye(theta.size)
15     if not H:#initialize H as an identity matrix
16         H=I
17     k=0
18     while k<max_iter and np.sum(np.abs(g))>eps:
19         d=-g*H
20         step=LineSearch(fun,dfun,x_pre,args,d,1,mode)
21         x_cur=x_pre+step*d
22         s=step*d
23         y=dfun(x_cur,args)-dfun(x_pre,args)
24         ys=np.sum(y*s.T)
25         if abs(ys)<eps:
26             return x_cur
27         change=(ys+np.sum(y*H*y.T))*(s.T*s)/(ys**2)-(H*y.T*s+s.T*y*H)/ys
28         H+=change
29         g=dfun(x_cur,args)
30         x_pre=x_cur
31         k+=1
32     return x_cur

下面我們分析如何構造下L-BFGS的算法[10,13]。假設我們現在處於優化過程的第k(k≥m)k(k≥m)次迭代，參數為xkxk，梯度gk=f′(xk)gk=f′(xk)，已經保存的mm條更新數據為sk=xk+1−xksk=xk+1−xk及yk=gk+1−gkyk=gk+1−gk。我們最終需要計算的是搜索方向dk=−Hkgkdk=−Hkgk，於是令Vk=(I−ρkyksTk)Vk=(I−ρkykskT)，ρk=1/(yTksk)ρk=1/(ykTsk)，將表1中BFGS的HkHk的更新規則展開，我們可以得到下式：

===HkgkVTk−1Hk−1Vk−1gk+sk−1ρk−1sTk−1gkVTk−1VTk−2Hk−2Vk−2Vk−1gk+VTk−1sk−2ρk−2sTk−2Vk−1gk+sk−1ρk−1sTk−1gk(VTk−1VTk−2⋯VTk−m)Hk−m(Vk−m⋯Vk−2Vk−1)gk+(VTk−1VTk−2⋯VTk−m+1)sk−mρk−msTk−m(Vk−m+1⋯Vk−1Vk)gk+(VTk−1VTk−2⋯VTk−m+2)sk−m+1ρk−m+1sTk−m+1(Vk−m+2⋯Vk−2Vk−1)gk+⋯+VTk−1sk−2ρk−2sTk−2Vk−1gk+sk−1ρk−1sTk−1gk(32)(32)Hkgk=Vk−1THk−1Vk−1gk+sk−1ρk−1sk−1Tgk=Vk−1TVk−2THk−2Vk−2Vk−1gk+Vk−1Tsk−2ρk−2sk−2TVk−1gk+sk−1ρk−1sk−1Tgk=(Vk−1TVk−2T⋯Vk−mT)Hk−m(Vk−m⋯Vk−2Vk−1)gk+(Vk−1TVk−2T⋯Vk−m+1T)sk−mρk−msk−mT(Vk−m+1⋯Vk−1Vk)gk+(Vk−1TVk−2T⋯Vk−m+2T)sk−m+1ρk−m+1sk−m+1T(Vk−m+2⋯Vk−2Vk−1)gk+⋯+Vk−1Tsk−2ρk−2sk−2TVk−1gk+sk−1ρk−1sk−1Tgk

上式非常有規律，這就為迭代求解奠定了很好的基礎。我們令q0=gkq0=gk,則當1≤i≤m1≤i≤m時有

qi=(Vk−i⋯Vk−2Vk−1)gk(33)(33)qi=(Vk−i⋯Vk−2Vk−1)gk

ai=ρk−isTk−iqi−1(34)(34)ai=ρk−isk−iTqi−1

那么可以得到如下的迭代規則：

qi=Vk−i+1qi−1=qi−1−ρk−i+1yk−i+1sTk−i+1qi−1=qi−1−ai−1yk−i+1(35)(35)qi=Vk−i+1qi−1=qi−1−ρk−i+1yk−i+1sk−i+1Tqi−1=qi−1−ai−1yk−i+1

到目前為止，我們已經可以求解出HkgkHkgk所有項的右半部分，那左半部分如何處理？在這里采用不斷提前最左端的公因式的方法完成迭代過程：

Hkgk=P1=VTk−1P2+sk−1a1(36)(36)Hkgk=P1=Vk−1TP2+sk−1a1

P2=VTk−2P3+sk−2a2(37)(37)P2=Vk−2TP3+sk−2a2

重復該過程，很快就可以發現規律：

Pi=VTk−iPi+1+sk−iai=Pi+1+sk−i(ai−ρk−iyTk−iPi+1)(38)(38)Pi=Vk−iTPi+1+sk−iai=Pi+1+sk−i(ai−ρk−iyk−iTPi+1)

其中Pm+1=Hk−mqmPm+1=Hk−mqm。 根據上述分析，我們可以得到L-BFGS的求解搜索方向的算法9。根據算法9的整個流程，可知通過兩個循環mm次的迭代運算即可出計算當前的搜索方向，需要存儲歷史數據{sk−i,yk−i|i=1,⋯,m}{sk−i,yk−i|i=1,⋯,m}和臨時數據{ak−i|i=1,⋯,m}{ak−i|i=1,⋯,m}，所以算法的時間和空間復雜度均為O(mn)O(mn)。如果目前處於迭代的初期，已有的歷史數據少於mm，那么就用這些已有的數據，在后續迭代過程中不斷新增歷史數據即可；若干當前的迭代次數不小於mm，那么在每次計算出搜索方向后，即可用sksk和ykyk替換sk−msk−m和yk−myk−m組成新的mm對歷史更新數據。

 

在算法9中，需要給出矩陣Hk−mHk−m。在第一次迭代時，Hk−mHk−m被初始化為單位陣，在隨后的迭代過程中Hk−m=γkIHk−m=γkI，其中

γk=yTk−1sk−1yTk−1yk−1(39)(39)γk=yk−1Tsk−1yk−1Tyk−1

另外，在內存受限的系統中存儲n×nn×n不是很現實的想法。用上述的方法，我們僅需存儲一個標量γkγk即可，這是一個簡單卻又高效的做法[13]。 最后，附上L-BFGS的Python版本代碼：

 

 1 def LBFGS(fun,dfun,theta,args,mode=0,eps=1e-12,max_iter=1e4):
 2     """
 3     #Functionality:find the minimum of objective function f(x) with LBFGS
 4     #@Parameters
 5     #fun:objective function f(x)
 6     #dfun:compute the gradient of f(x)
 7     #args:parameters needed by fun and dfun
 8     #theta:start vector of parameters of the model
 9     #H:initial inverse Hessian approximation
10     #mode:index of line search algorithm
11     """
12     x_pre=x_cur=theta
13     s_arr=[]
14     y_arr=[]
15     Hscale=1
16     k=0
17     while k<max_iter:
18         g=dfun(x_cur,args)
19         d=LBFGSSearchDirection(y_arr,s_arr,Hscale,-g)
20         step=LineSearch(fun,dfun,x_pre,args,d,1,mode)
21         s=step*d
22         x_cur=x_pre+s
23         y=dfun(x_cur,args)-dfun(x_pre,args)
24         ys=np.sum(y*s.T)
25         if np.sum(np.abs(s))<eps:
26             return x_cur
27         x_pre=x_cur
28         k+=1
29         y_arr,s_arr,Hscale=LBFGSUpdate(y,s,y_arr,s_arr)
30     return x_cur
31 
32     
33 def LBFGSSearchDirection(y_arr,s_arr,Hscale,g):
34     """
35     #Functionality:estimate search direction using with LBFGS
36     #@Parameters
37     #y_arr:m*dim matrix,where y_arr[i,:]=f'(x_{i+1})-f'(x_i)
38     #s_arr:m*dim matrix,where s_arr[i,:]=x_{k+1}-x_k
39     #Hscale:a scale to initilize the inverse of Hessian matrix
40     #g:a row vector representing -f'(x_{k})
41     """
42     histNum=len(s_arr)#number of update data stored
43     if not histNum:
44         return g
45     dim=s_arr[0].size
46     a_arr=[0 for i in range(histNum)]
47     rho=[0 for i in range(histNum)]
48     q=g
49     for i in range(1,histNum+1):
50         s=s_arr[histNum-i]
51         y=y_arr[histNum-i]
52         rho[histNum-i]=1/np.sum(s*y.T)
53         a_arr[i-1]=rho[histNum-i]*np.sum(s*q.T)
54         q-=(a_arr[i-1]*y)
55     P=Hscale*q
56     for i in range(histNum,0,-1):
57         y=y_arr[histNum-i]
58         s=s_arr[histNum-i]
59         beta=rho[histNum-i]*np.sum(y*P.T)
60         P+=s*(a_arr[i-1]-beta)
61     return P
62         
63 
64 def LBFGSUpdate(y,s,oldy,olds,m=1e2):
65     """
66     #Functionality:refresh the historical update data
67     #@Parameters
68     #y:f'(x_{k+1})-f'(x_k)
69     #s:x_{k+1}-x_k
70     #oldy:[y0,y1,...],which is a list
71     #olds:[s0,s1,...],which is a list
72     #m:number of historical data to store(default:100)
73     """
74     eps=1e-12
75     Hscale=np.sum(y*s.T/y*y.T) #a scale to initialize H_{k-m}
76     if Hscale<eps:#skip update
77         return oldy,olds,Hscale
78     
79     cur_m=len(oldy)
80     if cur_m>=m:
81         oldy.pop(0)
82         olds.pop(0)
83     oldy.append(copy.deepcopy(y))
84     olds.append(copy.deepcopy(s))
85     return oldy,olds,Hscale 

References

[1] Backtracking line search. http://en.wikipedia.org/wiki/Backtracking_line_search.

[2] Bisection method. http://en.wikipedia.org/wiki/Bisection_method.

[3] Gradient descent. http://en.wikipedia.org/wiki/Gradient_descent.

[4] Limited-memory bfgs. http://en.wikipedia.org/wiki/Limited-memory_BFGS.

[5] Line search methods. http://pages.cs.wisc.edu/~ferris/cs730/chap3.pdf.

[6] Line search methods:step length selection. http://terminus.sdsu.edu/SDSU/Math693a_f2013/Lectures/06/lecture.pdf.

[7] Math 408a line search methods. https://www.math.washington.edu/~burke/crs/408/lectures/L7-line-search.pdf.

[8] Newton’s method. http://en.wikipedia.org/wiki/Newton%27s_method.

[9] Nonlinear programming algorithms. http://www.math.bme.hu/~bog/GlobOpt/Chapter5.pdf.

[10] Oerview of quasi-newton optimization methods. https://homes.cs.washington.edu/~galen/files/quasi-newton-notes.pdf.

[11] Quasi-newton method. http://en.wikipedia.org/wiki/Quasi-Newton_method.

[12] Unconstrained minimization. http://www.ing.unitn.it/~bertolaz/2-teaching/2011-2012/AA-2011-2012-OPTIM/lezioni/slides-mND.pdf.

[13] Dong C Liu and Jorge Nocedal. On the limited memory bfgs method for large scale optimization. Mathematical programming, 45(1-3):503–528,1989.

  作者：JeromeWang 
  郵箱：[e-mail] 
  出處：http://www.cnblogs.com/jeromeblog/