拟牛顿法BFDS，matlab代码实现

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

deer(écho)

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

阅读量1.9k

点赞数

分类专栏： matlab 最优化算法文章标签： matlab 线性代数算法

本文链接：https://blog.csdn.net/deercho/article/details/121870728

版权

算法同时被 3 个专栏收录

6 篇文章 0 订阅

订阅专栏

matlab

2 篇文章 0 订阅

订阅专栏

最优化

2 篇文章 0 订阅

订阅专栏

该博客详细介绍了牛顿法和拟牛顿法，特别是BFGS校正策略在迭代优化过程中的应用。通过MATLAB代码展示了BFGS方法的实现，包括搜索方向的计算、Hesse矩阵近似的更新以及迭代终止条件。内容涵盖了梯度、目标函数和曲率近似的构建，讨论了方法的下降性质和正定性，以及如何通过秩一校正来更新Hesse逆矩阵的近似。

摘要由CSDN通过智能技术生成

牛顿法的关键就是利用了Hesse的曲率信息，但是Hesse的计算和存储都很困难，那可否用梯度和目标函数这些信息来构造曲率近似呢？

实际上就是用B_k代替牛顿法中的Hesse矩阵，H_k替牛顿法中的Hesse逆矩阵

	牛顿法	拟牛顿法
迭代公式	x_k+1= x_k- α_k(▽²f(x_k))^-1 ▽f(x_k))	x_k+1 = x_k - α_kB_k^-1 . ▽f(x_k)
d_k 搜索方向	- (▽²f(x_k))^-1 ▽f(x_k))	- B_k^-1 . ▽f(x_k)

所以就是搜索方向用 - B_k^-1 . ▽f(x_k) 迭代，
初始时B_k为单位矩阵，搜索方向为最速方向

性质

仅需一阶导数
B_k 正定即可使该方法有下降性质
搜索方向互相共轭

H_k更新

x_k+1 = x_k - α_kB_k^-1 . g_k

s_k = x_k+1 - x_k，y_k = g_k+1 - g_k

那么如何更新Hesse逆的近似——H_k呢？我们提出了两种修正方法。

DFP校正

在这里插入图片描述

BFGS校正

两次运用秩一校正得到关于B_k的迭代公式

$[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-VOZBr6HY-1639193111371)(C:\Users\chl\Desktop\最优化笔记.assets\image-20211204172908437.png)]$
matlab代码实现

function [best_x,best_fx,count]=BFGS(x0,ess) 

% ###########################
digits(5);
syms x1 x2 t;  
f=100*(x1^2-x2)^2+(x1-1)^2;
fx=diff(f,x1);%求表达式f对x1的一阶求导  
fy=diff(f,x2);%求表达式f对x2的一阶求导 
fi=[fx fy];%构造函数f的梯度函数 
%初始点的梯度和函数值  
g0=subs(fi,[x1 x2],x0); 
g0=vpa(g0,5);
f0=subs(f,[x1 x2],x0); 
B0=hessian(f,[x1,x2]); %B0是二阶梯度
B0=subs(B0,[x1,x2],x0);
x0
f0 
g0 
B0
xk=x0; 
fk=f0; 
gk=g0; 
Bk=B0; 
k=1;  
while(norm(gk)>ess)%迭代终止条件||gk||<=ess  
    xk0=xk;
%     X = sprintf('x%d=',k);
%         disp(X);
%         disp(xk0);
    disp('************************************************************')     
        disp(['第' num2str(k) '次寻优']) 
%确定搜索方向   
        pk=-inv(Bk)*gk'; 
%由步长找到下一点x(k+1)  
        xk=xk+t*pk';     
        f_t=subs(f,[x1 x2],xk); %构造一元搜索的一元函数φ(t) %由一维搜索找到最优步长  
        %更新原函数为关于步长t的函数，用以求解最优步长
        df_t=diff(f_t,t);%关于t求导    
        tk=solve(df_t);% 计算导数为0处t的值 
        X = sprintf('步长为 %f',tk);
        disp(X);

%计算下一点的函数值和梯度
        xk = subs(xk,t,tk(1));    
        fk=subs(f,[x1 x2],xk);    
        gk0=gk;     
        gk=subs(fi,[x1 x2],xk); 
%BFGS校正公式，找到修正矩阵    
        yk=gk'-gk0';    
        sk=xk'-xk0';
%       Bk=Bk+(yk*yk')/(sk'*yk); %修正公式    
        Bk=Bk+(yk*yk')/(sk'*yk)-(Bk*sk*sk'*Bk)/(sk'*Bk*sk); %修正公式 
        k=k+1; 
        disp(Bk);
xk
fk 
gk 
Bk
end

disp('结果如下：')
best_x=xk;%最优点 
best_fx=fk;%最优值 
count=k-1; 
end