Looking at PCA in Action

There are different approaches to perform a PCA operation. Software tools (such as Matlab and R) offer two common approaches that we are going to discuss here.

1. Using Eigen-decomposition

The first approach is based on what was written in the previous blog post on PCA. The step-by-step example of PCA below is in Matlab as it is easier to work with matrices. The first step is to mean-center the original data. Suppose you have 2D matrix called data, with [x₁, x₂, ..., x₁₀₀] and [y₁, y₂, ..., y₁₀₀] with 100 observations.

% Step-1: Mean-center the original data.
data2(:,1) = data(:,1) - mean(data(:,1));
data2(:,2) = data(:,2) - mean(data(:,2));

Next, we find the covariance matrix to capture all possible relationships of the variables in the dataset.

% Step-2: Calculate covariance matrix of the mean-centered data.
C = cov(data2);

Once the covariance matrix C is computed we then perform eigen-decomposition on it to find the eigenvectors and eigenvalues according to CV = VD, leading to C = VDV^T. (I think we have to flip the matrix V left and right).

% Step-3: Produce eigenvalues and eigenvectors: a diagonal matrix D of eigenvalues and matrix V
% whose columns are the corresponding eigenvectors so that C*V = V*D.
[V,D] = eig(C);

Remember now the eigenvalues measure the "amount of variation" retained by each principal component. So, we rank the eigenvalues according to the magnitude or "importance". For example, suppose the columns of matrix V and the diagonal elements of D contain the following (see right):

The diagonal elements contain two eigenvalues λ₁ and λ₂. After ranking these two, it is seen that  λ₁ > λ₂, which means that the eigenvector v₁ corresponds to the first component, and the one that corresponds to the second is v₂.
Finally, the last step is to rotate or re-orient the data by [new data] = [original data].[eigenvector]^T. The new dataset is called the component scores.

% Step-4: The new axes is formed using V, to rotate the mean-centered data.
newdata = data2 * V';
newdata = fliplr(newdata) 

What is the percentage of information explained or accounted for by each component? It turns out that this can be found by dividing the corresponding eigenvalue by the total sum of all eigenvalues.

% Step-5: Get the variance accounted for, use matrix D containing the eigenvalues.
explained = D/sum(D(:))

Look at Fig-1. The original dataset in the left panel shows correlated variables. The middle panel shows how the dataset can be represented in terms of two principal components. PC1 axis is the first principal direction along which the data show the biggest variation. PC2 axis is the second most important direction and it is orthogonal to the PC1 axis. The new dataset in the right panel is a "projection" onto the principal components, where projection onto PC1 (X-axis) is uncorrelated with projection onto PC2 (Y-axis). 

Fig-1: The process of rotating the original data into a new dimension of principal components.

Where is the dimensionality reduction here? Take a look at the right panel above. The transformed data have just as many dimensions as the original data. By visual inspection, most variability is along the PC1 (X-axis) so you can just throw away the data plotted on the Y-axis. In other words, our 2D data can be reduced to 1D by projecting each sample onto the first principal component. Do we lose anything? Well, a little bit, depending on the variation that the PC1 and PC2 explain. In data science, data variation means information. Based on the above example, we find that PC1 and PC2 carry 1.284/(0.049+1.284) = 96%; and 0.049/(0.049+1.284) = 4% of the variance of the data respectively.

2. Using Singular Value Decomposition
PCA is intimately related to another statistical analysis called the Singular Value Decomposition (SVD). Although notations may differ in different textbooks, SVD essentially involves the factorization of matrix X = UDV^T, where U is an n × n matrix, with columns as orthogonal unit vectors of length n (left singular vector of X); D is an n × p rectangular diagonal matrix, the singular values of X; V is a p × p, with columns as orthogonal unit vectors of length p (right singular vectors of X).

The new data can be obtained through Y = XV = UD.  The transpose of V is sometimes called the whitening and can be used as a preparation in ICA. Columns of V multiplied by the square root of corresponding eigenvalues, i.e. eigenvectors scaled up by the variances, are the coefficients or loadings.

Now PCA can also be done in Matlab directly using a function:
[coeff,newdata,latent,tsquared,explained,mu] =  pca(data)

This Matlab function automatically mean-centers the data. The coefficients are also known as loadings. The loadings can be understood as the weights for each original variable when calculating the new scores. Each column of coeff contains coefficients for each component, and the columns are in descending order of the variance. If you type coeff and V in the Matlab terminal, you may find similarity. Recall that the eigenvectors are like the coefficients to predict the new scores from the original data (see Step-4 above!). But eigenvectors scaled up by the variances are called loadings. The terms can be confusing because the Matlab command already does normalization by default.

The newdata refers to the scores, which is a 2 × 100 matrix that represents the transformed data in the principal component space. The maximum number of principal components equals the number of variables in X. The component variances or latents are actually the eigenvalues of the covariance matrix of X. The latent values are arranged in descending order according to the magnitude of variance, i.e. ranking the eigenvalues. The function also returns explained, i.e. the proportion of the total variance explained by each component; and mu, the estimated mean of each variable in X.

To recover the original data, simply multiply the component scores with the transpose of matrix V. Remember to put back the mean value to each of the variables.

% Recover the original data
origdata = fliplr(newdata) * V' + mu;

In R, PCA can be performed automatically using a function princomp (using eigen-decomposition) or prcomp (using the SVD operation, for better accuracy in the results) from the built-in stats package. We will discuss one example in R below.

3. Case-study Using PCA
One nice package to easily extract and visualize the results of PCA is by using functions from the "factoextra" package that is based on "ggplot2". As the demo dataset, I am using the decathlon2 dataset from the "factoextra". Briefly, it contains a list of 27 athletes who have performed several different types of sports such as javelin throw, long jump, and so on. We call the athletes as individuals and the sports performance as the numeric variables. The dataset can be divided into:
- Active individuals (row. 1 to 23).
- Active variables (col. 1 to 10), to be used in PCA.
- Supplementary individuals (row. 24 to 27).
- Supplementary variables (col. 11 to 13), used in later prediction.

We first perform the PCA on the Active individuals/variables with the scaling option. This is an important step when you have values of different scales, e.g. metres and feet. After that, we can use the spree plot to visualize the different principal components ranked according to the proportion of variance explained. Recall that the eigenvalues reflect the amount of variance retained by the components. Eigenvalues are large for the first PC and small for the subsequent PCs. The first PC corresponds to a direction with the maximum amount of variation in the data.

We can ask other questions arising from the PCA results:

(1) How do we know the variance accounted after PCA?
Get this by extracting using the summary() function, we can see that the PC1 accounts for 41% of the total variance in the data, PC2 accounts for 18% variance, and so on. This is also what fviz_eig() provides. There are up to 10 original variables so we get up to 10 principal components.

(2) How many principal components are useful?
PCA does not tell us # components are useful, so there is no magic answer for this. From the cumulative proportion as shown above (or the spree plot), we can safely say that the first six components (PC1 - PC6) are indeed sufficient as they account for more than 90% of the total variance in the data.

(3) How are the original axes look like in terms of the components?
Before PCA, we have in the dataset the original axes which are formed by the 10 variables. After PCA, we get rotated axes which are formed by the principal components. To extract the rotated variables, we write get_pca_var(res.pca). Note that Dim.xx and PCxx are the same thing!

We can then plot how the original variables look like in terms of the first two components, PC1 and PC2. To visualize can be achieved using fviz_pca_var(), which will output the graph as shown below, Fig-2a on the left panel. The plot shows a correlation plot or relationships between all variables based on the first two PCs. Negatively correlated variables are positioned on opposite sides of the quadrants. The distance between variables and the origin measures the quality of the variables on the new map. Variables that are away from the origin are well represented on the map.

Fig-2: (a) Left: A correlation circle with a radius of 1.0 shows how the original variables relate to one another in the new dimensions of the principal components. (b) Right: a graphical representation of how each variable contributes to the principal component. It can also be seen that the representations of the original variables to PC7 to PC10 are minimal.

(4) What is the representation of each variable to the principal components?

We can use the function from the "corrplot" package to highlight the most contributing variables for each component. The function takes in var$cos2 as the argument. This cos2 shows the quality of representation of the variables on the component map (square cosine, squared coordinates).

Refer to Fig-2b. There are 2 ways to interpret this figure. First, it indicates a degree of representation of the original variables (the rows) on each principal component (the columns). Both panels in Fig-2 are somehow related. If a variable has a high cos2 value, it will be seen as a large and dark blue circle in Fig-2b, the same variable is positioned closer to the circumference of the plot in Fig-2a. A small circle indicates that the variable is not perfectly represented by the PCs. In this case, the variable is close to the center of the circle. Another way is to say how each principal component contributes to the variable. So, for a given variable, the sum of the cos2 on all the principal components is equal to one.


For further reading, see below:
(1) PCA in 6 steps
(2) Tutorial on PCA
(3) Making sense of eigenvector/eigenvalue

Basic Math for Principal Component Analysis

To begin, Principle Component Analysis (PCA) is a statistical tool for exploratory analysis that is famous for reducing the dimensionality of the data. It is a technique that involves vector space transformation to represent data in another form. It is useful for feature extraction and feature elimination, removing redundant features. The outcome of a PCA operation is a set of orthogonal variables, called the principal components. To learn PCA from a technical angle, I recommend the work by J Shlens (ref.[1]).

1. Intro to Matrix Operations
Before going into the math, there are a few technical terms in matrix algebra we have to be familiar with. First, eigenvector and eigenvalue. An eigenvector v is of a linear transformation (or simply a square matrix A) is a non-zero vector that changes only by a scalar factor (i.e. by the eigenvalue, λ) when that transformation is applied to it. We can say, [matrix].[eigenvector] = [eigenvalue].[eigenvector]; or we can write: Av = λv.  If the eigenvalue is negative, the direction is reversed. Refer to the figure on the right for a geometric description.

An eigenvector, with a real nonzero eigenvalue, points in a direction that is stretched by the transformation and the eigenvalue is the factor by which it is stretched.

Another term is about an orthogonal matrix. Let a square matrix A is invertible (it can be inversed) or non-singular.

A square matrix is orthogonal with real entries if, and only if, its transpose (A^T) is equal to its inverse (A^-1), that is AA^T yields an identity matrix I.  The square matrix is also normal if it satisfies A^TA = AA^T.

Importantly, if the matrix is orthogonal then it is necessarily invertible, unitary, and normal. The determinant of an orthogonal matrix is either +1 or −1. Next, we learn about another special matrix called a diagonal matrix, i.e. a square matrix with non-zero diagonal elements and all-zero remaining (non-diagonal elements). When a square matrix A is diagonalizable, there exists an invertible matrix B (a matrix that can have an inverse) such that B^-1AB produces a diagonal matrix.

Matrix S is symmetric only if it has a characteristic of S^T = S. Such matrix is also orthogonally (orthonormally) diagonalisable.

An example of a square matrix related to PCA is called the covariance matrix. More specifically, in this matrix we have the variances on the diagonal entries of this matrix, and the covariances on the off-diagonal entries. Take a look below. Here, the red box denotes the variance of variable height, whereas the green box shows the covariance between the score and age.

Suppose we have a huge dataset which has too much information. Can we re-express the original dataset optimally in a new form? Now the word 'optimality' means we do not want correlated or redundant features. After a series of matrix operations, the new or transformed dataset will have a new coordinate system. The whole idea of PCA is to force as much of the variation as possible into fewer dimensions, you can throw away the rest without losing much information. Let's relate PCA with the previous ideas on matrices:

In PCA, each eigenvector is a unit vector pointing in the direction of a new coordinate axis. The axis with the highest eigenvalue is said to be the axis that explains the largest variation, meaning an eigenvalue represents variance. The directions with the largest variances are the most “important” or principal. 

2. Technical Definition

Suppose we have a dataset of m × n matrix X, where the n columns are the number of observations and the m rows are the variables. We wish to linearly transform this matrix X into another matrix Y, also of dimension m × n, so that for some m × m transformation matrix P, we will obtain Y = PX.

So, PCA attempts to express correlated m variables into uncorrelated new variables through some sort of transformation. In matrix algebra, all possible correlations of any variable-pairs can be captured using the covariance matrix C_X. How do we decouple the variables? It turns out, we can get this by maximizing the variances or diagonal elements, and minimizing the co-variances or off-diagonal entries of the covariance matrix. So our goal is to obtain a new data Y such that the covariance matrix of the transformed data C_Y is strictly diagonal. To put it in another way, if we can find the transformation matrix P in such a way that C_Y is diagonal, then our objective is met.

The solution is based on an important property of matrix algebra on C_X called the eigen-decomposition. The effort to make C_Y a diagonal matrix is essentially to find P = E^T. How to get this new matrix? It turns out that the orthonormal matrix E can be found by diagonalizing the covariance matrix of the original dataset, C_X = EDE^T , where E is the matrix of eigenvectors which applies the transformation on X; and matrix D has diagonal elements containing the eigenvalues in descending order. In short, the principal components are formed by:
- The “new directions” of our data as indicated by the eigenvectors;
- The “magnitude”, or importance of each direction, as indicated by the eigenvalues.

Often, the normalized eigenvectors by the variance are also called the coefficients or loadings of the principal components. Remember, these coefficients define how the new principal component is formed as a linear combination of the original variables. Each column in the eigenvector matrix contains coefficients for each component in the descending order of the variance.

Dimensionality reduction works as follows. Suppose our original data is an m × n matrix X with m dimension and n data points.

If we want to transform the points to k dimensions (where, k < m), then select the first k eigenvectors of the matrix C_X sorted decreasingly according to the eigenvalues, and form a matrix with them, then use them as k × m matrix P. The resulting matrix is Y with size (k × m)(m × n) = k × n.

3. Closing remarks

PCA does not simply do a rotation. Normalization is always involved prior to PCA, which can be done by mean-centering the data and divide them with the standard deviation. This is important if our data contain variables of different scales, e.g. kilogram and pounds.

PCA and ICA are two slightly different things. While PCA expresses the data into a set of orthogonal variables, ICA simply represents data into an independent set of variables (not necessarily orthogonal). In ICA, the order of importance is irrelevant. Prior to ICA, data whitening is performed and PCA is one method to whiten the data. The figure on the right is common on the internet. PCA is predominantly used as a dimensionality reduction technique, for example, in computer vision, facial recognition, and image compression. It is also used for finding patterns in a high dimensional dataset in the field of finance, data mining, biomechanics, bioinformatics, etc.

References:
(1) "A tutorial on Principal Components Analysis" by J. Shlens (2014). Very nice technical tutorial!
(2) https://towardsdatascience.com/a-step-by-step-explanation-of-principal-component-analysis-b836fb9c97e2
(3) https://stats.stackexchange.com/questions/612/is-pca-followed-by-a-rotation-such-as-varimax-still-pca