function [ W, H ] = NMFinitPCA( A, K )
%
% Corded by D. Kitamura (d-kitamura@nii.ac.jp) on 2 Jan. 2016.
%
% This function calculates initial basis and coefficient matrices
% of NMF based on principal component analysis.
% This initialization method has been proposed in the following literature:
% L. Zhao, G. Zhuang and X. Xu, "Facial expression recognition
% based on PCA and NMF," Proc. WCICA, pp. 6826?6829, 2008.
%
% see also
% http://d-kitamura.sakura.ne.jp/index.html
%
% [syntax]
%   [ W, H ] = NMFinitPCA( A )
%   [ W, H ] = NMFinitPCA( A, K )
%
% [inputs]
%      A: input nonnegative matrix (M x N)
%      K: number of bases (nonnegative scalar, default: N/10)
%
% [outputs]
%      W: initialized basis matrix (M x K)
%      H: initialized coefficient matrix (K x N)
%

% Check errors and set default values
[M,N] = size(A);
if min(min(A)) < 0
    error('All the entries of input matrix must be nonnegative.');
end
if (nargin < 2)
    K = N/10;
end

[H,W] = local_PCA(A,K); % PCA
W = max(abs(W.'),eps); % Make basis matrix nonnegative
H = max(abs(H),eps); % Make coefficient matrix nonnegative
end

%%% subfunction for PCA %%%%
function [ Y, Z ] = local_PCA( A, K )
[M,N] = size(A);
cX = A - repmat(mean(A,2),1,N); % Data centering
V = cX*cX.'/N; % Covariance matrix (M x M) (equivalent to V = cov(centerX.'))
[P,D] = eig(V); % Eigenvalue decomposition

% Sort eigenvalues in descending order
eig_val = diag(D);
[eig_val,idx] = sort(eig_val, 'descend');
D = D(idx,idx);
P = P(:,idx);

% Pick up top d eigenvalues and their eigenvectors
D2 = D(1:K,1:K);
P2 = P(:,1:K); 

Y = P2'*A; % Output matrix
Z = P2';  % Transformation matrix
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EOF %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%