2011-04-17 20:15:00 dgglx 阅读数 16487

注明:程序中调用的函数jintuifa.m golddiv.m我在之前的笔记中已贴出

DFP算法和BFGS算法不同在于H矩阵的修正公式不同

DFP算法

%拟牛顿法中DFP算法求解f = x1*x1+2*x2*x2-2*x1*x2-4*x1的最小值,起始点为x0=[1 1]  H0为二阶单位阵
%算法根据最优化方法(天津大学出版社)116页算法3.5.1编写
%v1.0 author: liuxi BIT

%format long
syms  x1 x2 alpha;
f = x1*x1+2*x2*x2-2*x1*x2-4*x1;%要最小化的函数
df=jacobian(f,[x1 x2]);%函数f的偏导
epsilon=1e-6;%0.000001
k=1;
x0=[1 1];%起始点
xk=x0;
gk=subs(df,[x1 x2],xk);%起始点的梯度
%gk=double(gk);
H0=[1 0;0 1];%初始矩阵为二阶单位阵
while(norm(gk)>epsilon)%迭代终止条件||gk||<=epsilon
     if k==1
        pk=-H0*gk';%负梯度方向
        Hk0=H0;%HK0代表HK(k-1)
     else
        pk=-Hk*gk';
        Hk0=Hk;%HK0代表HK(k-1)
    end
    f_alpha=subs(f,[x1 x2],xk+alpha*pk');%关于alpha的函数
    [left right] = jintuifa(f_alpha,alpha);%进退法求下单峰区间
    [best_alpha best_f_alpha]=golddiv(f_alpha,alpha,left,right);%黄金分割法求最优步长
    xk=xk+best_alpha*pk';
    gk0=gk;%gk0代表g(k-1)
    gk=subs(df,[x1 x2],xk);
    %gk=double(gk);
    yk=gk-gk0;
    sk=best_alpha*pk';%sk=x(k+1)-xk
    Hk=Hk0-Hk0*yk'*yk*Hk0/(yk*Hk0*yk')+sk'*sk/(yk*sk');%修正公式
    k=k+1;
end
best_x=xk;%最优点
best_fx=subs(f,[x1 x2],best_x);%最优值

BFGS算法

%拟牛顿法中BFGS算法求解f = x1*x1+2*x2*x2-2*x1*x2-4*x1的最小值,起始点为x0=[1 1]  H0为二阶单位阵
%算法根据最优化方法(天津大学出版社)122页编写
%v1.0 author: liuxi BIT

%format long
syms  x1 x2 alpha;
f = x1*x1+2*x2*x2-2*x1*x2-4*x1;%要最小化的函数
df=jacobian(f,[x1 x2]);%函数f的偏导
epsilon=1e-6;%0.000001
k=1;
x0=[1 1];%起始点
xk=x0;
gk=subs(df,[x1 x2],xk);%起始点的梯度
%gk=double(gk);
H0=[1 0;0 1];%初始矩阵为二阶单位阵
while(norm(gk)>epsilon)%迭代终止条件||gk||<=epsilon
     if k==1
        pk=-H0*gk';%负梯度方向
        Hk0=H0;%HK0代表HK(k-1)
     else
        pk=-Hk*gk';
        Hk0=Hk;%HK0代表HK(k-1)
    end
    f_alpha=subs(f,[x1 x2],xk+alpha*pk');%关于alpha的函数
    [left right] = jintuifa(f_alpha,alpha);%进退法求下单峰区间
    [best_alpha best_f_alpha]=golddiv(f_alpha,alpha,left,right);%黄金分割法求最优步长
    xk=xk+best_alpha*pk';
    gk0=gk;%gk0代表g(k-1)
    gk=subs(df,[x1 x2],xk);
    %gk=double(gk);
    yk=gk-gk0;
    sk=best_alpha*pk';%sk=x(k+1)-xk
    %====begin=============与DFP算法不同的地方==============
    wk=(yk*Hk*yk')^0.5*(sk'/(yk*sk')-Hk*yk'/(yk*Hk*yk'));
    Hk=Hk0-Hk0*yk'*yk*Hk0/(yk*Hk0*yk')+sk'*sk/(yk*sk')+wk*wk';%修正公式
    %====end===============与DFP算法不同的地方==============
    k=k+1;
end
best_x=xk%最优点
best_fx=subs(f,[x1 x2],best_x)%最优值

2019-11-24 15:49:13 weixin_41514525 阅读数 10

DFP算法是一种拟牛顿算法,它的好处是在于 继承了牛顿法的快速,又克服了牛顿法中每次求取G的逆矩阵的困难,且保证了Hk的正定,也就克服了牛顿法的G可能正定的缺点。

 

from 实用优化算法.helper import*   # 单独博客介绍
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


x1,x2 = symbols('x1,x2')

f = 100 * (x2 - x1 ** 2) ** 2 + (1 - x1) ** 2


def do_work(x0,h):
    k = 1
    step_x1 = [x0[0][0]]
    step_x2 = [x0[1][0]]
    while get_len_grad(get_grad(x0,f)) > 0.00001:
        d = -1 * np.dot(h, get_grad(x0,f))
        step = golden_search(0, 2, x0, d, f)
        x1 = x0 + step * d
        s = x1 - x0
        y = get_grad(x1,f) - get_grad(x0,f)
        t1 = np.dot(h,y)
        t2 = np.dot(t1,y.T)
        t3 = np.dot(t2,h)  # 分子
        t4 = np.dot(y.T,h)
        t5 = np.dot(t4,y)  # 分母
        t6 = np.dot(s,s.T)  # 分子
        t7 = np.dot(y.T, s)  # 分母
        h = h - t3*(1/t5[0][0]) + (1/t7[0][0]) * t6
        x0 = x1
        step_x1.append(x0[0][0])
        step_x2.append(x0[1][0])
        print(k,' ',x0)
        k += 1
    print('\n结束',x0)
    return step_x1,step_x2


if __name__ == '__main__':
    x0 = [[0],[0]]
    x0 = np.array(x0)
    h = [[1,0],[0,1]]
    h = np.array(h)
    step_x1,step_x2 = do_work(x0,h)

    fig = plt.figure()
    ax = Axes3D(fig)
    x = np.arange(-0.1, 1.5, 0.1)
    y = np.arange(-1, 2.1, 0.1)
    X1, X2 = np.meshgrid(x, y)
    plt.figure(figsize=(10, 6))
    Z = 100 * (X2 - X1 ** 2) ** 2 + (1 - X1) ** 2
    step_x1 = np.array(step_x1)
    step_x2 = np.array(step_x2)
    step_x3 = 100 * (step_x2 - step_x1 ** 2) ** 2 + (1 - step_x1) ** 2
    bar = plt.contourf(X1, X2, Z, 45, cmap=plt.cm.Blues)
    plt.plot(step_x1, step_x2, marker='*')
    ax.plot_surface(X1, X2, Z, rstride=1, cstride=1)
    ax.plot(step_x1, step_x2, step_x3, c='r',marker = '*')
    plt.colorbar(bar)
    plt.show()

 

2015-05-19 22:11:04 google19890102 阅读数 13151

一、牛顿法

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



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



上式两边对求导,即为:



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



求出其中的

从上式中发现,在牛顿法中要求Hesse矩阵是可逆的。  
    当时,上式为:



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

二、DFP拟牛顿法

1、DFP拟牛顿法简介

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



化简可得:



,可以得到:



DFP校正方法中,假设:



2、DFP校正方法的推导

    令:,其中均为的向量。
    则对于拟牛顿方程可以简化为:



代入上式:



代入上式:





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



代入上式:





,则:






则最终的DFP校正公式为:



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

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



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

4、求解具体的优化问题

    求解无约束优化问题



其中,

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()
    

5、实验结果







2017-04-25 17:00:43 philosophyatmath 阅读数 876

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/

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