1function [x,check,info] = trust_region(fcn,x0,j1,j2,jacobian_flag,gstep,tolf,tolx,maxiter,debug,varargin) 2% Solves systems of non linear equations of several variables, using a 3% trust-region method. 4% 5% INPUTS 6% fcn: name of the function to be solved 7% x0: guess values 8% j1: equations index for which the model is solved 9% j2: unknown variables index 10% jacobian_flag=true: jacobian given by the 'func' function 11% jacobian_flag=false: jacobian obtained numerically 12% gstep increment multiplier in numercial derivative 13% computation 14% tolf tolerance for residuals 15% tolx tolerance for solution variation 16% maxiter maximum number of iterations 17% debug debug flag 18% varargin: list of arguments following bad_cond_flag 19% 20% OUTPUTS 21% x: results 22% check=1: the model can not be solved 23% info: detailed exitcode 24% SPECIAL REQUIREMENTS 25% none 26 27% Copyright (C) 2008-2012 VZLU Prague, a.s. 28% Copyright (C) 2014-2019 Dynare Team 29% 30% This file is part of Dynare. 31% 32% Dynare is free software: you can redistribute it and/or modify 33% it under the terms of the GNU General Public License as published by 34% the Free Software Foundation, either version 3 of the License, or 35% (at your option) any later version. 36% 37% Dynare is distributed in the hope that it will be useful, 38% but WITHOUT ANY WARRANTY; without even the implied warranty of 39% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 40% GNU General Public License for more details. 41% 42% You should have received a copy of the GNU General Public License 43% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 44% 45% Initial author: Jaroslav Hajek <highegg@gmail.com>, for GNU Octave 46 47if (ischar (fcn)) 48 fcn = str2func (fcn); 49end 50 51n = length(j1); 52 53% These defaults are rather stringent. I think that normally, user 54% prefers accuracy to performance. 55 56macheps = eps (class (x0)); 57 58niter = 1; 59 60x = x0; 61info = 0; 62 63% Initial evaluation. 64% Handle arbitrary shapes of x and f and remember them. 65fvec = fcn (x, varargin{:}); 66fvec = fvec(j1); 67fn = norm (fvec); 68recompute_jacobian = true; 69 70% Outer loop. 71while (niter < maxiter && ~info) 72 73 % Calculate Jacobian (possibly via FD). 74 if recompute_jacobian 75 if jacobian_flag 76 [~, fjac] = fcn (x, varargin{:}); 77 fjac = fjac(j1,j2); 78 else 79 dh = max(abs(x(j2)),gstep(1)*ones(n,1))*eps^(1/3); 80 81 for j = 1:n 82 xdh = x ; 83 xdh(j2(j)) = xdh(j2(j))+dh(j) ; 84 t = fcn(xdh,varargin{:}); 85 fjac(:,j) = (t(j1) - fvec)./dh(j) ; 86 end 87 end 88 recompute_jacobian = false; 89 end 90 91 % Get column norms, use them as scaling factors. 92 jcn = sqrt(sum(fjac.*fjac))'; 93 if (niter == 1) 94 dg = jcn; 95 dg(dg == 0) = 1; 96 else 97 % Rescale adaptively. 98 % FIXME: the original minpack used the following rescaling strategy: 99 % dg = max (dg, jcn); 100 % but it seems not good if we start with a bad guess yielding Jacobian 101 % columns with large norms that later decrease, because the corresponding 102 % variable will still be overscaled. So instead, we only give the old 103 % scaling a small momentum, but do not honor it. 104 105 dg = max (0.1*dg, jcn); 106 end 107 108 if (niter == 1) 109 xn = norm (dg .* x(j2)); 110 % FIXME: something better? 111 delta = max (xn, 1); 112 end 113 114 % Get trust-region model (dogleg) minimizer. 115 s = - dogleg (fjac, fvec, dg, delta); 116 w = fvec + fjac * s; 117 118 sn = norm (dg .* s); 119 if (niter == 1) 120 delta = min (delta, sn); 121 end 122 123 x2 = x; 124 x2(j2) = x2(j2) + s; 125 fvec1 = fcn (x2, varargin{:}); 126 fvec1 = fvec1(j1); 127 fn1 = norm (fvec1); 128 129 if (fn1 < fn) 130 % Scaled actual reduction. 131 actred = 1 - (fn1/fn)^2; 132 else 133 actred = -1; 134 end 135 136 % Scaled predicted reduction, and ratio. 137 t = norm (w); 138 if (t < fn) 139 prered = 1 - (t/fn)^2; 140 ratio = actred / prered; 141 else 142 prered = 0; 143 ratio = 0; 144 end 145 146 % Update delta. 147 if (ratio < 0.1) 148 delta = 0.5*delta; 149 if (delta <= 1e1*macheps*xn) 150 % Trust region became uselessly small. 151 if (fn1 <= tolf) 152 info = 1; 153 else 154 info = -3; 155 end 156 break 157 end 158 elseif (abs (1-ratio) <= 0.1) 159 delta = 1.4142*sn; 160 elseif (ratio >= 0.5) 161 delta = max (delta, 1.4142*sn); 162 end 163 164 if (ratio >= 1e-4) 165 % Successful iteration. 166 x(j2) = x(j2) + s; 167 xn = norm (dg .* x(j2)); 168 fvec = fvec1; 169 fn = fn1; 170 recompute_jacobian = true; 171 end 172 173 niter = niter + 1; 174 175 176 % Tests for termination condition 177 if (fn <= tolf) 178 info = 1; 179 end 180end 181if info==1 182 check = 0; 183else 184 check = 1; 185end 186end 187 188 189% Solve the double dogleg trust-region least-squares problem: 190% Minimize norm(r*x-b) subject to the constraint norm(d.*x) <= delta, 191% x being a convex combination of the gauss-newton and scaled gradient. 192 193% TODO: error checks 194% TODO: handle singularity, or leave it up to mldivide? 195 196function x = dogleg (r, b, d, delta) 197% Get Gauss-Newton direction. 198x = r \ b; 199xn = norm (d .* x); 200if (xn > delta) 201 % GN is too big, get scaled gradient. 202 s = (r' * b) ./ d; 203 sn = norm (s); 204 if (sn > 0) 205 % Normalize and rescale. 206 s = (s / sn) ./ d; 207 % Get the line minimizer in s direction. 208 tn = norm (r*s); 209 snm = (sn / tn) / tn; 210 if (snm < delta) 211 % Get the dogleg path minimizer. 212 bn = norm (b); 213 dxn = delta/xn; snmd = snm/delta; 214 t = (bn/sn) * (bn/xn) * snmd; 215 t = t - dxn * snmd^2 + sqrt ((t-dxn)^2 + (1-dxn^2)*(1-snmd^2)); 216 alpha = dxn*(1-snmd^2) / t; 217 else 218 alpha = 0; 219 end 220 else 221 alpha = delta / xn; 222 snm = 0; 223 end 224 % Form the appropriate convex combination. 225 x = alpha * x + ((1-alpha) * min (snm, delta)) * s; 226end 227end 228