Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

5.4. Principal Components Analysis

Introduction

As we saw in Chapter 5.3, the singular value decomposition of a matrix XX into X=UΣVTX = U \Sigma V^T can be used to find the best low-rank approximation of XX, which has applications in image compression, for example. But, hiding in the singular value decomposition is another, more data-driven application: dimensionality reduction.

Consider, for example, the dataset below. Each row corresponds to a penguin.

Loading...

Each penguin has three categorical features – species, island, and sex – and four numerical features – bill_length_mm, bill_depth_mm, flipper_length_mm, and body_mass_g. If we restrict our attention to these numerical features, we have a 333×4333 \times 4 matrix, which we can think of as 333 points in R4\mathbb{R}^4. We, of course, can’t visualize in R4\mathbb{R}^4. One solution is to draw scatter plots of all pairs of numerical features. Points are colored by species.

Image produced in Jupyter

While we have four features (or “independent” variables) for each penguin, it seems that some of the features are strongly correlated, in the sense that the correlation coefficient

r=1n∑i=1n(xi−xˉσx)(yi−yˉσy)r = \frac{1}{n} \sum_{i=1}^n \left( \frac{x_i - \bar{x}}{\sigma_x} \right) \left( \frac{y_i - \bar{y}}{\sigma_y} \right)

between them is close to ±1\pm 1. (To be clear, the coloring of the points is by species, which does not impact the correlations between features.)

penguins[['bill_length_mm', 'bill_depth_mm', 'flipper_length_mm', 'body_mass_g']].corr()
Loading...

For instance, the correlation between flipper_length_mm and body_mass_g is about 0.87, which suggests that penguins with longer flippers tend to be heavier. Not only are correlated features redundant, but they can cause interpretability issues – that is, multicollinearity – when using them in a linear regression model, since the resulting optimal parameters w⃗∗\vec w^* will no longer be interpretable. (Remember, in h(x⃗i)=w⃗⋅Aug(x⃗i)h(\vec x_i) = \vec w \cdot \text{Aug}(\vec x_i), w⃗∗\vec w^* contains the optimal coefficients for each feature, which measure rates of change when the other features are held constant. But if two features are correlated, one can’t be held constant without also changing the other.)

Here’s the key idea: instead of having to work with d=4d = 4 features, perhaps we can instead construct some number k<dk < d of new features that:

  1. capture the “most important” information in the dataset,

  2. while being uncorrelated (i.e. linearly independent) from one another.

I’m not suggesting that we just select some of the original features and drop the others; rather, I’m proposing that we find 1 or 2 new features that are linear combinations of the original features. If everything works out correctly, we may be able to reduce the number of features we need to deal with, without losing too much information, and without the interpretability issues that come with multicollinearity.

To illustrate what I mean by constructing a new feature, let’s suppose we’re starting with a 2-dimensional dataset, which we’d like to reduce to 1 dimension. Consider the scatter plot below of flipper_length_mm vs. body_mass_g (with colors removed).

Image produced in Jupyter

Above, each penguin is represented by two numbers, which we can think of as a vector

x⃗i=[body massiflipper lengthi]\vec x_i = \begin{bmatrix} \text{body mass}_i \\ \text{flipper length}_i \end{bmatrix}

The fact that each point is a vector is not immediately intuitive given the plot above: I want you to imagine drawing arrows from the origin to each point.

Now, suppose these vectors x⃗i\vec x_i are in the rows of a matrix XX, i.e.

X=[body mass1flipper length1body mass2flipper length2⋮⋮body massnflipper lengthn]X = \begin{bmatrix} \text{body mass}_1 & \text{flipper length}_1 \\ \text{body mass}_2 & \text{flipper length}_2 \\ \vdots & \vdots \\ \text{body mass}_n & \text{flipper length}_n \end{bmatrix}

Our goal is to construct a single new feature, called a principal component, by taking a linear combination of both existing features. Principal component analysis, or PCA, is the process of finding the best principal components (new features) for a dataset.

new featurei=a⋅body massi+b⋅flipper lengthi\text{new feature}_i = a \cdot \text{body mass}_i + b \cdot \text{flipper length}_i

This is equivalent to approximating this 2-dimensional dataset with a single line. (The equation above is a re-arranged version of y=mx+by = mx + b, where xx is body massi\text{body mass}_i and yy is flipper lengthi\text{flipper length}_i.)

What is the best linear combination of the features, or equivalently, what is the best line to approximate this 2-dimensional dataset? To make this idea precise, we need to define what we mean by “best”. It may be tempting to think that we should find the best-fitting line that predicts flipper_length_mm from body_mass_g, but that’s not quite what we want, since this is not a prediction problem: it’s a representation problem. Regression lines in the way we’ve seen them before are fit by minimizing mean squared vertical error, yi−h(xi)y_i - h(x_i), but don’t at all capture horizontal errors, i.e. errors in the body_mass_g direction here.

So instead, we’d like to find the line such that, for our nn vectors (points) x⃗1,…,x⃗n\vec x_1, \ldots, \vec x_n, the mean squared distance to the closest vector (point) on the line is minimized. Crucially, this distance is calculated by computing the norm of the difference between the two vectors. If p⃗i\vec p_i is the closest vector (point) on the line to x⃗i\vec x_i, then we’d like to minimize

1n∑i=1n∥x⃗i−p⃗i∥2\frac{1}{n} \sum_{i=1}^n \lVert \vec x_i - \vec p_i \rVert^2

But, as we saw in Chapter 2.3, the vector (point) on a line that is closest to another vector (point) x⃗i\vec x_i is the orthogonal projection of x⃗i\vec x_i onto that line. The resulting error vectors are orthogonal to the line we’re projecting on. As such, I’ll say that we’re looking for the line that minimizes the mean squared orthogonal projection error.

Image produced in Jupyter

The difference between the two lines above is key. The top line reflects the supervised nature of linear regression, where the goal is to make accurate predictions. The bottom line reflects the fact that we’re looking for a new representation of the data that captures the most important information. These, in general, are not the same line, since they’re chosen to minimize different objective functions.

Recall from Chapter 1 that unsupervised learning is the process of finding structure in unlabeled data, without any specific target variable. This is precisely the problem we’re trying to solve here. To further distinguish our current task from the problem of finding regression lines, I’ll now refer to finding the line with minimal orthogonal projection error as finding the best direction.

Minimizing Orthogonal Error

The first section above told us that our goal is to find the best direction, which is the direction that minimizes the mean squared orthogonal error of projections onto it. How do we find this direction?

Fortunately, we spent all of Chapter 2.3 and Chapter 2.10 discussing orthogonal projections. But there, we discussed how to project one vector onto the span of one or more other given vectors (e.g. projecting the observation vector y⃗\vec y onto the span of the columns of a design matrix XX). That’s not quite we’re doing here, because here, we don’t know which vector we’re projecting onto: we’re trying to find the best vector to project our n=333n = 333 other vectors onto!

Suppose that in general, XX is an n×dn \times d matrix, whose rows x⃗i\vec x_i are data points in Rd\mathbb{R}^d. (In our current example, d=2d = 2, but everything we’re studying applies in higher dimensions.) Then,

  • The orthogonal projection of x⃗i\vec x_i onto v⃗\vec v is

    p⃗i=x⃗i⋅v⃗v⃗⋅v⃗v⃗\vec p_i = \frac{\vec x_i \cdot \vec v}{\vec v \cdot \vec v} \vec v

  • To keep the algebra slightly simpler, if we assume that v⃗\vec v is a unit vector (meaning ∥v⃗∥=v⃗⋅v⃗=1\lVert \vec v \rVert = \sqrt{\vec v \cdot \vec v} = 1), then the orthogonal projection of x⃗i\vec x_i onto v⃗\vec v is

    p⃗i=(x⃗i⋅v⃗)v⃗\vec p_i = (\vec x_i \cdot \vec v) \vec v

    This is a reasonable choice to make, since unit vectors are typically used to describe directions.

  • If v⃗∈Rd\vec v \in \mathbb{R}^d, then span({v⃗})\text{span}(\{\vec v\}) is a line through the origin in Rd\mathbb{R}^d.

Centering the Data

The issue is that our dataset – and most datasets in general, by default – aren’t centered at the origin, but the projections p⃗i=(x⃗i⋅v⃗)v⃗\vec p_i = (\vec x_i \cdot \vec v) \vec v all lie on the line through the origin spanned by v⃗\vec v. So, trying to approximate our dataset using a line through the origin is not a good idea.

Image produced in Jupyter
Image produced in Jupyter

So, the solution is to first center the data by subtracting the mean of each column from each entry in that column. This gives us a new matrix X~\tilde X.

X=[x1(1)x1(2)⋯x1(d)x2(1)x2(2)⋯x2(d)⋮⋮⋱⋮xn(1)xn(2)⋯xn(d)]→X~=[x1(1)−μ1x1(2)−μ2⋯x1(d)−μdx2(1)−μ1x2(2)−μ2⋯x2(d)−μd⋮⋮⋱⋮xn(1)−μ1xn(2)−μ2⋯xn(d)−μd]X = \begin{bmatrix} x_1^{(1)} & x_1^{(2)} & \cdots & x_1^{(d)} \\ x_2^{(1)} & x_2^{(2)} & \cdots & x_2^{(d)} \\ \vdots & \vdots & \ddots & \vdots \\ x_n^{(1)} & x_n^{(2)} & \cdots & x_n^{(d)} \end{bmatrix} \rightarrow \tilde X = \begin{bmatrix} x_1^{(1)} - \mu_1 & x_1^{(2)} - \mu_2 & \cdots & x_1^{(d)} - \mu_d \\ x_2^{(1)} - \mu_1 & x_2^{(2)} - \mu_2 & \cdots & x_2^{(d)} - \mu_d \\ \vdots & \vdots & \ddots & \vdots \\ x_n^{(1)} - \mu_1 & x_n^{(2)} - \mu_2 & \cdots & x_n^{(d)} - \mu_d \end{bmatrix}

In X~\tilde X, the mean of each column is 0, which means that it’s reasonable to approximate the dataset using a line through the origin. Importantly, this doesn’t change the shape of the dataset, meaning it doesn’t change the “best direction” we’re searching for.

Finding the Best Direction

Now, X~\tilde X is the n×dn \times d matrix of centered data, and we’re ready to return to our goal of finding the best direction. Let x~i∈Rd\tilde x_i \in \mathbb{R}^d refer to row ii of X~\tilde X. Note that I’ve dropped the vector symbol on top of x~i\tilde x_i, not because it’s not a vector, but just because drawing a tilde on top of a vector hat is a little funky: x~i⃗\vec{\tilde x_i}.

X~=[x1(1)−μ1x1(2)−μ2⋯x1(d)−μdx2(1)−μ1x2(2)−μ2⋯x2(d)−μd⋮⋮⋱⋮xn(1)−μ1xn(2)−μ2⋯xn(d)−μd]=[–x~1T––x~2T–⋮–x~nT–]\tilde X = \begin{bmatrix} x_1^{(1)} - \mu_1 & x_1^{(2)} - \mu_2 & \cdots & x_1^{(d)} - \mu_d \\ x_2^{(1)} - \mu_1 & x_2^{(2)} - \mu_2 & \cdots & x_2^{(d)} - \mu_d \\ \vdots & \vdots & \ddots & \vdots \\ x_n^{(1)} - \mu_1 & x_n^{(2)} - \mu_2 & \cdots & x_n^{(d)} - \mu_d \end{bmatrix} = \begin{bmatrix} – & \tilde x_1^T & – \\ – & \tilde x_2^T – \\ & \vdots & \\ – & \tilde x_n^T & – \end{bmatrix}

Our goal is to find the unit vector v⃗∈Rd\vec v \in \mathbb{R}^d that minimizes the mean squared orthogonal projection error. That is, if the orthogonal projection of x~i\tilde x_i onto v⃗\vec v is p⃗i\vec p_i, then we want to find the unit vector v⃗\vec v that minimizes

J(v⃗)=1n∑i=1n∥x~i−p⃗i∥2=1n∑i=1n∥x~i−(x~i⋅v⃗)v⃗∥2J(\vec v) = \frac{1}{n}\sum_{i = 1}^n \lVert \tilde x_i - \vec p_i \rVert^2 = \frac{1}{n}\sum_{i = 1}^n \lVert \tilde x_i - (\tilde x_i \cdot \vec v) \vec v \rVert^2

I’ve chosen the letter JJ since there is no standard notation for this quantity, and because JJ is a commonly used letter in optimization problems. It’s important to note that JJ is not a loss function, nor is it an empirical risk, since it has nothing to do with prediction.

How do we find the unit vector v⃗\vec v that minimizes J(v⃗)J(\vec v)? We can start by expanding the squared norm, using the fact that ∥a⃗−b⃗∥2=∥a⃗∥2−2a⃗⋅b⃗+∥b⃗∥2\lVert \vec a - \vec b \rVert^2 = \lVert \vec a \rVert^2 - 2 \vec a \cdot \vec b + \lVert \vec b \rVert^2:

J(v⃗)=1n∑i=1n∥x~i−(x~i⋅v⃗)v⃗∥2=1n∑i=1n(∥x~i∥2−2x~i⋅(x~i⋅v⃗)v⃗+((x~i⋅v⃗)v⃗)⋅((x~i⋅v⃗)v⃗))=1n∑i=1n(∥x~i∥2−2x~i⋅(x~i⋅v⃗)v⃗+(x~i⋅v⃗)2v⃗⋅v⃗⏟=1)=1n∑i=1n(∥x~i∥2−2(x~i⋅v⃗)2+(x~i⋅v⃗)2)=1n∑i=1n(∥x~i∥2−(x~i⋅v⃗)2)=1n∑i=1n∥x~i∥2−1n∑i=1n(x~i⋅v⃗)2\begin{align*} J(\vec v) &= \frac{1}{n}\sum_{i = 1}^n \lVert \tilde x_i - (\tilde x_i \cdot \vec v) \vec v \rVert^2 \\ &= \frac{1}{n} \sum_{i=1}^n \left( \lVert \tilde x_i \rVert^2 - 2 \tilde x_i \cdot (\tilde x_i \cdot \vec v) \vec v + \left( (\tilde x_i \cdot \vec v) \vec v \right) \cdot \left( (\tilde x_i \cdot \vec v) \vec v \right) \right) \\ &= \frac{1}{n} \sum_{i=1}^n \left( \lVert \tilde x_i \rVert^2 - 2 \tilde x_i \cdot (\tilde x_i \cdot \vec v) \vec v + (\tilde x_i \cdot \vec v)^2 \underbrace{\vec v \cdot \vec v}_{=1}\right) \\ &= \frac{1}{n} \sum_{i=1}^n \left( \lVert \tilde x_i \rVert^2 - 2 (\tilde x_i \cdot \vec v)^2 + (\tilde x_i \cdot \vec v)^2 \right) \\ &= \frac{1}{n} \sum_{i=1}^n \left( \lVert \tilde x_i \rVert^2 - (\tilde x_i \cdot \vec v)^2 \right) \\ &= \frac{1}{n} \sum_{i = 1}^n \lVert \tilde x_i \rVert^2 - \frac{1}{n} \sum_{i = 1}^n (\tilde x_i \cdot \vec v)^2 \end{align*}

The first term, 1n∑i=1n∥x~i∥2\frac{1}{n} \sum_{i = 1}^n \lVert \tilde x_i \rVert^2, is constant with respect to v⃗\vec v, meaning we can ignore it for the purposes of minimizing J(v⃗)J(\vec v). Our focus, then, is on minimizing the second term,

−1n∑i=1n(x~i⋅v⃗)2- \frac{1}{n} \sum_{i = 1}^n (\tilde x_i \cdot \vec v)^2

But, since the second term has a negative sign in front of it, minimizing it is equivalent to maximizing its positive counterpart,

1n∑i=1n(x~i⋅v⃗)2\frac{1}{n} \sum_{i = 1}^n (\tilde x_i \cdot \vec v)^2

Maximizing Variance

Finding the best direction is equivalent to finding the unit vector v⃗\vec v that maximizes the mean squared value of x~i⋅v⃗\tilde x_i \cdot \vec v. What is this quantity, though? Remember, the projection of x~i\tilde x_i onto v⃗\vec v is given by

p⃗i=(x~i⋅v⃗)v⃗\vec p_i = (\tilde x_i \cdot \vec v) \vec v

All nn of our original data points are projected onto v⃗\vec v, meaning that each x~i\tilde x_i is approximated by some scalar multiple (x~i⋅v⃗\tilde x_i \cdot \vec v) of v⃗\vec v. These scalar multiples are our new feature values – that is, the values of the first principal component (first because this is the “best” direction).

Image produced in Jupyter

The values x~1⋅v⃗,x~2⋅v⃗,…,x~n⋅v⃗\tilde x_1 \cdot \vec v, \tilde x_2 \cdot \vec v, \ldots, \tilde x_n \cdot \vec v – −2.62\color{orange}-2.62, −1.49\color{orange}-1.49, −1.14\color{orange}-1.14, −1.11\color{orange}-1.11, and 2.11\color{orange}2.11 above – describe the position of point x~i\tilde x_i on the line defined by v⃗\vec v. (I’ve highlighted just five of the points to make the picture less cluttered, but in reality every single blue point has a corresponding orange dot above!)

These coefficients on v⃗\vec v can be thought of as positions on a number line.

Image produced in Jupyter

Effectively, we’ve rotated the data so that the orange points are on the “new” xx-axis. (Rotations... where have we heard of those before? The connection is coming.)

So, we have a grasp of what the values x~1⋅v⃗,x~2⋅v⃗,…,x~n⋅v⃗{\color{orange}\tilde x_1 \cdot \vec v, \tilde x_2 \cdot \vec v, \ldots, \tilde x_n \cdot \vec v} describe. But why are they important?

The quantity we’re trying to maximize,

1n∑i=1n(x~i⋅v⃗)2\frac{1}{n} \sum_{i = 1}^n (\tilde x_i \cdot \vec v)^2

is the variance of the values x~1⋅v⃗,x~2⋅v⃗,…,x~n⋅v⃗\color{orange}\tilde x_1 \cdot \vec v, \tilde x_2 \cdot \vec v, \ldots, \tilde x_n \cdot \vec v. It’s not immediately obvious why this is the case. Remember, the variance of a list of numbers is

1n∑i=1n(valuei−mean)2\frac{1}{n} \sum_{i = 1}^n (\text{value}_i - \text{mean})^2

(I’m intentionally using vague notation to avoid conflating singular values with standard deviations.)

The quantity 1n∑i=1n(x~i⋅v⃗)2\frac{1}{n} \sum_{i = 1}^n (\tilde x_i \cdot \vec v)^2 resembles this definition of variance, but it doesn’t seem to involve a subtraction by a mean. This is because the mean of x~1⋅v⃗,x~2⋅v⃗,…,x~n⋅v⃗\tilde x_1 \cdot \vec v, \tilde x_2 \cdot \vec v, \ldots, \tilde x_n \cdot \vec v is 0! This seems to check out visually above: it seems that if we found the corresponding orange point for every single blue point, the orange points would be symmetrically distributed around 0.

Proof that the mean of x~1⋅v⃗,x~2⋅v⃗,…,x~n⋅v⃗\tilde x_1 \cdot \vec v, \tilde x_2 \cdot \vec v, \ldots, \tilde x_n \cdot \vec v is 0

Recall that each x~i\tilde x_i is a vector in Rd\mathbb{R}^d and corresponds to a row in X~\tilde X. The values x~1⋅v⃗,x~2⋅v⃗,…,x~n⋅v⃗\tilde x_1 \cdot \vec v, \tilde x_2 \cdot \vec v, \ldots, \tilde x_n \cdot \vec v can all be computed at once through the matrix-vector product X~v⃗\tilde X \vec v (remember that to multiply a matrix by a vector, we dot product the rows of the matrix with the vector).

The columns of X~\tilde X each have a mean of 0. As we discovered in Activity 2, this means that

X~T1⃗=[sum of column 1 of X~sum of column 2 of X~⋮sum of column d of X~]=[00⋮0]=0⃗d\tilde X^T \vec 1 = \begin{bmatrix} \text{sum of column 1 of } \tilde X \\ \text{sum of column 2 of } \tilde X \\ \vdots \\ \text{sum of column $d$ of } \tilde X\end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \end{bmatrix} = \vec 0_{d}

if 1⃗∈Rn\vec 1 \in \mathbb{R}^n is a vector of all ones.

So, what is the mean of the values x~1⋅v⃗,x~2⋅v⃗,…,x~n⋅v⃗\tilde x_1 \cdot \vec v, \tilde x_2 \cdot \vec v, \ldots, \tilde x_n \cdot \vec v? This mean, by definition, is 1n∑i=1nx~i⋅v⃗\frac{1}{n} \sum_{i = 1}^n \tilde x_i \cdot \vec v, but this is equalx to 1n(1⃗T(X~v⃗))\frac{1}{n} \left( \vec 1^T (\tilde X \vec v) \right), once again, using the fact that the dot product with the ones vector is equivalent to summing the values of the other vector.

1n∑i=1n(x~i⋅v⃗)=1n(1⃗T(X~v⃗))=1n(v⃗TX~T1⃗)=1n(v⃗T0⃗)=0\frac{1}{n} \sum_{i = 1}^n (\tilde x_i \cdot \vec v) = \frac{1}{n} \left( \vec 1^T (\tilde X \vec v) \right) = \frac{1}{n} \left( \vec v^T \tilde X^T \vec 1 \right) = \frac{1}{n} \left( \vec v^T \vec 0 \right) = 0

So, the mean of the values x~1⋅v⃗,x~2⋅v⃗,…,x~n⋅v⃗\tilde x_1 \cdot \vec v, \tilde x_2 \cdot \vec v, \ldots, \tilde x_n \cdot \vec v is 0.

Since the mean of the values x~1⋅v⃗,x~2⋅v⃗,…,x~n⋅v⃗\tilde x_1 \cdot \vec v, \tilde x_2 \cdot \vec v, \ldots, \tilde x_n \cdot \vec v is 0, it is true that the quantity we’re maximizing is the variance of the values x~1⋅v⃗,x~2⋅v⃗,…,x~n⋅v⃗\tilde x_1 \cdot \vec v, \tilde x_2 \cdot \vec v, \ldots, \tilde x_n \cdot \vec v, which I’ll refer to as

PV(v⃗)=1n∑i=1n(x~i⋅v⃗)2\text{PV}(\vec v) = \frac{1}{n} \sum_{i = 1}^n (\tilde x_i \cdot \vec v)^2

PV stands for “projected variance”; it is the variance of the scalar projections of the rows of X~\tilde X onto v⃗\vec v. This is non-standard notation, but I wanted to avoid calling this the “variance of v⃗\vec v”, since it’s not the variance of the components of v⃗\vec v.

The Duality

We’ve stumbled upon a beautiful duality – that is, two equivalent ways to think about the same problem. Given a centered n×dn \times d matrix X~\tilde X with rows x~i∈Rd\tilde x_i \in \mathbb{R}^d, the most important direction v⃗\vec v is the one that

  1. Minimizes the mean squared orthogonal error of projections onto it,

    J(v⃗)=1n∑i=1n∥x~i−p⃗i∥2J(\vec v) = \frac{1}{n} \sum_{i = 1}^n \lVert \tilde x_i - \vec p_i \rVert^2

  2. Equivalently, maximizes the variance of the projected values x~i⋅v⃗\tilde x_i \cdot \vec v, i.e. maximizes

PV(v⃗)=1n∑i=1n(x~i⋅v⃗)2\text{PV}(\vec v) = \frac{1}{n} \sum_{i = 1}^n (\tilde x_i \cdot \vec v)^2

Minimizing orthogonal error is equivalent to maximizing variance: don’t forget this. In the widget below, drag the slider to see this duality in action: the shorter the orthogonal errors tend to be, the more spread out the orange points tend to be (meaning their variance is large). When the orange points are close together, the orthogonal errors tend to be quite large.

Loading...

Here’s a related StackExchange post with a similar animation.

Minimizing mean squared orthogonal error is a useful goal: this brings our 1-dimensional approximation of the data as close as possible to the original data. Why is maximizing variance another useful goal?

Here’s an analogy: suppose I had a dataset with several characteristics about each student in this course, and I needed to come up with a single number to describe each student. Perhaps I’d take some weighted combination of the number of lectures they’ve attended, their Midterm 1 and Midterm 2 scores, and the number of office hours they’ve attended. I wouldn’t factor in the number of semesters they’ve been at Michigan, since this number is likely to be similar across all students, while these other characteristics vary more significantly.

The guiding philosophy here is that variance is information. I want to come up with a linear combination of features that remains as different as possible across different data points, as this will allow me to distinguish between different data points, i.e. most faithfully represent the data, while being in a lower-dimensional space.

Optimization Constraints

Okay, so we know that minimizing mean squared orthogonal error is equivalent to maximizing variance. Presumably, this interpretation will help us in our quest of finding the most important direction, v⃗\vec v.

Again, we want to maximize

PV(v⃗)=1n∑i=1n(x~i⋅v⃗)2\text{PV}(\vec v) = \frac{1}{n} \sum_{i = 1}^n (\tilde x_i \cdot \vec v)^2

PV(v⃗)\text{PV}(\vec v) is a vector-to-scalar function; it takes in vectors in Rd\mathbb{R}^d and returns a scalar. So, to maximize it, we’d need to take its gradient with respect to v⃗\vec v. (The same would have been true if we wanted to minimize the mean squared orthogonal error – we’d need to have taken the gradient of J(v⃗)J(\vec v). The point of highlighting the duality above, in addition to helping build intuition, was to address the fact that J(v⃗)J(\vec v) doesn’t have a simple gradient, but PV(v⃗)\text{PV}(\vec v) does.)

Is there a way we can write PV(v⃗)\text{PV}(\vec v) in a way that would allow us to apply the gradient rules we developed in Chapter 4? Using the observation that the values x~1⋅v⃗,x~2⋅v⃗,…,x~n⋅v⃗\tilde x_1 \cdot \vec v, \tilde x_2 \cdot \vec v, \ldots, \tilde x_n \cdot \vec v are all components of the matrix-vector product X~v⃗\tilde X \vec v, we can rewrite PV(v⃗)\text{PV}(\vec v) as

PV(v⃗)=1n∥X~v⃗∥2\text{PV}(\vec v) = \frac{1}{n} \lVert \tilde X \vec v \rVert^2

which can further be simplified to

PV(v⃗)=1nv⃗TX~TX~v⃗\text{PV}(\vec v) = \frac{1}{n} \vec v^T \tilde X^T \tilde X \vec v

Now, we can apply the gradient rules in Chapter 4 – specifically, the quadratic form rule – to find ∇PV(v⃗)\nabla \text{PV}(\vec v). But, there’s something we’ve forgotten. Let’s take a step back.

Remember, v⃗\vec v must be a unit vector! All of the optimization problems we’ve seen so far in this course have been unconstrained, meaning that there were never restrictions on what the optimal solutions could be. This is our first time encountering a constrained optimization problem. We can write this optimization problem as

max 1n∥X~v⃗∥2subject to ∥v⃗∥=1\begin{aligned} \text{max} \quad\ & \frac{1}{n} \lVert \tilde X \vec v \rVert^2 \\ \text{subject to} \quad\ & \lVert \vec v \rVert = 1 \end{aligned}

If we don’t constrain v⃗\vec v to be a unit vector, then PV(v⃗)\text{PV}(\vec v) is unbounded: the larger I make the components in v⃗\vec v, the larger PV(v⃗)\text{PV}(\vec v) becomes.

One general solution for solving constrained optimization problems is to use the method Lagrange multipliers, which are typically covered in a Calculus 3 course. We won’t cover this approach here. Instead, we’ll use another approach that works for this specific problem. Instead of maximizing

PV(v⃗)=1n∥X~v⃗∥2\text{PV}(\vec v) = \frac{1}{n} \lVert \tilde X \vec v \rVert^2

we’ll maximize a related function that automatically constrains v⃗\vec v to be a unit vector, without us having to write a constraint.

f(v⃗)=∥X~v⃗∥2∥v⃗∥2f(\vec v) = \frac{\lVert \tilde X \vec v \rVert^2}{\lVert \vec v \rVert^2}

Why is maximizing f(v⃗)f(\vec v) equivalent to maximizing PV(v⃗)\text{PV}(\vec v) with the unit vector constraint?

  • First, note that the factor of 1n\frac{1}{n} from PV(v⃗)\text{PV}(\vec v) has been dropped: maximizing 1n∥X~v⃗∥2\frac{1}{n} \lVert \tilde X \vec v \rVert^2 is equivalent to maximizing ∥X~v⃗∥2\lVert \tilde X \vec v \rVert^2. The same vector wins in both cases. So, for convenience, we’ll drop the 1n\frac{1}{n} for now.

  • Now, suppose v⃗\vec v is a unit vector. Then, ∥v⃗∥2=1\lVert \vec v \rVert^2 = 1, so f(v⃗)=∥X~v⃗∥2f(\vec v) = \lVert \tilde X \vec v \rVert^2.

  • If v⃗\vec v is not a unit vector, then v⃗=cu⃗\vec v = c \vec u, where u⃗\vec u is some unit vector. But then

    f(v⃗)=∥X~v⃗∥2∥v⃗∥2=∥X~(cu⃗)∥2∥cu⃗∥2=∥cX~u⃗∥2c2∥u⃗∥2=c2∥X~u⃗∥2c2∥u⃗∥2=∥X~u⃗∥2∥u⃗∥2=f(u⃗)f(\vec v) = \frac{\lVert \tilde X \vec v \rVert^2}{\lVert \vec v \rVert^2} = \frac{\lVert \tilde X (c \vec u) \rVert^2}{\lVert c \vec u \rVert^2} = \frac{\lVert c \tilde X \vec u \rVert^2}{c^2 \lVert \vec u \rVert^2} = \frac{c^2 \lVert \tilde X \vec u \rVert^2}{c^2 \lVert \vec u \rVert^2} = \frac{\lVert \tilde X \vec u \rVert^2}{\lVert \vec u \rVert^2} = f(\vec u)

    which means that f(v⃗)=f(u⃗)f(\vec v) = f(\vec u), i.e. f(v⃗)f(\vec v) is equal to ff applied to a unit vector with the same direction as v⃗\vec v.

Dividing by ∥v⃗∥2\lVert \vec v \rVert^2 effectively neutralizes the effect of the length of v⃗\vec v and instead focuses on the direction of v⃗\vec v, which is precisely what we’re looking for.

The Rayleigh Quotient

Now, our problem is to find the vector v⃗\vec v that maximizes

f(v⃗)=∥X~v⃗∥2∥v⃗∥2f(\vec v) = \frac{\lVert \tilde X \vec v \rVert^2}{\lVert \vec v \rVert^2}

Here’s the big idea: we’ve seen this function before! Specifically, in Homework 9, Problem 4, you studied the Rayleigh quotient, which is the function

g(v⃗)=v⃗TAv⃗v⃗Tv⃗g(\vec v) = \frac{\vec v^T A \vec v}{\vec v^T \vec v}

for a symmetric n×nn \times n matrix AA. There, you showed that the gradient of g(v⃗)g(\vec v) is given by

∇g(v⃗)=2v⃗Tv⃗(Av⃗−g(v⃗)v⃗)\nabla g(\vec v) = \frac{2}{\vec v^T \vec v} \left( A \vec v - g(\vec v) \vec v \right)

Which vector v⃗\vec v maximizes g(v⃗)g(\vec v)? If v⃗\vec v maximizes g(v⃗)g(\vec v), then ∇g(v⃗)=0⃗\nabla g(\vec v) = \vec 0.

2v⃗Tv⃗(Av⃗−g(v⃗)v⃗)=0⃗  ⟹  Av⃗−g(v⃗)v⃗=0⃗  ⟹  Av⃗=g(v⃗)v⃗\frac{2}{\vec v^T \vec v} \left( A \vec v - g(\vec v) \vec v \right) = \vec 0 \implies A \vec v - g(\vec v) \vec v = \vec 0 \implies A \vec v = g(\vec v) \vec v

What this is saying is that if ∇g(v⃗)=0⃗\nabla g(\vec v) = \vec 0, then v⃗\vec v is an eigenvector of AA!

As an example – and this is a total diversion from the main storyline of the lecture, but it’s an important one – let’s visualize the Rayleigh quotient for the symmetric 2×22 \times 2 matrix

A=[1551]A = \begin{bmatrix} 1 & 5 \\ 5 & 1 \end{bmatrix}

If we let v⃗=[xy]\vec v = \begin{bmatrix} x \\ y\end{bmatrix}, then this Rayleigh quotient is

g(v⃗)=v⃗TAv⃗v⃗Tv⃗=x2+10xy+y2x2+y2g(\vec v) = \frac{\vec v^T A \vec v}{\vec v^T \vec v} = \frac{x^2 + 10xy + y^2}{x^2 + y^2}

I have visualized it as a contour plot, a concept we first encountered in Chapter 4.1.

Loading...

Dark red regions correspond to the highest values of the Rayleigh quotient, and dark blue regions correspond to the lowest values. You’ll notice that the darkest red and darkest blue regions form lines: these lines are the directions of the eigenvectors of AA! The corresponding outputs along those lines are the eigenvalues of AA. Take a look at the scale on the colorbar: the graph maxes out at 6, which is the largest eigenvalue of AA, and its minimum is -4, which is the smallest eigenvalue of AA.

So, the fact that

∇g(v⃗)=0⃗  ⟹  Av⃗=g(v⃗)v⃗\nabla g(\vec v) = \vec 0 \implies A \vec v = g(\vec v) \vec v

means that:

  • the vector that maximizes g(v⃗)=v⃗TAv⃗v⃗Tv⃗g(\vec v) = \frac{\vec v^T A \vec v}{\vec v^T \vec v} is the eigenvector of AA corresponding to the largest eigenvalue of AA, and

  • the vector that minimizes g(v⃗)g(\vec v) is the eigenvector of AA corresponding to the smallest eigenvalue of AA!

How does this relate to our problem, of finding the vector that maximizes PV(v⃗)\text{PV}(\vec v), which is the same vector that maximizes

f(v⃗)=∥X~v⃗∥2∥v⃗∥2f(\vec v) = \frac{\lVert \tilde X \vec v \rVert^2}{\lVert \vec v \rVert^2}

?

Here’s the key: Our f(v⃗)f(\vec v) is a Rayleigh quotient, where A=X~TX~A = \tilde X^T \tilde X.

f(v⃗)=∥X~v⃗∥2∥v⃗∥2=v⃗TX~TX~v⃗v⃗Tv⃗f(\vec v) = \frac{\lVert \tilde X \vec v \rVert^2}{\lVert \vec v \rVert^2} = \frac{\vec v^T {\color{#3d81f6} \tilde X^T \tilde X} \vec v}{\vec v^T \vec v}

X~TX~\tilde X^T \tilde X is a symmetric matrix, so everything we know about the Rayleigh quotient applies. That is,

∇f(v⃗)=2v⃗Tv⃗(X~TX~v⃗−f(v⃗)v⃗)\nabla f(\vec v) = \frac{2}{\vec v^T \vec v} \left( \tilde X^T \tilde X \vec v - f(\vec v) \vec v \right)

and

∇f(v⃗)=0⃗  ⟹  X~TX~v⃗=f(v⃗)v⃗⏟if v⃗ is a critical point, it must be an eigenvector of X~TX~\underbrace{\nabla f(\vec v) = \vec 0 \implies \tilde X^T \tilde X \vec v = f(\vec v) \vec v}_{\text{if } \vec v \text{ is a critical point, it must be an eigenvector of } \tilde X^T \tilde X}

Using the same reasoning as before, the vector that maximizes f(v⃗)f(\vec v) is the eigenvector of X~TX~\tilde X^T \tilde X corresponding to the largest eigenvalue of X~TX~\tilde X^T \tilde X!

What does this really mean?

Image produced in Jupyter

It says that the line that maximizes the variance of the projected data – which is also the line that minimizes the mean squared orthogonal error – is spanned by the eigenvector of X~TX~\tilde X^T \tilde X corresponding to the largest eigenvalue of X~TX~\tilde X^T \tilde X.

Principal Components

We’ve discovered that the single best direction to project the data onto is the eigenvector of X~TX~\tilde X^T \tilde X corresponding to the largest eigenvalue of X~TX~\tilde X^T \tilde X. Let’s call this eigenvector v⃗1\vec v_1, as I did in the diagram above. More v⃗\vec v’s are coming.

PV(v⃗)=1n∥X~v⃗∥2  ⟹  v⃗∗=v⃗1\text{PV}(\vec v) = \frac{1}{n} \lVert \tilde X \vec v \rVert^2 \implies \vec v^* = \vec v_1

The values of the first principal component, PC1\text{PC}_1 (i.e. new feature 1), come from projecting each row of X~\tilde X onto v⃗1\vec v_1:

PC1=[x~1⋅v⃗1x~2⋅v⃗1⋮x~n⋅v⃗1]=X~v⃗1\text{PC}_1 = \begin{bmatrix} \tilde x_1 \cdot \vec v_1 \\ \tilde x_2 \cdot \vec v_1 \\ \vdots \\ \tilde x_n \cdot \vec v_1 \end{bmatrix} = \tilde X \vec v_1

This projection of our data onto the line spanned by v⃗1\vec v_1 is the linear combination of the original features that captures the largest possible amount of variance in the data.

(Whenever you see “principal component”, you should think “new feature”.)

Remember, one of our initial goals was to find multiple principal components that were uncorrelated with each other, meaning their correlation coefficient rr is 0. We’ve found the first principal component, PC1\text{PC}_1, which came from projecting onto the best direction. This is the single best linear combination of the original features we can come up with.

Let’s take a greedy approach. I now want to find the next best principal component, PC2\text{PC}_2, which should be uncorrelated with PC1\text{PC}_1. PC2\text{PC}_2 should capture all of the remaining variance in the data that PC1\text{PC}_1 couldn’t capture. Intuitively, since our dataset is 2-dimensional, together, PC1\text{PC}_1 and PC2\text{PC}_2 should contain the same information as the original two features.

PC2\text{PC}_2 comes from projecting onto the best direction among all directions that are orthogonal to v⃗1\vec v_1. It can be shown that this “second-best direction” is the eigenvector v⃗2\vec v_2 of X~TX~\tilde X^T \tilde X corresponding to the second largest eigenvalue of X~TX~\tilde X^T \tilde X.

PC2=X~v⃗2\text{PC}_2 = \tilde X \vec v_2

In other words, the vector v⃗\vec v that maximizes f(v⃗)=∥X~v⃗∥2∥v⃗∥2\displaystyle f(\vec v) = \frac{\lVert \tilde X \vec v \rVert^2}{\lVert \vec v \rVert^2} subject to the constraint that v⃗\vec v is orthogonal to v⃗1\vec v_1, is v⃗2\vec v_2. The proof of this is beyond the scope of what we’ll discuss here, as it involves some constrained optimization theory.

Image produced in Jupyter

Why is v⃗2\vec v_2 orthogonal to v⃗1\vec v_1? They are both eigenvectors of X~TX~\tilde X^T \tilde X corresponding to different eigenvalues, so they must be orthogonal, thanks to the spectral theorem (which applies because X~TX~\tilde X^T \tilde X is symmetric). Remember that while any vector on the line spanned by v⃗2\vec v_2 is also an eigenvector of X~TX~\tilde X^T \tilde X corresponding to the second largest eigenvalue, we pick the specific v⃗2\vec v_2 that is a unit vector.

v⃗2\vec v_2 and v⃗1\vec v_1 being orthogonal means that X~v⃗2\tilde X \vec v_2 and X~v⃗1\tilde X \vec v_1 are orthogonal, too. But because the columns of X~\tilde X are mean centered, we can show that this implies that the correlation between PC1\text{PC}_1 and PC2\text{PC}_2 is 0. I’ll save the algebra for now, but see if you can work this out yourself. You’ve done similar proofs in homeworks already.

The SVD Returns

The moment I said that v⃗1\vec v_1 is an eigenvector of X~TX~\tilde X^T \tilde X, something should have been ringing in your head: v⃗1\vec v_1 is a singular vector of X~\tilde X! Recall, if

X~=[∣⋯∣∣⋯∣u⃗1⋯u⃗ru⃗r+1⋯u⃗n∣⋯∣∣⋯∣]⏟U[σ10⋯0⋱⋮⋮σr0⋯00⋯00⋯0⋮⋮⋮⋮0⋯00⋯0]⏟Σ[—  v⃗1T  —⋮—  v⃗rT  ——  v⃗r+1T  —⋮—  v⃗dT  —]⏟VT\tilde X = \underbrace{\begin{bmatrix} {\color{#3d81f6} |} & {\color{#3d81f6} \cdots} & {\color{#3d81f6} |} & {\color{#aaaaaa} |} & {\color{#aaaaaa} \cdots} & {\color{#aaaaaa} |} \\ {\color{#3d81f6} \vec u_1} & {\color{#3d81f6} \cdots} & {\color{#3d81f6} \vec u_r} & {\color{#aaaaaa} \vec u_{r+1}} & {\color{#aaaaaa} \cdots} & {\color{#aaaaaa} \vec u_n} \\ {\color{#3d81f6} |} & {\color{#3d81f6} \cdots} & {\color{#3d81f6} |} & {\color{#aaaaaa} |} & {\color{#aaaaaa} \cdots} & {\color{#aaaaaa} |} \end{bmatrix}}_{U} \underbrace{\begin{bmatrix} {\color{orange} \sigma_1} & {\color{#aaaaaa} } & {\color{#aaaaaa} } & {\color{#aaaaaa} 0} & {\color{#aaaaaa} \cdots} & {\color{#aaaaaa} 0} \\ {\color{#aaaaaa} } & {\color{orange} \ddots} & {\color{#aaaaaa} } & {\color{#aaaaaa} \vdots} & {\color{#aaaaaa} } & {\color{#aaaaaa} \vdots} \\ {\color{#aaaaaa} } & {\color{#aaaaaa} } & {\color{orange} \sigma_r} & {\color{#aaaaaa} 0} & {\color{#aaaaaa} \cdots} & {\color{#aaaaaa} 0} \\[6pt] {\color{#aaaaaa} 0} & {\color{#aaaaaa} \cdots} & {\color{#aaaaaa} 0} & {\color{#aaaaaa} 0} & {\color{#aaaaaa} \cdots} & {\color{#aaaaaa} 0} \\ {\color{#aaaaaa} \vdots} & {\color{#aaaaaa} } & {\color{#aaaaaa} \vdots} & {\color{#aaaaaa} \vdots} & {\color{#aaaaaa} } & {\color{#aaaaaa} \vdots} \\ {\color{#aaaaaa} 0} & {\color{#aaaaaa} \cdots} & {\color{#aaaaaa} 0} & {\color{#aaaaaa} 0} & {\color{#aaaaaa} \cdots} & {\color{#aaaaaa} 0} \end{bmatrix}}_{\Sigma} \underbrace{\begin{bmatrix} {\color{#d81a60}\text{---} \:\: \vec v_1^T \:\: \text{---}} \\ {\color{#d81a60} \vdots} \\ {\color{#d81a60}\text{---} \:\: \vec v_r^T \:\: \text{---}} \\ \\ {\color{#aaaaaa} \text{---} \:\: \vec v_{r+1}^T \:\: \text{---}} \\ {\color{#aaaaaa} \vdots} \\ {\color{#aaaaaa} \text{---} \:\: \vec v_d^T \:\: \text{---}} \\ \end{bmatrix}}_{V^T}

is the singular value decomposition of X~\tilde X, then the columns of VV are the eigenvectors of X~TX~\tilde X^T \tilde X, all of which are orthogonal to each other and have unit length.

We arranged the components in the singular value decomposition in decreasing order of singular values of X~\tilde X, which are the square roots of the eigenvalues of X~TX~\tilde X^T \tilde X.

σi=λi,where λi is the i-th largest eigenvalue of X~TX~\sigma_i = \sqrt{\lambda_i}, \quad \text{where } \lambda_i \text{ is the } i\text{-th largest eigenvalue of } \tilde X^T \tilde X

So,

  • the first column of VV is v⃗1\vec v_1, the “best direction” to project the data onto,

  • the second column of VV is v⃗2\vec v_2, the “second-best direction” to project the data onto,

  • the third column of VV is v⃗3\vec v_3, the “third-best direction” to project the data onto,

  • and so on.

Let’s put this in context with a few examples.

Example: From R2\mathbb{R}^2 to R2\mathbb{R}^2

Let’s start with the 2-dimensional dataset of flipper length vs. body mass from the penguins dataset.

Image produced in Jupyter

The first thing you should notice is that while the original data points seem to have some positive correlation, the principal components are uncorrelated! This is a good thing; it’s what we wanted as a design goal. In effect, by converting from the original features to principal components, we’ve rotated the data to remove the correlation between the features.

I have picked two points to highlight, points 5 and 12 in the original data. The coloring in red and purple is meant to show you how original (centered) data points translate to points in the principal component space.

Notice the scale of the data: PC1\text{PC}_1 axis is much longer than the PC2\text{PC}_2 axis, since the first principal component captures much more variance than the second. We will make this notion – of the proportion of variance captured by each principal component – more precise soon.

(You’ll notice that the body mass and flipper length values on the left are much smaller than in the original datasets; in the original dataset, the body mass values were in the thousands, which distorted the scale of the plot and made it hard to see that the two eigenvectors are indeed orthogonal.)

Example: Starting from R3\mathbb{R}^3

The real power of PCA reveals itself when we start with high-dimensional data. Suppose we start with three of the features in the penguins dataset: bill_depth_mm, flipper_length_mm, and body_mass_g – and want to reduce the dimensionality of the data to 1 or 2. Points are colored by their species.

Loading...

Observe that penguins of the same species tend to be clustered together. This, alone, has nothing to do with PCA: we happen to have this information, so I’ve included it in the plot.

If XX is the 333×3333 \times 3 matrix of three features, and X~\tilde X is the mean-centered version, the best directions in which to project the data are the columns of VV in X~=UΣVT\tilde X = U \Sigma V^T.

X_three_features
Loading...
X_three_features.mean()
bill_depth_mm 17.164865 flipper_length_mm 200.966967 body_mass_g 4207.057057 dtype: float64
# This is the mean-centered version of X_three_features!
X_three_features - X_three_features.mean()
Loading...
X_three_features_centered = X_three_features - X_three_features.mean()
u, s, vt = np.linalg.svd(X_three_features_centered)

Now, the rows of VTV^T (vt), which are the columns of VV (vt.T), contain the best directions.

vt
array([[-1.15433983e-03, 1.51946036e-02, 9.99883889e-01], [ 1.02947493e-01, -9.94570148e-01, 1.52327042e-02], [ 9.94686122e-01, 1.02953123e-01, -4.16174416e-04]])

It’s important to remember that these “best directions” are nothing more than linear combinations of the original features. Since the first row of VTV^T is [−0.001150.015190.99988]\begin{bmatrix} -0.00115 & 0.01519 & 0.99988 \end{bmatrix}, the first principal component is

PC 1i=−0.00115⋅bill depthi+0.01519⋅flipper lengthi+0.99988⋅body massi\text{PC 1}_i = -0.00115 \cdot \text{bill depth}_i + 0.01519 \cdot \text{flipper length}_i + 0.99988 \cdot \text{body mass}_i

while the second is

PC 2i=0.10294⋅bill depthi−0.99457⋅flipper lengthi+0.01523⋅body massi\text{PC 2}_i = 0.10294 \cdot \text{bill depth}_i - 0.99457 \cdot \text{flipper length}_i + 0.01523 \cdot \text{body mass}_i

and third is

PC 3i=0.99468⋅bill depthi+0.10295⋅flipper lengthi−0.00042⋅body massi\text{PC 3}_i = 0.99468 \cdot \text{bill depth}_i + 0.10295 \cdot \text{flipper length}_i - 0.00042 \cdot \text{body mass}_i

where i=1,2,…,333i = 1, 2, \ldots, 333 is the index of the penguin.

To compute all three of these principal components at once, for every penguin, we just need to compute X~V\tilde X V.

pcs = X_three_features_centered @ vt.T
pcs
Loading...

Let’s plot the first two principal components: that is, for each of our 333 penguins, we’ll plot their value of X~v⃗1\tilde X \vec v_1 on the xx-axis and their value of X~v⃗2\tilde X \vec v_2 on the yy-axis.

Image produced in Jupyter

This is the best 2-dimensional projection of our 3-dimensional scatter plot. And here’s the kicker: penguins of the same species STILL tend to be clustered together in the principal component space!

What this tells us is that our technique for taking linear combinations of the original features is good at preserving the important information in the original dataset. We went from writing down 3 numbers per penguin to 2, but it seems that we didn’t lose much information per penguin.

What do I mean by “much important information”? Let’s make this idea more precise.

Explained Variance

The goal of PCA is to find new features – principal components – that capture as much of the variation in the data as possible, while being uncorrelated with each other.

PCA isn’t foolproof: it works better on some datasets than others. If the features we’re working with are already uncorrelated, PCA isn’t useful. And, even for datasets that are well suited for PCA, we need a systematic way to decide how many principal components to use. So far, we’ve often used 1 or 2 principal components for the purposes of visualization, but in general, we need a more systematic approach. What if we start with an n×50n \times 50 matrix and want to decide how many principal components to compute?

Let’s define the total variance of an n×dn \times d matrix XX as the sum of the variances of the columns of XX. If xi(j)x_i^{(j)} is the ii-th value of the jj-th feature, and μj\mu_j is the mean of the jj-th feature, then the total variance of XX is

total variance=∑j=1d(1n∑i=1n(xi(j)−μj)2)\text{total variance} = \sum_{j = 1}^d \left(\frac{1}{n} \sum_{i = 1}^n (x_i^{(j)} - \mu_j)^2 \right)

Let’s compute this for X_three_features.

X_three_features.var(ddof=0) # The variances of the individual columns of X.
bill_depth_mm 3.866243 flipper_length_mm 195.851762 body_mass_g 646425.423171 dtype: float64
X_three_features.var(ddof=0).sum()
646625.1411755901

So, the total variance of X_three_features is approximately 646625.

This is also equal to the sum of the variances of the columns of X~V\tilde XV, i.e. the sum of the variances of the principal components!

# These three numbers are DIFFERENT than the numbers above,
# but their sum is the same.
(X_three_features_centered @ vt.T).var(ddof=0)
0 646575.552578 1 47.055776 2 2.532822 dtype: float64
(X_three_features_centered @ vt.T).var(ddof=0).sum() # Same sum as before.
646625.1411755901

Why? If we create the same number of principal components as we have original features, we haven’t lost any information – we’ve just written the same data in a different basis. The goal, though, is to pick a number of principal components that is relatively small (smaller than the number of original features), but still captures most of the variance in the data.

Each new principal component – that is, each column of X~V\tilde X V – captures some amount of this total variance. What we’d like to measure is the proportion (that is, fraction, or percentage) of the total variance that each principal component captures. The first principal component captures the most variance, since it corresponds to the direction in which the data varies the most. The second principal component captures the second-most variance, and so on. But how much variance does the first, second, or jj-th principal component capture?

Recall from earlier that the variance of the data projected onto a vector v⃗\vec v is given by

PV(v⃗)=1n∥X~v⃗∥2\text{PV}(\vec v) = \frac{1}{n} \lVert \tilde X \vec v \rVert^2

v⃗1\vec v_1 maximizes this quantity, which is why the first principal component is the projection of the data onto v⃗1\vec v_1; v⃗2\vec v_2 maximizes the quantity subject to being orthogonal to v⃗1\vec v_1, and so on.

Then, the variance of PC1\text{PC}_1 is whatever we get back when we plug v⃗1\vec v_1 into the formula above, and in general, the variance of PCj\text{PC}_j is PV(v⃗j)\text{PV}(\vec v_j) (again, where v⃗j\vec v_j is the jj-th column of VV in the SVD of X~\tilde X). Observe that if v⃗j\vec v_j is the jj-th column of VV in X~=UΣVT\tilde X = U \Sigma V^T, then,

PV(v⃗j)=1n∥X~v⃗j∥2=1n∥σju⃗j∥2=1nσj2∥u⃗j∥2=σj2n\text{PV}(\vec v_j) = \frac{1}{n} \lVert \tilde X \vec v_j \rVert^2 = \frac{1}{n} \lVert \sigma_j \vec u_j \rVert^2 = \frac{1}{n} \sigma_j^2 \lVert \vec u_j \rVert^2 = \frac{\sigma_j^2}{n}

Here, we used the ever important fact that X~v⃗j=σju⃗j\tilde X \vec v_j = \sigma_j \vec u_j, where σj\sigma_j is the jj-th singular value of X~\tilde X and u⃗j\vec u_j is the jj-th column of UU in the SVD of X~\tilde X.

What this tells us is that the variance of the jj-th principal component is σj2n\frac{\sigma_j^2}{n}. This is a beautiful result – it tells us that the variance of the jj-th principal component is simply the square of the jj-th singular value of X~\tilde X, divided by nn. The σj\sigma_j’s in the SVD represent the amount that the data varies in the direction of v⃗j\vec v_j! We don’t need any other fancy information to compute the variance of the principal components; we don’t need to know the individual principal component values, or have access to a variance method in code.

u, s, vt = np.linalg.svd(X_three_features_centered)
s
array([14673.43378383, 125.1781673 , 29.04185933])
np.set_printoptions(precision=8)
n = X_three_features_centered.shape[0]
variances_of_pcs = s ** 2 / n
np.set_printoptions(suppress=True)
variances_of_pcs
array([646575.55257751, 47.05577648, 2.5328216 ])

Above, we see the exact same values as if we computed the variances of the principal components directly from the data.

(X_three_features_centered @ vt.T).var(ddof=0)
0 646575.552578 1 47.055776 2 2.532822 dtype: float64

Since the total variance of XX is the sum of the variances of its principal components, the total variance of XX is then the sum of σj2n\frac{\sigma_j^2}{n} over all j=1,2,…,rj = 1, 2, \ldots, r, where r=rank(X)r = \text{rank}(X). (Remember that if k>rk > r, then σk=0\sigma_k = 0.)

total variance=∑j=1rσj2n\text{total variance} = \sum_{j = 1}^r \frac{\sigma_j^2}{n}

So, if we single out just one principal component, how much of the overall variation in XX does it capture? The answer is given by the proportion of variance explained by the jj-th principal component:

proportion of variance explained by PC j=σj2∑k=1rσk2\text{proportion of variance explained by PC } j = \frac{\sigma_j^2}{\sum_{k=1}^r \sigma_k^2}

This is a number between 0 and 1, which we can interpret as a percentage.

s # The singular values of X_three_features.
array([14673.43378383, 125.1781673 , 29.04185933])
(s ** 2) / (s ** 2).sum() # The proportions of variance explained by each principal component.
array([0.99992331, 0.00007277, 0.00000392])

The above tells us that in X_three_features, the first principal component captures 99.99% of the variance in the data! There’s very little information lost in projecting this 3-dimensional dataset into 1-dimensional space.

Often, the proportions above are visualized in a scree plot, as you’ll see in Homework 11. Scree plots allow us to visually decide the number of principal components to keep, based on where it seems like we’ve captured most of the variation in the data. We’ll work on a related example in live lecture.

The PCA Recipe

Let’s briefly summarize what I’ll call the “PCA recipe”, which describes how to find the principal components (new features) of a dataset.

  1. Starting with an n×dn \times d matrix XX of nn data points in dd dimensions, mean-center the data by subtracting the mean of each column from itself. The new matrix is X~\tilde X.

  2. Compute the singular value decomposition of X~\tilde X: X~=UΣVT\tilde X = U \Sigma V^T. The columns of VV (rows of VTV^T) describe the directions of maximal variance in the data.

  3. Principal component jj comes from multiplying X~\tilde X by the jj-th column of VV.

    PCj=X~v⃗j=σju⃗j\text{PC}_j = \tilde X \vec v_j = \sigma_j \vec u_j

    The resulting principal components – which are the columns of X~V\tilde X V – are uncorrelated with one another.

  4. The variance of PCj\text{PC}_j is σj2n\frac{\sigma_j^2}{n}. The proportion of variance explained by PCj\text{PC}_j is

    proportion of variance explained by PC j=σj2∑k=1rσk2\text{proportion of variance explained by PC } j = \frac{\sigma_j^2}{\sum_{k=1}^r \sigma_k^2}
EECS 245 Course Notes
5.3. Singular Value Decomposition
Chapter 0: Math Review
Chapter 0.1. Summation Notation and the Mean