function [W,X,m] = pca(Y,c),
% The inputs of the function are the data matrix Y (with observations as 
% column vectors) and the number of principal components c. The outputs are 
% the matrix W, the principal components X and the bias vector m.
% 
% The principal components are normalized such that mean(X.^2,2) is a
% vector of ones. The eigenvalues can be computed as sum(W.^2,1).
% 
% See pca.pdf for more information.
%
% Alexander Ilin and Tapani Raiko 2010

check_convergence = 0;
[d,n] = size(Y);

% compute and remove the mean
m = mean(Y,2);
Y = bsxfun(@minus,Y,m);

% initialize W randomly
W = randn(d,c);

% find the principal subspace
if (check_convergence)
  err_old = 2e12;
  err = 1e12;
  while (err_old-err > 1e-12),
    X = pinv(W) * Y;
    W = Y * pinv(X);
    err_old = err;
    err = mean(mean((Y-W*X).^2));
  end
else
  for i=1:30,
    X = pinv(W) * Y;
    W = Y * pinv(X);
  end;
end;

% find the PCA axis in the principal subspace
[U, D_x] = eig(1/n*X*X');
D_x_half = diag(diag(D_x).^0.5);
[V, D_w] = eig(D_x_half*U'*W'*W*U*D_x_half);
W = W*U*D_x_half*V;
X = V'*diag(diag(D_x).^-0.5)*U'*X;

% sort principal components in decreasing importance
[temp,I] = sort(-sum(W.^2,1));
W = W(:,I);
X = X(I,:);



