wolffd@0: function [h, hdata] = rbfhess(net, x, t, hdata)
wolffd@0: %RBFHESS Evaluate the Hessian matrix for RBF network.
wolffd@0: %
wolffd@0: %	Description
wolffd@0: %	H = RBFHESS(NET, X, T) takes an RBF network data structure NET, a
wolffd@0: %	matrix X of input values, and a matrix T of target values and returns
wolffd@0: %	the full Hessian matrix H corresponding to the second derivatives of
wolffd@0: %	the negative log posterior distribution, evaluated for the current
wolffd@0: %	weight and bias values as defined by NET.  Currently, the
wolffd@0: %	implementation only computes the Hessian for the output layer
wolffd@0: %	weights.
wolffd@0: %
wolffd@0: %	[H, HDATA] = RBFHESS(NET, X, T) returns both the Hessian matrix H and
wolffd@0: %	the contribution HDATA arising from the data dependent term in the
wolffd@0: %	Hessian.
wolffd@0: %
wolffd@0: %	H = RBFHESS(NET, X, T, HDATA) takes a network data structure NET, a
wolffd@0: %	matrix X of input values, and a matrix T of  target values, together
wolffd@0: %	with the contribution HDATA arising from the data dependent term in
wolffd@0: %	the Hessian, and returns the full Hessian matrix H corresponding to
wolffd@0: %	the second derivatives of the negative log posterior distribution.
wolffd@0: %	This version saves computation time if HDATA has already been
wolffd@0: %	evaluated for the current weight and bias values.
wolffd@0: %
wolffd@0: %	See also
wolffd@0: %	MLPHESS, HESSCHEK, EVIDENCE
wolffd@0: %
wolffd@0: 
wolffd@0: %	Copyright (c) Ian T Nabney (1996-2001)
wolffd@0: 
wolffd@0: % Check arguments for consistency
wolffd@0: errstring = consist(net, 'rbf', x, t);
wolffd@0: if ~isempty(errstring);
wolffd@0:   error(errstring);
wolffd@0: end
wolffd@0: 
wolffd@0: if nargin == 3
wolffd@0:   % Data term in Hessian needs to be computed
wolffd@0:   [a, z] = rbffwd(net, x); 
wolffd@0:   hdata = datahess(net, z, t);
wolffd@0: end
wolffd@0: 
wolffd@0: % Add in effect of regularisation
wolffd@0: [h, hdata] = hbayes(net, hdata);
wolffd@0: 
wolffd@0: % Sub-function to compute data part of Hessian
wolffd@0: function hdata = datahess(net, z, t)
wolffd@0: 
wolffd@0: % Only works for output layer Hessian currently
wolffd@0: if (isfield(net, 'mask') & ~any(net.mask(...
wolffd@0:       1:(net.nwts - net.nout*(net.nhidden+1)))))
wolffd@0:   hdata = zeros(net.nwts);
wolffd@0:   ndata = size(z, 1);
wolffd@0:   out_hess = [z ones(ndata, 1)]'*[z ones(ndata, 1)];
wolffd@0:   for j = 1:net.nout
wolffd@0:     hdata = rearrange_hess(net, j, out_hess, hdata);
wolffd@0:   end
wolffd@0: else
wolffd@0:   error('Output layer Hessian only.');
wolffd@0: end
wolffd@0: return
wolffd@0: 
wolffd@0: % Sub-function to rearrange Hessian matrix
wolffd@0: function hdata = rearrange_hess(net, j, out_hess, hdata)
wolffd@0: 
wolffd@0: % Because all the biases come after all the input weights,
wolffd@0: % we have to rearrange the blocks that make up the network Hessian.
wolffd@0: % This function assumes that we are on the jth output and that all outputs
wolffd@0: % are independent.
wolffd@0: 
wolffd@0: % Start of bias weights block
wolffd@0: bb_start = net.nwts - net.nout + 1;
wolffd@0: % Start of weight block for jth output
wolffd@0: ob_start = net.nwts - net.nout*(net.nhidden+1) + (j-1)*net.nhidden...
wolffd@0:    + 1; 
wolffd@0: % End of weight block for jth output
wolffd@0: ob_end = ob_start + net.nhidden - 1; 
wolffd@0: % Index of bias weight
wolffd@0: b_index = bb_start+(j-1);   
wolffd@0: % Put input weight block in right place
wolffd@0: hdata(ob_start:ob_end, ob_start:ob_end) = out_hess(1:net.nhidden, ...
wolffd@0:    1:net.nhidden);
wolffd@0: % Put second derivative of bias weight in right place
wolffd@0: hdata(b_index, b_index) = out_hess(net.nhidden+1, net.nhidden+1);
wolffd@0: % Put cross terms (input weight v bias weight) in right place
wolffd@0: hdata(b_index, ob_start:ob_end) = out_hess(net.nhidden+1, ...
wolffd@0:    1:net.nhidden);
wolffd@0: hdata(ob_start:ob_end, b_index) = out_hess(1:net.nhidden, ...
wolffd@0:    net.nhidden+1);
wolffd@0: 
wolffd@0: return