无约束最优化方法——牛顿法、拟牛顿法、BFGS、LBFGS

转自http://blog.youkuaiyun.com/lansatiankongxxc/article/details/45873597

这是前一段时间写的博客，然后又重新整理了一下

【最速下降法】

无约束最优化方法不涉及约束条件，所以都是介绍如何寻找搜索方向以及搜索步长。 
无约束最优化问题的目标函数： 

minx∈Rnf(x)

感觉这latex还是有些别扭，稍不留意就直接当做字符处理了。 
还是首先介绍一下梯度下降，梯度下降学过优化的都很清楚，一般叫最速下降法，这个方法有两点，首先是 x 更新的方向是负梯度方向，第二个是沿着该方向搜索，找到该方向的最小值所对应的 x 就是下次更新的值。梯度下降是最简单的一种方法，但是很多情况下却并不使用这种方法，原因是收敛速率比较慢，问题出在第二步上，由于搜索搜索时一直打到该方向的最小值，那么很显然，继续沿着该方向搜索会使函数值变小，函数梯度与搜索方向夹角大于九十度，所以该点的梯度和搜索方向在此时正交，这样相邻搜索点的梯度就会呈现锯齿状，函数沿着锯齿状下降，严重降低目标函数的收敛速率。 
梯度下降的递推公式推导是根据函数的一阶泰勒展开近似得到的。将 f(x) 在 x(k) 附近进行一阶泰勒展开： 

f(x)≈f(x(k))+gTk(x−x(k))

这里, gk=g(x(k))=∇f(x(k)) 为 f(x) 在 x(k) 的梯度。 
那么第 k+1 次的迭代值就可以通过： 

x(k+1)←x(k)+λkpk

. 
其中 pk 是搜索方向，取负梯度方向 pk=−∇f(x(k)) 可以使函数下降最快， λk 是步长，并且取 λk 使得 

f(x(k)+λkpk)=minλ≥0f(x(k)+λpk)

最速下降法就是这样，不断地寻找搜索方向以及确定搜索步长，直到达到终止条件，相邻函数值相遇某个阈值或是 x(k) 和 x(k+1) 小于某个阈值。但是产生的问题就是最速下降在接近终点的时候收敛速度较慢，容易之字形收敛。当然步长也不必是取该方向下降尽头的值，可以取固定值，但是太大容易发散，太小收敛速率比较慢。 
关于随机梯度下降法与批量下降法，大多数用梯度下降是求无约束目标函数，例如求经验损失最小时函数的参数，含有大量的训练数据。批量下降法是同时使用所有数据对梯度进行更新，很显然需要好多次迭代。随机梯度下降是每次只使用一个数据对函数参数进行更新，这样往往只通过一部分数据更新参数就会收敛，但是由于每次根据一个数据跟新，容易造成噪音问题。

【牛顿法】

由于最速梯度下降收敛速度并不“最速”，局部搜索中最速下降的方向和全局的最小值方向并不一致，所以就有后来改进的方法，包括牛顿法以及拟牛顿法。 
牛顿法要求 f(x) 具有二阶连续可导性。 
仍然考虑无约束最优化问题的目标函数： 

minx∈Rf(x)

这里所不同的是进行二阶泰勒展开： 

f(x)≈f(x(k))+gTk(x−x(k))+12(x−x(k))TH(x(k))(x−x(k))

这里, gk=g(x(k))=∇f(x(k)) 为 f(x) 在 x(k) 的梯度。 H(x(k)) 是 f(x)  的海塞矩阵 

H(x)=[∂2f(x)∂xi∂xj]n×n

显然， f(x) 有极值的条件是在 xk 处的一阶导数为0， ∇f(x)=0 ,所以，当我们从 xk 处开始搜索时，搜索终止处 xk+1 应该满足 ∇f(x(k+1))=0 。所以我们对二阶近似求导。 

∇f(x)=gk+Hk(x−x(k))

所以 

gk+Hk(x−x(k))=0

then, 

x(k+1)=x(k)−H−1kgk

经典牛顿法虽然具有二次收敛性，但是要求初始点需要尽量靠近极小点，否则有可能不收敛。计算过程中需要计算目标函数的二阶偏导数，难度较大。更为复杂的是目标函数的Hesse矩阵无法保持正定，会导致算法产生的方向不能保证是f在Xk 处的下降方向，从而令牛顿法失效；特别的，如果Hesse矩阵奇异，牛顿方向可能根本是不存在的。

拟牛顿法

上面说了，虽然牛顿法能够具有二次收敛性，但是要求太高，个别情况下甚至无法求出牛顿法的迭代方向，所以就有了拟牛顿法，来对Hesse矩阵的逆进行近似。 
通过泰勒二阶近似可以得到：

∇f(xk+1)=∇f(xk)+Hk(x(k+1)−xk)

令， 

yk=∇f(xk+1)−∇f(xk),sk=x(k+1)−xk

then, 

yk=Hksk

或者说，

H−1kyk=sk

注意到，

sk=x(k+1)−x(k)=αdk

,所以拟牛顿法模拟了牛顿的方向。 
所以，拟牛顿法选取满足条件 Bksk=yk , Bk 作为Hesse矩阵 Hk 的近似，或者 sk=Gkyk   Gk 作为hesse矩阵逆的近似，而且要使得计算简便。当有了 Bk 之后，通过对 Bk 进行低秩修改得到 Bk+1 , 

Bk+1=Bk+Δk

使其仍满足近似条件。 
一般，最初始 Bk 都是使用单位矩阵或者随机初始化。

SR1

根据修改 Bk 方法的不同，衍生出很多不同的方法，最简单的就是给 Bk−1 加上一个秩为1的对称矩阵，由于秩为1的对称矩阵可以写成一个列向量和其转置相乘的形式，所以 Bk 的约束条件可以写成： 

(Bk−1+βkukuTk)sk=yk

展开得到： 

Bk−1sk+βkukuTksk=yk

注意到 uTksk 是个常数，所以， 

−Bk−1sk+yk=(βkuTksk)uk

所以我们可以选 βk 使其满足 βkuTksk=1  

uk=yk−Bk−1sk,βk=1uTksk=1sTkuk=1sTk(yk−Bk−1sk)

最后得到 Bk 的更新式子 

Bk=Bk−1+(yk−Bk−1sk)(yk−Bk−1sk)TsTk(yk−Bk−1sk)

当然，通过 Gk 也能得到类似的式子，

BFGS

BFGS方法是一种秩2近似，至于为什么使用秩2近似这个暂时还不得而知。先讲一下是如何推导的。 
BFGS是近似海瑟矩阵 H ,首先，相应的牛顿条件是 

Bk+1sk=yk

使用秩2近似， 

Bk+1=Bk+Pk+Qk=Bk+αukuTk+βvkvTk

所以， 

Bk+1sk=(Bk+Pk+Qk)sk=Bksk+αukuTksk+βvkvTksk=yk

Bk+1sk=Bksk+(αuTksk)uk+(βvTksk)vk=yk

由于满足条件的 α,β,uk,vk 相当多，所以可以这样设置， 

αuTksk=1,βvTksk=1

α=1uTksk,β=1vTksk

这样式子就成了 

Bk+1sk=Bksk+uk+vk=yk

令 uk=yk,Bksk+vk=0,vk=−Bksk  
所以( Bk 是对称的) 

Bk=Bk+αukuTk+βvkvTk

=Bk+ykyTkyTksk−BksksTkBksTkBksk

我们使用的 Bk 的逆，所以这里还需要使用Sherman-Morrison公式，假设A是n阶可逆矩阵， u,v 是n维向量，且 A+uvT 也是可逆矩阵，则 

(A+uvT)−1=A−1−A−1uvTA−11+vTA−1u

得到 

B−1k+1=(I−skyTkyTksk)B−1k(I−yksTkyTksk)+sksTkyTksk

或者说使用Sherman–Morrison–Woodbury formula 进行一步变换【7】 

(A+UVT)−1=A−1−A−1U(I+VTA−1U)−1VTA−1

由这个式子就很容推了，上面式子可以写成 

(A+∑i=1kuivTi)−1=A−1−A−1U(C)−1VTA−1

Cij=δij+vTiA−1uj.i,j=1,2,...k

很明显，BFGS是对于Sherman–Morrison–Woodbury k=2的情况， 
我们可以令 

u1=v1=yk(sTkyk)1/2,u2=−v2=Bksk(sTkBksk)1/2

我们可以令 Hk=B−1k  

C11=1+vT1A−1u1=1+yTkHkyksTkyk

C22=1+vT2A−1u2=1−sTkBkHkBksksTkBksk=1−1=0

C12=vT1u2=yTkBksk(sTkyk)1/2)(sTkBksk)1/2=(sTkBksk)1/2(sTkyk)1/2

C21=vT2u1=−C11

回想一下2x2矩阵的逆 

C={β−αα0}

C−1=1α2{0α−αβ}

β=C11=1+yTkHkyksTkyk,α=C12=(sTkBksk)1/2(sTkyk)1/2

然后就是代入了,可以令 

U˜=HkU,V˜=HkV

这样，对于每一维 

u˜i=Hkui,v˜i=Hkvi,i=1,2

Hk+1←Hk−HkUC−1VTHk=Hk−U˜C−1V˜T=Hk+1α(−u˜1v˜T1)−βα2u˜2vT2=Hk−HkyksTk+skyTkHksTkyk+sksTksTkyk(1+yTkHkyksTkyk)

整理就得到BFGS的一般式了

DFP

DFP推导方法和BFGS类似，只不过是对hesse矩阵的逆进行近似,略。

LBFGS

关于LBFGS的推导，可以参考【3】和【4】，主要是通过BFGS的最后目标式子，不再保留完整的矩阵B_k^{-1}，因为当维度很大的时候（n>10^4），需要的空间非常大，所以保留了一些计算 B−1k 需要的 sk,yk 序列，而且只保存最近的m个序列。 
这里不妨用 Hk 表示 B−1k ,非hesse矩阵. 

Hk+1=(I−skyTkyTksk)Hk(I−yksTkyTksk)+sksTkyTksk

define: 
ρk=1yTksk , Vk=I−ρkyksTk ,then the above formulation can be rewritten as: 

Hk+1=VTkHkVk+ρksksTk

Then,recursively 

H1=VT0H0V0+ρ0s0sT0

H2===VT1H1V1+ρ1s1sT1VT1(VT0H0V0+ρ0s0sT0)V1+ρ1s1sT1VT1VT0H0V0V1+VT1ρ0s0sT0)V1+ρ1s1sT1

所以就有了这个公式：

Hk+1=++++(VTkVTk−1...VT1VT0)H0(V0V1...Vk−1Vk)(VTkVTk−1...VT1)ρ1s1sT1(V1...Vk−1Vk)...(VTk)ρk−1sk−1sTk−1(Vk)ρksksTk

然后为了算这个式子，需要不断迭代LBFGS原著中给了一个两层的递推程序求这个式子，只保留最近m步：

Hk+1=++...++(VTkVTk−1...VTk−m)H0(Vk−m...Vk−1Vk)(VTkVTk−1...VTk−m+1)ρk−msk−msTk−m(Vk−m+1...Vk−1Vk)(VTk)ρk−1sk−1sTk−1(Vk)ρksksTk

更新的方向：

Hk+1∇f(x)=++++(VTkVTk−1...VTk−m)H0(Vk−m...Vk−1Vk)∇f(x)(VTkVTk−1...VTk−m+1)ρk−msk−msTk−m(Vk−m+1...Vk−1Vk)∇f(x)...(VTk)ρk−1sk−1sTk−1(Vk)∇f(x)ρksksTk∇f(x)

所谓的Two-loop算法：

qk←∇f(xk)  
对 
i=k−1  to  k−m  
αi=ρisTiqi+1  
qi=qi+1−αiyi  
然后第二次循环， 
根据 wiki LBFGS 【5】 
H0=yTk−1sk−1yTk−1yk−1  
初始化： rk−m−1=H0qk−m  
对于  i=k−m,k−m+1  to  k−1  
βi=ρiyTiri−1  
ri=ri−1+si(αi−βi)  
最后得到的 r 即为所求。上面的q以及 r都只有最后一步结果，中间结果的可以用一个变量代替。

参考： 
【1】http://blog.youkuaiyun.com/lilyth_lilyth/article/details/8973972 
【2】统计学习方法 
【3】http://blog.youkuaiyun.com/lansatiankongxxc/article/details/38801863 
【4】http://blog.youkuaiyun.com/zhirom/article/details/38332111 
【5】http://en.wikipedia.org/wiki/Limited-memory_BFGS 
【6】http://en.wikipedia.org/wiki/Woodbury_matrix_identity 
【7】http://www.ing.unitn.it/~bertolaz/2-teaching/2004-2005/AA-2004-2005-PHD/lucidi/slides-mQN-1x2.pdf 
【8】http://www.iaeng.org/publication/WCE2012/WCE2012_pp1-5.pdf