拟牛顿法之DFP算法

DFP算法(Davidon-Fletcher-Powell algorithm)一种秩2拟牛顿法.由戴维登(Davidon, W. D.)于1959年导出,并由弗莱彻(Fletcher,R.)和鲍威尔(Powe11,M. J. D.)于1963年进行了改善.

对函数f(x)x=xk+1处进行泰勒展开到二阶:

f(x)=f(xk+1)+f(xk+1)(xxk+1)+12f′′(xk+1)(xxk+1)2+Rn(x)f(xk+1)+f(xk+1)(xxk+1)+12f′′(xk+1)(xxk+1)2

对上式求导并令其为0,由于f(x)中的x是一个向量,f(x)x求导意味着对x向量中的每个值求偏导。即,f(x)x的一阶导数为一个向量,对x的二阶导数为一个nn的矩阵
f(x)=(f(x)x1f(x)x2,,f(x)xn)f′′(x)=[2f(x)xixj]nn

gTk+1=f(xk+1)=f(xk+1),表示的是目标函数在xk+1的梯度,是一个向量。Gk+1=2f(xk+1)=f′′(xk+1)(xxk+1)2表示的是目标函数在xk+1处的Hesse矩阵。
求导后得:
f(x)=f(xk+1)+Gk+1(xxk+1)

在基本牛顿法中,取得最值的点处的导数值为0,即上式左侧为0。则:
f(xk+1)+Gk+1(xxk+1)=0

求出其中的x
x=xk+1G1k+1f(xk+1)

从上式中发现,在牛顿法中要求Hesse矩阵是可逆的。
x=xk时,上式为(拟牛顿方程):

f(k)=f(xk+1)+Gk+1(xkxk+1)

此时,是否可以通过xk,xk+1,f(xk),f(xk+1)模拟出Hesse矩阵的构造过程?此方法便称为拟牛顿法(QuasiNewton),上式称为拟牛顿方程。

DFP拟牛顿法

DFP拟牛顿法简介

对于拟牛顿方程:

f(k)=f(xk+1)+Gk+1(xkxk+1)

化简可得:
G1k+1[f(xk+1)f(xk)]=xk+1xk

Hk+1G1k+1,可以得到:
Hk+1[f(xk+1)f(xk)]=xk+1xk

在DFP校正方法中,假设:
Hk+1=Hk+Ek

DFP校正方法的推导

令:Ek=αukuTk+βvkvTk,其中uk,vk均为n×1的向量。其中:

yksk=f(xk+1)f(xk),=xk+1xk

则对于Hk+1[f(xk+1)f(xk)]=xk+1xk可以简化为:

Hk+1yk=sk

Hk+1=Hk+Ek代入上式:
(Hk+Ek)yk=sk

Ek=αukuTk+βvkvTk代入上式:
(Hk+αukuTk+βvkvTk)ykα(uTkyk)uk+β(vTkyk)vk=sk=skHkyk

已知:uTkyk,vTkyk为实数,skHkykn×1的向量。上式中,参数α,β解的可能性有很多,我们取特殊的情况,假设uk=rHkyk,vk=θsk。则:
Ek=αr2HkykyTkHk+βθ2sksTk

代入上式:
α[(rHkyk)Tkyk)](rHkyk)+β[(θsk)Tkyk)(θsk)=skHkyk[αr2(yTkHkyk)+1](Hkyk)+[βθ2(sTkyk)1](xk)=0

αr2(yTkHkyk)+1=0,βθ2(sTkyk)1=0,则:
αr2βθ2=1yTkHkyk=1sTkyk

则最终的DFP校正公式为:
Hk+1=HkHkykyTkHkyTkHkyk+sksTksTkyk

DFP拟牛顿法的算法流程

Hk对称正定,Hk+1由上述的DFP校正公式确定,那么Hk+1对称正定的充要条件是sTkyk>0
Armijo搜索准则,搜索准则的目的是为了帮助我们确定学习率,还有其他的一些准则,如Wolfe准则以及精确线搜索等。在利用Armijo搜索准则时并不是都满足上述的充要条件,此时可以对DFP校正公式做些许改变:

Hk+1=Hk,HkHkykyTkHkyTkHkyk+sksTksTksk,ifsTkyk0ifsTkyk>0

DFP拟牛顿法的算法流程如下:
1. 给定参数δ(0,1),σ(0,0.5),初始化点x0Rn,终止误差0ϵ1,初始化对称正定阵H0,通常取为G(xo)1或单位阵In;令k=0
2. 计算gk=f(xk),若gkϵ,终止,输出xk作为近似极小点。
3. 计算搜索方向 :dk=Hkgk.
4. 设mk是满足下列不等式的最小非负整数m:
f(xk+δmdk)f(xk)+σδmgTkdk

αk=δmk,xk+1=xk+αkdk.
5. 由校正公式确定Hk+1
6. 令k=k+1,转向步骤“2”

代码:

dfp.py

#coding:UTF-8
'''
Created on 2017年4月25日

@author: zhangdapeng
'''
from numpy import *
from function import *

def dfp(fun, gfun, x0):
    result = []
    maxk = 500
    delta = 0.55
    sigma = 0.4
    m = shape(x0)[0]
    Hk = eye(m)
    k = 0
    epsilon=1e-10

    while (k < maxk):
        gk = mat(gfun(x0))#计算梯度
        if linalg.norm(gk,1)<epsilon:
            break
        dk = -mat(Hk)*gk
        m = 0
        mk = 0
        while (m < 20):
            newf = fun(x0 + delta ** m * dk)
            oldf = fun(x0)
            if (newf < oldf + sigma * (delta ** m) * (gk.T * dk)[0,0]):
                mk = m
                break
            m = m + 1

        #DFP校正
        x = x0 + delta ** 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

function.py

#coding:UTF-8
'''
Created on 2017年4月24日

@author: zhangdapeng
'''

from numpy import *

#fun  原始函数
def fun(x):  
    return 100 * (x[0,0] ** 2 - x[1,0]) ** 2 + (x[0,0] - 1) ** 2  

#对x1,x2求导后的函数  
def gfun(x):  
    result = zeros((2, 1))  
#     对x1求导
    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])  #对x2求导
    return result

testDFP.py

#coding:UTF-8
'''
Created on 2017年4月25日

@author: zhangdapeng
'''
from dfp import *

import matplotlib.pyplot as plt  

x0 = mat([[-1.2], [1]])
result = dfp(fun, gfun, x0)
print(result[-1])
n = len(result)
ax = plt.figure().add_subplot(111)
x = arange(0, n, 1)
y = result
ax.plot(x,y)

plt.show()

输出结果:

1.54074395551e-28

输出图片:

这里写图片描述

http://blog.csdn.net/google19890102/article/details/45848439
http://blog.csdn.net/u012102306/article/details/50534086
http://www.codelast.com/%E5%8E%9F%E5%88%9B%E6%8B%9F%E7%89%9B%E9%A1%BF%E6%B3%95quasi-newton%EF%BC%8Cdfp%E7%AE%97%E6%B3%95davidon-fletcher-powell%EF%BC%8C%E5%8F%8Abfgs%E7%AE%97%E6%B3%95broyden-fletcher-goldfarb-shanno/

发布了203 篇原创文章 · 获赞 68 · 访问量 73万+
展开阅读全文

没有更多推荐了,返回首页

©️2019 CSDN 皮肤主题: 大白 设计师: CSDN官方博客

分享到微信朋友圈

×

扫一扫,手机浏览