Add trust-region nonlinear solver.

Closes #260

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

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

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

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
+

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 \

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);