Let's use the Iris dataset to understand how to use PCA efficiently in reducing the dimension of the dataset. The Iris dataset contains measurements for 150 iris flowers from three different species.

The three classes in the Iris dataset are as follows:

The following are the four features in the Iris dataset:

Can we use, say, two columns instead of all the four to express most of the variations in the data? Our quest is to reduce the dimension of the data. In this case, our instances have four columns. Let's say that we are building a classifier to predict the type of flower with a new instance; can we do this task using instances in the reduced dimension space? Can we reduce the number of columns from four to two and still achieve a good accuracy for our classifier?

PCA is done using the following steps:

Let's load the necessary libraries and call the utility function load_iris from scikit-learn to get the Iris dataset:

import numpy as np
from sklearn.datasets import load_iris
from sklearn.preprocessing import scale
import scipy
import matplotlib.pyplot as plt

# Load Iris data
data = load_iris()
x = data['data']
y = data['target']

# Since PCA is an unsupervised method, we will not be using the target variable y
# scale the data such that mean = 0 and standard deviation = 1
x_s = scale(x,with_mean=True,with_std=True,axis=0)

# Calculate correlation matrix
x_c = np.corrcoef(x_s.T)

# Find eigen value and eigen vector from correlation matrix
eig_val,r_eig_vec = scipy.linalg.eig(x_c)
print 'Eigen values 
%s'%(eig_val)
print '
 Eigen vectors 
%s'%(r_eig_vec)


# Select the first two eigen vectors.
w = r_eig_vec[:,0:2]

# # Project the dataset in to the dimension
# from 4 dimension to 2 using the right eignen vector
x_rd = x_s.dot(w)

# Scatter plot the new two dimensions
plt.figure(1)
plt.scatter(x_rd[:,0],x_rd[:,1],c=y)
plt.xlabel("Component 1")
plt.ylabel("Component 2")

Now, we will proceed to Standardize this data, with a zero mean and standard deviation of one, we will leverage the numpyscorr_coef function to find the correlation matrix:

x_s = scale(x,with_mean=True,with_std=True,axis=0)
x_c = np.corrcoef(x_s.T)

We will then do the Eigenvalue decomposition and project our Iris data on the first two principal Eigenvectors. Finally, we will plot the dataset in the reduced space:

eig_val,r_eig_vec = scipy.linalg.eig(x_c)
print 'Eigen values 
%s'%(eig_val)
print '
 Eigen vectors 
%s'%(r_eig_vec)
# Select the first two eigen vectors.
w = r_eig_vec[:,0:2]

# # Project the dataset in to the dimension
# from 4 dimension to 2 using the right eignen vector
x_rd = x_s.dot(w)

# Scatter plot the new two dimensions
plt.figure(1)
plt.scatter(x_rd[:,0],x_rd[:,1],c=y)
plt.xlabel("Component 1")
plt.ylabel("Component 2")

Using function scale. The scale function can perform centering, scaling and standardization. Centering is subtracting the mean value from individual values, Scaling is dividing each value by the variable's standard deviation and finally Standardization is performing centering followed by scaling. Using variables with_mean and with_std function scale can be used to perform all three normalization techniques.

The Iris dataset has four columns. Though there are not many columns, it will serve our purpose. We intend to reduce the dimensionality of the Iris dataset to two from four and still retain all the information about the data.

We will load the Iris data to the x and y variables using the convenient load_iris function from scikit-learn. The x variable is our data matrix and we can inspect its shape as follows:

>>>x.shape
(150, 4)
>>>

We will scale the data matrix x to have zero mean and unit standard deviation. The rule of thumb is that if all your columns are measured in the same scale in your data and have the same unit of measurement, you don't have to scale the data. This will allow PCA to capture these basic units with the maximum variation:

x_s = scale(x,with_mean=True,with_std=True,axis=0)

We will proceed to build the correlation matrix of our input data:

The correlation matrix of n random variables X1, ..., Xn is then × n matrix whosei, jentry is corr (Xi, Xj), Wikipedia.

We will then use the SciPy library to calculate the Eigenvalues and Eigenvectors of the matrix.Let's look at our Eigenvalues and Eigenvectors:

print Eigen values 
%s%(eig_val)
print 
 Eigen vectors 
%s%(r_eig_vec)

The output looks as follows:

In our case, the Eigenvalues are printed in a descending order. A key question is how many components should we choose? In the next section, we will explain a few ways of choosing the number of components.

You can see that we selected only the first two columns of our right-hand side Eigenvectors. The discrimination capability of the retained components on the y variable is a good test of how much information or variation is retained in the data.

We will project the data to the new reduced dimension.

Finally, we will plot the components in the x and y axes and color them by the target variable:

You can see that components 1 and 2 are able to discriminate the three classes of the iris flowers. Thus we have effectively used PCA in reducing the dimension to two from four and still able to discriminate the instances belonging to different classes of Iris flower.

In the previous section, we said that we would outline a couple of ways to help us select how many components should we include. In our recipe, we included only two. The following are a list of ways to select the components more empirically:

For each component, we printed the Eigenvalue, percentage of the variance explained by that component, and cumulative percentage value of the variance explained. For example, component 1 has an Eigenvalue of 2.91; 2.91/4 gives the percentage of the variance explained, which is 72.80%. Now, if we include the first two components, then we can explain 95.80% of the variance in the data.

The decomposition of a correlation matrix into its Eigenvectors and values is a general technique that can be applied to any matrix. In this case, we will apply it to a correlation matrix in order to understand the principal axes of data distribution, that is, axes through which the maximum variation in the data is observed.

PCA can be used either as an exploratory technique or as a data preparation technique for a downstream algorithm. Document classification dataset problems typically have very large dimensional feature vectors. PCA can be used to reduce the dimension of the dataset in order to include only the most relevant features before feeding the data to a classification algorithm.

A drawback of PCA worth mentioning here is that it is computationally expensive operation. Finally a point about numpy's corrcoeff function. The corrcoeff function will standardize your data internally as a part of its calculation. But since we want to explicitly state the reason for scaling, we have included it in our recipe.