% Principal Component Analysis without Tears
%
% This is a lecture on PCA with Matlab code
% Xiaojin Zhu, jerryzhu@cs.wisc.edu
% Oct. 2017

% read vocabulary
vocab=fileread('vocabulary_stopword_removal.txt');
vocab=string(vocab);
vocab=splitlines(vocab);

% read BOW from the WARC201709 corpus
load WARC201709BOW.txt
x=sparse(WARC201709BOW(:,1), WARC201709BOW(:,2), WARC201709BOW(:,3))';
[d,n]=size(x)

% center x
mu=mean(x,2);
x=x-repmat(mu,1,n);


% PCA goal: find w to maximize the variance (preserve differences)
% max_w 1/(n-1) sum_i (w'x_i)^2 
%       = w' S w
% s.t. w'w=1
% where S=1/(n-1) sum x x' is the covariance matrix of the centered x
S=x*x'/(n-1);

% Lagrangian w'Sw + lambda (1-w'w)
% grad = 1/2 Sw - 1/2 lambda w = 0
% Sw=lambda w, w is an eigenvector of S!
% w'Sw = lambda w'w = lambda => want largest eigenvalue.
[evec, eval]=eigs(S,1);

% if instead of a single projection direction, we want a M dimensional projection space,
% these are the largest M eigenvectors/values of S.
% S \approx U Lambda U'
M=2;
[U, Lambda]=eigs(S,M);
var(U(:,1)'*x, 1)
var(U(:,2)'*x, 1)
Lambda
% why variance=LAMBDA? var=w'Sw = lambda w'w = lambda

% what projection directions do we get?
dim=1
[tmpy tmpi]=sort(abs(U(:,dim)), 1, 'descend');
[U(tmpi(1:20), dim) vocab(tmpi(1:20))]
dim=2
[tmpy tmpi]=sort(abs(U(:,dim)), 1, 'descend');
[U(tmpi(1:20), dim) vocab(tmpi(1:20))]


% geometric view: M major axes
% M can at most be D, then M is a *rotated* coordinate system aligned with major axes

% for any orthonomal basis u1...uD, yi can be written as
% 	yi = sum_{j=1}^D alpha_ij uj
% where
%	alpha_ij = uj' yi.

% 3. For xi, we project it to the M dimensional space as alpha_i{1..M}
% where
%	alpha_ij = uj' xi.
% This is dimensionality reduction if M<D
alpha=U'*x;

% 4. To recover the original xi, do
%    xi \approx sum_{j=1}^M alpha_ij uj 
xapprox = U*alpha;

%% usage of PCA: visualization
figure; hold on
for i=1:n,
  text(alpha(1,i), alpha(2,i), sprintf('%03d',i));
end
axis([ min(alpha(1,:)), max(alpha(1,:)), min(alpha(2,:)), max(alpha(2,:)) ])
xlabel('\alpha1')
ylabel('\alpha2')

% PCA is an *unsupervised* technique.
% When it can go wrong: classification (supervised learning).
% In that case, the corresponding method is Fisher's Linear Discriminant Analysis.
