2021-07-22 16:46:10 +02:00
|
|
|
function [x,check,info] = trust_region(fcn,x0,j1,j2,jacobian_flag,gstep,tolf,tolx,maxiter,factor,debug,varargin)
|
2014-02-04 17:55:55 +01:00
|
|
|
% Solves systems of non linear equations of several variables, using a
|
|
|
|
% trust-region method.
|
|
|
|
%
|
|
|
|
% INPUTS
|
|
|
|
% fcn: name of the function to be solved
|
|
|
|
% x0: guess values
|
|
|
|
% j1: equations index for which the model is solved
|
|
|
|
% j2: unknown variables index
|
2019-03-19 14:26:16 +01:00
|
|
|
% jacobian_flag=true: jacobian given by the 'func' function
|
|
|
|
% jacobian_flag=false: jacobian obtained numerically
|
2014-02-04 17:55:55 +01:00
|
|
|
% gstep increment multiplier in numercial derivative
|
|
|
|
% computation
|
|
|
|
% tolf tolerance for residuals
|
|
|
|
% tolx tolerance for solution variation
|
|
|
|
% maxiter maximum number of iterations
|
2021-07-22 16:46:10 +02:00
|
|
|
% factor real scalar, determines the initial step bound
|
2014-02-04 17:55:55 +01:00
|
|
|
% debug debug flag
|
|
|
|
% varargin: list of arguments following bad_cond_flag
|
|
|
|
%
|
|
|
|
% OUTPUTS
|
|
|
|
% x: results
|
|
|
|
% check=1: the model can not be solved
|
2016-03-12 19:48:05 +01:00
|
|
|
% info: detailed exitcode
|
2014-02-04 17:55:55 +01:00
|
|
|
% SPECIAL REQUIREMENTS
|
|
|
|
% none
|
|
|
|
|
|
|
|
% Copyright (C) 2008-2012 VZLU Prague, a.s.
|
2021-01-15 17:17:48 +01:00
|
|
|
% Copyright (C) 2014-2021 Dynare Team
|
2014-02-04 17:55:55 +01:00
|
|
|
%
|
|
|
|
% This file is part of Dynare.
|
|
|
|
%
|
|
|
|
% Dynare is free software: you can redistribute it and/or modify
|
|
|
|
% it under the terms of the GNU General Public License as published by
|
|
|
|
% the Free Software Foundation, either version 3 of the License, or
|
|
|
|
% (at your option) any later version.
|
|
|
|
%
|
|
|
|
% Dynare is distributed in the hope that it will be useful,
|
|
|
|
% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
|
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
|
|
% GNU General Public License for more details.
|
|
|
|
%
|
|
|
|
% You should have received a copy of the GNU General Public License
|
2021-06-09 17:33:48 +02:00
|
|
|
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
2014-02-04 17:55:55 +01:00
|
|
|
%
|
|
|
|
% Initial author: Jaroslav Hajek <highegg@gmail.com>, for GNU Octave
|
|
|
|
|
|
|
|
if (ischar (fcn))
|
|
|
|
fcn = str2func (fcn);
|
|
|
|
end
|
|
|
|
|
|
|
|
n = length(j1);
|
|
|
|
|
|
|
|
% These defaults are rather stringent. I think that normally, user
|
|
|
|
% prefers accuracy to performance.
|
|
|
|
|
|
|
|
macheps = eps (class (x0));
|
|
|
|
|
|
|
|
niter = 1;
|
|
|
|
|
|
|
|
x = x0;
|
|
|
|
info = 0;
|
|
|
|
|
|
|
|
% Initial evaluation.
|
|
|
|
% Handle arbitrary shapes of x and f and remember them.
|
|
|
|
fvec = fcn (x, varargin{:});
|
|
|
|
fvec = fvec(j1);
|
|
|
|
fn = norm (fvec);
|
2019-11-14 16:03:50 +01:00
|
|
|
recompute_jacobian = true;
|
2014-02-04 17:55:55 +01:00
|
|
|
|
|
|
|
% Outer loop.
|
|
|
|
while (niter < maxiter && ~info)
|
|
|
|
|
2019-11-14 16:03:50 +01:00
|
|
|
% Calculate Jacobian (possibly via FD).
|
|
|
|
if recompute_jacobian
|
|
|
|
if jacobian_flag
|
|
|
|
[~, fjac] = fcn (x, varargin{:});
|
|
|
|
fjac = fjac(j1,j2);
|
|
|
|
else
|
|
|
|
dh = max(abs(x(j2)),gstep(1)*ones(n,1))*eps^(1/3);
|
|
|
|
|
|
|
|
for j = 1:n
|
|
|
|
xdh = x ;
|
|
|
|
xdh(j2(j)) = xdh(j2(j))+dh(j) ;
|
|
|
|
t = fcn(xdh,varargin{:});
|
|
|
|
fjac(:,j) = (t(j1) - fvec)./dh(j) ;
|
|
|
|
end
|
2014-02-04 17:55:55 +01:00
|
|
|
end
|
2019-11-14 16:03:50 +01:00
|
|
|
recompute_jacobian = false;
|
2014-02-04 17:55:55 +01:00
|
|
|
end
|
|
|
|
|
|
|
|
% Get column norms, use them as scaling factors.
|
2014-05-11 14:44:44 +02:00
|
|
|
jcn = sqrt(sum(fjac.*fjac))';
|
2014-02-04 17:55:55 +01:00
|
|
|
if (niter == 1)
|
|
|
|
dg = jcn;
|
|
|
|
dg(dg == 0) = 1;
|
|
|
|
else
|
|
|
|
% Rescale adaptively.
|
2020-02-20 09:34:13 +01:00
|
|
|
dg = max (dg, jcn);
|
2014-02-04 17:55:55 +01:00
|
|
|
end
|
|
|
|
|
|
|
|
if (niter == 1)
|
|
|
|
xn = norm (dg .* x(j2));
|
|
|
|
% FIXME: something better?
|
2021-07-22 16:46:10 +02:00
|
|
|
delta = max (xn, 1)*factor;
|
2014-02-04 17:55:55 +01:00
|
|
|
end
|
|
|
|
|
|
|
|
% Get trust-region model (dogleg) minimizer.
|
|
|
|
s = - dogleg (fjac, fvec, dg, delta);
|
|
|
|
w = fvec + fjac * s;
|
|
|
|
|
|
|
|
sn = norm (dg .* s);
|
|
|
|
if (niter == 1)
|
|
|
|
delta = min (delta, sn);
|
|
|
|
end
|
|
|
|
|
|
|
|
x2 = x;
|
|
|
|
x2(j2) = x2(j2) + s;
|
|
|
|
fvec1 = fcn (x2, varargin{:});
|
|
|
|
fvec1 = fvec1(j1);
|
|
|
|
fn1 = norm (fvec1);
|
|
|
|
|
|
|
|
if (fn1 < fn)
|
|
|
|
% Scaled actual reduction.
|
|
|
|
actred = 1 - (fn1/fn)^2;
|
|
|
|
else
|
|
|
|
actred = -1;
|
|
|
|
end
|
|
|
|
|
|
|
|
% Scaled predicted reduction, and ratio.
|
|
|
|
t = norm (w);
|
|
|
|
if (t < fn)
|
|
|
|
prered = 1 - (t/fn)^2;
|
|
|
|
ratio = actred / prered;
|
|
|
|
else
|
|
|
|
prered = 0;
|
|
|
|
ratio = 0;
|
|
|
|
end
|
|
|
|
|
|
|
|
% Update delta.
|
|
|
|
if (ratio < 0.1)
|
|
|
|
delta = 0.5*delta;
|
|
|
|
if (delta <= 1e1*macheps*xn)
|
|
|
|
% Trust region became uselessly small.
|
2015-10-06 11:09:41 +02:00
|
|
|
if (fn1 <= tolf)
|
|
|
|
info = 1;
|
|
|
|
else
|
|
|
|
info = -3;
|
|
|
|
end
|
2017-05-16 12:42:01 +02:00
|
|
|
break
|
2014-02-04 17:55:55 +01:00
|
|
|
end
|
|
|
|
elseif (abs (1-ratio) <= 0.1)
|
|
|
|
delta = 1.4142*sn;
|
|
|
|
elseif (ratio >= 0.5)
|
|
|
|
delta = max (delta, 1.4142*sn);
|
|
|
|
end
|
|
|
|
|
|
|
|
if (ratio >= 1e-4)
|
|
|
|
% Successful iteration.
|
|
|
|
x(j2) = x(j2) + s;
|
|
|
|
xn = norm (dg .* x(j2));
|
|
|
|
fvec = fvec1;
|
|
|
|
fn = fn1;
|
2019-11-14 16:03:50 +01:00
|
|
|
recompute_jacobian = true;
|
2014-02-04 17:55:55 +01:00
|
|
|
end
|
|
|
|
|
|
|
|
niter = niter + 1;
|
|
|
|
|
|
|
|
|
2019-11-14 16:03:50 +01:00
|
|
|
% Tests for termination condition
|
2014-07-01 18:43:16 +02:00
|
|
|
if (fn <= tolf)
|
2014-02-04 17:55:55 +01:00
|
|
|
info = 1;
|
|
|
|
end
|
|
|
|
end
|
2015-10-06 11:09:41 +02:00
|
|
|
if info==1
|
|
|
|
check = 0;
|
|
|
|
else
|
|
|
|
check = 1;
|
|
|
|
end
|
2014-02-04 17:55:55 +01:00
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
% Solve the double dogleg trust-region least-squares problem:
|
|
|
|
% Minimize norm(r*x-b) subject to the constraint norm(d.*x) <= delta,
|
|
|
|
% x being a convex combination of the gauss-newton and scaled gradient.
|
|
|
|
|
|
|
|
% TODO: error checks
|
|
|
|
% TODO: handle singularity, or leave it up to mldivide?
|
|
|
|
|
|
|
|
function x = dogleg (r, b, d, delta)
|
|
|
|
% Get Gauss-Newton direction.
|
2021-01-15 17:17:48 +01:00
|
|
|
if isoctave || matlab_ver_less_than('9.3')
|
|
|
|
% The decomposition() function does not exist in Octave and MATLAB < R2017b
|
|
|
|
x = r \ b;
|
|
|
|
else
|
|
|
|
x = decomposition(r, 'CheckCondition', false) \ b;
|
|
|
|
end
|
2014-02-04 17:55:55 +01:00
|
|
|
xn = norm (d .* x);
|
|
|
|
if (xn > delta)
|
|
|
|
% GN is too big, get scaled gradient.
|
|
|
|
s = (r' * b) ./ d;
|
|
|
|
sn = norm (s);
|
|
|
|
if (sn > 0)
|
|
|
|
% Normalize and rescale.
|
|
|
|
s = (s / sn) ./ d;
|
|
|
|
% Get the line minimizer in s direction.
|
|
|
|
tn = norm (r*s);
|
|
|
|
snm = (sn / tn) / tn;
|
|
|
|
if (snm < delta)
|
|
|
|
% Get the dogleg path minimizer.
|
|
|
|
bn = norm (b);
|
|
|
|
dxn = delta/xn; snmd = snm/delta;
|
|
|
|
t = (bn/sn) * (bn/xn) * snmd;
|
2014-05-07 13:34:19 +02:00
|
|
|
t = t - dxn * snmd^2 + sqrt ((t-dxn)^2 + (1-dxn^2)*(1-snmd^2));
|
2014-02-04 17:55:55 +01:00
|
|
|
alpha = dxn*(1-snmd^2) / t;
|
|
|
|
else
|
|
|
|
alpha = 0;
|
|
|
|
end
|
|
|
|
else
|
|
|
|
alpha = delta / xn;
|
|
|
|
snm = 0;
|
|
|
|
end
|
|
|
|
% Form the appropriate convex combination.
|
|
|
|
x = alpha * x + ((1-alpha) * min (snm, delta)) * s;
|
|
|
|
end
|
|
|
|
end
|
2020-03-11 12:47:23 +01:00
|
|
|
|