General Matlab help file for solving Principal Component Analysis / Proper Orthogonal Decomposition in detail
Contents
- Clear all variables, close all open windows (figures), clear screen
- File for initializing parameters
- Create the 2D grid
- Create snapshot matrix
- Generate the snapshots and store to snapshot matrix S
- Calculate mean centered data
- Proper Orthogonal Decomposition
- Calculate relative energies
- Proportion of variance explained
- Calculate the PCA modes
- Mean image
- PCA modes
- PCA reconstruction : Random snapshot
- PCA reconstruction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 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));