信赖域狗腿（dogleg）方法

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

置顶陆嵩

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

阅读量2.2w

点赞数 17

分类专栏：计算数学数据科学与人工智能 MATLAB与Mathematica 数学原理文章标签：最优化非线性规划无约束优化信赖域方法狗腿法

本文链接：https://blog.csdn.net/lusongno1/article/details/80558684

版权

计算数学同时被 3 个专栏收录

167 篇文章 153 订阅

订阅专栏

MATLAB与Mathematica

141 篇文章 103 订阅

订阅专栏

数学原理

116 篇文章 17 订阅

订阅专栏

信赖域狗腿方法

信赖域方法（Trust-region methods）又称为TR方法，它是一种最优化方法，能够保证最优化方法总体收敛。现今，信赖域算法广泛应用于应用数学、物理、化学、工程学、计算机科学、生物学与医学等学科。相信在不远将来，信赖域方法会在更广泛多样的领域有着更深远的的发展。本文简单介绍信赖域狗腿（dogleg）方法，希望能为其他专业诸如机器学习的同学在使用优化工具时提供一个参考。

原理简述

考虑模型问题：
这里写图片描述

它具有全局最优解：
这里写图片描述

该问题沿着负梯度方向，具有全局最优解：
这里写图片描述

由于信赖域的限制，全局最优解可能落在信赖域之外。我们依据两个全局最优点是否落于信赖域中，来决定下降方向。

当B点在信赖域内部时，取方向为 $p^U$ 方向， $x_{k+1}$ 点为B点：
当B点和U点都在信赖域外时，取下一个迭代点为 $p^U$ 和信赖域边界的交点：
当B点在外，U点在内时，去BU连线和信赖域边界的交点：

我们可以用一个参数 $\tau$ 将这三种情况统一为一个表达式，即下降取值为：
这里写图片描述

这里的 $\tau$ 是如何取的呢？通过简单的数学推导，我们知道它可以这样来取，刚好满足上面的三种情况。
这里写图片描述

简单地说，第一种情况， $\tau$ 取2，第二种情况， $\tau$ 取为信赖域半径和 $p^U$ 长的比值，第三种情况，要求那个交点，我们就要求一元二次方程的一个根，这可以用求根公式算得。

算法步骤

总的算法框架如下：
这里写图片描述

$\tau$ 和 $s_k = x_{k+1} - x_k$ 的计算就是使用上面提到的狗腿方法。

下降比 $\rho_k$ 的计算表达式为：
这里写图片描述

所以，到这事情就很明了了。

代码

function [x,k] = trust_region1(varargin)
%帮助信息：
%clc
%clear
%输入参数为options = struct('f',f,'x0',[0;0],'Delta',1,'Delta0',0.1,'eta',1/4);
if nargin < 1, help(mfilename),
    f = @(x1,x2)100*(x2-x1^2)^2 + (1-x1)^2;
    options = struct('f',f,'x0',[0;0],'Delta_hat',10,'Delta0',0.1,'eta',0);
    fprintf('使用默认参数…… \n \n');
    addpath(genpath(pwd));
end

f = options.f;
syms x1 x2;
f_s = f(x1,x2);
x0 = options.x0;
Delta_hat = options.Delta_hat;
Delta0 = options.Delta0;
eta = options.eta;
df = diff_handle(f_s);
B = hessian(f_s);
Delta = Delta0;
x = x0;
step = 100;
for k = 0:step-1
    x
    fk = f(x(1),x(2));
    dfk = df(x(1),x(2));
    Bk = subs(B,{x1,x2},{x(1),x(2)});
    Bk = double(Bk);
    sk = cal_sk(dfk,Bk,Delta);
    m = @(p) fk + dfk'*p + 1/2*p'*Bk*p;
    rho_k = cal_rho_k(f,m,x,sk);
    if rho_k < 1/4
        Delta = 1/4*Delta;
    elseif rho_k > 3/4 && abs(norm(sk,2) - Delta)<1.0e-10;
        Delta = min(2*Delta,Delta_hat);
    else
        %不做任何操作；
    end

    if rho_k > eta
        x = x + sk;
    else
        %不做任何操作；
    end
end

end








function df = diff_handle(f_s)
syms x1 x2;
df = [diff(f_s,x1); diff(f_s,x2)];
df = matlabFunction(df);
end



function tau = cal_tau(pB,pU,Delta)
npB = sqrt(pB'*pB);
npU = sqrt(pU'*pU);
if npB <= Delta
    tau = 2;
elseif npU >= Delta
    tau = Delta/npU;
else
    pB_U = pB-pU;
    tau = (-pU'*pB_U+sqrt((pU'*pB_U)^2-pB_U'*pB_U*(pU'*pU-Delta^2)))/(pB_U'*pB_U);
    tau = tau + 1;
end
end
% temp = conj(dfk')*Bk*dfk;
% if  temp<= 0
%     tau = 1;
% else
%     tau = min(norm(dfk,2)/temp,1);




function sk = cal_sk(dfk,Bk,Delta)
pU = -(conj(dfk')*dfk)/(conj(dfk')*Bk*dfk).*dfk;
pB = -Bk^(-1)*dfk;
tau = cal_tau(pB,pU,Delta);
if tau >=0 && tau <=1
    sk = tau*pU;
elseif tau >= 1 && tau <=2
    sk = pU + (tau-1)*(pB-pU);
else
    error('tau的值不能为%f',tau);
end

end


function rho_k = cal_rho_k(f,m,x,sk)
rho_k = (f(x(1),x(2)) - f(x(1)+sk(1),x(2)+sk(2)))/((m([0;0])-m(sk)));
end