From 0cdfe4d6c7d773059e6346bc75f80c5ae179267d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Fri, 20 Sep 2013 23:15:00 +0200 Subject: [PATCH] Added new routines for computing gradients (options 13 and 15). --- matlab/csminwel1.m | 20 ++++++++++++++ matlab/numgrad3_.m | 67 ++++++++++++++++++++++++++++++++++++++++++++++ matlab/numgrad5_.m | 64 +++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 151 insertions(+) create mode 100644 matlab/numgrad3_.m create mode 100644 matlab/numgrad5_.m diff --git a/matlab/csminwel1.m b/matlab/csminwel1.m index 6b0ea0166..403b9a58e 100644 --- a/matlab/csminwel1.m +++ b/matlab/csminwel1.m @@ -81,6 +81,10 @@ if NumGrad [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 @@ -135,6 +139,10 @@ while ~done [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 @@ -174,6 +182,10 @@ while ~done [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 @@ -214,6 +226,10 @@ while ~done [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 @@ -283,6 +299,10 @@ while ~done [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 diff --git a/matlab/numgrad3_.m b/matlab/numgrad3_.m new file mode 100644 index 000000000..a4b4e56a6 --- /dev/null +++ b/matlab/numgrad3_.m @@ -0,0 +1,67 @@ +function [g, badg, f0, f1, f2] = numgrad3_(fcn,f0,x,epsilon,varargin) +% Computes the gradient of the objective function fcn using a three points +% formula if possible. +% +% Adapted from Sims' numgrad routine. +% +% See section 25.3.4 in Abramovitz and Stegun (1972, Tenth Printing, December) Handbook of Mathematical Functions. +% http://www.math.sfu.ca/~cbm/aands/ + +% Original file downloaded from: +% http://sims.princeton.edu/yftp/optimize/mfiles/numgrad.m + +% Copyright (C) 1993-2007 Christopher Sims +% Copyright (C) 2008-2012 Dynare Team +% +% 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 . + +delta = epsilon; +n = length(x); +g = zeros(n,1); + +badg = 0; + +scale = []; % one(n,1); + +for i=1:n + xiold = x(i); + h = step_length_correction(xiold,scale,i)*delta; + x(i) = xiold + h; + [f1,junk1,junk2,cost_flag1] = feval(fcn, x, varargin{:}); + if ~cost_flag1 + fprintf('Gradient w.r.t. parameter number %3d (x=%16.8f,+h=%16.8f,f0=%16.8f,f1=%16.8f,f2=%16.8f,g0=%16.8f): penalty on the right!\n',i,xiold,h,f0,f1,f2,(f1 - f2) / (2*h)) + end + x(i) = xiold - h; + [f2,junk2,junk2,cost_flag2] = feval(fcn, x, varargin{:}); + if ~cost_flag2 + fprintf('Gradient w.r.t. parameter number %3d (x=%16.8f,+h=%16.8f,f0=%16.8f,f1=%16.8f,f2=%16.8f,g0=%16.8f): penalty on the left!\n',i,xiold,h,f0,f1,f2,(f1 - f2) / (2*h)) + end + if f0. + +delta = epsilon; +n = length(x); +g = zeros(n,1); + +badg = 0; +scale = []; % ones(n,1); + +for i=1:n + xiold = x(i); + h = step_length_correction(xiold,scale,i)*delta; + x(i) = xiold+h; + [f1,junk1,junk2,cost_flag1] = feval(fcn, x, varargin{:}); + x(i) = xiold-h; + [f2,junk1,junk2,cost_flag2] = feval(fcn, x, varargin{:}); + x(i) = xiold+2*h; + [f3,junk1,junk2,cost_flag3] = feval(fcn, x, varargin{:}); + x(i) = xiold-2*h; + [f4,junk1,junk2,cost_flag4] = feval(fcn, x, varargin{:}); + if f0