function ddEW = lyngby_nn_qddews(H, W);

% lyngby_nn_qddews     - 2nd. order der., quadratic, output, gauss approx.
%
%	function ddEW = lyngby_nn_qddews(H, W);
%
%	Input:	H      Hidden layer
%               W      Output weights
%
%	Output:	ddEW   Derivative, Hessian matrix 
%
%       Gauss-Newton approximation for the second order derivative of
%       the quadratic errofr function for the output weight in the
%       neural network.

% cvs : $Id: lyngby_nn_qddews.m,v 1.3 1998/02/12 15:17:38 fnielsen Exp $
%       $Revision: 1.3 $


    [Np, Nh] = size(H);
    [Nhh, No] = size(W);
    if Nhh > Nh
      H = [H ones(Np,1)];
    end

    ddEW = 1/(Np*No) * kron(ones(No,No), H'*H); 

    % Numerical errors? 
    indices = find(W==0);
    if indices
      ddEW(indices,:) = 0;
      ddEW(:,indices) = 0;
    end




