New implementation of the trust region algorithm.
Main difference is the presence of traps for NaN/Inf/Complex numbers in residuals or the Jacobian matrix. Also added new unit tests.trustregion
parent
1b2b3db380
commit
aa8439d4cc
GPG Key ID:
295C1FE89E17EB3C
|
|
@ -118,11 +118,6 @@ Copyright: 2011 Lawrence J. Christiano, Mathias Trabandt and Karl Walentin
|
|||
2013-2017 Dynare Team
|
||||
License: GPL-3+
|
||||
|
||||
Files: matlab/trust_region.m
|
||||
Copyright: 2008-2012 VZLU Prague, a.s.
|
||||
2014-2021 Dynare Team
|
||||
License: GPL-3+
|
||||
|
||||
Files: matlab/one_sided_hp_filter.m
|
||||
Copyright: 2010-2015 Alexander Meyer-Gohde
|
||||
2015-2017 Dynare Team
|
||||
|
|
|
|||
|
|
@ -14,6 +14,7 @@ function [x, errorflag, fvec, fjac, exitflag] = dynare_solve(f, x, options, vara
|
|||
% - errorflag [logical] scalar, true iff the model can not be solved.
|
||||
% - fvec [double] n×1 vector, function value at x (f(x), used for debugging when errorflag is true).
|
||||
% - fjac [double] n×n matrix, Jacobian value at x (J(x), used for debugging when errorflag is true).
|
||||
% - exitflag [integer] scalar,
|
||||
|
||||
% Copyright © 2001-2022 Dynare Team
|
||||
%
|
||||
|
|
@ -94,7 +95,7 @@ if jacobian_flag
|
|||
if ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:)))) ...
|
||||
|| any(~isreal(fvec)) || any(~isreal(fjac(:)))
|
||||
if max(abs(fvec)) < tolf %return if initial value solves problem
|
||||
info = 0;
|
||||
exitflag = -1;
|
||||
return;
|
||||
end
|
||||
disp_verbose('Randomize initial guess...',options.verbosity)
|
||||
|
|
@ -231,10 +232,7 @@ if options.solve_algo == 0
|
|||
elseif options.solve_algo==1
|
||||
[x, errorflag, exitflag] = solve1(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, tolf, tolx, maxit, [], options.debug, arguments{:});
|
||||
elseif options.solve_algo==9
|
||||
[x, errorflag] = trust_region(f,x, 1:nn, 1:nn, jacobian_flag, options.gstep, ...
|
||||
tolf, tolx, maxit, ...
|
||||
options.trust_region_initial_step_bound_factor, ...
|
||||
options.debug, arguments{:});
|
||||
[x, errorflag, exitflag] = trust_region(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, tolf, tolx, maxit, options.trust_region_initial_step_bound_factor, options.debug, arguments{:});
|
||||
elseif ismember(options.solve_algo, [2, 12, 4])
|
||||
if ismember(options.solve_algo, [2, 12])
|
||||
solver = @solve1;
|
||||
|
|
@ -310,11 +308,11 @@ elseif ismember(options.solve_algo, [2, 12, 4])
|
|||
dprintf('DYNARE_SOLVE (solve_algo=2|4|12): solving block %u with trust_region routine.', i);
|
||||
end
|
||||
end
|
||||
[x, errorflag] = solver(f, x, j1(j), j2(j), jacobian_flag, ...
|
||||
options.gstep, ...
|
||||
tolf, options.solve_tolx, maxit, ...
|
||||
options.trust_region_initial_step_bound_factor, ...
|
||||
options.debug, arguments{:});
|
||||
[x, errorflag, exitflag] = solver(f, x, j1(j), j2(j), jacobian_flag, ...
|
||||
options.gstep, ...
|
||||
tolf, options.solve_tolx, maxit, ...
|
||||
options.trust_region_initial_step_bound_factor, ...
|
||||
options.debug, arguments{:});
|
||||
fre = true;
|
||||
if errorflag
|
||||
return
|
||||
|
|
@ -323,10 +321,10 @@ elseif ismember(options.solve_algo, [2, 12, 4])
|
|||
fvec = feval(f, x, arguments{:});
|
||||
if max(abs(fvec))>tolf
|
||||
disp_verbose('Call solver on the full nonlinear problem.',options.verbosity)
|
||||
[x, errorflag] = solver(f, x, 1:nn, 1:nn, jacobian_flag, ...
|
||||
options.gstep, tolf, options.solve_tolx, maxit, ...
|
||||
options.trust_region_initial_step_bound_factor, ...
|
||||
options.debug, arguments{:});
|
||||
[x, errorflag, exitflag] = solver(f, x, 1:nn, 1:nn, jacobian_flag, ...
|
||||
options.gstep, tolf, options.solve_tolx, maxit, ...
|
||||
options.trust_region_initial_step_bound_factor, ...
|
||||
options.debug, arguments{:});
|
||||
end
|
||||
elseif options.solve_algo==3
|
||||
if jacobian_flag
|
||||
|
|
|
|||
|
|
@ -1,32 +1,42 @@
|
|||
function [x,check,info] = trust_region(fcn,x0,j1,j2,jacobian_flag,gstep,tolf,tolx,maxiter,factor,debug,varargin)
|
||||
function [x, errorflag, info] = trust_region(objfun, x, j1, j2, jacobianflag, gstep, tolf, tolx, maxiter, factor, debug, varargin)
|
||||
|
||||
% 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
|
||||
% jacobian_flag=true: jacobian given by the 'func' function
|
||||
% jacobian_flag=false: jacobian obtained numerically
|
||||
% gstep increment multiplier in numercial derivative
|
||||
% computation
|
||||
% tolf tolerance for residuals
|
||||
% tolx tolerance for solution variation
|
||||
% maxiter maximum number of iterations
|
||||
% factor real scalar, determines the initial step bound
|
||||
% debug debug flag
|
||||
% varargin: list of arguments following bad_cond_flag
|
||||
% - objfun [function handle, char] name of the routine evaluating the system of nonlinear equations (and possibly jacobian).
|
||||
% - x [double] n×1 vector, initial guess for the solution.
|
||||
% - j1 [integer] vector, equation indices defining a subproblem to be solved.
|
||||
% - j2 [integer] vector, unknown variable indices to be solved for in the subproblem.
|
||||
% - jacobianflag [logical] scalar, if true the jacobian matrix is expected to be returned as a second output argument when calling objfun, otherwise
|
||||
% the jacobian is computed numerically.
|
||||
% - gstep [double] scalar, increment multiplier in numerical derivative computation (only used if jacobianflag value is false).
|
||||
% - tolf [double] scalar, tolerance for residuals.
|
||||
% - tolx [double] scalar, tolerance for solution variation.
|
||||
% - maxiter [integer] scalar, maximum number of iterations;
|
||||
% - factor [double] scalar, determines the initial step bound.
|
||||
% - debug [logical] scalar, dummy argument.
|
||||
% - varargin: [cell] list of additional arguments to be passed to objfun.
|
||||
%
|
||||
% OUTPUTS
|
||||
% x: results
|
||||
% check=1: the model can not be solved
|
||||
% info: detailed exitcode
|
||||
% SPECIAL REQUIREMENTS
|
||||
% none
|
||||
% - x [double] n⨱1 vector, solution of the nonlinear system of equations.
|
||||
% - errorflag [logical] scalar, false iff nonlinear solver is successful.
|
||||
% - info [integer] scalar, information about the failure.
|
||||
%
|
||||
% REMARKS
|
||||
% [1] j1 and j2 muyst have the same number of elements.
|
||||
% [2] debug is here for compatibility purpose (see solve1), it does not affect the output.
|
||||
% [3] Possible values for info are:
|
||||
%
|
||||
% -1 if the initial guess is a solution of the nonlinear system of equations.
|
||||
% 0 if the nonlinear solver failed because the problem is ill behaved at the initial guess.
|
||||
% 1 if the nonlinear solver found a solution.
|
||||
% 2 if the maximum number of iterations has been reached.
|
||||
% 3 if spurious convergence (trust region radius is too small).
|
||||
% 4 if iteration is not making good progress, as measured by the improvement from the last 15 iterations.
|
||||
% 5 if no further improvement in the approximate solution x is possible (xtol is too small).
|
||||
|
||||
% Copyright (C) 2008-2012 VZLU Prague, a.s.
|
||||
% Copyright (C) 2014-2021 Dynare Team
|
||||
% Copyright © 2014-2022 Dynare Team
|
||||
%
|
||||
% This file is part of Dynare.
|
||||
%
|
||||
|
|
@ -42,186 +52,275 @@ function [x,check,info] = trust_region(fcn,x0,j1,j2,jacobian_flag,gstep,tolf,tol
|
|||
%
|
||||
% You should have received a copy of the GNU General Public License
|
||||
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
||||
%
|
||||
% Initial author: Jaroslav Hajek <highegg@gmail.com>, for GNU Octave
|
||||
|
||||
if (ischar (fcn))
|
||||
fcn = str2func (fcn);
|
||||
% Convert to function handle if necessary
|
||||
if ischar(objfun)
|
||||
objfun = str2func(objfun);
|
||||
end
|
||||
|
||||
%
|
||||
% Set constants
|
||||
%
|
||||
|
||||
radiusfactor = .5;
|
||||
actualreductionthreshold = .001;
|
||||
relativereductionthreshold = .1;
|
||||
|
||||
% number of equations
|
||||
n = length(j1);
|
||||
|
||||
% These defaults are rather stringent. I think that normally, user
|
||||
% prefers accuracy to performance.
|
||||
%
|
||||
% Initialization
|
||||
%
|
||||
|
||||
macheps = eps (class (x0));
|
||||
errorflag = true;
|
||||
|
||||
niter = 1;
|
||||
|
||||
x = x0;
|
||||
iter = 1;
|
||||
info = 0;
|
||||
delta = 0.0;
|
||||
|
||||
% Initial evaluation.
|
||||
% Handle arbitrary shapes of x and f and remember them.
|
||||
fvec = fcn (x, varargin{:});
|
||||
fvec = fvec(j1);
|
||||
fn = norm (fvec);
|
||||
recompute_jacobian = true;
|
||||
ncsucc = 0; ncslow = 0;
|
||||
|
||||
% Outer loop.
|
||||
while (niter < maxiter && ~info)
|
||||
%
|
||||
% Attempt to evaluate the residuals and jacobian matrix on the initial guess
|
||||
%
|
||||
|
||||
% 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
|
||||
try
|
||||
if jacobianflag
|
||||
[fval, fjac] = objfun(x, varargin{:});
|
||||
fval = fval(j1);
|
||||
fjac = fjac(j1,j2);
|
||||
else
|
||||
fval = objfun(x, varargin{:});
|
||||
fval = fval(j1);
|
||||
dh = max(abs(x(j2)), gstep(1)*ones(n,1))*eps^(1/3);
|
||||
fjac = zeros(n);
|
||||
for j = 1:n
|
||||
xdh = x;
|
||||
xdh(j2(j)) = xdh(j2(j))+dh(j);
|
||||
Fval = objfun(xdh, varargin{:});
|
||||
fjac(:,j) = (Fval(j1)-fval)./dh(j);
|
||||
end
|
||||
recompute_jacobian = false;
|
||||
end
|
||||
fnorm = norm(fval);
|
||||
catch
|
||||
% System of equation cannot be evaluated at the initial guess.
|
||||
return
|
||||
end
|
||||
|
||||
% Get column norms, use them as scaling factors.
|
||||
jcn = sqrt(sum(fjac.*fjac))';
|
||||
if (niter == 1)
|
||||
dg = jcn;
|
||||
dg(dg == 0) = 1;
|
||||
if any(isnan(fval)) || any(isinf(fval)) || any(~isreal(fval)) || any(isnan(fjac(:))) || any(isinf(fjac(:))) || any(~isreal(fjac(:)))
|
||||
% System of equations is ill-behaved at the initial guess.
|
||||
return
|
||||
end
|
||||
|
||||
if fnorm<tolf
|
||||
% Initial guess is a solution.
|
||||
errorflag = false;
|
||||
info = -1;
|
||||
return
|
||||
end
|
||||
|
||||
%
|
||||
% Main loop
|
||||
%
|
||||
|
||||
while iter<=maxiter && ~info
|
||||
% Compute the columns norm ofr the Jacobian matrix.
|
||||
fjacnorm = transpose(sqrt(sum((fjac.*fjac))));
|
||||
if iter==1
|
||||
% On the first iteration, calculate the norm of the scaled vector of unknowns x
|
||||
% and initialize the step bound delta. Scaling is done according to the norms of
|
||||
% the columns of the initial jacobian.
|
||||
fjacnorm__ = fjacnorm;
|
||||
fjacnorm__(fjacnorm<eps(1.0)) = 1.0;
|
||||
xnorm = norm(fjacnorm__.*x(j2));
|
||||
if xnorm>0
|
||||
delta = xnorm*factor;
|
||||
else
|
||||
delta = factor;
|
||||
end
|
||||
else
|
||||
% Rescale adaptively.
|
||||
dg = max (dg, jcn);
|
||||
fjacnorm__ = max(.1*fjacnorm__, fjacnorm);
|
||||
xnorm = norm(fjacnorm__.*x(j2));
|
||||
end
|
||||
|
||||
if (niter == 1)
|
||||
xn = norm (dg .* x(j2));
|
||||
% FIXME: something better?
|
||||
delta = max (xn, 1)*factor;
|
||||
% Determine the direction p (with trust region model defined in dogleg routine).
|
||||
p = dogleg(fjac, fval, fjacnorm__, delta);
|
||||
% Compute the norm of p.
|
||||
pnorm = norm(fjacnorm__.*p);
|
||||
xx = x;
|
||||
x0 = x;
|
||||
% Move along the direction p. Set a candidate value for x and predicted improvement for f.
|
||||
xx(j2) = x(j2) - p;
|
||||
ww = fval - fjac*p;
|
||||
% Evaluate the function at xx and calculate its norm.
|
||||
try
|
||||
fval1 = objfun(xx, varargin{:});
|
||||
fval1 = fval1(j1);
|
||||
catch
|
||||
% If evaluation of the residuals returns an error, then restart but with a smaller radius of the trust region.
|
||||
delta = delta*radiusfactor;
|
||||
iter = iter+1;
|
||||
continue
|
||||
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);
|
||||
if any(isnan(fval1)) || any(isinf(fval1)) || any(~isreal(fval1))
|
||||
% If evaluation of the residuals returns a NaN, an infinite number or a complex number, then restart but with a smaller radius of the trust region.
|
||||
delta = delta*radiusfactor;
|
||||
iter = iter+1;
|
||||
continue
|
||||
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;
|
||||
fnorm1 = norm(fval1);
|
||||
if fnorm1<tolf
|
||||
x = xx;
|
||||
errorflag = false;
|
||||
info = 1;
|
||||
continue
|
||||
end
|
||||
% Compute the scaled actual reduction.
|
||||
if fnorm1<fnorm
|
||||
actualreduction = 1.0-(fnorm1/fnorm)^2;
|
||||
else
|
||||
actred = -1;
|
||||
actualreduction = -1.0;
|
||||
end
|
||||
|
||||
% Scaled predicted reduction, and ratio.
|
||||
t = norm (w);
|
||||
if (t < fn)
|
||||
prered = 1 - (t/fn)^2;
|
||||
ratio = actred / prered;
|
||||
% Compute the scaled predicted reduction and the ratio of the actual to the
|
||||
% predicted reduction.
|
||||
tt = norm(ww);
|
||||
if tt<fnorm
|
||||
predictedreduction = 1.0 - (tt/fnorm)^2;
|
||||
else
|
||||
predictedreduction = 0.0;
|
||||
end
|
||||
if predictedreduction>0
|
||||
ratio = actualreduction/predictedreduction;
|
||||
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.
|
||||
if (fn1 <= tolf)
|
||||
info = 1;
|
||||
else
|
||||
info = -3;
|
||||
end
|
||||
break
|
||||
% Update the radius of the trust region if need be.
|
||||
if iter==1
|
||||
% On first iteration adjust the initial step bound.
|
||||
delta = min(delta, pnorm);
|
||||
end
|
||||
if ratio<relativereductionthreshold
|
||||
% Reduction is much smaller than predicted… Reduce the radius of the trust region.
|
||||
ncsucc = 0;
|
||||
delta = delta*radiusfactor;
|
||||
else
|
||||
ncsucc = ncsucc + 1;
|
||||
if abs(ratio-1.0)<relativereductionthreshold
|
||||
delta = pnorm/radiusfactor;
|
||||
elseif ratio>=radiusfactor || ncsucc>1
|
||||
delta = max(delta, pnorm/radiusfactor);
|
||||
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;
|
||||
recompute_jacobian = true;
|
||||
if ratio>1e-4
|
||||
% Successful iteration. Update x, xnorm, fval, fnorm and fjac.
|
||||
x = xx;
|
||||
fval = fval1;
|
||||
xnorm = norm(fjacnorm__.*x(j2));
|
||||
fnorm = fnorm1;
|
||||
end
|
||||
|
||||
niter = niter + 1;
|
||||
|
||||
|
||||
% Tests for termination condition
|
||||
if (fn <= tolf)
|
||||
info = 1;
|
||||
% Determine the progress of the iteration.
|
||||
ncslow = ncslow+1;
|
||||
if actualreduction>=actualreductionthreshold
|
||||
ncslow = 0;
|
||||
end
|
||||
iter = iter+1;
|
||||
if iter==maxiter
|
||||
info = 2;
|
||||
x(:) = inf;
|
||||
continue
|
||||
end
|
||||
if delta<tolx*xnorm
|
||||
info = 3;
|
||||
x(:) = inf;
|
||||
errorflag = true;
|
||||
continue
|
||||
end
|
||||
% Tests for termination and stringent tolerances.
|
||||
if max(.1*delta, pnorm)<=10*eps(xnorm)*xnorm
|
||||
% xtol is too small. no further improvement in
|
||||
% the approximate solution x is possible.
|
||||
info = 5;
|
||||
x(:) = inf;
|
||||
errorflag = true;
|
||||
continue
|
||||
end
|
||||
if ncslow==15
|
||||
info = 4;
|
||||
x(:) = inf;
|
||||
errorflag = true;
|
||||
continue
|
||||
end
|
||||
% Compute the jacobian for the next iteration.
|
||||
if jacobianflag
|
||||
try
|
||||
[~, fjac] = objfun(x, varargin{:});
|
||||
fjac = fjac(j1,j2);
|
||||
catch
|
||||
% If evaluation of the Jacobian matrix returns an error, then restart but with a smaller radius of the trust region.
|
||||
x = x0;
|
||||
delta = delta*radiusfactor;
|
||||
continue
|
||||
end
|
||||
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);
|
||||
Fval = objfun(xdh, varargin{:});
|
||||
fjac(:,j) = (Fval(j1)-fval)./dh(j);
|
||||
end
|
||||
end
|
||||
if any(isnan(fjac(:))) || any(isinf(fjac(:))) || any(~isreal(fjac(:)))
|
||||
% If evaluation of the Jacobian matrix returns NaNs, an infinite numbers or a complex numbers, then restart but with a smaller radius of the trust region.
|
||||
x = x0;
|
||||
delta = delta*radiusfactor;
|
||||
end
|
||||
end
|
||||
if info==1
|
||||
check = 0;
|
||||
else
|
||||
check = 1;
|
||||
end
|
||||
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.
|
||||
% Compute the Gauss-Newton direction.
|
||||
if isoctave || matlab_ver_less_than('9.3')
|
||||
% The decomposition() function does not exist in Octave and MATLAB < R2017b
|
||||
% The decomposition() function does not exist in Octave and MATLAB < R2017b
|
||||
x = r \ b;
|
||||
else
|
||||
x = decomposition(r, 'CheckCondition', false) \ b;
|
||||
end
|
||||
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;
|
||||
% Compute norm of scaled x
|
||||
qnorm = norm(d.*x);
|
||||
if qnorm<=delta
|
||||
% Gauss-Newton direction is acceptable. There is nothing to do here.
|
||||
else
|
||||
% Gauss-Newton direction is not acceptable…
|
||||
% Compute the scale gradient direction and its norm
|
||||
s = (r'*b)./d;
|
||||
gnorm = norm(s);
|
||||
if gnorm>0
|
||||
% Normalize and rescale → gradient direction.
|
||||
s = (s/gnorm)./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;
|
||||
t = t - dxn * snmd^2 + sqrt ((t-dxn)^2 + (1-dxn^2)*(1-snmd^2));
|
||||
alpha = dxn*(1-snmd^2) / t;
|
||||
temp0 = norm(r*s);
|
||||
sgnorm = gnorm/(temp0*temp0);
|
||||
if sgnorm<delta
|
||||
% The scaled gradient direction is not acceptable…
|
||||
% Compute the point along the dogleg at which the
|
||||
% quadratic is minimized.
|
||||
bnorm = norm(b);
|
||||
temp1 = delta/qnorm;
|
||||
temp2 = sgnorm/delta;
|
||||
temp0 = bnorm*bnorm*temp2/(gnorm*qnorm);
|
||||
temp0 = temp0 - temp1*temp2*temp2 + sqrt((temp0-temp1)^2+(1.0-temp1*temp1)*(1.0-temp2*temp2));
|
||||
alpha = temp1*(1.0-temp2*temp2)/temp0;
|
||||
else
|
||||
alpha = 0;
|
||||
% The scaled gradient direction is acceptable.
|
||||
alpha = 0.0;
|
||||
end
|
||||
else
|
||||
alpha = delta / xn;
|
||||
snm = 0;
|
||||
% If the norm of the scaled gradient direction is zero.
|
||||
alpha = delta/qnorm;
|
||||
sgnorm = 0.0;
|
||||
end
|
||||
% Form the appropriate convex combination.
|
||||
x = alpha * x + ((1-alpha) * min (snm, delta)) * s;
|
||||
end
|
||||
end
|
||||
|
||||
% Form the appropriate convex combination of the Gauss-Newton direction and the
|
||||
% scaled gradient direction.
|
||||
x = alpha*x + (1.0-alpha)*min(sgnorm, delta)*s;
|
||||
end
|
||||
|
|
@ -33,176 +33,120 @@ factor = 10;
|
|||
|
||||
auxstruct = struct();
|
||||
|
||||
% List of function handles
|
||||
objfun = { @rosenbrock,
|
||||
@powell1,
|
||||
@powell2,
|
||||
@wood,
|
||||
@helicalvalley,
|
||||
@watson,
|
||||
@chebyquad,
|
||||
@brown,
|
||||
@discreteboundaryvalue,
|
||||
@discreteintegralequation,
|
||||
@trigonometric,
|
||||
@variablydimensioned,
|
||||
@broydentridiagonal,
|
||||
@broydenbanded };
|
||||
|
||||
% FIXME block_trust_region (mex) or trust_region (matlab) do not work for all n (not sure we can fix that).
|
||||
% FIXME block_trust_region (mex) and trust_region (matlab) do not behave the same (spurious convergence for powell2 and trigonometric with block_trust_region).
|
||||
|
||||
%
|
||||
% Test mex routine
|
||||
%
|
||||
|
||||
t0 = clock;
|
||||
|
||||
try
|
||||
for i=1:length(objfun)
|
||||
NumberOfTests = NumberOfTests+1;
|
||||
x = rosenbrock();
|
||||
[x, errorflag] = block_trust_region(@rosenbrock, x, tolf, tolx, maxit, factor, false, auxstruct);
|
||||
if errorflag
|
||||
testFailed = testFailed+1;
|
||||
switch func2str(objfun{i})
|
||||
case 'helicalvalley'
|
||||
% FIXME block_trust_region is diverging if x(1)<0.
|
||||
x = helicalvalley();
|
||||
x(1) = 5;
|
||||
case 'chebyquad'
|
||||
% Fails with a system of 10 equations.
|
||||
x = objfun{i}(nan(9,1));
|
||||
case {'watson', 'brown', 'discreteintegralequation', 'discreteboundaryvalue', 'chebyquad', 'trigonometric', 'variablydimensioned', 'broydenbanded', 'broydentridiagonal'}
|
||||
x = objfun{i}(nan(4,1));
|
||||
otherwise
|
||||
x = objfun{i}();
|
||||
end
|
||||
try
|
||||
[x, errorflag] = block_trust_region(objfun{i}, x, tolf, tolx, maxit, factor, false, auxstruct);
|
||||
if isequal(func2str(objfun{i}), 'powell2')
|
||||
if ~errorflag
|
||||
testFailed = testFailed+1;
|
||||
if debug
|
||||
dprintf('Nonlinear solver is expected to fail on %s function but did not return an error.', func2str(objfun{i}))
|
||||
end
|
||||
end
|
||||
else
|
||||
if errorflag || norm(objfun{i}(x))>tolf
|
||||
testFailed = testFailed+1;
|
||||
if debug
|
||||
dprintf('Nonlinear solver (mex) failed on %s function (norm(f(x))=%s).', func2str(objfun{i}), num2str(norm(objfun{i}(x))))
|
||||
end
|
||||
end
|
||||
end
|
||||
catch
|
||||
testFailed = testFailed+1;
|
||||
if debug
|
||||
dprintf('Nonlinear solver (mex) failed on %s function.', func2str(objfun{i}))
|
||||
end
|
||||
end
|
||||
catch
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
|
||||
try
|
||||
t1 = clock; etime(t1, t0)
|
||||
|
||||
%
|
||||
% Test matlab routine
|
||||
%
|
||||
|
||||
for i=1:length(objfun)
|
||||
NumberOfTests = NumberOfTests+1;
|
||||
x = powell1();
|
||||
[x, errorflag] = block_trust_region(@powell1, x, tolf, tolx, maxit, factor, false, auxstruct);
|
||||
if errorflag
|
||||
testFailed = testFailed+1;
|
||||
switch func2str(objfun{i})
|
||||
case 'chebyquad'
|
||||
% Fails with a system of 10 equations.
|
||||
x = objfun{i}(nan(9,1));
|
||||
case {'watson', 'brown', 'discreteintegralequation', 'discreteboundaryvalue', 'trigonometric', 'variablydimensioned', 'broydenbanded', 'broydentridiagonal'}
|
||||
x = objfun{i}(nan(10,1));
|
||||
otherwise
|
||||
x = objfun{i}();
|
||||
end
|
||||
try
|
||||
[x, errorflag, info] = trust_region(objfun{i}, x, 1:length(x), 1:length(x), true, [], tolf, tolx, maxit, factor);
|
||||
if isequal(func2str(objfun{i}), 'powell2')
|
||||
if ~errorflag
|
||||
testFailed = testFailed+1;
|
||||
if debug
|
||||
dprintf('Nonlinear solver is expected to fail on %s function but did not return an error.', func2str(objfun{i}))
|
||||
end
|
||||
end
|
||||
if info~=3
|
||||
testFailed = testFailed+1;
|
||||
if debug
|
||||
dprintf('Nonlinear solver is expected to fail on %s function with info==3 but did not the correct value of info.', func2str(objfun{i}))
|
||||
end
|
||||
end
|
||||
else
|
||||
if errorflag
|
||||
testFailed = testFailed+1;
|
||||
if debug
|
||||
dprintf('Nonlinear solver failed on %s function (info=%s).', func2str(objfun{i}), int2str(info))
|
||||
end
|
||||
end
|
||||
end
|
||||
catch
|
||||
testFailed = testFailed+1;
|
||||
if debug
|
||||
dprintf('Nonlinear solver failed on %s function.', func2str(objfun{i}))
|
||||
end
|
||||
end
|
||||
catch
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
|
||||
try
|
||||
NumberOfTests = NumberOfTests+1;
|
||||
x = powell2();
|
||||
[x, errorflag] = block_trust_region(@powell2, x, tolf, tolx, maxit, factor, false, auxstruct);
|
||||
if errorflag
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
catch
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
|
||||
try
|
||||
NumberOfTests = NumberOfTests+1;
|
||||
x = wood();
|
||||
[x, errorflag] = block_trust_region(@wood, x, tolf, tolx, maxit, factor, false, auxstruct);
|
||||
if errorflag
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
catch
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
|
||||
try
|
||||
NumberOfTests = NumberOfTests+1;
|
||||
% FIXME block_trust_region is diverging if x(1)<0. Note that trust_region is not finding the
|
||||
% solution for the same initial conditions.
|
||||
x = helicalvalley();
|
||||
x(1) = 5;
|
||||
[x, errorflag] = block_trust_region(@helicalvalley, x, tolf, tolx, maxit, factor, false, auxstruct);
|
||||
if errorflag
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
catch
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
|
||||
try
|
||||
NumberOfTests = NumberOfTests+1;
|
||||
n = 10;
|
||||
x = watson(nan(n,1));
|
||||
[x, errorflag] = block_trust_region(@watson, x, tolf, tolx, maxit, factor, false, auxstruct);
|
||||
if errorflag
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
catch
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
|
||||
try
|
||||
NumberOfTests = NumberOfTests+1;
|
||||
% FIXME block_trust_region does not work for all n.
|
||||
n = 9;
|
||||
x = chebyquad(nan(n,1));
|
||||
[x, errorflag] = block_trust_region(@chebyquad, x, tolf, tolx, maxit, factor, false, auxstruct);
|
||||
if errorflag
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
catch
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
|
||||
try
|
||||
NumberOfTests = NumberOfTests+1;
|
||||
n = 10;
|
||||
x = brown(nan(n,1));
|
||||
[x, errorflag] = block_trust_region(@brown, x, tolf, tolx, maxit, factor, false, auxstruct);
|
||||
if errorflag
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
catch
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
|
||||
try
|
||||
NumberOfTests = NumberOfTests+1;
|
||||
n = 10;
|
||||
x = discreteboundaryvalue(nan(n,1));
|
||||
[x, errorflag] = block_trust_region(@discreteboundaryvalue, x, tolf, tolx, maxit, factor, false, auxstruct);
|
||||
if errorflag
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
catch
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
|
||||
try
|
||||
NumberOfTests = NumberOfTests+1;
|
||||
n = 10;
|
||||
x = discreteintegralequation(nan(n,1));
|
||||
[x, errorflag] = block_trust_region(@discreteintegralequation, x, tolf, tolx, maxit, factor, false, auxstruct);
|
||||
if errorflag
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
catch
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
|
||||
try
|
||||
NumberOfTests = NumberOfTests+1;
|
||||
n = 10;
|
||||
x = trigonometric(nan(n,1));
|
||||
[x, errorflag] = block_trust_region(@trigonometric, x, tolf, tolx, maxit, factor, false, auxstruct);
|
||||
if errorflag
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
catch
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
|
||||
try
|
||||
NumberOfTests = NumberOfTests+1;
|
||||
n = 10;
|
||||
x = variablydimensioned(nan(n,1));
|
||||
[x, errorflag] = block_trust_region(@variablydimensioned, x, tolf, tolx, maxit, factor, false, auxstruct);
|
||||
if errorflag
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
catch
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
|
||||
try
|
||||
NumberOfTests = NumberOfTests+1;
|
||||
n = 10;
|
||||
x = broydentridiagonal(nan(n,1));
|
||||
[x, errorflag] = block_trust_region(@broydentridiagonal, x, tolf, tolx, maxit, factor, false, auxstruct);
|
||||
if errorflag
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
catch
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
|
||||
try
|
||||
NumberOfTests = NumberOfTests+1;
|
||||
n = 10;
|
||||
x = broydenbanded(nan(n,1));
|
||||
[x, errorflag] = block_trust_region(@broydenbanded, x, tolf, tolx, maxit, factor, false, auxstruct);
|
||||
if errorflag
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
catch
|
||||
testFailed = testFailed+1;
|
||||
end
|
||||
|
||||
t1 = clock;
|
||||
t2 = clock; etime(t2, t1)
|
||||
|
||||
if ~debug
|
||||
cd(getenv('TOP_TEST_DIR'));
|
||||
|
|
@ -216,14 +160,14 @@ else
|
|||
fid = fopen('nonlinearsolvers.m.trs', 'w+');
|
||||
end
|
||||
if testFailed
|
||||
fprintf(fid,':test-result: FAIL\n');
|
||||
fprintf(fid,':test-result: FAIL\n');
|
||||
else
|
||||
fprintf(fid,':test-result: PASS\n');
|
||||
fprintf(fid,':test-result: PASS\n');
|
||||
end
|
||||
fprintf(fid,':number-tests: %i\n', NumberOfTests);
|
||||
fprintf(fid,':number-failed-tests: %i\n', testFailed);
|
||||
fprintf(fid,':list-of-passed-tests: nonlinearsolvers.m\n');
|
||||
fprintf(fid,':elapsed-time: %f\n', etime(t1, t0));
|
||||
fprintf(fid,':elapsed-time: %f\n', etime(t2, t0));
|
||||
fclose(fid);
|
||||
|
||||
if ~debug
|
||||
|
|
|
|||
Loading…
Reference in New Issue