csminwel1.m 15.8 KB
Newer Older
1
2
function [fh,xh,gh,H,itct,fcount,retcodeh] = csminwel1(fcn,x0,H0,grad,crit,nit,method,epsilon,varargin)
%[fhat,xhat,ghat,Hhat,itct,fcount,retcodehat] = csminwel1(fcn,x0,H0,grad,crit,nit,method,epsilon,varargin)
3
4
5
6
7
8
9
10
11
12
% fcn:   string naming the objective function to be minimized
% x0:    initial value of the parameter vector
% H0:    initial value for the inverse Hessian.  Must be positive definite.
% grad:  Either a string naming a function that calculates the gradient, or the null matrix.
%        If it's null, the program calculates a numerical gradient.  In this case fcn must
%        be written so that it can take a matrix argument and produce a row vector of values.
% crit:  Convergence criterion.  Iteration will cease when it proves impossible to improve the
%        function value by more than crit.
% nit:   Maximum number of iterations.
% method: integer scalar, 2, 3 or 5 points formula.
13
% epsilon: scalar double, numerical differentiation increment
14
15
16
17
18
19
20
21
22
23
24
25
% varargin: A list of optional length of additional parameters that get handed off to fcn each
%        time it is called.
%        Note that if the program ends abnormally, it is possible to retrieve the current x,
%        f, and H from the files g1.mat and H.mat that are written at each iteration and at each
%        hessian update, respectively.  (When the routine hits certain kinds of difficulty, it
%        write g2.mat and g3.mat as well.  If all were written at about the same time, any of them
%        may be a decent starting point.  One can also start from the one with best function value.)

% Original file downloaded from:
% http://sims.princeton.edu/yftp/optimize/mfiles/csminwel.m

% Copyright (C) 1993-2007 Christopher Sims
Sébastien Villemot's avatar
Sébastien Villemot committed
26
% Copyright (C) 2006-2012 Dynare Team
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
%
% 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/>.

43
global objective_function_penalty_base
44

45
46
47
48
49
50
51
52
53
54
fh = [];
xh = [];
[nx,no]=size(x0);
nx=max(nx,no);
Verbose=1;
NumGrad= isempty(grad);
done=0;
itct=0;
fcount=0;
snit=100;
55
56
57
gh = [];
H = [];
retcodeh = [];
58
59
60
61
62
63
64
65
66
67
68
%tailstr = ')';
%stailstr = [];
% Lines below make the number of Pi's optional.  This is inefficient, though, and precludes
% use of the matlab compiler.  Without them, we use feval and the number of Pi's must be
% changed with the editor for each application.  Places where this is required are marked
% with ARGLIST comments
%for i=nargin-6:-1:1
%   tailstr=[ ',P' num2str(i)  tailstr];
%   stailstr=[' P' num2str(i) stailstr];
%end

69
[f0,junk1,junk2,cost_flag] = feval(fcn,x0,varargin{:});
70

71
if ~cost_flag
Stéphane Adjemian's avatar
Stéphane Adjemian committed
72
73
    disp('Bad initial parameter.')
    return
74
75
76
77
78
end

if NumGrad
    switch method
      case 2
79
        [g,badg] = numgrad2(fcn, f0, x0, epsilon, varargin{:});
80
      case 3
81
        [g,badg] = numgrad3(fcn, f0, x0, epsilon, varargin{:});
82
      case 5
Stéphane Adjemian's avatar
Stéphane Adjemian committed
83
84
85
86
87
88
89
90
91
        [g,badg] = numgrad5(fcn, f0, x0, epsilon, varargin{:});
      case {12, 22}
        [g,badg] = numgrad2_(fcn, f0, x0, epsilon, [], varargin{:});
      case {13, 23}
        [g,badg] = numgrad3_(fcn, f0, x0, epsilon, [], varargin{:});
      case {15, 25}
        [g,badg] = numgrad5_(fcn, f0, x0, epsilon, [], varargin{:});
      otherwise
        error('csminwel1: Unknown method for gradient evaluation!')
92
    end
93
elseif ischar(grad)
94
    [g,badg] = feval(grad,x0,varargin{:});
95
96
97
else
    g=junk1;
    badg=0;
98
99
100
101
102
103
104
end
retcode3=101;
x=x0;
f=f0;
H=H0;
cliff=0;
while ~done
105
106
    % penalty for dsge_likelihood and DsgeVarLikelihood
    objective_function_penalty_base = f;
Stéphane Adjemian's avatar
Stéphane Adjemian committed
107

108
109
110
111
112
113
114
115
116
117
    g1=[]; g2=[]; g3=[];
    %addition fj. 7/6/94 for control
    disp('-----------------')
    disp('-----------------')
    %disp('f and x at the beginning of new iteration')
    disp(sprintf('f at the beginning of new iteration, %20.10f',f))
    %-----------Comment out this line if the x vector is long----------------
    %   disp([sprintf('x = ') sprintf('%15.8g %15.8g %15.8g %15.8g\n',x)]);
    %-------------------------
    itct=itct+1;
118
    [f1 x1 fc retcode1] = csminit1(fcn,x,f,g,badg,H,varargin{:});
119
120
121
122
123
124
    %ARGLIST
    %[f1 x1 fc retcode1] = csminit(fcn,x,f,g,badg,H,P1,P2,P3,P4,P5,P6,P7,...
    %           P8,P9,P10,P11,P12,P13);
    % itct=itct+1;
    fcount = fcount+fc;
    % erased on 8/4/94
125
    % if (retcode == 1) || (abs(f1-f) < crit)
126
127
128
129
130
131
132
    %    done=1;
    % end
    % if itct > nit
    %    done = 1;
    %    retcode = -retcode;
    % end
    if retcode1 ~= 1
133
        if retcode1==2 || retcode1==4
134
135
136
            wall1=1; badg1=1;
        else
            if NumGrad
Stéphane Adjemian's avatar
Stéphane Adjemian committed
137
                switch method
138
                  case 2
139
                    [g1 badg1] = numgrad2(fcn, f1, x1, epsilon, varargin{:});
140
                  case 3
141
                    [g1 badg1] = numgrad3(fcn, f1, x1, epsilon, varargin{:});
142
                  case 5
Stéphane Adjemian's avatar
Stéphane Adjemian committed
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
                    [g1,badg1] = numgrad5(fcn, f1, x1, epsilon, varargin{:});
                  case 12
                    [g1,badg1] = numgrad2_(fcn, f1, x1, epsilon, [], varargin{:});
                  case 13
                    [g1,badg1] = numgrad3_(fcn, f1, x1, epsilon, [], varargin{:});
                  case 15
                    [g1,badg1] = numgrad5_(fcn, f1, x1, epsilon, [], varargin{:});
                  case 22
                    [g1,badg1] = numgrad2_(fcn, f1, x1, epsilon, abs(diag(H)), varargin{:});
                  case 23
                    [g1,badg1] = numgrad3_(fcn, f1, x1, epsilon, abs(diag(H)), varargin{:});
                  case 25
                    [g1,badg1] = numgrad5_(fcn, f1, x1, epsilon, abs(diag(H)), varargin{:});
                  otherwise
                    error('csminwel1: Unknown method for gradient evaluation!')
158
                end
159
            elseif ischar(grad),
160
                [g1 badg1] = feval(grad,x1,varargin{:});
161
162
163
            else
                [junk1,g1,junk2, cost_flag] = feval(fcn,x1,varargin{:});
                badg1 = ~cost_flag;
164
165
166
167
168
169
170
            end
            wall1=badg1;
            % g1
            save g1.mat g1 x1 f1 varargin;
            %ARGLIST
            %save g1 g1 x1 f1 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13;
        end
171
        if wall1 % && (~done) by Jinill
172
173
174
175
176
177
                 % Bad gradient or back and forth on step length.  Possibly at
                 % cliff edge.  Try perturbing search direction.
                 %
                 %fcliff=fh;xcliff=xh;
            Hcliff=H+diag(diag(H).*rand(nx,1));
            disp('Cliff.  Perturbing search direction.')
178
            [f2 x2 fc retcode2] = csminit1(fcn,x,f,g,badg,Hcliff,varargin{:});
179
180
181
182
183
            %ARGLIST
            %[f2 x2 fc retcode2] = csminit(fcn,x,f,g,badg,Hcliff,P1,P2,P3,P4,...
            %     P5,P6,P7,P8,P9,P10,P11,P12,P13);
            fcount = fcount+fc; % put by Jinill
            if  f2 < f
184
                if retcode2==2 || retcode2==4
185
186
187
188
189
                    wall2=1; badg2=1;
                else
                    if NumGrad
                        switch method
                          case 2
190
                            [g2 badg2] = numgrad2(fcn, f2, x2, epsilon, varargin{:});
191
                          case 3
192
                            [g2 badg2] = numgrad3(fcn, f2, x2, epsilon, varargin{:});
193
                          case 5
Stéphane Adjemian's avatar
Stéphane Adjemian committed
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
                            [g2,badg2] = numgrad5(fcn, f2, x2, epsilon, varargin{:});
                          case 12
                            [g2,badg2] = numgrad2_(fcn, f2, x2, epsilon, [], varargin{:});
                          case 13
                            [g2,badg2] = numgrad3_(fcn, f2, x2, epsilon, [], varargin{:});
                          case 15
                            [g2,badg2] = numgrad5_(fcn, f2, x2, epsilon, [], varargin{:});
                          case 22
                            [g2,badg2] = numgrad2_(fcn, f2, x2, epsilon, abs(diag(H)), varargin{:});
                          case 23
                            [g2,badg2] = numgrad3_(fcn, f2, x2, epsilon, abs(diag(H)), varargin{:});
                          case 25
                            [g2,badg2] = numgrad5_(fcn, f2, x2, epsilon, abs(diag(H)), varargin{:});
                          otherwise
                            error('csminwel1: Unknown method for gradient evaluation!')
209
                        end
210
                    elseif ischar(grad),
211
                        [g2 badg2] = feval(grad,x2,varargin{:});
212
213
214
                    else
                        [junk1,g2,junk2, cost_flag] = feval(fcn,x1,varargin{:});
                        badg2 = ~cost_flag;
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
                    end
                    wall2=badg2;
                    % g2
                    badg2
                    save g2.mat g2 x2 f2 varargin
                    %ARGLIST
                    %save g2 g2 x2 f2 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13;
                end
                if wall2
                    disp('Cliff again.  Try traversing')
                    if norm(x2-x1) < 1e-13
                        f3=f; x3=x; badg3=1;retcode3=101;
                    else
                        gcliff=((f2-f1)/((norm(x2-x1))^2))*(x2-x1);
                        if(size(x0,2)>1), gcliff=gcliff', end
230
                        [f3 x3 fc retcode3] = csminit1(fcn,x,f,gcliff,0,eye(nx),varargin{:});
231
232
233
234
235
                        %ARGLIST
                        %[f3 x3 fc retcode3] = csminit(fcn,x,f,gcliff,0,eye(nx),P1,P2,P3,...
                        %         P4,P5,P6,P7,P8,...
                        %      P9,P10,P11,P12,P13);
                        fcount = fcount+fc; % put by Jinill
236
                        if retcode3==2 || retcode3==4
237
238
239
240
241
                            wall3=1; badg3=1;
                        else
                            if NumGrad
                                switch method
                                  case 2
242
                                    [g3 badg3] = numgrad2(fcn, f3, x3, epsilon, varargin{:});
243
                                  case 3
244
                                    [g3 badg3] = numgrad3(fcn, f3, x3, epsilon, varargin{:});
245
                                  case 5
Stéphane Adjemian's avatar
Stéphane Adjemian committed
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
                                    [g3,badg3] = numgrad5(fcn, f3, x3, epsilon, varargin{:});
                                  case 12
                                    [g3,badg3] = numgrad2_(fcn, f3, x3, epsilon, [], varargin{:});
                                  case 13
                                    [g3,badg3] = numgrad3_(fcn, f3, x3, epsilon, [], varargin{:});
                                  case 15
                                    [g3,badg3] = numgrad5_(fcn, f3, x3, epsilon, [], varargin{:});
                                  case 22
                                    [g3,badg3] = numgrad2_(fcn, f3, x3, epsilon, abs(diag(H)), varargin{:});
                                  case 23
                                    [g3,badg3] = numgrad3_(fcn, f3, x3, epsilon, abs(diag(H)), varargin{:});
                                  case 25
                                    [g3,badg3] = numgrad5_(fcn, f3, x3, epsilon, abs(diag(H)), varargin{:});
                                  otherwise
                                    error('csminwel1: Unknown method for gradient evaluation!')
261
                                end
262
                            elseif ischar(grad),
263
                                [g3 badg3] = feval(grad,x3,varargin{:});
264
265
266
                            else
                                [junk1,g3,junk2, cost_flag] = feval(fcn,x1,varargin{:});
                                badg3 = ~cost_flag;
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
                            end
                            wall3=badg3;
                            % g3
                            save g3.mat g3 x3 f3 varargin;
                            %ARGLIST
                            %save g3 g3 x3 f3 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13;
                        end
                    end
                else
                    f3=f; x3=x; badg3=1; retcode3=101;
                end
            else
                f3=f; x3=x; badg3=1;retcode3=101;
            end
        else
            % normal iteration, no walls, or else we're finished here.
            f2=f; f3=f; badg2=1; badg3=1; retcode2=101; retcode3=101;
        end
Stéphane Adjemian's avatar
Stéphane Adjemian committed
285
    else
286
287
288
        f2=f;f3=f;f1=f;retcode2=retcode1;retcode3=retcode1;
    end
    %how to pick gh and xh
289
    if f3 < f - crit && badg3==0 && f3 < f2 && f3 < f1
290
291
        ih=3;
        fh=f3;xh=x3;gh=g3;badgh=badg3;retcodeh=retcode3;
292
    elseif f2 < f - crit && badg2==0 && f2 < f1
293
294
        ih=2;
        fh=f2;xh=x2;gh=g2;badgh=badg2;retcodeh=retcode2;
295
    elseif f1 < f - crit && badg1==0
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
        ih=1;
        fh=f1;xh=x1;gh=g1;badgh=badg1;retcodeh=retcode1;
    else
        [fh,ih] = min([f1,f2,f3]);
        %disp(sprintf('ih = %d',ih))
        %eval(['xh=x' num2str(ih) ';'])
        switch ih
          case 1
            xh=x1;
          case 2
            xh=x2;
          case 3
            xh=x3;
        end %case
            %eval(['gh=g' num2str(ih) ';'])
            %eval(['retcodeh=retcode' num2str(ih) ';'])
        retcodei=[retcode1,retcode2,retcode3];
        retcodeh=retcodei(ih);
        if exist('gh')
            nogh=isempty(gh);
        else
            nogh=1;
        end
        if nogh
            if NumGrad
                switch method
                  case 2
323
                    [gh,badgh] = numgrad2(fcn, fh, xh, epsilon, varargin{:});
324
                  case 3
325
                    [gh,badgh] = numgrad3(fcn, fh, xh, epsilon, varargin{:});
326
                  case 5
327
                    [gh,badgh] = numgrad5(fcn, fh, xh, epsilon, varargin{:});
Stéphane Adjemian's avatar
Stéphane Adjemian committed
328
329
330
331
332
333
334
335
336
337
338
339
340
341
                  case 12
                    [gh,badgh] = numgrad2_(fcn, fh, xh, epsilon, [], varargin{:});
                  case 13
                    [gh,badgh] = numgrad3_(fcn, fh, xh, epsilon, [], varargin{:});
                  case 15
                    [gh,badgh] = numgrad5_(fcn, fh, xh, epsilon, [], varargin{:});
                  case 22
                    [gh,badgh] = numgrad2_(fcn, fh, xh, epsilon, abs(diag(H)), varargin{:});
                  case 23
                    [gh,badgh] = numgrad3_(fcn, fh, xh, epsilon, abs(diag(H)), varargin{:});
                  case 25
                    [gh,badgh] = numgrad5_(fcn, fh, xh, epsilon, abs(diag(H)), varargin{:});
                  otherwise
                    error('csminwel1: Unknown method for gradient evaluation!')
342
                end
343
            elseif ischar(grad),
344
                [gh badgh] = feval(grad, xh,varargin{:});
345
346
347
            else
                [junk1,gh,junk2, cost_flag] = feval(fcn,x1,varargin{:});
                badgh = ~cost_flag;
348
349
350
351
352
353
354
355
356
357
358
            end
        end
        badgh=1;
    end
    %end of picking
    %ih
    %fh
    %xh
    %gh
    %badgh
    stuck = (abs(fh-f) < crit);
359
    if (~badg) && (~badgh) && (~stuck)
Marco Ratto's avatar
Marco Ratto committed
360
        H = bfgsi1(H,gh-g,xh-x);
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
    end
    if Verbose
        disp('----')
        disp(sprintf('Improvement on iteration %d = %18.9f',itct,f-fh))
    end
    % if Verbose
    if itct > nit
        disp('iteration count termination')
        done = 1;
    elseif stuck
        disp('improvement < crit termination')
        done = 1;
    end
    rc=retcodeh;
    if rc == 1
        disp('zero gradient')
    elseif rc == 6
        disp('smallest step still improving too slow, reversed gradient')
    elseif rc == 5
        disp('largest step still improving too fast')
381
    elseif (rc == 4) || (rc==2)
382
383
384
385
386
387
388
389
390
391
392
393
394
        disp('back and forth on step length never finished')
    elseif rc == 3
        disp('smallest step still improving too slow')
    elseif rc == 7
        disp('warning: possible inaccuracy in H matrix')
    end
    % end
    f=fh;
    x=xh;
    g=gh;
    badg=badgh;
end
% what about making an m-file of 10 lines including numgrad.m
395
% since it appears three times in csminwel.m