数学软件 之 基于MATLAB的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

发布了0 篇原创文章 · 获赞 12 · 访问量 5万+
展开阅读全文
评论将由博主筛选后显示,对所有人可见 | 还能输入1000个字符

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

©️2019 CSDN 皮肤主题: 编程工作室 设计师: CSDN官方博客

分享到微信朋友圈

×

扫一扫,手机浏览