diff --git a/matlab/optimization/csminwel1.m b/matlab/optimization/csminwel1.m index b6e38fedc..dcc579f50 100644 --- a/matlab/optimization/csminwel1.m +++ b/matlab/optimization/csminwel1.m @@ -1,29 +1,48 @@ 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) -% 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. -% epsilon: scalar double, numerical differentiation increment -% varargin: A list of optional length of additional parameters that get handed off to fcn each -% time it is called. +%[fhat,xhat,ghat,Hhat,itct,fcount,retcodeh] = csminwel1(fcn,x0,H0,grad,crit,nit,method,epsilon,varargin) +% Inputs: +% fcn: [string] string naming the objective function to be minimized +% x0: [npar by 1] initial value of the parameter vector +% H0: [npar by npar] initial value for the inverse Hessian. Must be positive definite. +% grad: [string or empty matrix] 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: [scalar] Convergence criterion. Iteration will cease when it proves impossible to improve the +% function value by more than crit. +% nit: [scalar] Maximum number of iterations. +% method: [scalar] integer scalar for selecting gradient method: 2, 3 or 5 points formula. +% epsilon: [scalar] scalar double, numerical differentiation increment +% varargin: Optional additional inputs 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.) - +% writes 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.) +% +% Outputs: +% fh: [scalar] function value at minimum +% xh: [npar by 1] parameter vector at minimum +% gh [npar by 1] gradient vector +% H [npar by npar] inverse of the Hessian matrix +% itct [scalar] iteration count upon termination +% fcount [scalar] function iteration count upon termination +% retcodeh [scalar] return code: +% 0: normal step +% 1: zero gradient +% 2: back and forth on step length never finished +% 3: smallest step still improving too slow +% 4: back and forth on step length never finished +% 5: largest step still improving too fast +% 6: smallest step still improving too slow, reversed gradient +% 7: warning: possible inaccuracy in H matrix +% % Original file downloaded from: % http://sims.princeton.edu/yftp/optimize/mfiles/csminwel.m - +% % Copyright (C) 1993-2007 Christopher Sims -% Copyright (C) 2006-2012 Dynare Team +% Copyright (C) 2006-2015 Dynare Team % % This file is part of Dynare. % @@ -51,7 +70,6 @@ NumGrad= isempty(grad); done=0; itct=0; fcount=0; -snit=100; gh = []; H = []; retcodeh = []; @@ -74,20 +92,7 @@ if ~cost_flag end if NumGrad - switch method - case 2 - [g,badg] = numgrad2(fcn, f0, x0, epsilon, varargin{:}); - case 3 - [g,badg] = numgrad3(fcn, f0, x0, epsilon, varargin{:}); - case 5 - [g,badg] = numgrad5(fcn, f0, x0, epsilon, varargin{:}); - case 13 - [g,badg] = numgrad3_(fcn, f0, x0, epsilon, varargin{:}); - case 15 - [g,badg] = numgrad5_(fcn, f0, x0, epsilon, varargin{:}); - otherwise - error('csminwel1: Unknown method for gradient evaluation!') - end + [g, badg]=get_num_grad(method,fcn,f0,x0,epsilon,varargin{:}); elseif ischar(grad) [g,badg] = feval(grad,x0,varargin{:}); else @@ -105,19 +110,15 @@ while ~done 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)) + if Verbose + disp('-----------------') + disp(sprintf('f at the beginning of new iteration, %20.10f',f)) + end %-----------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; - [f1 x1 fc retcode1] = csminit1(fcn,x,f,g,badg,H,varargin{:}); - %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; + [f1, x1, fc, retcode1] = csminit1(fcn,x,f,g,badg,H,varargin{:}); fcount = fcount+fc; % erased on 8/4/94 % if (retcode == 1) || (abs(f1-f) < crit) @@ -132,22 +133,9 @@ while ~done wall1=1; badg1=1; else if NumGrad - switch method - case 2 - [g1 badg1] = numgrad2(fcn, f1, x1, epsilon, varargin{:}); - case 3 - [g1 badg1] = numgrad3(fcn, f1, x1, epsilon, varargin{:}); - case 5 - [g1,badg1] = numgrad5(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{:}); - otherwise - error('csminwel1: Unknown method for gradient evaluation!') - end + [g1, badg1]=get_num_grad(method,fcn,f1,x1,epsilon,varargin{:}); elseif ischar(grad), - [g1 badg1] = feval(grad,x1,varargin{:}); + [g1, badg1] = feval(grad,x1,varargin{:}); else [junk1,g1,junk2, cost_flag] = feval(fcn,x1,varargin{:}); badg1 = ~cost_flag; @@ -155,8 +143,6 @@ while ~done 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 if wall1 % && (~done) by Jinill % Bad gradient or back and forth on step length. Possibly at @@ -164,33 +150,19 @@ while ~done % %fcliff=fh;xcliff=xh; Hcliff=H+diag(diag(H).*rand(nx,1)); - disp('Cliff. Perturbing search direction.') - [f2 x2 fc retcode2] = csminit1(fcn,x,f,g,badg,Hcliff,varargin{:}); - %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); + if Verbose + disp('Cliff. Perturbing search direction.') + end + [f2, x2, fc, retcode2] = csminit1(fcn,x,f,g,badg,Hcliff,varargin{:}); fcount = fcount+fc; % put by Jinill if f2 < f if retcode2==2 || retcode2==4 wall2=1; badg2=1; else if NumGrad - switch method - case 2 - [g2 badg2] = numgrad2(fcn, f2, x2, epsilon, varargin{:}); - case 3 - [g2 badg2] = numgrad3(fcn, f2, x2, epsilon, varargin{:}); - case 5 - [g2,badg2] = numgrad5(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{:}); - otherwise - error('csminwel1: Unknown method for gradient evaluation!') - end + [g2, badg2]=get_num_grad(method,fcn,f2,x2,epsilon,varargin{:}); elseif ischar(grad), - [g2 badg2] = feval(grad,x2,varargin{:}); + [g2, badg2] = feval(grad,x2,varargin{:}); else [junk1,g2,junk2, cost_flag] = feval(fcn,x1,varargin{:}); badg2 = ~cost_flag; @@ -199,8 +171,6 @@ while ~done % 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') @@ -208,33 +178,19 @@ while ~done 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 - [f3 x3 fc retcode3] = csminit1(fcn,x,f,gcliff,0,eye(nx),varargin{:}); - %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); + if(size(x0,2)>1) + gcliff=gcliff'; + end + [f3, x3, fc, retcode3] = csminit1(fcn,x,f,gcliff,0,eye(nx),varargin{:}); fcount = fcount+fc; % put by Jinill if retcode3==2 || retcode3==4 - wall3=1; badg3=1; + wall3=1; + badg3=1; else if NumGrad - switch method - case 2 - [g3 badg3] = numgrad2(fcn, f3, x3, epsilon, varargin{:}); - case 3 - [g3 badg3] = numgrad3(fcn, f3, x3, epsilon, varargin{:}); - case 5 - [g3,badg3] = numgrad5(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{:}); - otherwise - error('csminwel1: Unknown method for gradient evaluation!') - end + [g3, badg3]=get_num_grad(method,fcn,f3,x3,epsilon,varargin{:}); elseif ischar(grad), - [g3 badg3] = feval(grad,x3,varargin{:}); + [g3, badg3] = feval(grad,x3,varargin{:}); else [junk1,g3,junk2, cost_flag] = feval(fcn,x1,varargin{:}); badg3 = ~cost_flag; @@ -242,8 +198,6 @@ while ~done 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 @@ -292,22 +246,9 @@ while ~done end if nogh if NumGrad - switch method - case 2 - [gh,badgh] = numgrad2(fcn, fh, xh, epsilon, varargin{:}); - case 3 - [gh,badgh] = numgrad3(fcn, fh, xh, epsilon, varargin{:}); - case 5 - [gh,badgh] = numgrad5(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{:}); - otherwise - error('csminwel1: Unknown method for gradient evaluation!') - end + [gh, badgh]=get_num_grad(method,fcn,fh,xh,epsilon,varargin{:}); elseif ischar(grad), - [gh badgh] = feval(grad, xh,varargin{:}); + [gh, badgh] = feval(grad, xh,varargin{:}); else [junk1,gh,junk2, cost_flag] = feval(fcn,x1,varargin{:}); badgh = ~cost_flag; @@ -316,11 +257,6 @@ while ~done badgh=1; end %end of picking - %ih - %fh - %xh - %gh - %badgh stuck = (abs(fh-f) < crit); if (~badg) && (~badgh) && (~stuck) H = bfgsi1(H,gh-g,xh-x); @@ -338,24 +274,47 @@ while ~done 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') - elseif (rc == 4) || (rc==2) - 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') + if Verbose || done + if rc ==0 + %do nothing, just a normal step + elseif 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') + elseif (rc == 4) || (rc==2) + 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') + else + error('Unaccounted Case, please contact the developers') + end end - % end + f=fh; x=xh; g=gh; badg=badgh; end -% what about making an m-file of 10 lines including numgrad.m -% since it appears three times in csminwel.m \ No newline at end of file + +end + +function [g, badg]=get_num_grad(method,fcn,f0,x0,epsilon,varargin) + switch method + case 2 + [g,badg] = numgrad2(fcn, f0, x0, epsilon, varargin{:}); + case 3 + [g,badg] = numgrad3(fcn, f0, x0, epsilon, varargin{:}); + case 5 + [g,badg] = numgrad5(fcn, f0, x0, epsilon, varargin{:}); + case 13 + [g,badg] = numgrad3_(fcn, f0, x0, epsilon, varargin{:}); + case 15 + [g,badg] = numgrad5_(fcn, f0, x0, epsilon, varargin{:}); + otherwise + error('csminwel1: Unknown method for gradient evaluation!') + end +end \ No newline at end of file