一、牛顿法

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

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

上式两边对求导，即为：

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

求出其中的：

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

当时，上式为：

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

二、DFP拟牛顿法

1、DFP拟牛顿法简介

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

对于拟牛顿方程：

化简可得：

令，可以得到：

在DFP校正方法中，假设：

2、DFP校正方法的推导

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

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

将代入上式：

将代入上式：

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

代入上式：

令，，则：

则最终的DFP校正公式为：

3、求解具体的优化问题

求解无约束优化问题

其中，。

python程序实现：

1. 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
   

2. dfp.py

   #coding:UTF-8
   '''
   Created on 2015年5月19日
   
   @author: zhaozhiyong
   '''
   
   from numpy import *
   from function import *
   
   def dfp(fun, gfun, x0):
       result = []
       maxk = 500
       rho = 0.55
       sigma = 0.4
       m = shape(x0)[0]
       Hk = eye(m)
       k = 0
       while (k < maxk):
           gk = mat(gfun(x0))#计算梯度
           dk = -mat(Hk)*gk
           m = 0
           mk = 0
           while (m < 20):
               newf = fun(x0 + rho ** m * dk)
               oldf = fun(x0)
               if (newf < oldf + sigma * (rho ** m) * (gk.T * dk)[0,0]):
                   mk = m
                   break
               m = m + 1
   
           #DFP校正
           x = x0 + rho ** mk * dk
           sk = x - x0
           yk = gfun(x) - gk
           if (sk.T * yk > 0):
               Hk = Hk - (Hk * yk * yk.T * Hk) / (yk.T * Hk * yk) + (sk * sk.T) / (sk.T * yk)
   
           k = k + 1
           x0 = x
           result.append(fun(x0))
   
       return result
   

3. testDFP.py

   #coding:UTF-8
   '''
   Created on 2015年5月19日
   
   @author: zhaozhiyong
   '''
   
   from bfgs import *
   from dfp import dfp
   
   import matplotlib.pyplot as plt  
   
   x0 = mat([[-1.2], [1]])
   result = dfp(fun, gfun, x0)
   
   n = len(result)
   ax = plt.figure().add_subplot(111)
   x = arange(0, n, 1)
   y = result
   ax.plot(x,y)
   
   plt.show()
   

4、实验结果