Principal Component Analysis - A Geometric Approach

Overview

Principal Component Analysis is a neat statistical trick, with a simple bit of linear algebra backing it up.  Imagine you've got a dependent, or response variable, $Y$ and a large number of independent, a.k.a. explanatory, variables $X_1, ... X_p$.  You also have $n$ measurements of each, giving you a matrix

$$
\begin{align}
X = &
\begin{pmatrix}
x_{11} & x_{12} & ... & x_{1n}\\
\vdots & & & \\
x_{p1} & x_{p2} & ... & x_{pn}\\
\end{pmatrix} \\
= &
\begin{pmatrix}
-\mathbb{x_1}- \\
\vdots \\
-\mathbb{x_p}- \\
\end{pmatrix} \\
\end{align}
$$

You want to build a model explaining how $Y$ depends on the $X_i$, perhaps a linear model like

$$
Y = \beta_0 + \sum_{i=1}^p \beta_i X_i
$$
but first you need to check that the $X_i$ are more or less independent of each other, otherwise there's no way of uniquely setting the $\beta_i$ values and the stats package is likely to produce an unreliable output.  In general the $X_i$ are not independent, which presents a problem.  This is the problem that PCA seeks to solve, and it does it by rotating the $X_i$ variables in p-space to produce new variables $X_i'$ which are completely independent of eachother.  In the process it provides a mechanism for reducing the number of explanatory variables, enabling a more parsimonious model.

Variance

From here on let's assume that the $X_i$ have been scaled and centered as it will make everything that follows simpler.  We can think of each observation as being a geometric vector in $p$ dimensional space.  These vectors can be squared and summed to give a variance $\sigma^2$ i.e.
$$
\sigma^2 = \frac{1}{n}\sum_{j=1}^n{\left|
\begin{pmatrix}
x_{1j} \\
\vdots \\
x_{pj}\\
\end{pmatrix}
\right|^2
}
$$
This can be rewritten using the per variable variances to get
$$
\begin{align}
\sigma^2 &= \frac{1}{n}\sum_{i,j}{x_{ij}^2} \\
&= \sigma_1^2 + ... + \sigma_p^2 \\
&= p
\end{align}
$$

Covariance

Since the $X_i$ are centered, the requirement that distinct variables $X_i$ and $X_j$ are independent is the same as saying that the covariance matrix $C = \frac{1}{n}XX^{T}$ is diagonal.  In general this is not the case, but perhaps we can find linear combinations $X_i'$ of the original $X_i$ such that it is!

Rotation

We know that $C$ is symmetric and real, and so it follows from basic linear algebra that it has $p$ linearly independent and orthonormal eigenvectors $\mathbb{e_1}, ...,\mathbb{e_p}$ corresponding to real eigenvalues $\lambda_1, ... \lambda_p$.  I.e. we can find column vectors $\mathbb{e_i}$ and scalars $\lambda_i$ such that
$$
C \mathbb{e_i} = \lambda_i \mathbb{e_i}
$$
from which it follows that whenever $i \ne j$
$$
\mathbb{e_j}^TXX^T \mathbb{e_i} = 0
$$
But the LHS above is just the dot product of $\mathbb{x_j}'$ and $\mathbb{x_i}'$ defined by the rule
$$
\mathbb{x_i}' = \mathbb{e_i}^T X
$$
Hey presto! We've defined a rotation of the $X_i$ which results in what we wanted: independent explanatory variables.  (I'm calling it a rotation because just like with rotations in 3D space, the matrix defined by the $\mathbb{e_i}$ is inverse to its transpose.)  These transformed variables are referred to as Principal Components.

Parsimony

Typically we want to build a model with as few independent variables as possible. If a model with lots of explanatory variables explains most of the variation in the response variable, it could well be just that there's enough degrees of freedom to explain any response!  So it's much better to get rid of a few variables, as long as those remaining can still explain most of the variation.  With the original $X_i$ there's no obvious way to decide which to ditch but, as we shall see, with the transformed $X_i'$ variables the eigenvalues can help us find the subset that is most "important".
 
The original variance $\sigma^2$ is just the average of the squares of $n$ p-vectors.  Since our transformation merely rotates these p-vectors their sizes do not change, and consequently
$$
\sigma'^2 = \sigma^2 = p
$$
It turns out that the $\lambda_i$ are the same as the per-principal component variances since
$$
\begin{align}
\sigma_i'^2 &= \frac{1}{n}\mathbb{x_i'}\mathbb{x_i'}^T \\
&= \frac{1}{n} \mathbb{e_i}^T X X^T\mathbb{e_i} \\
&=\mathbb{e_i}^T C\mathbb{e_i} \\
&=\lambda_i \\
\end{align}
$$
This tells us that the eigenvalues (variances) must all be positive and sum to $p$. Remember, we're looking at exactly the same data, just rotated in p-space.  Thus, if any eigenvalues happen to be near zero, this says that we can ditch the corresponding principal component.  In general we arrange the eigenvalues in decreasing order and name the corresponding principal components PC1, PC2, and so on, and then we ditch the ones at the end that contribute less than 5% of the variance between them.  This is equivalent to projecting the observations onto a hypersurface in such a way as to lose the least amount of information.  Very often just two principal components suffice to cover most of the variance enabling the creation of a simple scatterplot which can help identify modal groups.

Example

The following short R script shows PCA at work for an imagined experiment.  Photons are fired from a torch and the $x,y,z$ distances travelled are measured, along with the time taken $t$.  Obviously $t$ depends wholly on the other components, so let's see if PCA spots that we can reduce these 4 variables to just 3.

# Print real numbers, not scientific notation.
options(scipen=999)

# Produce 1000 x,y,z,t displacement observations for photons emitted by a torch
n <- 1000
x <- rnorm(n, mean=10, sd=1)
y <- rnorm(n, mean=8,  sd=0.5)
z <- rnorm(n, mean=12, sd=1.5)
t <- (x^2 + y^2 + z^2)^(1/2)

# Combine into dataframe
df <- data.frame(x=x, y=y, z=z, t=t)

# Show how t is not independent of x,y,z
library(car)
scatterplotMatrix(df)

# Do a PCA
pca <- prcomp(df, scale=T, center=T)
summary(pca)
print(pca$rotation)

The dependency between the original variables can be seen by the scatterplot matrix the script generates:

The text output shows that the PCA does in fact detect this

> summary(pca)
Importance of components:
                          PC1    PC2    PC3     PC4
Standard deviation     1.3934 1.0327 0.9950 0.04177
Proportion of Variance 0.4854 0.2666 0.2475 0.00044
Cumulative Proportion  0.4854 0.7520 0.9996 1.00000

The above output shows that half the variation is covered by PC1 on it's own (corresponding to an eigenvalue of 2 where $p =4$). This is because R has noticed that if you scale each of $x, y$, and $z$ and add them together you get a first order approximation for $t$ as well.  The final PC adds almost nothing to the variance, as expected, and can be dropped.

Finally, R lets you see exactly how each PC is defined in terms of the original variables, and you can check that these form an orthonormal set.

> print(pca$rotation)
        PC1         PC2         PC3        PC4
x 0.2746490  0.89112242 -0.08053689 -0.3521116
y 0.1857603 -0.02987745  0.97020834 -0.1526310
z 0.6142378 -0.44980563 -0.22700113 -0.6073362
t 0.7160817  0.05179825 -0.02607824  0.6956033

Comments

Post a Comment

Popular posts from this blog

How To Make ASCII Diagrams Beautifuller

Image
I've discovered an excellent tool in asciiflow.com .  The website makes it really easy to create ASCII box diagrams like this They put these things in fruit machines you know! This is ideal for source code banners, which I think should contain helpful documentation - but most programmers think it's a good place for the COPYRIGHT information and nothing else. But wait! we can make it beautifuller... and easier to read... by replacing some of the ASCII characters with ones available in UTF-8: There!  Isn't that better?  (Although some purists may object to non-ASCII characters in your code base.) SOURCE CODE: #!/usr/bin/python2 # coding: utf-8 # + gets converted in different ways depending on it's 4 neighbours # # . N . { nsew(N,S,E,W) has bit 3 set if N in "+|<>" # W + E { nsew(N,S,E,W) has bit 2 set if S in "+|<>" # . S . { nsew(N,S,E,W) has bit 1 set if W in "+-^v" # { nsew(N,S,E,W) has bit 0 set if E ...
Read more

Why growth is falling in all developed countries (as a long term trend)

Image
In this post I argue that economic growth is not the normal state of affairs, it is a blip caused by the discovery of a one-off resource of non-renewable energy.   The logistic equation $$ \frac{dy}{dt} = y(1-y) $$ Here's a simple model.  Humans discover an exploitable but limited resource and start consuming it.  The amount consumed, $y(t)$, is a function of time.  The general form of the equation is $\frac{dy}{dt} = \alpha y (\beta - y)$, but if you choose the right units $\alpha$ and $\beta$ both become $1$. Why should this work, in principle? Early on, the factor $1-y$ is approximately $1$ and can be ignored.  So the model states that annual consumption $dy/dt$ starts off proportional to $y$.  In other words, $y$ grows exponentially at first.  This could happen if exploiting the resource enables further exploitation of the resource.  For example, suppose a few humans are shipwrecked on an island with 1000 trees: they take ages to cut down the ...
Read more

Three ways to look at the Bell/GHZ experiment

Image
In 1964 John Stewart Bell proposed an experiment to determine whether the results of quantum measurements were truly random, or governed by hidden variables, i.e. state that exists prior to the measurement, but which we don’t have access to. The experiment involved creating a large number of EPR pairs, and firing them at two observers, Alice and Bob, who measure their photon’s polarisation, choosing the $\updownarrow$ direction or the $\nearrow\llap\swarrow$ direction at random. Determining the result of the experiment involves doing a complex statistical calculation to see if something called Bell’s inequality is satisfied or violated. The Bell experiment was first performed by in 1982 by Alain Aspect, and the result, as most commonly interpreted, is that hidden variables can only exist if Quantum Mechanics is non-local, i.e. if it supports faster-than-light causality! Some time after Bell proposed his experiment, Greenberger, Horne, and Zeilinger suggested an alternativ...
Read more