优化算法——拟牛顿法之L-BFGS算法

一、BFGS算法

    在“优化算法——拟牛顿法之BFGS算法”中，我们得到了BFGS算法的校正公式：

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

令，则得到：

二、BGFS算法存在的问题

    在BFGS算法中，每次都要存储近似Hesse矩阵，在高维数据时，存储浪费很多的存储空间，而在实际的运算过程中，我们需要的是搜索方向，因此出现了L-BFGS算法，是对BFGS算法的一种改进算法。在L-BFGS算法中，只保存最近的次迭代信息，以降低数据的存储空间。

三、L-BFGS算法思路

    令，，则BFGS算法中的可以表示为：

若在初始时，假定初始的矩阵，则我们可以得到：

若此时，只保留最近的步：

这样在L-BFGS算法中，不再保存完整的，而是存储向量序列和，需要矩阵时，使用向量序列和计算就可以得到，而向量序列和也不是所有都要保存，只要保存最新的步向量即可。

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

五、实验仿真

lbfgs.py

[python]  view plain  copy

1. #coding:UTF-8  
2.   
3. from numpy import *  
4. from function import *  
5.   
6. def lbfgs(fun, gfun, x0):  
7.     result = []#保留最终的结果  
8.     maxk = 500#最大的迭代次数  
9.     rho = 0.55  
10.     sigma = 0.4  
11.       
12.     H0 = eye(shape(x0)[0])  
13.       
14.     #s和y用于保存最近m个，这里m取6  
15.     s = []  
16.     y = []  
17.     m = 6  
18.       
19.     k = 1  
20.     gk = mat(gfun(x0))#计算梯度  
21.     dk = -H0 * gk  
22.     while (k < maxk):               
23.         n = 0  
24.         mk = 0  
25.         gk = mat(gfun(x0))#计算梯度  
26.         while (n < 20):  
27.             newf = fun(x0 + rho ** n * dk)  
28.             oldf = fun(x0)  
29.             if (newf < oldf + sigma * (rho ** n) * (gk.T * dk)[0, 0]):  
30.                 mk = n  
31.                 break  
32.             n = n + 1  
33.           
34.         #LBFGS校正  
35.         x = x0 + rho ** mk * dk  
36.         #print x  
37.           
38.         #保留m个  
39.         if k > m:  
40.             s.pop(0)  
41.             y.pop(0)  
42.               
43.         #计算最新的  
44.         sk = x - x0  
45.         yk = gfun(x) - gk  
46.           
47.         s.append(sk)  
48.         y.append(yk)  
49.           
50.         #two-loop的过程  
51.         t = len(s)  
52.         qk = gfun(x)  
53.         a = []  
54.         for i in xrange(t):  
55.             alpha = (s[t - i - 1].T * qk) / (y[t - i - 1].T * s[t - i - 1])  
56.             qk = qk - alpha[0, 0] * y[t - i - 1]  
57.             a.append(alpha[0, 0])  
58.         r = H0 * qk  
59.               
60.         for i in xrange(t):  
61.             beta = (y[i].T * r) / (y[i].T * s[i])  
62.             r = r + s[i] * (a[t - i - 1] - beta[0, 0])  
63.   
64.               
65.         if (yk.T * sk > 0):  
66.             dk = -r              
67.           
68.         k = k + 1  
69.         x0 = x  
70.         result.append(fun(x0))  
71.       
72.     return result  

function.py

[python]  view plain  copy

1. #coding:UTF-8  
2. ''''' 
3. Created on 2015年5月19日 
4.  
5. @author: zhaozhiyong 
6. '''  
7.   
8. from numpy import *  
9.   
10. #fun  
11. def fun(x):  
12.     return 100 * (x[0,0] ** 2 - x[1,0]) ** 2 + (x[0,0] - 1) ** 2  
13.   
14. #gfun  
15. def gfun(x):  
16.     result = zeros((2, 1))  
17.     result[0, 0] = 400 * x[0,0] * (x[0,0] ** 2 - x[1,0]) + 2 * (x[0,0] - 1)  
18.     result[1, 0] = -200 * (x[0,0] ** 2 - x[1,0])  
19.     return result  

testLBFGS.py

[python]  view plain  copy

1. #coding:UTF-8  
2. ''''' 
3. Created on 2015年6月6日 
4.  
5. @author: zhaozhiyong 
6. '''  
7.   
8. from lbfgs import *  
9.   
10. import matplotlib.pyplot as plt    
11.   
12. x0 = mat([[-1.2], [1]])  
13. result = lbfgs(fun, gfun, x0)  
14. print result  
15.   
16. n = len(result)  
17. ax = plt.figure().add_subplot(111)  
18. x = arange(0, n, 1)  
19. y = result  
20. ax.plot(x,y)  
21.   
22. plt.show()  

实验结果