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