Principal Component Analysis (PCA) - An Easy Tutorial with Python

Figure 1: MNIST Database: Handwritten digits

In a typical machine learning or data analysis problem, the total number of variables or features are usually large and that makes the visualization of the data (which is often the first step in any data analysis) is somewhat challenging. There are many visualizing tools to explore univaruiate or bivariate data, but the plots to reveal the association between many variables in higher dimension is complicated.

To understand better, suppose there are 𝑝 features 𝑋1, 𝑋2, … ,𝑋𝑝. Then we can examine all pairwise scatterplots of the data, each of which contains the plots of 𝑛 observations on two of the features. But that will result in a 𝑝(𝑝−1)/2 scatter plots. For 𝑝 = 10, this means we need to check 45 scatter plots. Now consider an image data with only 28x28 pixels. Then an image has 784 features variables, or 𝑝 = 784 and we need to explore 306,936 scatter plots. This is an impossible task and more importantly, most of them will be very uninformative as they contain only a very small fraction of the total information hidden in the data. Clearly, a better method is required to visualize the 𝑛 observations when 𝑝 is large. In particular, we would like to find a low-dimensional representation of the data that captures as much of the information as possible. For instance, if we can obtain a two-dimensional representation of the data that captures most of the information, then we can plot the observations in this low-dimensional space.

In this article, we discuss some powerdul techniques in unsupervised learning known as dimension reduction. The general idea is to reduce the dimension of the dataset while preserving important characteristics, such as the distance between features or observations. With fewer dimensions, visualization then becomes more feasible. We will restrict our discussion to a linear dimesion reduction tool, Principal component analysis (PCA).

In the following, we first motivate PCA with a simple example.

Preserving the Distance

Let us consider an artificial example of marks of 100 students in Maths and English tests. Some of these students are from a priviledged background and some are from a less priviledged background. We simulate 100 two dimensionals points that represent the number of standard deviations each individual is from the mean score. Each points represents the scores in Maths and English for a student.

#Load the libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

np.random.seed(42)

Sigma = np.array([[9.0, 9 * 0.9], [9 * 0.9, 9.0]])
mu1 = [69, 69]
mu2 = [55, 55]
n = 100
x1,y1 = np.random.multivariate_normal(mu1, Sigma, int(n/2)).T
x2, y2 = np.random.multivariate_normal(mu2, Sigma, int(n/2)).T

x = np.append(x1,x2)
y = np.append(y1,y2)

X = np.vstack((x,y)).T

A scatterplot quickly reveals that the correlation between x and y is high and that there are two groups of students, the more priviledged (upper right points) and the less proviledged (lower left points):

plt.plot(x,y,'.', ms=10.0)
plt.xlabel('Maths Scores')
plt.ylabel('English Scores')
plt.show()

Figure 2: Scatter plot of Simulated Maths and English Scores

In this example, we have only 2 features. Let us consider it is very challenging to visualise even a two-diumensional data and we wish to visualise only in one dimension. We therefore want to reduce the dimensions from two to one, but still be able to understand the important characteristics of the data, for example that the observations cluster into two groups.

Suppose we want a one-dimensional summary of these two features from which we can approximate the distance between any two observations. In the plot below we show the distance between observation 1 and 2 (magenta), and observation 1 and 51 (red). Note that the blue line is shorter, which implies 1 and 2 are closer.

plt.plot(x,y,'.', ms=10.0)
plt.plot(x[:2], y[:2], 'm-', lw=2.0)
plt.plot(x[[0,50]], y[[0,50]], 'r-', lw=2.0)
plt.xlabel('Maths Scores')
plt.ylabel('English Scores')
plt.show()

Figure 3: Scatterplot with distances between two points

We can compute these distances using pdist

from scipy.spatial.distance import pdist, squareform

dist = pdist(X)
distm = squareform(dist) #In a matrix form
print('Distance between obs 1 and 2:', distm[0,1])
print('Distance between obs 1 and 51:', distm[0,50])

Distance between obs 1 and 2: 1.6951912225188088
Distance between obs 1 and 51: 11.895125409898515

This distance is based on two dimensions and we need a distance approximation based on just one.

Let us start with the naive approach of simply removing one of the two dimensions. Let us compare the actual distances to the distance computed with just the first dimension of Maths scores only:

distx = pdist(np.asmatrix(x).T)
plt.plot(dist, distx, '.', ms=5.0)
plt.plot([0,25],[0,25],'-r')
plt.xlabel('Distances Based on Two Dimensions')
plt.ylabel('Distances Based on One Dimension')
plt.show()

Figure 4: Ditances between the points based on two dimensions and one dimension

Here are the approximate distances versus the original distances. The plot looks almost the same if we use English scores instead of the Maths scores. In geneeral, we underestimate the actual distances with one dimension. This is to be expected because we are adding more positive quantities in the distance calculation as we increase the number of dimensions.

However, if we divide the two dimensional distances by √2, this underestimation goes away.

plt.plot(dist/np.sqrt(2), distx, '.', ms=5.0)
plt.plot([0,25],[0,25],'-r')
plt.xlabel('Distances Based on Two Dimensions/sqrt(2)')
plt.ylabel('Distances Based on One Dimension')
plt.show()

Figure 5: Scaled distances between the points in one dimension and two dimensions

Now, the question is can we make it even better?

If we look back at the previous scatterplot and visualize a line between any pair of points, the length of this line is the distance between the two points. These lines tend to go along the direction of the diagonal. Notice that if we instead plot the difference versus the average:

Z = np.vstack(((y + x)/2,  (y-x))).T

plt.plot(Z[:,0],Z[:,1],'.', ms=10.0)
plt.plot(Z[:2,0], Z[:2,1], 'm-', lw=2.0)
plt.plot(Z[[0,50],0], Z[[0,50],1], 'r-', lw=2.0)
plt.xlabel('Average of two scores')
plt.ylabel('Difference between two scoress')
plt.ylim([-15,15])
plt.show()

Figure 6: Scatter plot of transformed points

We can see how the distance between points is mostly explained by the first dimension: the average. This means that we can ignore the second dimension and not lose too much information. If the line is completely horizontal, we do not lose any information at all. Using the first dimension of this transformed data we obtain an even better approximation:

distz = pdist(np.asmatrix(Z[:,0]).T)
plt.plot(dist/np.sqrt(2), distz, '.', ms=5.0)
plt.plot([0,28],[0,28],'-r')
plt.xlabel('Distances Based on Two Dimensions/sqrt(2)')
plt.ylabel('Distances Based on Averages')
plt.show()

Figure 7: Distances between the points in the transformed coordinates

Linear Transformations

To understand the Mathematical principles, let us check what is heppening here. Observe that the rows of the transformed data matrix 𝑍 is obtained using a linear transformation of the rows of 𝑋.

We can also use linear transformation of 𝑍 to get back 𝑋 using

If we define a matrix

Then 𝑍=𝑋𝐴 and transforming back we can write

Dimension reduction can often be described as applying a transformation 𝐴 to a data matrix 𝑋 with many columns that moves the information contained in 𝑋 to the first few columns of 𝑍=𝑋𝐴, then we keep only these few informative columns, which reduces the dimension of the vectors contained in the rows.

Orthogonal transformations

Earlier we had to divide by √2 to account for the differences in dimensions when comparing a 2 dimensional distance to an one dimensional distance. To ensure that the distance scales remain the same in all dimensions, we can re-scale the columns of the transformation matrix 𝐴, so that the length of each column of 𝐴 is 1. That is

Further, if the columns of the transfomed data 𝑍 are uncorrelated, then the ignored columns do not contain any further information about the columns included and that makes our dimension reduction useful. To ensure correlation to be zero we must have

In our example, to achieve orthogonality, we multiply the first set of coefficients (first column of 𝐴) by √2 and the second by 1/√2, then we get the same exact distance if we use both dimensions:

Z = np.vstack(((y + x)/np.sqrt(2),  (y-x)/np.sqrt(2))).T

distz = pdist(Z)

np.std(dist - distz)

8.642993477664565e-15

If we use only the first column of 𝑍, then standard deviation of the differences between the distances based on 2 dimensions and the one dimension is:

distz1 = pdist(np.asmatrix(Z[:,0]).T)
np.std(dist - distz1)

0.3831326110045057

This gives us an improved approximation of the distances.

In this case, 𝑍 is called an orthogonal rotation of 𝑋 and it preserves the distances between rows.

Note that by using the transformation above we can summarize the distance between any two pairs of test scores with just one dimension. For example, one-dimensional data exploration of the first dimension of 𝑍 clearly shows that there are two groups:

plt.hist(Z[:,0], bins=20)
plt.xlabel('Transformed Dimension 1')
plt.ylabel('Count')
plt.show()

Figure 8: Histogram of the transformed observations in one dimension

Observe that we have here reduced the number of dimensions from two to one with very little loss of information.

The reason we were able to do this is because the columns of 𝑋 were highly correlated:

np.corrcoef(X[:,0], X[:,1])[0,1]

0.9844468428261011

The orthogonal transformation produces uncorrelated columns of 𝑍. (In this case, the correlation is not exactly 0, but very small).

np.corrcoef(Z[:,0], Z[:,1])[0,1]

-0.1105356767780309

Principal Component Analysis

Now, we formalize the dimesion rediction using principal components.

Let us define the total variability in our data as the sum of sum of squares of each of the columns. If each of the colmns of 𝑋 are centered by its average, then the total variablity is 𝑣1+𝑣2, where

print(np.mean(X**2, axis=0))
np.sum(np.mean(X**2, axis=0))

[3940.74842046 3944.23285603]
7884.981276482109

Therefore, we have 𝑣1 = 3940.75 and 𝑣2 = 3944.233 with 𝑣1+𝑣2 = 7884.98. For the transformed data 𝑍

print(np.mean(Z**2, axis=0))
print(np.sum(np.mean(Z**2, axis=0)))

[7.88409127e+03 8.90008876e-01]
7884.98127648211

Therefore, the total variability in 𝑍 renains exactly the same as the total variability in 𝑋. However, 99.9% of the total variability is contained in the first column of 𝑍.

np.mean(Z**2, axis=0)/np.sum(np.mean(Z**2, axis=0))

array([9.99887126e-01, 1.12873937e-04])

Principal Component Analysis finds a low-dimensional representation of a data set that contains as much as possible of the variation. The first principal component of a set of features 𝑋1,𝑋2,…,𝑋𝑝 is the normalized linear combination of the features

that has the largest variance. Normalization means

The elements 𝑎11,𝑎21,…,𝑎𝑝1 are known as the loadings of the first principal component.

from sklearn.decomposition import PCA
pca_X = PCA(n_components=2)
prcomps_X= pca_X.fit_transform(X)
pca_X.components_

array([[ 0.71408813,  0.70005581],
       [-0.70005581,  0.71408813]])

Note that, in sklearn, the PCA class object returns the principal components in the rows instead of columns. That is, the first principal component loadings are given by the first row. Observe that, they are very close to our earlier transformation of (1/√2,1/√2).

The transformed data matrix 𝑍, which are known as the prinicipal components scores are given by the object prcomps_X inthe above code. Check that the columns of that matrix are uncorrelated.

np.corrcoef(prcomps_X[:,0], prcomps_X[:,1])[0,1]

-4.1186897937849173e-17

Example with Handwritten Digits Data (MNIST)

Let us now look into an example with a large datset.

The MNIST dataset is a dataset of 70,000 small square 28×28 pixel grayscale images of handwritten single digits between 0 and 9. As each of the image has 28x20 pixels, there are 784 features for each image. Is there any room for data reduction? Can we create simple machine learning algorithms using fewer features?

Let us load the data first.

from sklearn.datasets import fetch_openml
digits = fetch_openml('mnist_784', parser='auto')

digits.data.shape

(70000, 784)

First we split the data into a training and test set with 60000 observations in the training set and 10,000 images in the test set. Then to have better numerical accuracy, divide the pixel values of each image by 255 so that all features are between 0 and 1.

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(digits.data, digits.target, test_size=10000, random_state=42)

X_train = X_train/255.0

Let us try PCA and explore the variance of the principal components.

PCA_digits = PCA()
PCA_digits.fit(X_train)

pcvar = PCA_digits.explained_variance_
plt.plot(range(X_train.shape[1]), pcvar, '.')
plt.xlabel('Principla Components')
plt.ylabel('Explained Variance')
plt.show()

Figure 9: Explained variances of the principal components

We observe that a large amount of variation is explained by only first few principal components. Let us check the proprtions of explained variances for the first 6 components.

# Proportion of Explained Vraiances
print(PCA_digits.explained_variance_ratio_[:6])

#Cumulative proportion of explained variances
print(np.cumsum(PCA_digits.explained_variance_ratio_[:6]))

[0.09736019 0.07162769 0.06157279 0.05407583 0.04894241 0.04314663]
[0.09736019 0.16898788 0.23056067 0.2846365  0.33357892 0.37672555]

So, the first 6 principal components explain 37.67% of the total variation in the data with 784 features.

Let us look into the first two principal components of 2000 randomly selected images with the class information.

sample = np.random.choice(X_train.shape[0], 2000, replace=False)

PC12 = PCA_digits.transform(X_train)[sample,:2]
PClabels = y_train.reset_index(drop=True)[sample]

df = pd.DataFrame({'PC1': PC12[:,0], 'PC2': PC12[:,1], 'digits':PClabels})
sns.scatterplot(data=df, x='PC1', y='PC2', hue='digits')
plt.show()

Figure 10: First two principal components

We can also see the loadings of the first 6 principal components to understand which features are getting more weights.

pccomps = PCA_digits.components_[:6,:]

fig, ax = plt.subplots(2, 3, sharex=True, sharey=True)

for i in range(2):
    for j in range(3):
        icol = i*3+j
        ax[i,j].imshow(np.reshape(pccomps[icol,:], (28,28)), cmap='pink')
        ax[i,j].set_title('PCA'+str(icol+1))

Figure 11: Loadings of the first 6 principal components

We can also look into the principal components with lowest variances.

Figure 12: Loading of the last 6 principal components

From the above plots of the principal component loadings, we observe that the first few principal components with high variances weights very highly the central pixels, which are very important for the idenification of the digits. On the other hand, the last few principal components puts larger weights to the unimportant corner pixels of the images.

Let us train a k-nearest neighbour classifier with 𝑘=3 nearest neighbours using only the first 35 prinicipal components of the training data.

traindata = PCA_digits.transform(X_train)[:,:35]

from sklearn.neighbors import  KNeighborsClassifier

clf = KNeighborsClassifier(3)
clf.fit(traindata, y_train)

X_test = X_test/255.0
testdata = PCA_digits.transform(X_test)[:,:35]

yhat = clf.predict(testdata)

from sklearn.metrics import accuracy_score

print('Accuracy of the test data:',accuracy_score(yhat, y_test))

Accuracy of the test data: 0.9752

The accuracy is more than 97.52% with only 35 principal components.

Hope, you have enjoyed the article!!

For consulting on any data science problems, contact biman.pph@gmail.com