信赖域方法解优化问题

考虑无约束光滑优化问题：

$\min_x \;\; f(x).$

在第 $k$ 步迭代，信赖域方法对目标函数 $f(x)$ 在 $x=x^k$ 附近用如下二次函数逼近： $\displaystyle m_k(d)=f(x^k)+\nabla f(x^k)^\top d+\frac{1}{2}d^\top B^k d$ ，其中 $B^k$ 是一个对称矩阵，可以是精确海瑟矩阵或者海瑟矩阵的近似。同时，令 $\|d\|\le\Delta_k$ 限制区域大小以保证该二次函数的逼近质量。结合两者，信赖域子问题可以写为：

$d^k=\arg\min_{d}m_k(d),\quad \mathrm{s.t.}\quad \|d\|\le\Delta_k.$

在得到 $d^k$ 之后，计算 $\hat{x}^{k+1} = x^k+d^k$ 和比值 $\rho_k$ : $\displaystyle\rho_k=\frac{f(x^k)-f(\hat{x}^{k+1})}{m_k(0)-m_k(d^k)}$ 来确定是否更新当前迭代。如果 $\rho_k$ 较大，我们认为 $m_k$ 在 $x^k$ 附近逼近 $f(x)$ 的质量较好，更新 $x^{k+1} = \hat{x}^{k+1}$ ，并增大信赖域半径 $\Delta_k$ （如乘以一个大于1的常数）。否则，则拒绝更新，即令 $x^{k+1} = x^k$ ，并减小 $\Delta_k$ （如乘以一个小于1的数）。

这里，对于信赖域子问题的求解，我们采用截断共轭梯度法，见截断共轭梯度法。

初始化和迭代准备

输入信息：初始迭代点 x ，目标函数值和梯度计算函数 fun ，海瑟矩阵 hess 以及算法求解相关参数的结构体 opts 。

输出信息：迭代得到的解 x 和迭代信息结构体 out 。

out.msg ：标记停机时算法的求解状态
out.nrmG ：记录迭代过程的梯度范数
out.nrmg ：记录迭代退出时的梯度范数
out.iter ：记录迭代退出时的迭代步数
out.f ：记录迭代退出时的函数值
out.nfe ：调用原函数的次数

function [x, out] = fminTR(x,fun, hess, opts, varargin)

检查输入信息，MATLAB 以 nargin 表示函数输入的参数个数。当参数量不足 3 时报错，等于 3 时认为 opts 为空结构体，即全部采用默认参数。

从输入的结构体中读取参数或采取默认参数。

opts.ftol ：针对函数值的停机准则
opts.gtol ：针对梯度范数的停机准则
opts.eta1 ： $\rho^k$ 的下界（当超出此界时意味着信赖域半径需要缩小并拒绝更新）
opts.eta2 ： $\rho^k$ 的上界（当超出此界时意味着信赖域半径需要增大并接受更新）
opts.gamma1 ：每次调整信赖域半径缩小的比例
opts.gamma2 ：每次调整信赖域半径增大的比例
opts.maxit ：主循环的最大迭代次数
opts.record ：是否需要打印迭代信息
opts.verbose ：是否需要打印半径调整信息
opts.itPrint ：每隔几步输出一次迭代信息

if (nargin < 3); error('[x, out] = fminTR(fun, x, opts)'); end
if (nargin < 4); opts = []; end
if ~isfield(opts,'gtol');     opts.gtol = 1e-6;   end
if ~isfield(opts,'ftol');     opts.ftol = 1e-12;  end
if ~isfield(opts,'eta1');     opts.eta1 = 1e-2;   end
if ~isfield(opts,'eta2');     opts.eta2 = 0.9;    end
if ~isfield(opts,'gamma1');   opts.gamma1 = .25;  end
if ~isfield(opts,'gamma2');   opts.gamma2 = 10;   end
if ~isfield(opts,'maxit');    opts.maxit = 200;   end
if ~isfield(opts,'record');   opts.record = 0;    end
if ~isfield(opts,'itPrint');  opts.itPrint = 1;   end
if ~isfield(opts,'verbose');  opts.verbose = 0;    end

从结构体中复制参数。

maxit   = opts.maxit;   record  = opts.record;  itPrint = opts.itPrint;
gtol    = opts.gtol;    ftol    = opts.ftol;
eta1    = opts.eta1;    eta2    = opts.eta2;
gamma1  = opts.gamma1;  gamma2  = opts.gamma2;

迭代准备，计算初始点 $x$ 处的函数值和梯度。

out = struct();
[f,g] = fun(x);
out.nfe = 1;
nrmg = norm(g,2);
out.nrmG = nrmg;
fp = f; gp = g;

信赖域子问题利用截断共轭梯度法求解，利用结构体 opts_tCG 为其提供参数。

opts_tCG = struct();

初始化信赖域半径。初始值设定为提供的参数值或默认值 $\sqrt{\mathrm{len}(x)}/8$ 并令 $\Delta$ 的上界 $\bar{\Delta}$ 为 $\sqrt{\mathrm{len}(x)}$ 。

Delta_bar = sqrt(length(x));
Delta0 = Delta_bar/8;
if isfield(opts,'Delta')
    Delta = opts.Delta;
else
    Delta = Delta0;
end

两个参数分别记录信赖域半径连续增大或减小的次数，以方便初始值的调整。

consecutive_TRplus = 0;
consecutive_TRminus = 0;

将 tCG.m 存放的目录加入工作目录。

addpath('../newton');

当需要详细输出时，设定输出格式。

if record >= 1
    if ispc; str1 = '  %10s'; str2 = '  %8s';
    else     str1 = '  %10s'; str2 = '  %8s'; end
    stra = ['%5s', str1, str2, str2, str2, str2, str2, str2, '\n'];
    str_head = sprintf(stra, ...
        'iter', 'F', 'fdiff', 'mdiff', 'redf', 'ratio', 'Delta', 'nrmg');
    fprintf('%s', str_head);
    str_num = ['%4d  %+8.7e  %+2.1e  %+2.1e  %+2.1e  %+2.1e  %2.1e  %2.1e'];
end

迭代主循环

以 maxit 为最大迭代次数。

for iter = 1:maxit

截断共轭梯度法中用于判断收敛的参数。

    opts_tCG.kappa = 0.1;
    opts_tCG.theta = 1;

调用 tCG 函数，利用截断共轭梯度法求解信赖域子问题，得到迭代方向 $d$ 。 stop_tCG 表示截断共轭梯度法的退出原因。

    hess_tCG = @(d) hess(x,d);
    [d, Hd, out_tCG] = tCG(x, gp, hess_tCG, Delta, opts_tCG);
    stop_tCG = out_tCG.stop_tCG;

计算比值

$\displaystyle\rho_k=\frac{f(x^k)-f(x^k+d^k)}{m_k(0)-m_k(d^k)},$

以确定是否更新迭代和修正信赖域半径。首先计算 $m_k(0)-m_k(d^k)=-\nabla f(x^k)^\top d^k + \frac{1}{2}(d^k)^\top B^k d^k$ ，为了保证数值稳定性，增加一个小常数 rreg 。

    mdiff = g'*d + .5*d'*Hd;
    rreg = 10*max(1, abs(fp)) * eps;
    mdiff = -mdiff + rreg;
    model_decreased = mdiff > 0;

构造试验点 $\hat{x}^{k+1}=x^k+d^k$ ，并计算该点处计算函数值和梯度。

    xnew = x + d;
    [f,g] = fun(xnew);
    nrmg = norm(g,2);
    out.nfe = out.nfe + 1;

计算 $f(x^k)-f(x^k+d^k)$ ，同样地增加一个小常数 rreg。然后计算 $\rho^k$ 作为修正信赖域半径和判断是否更新的依据。

    redf = fp - f + rreg;
    ratio = redf/mdiff;

当 $\rho^k>\eta_1$ 时，接受此次更新，并记录上一步的函数值和梯度。

    if ratio >= eta1
        x = xnew;   fp = f;  gp = g;
        out.nrmG = [out.nrmG; nrmg];
    end

计算函数值相对变化。

    fdiff = abs(redf/(abs(fp)+1));

停机准则：当满足（1）梯度范数小于阈值或（2）函数值的相对变化小于阈值且 $\rho^k>0$ 时，停止迭代。

    cstop = nrmg <= gtol || (abs(fdiff) <= ftol && ratio > 0);

    % 当需要详细输出时，在（1）开始迭代时（2）达到收敛时（3）达到最大迭代次数退出迭代时
    % （4）每若干步时，打印详细结果。
    if record>=1 && (cstop || ...
            iter == 1 || iter==maxit || mod(iter,itPrint)==0)
        if mod(iter,20*itPrint) == 0 && iter ~= maxit && ~cstop
            fprintf('\n%s', str_head);
        end

stop_tCG 记录内层的截断共轭梯度法退出的原因，分别对应输出如下：

        switch stop_tCG
            case{1}
                str_tCG = ' [negative curvature]\n';
            case{2}
                str_tCG = ' [exceeded trust region]\n';
            case{3}
                str_tCG = ' [linear convergence]\n';
            case{4}
                str_tCG = ' [superlinear convergence]\n';
            case{5}
                str_tCG = ' [maximal iteration number reached]\n';
            case{6}
                str_tCG = ' [model did not decrease]\n';
        end

        fprintf(strcat(str_num,str_tCG), ...
            iter, f, fdiff, mdiff, redf, ratio, Delta, nrmg);
    end

当满足停机准则时，记录达到最优值，退出循环。

    if cstop
        out.msg = 'optimal';
        break;
    end

信赖域半径的调整

如果 $\rho^k<\eta_1$ 或者算法非降（或者当 $\rho^k$ 计算时出现分母约为 0 ，即 ${m_k(0)-m_k(d^k)}\approx 0$ ）时，不接受当前步迭代。这种情况下，需要对信赖域半径进行缩减。

具体而言，令 $\Delta \leftarrow \gamma_1\Delta$ ，将信赖域半径的连续增大次数置零，缩小次数加一。

    if ratio < eta1 || ~model_decreased || isnan(ratio)
        Delta = Delta*gamma1;
        consecutive_TRplus = 0;
        consecutive_TRminus = consecutive_TRminus + 1;

当信赖域半径连续 5 次减小时，认为当前的信赖域半径过大，并输出相应的提示信息。

        if consecutive_TRminus >= 5 && opts.verbose >= 1
            consecutive_TRminus = -inf;
            fprintf(' +++ Detected many consecutive TR- (radius decreases).\n');
            fprintf(' +++ Consider decreasing options.Delta_bar by an order of magnitude.\n');
            fprintf(' +++ Current values: options.Delta_bar = %g and options.Delta0 = %g.\n', options.Delta_bar, options.Delta0);
        end

当 $\rho_k>\eta_2$ 且 $\|d_k\|=\Delta$ （对应为截断共轭梯度法因遇到负曲率或者超出信赖域半径而终止），增大信赖域半径。 $\Delta \leftarrow \min\{\gamma_2\Delta, \bar{\Delta}\}$ （ $\bar{\Delta}$ 为信赖域半径的上界，设置为 $\sqrt{\mathrm{len}(x)}$ ）。将信赖域半径的连续减小次数置零，增大次数加一。

    elseif ratio > eta2 && (stop_tCG == 1 || stop_tCG == 2)
        Delta = min(gamma2*Delta, Delta_bar);
        consecutive_TRminus = 0;
        consecutive_TRplus = consecutive_TRplus + 1;

考虑当信赖域半径连续 5 次增大时，认为当前的信赖域半径过小，输出相应的提示信息。

        if consecutive_TRplus >= 5 && opts.verbose >= 1
            consecutive_TRplus = -inf;
            fprintf(' +++ Detected many consecutive TR+ (radius increases).\n');
            fprintf(' +++ Consider increasing options.Delta_bar by an order of magnitude.\n');
            fprintf(' +++ Current values: options.Delta_bar = %g and options.Delta0 = %g.\n', Delta_bar, Delta0);
        end

除了以上两种情况，不需要对信赖域半径进行调整，将其连续增大和减小次数都置零。

    else
        consecutive_TRplus = 0;
        consecutive_TRminus = 0;
    end

end

当从外层迭代退出时，记录退出信息。

out.iter = iter;
out.f    = f;
out.nrmg = nrmg;

end

参考页面

我们在页面实例：信赖域算法解逻辑回归问题中展示了该算法的一个应用。另外，该算法调用截断共轭梯度法 tCG 求解信赖域子问题，关于截断共轭梯度法参考截断共轭梯度法。

此页面的源代码请见： fminTR.m。

版权声明

此页面为《最优化：建模、算法与理论》、《最优化计算方法》配套代码。代码作者：文再文、刘浩洋、户将，代码整理与页面制作：杨昊桐。