一、BFGS算法

BFGS算法的校正公式：

利用Sherman-Morrison公式可对上式进行变换，得到

令

，则得到：

二、BGFS算法存在的问题

在BFGS算法中。每次都要存储近似Hesse矩阵

，在高维数据时，存储

浪费非常多的存储空间，而在实际的运算过程中。我们须要的是搜索方向。因此出现了L-BFGS算法。是对BFGS算法的一种改进算法。

在L-BFGS算法中。仅仅保存近期的

次迭代信息。以减少数据的存储空间。

三、L-BFGS算法思路

令

。

，则BFGS算法中的

能够表示为：

若在初始时，假定初始的矩阵

，则我们能够得到：

若此时。仅仅保留近期的

步：

这样在L-BFGS算法中。不再保存完整的

。而是存储向量序列

和

。须要矩阵

时，使用向量序列

和

计算就能够得到。而向量序列

和

也不是全部都要保存，仅仅要保存最新的

步向量就可以。

四、L-BFGS算法中的方向的计算方法

五、实验仿真

lbfgs.py

#coding:UTF-8

from numpy import *

from function import *

def lbfgs(fun, gfun, x0):

result = []#保留终于的结果

maxk = 500#最大的迭代次数

rho = 0.55

sigma = 0.4

H0 = eye(shape(x0)[0])

#s和y用于保存近期m个，这里m取6

s = []

y = []

m = 6

k = 1

gk = mat(gfun(x0))#计算梯度

dk = -H0 * gk

while (k < maxk):

n = 0

mk = 0

gk = mat(gfun(x0))#计算梯度

while (n < 20):

newf = fun(x0 + rho ** n * dk)

oldf = fun(x0)

if (newf < oldf + sigma * (rho ** n) * (gk.T * dk)[0, 0]):

mk = n

break

n = n + 1

#LBFGS校正

x = x0 + rho ** mk * dk

#print x

#保留m个

if k > m:

s.pop(0)

y.pop(0)

#计算最新的

sk = x - x0

yk = gfun(x) - gk

s.append(sk)

y.append(yk)

#two-loop的过程

t = len(s)

qk = gfun(x)

a = []

for i in xrange(t):

alpha = (s[t - i - 1].T * qk) / (y[t - i - 1].T * s[t - i - 1])

qk = qk - alpha[0, 0] * y[t - i - 1]

a.append(alpha[0, 0])

r = H0 * qk

for i in xrange(t):

beta = (y[i].T * r) / (y[i].T * s[i])

r = r + s[i] * (a[t - i - 1] - beta[0, 0])

if (yk.T * sk > 0):

dk = -r

k = k + 1

x0 = x

result.append(fun(x0))

return result

function.py

#coding:UTF-8

'''

Created on 2015年5月19日

@author: zhaozhiyong

'''

from numpy import *

#fun

def fun(x):

return 100 * (x[0,0] ** 2 - x[1,0]) ** 2 + (x[0,0] - 1) ** 2

#gfun

def gfun(x):

result = zeros((2, 1))

result[0, 0] = 400 * x[0,0] * (x[0,0] ** 2 - x[1,0]) + 2 * (x[0,0] - 1)

result[1, 0] = -200 * (x[0,0] ** 2 - x[1,0])

return result

testLBFGS.py

#coding:UTF-8

'''

Created on 2015年6月6日

@author: zhaozhiyong

'''

from lbfgs import *

import matplotlib.pyplot as plt

x0 = mat([[-1.2], [1]])

result = lbfgs(fun, gfun, x0)

print result

n = len(result)

ax = plt.figure().add_subplot(111)

x = arange(0, n, 1)

y = result

ax.plot(x,y)

plt.show()

实验结果

參考文献