General Matlab help file for solving Principal Component Analysis / Proper Orthogonal Decomposition in detail

Contents

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FILE:   PCA_Gaussian.m
% AUTHOR: Santosh Tirunagari, Feb 25/2012
% DESCRIPTION:
% General matlab help file for solving Principal Component Analysis / proper orthogonal decomposition
% with images of Gaussian processes.
% Also contains usage of quiver plot to view velocity vectors,
% Eigen decompostion method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Clear all variables, close all open windows (figures), clear screen

clear all
close all
clc

File for initializing parameters

N  = 2^7;     % number of cells in a coordinate direction; 2^n -> best efficiency
L  = 2*pi;    % domain side length, do not change
dx = L/N;     % grid spacing, do not change
dy = L/N;     % always dx=dy, do not change

nos = 300;    % Number of snapshots

Create the 2D grid

x = (L-dx)*(0:(N-1))/(N-1); % to be precise, dx is subtracted here
y = (L-dx)*(0:(N-1))/(N-1); % otherwise the derivative of e.g. sin(x) wouldnt be differentiable on this periodic grid
[X,Y] = meshgrid(x,y);      % basically, these grid points do not really affect the solution algorithm

Create snapshot matrix

[m,n] = size(X); % take the size of the meshrid to reshape mxn matrix to m*nx1 column vector
S = []; % S is the snapshot matrix, to store the snapshots or the velocity components

Generate the snapshots and store to snapshot matrix S

 for i = 1:1:nos

    alpha = (rand)*(pi./2); X1=X-pi; Y1 = Y - pi -alpha;
    u1 = (Y1.*exp(-((X1).^2+(Y1).^2)/(2*(pi/6)^2)));
    v1 = -(X1.*exp(-((X1).^2+(Y1).^2)/(2*(pi/6)^2)));

     alpha = (rand)*(pi./2); X1=X-pi-alpha; Y1 = Y - pi;
    u2 = (Y1.*exp(-((X1).^2+(Y1).^2)/(2*(pi/6)^2)));
    v2 = -(X1.*exp(-((X1).^2+(Y1).^2)/(2*(pi/6)^2)));

    alpha = (rand)*(pi./2); X1=X-pi; Y1 = Y - pi + alpha;
    u3 = (Y1.*exp(-((X1).^2+(Y1).^2)/(2*(pi/6)^2)));
    v3 = -(X1.*exp(-((X1).^2+(Y1).^2)/(2*(pi/6)^2)));

    alpha = (rand)*(pi./2); X1=X-pi+alpha; Y1 = Y - pi;
    u4 = (Y1.*exp(-((X1).^2+(Y1).^2)/(2*(pi/6)^2)));
    v4 = -(X1.*exp(-((X1).^2+(Y1).^2)/(2*(pi/6)^2)));


    u = u1+u2+u3+u4; v = v1+v2+v3+v4;

    U = reshape(u,m*n,1);
    V = reshape(v,m*n,1);
   % fprintf('%d\n',i);
    temp =[U;V];
    S = [S temp];
end

clc

Calculate mean centered data

mean_all=mean(S,2);
dbx=[];
dbx=S -repmat(mean_all,1,nos);
dbx(isnan(dbx)) = 0;

Proper Orthogonal Decomposition

R=(dbx'*dbx)./nos; % calculate the covariance matrix
[eV,D]=eig(R);
% solve: eV is eigenvectors, D is eigenvalues in diagonal matrix
[L,I]=sort(diag(D));
%sort eigenvalues in ascending order - I is sorted index vector
for i=1:length(D)
eValue(length(D)+1-i)=L(i);
% Eigenvalues sorted in descending order
eVec(:,length(D)+1-i)=eV(:,I(i)); % Eigenvectors sorted in the same order
end;
eValue(length(eValue))=0;
%last eigenvalue should be zero

Calculate relative energies

menergy=eValue/sum(eValue); % relative energy associated with mode m
fh = figure();
plot(menergy,'r-*')
xlabel('Number of PCA eigenvalues')
ylabel('Energy captured')
xlim([1 18]), ylim([0 0.6])
grid on;

Proportion of variance explained

normalised_evalues = eValue ./ sum(eValue);
figure, plot(cumsum(normalised_evalues),'k-*');
xlabel('No. of eigenvectors');
ylabel('Variance accounted for');
xlim([1 18]), ylim([0 1]), grid on;

Calculate the PCA modes

for i=1:50
tmp=dbx*eVec(:,i);
% find mode
phi(:,i)=tmp/norm(tmp);
% normalize mode
end;

Mean image

k = 1;
d = size(mean_all);
d1 = d/2;
d4= mean_all(1:d1,1);
 d5 = mean_all((d1+1):end,1);
 vel = d4.^2 + d5.^2;
 s = reshape(vel,m,n);
  imagesc(s)

PCA modes

for k = 1:1:18
  d = size(phi(:,k));
 d1 = d/2;
 d4= phi(1:d1,k);
 d5 = phi((d1+1):end,k);
 vel = d4.^2 + d5.^2;
 %
  s = reshape(vel,n,m);
  imagesc(s)
  pause(2);
 end

PCA reconstruction : Random snapshot

NormImage = S(:,15); % random snapshot
% visualize random snapshot
d4= NormImage(1:d1,1);
d5 = NormImage((d1+1):end,1);
n_vel = d4.^2 + d5.^2;
figure(1), imagesc(reshape(n_vel,m,n));

PCA reconstruction

P is the vector for storing dotproducts between randomsnapshot and poms

p = [];

k = 10; % K is the number of PCA modes we are interested
for i = 1:1:k   %size(u,2)


    pare = dot(NormImage,phi(:,i));
    p = [p; pare];

end
ReshapedImage =  (phi(:,1:k)*p)+mean_all;
d4= ReshapedImage(1:d1,1);
d5 = ReshapedImage((d1+1):end,1);
r_vel = d4.^2 + d5.^2;
figure(2), imagesc(reshape(r_vel,m,n));

matlabPCAPODPrinciple component analysis