wolffd@0: function [gradient, delta] = gradchek(w, func, grad, varargin)
wolffd@0: %GRADCHEK Checks a user-defined gradient function using finite differences.
wolffd@0: %
wolffd@0: %	Description
wolffd@0: %	This function is intended as a utility for other netlab functions
wolffd@0: %	(particularly optimisation functions) to use.  It enables the user to
wolffd@0: %	check whether a gradient calculation has been correctly implmented
wolffd@0: %	for a given function. GRADCHEK(W, FUNC, GRAD) checks how accurate the
wolffd@0: %	gradient  GRAD of a function FUNC is at a parameter vector X.   A
wolffd@0: %	central difference formula with step size 1.0e-6 is used, and the
wolffd@0: %	results for both gradient function and finite difference
wolffd@0: %	approximation are printed. The optional return value GRADIENT is the
wolffd@0: %	gradient calculated using the function GRAD and the return value
wolffd@0: %	DELTA is the difference between the functional and finite difference
wolffd@0: %	methods of calculating the graident.
wolffd@0: %
wolffd@0: %	GRADCHEK(X, FUNC, GRAD, P1, P2, ...) allows additional arguments to
wolffd@0: %	be passed to FUNC and GRAD.
wolffd@0: %
wolffd@0: %	See also
wolffd@0: %	CONJGRAD, GRADDESC, HMC, OLGD, QUASINEW, SCG
wolffd@0: %
wolffd@0: 
wolffd@0: %	Copyright (c) Ian T Nabney (1996-2001)
wolffd@0: 
wolffd@0: % Reasonable value for step size
wolffd@0: epsilon = 1.0e-6;
wolffd@0: 
wolffd@0: func = fcnchk(func, length(varargin));
wolffd@0: grad = fcnchk(grad, length(varargin));
wolffd@0: 
wolffd@0: % Treat
wolffd@0: nparams = length(w);
wolffd@0: deltaf = zeros(1, nparams);
wolffd@0: step = zeros(1, nparams);
wolffd@0: for i = 1:nparams
wolffd@0:   % Move a small way in the ith coordinate of w
wolffd@0:   step(i) = 1.0;
wolffd@0:   fplus  = feval('linef', epsilon, func, w, step, varargin{:});
wolffd@0:   fminus = feval('linef', -epsilon, func, w, step, varargin{:});
wolffd@0:   % Use central difference formula for approximation
wolffd@0:   deltaf(i) = 0.5*(fplus - fminus)/epsilon;
wolffd@0:   step(i) = 0.0;
wolffd@0: end
wolffd@0: gradient = feval(grad, w, varargin{:});
wolffd@0: fprintf(1, 'Checking gradient ...\n\n');
wolffd@0: delta = gradient - deltaf;
wolffd@0: fprintf(1, '   analytic   diffs     delta\n\n');
wolffd@0: disp([gradient', deltaf', delta'])