% PCA demo

p = 10; % number of training points
dimx = 3; % dimension of the data
xdata = randn(dimx,p) + 2*ones(dimx,p); 
A = rand(dimx); xdata = A*xdata; % generate some training data, xdata

% perform PCA :
m = mean(xdata,2); x = xdata - repmat(m,1,p);  % subtract the mean
S = cov(x'); % covariance of the data

[Evec,Evalm] = eig(S); Eval = diag(Evalm); % find the e-vals and e-vecs 
[evals,index]=sort(Eval); index=flipud(index); % find the largest e-vals
num_retain = 1; % number of eigenvectors to retain
Evec_retained = Evec(:,index(1:num_retain)); % get the largest e-vecs

x_lowerdim = Evec_retained'*x; % lower dimensional representation 
x_reconstruction = Evec_retained*x_lowerdim + repmat(m,1,p); % reconstruction of x

subplot(1,2,1); imagesc(xdata); subplot(1,2,2); imagesc(x_reconstruction);
