精华内容
下载资源
问答
  • BFGS和DFP公式的拟牛顿pptimization算法的数学优化代码 您将需要使用Octave / Matlab。 简短的教程 dim = [2,1] vnn(@convex,dim,@DFP,1000)vnn(@convex,dim,@BFGS,1000)vnn(@rossenbrock,dim,@DFP,5...
  • DFP算法

    2020-12-31 23:34:10
    一、算法公式 二、步骤 1. 给定,误差 2. 令,算出,令 3 令

    一、算法公式

    H_{k+1}=H_{k}+\Delta H_{k}=H_{k}+\frac{\Delta x^{k}\left ( \Delta x^{k} \right )^{T}}{\left ( \Delta x^{k} \right )^{T}\Delta g^{k}}-\frac{H_{k} \Delta g^{k}\left ( \Delta g^{k} \right )^{T}{H_{k}}}{\left ( \Delta g^{k} \right )^{T}H_{k} \Delta g^{k}}

    二、步骤

    1. 给定x^{0}\epsilon R^n,,误差\varepsilon >0

    2. 令H_{0}:=E_{n},算出g^{0}:=\bigtriangledown f\left ( x^0 \right ),令k:=0

    3  令p^k=-H_{k} g^{k}}

     

    展开全文
  • 3.1 DFP公式推导 3.2 要求解的问题 3.3 python实现 1.拟牛顿法 1.1拟牛顿法的导出与优点 在上一文中(牛顿法公式推导与python实现),谈到说牛顿法需要计算一个Hessian矩阵的逆,才能够迭代,但在实际工程中...

    目录

    1. 拟牛顿法 
      1.1拟牛顿法的导出与优点 
      1.2 算法步骤与特点
    2. 对称秩一校正公式
    3. DFP算法 
      3.1 DFP公式推导 
      3.2 要求解的问题 
      3.3 python实现

    1.拟牛顿法

    1.1拟牛顿法的导出与优点

    在上一文中(牛顿法公式推导与python实现),谈到说牛顿法需要计算一个Hessian矩阵的逆,才能够迭代,但在实际工程中,计算如此大型的矩阵需要很大的计算资源,因此,有人提出能否不计算Hessian矩阵,在迭代过程中,仅仅利用相邻两个迭代点以及梯度信息,产生一个对称正定矩阵,使之逐步逼近目标函数Hessian矩阵的逆阵。

    其实这就是你牛顿法的基本思想,这样做,既能保存Hessian矩阵的大部分信息(曲率),也能极大的减小计算量。

    考虑无约束极小化问题。假设目标函数f:Rn→Rf:Rn→R是二次连续可微的,那么∇f∇f在xk+1xk+1处的泰勒展开为:

    ∇f(x)=∇f(xk+1)+∇2f(xk+1)(x−xk+1)+o||x−xk+1||∇f(x)=∇f(xk+1)+∇2f(xk+1)(x−xk+1)+o||x−xk+1||

    ,取x:=xkx:=xk.当xk与xk+1xk与xk+1充分接近时,有:

    ∇2f(xk+1)(xk+1−xk)≈∇f(xk+1)−∇f(xk)∇2f(xk+1)(xk+1−xk)≈∇f(xk+1)−∇f(xk)

    ∇2f(xk+1)∇2f(xk+1)就是f(x)f(x)xk+1xk+1处的Hessian矩阵,那么我们可以用它的近似矩阵Bk+1Bk+1来代替它,得到如下等式:

    Bk+1(xk+1−xk)=∇f(xk+1)−∇f(xk)(1)(1)Bk+1(xk+1−xk)=∇f(xk+1)−∇f(xk)

    如该矩阵存在逆矩阵有:

    Hk+1(∇f(xk+1)−∇f(xk))=xk+1−xk(2)(2)Hk+1(∇f(xk+1)−∇f(xk))=xk+1−xk

    以上两个方程成为拟牛顿方程(条件)。其中Hk+1=∇2f(xk+1)−1Hk+1=∇2f(xk+1)−1,为Hessian的逆阵。

     

    1.2 算法步骤与特点

    拟牛顿法的算法步骤如下:

    1. 给出x0∈Rn,H0∈Rnxn,0≤ϵ<1,k:=0x0∈Rn,H0∈Rnxn,0≤ϵ<1,k:=0;

    2. 若|∇f(xk)|≤ϵ|∇f(xk)|≤ϵ,迭代停止;否则求方向:dk=−Hk∇f(xk)dk=−Hk∇f(xk)

    3. 沿着方向做线性搜索αk>0αk>0,令xk+1=xk+αkdkxk+1=xk+αkdk

    4. 校正Hk产生Hk+1,使得牛顿条件(2)Hk产生Hk+1,使得牛顿条件(2)依然成立

    5. k:=k+1,转至第二步

    总结一下拟牛顿法的特点:

    • 这种情况下,Hk+1≈∇2f(xk+1)−1Hk+1≈∇2f(xk+1)−1,使得算法产生的方向近似于牛顿方向,从而确保算法具有比较好的收敛性。
    • 对任意的k,近似矩阵Bk+1Bk+1都是正定的,使得算法选取的方向(dk=−Hk∇f(xk)dk=−Hk∇f(xk))都是下降方向。
    • 仅需一阶导数,就能完整整个迭代过程
    • 需要校正HkHk产生Hk+1Hk+1

    2.对称秩一校正公式

    前面我们说过要用Hk+1Hk+1来近似Hessian的逆阵,但不可能说一次取值,就能得到最优的Hk+1Hk+1,所以我们接下来讨论一下,如何通过迭代,不断的校正这个近似矩阵,使得:

    Hk+1=Hk+Ek(3)(3)Hk+1=Hk+Ek

    在秩一校正情形下,有:

    Hk+1=Hk+uvT(4)(4)Hk+1=Hk+uvT

    其中rank(uvT)=1(秩为1)。rank(uvT)=1(秩为1)。

     

    它的想法是希望通过以上这个迭代公式,将u,vTu,vT换成我们可以求得的xk,∇f(x)xk,∇f(x)等,达到的迭代的效果。

    令sk=xk+1−xk,yk=∇f(xk+1)−∇f(xk)sk=xk+1−xk,yk=∇f(xk+1)−∇f(xk)将Hk+1Hk+1代入(2)有:

    Hk+1yk=(Hk+uvT)yk=sk(5)(5)Hk+1yk=(Hk+uvT)yk=sk

    (vTyk)u=sk−Hkyk(6)(6)(vTyk)u=sk−Hkyk

    故u必定在sk−Hkyksk−Hkyk方向上,且sk−Hkyk≠0sk−Hkyk≠0(如果等于0,则已经满足拟牛顿条件了),则u=sk−HkykvTyku=sk−HkykvTyk,代入(4),我们有:

    Hk+1=Hk+(sk−Hkyk)vTvTyk(7)(7)Hk+1=Hk+(sk−Hkyk)vTvTyk

    由于要求Hess矩阵对称,故其逆也必定对称,故v=sk−Hkykv=sk−Hkyk,有

    Hk+1=Hk+(sk−Hkyk)(sk−Hkyk)TvTyk(8)(8)Hk+1=Hk+(sk−Hkyk)(sk−Hkyk)TvTyk

    该公式被称为对称秩一校正公式,可以用它来校正我们要校正的HkHk.

     

    3.DFP算法

    3.1 DFP公式推导

    由前面的对称秩一校正公式的导出,我们发现把末尾的未知参数用已知参数代替后,就能完成校正的功能,但对称秩一校正的效果并不是太好,我们可以再加一个校正,让他们协调一下,就有了DFP算法。

    DFP算法是设出一个对称秩二校正:

    Hk+1=Hk+auuT+bvvT(9)(9)Hk+1=Hk+auuT+bvvT

    在满足你牛顿条件的情况下,将式中所有的未知参数a,u,b,va,u,b,v都用已知条件代替,得到一个迭代公式,校正HKHK

     

    用同样的思想,我们有:

    Hkyk+auuTyk+bvvTyk=sk(10)(10)Hkyk+auuTyk+bvvTyk=sk

    这里的u,v都不是唯一确定的,但很明显,如果要让等式成立,有:

    u=sk,v=Hkyk(11)(11)u=sk,v=Hkyk

    与(10)联立,可得:

    auTyk=1,bvTyk=−1auTyk=1,bvTyk=−1

    确定出:

    a=1uTyk=1sTkyk,b=−1yTkHkyka=1uTyk=1skTyk,b=−1ykTHkyk

    最后得到DFP公式:

    Hk+1=Hk+sksTksTkyk−HkykyTkHkyTkHkykHk+1=Hk+skskTskTyk−HkykykTHkykTHkyk


    注意式中的分数结构,分子sTkyk,yTkHkykskTyk,ykTHkyk都是标量,分母sksTk,HkykyTkHkskskT,HkykykTHk则是与HkHk同型的矩阵,且都是正定矩阵。若Hk为Hk为正定矩阵,且sTkyk>0skTyk>0,则Hk+1Hk+1也正定。

     

    当采用精确线搜索时,矩阵序列HkHk的正定新条件sTkyk>0skTyk>0可以被满足。但对于Armijo搜索准则来说,不能满足这一条件,需要做如下修正: 

    Hk+1=⎧⎩⎨HkHk−HkykyTkHkyTkHkyk+sksTksTkyksTkyk≤0sTkyk>0⎫⎭⎬Hk+1={HkskTyk≤0Hk−HkykykTHkykTHkyk+skskTskTykskTyk>0}

     

    3.2 要求解的问题

    求解无约束线性优化问题

    minx∈R2f(x)=100(x21−x2)2+(x1−1)2minx∈R2f(x)=100(x12−x2)2+(x1−1)2

    该问题有精确解x∗=(1,1)T,f(x∗)=0x∗=(1,1)T,f(x∗)=0其梯度为

    g(x)=(400x1(x21−x2)+2(x1−1),−200(x21−x2))g(x)=(400x1(x12−x2)+2(x1−1),−200(x12−x2))

    其Hessian矩阵为

    G(x)=[1200x21−400x2+2−400x1−400x1200]G(x)=[1200x12−400x2+2−400x1−400x1200]

     

    3.3 python实现

    由1.2的算法步骤,可得:

    import numpy as np
    
    #函数表达式
    fun = lambda x:100*(x[0]**2 - x[1]**2)**2 +(x[0] - 1)**2
    
    #梯度向量
    gfun = lambda x:np.array([400*x[0]*(x[0]**2 - x[1]) + 2*(x[0] - 1),-200*(x[0]**2 - x[1])])
    
    #Hessian矩阵
    hess = lambda x:np.array([[1200*x[0]**2 - 400*x[1] + 2,-400*x[0]],[-400*x[0],200]])
    
    def dfp(fun,gfun,hess,x0):
        #功能:用DFP算法求解无约束问题:min fun(x)
        #输入:x0式初始点,fun,gfun,hess分别是目标函数和梯度,Hessian矩阵格式
        #输出:x,val分别是近似最优点,最优解,k是迭代次数
        maxk = 1e5
        rho = 0.05
        sigma = 0.4
        epsilon = 1e-5 #迭代停止条件
        k = 0
        n = np.shape(x0)[0]
        #将Hessian矩阵初始化为单位矩阵
        Hk = np.linalg.inv(hess(x0))
    
        while k < maxk:
            gk = gfun(x0)
            if np.linalg.norm(gk) < epsilon:
                break
            dk = -1.0*np.dot(Hk,gk)
    #         print dk
    
            m = 0;
            mk = 0
            while m < 20:#用Armijo搜索步长
                if fun(x0 + rho**m*dk) < fun(x0) + sigma*rho**m*np.dot(gk,dk):
                    mk = m
                    break
                m += 1
            #print mk
            #DFP校正
            x = x0 + rho**mk*dk
            print "第"+str(k)+"次的迭代结果为:"+str(x)
            sk = x - x0
            yk = gfun(x) - gk
    
            if np.dot(sk,yk) > 0:
                Hy = np.dot(Hk,yk)
                sy = np.dot(sk,yk) #向量的点积
                yHy = np.dot(np.dot(yk,Hk),yk) #yHy是标量
                Hk = Hk - 1.0*Hy.reshape((n,1))*Hy/yHy + 1.0*sk.reshape((n,1))*sk/sy
    
            k += 1
            x0 = x
        return x0,fun(x0),k
    
    x0 ,fun0 ,k = dfp(fun,gfun,hess,np.array([0,0]))
    print x0,fun0,k
    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21
    • 22
    • 23
    • 24
    • 25
    • 26
    • 27
    • 28
    • 29
    • 30
    • 31
    • 32
    • 33
    • 34
    • 35
    • 36
    • 37
    • 38
    • 39
    • 40
    • 41
    • 42
    • 43
    • 44
    • 45
    • 46
    • 47
    • 48
    • 49
    • 50
    • 51
    • 52
    • 53
    • 54
    • 55
    • 56
    • 57

    输出:

    第0次的迭代结果为:[ 0.05  0.  ]
    第1次的迭代结果为:[ 0.08583333  0.0015    ]
    第2次的迭代结果为:[ 0.10536555  0.00351201]
    -----
    第53次的迭代结果为:[ 1.00007963  1.00015789]
    第54次的迭代结果为:[ 1.00000251  1.00000578]
    第55次的迭代结果为:[ 1.00000079  1.00000187]
    第56次的迭代结果为:[ 1.  1.]
    [ 1.  1.] 7.69713624862e-16 57
    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9

    迭代57次后得到解(1,1)

    展开全文
  • 3.1 DFP公式推导 3.2 要求解的问题 3.3 python实现 1.拟牛顿法 1.1拟牛顿法的导出与优点 在上一文中(牛顿法公式推导与python实现),谈到说牛顿法需要计算一个Hessian矩阵的逆,才能够迭代,但在实际...

    目录

    1. 拟牛顿法
      1.1拟牛顿法的导出与优点
      1.2 算法步骤与特点
    2. 对称秩一校正公式
    3. DFP算法
      3.1 DFP公式推导
      3.2 要求解的问题
      3.3 python实现

    1.拟牛顿法

    1.1拟牛顿法的导出与优点

    在上一文中(牛顿法公式推导与python实现),谈到说牛顿法需要计算一个Hessian矩阵的逆,才能够迭代,但在实际工程中,计算如此大型的矩阵需要很大的计算资源,因此,有人提出能否不计算Hessian矩阵,在迭代过程中,仅仅利用相邻两个迭代点以及梯度信息,产生一个对称正定矩阵,使之逐步逼近目标函数Hessian矩阵的逆阵。

    其实这就是你牛顿法的基本思想,这样做,既能保存Hessian矩阵的大部分信息(曲率),也能极大的减小计算量。

    考虑无约束极小化问题。假设目标函数f:RnR是二次连续可微的,那么fxk+1处的泰勒展开为:

    f(x)=f(xk+1)+2f(xk+1)(xxk+1)+o||xxk+1||
    ,取x:=xk.当xkxk+1充分接近时,有:
    2f(xk+1)(xk+1xk)f(xk+1)f(xk)
    2f(xk+1)就是f(x)xk+1处的Hessian矩阵,那么我们可以用它的近似矩阵Bk+1来代替它,得到如下等式:
    (1)Bk+1(xk+1xk)=f(xk+1)f(xk)
    如该矩阵存在逆矩阵有:
    (2)Hk+1(f(xk+1)f(xk))=xk+1xk
    以上两个方程成为拟牛顿方程(条件)。其中Hk+1=2f(xk+1)1,为Hessian的逆阵。

    1.2 算法步骤与特点

    拟牛顿法的算法步骤如下:

    1. 给出x0Rn,H0Rnxn,0ϵ<1,k:=0;

    2. |f(xk)|ϵ,迭代停止;否则求方向:dk=Hkf(xk)

    3. 沿着方向做线性搜索αk>0,令xk+1=xk+αkdk

    4. 校正HkHk+1,使2依然成立

    5. k:=k+1,转至第二步

    总结一下拟牛顿法的特点:

    • 这种情况下,Hk+12f(xk+1)1,使得算法产生的方向近似于牛顿方向,从而确保算法具有比较好的收敛性。
    • 对任意的k,近似矩阵Bk+1都是正定的,使得算法选取的方向(dk=Hkf(xk))都是下降方向。
    • 仅需一阶导数,就能完整整个迭代过程
    • 需要校正Hk产生Hk+1

    2.对称秩一校正公式

    前面我们说过要用Hk+1来近似Hessian的逆阵,但不可能说一次取值,就能得到最优的Hk+1,所以我们接下来讨论一下,如何通过迭代,不断的校正这个近似矩阵,使得:

    (3)Hk+1=Hk+Ek
    在秩一校正情形下,有:
    (4)Hk+1=Hk+uvT
    其中rank(uvT)=1(1)

    它的想法是希望通过以上这个迭代公式,将u,vT换成我们可以求得的xk,f(x)等,达到的迭代的效果。

    sk=xk+1xk,yk=f(xk+1)f(xk)Hk+1代入(2)有:

    (5)Hk+1yk=(Hk+uvT)yk=sk
    (6)(vTyk)u=skHkyk
    故u必定在skHkyk方向上,且skHkyk0(如果等于0,则已经满足拟牛顿条件了),则u=skHkykvTyk,代入(4),我们有:
    (7)Hk+1=Hk+(skHkyk)vTvTyk
    由于要求Hess矩阵对称,故其逆也必定对称,故v=skHkyk,有
    (8)Hk+1=Hk+(skHkyk)(skHkyk)TvTyk
    该公式被称为对称秩一校正公式,可以用它来校正我们要校正的Hk.

    3.DFP算法

    3.1 DFP公式推导

    由前面的对称秩一校正公式的导出,我们发现把末尾的未知参数用已知参数代替后,就能完成校正的功能,但对称秩一校正的效果并不是太好,我们可以再加一个校正,让他们协调一下,就有了DFP算法。

    DFP算法是设出一个对称秩二校正:

    (9)Hk+1=Hk+auuT+bvvT
    在满足你牛顿条件的情况下,将式中所有的未知参数a,u,b,v都用已知条件代替,得到一个迭代公式,校正HK

    用同样的思想,我们有:

    (10)Hkyk+auuTyk+bvvTyk=sk
    这里的u,v都不是唯一确定的,但很明显,如果要让等式成立,有:
    (11)u=sk,v=Hkyk
    与(10)联立,可得:
    auTyk=1,bvTyk=1
    确定出:
    a=1uTyk=1skTyk,b=1ykTHkyk
    最后得到DFP公式:
    Hk+1=Hk+skskTskTykHkykykTHkykTHkyk

    注意式中的分数结构,分子skTyk,ykTHkyk都是标量,分母skskT,HkykykTHk则是与Hk同型的矩阵,且都是正定矩阵。若Hk正定矩阵,且skTyk>0,则Hk+1也正定。

    当采用精确线搜索时,矩阵序列Hk的正定新条件skTyk>0可以被满足。但对于Armijo搜索准则来说,不能满足这一条件,需要做如下修正:

    (1)Hk+1={HkskTyk0HkHkykykTHkykTHkyk+skskTskTykskTyk>0}

    3.2 要求解的问题

    求解无约束线性优化问题

    minxR2f(x)=100(x12x2)2+(x11)2
    该问题有精确解x=(1,1)T,f(x)=0其梯度为
    g(x)=(400x1(x12x2)+2(x11),200(x12x2))
    其Hessian矩阵为
    G(x)=[1200x12400x2+2400x1400x1200]

    3.3 python实现

    由1.2的算法步骤,可得:

    import numpy as np
    
    #函数表达式
    fun = lambda x:100*(x[0]**2 - x[1]**2)**2 +(x[0] - 1)**2
    
    #梯度向量
    gfun = lambda x:np.array([400*x[0]*(x[0]**2 - x[1]) + 2*(x[0] - 1),-200*(x[0]**2 - x[1])])
    
    #Hessian矩阵
    hess = lambda x:np.array([[1200*x[0]**2 - 400*x[1] + 2,-400*x[0]],[-400*x[0],200]])
    
    def dfp(fun,gfun,hess,x0):
        #功能:用DFP算法求解无约束问题:min fun(x)
        #输入:x0式初始点,fun,gfun,hess分别是目标函数和梯度,Hessian矩阵格式
        #输出:x,val分别是近似最优点,最优解,k是迭代次数
        maxk = 1e5
        rho = 0.05
        sigma = 0.4
        epsilon = 1e-5 #迭代停止条件
        k = 0
        n = np.shape(x0)[0]
        #将Hessian矩阵初始化为单位矩阵
        Hk = np.linalg.inv(hess(x0))
    
        while k < maxk:
            gk = gfun(x0)
            if np.linalg.norm(gk) < epsilon:
                break
            dk = -1.0*np.dot(Hk,gk)
    #         print dk
    
            m = 0;
            mk = 0
            while m < 20:#用Armijo搜索步长
                if fun(x0 + rho**m*dk) < fun(x0) + sigma*rho**m*np.dot(gk,dk):
                    mk = m
                    break
                m += 1
            #print mk
            #DFP校正
            x = x0 + rho**mk*dk
            print "第"+str(k)+"次的迭代结果为:"+str(x)
            sk = x - x0
            yk = gfun(x) - gk
    
            if np.dot(sk,yk) > 0:
                Hy = np.dot(Hk,yk)
                sy = np.dot(sk,yk) #向量的点积
                yHy = np.dot(np.dot(yk,Hk),yk) #yHy是标量
                Hk = Hk - 1.0*Hy.reshape((n,1))*Hy/yHy + 1.0*sk.reshape((n,1))*sk/sy
    
            k += 1
            x0 = x
        return x0,fun(x0),k
    
    x0 ,fun0 ,k = dfp(fun,gfun,hess,np.array([0,0]))
    print x0,fun0,k

    输出:

    第0次的迭代结果为:[ 0.05  0.  ]
    第1次的迭代结果为:[ 0.08583333  0.0015    ]
    第2次的迭代结果为:[ 0.10536555  0.00351201]
    -----
    第53次的迭代结果为:[ 1.00007963  1.00015789]
    第54次的迭代结果为:[ 1.00000251  1.00000578]
    第55次的迭代结果为:[ 1.00000079  1.00000187]
    第56次的迭代结果为:[ 1.  1.]
    [ 1.  1.] 7.69713624862e-16 57

    迭代57次后得到解(1,1)

    reference
    梯度-牛顿-拟牛顿优化算法和实现
    牛顿法与拟牛顿法学习笔记(四)BFGS 算法
    最优化理论与方法-袁亚湘
    我的课本-最优化选讲

    展开全文
  • 作业二 用DFP算法求解取 一求解 求迭代点x1 令得的极小值点 所以得 于是由DFP修正公式有 下一个搜索方向为 求迭代点x2 令得的极小值点 于是得所以 因Hesse阵为正定阵为严格凸函数所以为整体 极小点 二DFP算法迭代...
  • DFP算法及Matlab程序 精品文档 精品文档 收集于网络如有侵权...所以得 于是由DFP修正公式有 下一个搜索方向为 求迭代点x2 令得的极小值点 于是得所以 因Hesse阵为正定阵为严格凸函数所以为整体 极小点 二DFP算法迭代步
  • 数学软件 之 基于MATLAB的DFP算法

    千次阅读 2019-10-08 19:24:56
    也是无约束最优化方法中非常重要的两个拟Newton算法之一,上一周写了一周的数学软件课程论文,姑且将DFP算法的实现细节贴出来分享给学弟学妹参考吧,由于博客不支持数学公式,所以就不累述算法原理及推导公式了。...

      DFP算法是本科数学系中最优化方法的知识,也是无约束最优化方法中非常重要的两个拟Newton算法之一,上一周写了一周的数学软件课程论文,姑且将DFP算法的实现细节贴出来分享给学弟学妹参考吧,由于博客不支持数学公式,所以就不累述算法原理及推导公式了。

     


     

     

     

     

    DFP算法流程图

    先给出DFP算法迭代流程图,总体上是拟Newton方法的通用迭代步骤,唯独在校正公式的地方有所区别。

     

     

     

    MATLAB实现DFP

      基于此图便可以设计DFP算法的MATLAB程序:

      对分法及加步探索法的实现

      首先由于DFP算法中需要利用一维搜索得到最优步长,因此需要先设计一个一维搜索函数,博主选用的是简单的对分法(二分法):

      

    %本函数利用二分法求解X = ls(Xk,Pk)问题
    %目标函数:f
    %符号参数:var
    %终止限:eps
    function x = dichotomy(f,var,eps)
        g = diff(f,var);
        [a, b] = search(f,var);
        x = (a + b)/2;      %防止eps过大导致x无值
        while b - a > eps
            x = (a+b)/2;
            gx = subs(g, var, x);
            if gx > 0
                b = x;
            elseif gx < 0
                a = x;
            else break;
            end
        end
        %加步搜索法-确定搜索区间
        function [a, b] = search(g,var)
            gt = matlabFunction(g);
            X = 0;  tmp = X;
            h = 1;  k = 0;
            while 1
                Xk = X + h; k = k+1;
                Y = subs(gt,var,X);
                Yk = subs(gt,var,Xk);
                if Y > Yk   %加大步长搜索
                    h = 2 * h;
                    tmp = X;
                    X = Xk;
                elseif Y == Yk  %缩小步长搜索
                    h = h/2;
                elseif k == 1
                    h = -h;     %反向搜索
                else break;
                end
            end
            a = min(tmp, Xk);
            b = max(tmp, Xk);
        end
    end
    

     

      DFP算法的实现

      有了一维搜索函数,那么实现DFP算法也就能依照算法流程图来设计了:

      

    %DFP算法主程序
    %目标函数:f
    %初始点:X0
    %参数:var
    %终止限:eps
    function DFP(f, X0, var, eps)
    %初始化符号函数,梯度,维数等
    syms var t;
    g = jacobian(f)';   %Jacobian转置->Grad
    fx = matlabFunction(f); %符号函数->函数句柄(R2009以上支持)
    gx = matlabFunction(g);
    n = length(var);    %维数
    X = X0; Xk = X0;
    while 1
        fx0 = fx(X(1),X(2));   gx0 = gx(X(1),X(2));
        Hk = eye(2);    Pk = -gx0;  %初始方向
        k = 0;  %迭代次数
        while 1
            Y = Xk + t*Pk;
            y = fx(Y(1),Y(2));
            tk = dichotomy(y, t, eps);  %一维搜索
            Xk = Xk + tk*Pk;
            fx1 = fx(Xk(1),Xk(2));
            gx1 = gx(Xk(1),Xk(2));
            if norm(gx1) < eps || k == n
                X = Xk; fx0 = fx1;
                break;
            end
            Sk = Xk - X;    Yk = gx1 - gx0;
            Hk = Hk + Sk*Sk'/(Sk'*Yk) - Hk*(Yk)*Yk'*Hk/(Yk'*Hk*Yk);
            Pk = -Hk*gx1;   %校正方向
            k = k+1;
        end
        if norm(gx1) < eps
            disp('X(k+1) = ');  disp(Xk);
            disp('F(K+1) = ');  disp(fx0);
            break;
        end
    end
    

     


     

     

    实例验证

     

      有了DFP算法的实现函数,那么应用于实例也就难了。

      可以在命令文件下输入如下代码就能得到目标函数极值点及极值

    clear; clc; format long;
    syms x1 x2;
    f = 4*(x1-5)^2 + (x2-6)^2;
    tic;    %初始时间
    DFP(f, [8;9], [x1, x2], 0.00000001);
    toc;    %结束时间
    

     

    输出结果如下:

    X(k+1) =

       4.999995811278565

       5.999767686222325

    F(K+1) =

         5.403987284687523e-08

    Elapsed time is 8.229108 seconds.

     

      算法时间度分析:

    由此可知,函数 [8,9]附近的点[5.00,6.00]处取得局部最小值,其中局部极值点约为5.40e-8.

    此算法运行时间约为8.23s,并且我们在降低终止限eps后,针对本题,算法运行时间增长较快,例如若eps = 1e-3,耗时11.6s,若eps = 1e-5,耗时22.94s,而eps = 1e-7,耗时甚至超过15分钟.这说明DFP算法在求解高精度运算时的运行效率表现得并不是那么好,甚至有可能无法得出最优解.

      

      实例搜索图

      基于该实例,对算法的迭代过程进行绘图,得到如下搜索图

      

    可以由以上两个搜索图像得出一个结论:DFP算法的实质是在每一次迭代过程中调整自己的搜索方向,以使得该方向能够尽可能接近极值点,这也正是几乎所有拟Newton算法中校正矩阵的作用.

     

     


     

     

    转载于:https://www.cnblogs.com/Inkblots/p/5432731.html

    展开全文
  • MATLAB拟牛顿法之DFP与BFGS算法

    千次阅读 多人点赞 2020-04-22 10:39:44
    由于博主使用WPS编辑的文本,公式无法赋值粘贴,这里以截图的方法给出了推导过程。博主会上传该DOC文档。 matlab代码 syms x1 x2 f=@(x1,x2) x1.^2+x2.^2-x1*x2-10*x1-4*x2+60; X=DFP(f,[0 0],1e-8,100) ...
  • 对非凸目标函数Broyden变尺度算法的收敛性是一个没有...针对DFP修正公式证明在不假定精确线搜索条件下,对光滑的目标函数,当DFP算法得到的点列收敛时,该点列一定趋向于稳定点。指出对于其他Broyden算法结论都是成立的。
  • DFP和BFGS迭代公式的证明
  • 研究了一种基于可变收敛因子的Davidon-Fletcher-Powell(DFP)自适应算法,给出算法中自相关逆矩阵估计的递归更新公式。对DFP算法中参数τ(n)的作用及算法的计算复杂度进行了分析。当分别输入正弦信号和高斯白噪声时,对...
  • 由于博主使用WPS编辑的文本,公式无法赋值粘贴,这里以截图的方法给出了推导过程。博主会上传该DOC文档。该资源为博客配套讲义资源。
  • 拟牛顿法中的DFP算法和BFGS算法

    万次阅读 2011-04-17 20:15:00
    注明:程序中调用的函数jintuifa.m golddiv.m我在之前的笔记中已贴出DFP算法和BFGS算法不同在于H矩阵的修正公式不同DFP算法%拟牛顿法中DFP算法求解f = x1*x1+2*x2*x2-2*x1*x2-4*x1的最小值,起始点为x0=[1 1] H0为...
  • DFP(x0,ess)  % ########################### syms x1 x2 t; f=0.5*2000*sqrt((3-x1)^2+(8-x2)^2)+0.5*3000*sqrt((8-x1)^2+(2-x2)^2)+0.75*2500*sqrt((2-x1)^2+(5-x2)^...
  • 其中E矩阵的求解公式为: 2.MATLAB可执行程序 function [k,x_min,f_min]= variable_metric_algorithm(f,x,exp) %% 程序说明 %{ 该程序应用于多维无条件优化问题中的变尺度法 变量说明: 输入值部分 f为目标函数,...
  • 牛顿法的特点就是收敛快。但是运用牛顿法需要计算二阶偏导数,而且目标函数的Hesse矩阵可能非正定。...牛顿法的迭代公式 x(k+1)=x(k)+λd(k)x^{(k+1)}=x^{(k)}+\lambda d^{(k)}x(k+1)=x(k)+λd(k)d(k)=−▽2f(x(k))...
  • BFGS算法与DFP步骤基本相同,区别在于更新公式的差异 ''' def bfgs(fun,gfun,hess,x0): #功能:用BFGS族算法求解无约束问题:min fun(x) #输入:x0是初始点,fun,gfun分别是目标函数和梯度...
  • 在Broyden凸族建立了Hesse近似矩阵关于目标函数梯度向量等内积分解矩阵的校正公式,从而把由校正矩阵的等内积分解矩阵确定搜索方向的DFP和BFGS算法推广到Broyden凸族.
  • 一直记不住这些算法的推导,所以打算详细点写到博客中以后不记得就翻阅自己的笔记。...上述的是当x是标量时的展开式,当x是多元时可以根据以下公式进行推导: 舍去二阶项以上的项可以得到: 参考文...
  • 牛顿法Taylor公式 – Maclaurin公式:1.求根问题:牛顿法的最初提出是用来求解方程的根的。我们假设点x∗为函数f(x)的根,那么有f(x∗)=0。现在我们把函数f(x)在点xk处一阶泰勒展开有: 那么假设点xk+1为该方程的根,...
  • BFGS是可以认为是由DFP算法推导出来的,上篇文章有详细的推导:(拟牛顿法公式推导以及python代码实现(一))目前BFGS被证明是最有效的拟牛顿优化方法。它的思想是根据我们已知的两个拟牛顿条件:关于Hess近似的:...
  • 优化算法——拟牛顿法之BFGS算法

    万次阅读 2015-05-20 11:31:14
    一、BFGS算法简介  BFGS算法是使用较多的... 同DFP校正的推导公式一样,DFP校正见博文“优化算法——拟牛顿法之DFP算法”。对于拟牛顿方程: 可以化简为: 令,则可得: 在B
  • 本节介绍: hessian matrix 近似 DFP算法 bfgs算法 hessian matrix 近似牛顿法的基本思路是用二次函数来局部逼近目标函数 ff 并解近似函数的极小点作为下一个迭代点,迭代公式 但是牛顿法的缺陷是需要
  • 在此基础上,提出了求解矩阵最大特征值及其对应特征向量的拟Newton法,给出求解矩阵最大特征值及其单位化向量重新整理后的Broyden方法公式、BFS方法公式DFP方法公式及其对应的Broyden算法,BFS算法,DFP算法。...
  • DFP BFGS 一.Taylor公式-Maclaurin公式 1.1 应用1:函数值计算 数值计算:初等函数值的计算(在原点展开) 计算: 求整数k和小数r,使得: x=kln2+r,|r|≤0.5ln2 1.2 应用2:解释Gini系数...
  • DFP – BFGS 一 Taylor 公式-Maclaurin公式 泰勒展开式可以在任意一点展开,即第一个式子。当在x=0处展开时称为maclaurin(麦克劳林)公式,即第二个式子。 第三步的约等于号变成等于号是为了方便计算。对于...
  • 机器学习中的数学 觉得有用的话,欢迎一起讨论相互学习~Follow Me 原创文章,如需转载请保留出处 本博客为七月在线邹博老师机器...DFP BFGS Taylor公式 如果函数在x0点可以计算n阶导数,则有Taylor展开 如果取x0...

空空如也

空空如也

1 2 3
收藏数 46
精华内容 18
关键字:

dfp公式