Daniel@0: function [h, hdata] = mlphess(net, x, t, hdata)
Daniel@0: %MLPHESS Evaluate the Hessian matrix for a multi-layer perceptron network.
Daniel@0: %
Daniel@0: %	Description
Daniel@0: %	H = MLPHESS(NET, X, T) takes an MLP network data structure NET, a
Daniel@0: %	matrix X of input values, and a matrix T of target values and returns
Daniel@0: %	the full Hessian matrix H corresponding to the second derivatives of
Daniel@0: %	the negative log posterior distribution, evaluated for the current
Daniel@0: %	weight and bias values as defined by NET.
Daniel@0: %
Daniel@0: %	[H, HDATA] = MLPHESS(NET, X, T) returns both the Hessian matrix H and
Daniel@0: %	the contribution HDATA arising from the data dependent term in the
Daniel@0: %	Hessian.
Daniel@0: %
Daniel@0: %	H = MLPHESS(NET, X, T, HDATA) takes a network data structure NET, a
Daniel@0: %	matrix X of input values, and a matrix T of  target values, together
Daniel@0: %	with the contribution HDATA arising from the data dependent term in
Daniel@0: %	the Hessian, and returns the full Hessian matrix H corresponding to
Daniel@0: %	the second derivatives of the negative log posterior distribution.
Daniel@0: %	This version saves computation time if HDATA has already been
Daniel@0: %	evaluated for the current weight and bias values.
Daniel@0: %
Daniel@0: %	See also
Daniel@0: %	MLP, HESSCHEK, MLPHDOTV, EVIDENCE
Daniel@0: %
Daniel@0: 
Daniel@0: %	Copyright (c) Ian T Nabney (1996-2001)
Daniel@0: 
Daniel@0: % Check arguments for consistency
Daniel@0: errstring = consist(net, 'mlp', x, t);
Daniel@0: if ~isempty(errstring);
Daniel@0:   error(errstring);
Daniel@0: end
Daniel@0: 
Daniel@0: if nargin == 3
Daniel@0:   % Data term in Hessian needs to be computed
Daniel@0:   hdata = datahess(net, x, t);
Daniel@0: end
Daniel@0: 
Daniel@0: [h, hdata] = hbayes(net, hdata);
Daniel@0: 
Daniel@0: % Sub-function to compute data part of Hessian
Daniel@0: function hdata = datahess(net, x, t)
Daniel@0: 
Daniel@0: hdata = zeros(net.nwts, net.nwts);
Daniel@0: 
Daniel@0: for v = eye(net.nwts);
Daniel@0:   hdata(find(v),:) = mlphdotv(net, x, t, v);
Daniel@0: end
Daniel@0: 
Daniel@0: return