diff --git a/doc/dynare.texi b/doc/dynare.texi
index e1c151facdd1b9f614a779fd9aee7407994eaf3e..0ceb22d9f0be51489fd744cd12044a678ba5e98c 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -2649,8 +2649,7 @@ using the same solver as value @code{1}
Use Chris Sims' solver
@item 4
-Same as value @code{2}, except that it does not try to adapt the
-search direction when the Jacobian is nearly singular
+Trust-region solver with autoscaling.
@item 5
Newton algorithm with a sparse Gaussian elimination (SPE) (requires
diff --git a/license.txt b/license.txt
index 8c5be8f0bf43b8d1e75f51f300cbe517c1c4d620..57fc483ebb3d3533d5e7d64d5b50050fae47b234 100644
--- a/license.txt
+++ b/license.txt
@@ -58,6 +58,11 @@ Copyright: 2011 Lawrence J. Christiano, Mathias Trabandt and Karl Walentin
2013 Dynare Team
License: GPL-3+
+Files: matlab/trust_region.m
+Copyright: 2008-2012 VZLU Prague, a.s.
+ 2014 Dynare Team
+License: GPL-3+
+
Files: matlab/simpsa.m matlab/simpsaget.m matlab/simpsaset.m
Copyright: 2005 Henning Schmidt, FCC, henning@fcc.chalmers.se
2006 Brecht Donckels, BIOMATH, brecht.donckels@ugent.be
diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index ba5a676f962376c702b49e6be5cc2b140505553d..b1ab7aefd9cb1232a448131c9caf7e53b99b5e13 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -124,6 +124,13 @@ elseif options_.solve_algo == 1
tolf,options_.solve_tolx, ...
options_.steady.maxit,options_.debug,varargin{:});
elseif options_.solve_algo == 2 || options_.solve_algo == 4
+
+ if options_.solve_algo == 2
+ solver = @solve1;
+ else
+ solver = @trust_region;
+ end
+
if ~jacobian_flag
fjac = zeros(nn,nn) ;
dh = max(abs(x),options_.gstep(1)*ones(nn,1))*eps^(1/3);
@@ -144,7 +151,7 @@ elseif options_.solve_algo == 2 || options_.solve_algo == 4
if options_.debug
disp(['DYNARE_SOLVE (solve_algo=2|4): solving block ' num2str(i) ', of size ' num2str(r(i+1)-r(i)) ]);
end
- [x,info]=solve1(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),jacobian_flag, ...
+ [x,info]=solver(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),jacobian_flag, ...
options_.gstep, ...
tolf,options_.solve_tolx, ...
options_.steady.maxit,options_.debug,varargin{:});
@@ -154,7 +161,7 @@ elseif options_.solve_algo == 2 || options_.solve_algo == 4
end
fvec = feval(func,x,varargin{:});
if max(abs(fvec)) > tolf
- [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag, ...
+ [x,info]=solver(func,x,1:nn,1:nn,jacobian_flag, ...
options_.gstep, tolf,options_.solve_tolx, ...
options_.steady.maxit,options_.debug,varargin{:});
end
diff --git a/matlab/trust_region.m b/matlab/trust_region.m
new file mode 100644
index 0000000000000000000000000000000000000000..10e8d098e6cd633d737216b63a354c9cdc96a819
--- /dev/null
+++ b/matlab/trust_region.m
@@ -0,0 +1,244 @@
+function [x,check] = trust_region(fcn,x0,j1,j2,jacobian_flag,gstep,tolf,tolx,maxiter,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=1: jacobian given by the 'func' function
+% jacobian_flag=0: jacobian obtained numerically
+% gstep increment multiplier in numercial derivative
+% computation
+% tolf tolerance for residuals
+% tolx tolerance for solution variation
+% maxiter maximum number of iterations
+% debug debug flag
+% varargin: list of arguments following bad_cond_flag
+%
+% OUTPUTS
+% x: results
+% check=1: the model can not be solved
+%
+% SPECIAL REQUIREMENTS
+% none
+
+% Copyright (C) 2008-2012 VZLU Prague, a.s.
+% Copyright (C) 2014 Dynare Team
+%
+% 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
+% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+%
+% 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);
+jcn = nan(n, 1);
+
+% Outer loop.
+while (niter < maxiter && ~info)
+
+ % Calculate function value and Jacobian (possibly via FD).
+ if jacobian_flag
+ [fvec, fjac] = fcn (x, varargin{:});
+ fvec = fvec(j1);
+ fjac = fjac(j1,j2);
+ else
+ dh = max(abs(x(j2)),gstep(1)*ones(nn,1))*eps^(1/3);
+
+ for j = 1:nn
+ xdh = x ;
+ xdh(j2(j)) = xdh(j2(j))+dh(j) ;
+ t = fcn(xdh,varargin{:});
+ fjac(:,j) = (t(j1) - fvec)./dh(j) ;
+ g(j) = fvec'*fjac(:,j) ;
+ end
+ end
+
+ % Get column norms, use them as scaling factors.
+ for j = 1:n
+ jcn(j) = norm(fjac(:,j));
+ end
+ if (niter == 1)
+ dg = jcn;
+ dg(dg == 0) = 1;
+ else
+ % Rescale adaptively.
+ % FIXME: the original minpack used the following rescaling strategy:
+ % dg = max (dg, jcn);
+ % but it seems not good if we start with a bad guess yielding Jacobian
+ % columns with large norms that later decrease, because the corresponding
+ % variable will still be overscaled. So instead, we only give the old
+ % scaling a small momentum, but do not honor it.
+
+ dg = max (0.1*dg, jcn);
+ end
+
+ if (niter == 1)
+ xn = norm (dg .* x(j2));
+ % FIXME: something better?
+ delta = max (xn, 1);
+ 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.
+ info = -3;
+ break;
+ 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;
+ end
+
+ niter = niter + 1;
+
+
+ % Tests for termination conditions. A mysterious place, anything
+ % can happen if you change something here...
+
+ % The rule of thumb (which I'm not sure M*b is quite following)
+ % is that for a tolerance that depends on scaling, only 0 makes
+ % sense as a default value. But 0 usually means uselessly long
+ % iterations, so we need scaling-independent tolerances wherever
+ % possible.
+
+ % FIXME -- why tolf*n*xn? If abs (e) ~ abs(x) * eps is a vector
+ % of perturbations of x, then norm (fjac*e) <= eps*n*xn, i.e. by
+ % tolf ~ eps we demand as much accuracy as we can expect.
+ if (fn <= tolf*n*xn)
+ info = 1;
+ % The following tests done only after successful step.
+ elseif (ratio >= 1e-4)
+ % This one is classic. Note that we use scaled variables again,
+ % but compare to scaled step, so nothing bad.
+ if (sn <= tolx*xn)
+ info = 2;
+ % Again a classic one. It seems weird to use the same tolf
+ % for two different tests, but that's what M*b manual appears
+ % to say.
+ elseif (actred < tolf)
+ info = 3;
+ end
+ end
+end
+
+check = ~info;
+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.
+x = r \ b;
+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;
+ t = t - dxn * snmd^2 - sqrt ((t-dxn)^2 + (1-dxn^2)*(1-snmd^2));
+ 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
+
diff --git a/tests/Makefile.am b/tests/Makefile.am
index a02dfb27e42b9c995eea00c4ab0dc152343501fa..625979fdcb679075c30d32c72d147dfd87ca7540 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -45,6 +45,7 @@ MODFILES = \
steady_state/walsh1_ssm.mod \
steady_state/walsh1_ssm_block.mod \
steady_state/multi_leads.mod \
+ steady_state/example1_trust_region.mod \
steady_state_operator/standard.mod \
steady_state_operator/use_dll.mod \
steady_state_operator/block.mod \
diff --git a/tests/steady_state/example1_trust_region.mod b/tests/steady_state/example1_trust_region.mod
new file mode 100644
index 0000000000000000000000000000000000000000..a6ed9395b70be898bf12eded85ed65735b39c439
--- /dev/null
+++ b/tests/steady_state/example1_trust_region.mod
@@ -0,0 +1,46 @@
+// Test trust region nonlinear solver (solve_algo=4)
+var y, c, k, a, h, b;
+varexo e, u;
+
+parameters beta, rho, alpha, delta, theta, psi, tau;
+
+alpha = 0.36;
+rho = 0.95;
+tau = 0.025;
+beta = 0.99;
+delta = 0.025;
+psi = 0;
+theta = 2.95;
+
+phi = 0.1;
+
+model;
+c*theta*h^(1+psi)=(1-alpha)*y;
+k = beta*(((exp(b)*c)/(exp(b(+1))*c(+1)))
+ *(exp(b(+1))*alpha*y(+1)+(1-delta)*k));
+y = exp(a)*(k(-1)^alpha)*(h^(1-alpha));
+k = exp(b)*(y-c)+(1-delta)*k(-1);
+a = rho*a(-1)+tau*b(-1) + e;
+b = tau*a(-1)+rho*b(-1) + u;
+end;
+
+initval;
+y = 1;
+c = 0.8;
+h = 0.3;
+k = 10;
+a = 0;
+b = 0;
+e = 0;
+u = 0;
+end;
+
+steady(solve_algo=4);
+
+shocks;
+var e; stderr 0.009;
+var u; stderr 0.009;
+var e, u = phi*0.009*0.009;
+end;
+
+stoch_simul(order=1,nomoments,irf=0);