From aa8439d4ccbfcf3761c7e168a8423adc8cfce204 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Ry=C3=BBk=29?= Date: Fri, 11 Mar 2022 18:39:40 +0100 Subject: [PATCH] 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. --- license.txt | 5 - matlab/dynare_solve.m | 26 ++- matlab/trust_region.m | 427 ++++++++++++++++++++++++--------------- tests/nonlinearsolvers.m | 270 ++++++++++--------------- 4 files changed, 382 insertions(+), 346 deletions(-) diff --git a/license.txt b/license.txt index 61a2f6a5d..a479aeeb7 100644 --- a/license.txt +++ b/license.txt @@ -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 diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m index 024485500..f68d14cce 100644 --- a/matlab/dynare_solve.m +++ b/matlab/dynare_solve.m @@ -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 diff --git a/matlab/trust_region.m b/matlab/trust_region.m index ab5d6a908..8fcaaaa66 100644 --- a/matlab/trust_region.m +++ b/matlab/trust_region.m @@ -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 . -% -% Initial author: Jaroslav Hajek , 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 fnorm0 + 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 fnorm10 + 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=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 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 sgnormtolf + 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