拟牛顿法-DFP算法举例与matlab代码实现(转载＋整理)

最新推荐文章于 2024-07-17 02:05:32 发布

微电子学与固体电子学-俞驰

最新推荐文章于 2024-07-17 02:05:32 发布

阅读量1.8w

点赞数 24

分类专栏：机器学习算法

本文链接：https://blog.csdn.net/appleyuchi/article/details/97395358

版权

机器学习算法专栏收录该内容

87 篇文章 7 订阅

订阅专栏

算法来自于[1],如下：在这里插入图片描述
值得一提的是,[1]中的python代码实现了对Rosenbrock函数的求极值测试.
例子来自于[2]:

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

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

解答:
$g(x)=(2x_1-2x_2-4,-2x_1+4x_2)^T$
$g_0=(-4,2)^T$
$p_0=-H_0g_0=(4,-2)^T$ ,
(i)求迭代点 $x_1$ ,令
$\phi_0(\alpha)=f(x_0+\alpha p_0)=40\alpha^2-20\alpha-3$ ,
得到 $\phi(\alpha)$ 的极小值点为 $\alpha_0=\frac{1}{4}$ ,
所以得：
$x_1=x_0+\alpha_0p_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^TH_0}{y_0^TH_0y_0}+\frac{s_0s_0^T}{y_0^Ts_0}=\frac{1}{100}$ $\left[ \begin{matrix} 84 & 38 \\ 38&41 \\ \end{matrix} \right]$

所以下一个搜索方向为 $p_1=-H_1g_1=\frac{1}{5}(8,6)^T$

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

得到 $\phi(\alpha)$ 的极小值点 $\alpha_1=\frac{5}{4}$
于是得:
$x_2=x_1+\alpha_1p_1=(4,2)^T,g_2=(0,0)^T,所以:x^*=x_2=(4,2)^T,f^*=-8$
因为Hessian矩阵G(x)=G=
$\left[ \begin{matrix} 2 & -2 \\ -2&4 \\ \end{matrix} \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); 
else
    break; 
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; 
end

disp('结果如下：')
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     1

 
f0 =
 
-3
 
 
g0 =
 
[ -4, 2]
 
************************************************************
第1次寻优
 
xk =
 
[ 2, 1/2]
 
 
fk =
 
-11/2
 
 
gk =
 
[ -1, -2]
 
 
Hk =
 
[ 21/25,  19/50]
[ 19/50, 41/100]
 
************************************************************
第2次寻优
 
xk =
 
[ 4, 2]
 
 
fk =
 
-8
 
 
gk =
 
[ 0, 0]
 
 
Hk =
 
[   1, 1/2]
[ 1/2, 1/2]
 
结果如下：
 
best_x =
 
[ 4, 2]
 
 
best_fx =
 
-8
 

count =

     2

>> x0=[1 1];
>> ess=1e-6

ess =

   1.0000e-06

Reference:

[1]优化算法——拟牛顿法之DFP算法
[2]拟牛顿法－最优化方法－百度文库
[3]DFP算法及Matlab程序
[4]DFP算法