信赖域狗腿方法

信赖域方法（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

例子

使用狗腿方法来求解信赖域子问题，应用这个方法来最小化下述函数。这里的$B_k = \nabla^2 f(x_k)$。使用信赖域方法的更新规则更新。这里给出前两步迭代结果。

选定的参数如下所示：

这里的$Delta\_hat$表示信赖域上界，$Delta0$表示信赖域初始半径，$eta$即为$\eta$的值。

在此参数下，迭代的前两个$x$值为：