function ddEV = lyngby_nn_qddevs(X, V, W, H12)

% lyngby_nn_qddevs     - 2nd der., quadratic, input, gauss approx.
%
%       function ddEV = lyngby_nn_qddevs(X, V, W, H12)
%
%	Input:  X     Neural network input
%               V     Input weights
%               W     Output weights
%               H12   (1-tanh(XV).^2) 
%
%       Output: ddEV  Derivative as a Hessian matrix
%
%       Gauss-Newton (symmetric) approximation to the second order
%       derivative of the qaudratic error function with respect to the
%       input weights in the neural network.

% cvs : $Id: lyngby_nn_qddevs.m,v 1.3 1998/02/12 14:08:10 fnielsen Exp $
%       $Revision: 1.3 $


    [Np, Nii]  = size(X);
    [Nii, Nh] = size(V);
    [Nhh, No] = size(W);

    if No == 1
      W1 = W(1:Nh, :)';
    else
      W1 = sum(W(1:Nh, :)');
    end

    A = kron(ones(1,Nh), X) .* kron(H12 .* (ones(Np,1) * W1), ones(1, Nii));
    ddEV = A'*A / (Np*No);
    
    if 1==0
      ddEVtest = zeros(Nii*Nh);
      for p=1:Np
	for o=1:No
	  for h1=1:Nh
	    for h2 = 1:Nh
	      for i1 = 1:Nii
		for i2 = 1:Nii
		  n1 = i1+Nii*(h1-1);
		  n2 = i2+Nii*(h2-1);
		  ddEVtest(n1,n2) = ddEVtest(n1,n2) + ...
		      1/(Np*No) * W(h1,o)*W(h2,o) * H12(p,h1)*H12(p,h2) ...
		      * X(p,i1)*X(p,i2);
		end
	      end
	    end
	  end
	end
      end
    end

