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