# Problems with implementing PCA in Matlab

I'm implementing PCA using eigenvalue decomposition in Matlab. I know Matlab has PCA implemented, but it helps me understand all the technicalities when I write code. I've been following the guidance from here, but I'm getting different results in comparison to built-in function princomp.

Could anybody look at it and point me in the right direction.

Here's the code:

function [mu, Ev, Val ] = pca(data)

% mu - mean image
% Ev - matrix whose columns are the eigenvectors corresponding to the eigen
% values Val
% Val - eigenvalues

if nargin ~= 1
error ('usage: [mu,E,Values] = pca_q1(data)');
end

mu = mean(data)';

nimages = size(data,2);

for i = 1:nimages
data(:,i) = data(:,i)-mu(i);
end

L = data'*data;
[Ev, Vals]  = eig(L);
[Ev,Vals] = sort(Ev,Vals);

% computing eigenvector of the real covariance matrix
Ev = data * Ev;

Val = diag(Vals);
Vals = Vals / (nimages - 1);

% normalize Ev to unit length
proper = 0;
for i = 1:nimages
Ev(:,i) = Ev(:,1)/norm(Ev(:,i));
if Vals(i) < 0.00001
Ev(:,i) = zeros(size(Ev,1),1);
else
proper = proper+1;
end;
end;

Ev = Ev(:,1:nimages);


The line [Ev,Vals] = sort(Ev,Vals); probably does not do what you think it should. The sort command would have to be applied separately. Moreover, eig returns a diagonal matrix of eigenvalues, which you have to strip out to a vector. Probably you want to do the following:

[Ev, Vals]  = eig(L);
Vals = diag(Vals);        %strip out the eigenvalues;
[Vals,sidx] = sort(Vals); %do the sort
Ev = Ev(:,sidx);          %apply the sort to the eigenvectors as well.


I would guess this is actually unecessary, because eig returns the eigenvectors, -values, in ascending order anyway.

On a stylistic note, you can subtract out the mean using bsxfun as follows:

X = bsxfun(@minus,X,mu);


(This was supposed to be a comment, but it got too long)

I'd say

L = data'*data;
[Ev, Vals]  = eig(L);
[Ev,Vals] = sort(Ev,Vals);


is pretty dodgy; it's a common notion in numerical linear algebra that forming matrices like $\mathbf A^T\mathbf A$ worsens the so-called "conditioning" of your problem, in that merely forming the matrix can cause the loss of significant figures, and any further operations with it, like trying to obtain the required eigenpairs, will result in a complete mess.

One possible way to show how "squaring" ruins things is to consider the so-called Lauchli matrix, which is returned by the MATLAB function gallery('lauchli',n), which returns an $(n+1)\times n$ matrix. If you attempt to square that and give that to an eigenroutine like eig()... disaster! That's just for a small matrix; the larger the matrix you have, the more likely the conditioning of your matrix will be rather bad.

The best way of doing things, then, is to use the singular value decomposition. I discussed SVD and the Lauchli example in this m.SE answer, but for your purposes, here's what you might try replacing that snippet I quoted with (following shabbychef's answer):

[U, Vals, Ev]  = svd(data, 'econ');
Vals = diag(Vals);
Vals = Vals.^2
[Vals, sidx] = sort(Vals);
Ev = Ev(:,sidx);


On the other hand, svd() returns the singular values in decreasing order, so the sorting might no longer be required...

• If the MATLAB looks wonky, I apologize; I haven't touched MATLAB in years... XD – J. M. is not a statistician Dec 12 '10 at 13:08