SANTTUcurriculum vitae
03 Sep 2018

## 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));
```

machine learning • matlab