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

转自http://blog.csdn.net/lansatiankongxxc/article/details/45873597

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

【最速下降法】

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

minxRnf(x)

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

minxRf(x)

这里所不同的是进行二阶泰勒展开: 
f(x)f(x(k))+gTk(xx(k))+12(xx(k))TH(x(k))(xx(k))

这里, gk=g(x(k))=f(x(k)) f(x) x(k) 的梯度。 H(x(k)) f(x)  的海塞矩阵 
H(x)=[2f(x)xixj]n×n

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

所以 
gk+Hk(xx(k))=0

then, 
x(k+1)=x(k)H1kgk

经典牛顿法虽然具有二次收敛性,但是要求初始点需要尽量靠近极小点,否则有可能不收敛。计算过程中需要计算目标函数的二阶偏导数,难度较大。更为复杂的是目标函数的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

或者说,
H1kyk=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 方法的不同,衍生出很多不同的方法,最简单的就是给 Bk1 加上一个秩为1的对称矩阵,由于秩为1的对称矩阵可以写成一个列向量和其转置相乘的形式,所以 Bk 的约束条件可以写成: 

(Bk1+βkukuTk)sk=yk

展开得到: 
Bk1sk+βkukuTksk=yk

注意到 uTksk 是个常数,所以, 
Bk1sk+yk=(βkuTksk)uk

所以我们可以选 βk 使其满足 βkuTksk=1  
uk=ykBk1sk,βk=1uTksk=1sTkuk=1sTk(ykBk1sk)

最后得到 Bk 的更新式子 
Bk=Bk1+(ykBk1sk)(ykBk1sk)TsTk(ykBk1sk)

当然,通过 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+ykyTkyTkskBksksTkBksTkBksk

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

得到 
B1k+1=(IskyTkyTksk)B1k(IyksTkyTksk)+sksTkyTksk

或者说使用Sherman–Morrison–Woodbury formula 进行一步变换【7】 
(A+UVT)1=A1A1U(I+VTA1U)1VTA1

由这个式子就很容推了,上面式子可以写成 
(A+i=1kuivTi)1=A1A1U(C)1VTA1

Cij=δij+vTiA1uj.i,j=1,2,...k

很明显,BFGS是对于Sherman–Morrison–Woodbury k=2的情况, 
我们可以令 
u1=v1=yk(sTkyk)1/2,u2=v2=Bksk(sTkBksk)1/2

我们可以令 Hk=B1k  
C11=1+vT1A1u1=1+yTkHkyksTkyk

C22=1+vT2A1u2=1sTkBkHkBksksTkBksk=11=0

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

C21=vT2u1=C11

回想一下2x2矩阵的逆 
C={βαα0}

C1=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+1HkHkUC1VTHk=HkU˜C1V˜T=Hk+1α(u˜1v˜T1)βα2u˜2vT2=HkHkyksTk+skyTkHksTkyk+sksTksTkyk(1+yTkHkyksTkyk)

整理就得到BFGS的一般式了

DFP

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

LBFGS

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

Hk+1=(IskyTkyTksk)Hk(IyksTkyTksk)+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=++++(VTkVTk1...VT1VT0)H0(V0V1...Vk1Vk)(VTkVTk1...VT1)ρ1s1sT1(V1...Vk1Vk)...(VTk)ρk1sk1sTk1(Vk)ρksksTk

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

Hk+1=++...++(VTkVTk1...VTkm)H0(Vkm...Vk1Vk)(VTkVTk1...VTkm+1)ρkmskmsTkm(Vkm+1...Vk1Vk)(VTk)ρk1sk1sTk1(Vk)ρksksTk

更新的方向:

Hk+1f(x)=++++(VTkVTk1...VTkm)H0(Vkm...Vk1Vk)f(x)(VTkVTk1...VTkm+1)ρkmskmsTkm(Vkm+1...Vk1Vk)f(x)...(VTk)ρk1sk1sTk1(Vk)f(x)ρksksTkf(x)

所谓的Two-loop算法:

qkf(xk)  
对 
i=k1  to  km  
αi=ρisTiqi+1  
qi=qi+1αiyi  
然后第二次循环, 
根据 wiki LBFGS 【5】 
H0=yTk1sk1yTk1yk1  
初始化: rkm1=H0qkm  
对于  i=km,km+1  to  k1  
βi=ρiyTiri1  
ri=ri1+si(αiβi)  
最后得到的 r 即为所求。上面的q以及 r都只有最后一步结果,中间结果的可以用一个变量代替。

参考: 
【1】http://blog.csdn.net/lilyth_lilyth/article/details/8973972 
【2】统计学习方法 
【3】http://blog.csdn.net/lansatiankongxxc/article/details/38801863 
【4】http://blog.csdn.net/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

  • 1
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值