优化算法——拟牛顿法之DFP

转载地址：http://blog.csdn.net/google19890102/article/details/45848439

一、牛顿法

    在博文“优化算法——牛顿法(Newton Method)”中介绍了牛顿法的思路，牛顿法具有二阶收敛性，相比较最速下降法，收敛的速度更快。在牛顿法中使用到了函数的二阶导数的信息，对于函数，其中表示向量。在牛顿法的求解过程中，首先是将函数在处展开，展开式为：

其中，，表示的是目标函数在的梯度，是一个向量。，表示的是目标函数在处的 Hesse 矩阵。省略掉最后面的高阶无穷小项，即为：

上式两边对求导，即为：

    在基本牛顿法中，取得最值的点处的导数值为，即上式左侧为。则：

求出其中的：

从上式中发现，在牛顿法中要求 Hesse 矩阵是可逆的。  

    当时，上式为：

此时，是否可以通过，，和模拟出 Hesse 矩阵的构造过程？此方法便称为拟牛顿法 (QuasiNewton)，上式称为拟牛顿方程。在拟牛顿法中，主要包括DFP拟牛顿法，BFGS拟牛顿法。

二、DFP拟牛顿法

1、DFP拟牛顿法简介

        DFP拟牛顿法也称为DFP校正方法，DFP校正方法是第一个拟牛顿法，是有Davidon最早提出，后经Fletcher和Powell解释和改进，在命名时以三个人名字的首字母命名。

对于拟牛顿方程：

化简可得：

令，可以得到：

在DFP校正方法中，假设：

2、DFP校正方法的推导

    令：，其中均为的向量。，。

    则对于拟牛顿方程可以简化为：

将代入上式：

将代入上式：

已知：为实数，为的向量。上式中，参数和解的可能性有很多，我们取特殊的情况，假设，。则：

代入上式：

令，，则：

则最终的 DFP 校正公式为：

3、DFP拟牛顿法的算法流程

    设对称正定，由上述的 DFP 校正公式确定，那么对称正定的充要条件是。

    在博文“优化算法——牛顿法(Newton Method)”中介绍了非精确的线搜索准则：Armijo搜索准则，搜索准则的目的是为了帮助我们确定学习率，还有其他的一些准则，如Wolfe准则以及精确线搜索等。在利用Armijo搜索准则时并不是都满足上述的充要条件，此时可以对DFP校正公式做些许改变：

       DFP 拟牛顿法的算法流程如下：

4、求解具体的优化问题

    求解无约束优化问题

其中，。

python 程序实现：

1. 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  
   
   
2. dfp.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. from function import *  
   10.   
   11. def dfp(fun, gfun, x0):  
   12.     result = []  
   13.     maxk = 500  
   14.     rho = 0.55  
   15.     sigma = 0.4  
   16.     m = shape(x0)[0]  
   17.     Hk = eye(m)  
   18.     k = 0  
   19.     while (k < maxk):  
   20.         gk = mat(gfun(x0))#计算梯度  
   21.         dk = -mat(Hk)*gk  
   22.         m = 0  
   23.         mk = 0  
   24.         while (m < 20):  
   25.             newf = fun(x0 + rho ** m * dk)  
   26.             oldf = fun(x0)  
   27.             if (newf < oldf + sigma * (rho ** m) * (gk.T * dk)[0,0]):  
   28.                 mk = m  
   29.                 break  
   30.             m = m + 1  
   31.           
   32.         #DFP校正  
   33.         x = x0 + rho ** mk * dk  
   34.         sk = x - x0  
   35.         yk = gfun(x) - gk  
   36.         if (sk.T * yk > 0):  
   37.             Hk = Hk - (Hk * yk * yk.T * Hk) / (yk.T * Hk * yk) + (sk * sk.T) / (sk.T * yk)  
   38.           
   39.         k = k + 1  
   40.         x0 = x  
   41.         result.append(fun(x0))  
   42.       
   43.     return result  
   
   
3. testDFP.py
   [python]  view plain  copy

    

   1. #coding:UTF-8  
   2. ''''' 
   3. Created on 2015年5月19日 
   4.  
   5. @author: zhaozhiyong 
   6. '''  
   7.   
   8. from bfgs import *  
   9. from dfp import dfp  
   10.   
   11. import matplotlib.pyplot as plt    
   12.   
   13. x0 = mat([[-1.2], [1]])  
   14. result = dfp(fun, gfun, x0)  
   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()