function ddEW = lyngby_nn_eddewds(H, W);

% lyngby_nn_qddewds    - 2nd order der, quadratic, output, diag gauss
%
%       function ddEW = lyngby_nn_eddewds(H, W);
%
%	Input:	H      Hidden layer
%               W      Output weights
%
%	Output:	ddEW   Derivative, Matrix with the hessian diagonal
%
%       Diagonal Gauss-Newton approximation for the second order
%       derivative of the quadratic error function with respect to the
%       output weights in the neural network.

% cvs : $Id: lyngby_nn_qddewds.m,v 1.1 1997/08/22 16:43:12 fnielsen Exp $
%       $Revision: 1.1 $

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

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

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


