wolffd@0: function hdv = mlphdotv(net, x, t, v)
wolffd@0: %MLPHDOTV Evaluate the product of the data Hessian with a vector. 
wolffd@0: %
wolffd@0: %	Description
wolffd@0: %
wolffd@0: %	HDV = MLPHDOTV(NET, X, T, V) takes an MLP network data structure NET,
wolffd@0: %	together with the matrix X of input vectors, the matrix T of target
wolffd@0: %	vectors and an arbitrary row vector V whose length equals the number
wolffd@0: %	of parameters in the network, and returns the product of the data-
wolffd@0: %	dependent contribution to the Hessian matrix with V. The
wolffd@0: %	implementation is based on the R-propagation algorithm of
wolffd@0: %	Pearlmutter.
wolffd@0: %
wolffd@0: %	See also
wolffd@0: %	MLP, MLPHESS, HESSCHEK
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, 'mlp', x, t);
wolffd@0: if ~isempty(errstring);
wolffd@0:   error(errstring);
wolffd@0: end
wolffd@0: 
wolffd@0: ndata = size(x, 1);
wolffd@0: 
wolffd@0: [y, z] = mlpfwd(net, x);		% Standard forward propagation.
wolffd@0: zprime = (1 - z.*z);			% Hidden unit first derivatives.
wolffd@0: zpprime = -2.0*z.*zprime;		% Hidden unit second derivatives.
wolffd@0: 
wolffd@0: vnet = mlpunpak(net, v);	% 		Unpack the v vector.
wolffd@0: 
wolffd@0: % Do the R-forward propagation.
wolffd@0: 
wolffd@0: ra1 = x*vnet.w1 + ones(ndata, 1)*vnet.b1;
wolffd@0: rz = zprime.*ra1;
wolffd@0: ra2 = rz*net.w2 + z*vnet.w2 + ones(ndata, 1)*vnet.b2;
wolffd@0: 
wolffd@0: switch net.outfn
wolffd@0: 
wolffd@0:   case 'linear'      % Linear outputs
wolffd@0: 
wolffd@0:     ry = ra2;
wolffd@0: 
wolffd@0:   case 'logistic'    % Logistic outputs
wolffd@0: 
wolffd@0:     ry = y.*(1 - y).*ra2;
wolffd@0: 
wolffd@0:   case 'softmax'     % Softmax outputs
wolffd@0:   
wolffd@0:     nout = size(t, 2);
wolffd@0:     ry = y.*ra2 - y.*(sum(y.*ra2, 2)*ones(1, nout));
wolffd@0: 
wolffd@0:   otherwise
wolffd@0:     error(['Unknown activation function ', net.outfn]);  
wolffd@0: end
wolffd@0: 
wolffd@0: % Evaluate delta for the output units.
wolffd@0: 
wolffd@0: delout = y - t;
wolffd@0: 
wolffd@0: % Do the standard backpropagation.
wolffd@0: 
wolffd@0: delhid = zprime.*(delout*net.w2');
wolffd@0: 
wolffd@0: % Now do the R-backpropagation.
wolffd@0: 
wolffd@0: rdelhid = zpprime.*ra1.*(delout*net.w2') + zprime.*(delout*vnet.w2') + ...
wolffd@0:           zprime.*(ry*net.w2');
wolffd@0: 
wolffd@0: % Finally, evaluate the components of hdv and then merge into long vector.
wolffd@0: 
wolffd@0: hw1 = x'*rdelhid;
wolffd@0: hb1 = sum(rdelhid, 1);
wolffd@0: hw2 = z'*ry + rz'*delout;
wolffd@0: hb2 = sum(ry, 1);
wolffd@0: 
wolffd@0: hdv = [hw1(:)', hb1, hw2(:)', hb2];