function [new_data, eigvector, eigvalue] = pca(data, eigvector_n)
%PCA Principal component analysis
% Usage:
% [NEW_DATA, EIGVECTOR, EIGVALUE] = PCA(DATA, EIG_N)
%
% DATA: Rows of vectors of (zero-mean) data points
% EIG_N: No. of selected eigenvectors
% NEW_DATA: The data after projection
% EIGVECTOR: Each column of this matrix is a eigenvector of DATA'*DATA
% EIGVALUE: Eigenvalues of DATA'*DATA
%
% Note that DATA must be zero-mean'ed before calling this function.
%
% Type "pca" for a self-demo.
% Roger Jang, 970406, 990612, 991215
if nargin == 0, selfdemo; return; end
if nargin < 2, eigvector_n = min(size(data)); end
if size(data, 1) >= size(data, 2),
[eigvector, d] = eig(data'*data);
eigvalue = diag(d);
% ====== Sort based on descending order
[junk, index] = sort(-eigvalue);
eigvalue = eigvalue(index);
eigvector = eigvector(:, index);
new_data = data*eigvector;
if eigvector_n < size(data, 2),
eigvalue = eigvalue(1:eigvector_n);
eigvector = eigvector(:, 1:eigvector_n);
new_data = new_data(:, 1:eigvector_n);
end
else % This is an efficient method which computes the eigvectors of
% of A*A^T (instead of A^T*A) first, and then convert them back to
% the eigenvectors of A^T*A.
[eigvector1, d] = eig(data*data');
eigvalue = diag(d);
% ====== Sort based on descending order
[junk, index] = sort(-eigvalue);
eigvalue = eigvalue(index);
eigvector1 = eigvector1(:, index); % Eigenvectors of A*A^T
eigvector = data'*eigvector1; % Eigenvectors of A^T*A
eigvector = eigvector*diag(1./(sum(eigvector.^2).^0.5)); % Normalization
new_data = data*eigvector;
if eigvector_n < size(data, 1),
eigvalue = eigvalue(1:eigvector_n);
eigvector = eigvector(:, 1:eigvector_n);
new_data = new_data(:, 1:eigvector_n);
end
end
function selfdemo
% ====== Demo for 2D data
data_n = 1000;
x = 5+randn(data_n, 1);
y = 10+randn(data_n, 1)/2;
x_mean = mean(x);
y_mean = mean(y);
theta = pi/6;
z = ((x-x_mean) + i*(y-y_mean))*exp(i*theta)+x_mean+i*y_mean;
x = real(z);
y = imag(z);
plot(x, y, '.');
axis image;
t = [x-x_mean y-y_mean];
[new_data, v, d] = feval(mfilename, t);
v1 = -v(:, 1);
v2 = -v(:, 2);
arrow_x = [-1 0 nan -0.1 0 -0.1]+1;
arrow_y = [0 0 nan 0.1 0 -0.1];
arrow = arrow_x + j*arrow_y;
arrow1 = 2*arrow*(v1(1)+i*v1(2))*d(1)/data_n+x_mean+i*y_mean;
arrow2 = 2*arrow*(v2(1)+i*v2(2))*d(2)/data_n+x_mean+i*y_mean;
line(real(arrow1), imag(arrow1), 'color', 'c', 'linewidth', 2);
line(real(arrow2), imag(arrow2), 'color', 'g', 'linewidth', 2);
title('Axes for PCA');
% ====== Demo for Iris data
load iris.dat
data = iris(:, 1:4); % Take only the input part
data_n = size(data, 1);
data = data - ones(data_n, 1)*mean(data); % Make data zero-mean
[new_data, eigvector, eigvalue] = feval(mfilename, data);
index1 = find(iris(:,5)==1);
index2 = find(iris(:,5)==2);
index3 = find(iris(:,5)==3);
figure;
plot( new_data(index1, 1), new_data(index1, 2), '*', ...
new_data(index2, 1), new_data(index2, 2), 'o', ...
new_data(index3, 1), new_data(index3, 2), 'x');
legend('Class 1', 'Class 2', 'Class 3');
title('PCA projection of IRIS data onto the first 2 eigenvectors');
axis equal; axis tight;
figure;
plot( new_data(index1, 3), new_data(index1, 4), '*', ...
new_data(index2, 3), new_data(index2, 4), 'o', ...
new_data(index3, 3), new_data(index3, 4), 'x');
legend('Class 1', 'Class 2', 'Class 3');
title('PCA projection of IRIS data onto the last 2 eigenvectors');
axis equal; axis tight;
