stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Trust region solver: recompute Jacobian only when necessary 

Previously, the solver would recompute the Jacobian at every iteration. But, if
an iteration fails, the current point is not moved (only the radius of the
trust region changes), hence it is not necessary to recompute the Jacobian in
that case. This commit implements this optimization.
time-shift
Sébastien Villemot 2019-11-14 16:03:50 +01:00
parent 49a17e75df
commit caf0c8e1f8
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
1 changed files with 18 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

40
matlab/trust_region.m
View File

@ -25,7 +25,7 @@ function [x,check,info] = trust_region(fcn,x0,j1,j2,jacobian_flag,gstep,tolf,tol
% none
% Copyright (C) 2008-2012 VZLU Prague, a.s.
% Copyright (C) 2014-2017 Dynare Team
% Copyright (C) 2014-2019 Dynare Team
%
% This file is part of Dynare.
%
@ -65,24 +65,27 @@ info = 0;
fvec = fcn (x, varargin{:});
fvec = fvec(j1);
fn = norm (fvec);
recompute_jacobian = true;
% 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(n,1))*eps^(1/3);
% 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) ;
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
end
recompute_jacobian = false;
end
% Get column norms, use them as scaling factors.
@ -164,20 +167,13 @@ while (niter < maxiter && ~info)
xn = norm (dg .* x(j2));
fvec = fvec1;
fn = fn1;
recompute_jacobian = true;
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.
% Tests for termination condition
if (fn <= tolf)
info = 1;
end