算法來自于[1],如下：
值得一提的是,[1]中的python代碼實現了對Rosenbrock函數的求極值測試.
例子來自于[2]:

----------------------------------------------------
用DFP算法求解：
$minf(x)=x_1^2+2x_2^2?2x_1{x}_2?4x_1$
取$x_0=(1,1{)}^T,H_0=$
$$\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$$

----------------------------------------------------

解答:
$g(x)=(2x_1?2x_2?4,?2x_1+4x_2{)}^T$
$g_0=(?4,2{)}^T$
$p_0=?H_0{g}_0=(4,?2{)}^T$,
(i)求迭代點$x_1$,令
${?}_0(\alpha )=f(x_0+\alpha p_0)=40{\alpha }^2?20\alpha ?3$,
得到$?(\alpha )$的極小值點為${\alpha }_0=\frac{1}{4}$,
所以得：
$x_1=x_0+{\alpha }_0{p}_0=(2,0.5{)}^T,g_1=(?1,?2{)}^T,$
$s_0=x_1?x_0=(1,?0.5{)}^T,y_0=g_1?g_0=(3,?4{)}^T$
這里的 $s_0$是因為需要滿足一個擬Newton條件，可以參考[4]

于是,
由DFP修正公式有$H_1=H_0?\frac{H_0{y}_0{y}_0^T{H}_0}{y_0^T{H}_0{y}_0}+\frac{s_0{s}_0^T}{y_0^T{s}_0}=\frac{1}{100}$$\left[\begin{array}{cc}84& 38\\ 38& 41\end{array}\right]$

所以下一個搜索方向為$p_1=?H_1{g}_1=\frac{1}{5}(8,6{)}^T$

(2)求迭代點x2
令
${?}_1(\alpha )=f(x_1+\alpha p_1)=\frac{8}{5}{\alpha }^2?4\alpha ?5.5$,

得到$?(\alpha )$的極小值點${\alpha }_1=\frac{5}{4}$
于是得:
$x_2=x_1+{\alpha }_1{p}_1=(4,2{)}^T,g_2=(0,0{)}^T,所以:x^{?}=x_2=(4,2{)}^T,f^{?}=?8$
因為Hessian矩陣G(x)=G=
${\left[\begin{array}{cc}2& ?2\\ ?2& 4\end{array}\right]}^T$
為正定矩陣，$f(x)$為嚴格凸函數，所以$x?$為整體極小點

[3]提供了matlab代碼,建立一個文件DFP.m(文件名必須和代碼中的函數名保持一致)，代碼如下：

function [best_x,best_fx,count]=DFP(x0,ess) colormap Jet % ########################### syms x1 x2 t; f=x1*x1+2*x2*x2-2*x1*x2-4*x1; fx=diff(f,x1);%求表達式f對x1的一階求導 fy=diff(f,x2);%求表達式f對x2的一階求導 fi=[fx fy];%構造函數f的梯度函數 %初始點的梯度和函數值 g0=subs(fi,[x1 x2],x0); f0=subs(f,[x1 x2],x0); H0=eye(2); %輸出x0,f0,g0 x0 f0 g0 xk=x0; fk=f0; gk=g0; Hk=H0; k=1; while(norm(gk)>ess)%迭代終止條件||gk||<=ess disp('************************************************************') disp(['第' num2str(k) '次尋優']) %確定搜索方向 pk=-Hk*gk'; %由步長找到下一點x(k+1) xk=xk+t*pk'; f_t=subs(f,[x1 x2],xk); %構造一元搜索的一元函數φ(t) %由一維搜索找到最優步長 df_t=diff(f_t,t); tk=solve(df_t); if tk~=0 tk=double(tk); elsebreak; end %計算下一點的函數值和梯度xk = subs(xk,t,tk) fk=subs(f,[x1 x2],xk) gk0=gk; gk=subs(fi,[x1 x2],xk) %DPF校正公式，找到修正矩陣 yk=gk-gk0; sk=tk*pk';Hk=Hk-(Hk*yk'*yk*Hk)/(yk*Hk*yk')+sk'*sk/(yk*sk')%修正公式 k=k+1; enddisp('結果如下：') best_x=xk;%最優點 best_fx=fk;%最優值 count=k-1; end

matlab終端運行方法如下:
>> x0=[1 1];
>> ess=1e-6
>> [best_x,best_fx,count]=DFP(x0,ess)

輸出如下：

x0 =1 1f0 =-3g0 =[ -4, 2]************************************************************ 第1次尋優xk =[ 2, 1/2]fk =-11/2gk =[ -1, -2]Hk =[ 21/25, 19/50] [ 19/50, 41/100]************************************************************ 第2次尋優xk =[ 4, 2]fk =-8gk =[ 0, 0]Hk =[ 1, 1/2] [ 1/2, 1/2]結果如下：best_x =[ 4, 2]best_fx =-8count =2>> x0=[1 1]; >> ess=1e-6ess =1.0000e-06

Reference:

[1]優化算法——擬牛頓法之DFP算法
[2]擬牛頓法－最優化方法－百度文庫
[3]DFP算法及Matlab程序
[4]DFP算法

總結

以上是生活随笔為你收集整理的拟牛顿法-DFP算法举例与matlab代码实现(转载＋整理)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。