【优化算法】信赖域折线（狗腿）法（Trust Region Dogleg Method）

最新推荐文章于 2024-01-18 23:14:24 发布

HeroOreo

最新推荐文章于 2024-01-18 23:14:24 发布

阅读量6.4k

点赞数 1

文章标签：算法 matlab 机器学习

本文链接：https://blog.csdn.net/weixin_43772929/article/details/121960254

版权

出发点

刘红英老师的最优化课程需要用到信赖域方法，求解问题并画图，网上找到的中文内容的矩阵实验室（MATLAB）信赖域折线（狗腿）法求解代码，都很繁琐，而且心理感觉不可靠（手动狗头），而且我没心思一点一点调试，最后不如自己写一个靠谱。

特点

信赖域法相比线搜索而言，最大的特点在于用少量的步骤，在一定约束范围内（通常是限定迭代步长），求出一个局部优化的值。

折线法算法

其中，函数dogleg就是用折线法求信赖域子问题，即二次近似函数步长约束最优问题，即

$\min q(s)=g^Ts+s^TBs,s.t.|s|≤\Delta$
$g, B$ 分别为 $x_k$ 点的梯度和（近似）二阶导， $B$ 不要求正定
$s$ 是下一步的最优迭代步长（下面用 $p$ 表示）

全局最优点 $p_B=-B^{-1}g$
沿负梯度方向的全局最优点 $p_U=-\cfrac{g^Tg}{g^TBg}g$
$若p_B在\Delta范围内，取p=p_B$
$若p_B和p_U都在\Delta范围外，取p=-\cfrac{\Delta}{|g|}g$
$若p_B在\Delta范围外，p_U在\Delta范围内，p的计算如下：$

$p_{BU}=p_B-p_U$
$\tau-1=\cfrac{\sqrt{(p_U^Tp_{BU})^2-||p_{BU}||^2(||p_U||^2-\Delta^2))}-p_U^Tp_{BU}}{||p_{BU}||^2}$
$p=p_U+(\tau-1)p_{BU}$

二次近似函数步长约束最优问题示例

备注

如果要优化自己的目标函数，只要把代码中的函数myf(x)的内容改掉就行，本代码理论上适用于任意维度的自变量
本程序包括了二阶导阵（海森阵Hessian）的模式选择，可以选择中心差分方法，也可以选BFGS矫正法。
函数Show中，默认不显示等高线图。如果需要，就自己提前准备好，放在同一文件夹里。

MATLAB代码

%%  main function Version 8信赖域法
%   创建：2021年12月14日11:23:22
%   末改：2021年12月15日19:40:54
%   继承：mainV6,BFGS拟牛顿法
%   修改：信赖域法
%   总结：别人写的东西都太复杂太不靠谱了

%%  Initialize
clear;clc;close all;tic
% 函数调用次数统计：
global kFcn % 全局变量，统计函数@myf的被调用次数
kFcn=0;

% 参数设置：
err=1e-6;h=1e-6;MaxIter=1000;
Dta=0.1;  Dtamax=1.0;
eta1=0.1; eta2=0.75;
gma1=0.5; gma2=2.0;
% HessianMode="Central finite difference";
% HessianMode="Diff Positive Definite";
HessianMode="BFGS";

% 初始点:
% x0=[0.7395, 0.3144]';
% x0=[0.8, 0.4]';
x0=[1,0.5]';
len=length(x0);
k=1;xk=x0;
x_store=zeros(len,MaxIter);
x_store(:,k)=xk;
fprintf("x%d=%f,%f\n",k,xk);
    
%计算函数值和梯度：
fk=myf(xk);
gk=grd(xk,h);
% Hk为海森阵（二阶导阵）Gk的逆：
Gk=Hess(xk,h);Bk=Gk;Hk=inv(Bk);

%% Iter进入迭代循环
while(k<MaxIter)    %不能使用for循环计数因为中间有的循环被舍弃了
    
    if (norm(gk)<err);fprintf("uhu!!reach the end!\n");break;end
    
    % 折线法（Dogleg法，狗腿法）求子问题：
    % 二次近似函数圆约束最优问题，即
    % min q(s)=g·s+½s·B·s, s.t.|s|≤△
    sk=dogleg(gk,Bk,Dta);
    
    % 计算下降效果ρk=δf/δq：
    deltaq=-(gk'*sk+.5*sk'*Bk*sk);
    deltaf=myf(xk)-myf(xk+sk);
    rhok=deltaf/deltaq;
    
    % 根据下降效果ρk，进行信赖域分析：
    if rhok<=0||imag(rhok);Dta=gma1*Dta;continue;end %不符，缩步长
    if rhok<=eta1;Dta=gma1*Dta;end                   %凑合，缩步长
    if rhok>=eta1;Dta=gma2*Dta;Dta=min(Dta,Dtamax);end %好，扩步长
    
    % 更新迭代点并储存：
    xk=xk+sk;
    k=k+1;
    x_store(:,k)=xk;
    fprintf("x%d=%f,%f, ρ%d=%f, Δ=%f\n",k,xk,k,rhok,Dta);
    
    % 计算函数值和梯度：
    fk=myf(xk);
    gk1=grd(xk,h);
    yk=gk1-gk;
    gk=gk1;
    %% 计算二阶导阵（海森阵Gk）的近似阵Bk：
    if HessianMode=="central difference"
        Gk=Hess(xk,h);
        Bk=Gk;
    elseif HessianMode=="Diff Positive Definite"
        Gk=Hess(xk,h);
        while(min(eig(Gk))<0);Gk=Gk+2*norm(eig(Gk))*eye(2);end%正定化
        Bk=Gk;
    elseif HessianMode=="BFGS"
        Hk=Hk+(1/(sk'*yk)+yk'*Hk*yk/(sk'*yk)^2)*(sk*sk')...
            -1/(sk'*yk)*(Hk*yk*sk'+sk*yk'*Hk);
        Bk=Bk-Bk*(sk*sk')*Bk/(sk'*Bk*sk)+yk*yk'/(yk'*sk);
    end
end

toc

%% show
Show(x_store,k);

%% Functions
function f=myf(x)       % 目标函数
global kFcn;kFcn=kFcn+1;% 全局变量，统计函数@myf的被调用次数
x1=x(1);x2=x(2);
a=[0.3;0.6;0.2];b=[5;26;3];c=[40,1,10];
u=x1-0.8;
v=x2-(a(1)+a(2)*u^2*(1-u)^.5-a(3)*u);
alpha=-b(1)+b(2)*u^2*(1+u)^.5+b(3)*u;
beta=c(1)*v^2*(1-c(2)*v)/(1+c(3)*u^2);
f=alpha*exp(-beta);
end

function gf=grd(x,h)    % 中心查分求梯度
len_x=length(x);
g=zeros(len_x,1);
for i_x=1:len_x
    xp=x;xm=x;
    xp(i_x)=xp(i_x)+h;
    xm(i_x)=xm(i_x)-h;
    g(i_x)=(myf(xp)-myf(xm))/h/2;
end
gf=g;
end

function Hs=Hess(x,h)   % 中心查分求二阶导阵
% len_x=2;
H=zeros(length(x));
for i_x=1:length(x)          % Form the F
    xp=x;xm=x;
    xp(i_x)=xp(i_x)+h;
    xm(i_x)=xm(i_x)-h;
    % fprintf("dj=\t%.8f\t%.8f\t\n",dj);
    resp=grd(xp,h);resm=grd(xm,h);
    H(:,i_x)=(resp-resm)/2/h;
end
Hs=(H+H')/2;
% eig_H=eig(Hs);
% while(eig_H(1)<0||eig_H(2)<0)
%     Hs=Hs+0.5*eye(2);eig_H=eig(Hs);
% end
end

function sk=dogleg(gk,Bk,Dta)
% 折线法（Dogleg法，狗腿法）求子问题：
% 二次近似函数圆约束最优问题，即
% min q(s)=g·s+½s·B·s, s.t.|s|≤△
pB=-Bk\gk;                  % 全局最优点
pU=-gk'*gk/(gk'*Bk*gk)*gk;  % 最降向全局最优点

if norm(pB)<=Dta;sk=pB;     %fprintf("pB ");
elseif norm(pU)>Dta;sk=-Dta/norm(gk)*gk;%fprintf("pU ");
else; pBU=pB-pU;        dotp=pU'*pBU;
    tau_1=((dotp^2-pBU'*pBU*(pU'*pU-Dta^2))^.5-dotp)/(pBU'*pBU);
    sk=pU+tau_1*pBU;        %fprintf("dogleg ");
end

%     fprintf("|pB|=%f,|pU|=%f,△=%f\n",norm(pB),norm(pU),Dta);
end 

function  Show(x_store,k)
global kFcn;
fprintf('x%d=%f,%f\n',k,x_store(:,k));
fprintf('总调用函数次数：%d\n',kFcn);

% Point Sequence
figure
x_show=x_store(:,1:k);
k1_show=x_show(1,:)';k2_show=x_show(2,:)';% column vector
plot(k1_show,k2_show,'o-')
PSeqTxt=num2str((1:length(x_show))');
text(k1_show,k2_show,PSeqTxt);
axis equal
xlabel('k1');ylabel('k2');

% Contour
% load ContourData
% hold on;
% contour(x,y,z,20,'ShowText','off')


% toc;
end

结果：

x1=1.000000,0.500000
x2=0.982381,0.401564, ρ2=0.984471, Δ=0.200000
x3=0.893268,0.222515, ρ3=0.305689, Δ=0.400000
x4=0.784681,0.325616, ρ4=0.907275, Δ=0.800000
x5=0.742343,0.299018, ρ5=0.457107, Δ=0.100000
x6=0.738332,0.315455, ρ6=0.939598, Δ=0.200000
x7=0.739620,0.314362, ρ7=1.011037, Δ=0.400000
x8=0.739479,0.314362, ρ8=0.792868, Δ=0.800000
x9=0.739505,0.314360, ρ9=0.999604, Δ=1.000000
x10=0.739505,0.314360, ρ10=0.987830, Δ=1.000000
uhu!!reach the end!
时间已过 0.057087 秒。
x10=0.739505,0.314360
总调用函数次数：92

请添加图片描述

参考

刘红英，夏勇，周水生，数学规划基础（研究生教材），北京航空航天大学出版社，2012.10

信赖域算法-The Dogleg Method（含例题及Python实现）-CSDN等待戈多

其他

引用请标注来源。

HeroOreo

关注

1
点赞
踩
35

收藏

觉得还不错? 一键收藏
5
评论
【优化算法】信赖域折线（狗腿）法（Trust Region Dogleg Method）

刘红英老师的最优化课程需要用到信赖域方法，求解问题并画图，网上找到的中文内容的矩阵实验室（MATLAB）信赖域折线（狗腿）法求解代码，都很繁琐，而且心理感觉不可靠（手动狗头），而且我没心思一点一点调试，最后不如自己写一个靠谱。
复制链接

扫一扫