數學軟件 之 基於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算法中校正矩陣的作用.

 

 

---