Dimensionality Reduction(Eigen Vector Space): Principal Components Analysis (PCA)

PCA is designed to model linear variabilities in high-dimensional data.

Principal components analysis (PCA) is a very popular technique for dimensionality reduction. Given a set of data on n dimensions, PCA aims to find a linear subspace of dimension d lower than n such that the data points lie mainly on this linear subspace. Such a reduced subspace attempts to maintain most of the variability of the data.

The linear subspace can be specified by d orthogonal vectors that form a new coordinate system, called the ‘principal components’. The principal components are orthogonal, linear transformations of the original data points, so there can be no more than n of them.However, the hope is that only d < n principal components are needed to approximate the space spanned by the n original axes.

It is one of the most popular techniques for dimensionality reduction. However, its effectiveness is limited by its global linearity.

library("e1071")
library("RSpectra")
library("dplyr")
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# load original dataset
# image data
o_df <- read.csv("C:/Python36/my_project/project_test/olivetti_X.csv",header = F) %>% as.matrix()
# label data
y_df <- read.csv("C:/Python36/my_project/project_test/olivetti_y.csv",header = F) 
# Function to plot image data
plt_img <- function(x){ image(x, col=grey(seq(0, 1, length=256)))}

par(mfrow=c(2,2))
par(mar=c(1,1,1,1))

Display Person 10,20,30,40

# Display Person 2,3,4,5
par(mfrow=c(2,2))
par(mar=c(1,1,1,1))

c1 <- t(apply(matrix(as.numeric(o_df[100, ]), nrow=64, byrow=T), 2, rev))
plt_img(c1)
c2 <- t(apply(matrix(as.numeric(o_df[200, ]), nrow=64, byrow=T), 2, rev))
plt_img(c2)
c3 <- t(apply(matrix(as.numeric(o_df[300, ]), nrow=64, byrow=T), 2, rev))
plt_img(c3)
c4 <- t(apply(matrix(as.numeric(o_df[400, ]), nrow=64, byrow=T), 2, rev))
plt_img(c4)

Compute First 40 Eigenvalues and Eigenvectors

df <- read.csv("C:/Python36/my_project/project_test/train_faces.csv",header = T) 
D <- data.matrix(df)

# Perform PCA on the data

# Step 1: scale data
# Scale as follows: mean equal to 0, stdv equal to 1
D <- scale(D)

# Step 2: calculate covariance matrix
A <- cov(D)


# Step 3: calculate eigenvalues and eigenvectors

# Calculate the largest 40 eigenvalues and corresponding eigenvectors
eigs <- eigs(A, 40, which = "LM")
# Eigenvalues
eigenvalues <- eigs$values
# Eigenvectors (also called loadings or "rotation" in R prcomp function: i.e. prcomp(A)$rotation)
eigenvectors <- eigs$vectors

Display first 4 Eigenfaces (Eigenvectors)

# Plot of first 6 eigenfaces

par(mfrow=c(2,2))
par(mar=c(1,1,1,1))
for (i in 1:4){
  plt_img(matrix(as.numeric(eigenvectors[, i]),nrow=64,byrow=T))
}

Projection onto Eigenvector Space, reduce dimension from 4096 to 40

par(mfrow=c(2,2))
par(mar=c(2,2,2,2))

# projection of 1st photo into eigen space
# reduce the dimension from 4096 to 40
PF1 <- data.matrix(df[15,]) %*% eigenvectors
barplot(PF1,main="projection coefficients in eigen space", col="blue",ylim = c(-40,10))
legend("topright", legend = "21st photo")

# projection of 11th photo into eigen space
# reduce the dimension from 4096 to 40
PF2 <- data.matrix(df[25,]) %*% eigenvectors
barplot(PF2,main="projection coefficients in eigen space", col="green",ylim = c(-40,10))
legend("topright", legend = "31st photo")

# projection of 1st photo into eigen space
# reduce the dimension from 4096 to 40
PF3 <- data.matrix(df[35,]) %*% eigenvectors
barplot(PF3,main="projection coefficients in eigen space", col="red",ylim = c(-40,10))
legend("topright", legend = "41st photo")

# projection of 11th photo into eigen space
# reduce the dimension from 4096 to 40
PF4 <- data.matrix(df[45,]) %*% eigenvectors
barplot(PF4,main="projection coefficients in eigen space", col="purple",ylim = c(-40,10))
legend("topright", legend = "51st photo")

Reconstruct images from the reduced dimension eigen space

par(mfrow=c(2,2))
par(mar=c(1,1,1,1))

# 21st person 1st photo
# average face
average_face=colMeans(df)
AVF=matrix(average_face,nrow=1,byrow=T)

# project into eigen space and back
PF1 <- data.matrix(df[21,]) %*% eigenvectors
RE1 <- PF1 %*% t(eigenvectors)

# add the average face
RE1AVF=RE1+AVF
plt_img(matrix(as.numeric(RE1AVF),nrow=64,byrow=T))


# project into eigen space and back
PF1 <- data.matrix(df[31,]) %*% eigenvectors
RE1 <- PF1 %*% t(eigenvectors)

# add the average face
RE1AVF=RE1+AVF
plt_img(matrix(as.numeric(RE1AVF),nrow=64,byrow=T))


# project into eigen space and back
PF1 <- data.matrix(df[41,]) %*% eigenvectors
RE1 <- PF1 %*% t(eigenvectors)

# add the average face
RE1AVF=RE1+AVF
plt_img(matrix(as.numeric(RE1AVF),nrow=64,byrow=T))

# project into eigen space and back
PF1 <- data.matrix(df[51,]) %*% eigenvectors
RE1 <- PF1 %*% t(eigenvectors)

# add the average face
RE1AVF=RE1+AVF
plt_img(matrix(as.numeric(RE1AVF),nrow=64,byrow=T))

Transform to Higher Dimension(Suport Vector Space) to find Hyperplan: Support Vector Machines(SVM)

A Support Vector Machine (SVM) is a discriminative classifier formally defined by a separating hyperplane. In other words, given labeled training data (supervised learning), the algorithm outputs an optimal hyperplane which categorizes new examples. They have their roots in Statistical Learning Theory and have gained prominence because they are robust, accurate and are effective even when using a small training sample.

SVMs function by nonlinearly projecting the training data in the input space to a feature space of higher (infinite) dimension by use of a kernel function. This results in a linearly separable dataset that can be separated by a linear classifier.

Load the low dimension features for classification

x_df <- read.csv('C:/Python36/my_project/project_test/proj_faces.csv',header = T)

y_df <- lapply(y_df, factor)
dataset=cbind(y_df,x_df)
names(dataset)[1]<-paste("y")

test_dataset = dataset[seq(1, nrow(dataset), 10), ]
train_dataset = dataset[-(seq(1,to=nrow(dataset),by=10)),]

SVM project the features to higher dimension(support vectors) to build the classifier

Myclassifier = svm(formula = train_dataset$y ~ ., data = train_dataset, type = 'C-classification', kernel = 'radial') 
summary(Myclassifier)
## 
## Call:
## svm(formula = train_dataset$y ~ ., data = train_dataset, type = "C-classification", 
##     kernel = "radial")
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  1 
##       gamma:  0.025 
## 
## Number of Support Vectors:  343
## 
##  ( 9 8 9 9 9 8 9 9 8 9 9 9 9 8 9 9 8 9 7 9 9 8 9 9 8 9 8 8 9 9 8 9 7 9 9 9 8 7 9 9 )
## 
## 
## Number of Classes:  40 
## 
## Levels: 
##  0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39

Test classifier on test data

pred <- predict(Myclassifier,test_dataset[3:6,])
pred
## 21 31 41 51 
##  2  3  4  5 
## 40 Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 ... 39
test_df <- select(test_dataset,-1)
par(mfrow=c(2,2))
par(mar=c(1,1,1,1))

# 21st person 1st photo
# average face
average_face=colMeans(df)
AVF=matrix(average_face,nrow=1,byrow=T)

# project back to original dimension(4096)
RE1 <- data.matrix(test_df[3,]) %*% t(eigenvectors)

# add the average face
RE1AVF=RE1+AVF
plt_img(matrix(as.numeric(RE1AVF),nrow=64,byrow=T))

# project back to original dimension(4096)
RE1 <- data.matrix(test_df[4,]) %*% t(eigenvectors)

# add the average face
RE1AVF=RE1+AVF
plt_img(matrix(as.numeric(RE1AVF),nrow=64,byrow=T))

# project back to original dimension(4096)
RE1 <- data.matrix(test_df[5,]) %*% t(eigenvectors)

# add the average face
RE1AVF=RE1+AVF
plt_img(matrix(as.numeric(RE1AVF),nrow=64,byrow=T))

# project back to original dimension(4096)
RE1 <- data.matrix(test_df[6,]) %*% t(eigenvectors)

# add the average face
RE1AVF=RE1+AVF
plt_img(matrix(as.numeric(RE1AVF),nrow=64,byrow=T))

# system.time(pred <- predict(Myclassifier,test_dataset))
# svm_tune <- tune(svm, train.x=x, train.y=y, kernel="radial", ranges=list(cost=10^(-1:2), gamma=c(.5,1,2)))
# print(svm_tune)
# svm_model_after_tune <- svm(Species ~ ., data=iris, kernel="radial", cost=1, gamma=0.5)
# summary(svm_model_after_tune)
# pred <- predict(svm_model_after_tune,x)
# system.time(predict(svm_model_after_tune,x))

Being Idenfied as Person 2,3,4,5

# Display Person 2,3,4,5
par(mfrow=c(2,2))
par(mar=c(1,1,1,1))

c1 <- t(apply(matrix(as.numeric(o_df[30, ]), nrow=64, byrow=T), 2, rev))
plt_img(c1)
c2 <- t(apply(matrix(as.numeric(o_df[40, ]), nrow=64, byrow=T), 2, rev))
plt_img(c2)
c3 <- t(apply(matrix(as.numeric(o_df[50, ]), nrow=64, byrow=T), 2, rev))
plt_img(c3)
c4 <- t(apply(matrix(as.numeric(o_df[60, ]), nrow=64, byrow=T), 2, rev))
plt_img(c4)