Description of csminit1

0001 function [fhat,xhat,fcount,retcode] = csminit1(fcn,x0,f0,g0,badg,H0,varargin)
0002 % [fhat,xhat,fcount,retcode] = csminit1(fcn,x0,f0,g0,badg,H0,...
0003 %                                       P1,P2,P3,P4,P5,P6,P7,P8)
0004 % retcodes: 0, normal step.  5, largest step still improves too fast.
0005 % 4,2 back and forth adjustment of stepsize didn't finish.  3, smallest
0006 % stepsize still improves too slow.  6, no improvement found.  1, zero
0007 % gradient.
0008 %---------------------
0009 % Modified 7/22/96 to omit variable-length P list, for efficiency and compilation.
0010 % Places where the number of P's need to be altered or the code could be returned to
0011 % its old form are marked with ARGLIST comments.
0012 %
0013 % Fixed 7/17/93 to use inverse-hessian instead of hessian itself in bfgs
0014 % update.
0015 %
0016 % Fixed 7/19/93 to flip eigenvalues of H to get better performance when
0017 % it's not psd.
0018 
0019 % Original file downloaded from:
0020 % http://sims.princeton.edu/yftp/optimize/mfiles/csminit.m
0021 
0022 % Copyright (C) 1993-2007 Christopher Sims
0023 % Copyright (C) 2008-2011 Dynare Team
0024 %
0025 % This file is part of Dynare.
0026 %
0027 % Dynare is free software: you can redistribute it and/or modify
0028 % it under the terms of the GNU General Public License as published by
0029 % the Free Software Foundation, either version 3 of the License, or
0030 % (at your option) any later version.
0031 %
0032 % Dynare is distributed in the hope that it will be useful,
0033 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0034 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0035 % GNU General Public License for more details.
0036 %
0037 % You should have received a copy of the GNU General Public License
0038 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0039 
0040 %tailstr = ')';
0041 %for i=nargin-6:-1:1
0042 %   tailstr=[ ',P' num2str(i)  tailstr];
0043 %end
0044 %ANGLE = .03;
0045 ANGLE = .005;
0046 %THETA = .03;
0047 THETA = .3; %(0<THETA<.5) THETA near .5 makes long line searches, possibly fewer iterations.
0048 FCHANGE = 1000;
0049 MINLAMB = 1e-9;
0050 % fixed 7/15/94
0051 % MINDX = .0001;
0052 % MINDX = 1e-6;
0053 MINDFAC = .01;
0054 fcount=0;
0055 lambda=1;
0056 xhat=x0;
0057 f=f0;
0058 fhat=f0;
0059 g = g0;
0060 gnorm = norm(g);
0061 %
0062 if (gnorm < 1.e-12) && ~badg % put ~badg 8/4/94
0063     retcode =1;
0064     dxnorm=0;
0065     % gradient convergence
0066 else
0067     % with badg true, we don't try to match rate of improvement to directional
0068     % derivative.  We're satisfied just to get some improvement in f.
0069     %
0070     %if(badg)
0071     %   dx = -g*FCHANGE/(gnorm*gnorm);
0072     %  dxnorm = norm(dx);
0073     %  if dxnorm > 1e12
0074     %     disp('Bad, small gradient problem.')
0075     %     dx = dx*FCHANGE/dxnorm;
0076     %   end
0077     %else
0078     % Gauss-Newton step;
0079     %---------- Start of 7/19/93 mod ---------------
0080     %[v d] = eig(H0);
0081     %toc
0082     %d=max(1e-10,abs(diag(d)));
0083     %d=abs(diag(d));
0084     %dx = -(v.*(ones(size(v,1),1)*d'))*(v'*g);
0085     %      toc
0086     dx = -H0*g;
0087     %      toc
0088     dxnorm = norm(dx);
0089     if dxnorm > 1e12
0090         disp('Near-singular H problem.')
0091         dx = dx*FCHANGE/dxnorm;
0092     end
0093     dfhat = dx'*g0;
0094     %end
0095     %
0096     %
0097     if ~badg
0098         % test for alignment of dx with gradient and fix if necessary
0099         a = -dfhat/(gnorm*dxnorm);
0100         if a<ANGLE
0101             dx = dx - (ANGLE*dxnorm/gnorm+dfhat/(gnorm*gnorm))*g;
0102             dfhat = dx'*g;
0103             dxnorm = norm(dx);
0104             disp(sprintf('Correct for low angle: %g',a))
0105         end
0106     end
0107     disp(sprintf('Predicted improvement: %18.9f',-dfhat/2))
0108     %
0109     % Have OK dx, now adjust length of step (lambda) until min and
0110     % max improvement rate criteria are met.
0111     done=0;
0112     factor=3;
0113     shrink=1;
0114     lambdaMin=0;
0115     lambdaMax=inf;
0116     lambdaPeak=0;
0117     fPeak=f0;
0118     lambdahat=0;
0119     while ~done
0120         if size(x0,2)>1
0121             dxtest=x0+dx'*lambda;
0122         else
0123             dxtest=x0+dx*lambda;
0124         end
0125         % home
0126         f = feval(fcn,dxtest,varargin{:});
0127         %ARGLIST
0128         %f = feval(fcn,dxtest,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13);
0129         % f = feval(fcn,x0+dx*lambda,P1,P2,P3,P4,P5,P6,P7,P8);
0130         disp(sprintf('lambda = %10.5g; f = %20.7f',lambda,f ))
0131         %debug
0132         %disp(sprintf('Improvement too great? f0-f: %g, criterion: %g',f0-f,-(1-THETA)*dfhat*lambda))
0133         if f<fhat
0134             fhat=f;
0135             xhat=dxtest;
0136             lambdahat = lambda;
0137         end
0138         fcount=fcount+1;
0139         shrinkSignal = (~badg & (f0-f < max([-THETA*dfhat*lambda 0]))) | (badg & (f0-f) < 0) ;
0140         growSignal = ~badg & ( (lambda > 0)  &  (f0-f > -(1-THETA)*dfhat*lambda) );
0141         if  shrinkSignal  &&   ( (lambda>lambdaPeak) || (lambda<0) )
0142             if (lambda>0) && ((~shrink) || (lambda/factor <= lambdaPeak))
0143                 shrink=1;
0144                 factor=factor^.6;
0145                 while lambda/factor <= lambdaPeak
0146                     factor=factor^.6;
0147                 end
0148                 %if (abs(lambda)*(factor-1)*dxnorm < MINDX) || (abs(lambda)*(factor-1) < MINLAMB)
0149                 if abs(factor-1)<MINDFAC
0150                     if abs(lambda)<4
0151                         retcode=2;
0152                     else
0153                         retcode=7;
0154                     end
0155                     done=1;
0156                 end
0157             end
0158             if (lambda<lambdaMax) && (lambda>lambdaPeak)
0159                 lambdaMax=lambda;
0160             end
0161             lambda=lambda/factor;
0162             if abs(lambda) < MINLAMB
0163                 if (lambda > 0) && (f0 <= fhat)
0164                     % try going against gradient, which may be inaccurate
0165                     lambda = -lambda*factor^6
0166                 else
0167                     if lambda < 0
0168                         retcode = 6;
0169                     else
0170                         retcode = 3;
0171                     end
0172                     done = 1;
0173                 end
0174             end
0175         elseif  (growSignal && lambda>0) ||  (shrinkSignal && ((lambda <= lambdaPeak) && (lambda>0)))
0176             if shrink
0177                 shrink=0;
0178                 factor = factor^.6;
0179                 %if ( abs(lambda)*(factor-1)*dxnorm< MINDX ) || ( abs(lambda)*(factor-1)< MINLAMB)
0180                 if abs(factor-1)<MINDFAC
0181                     if abs(lambda)<4
0182                         retcode=4;
0183                     else
0184                         retcode=7;
0185                     end
0186                     done=1;
0187                 end
0188             end
0189             if ( f<fPeak ) && (lambda>0)
0190                 fPeak=f;
0191                 lambdaPeak=lambda;
0192                 if lambdaMax<=lambdaPeak
0193                     lambdaMax=lambdaPeak*factor*factor;
0194                 end
0195             end
0196             lambda=lambda*factor;
0197             if abs(lambda) > 1e20;
0198                 retcode = 5;
0199                 done =1;
0200             end
0201         else
0202             done=1;
0203             if factor < 1.2
0204                 retcode=7;
0205             else
0206                 retcode=0;
0207             end
0208         end
0209     end
0210 end
0211 disp(sprintf('Norm of dx %10.5g', dxnorm))
csminit1

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
