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

一、牛顿法

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

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

当时，上式为：

令，可以得到：

令：，其中的向量。。

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

将代入上式：

将代入上式：

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

令则：

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

其中，。

[python] view plain copy 在CODE上查看代码片派生到我的代码片
#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  

dfp.py
[python] view plain copy 在CODE上查看代码片派生到我的代码片
#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  

testDFP.py
[python] view plain copy 在CODE上查看代码片派生到我的代码片
#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()