Initialization for symmetric NMF using Eigendecomposition based on NNDSVD technique by Boutsidis & Gallopoulos (still experimental as the math may not be correct)

NNDEIG.m

function [W] = NNDEIG(A,k,flag);
%
% This function implements the NNDSVD algorithm described in [1] for
% initialization of Nonnegative Matrix Factorization Algorithms
% for symmetric NMF so uses Eigendecomposition
%
% [W] = nndeig(A,k,flag);
%
% INPUT
% ------------
%
% A    : the input nonnegative n x n matrix A
% k    : the rank of the computed factors W
% flag : indicates the variant of the NNDSVD Algorithm
%        flag = 0 --> NNDEIG
%        flag = 1 --> NNDEIGa
%        flag = 2 --> NNDEIGar
%
% OUTPUT
% -------------
%   
% W   : nonnegative n x k matrix
%
% 
% References:
% 
% [1] C. Boutsidis and E. Gallopoulos, SVD-based initialization: A head
%     start for nonnegative matrix factorization, Pattern Recognition,
%     Elsevier
%
% This original NNDSVD code was kindly provided by the authors for research porpuses.
% - Efstratios Gallopoulos (stratis@ceid.upatras.gr)
% - Christos Boutsidis (boutsc@cs.rpi.edu)     
% 
%--------------------------------------------------------------------------


%----------------------check the input matrix------------------------------
if numel(find(A<0)) > 0
    error('The input matrix contains negative elements !')
end
%--------------------------------------------------------------------------

%size of the input matrix
[n,m] = size(A);

%the matrices of the factorization
W = zeros(n,k);

% 1st SVD --> partial SVD rank-k to the input matrix A. 
%[U,S,V] = svds(A,k);
[U,S] = eigs(A) % original [V,D] changed to [U,S]

%choose the first singular triplet to be nonnegative
W(:,1)     =  sqrt(S(1,1)) * abs(U(:,1) );         

% 2nd SVD for the other factors (see table 1 in our paper)
for i=2:k
    uu = U(:,i); 
    %vv = V(:,i);
    uup = pos(uu); uun = neg(uu) ;
    %vvp = pos(vv); vvn = neg(vv);
    n_uup = norm(uup);
    %n_vvp = norm(vvp) ;
    n_uun = norm(uun) ;
    %n_vvn = norm(vvn) ;
    %termp = n_uup*n_vvp; termn = n_uun*n_vvn;
    termp = n_uup*n_uup; 
    % not sure what to do here
    %if (termp >= termn)
    %    W(:,i) = sqrt(S(i,i)*termp)*uup/n_uup; 
    %    H(i,:) = sqrt(S(i,i)*termp)*vvp'/n_vvp;
    %else
    %   W(:,i) = sqrt(S(i,i)*termn)*uun/n_uun; 
    %    H(i,:) = sqrt(S(i,i)*termn)*vvn'/n_vvn;
    %end
    W(:,i) = sqrt(S(i,i)*termp)*uup/n_uup; 
end
%------------------------------------------------------------

%actually these numbers are zeros
W(find(W<0.0000000001))=0.1;
%H(find(H<0.0000000001))=0.1;

% NNDSVDa: fill in the zero elements with the average 
if flag==1
   ind1      =  find(W==0) ;
   %ind2      =  find(H==0) ;
   average   =  mean(A(:)) ; 
   W( ind1 ) =  average    ; 
   %H( ind2 ) =  average    ;

% NNDSVDar: fill in the zero elements with random values in the space [0:average/100]
elseif flag==2
   ind1      =  find(W==0) ;
   %ind2      =  find(H==0) ;
   n1        =  numel(ind1);
   %n2        =  numel(ind2);
   
   average   =  mean(A(:))       ;
   W( ind1 ) =  (average*rand(n1,1)./100)  ; 
   %H( ind2 ) =  (average*rand(n2,1)./100)  ;   
end

end
%--------------------------------------------------------------------------
%end of the nndsvd function



%This function sets to zero the negative elements of a matrix
%--------------------------------------------------------------------------
function [Ap] = pos(A)
Ap = (A>=0).*A;
end
%--------------------------------------------------------------------------

%This functions sets to zero the positive elements of a matrix and takes
%the absolute value of the negative elements
%--------------------------------------------------------------------------
function [Am] = neg(A);
Am = (A<0).*(-A);
end
%-------------------------------------------------------------------------

It seems it can't find W matrix if k > 6. Do you have any suggestions on how to fix this?

I am getting the following error for k > 6.

Index in position 2 exceeds array bounds (must not exceed 6).

Error in NNDEIG (line 59)
uu = U(:,i);