数学软件之基于MATLAB的DFP算法

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

DFP算法流程图

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

数学软件之基于MATLAB的DFP算法

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;    %结束时间

Normal
0

7.8 磅
0
2

false
false
false

EN-US
ZH-CN
X-NONE

/* Style Definitions */
table.MsoNormalTable
{mso-style-name:普通表格;
mso-tstyle-rowband-size:0;
mso-tstyle-colband-size:0;
mso-style-noshow:yes;
mso-style-priority:99;
mso-style-parent:"";
mso-padding-alt:0cm 5.4pt 0cm 5.4pt;
mso-para-margin:0cm;
mso-para-margin-bottom:.0001pt;
mso-pagination:widow-orphan;
font-size:10.0pt;
font-family:"Times New Roman",serif;}

输出结果如下：

X(k+1) =

4.999995811278565

5.999767686222325

F(K+1) =

5.403987284687523e-08

Elapsed time is 8.229108 seconds.

　　算法时间度分析：

Normal
0

7.8 磅
0
2

false
false
false

EN-US
ZH-CN
X-NONE

由此可知，函数在[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算法在求解高精度运算时的运行效率表现得并不是那么好，甚至有可能无法得出最优解.

　　实例搜索图

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

　　数学软件之基于MATLAB的DFP算法

Normal
0

7.8 磅
0
2

false
false
false

EN-US
ZH-CN
X-NONE

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

Normal
0

7.8 磅
0
2

false
false
false

EN-US
ZH-CN
X-NONE

秒客网

数学软件之基于MATLAB的DFP算法

相关文章

数学软件 之 基于MATLAB的DFP算法

相关文章

数学软件之基于MATLAB的DFP算法