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

来源：互联网发布：在淘宝买摩托车可靠吗编辑：程序博客网时间：2024/04/30 07:35

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

最速下降法
牛顿法
拟牛顿法
SR1
BFGS
DFP
LBFGS

【最速下降法】

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

min x \in R n f (x)

感觉这latex还是有些别扭，稍不留意就直接当做字符处理了。
还是首先介绍一下梯度下降，梯度下降学过优化的都很清楚，一般叫最速下降法，这个方法有两点，首先是

x更新的方向是负梯度方向，第二个是沿着该方向搜索，找到该方向的最小值所对应的

x就是下次更新的值。梯度下降是最简单的一种方法，但是很多情况下却并不使用这种方法，原因是收敛速率比较慢，问题出在第二步上，由于搜索搜索时一直打到该方向的最小值，那么很显然，继续沿着该方向搜索会使函数值变小，函数梯度与搜索方向夹角大于九十度，所以该点的梯度和搜索方向在此时正交，这样相邻搜索点的梯度就会呈现锯齿状，函数沿着锯齿状下降，严重降低目标函数的收敛速率。
梯度下降的递推公式推导是根据函数的一阶泰勒展开近似得到的。将

f(x)在

x(k)附近进行一阶泰勒展开：

f (x) \approx f (x (k)) + g T k (x - x (k))

这里,

gk=g(x(k))=∇f(x(k))为

f(x)在

x(k)的梯度。
那么第

k+1次的迭代值就可以通过：

x (k + 1) \leftarrow x (k) + λ k p k

.
其中

pk是搜索方向，取负梯度方向

pk=−∇f(x(k))可以使函数下降最快，

λk是步长，并且取

λk使得

f (x (k) + λ k p k) = min λ \geq 0 f (x (k) + λ p k)

最速下降法就是这样，不断地寻找搜索方向以及确定搜索步长，直到达到终止条件，相邻函数值相遇某个阈值或是

x(k)和

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

【牛顿法】

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

min x \in R f (x)

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

f (x) \approx f (x (k)) + g T k (x - x (k)) + 1 2 (x - x (k)) T H (x (k)) (x - x (k))

这里,

gk=g(x(k))=∇f(x(k))为

f(x)在

x(k)的梯度。

H(x(k))是

f(x) 的海塞矩阵

H (x) = [\partial 2 f ( x ) \partial x i \partial x j] n \times n

显然，

f(x)有极值的条件是在

xk处的一阶导数为0，

∇f(x)=0,所以，当我们从

xk处开始搜索时，搜索终止处

xk+1应该满足

∇f(x(k+1))=0。所以我们对二阶近似求导。

\nabla f (x) = g k + H k (x - x (k))

所以

g k + H k (x - x (k)) = 0

then,

x (k + 1) = x (k) - H - 1 k g k

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

拟牛顿法

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

\nabla f (x k + 1) = \nabla f (x k) + H k (x (k + 1) - x k)

令，

y k = \nabla f (x k + 1) - \nabla f (x k), s k = x (k + 1) - x k

then,

y k = H k s k

或者说，

H - 1 k y k = s k

注意到，

s k = x (k + 1) - x (k) = α d k

,所以拟牛顿法模拟了牛顿的方向。
所以，拟牛顿法选取满足条件

Bksk=yk,

Bk作为Hesse矩阵

Hk的近似，或者

sk=Gkyk

Gk作为hesse矩阵逆的近似，而且要使得计算简便。当有了

Bk之后，通过对

Bk进行低秩修改得到

Bk+1,

B k + 1 = B k + Δ k

使其仍满足近似条件。
一般，最初始

Bk都是使用单位矩阵或者随机初始化。

SR1

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

(B k - 1 + β k u k u T k) s k = y k

展开得到：

B k - 1 s k + β k u k u T k s k = y k

注意到

uTksk是个常数，所以，

- B k - 1 s k + y k = (β k u T k s k) u k

所以我们可以选

βk使其满足

βkuTksk=1

u k = y k - B k - 1 s k, β k = 1 u T k s k = 1 s T k u k = 1 s T k ( y k - B k - 1 s k )

最后得到

Bk的更新式子

B k = B k - 1 + ( y k - B k - 1 s k ) ( y k - B k - 1 s k ) T s T k ( y k - B k - 1 s k )

当然，通过

Gk也能得到类似的式子，

BFGS

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

B k + 1 s k = y k

使用秩2近似，

B k + 1 = B k + P k + Q k = B k + α u k u T k + β v k v T k

所以，

B k + 1 s k = (B k + P k + Q k) s k = B k s k + α u k u T k s k + β v k v T k s k = y k

B k + 1 s k = B k s k + (α u T k s k) u k + (β v T k s k) v k = y k

由于满足条件的

α,β,uk,vk相当多，所以可以这样设置，

α u T k s k = 1, β v T k s k = 1

α = 1 u T k s k, β = 1 v T k s k

这样式子就成了

B k + 1 s k = B k s k + u k + v k = y k

令

uk=yk,Bksk+vk=0,vk=−Bksk
所以(

Bk是对称的)

B k = B k + α u k u T k + β v k v T k

= B k + y k y T k y T k s k - B k s k s T k B k s T k B k s k

我们使用的

Bk的逆，所以这里还需要使用Sherman-Morrison公式，假设A是n阶可逆矩阵，

u,v是n维向量，且

A+uvT也是可逆矩阵，则

(A + u v T) - 1 = A - 1 - A - 1 u v T A - 1 1 + v T A - 1 u

得到

B - 1 k + 1 = (I - s k y T k y T k s k) B - 1 k (I - y k s T k y T k s k) + s k s T k y T k s k

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

(A + U V T) - 1 = A - 1 - A - 1 U (I + V T A - 1 U) - 1 V T A - 1

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

(A + \sum i = 1 k u i v T i) - 1 = A - 1 - A - 1 U (C) - 1 V T A - 1

C i j = δ i j + v T i A - 1 u j . i, j = 1, 2, . . . k

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

u 1 = v 1 = y k ( s T k y k ) 1 / 2, u 2 = - v 2 = B k s k ( s T k B k s k ) 1 / 2

我们可以令

Hk=B−1k

C 11 = 1 + v T 1 A - 1 u 1 = 1 + y T k H k y k s T k y k

C 22 = 1 + v T 2 A - 1 u 2 = 1 - s T k B k H k B k s k s T k B k s k = 1 - 1 = 0

C 12 = v T 1 u 2 = y T k B k s k ( s T k y k ) 1 / 2 ) ( s T k B k s k ) 1 / 2 = ( s T k B k s k ) 1 / 2 ( s T k y k ) 1 / 2

C 21 = v T 2 u 1 = - C 11

回想一下2x2矩阵的逆

C = {β - α α 0}

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

β = C 11 = 1 + y T k H k y k s T k y k, α = C 12 = ( s T k B k s k ) 1 / 2 ( s T k y k ) 1 / 2

然后就是代入了,可以令

U ˜ = H k U, V ˜ = H k V

这样，对于每一维

u ˜ i = H k u i, v ˜ i = H k v i, i = 1, 2

H k + 1 \leftarrow H k - H k U C - 1 V T H k = H k - U ˜ C - 1 V ˜ T = H k + 1 α (- u ˜ 1 v ˜ T 1) - β α 2 u ˜ 2 v T 2 = H k - H k y k s T k + s k y T k H k s T k y k + s k s T k s T k y k (1 + y T k H k y k s T k y k)

整理就得到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矩阵.

H k + 1 = (I - s k y T k y T k s k) H k (I - y k s T k y T k s k) + s k s T k y T k s k

define:

ρk=1yTksk,

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

H k + 1 = V T k H k V k + ρ k s k s T k

Then,recursively

H 1 = V T 0 H 0 V 0 + ρ 0 s 0 s T 0

H 2 = = = V T 1 H 1 V 1 + ρ 1 s 1 s T 1 V T 1 (V T 0 H 0 V 0 + ρ 0 s 0 s T 0) V 1 + ρ 1 s 1 s T 1 V T 1 V T 0 H 0 V 0 V 1 + V T 1 ρ 0 s 0 s T 0) V 1 + ρ 1 s 1 s T 1

所以就有了这个公式：

H k + 1 = + + + + (V T k V T k - 1 . . . V T 1 V T 0) H 0 (V 0 V 1 . . . V k - 1 V k) (V T k V T k - 1 . . . V T 1) ρ 1 s 1 s T 1 (V 1 . . . V k - 1 V k) . . . (V T k) ρ k - 1 s k - 1 s T k - 1 (V k) ρ k s k s T k

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

H k + 1 = + + . . . + + (V T k V T k - 1 . . . V T k - m) H 0 (V k - m . . . V k - 1 V k) (V T k V T k - 1 . . . V T k - m + 1) ρ k - m s k - m s T k - m (V k - m + 1 . . . V k - 1 V k) (V T k) ρ k - 1 s k - 1 s T k - 1 (V k) ρ k s k s T k

更新的方向：

H k + 1 \nabla f (x) = + + + + (V T k V T k - 1 . . . V T k - m) H 0 (V k - m . . . V k - 1 V k) \nabla f (x) (V T k V T k - 1 . . . V T k - m + 1) ρ k - m s k - m s T k - m (V k - m + 1 . . . V k - 1 V k) \nabla f (x) . . . (V T k) ρ k - 1 s k - 1 s T k - 1 (V k) \nabla f (x) ρ k s k s T k \nabla 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.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

2 0