Principal component analysis

Examples

Directions of largest variance

The first principal component is the direction in which the projected data has the largest variance.

Given data \(x_{ij}\) with column means equal to zero, the score of point \(i\) on the first component is

\[ z_{i1} = \langle \phi_1, x_i \rangle = \sum_{j=1}^{p} x_{ij}\phi_{j1}. \]

We look for the loadings \(\phi_{11}, \dots, \phi_{p1}\) that maximize \(\frac{1}{n}\sum_i z_{i1}^2\) subject to \(\sum_j \phi_{j1}^2 = 1\).

The \(k\)-th principal component is the direction of largest variance among all directions orthogonal to the first \(k-1\).

Three dimensional data, shown as two of its axes at a time. Each line runs from a point to the plane through the origin perpendicular to the chosen component, so its length is that point’s score. The first component has the longest lines.

The solution is that \(\phi_k\) is the eigenvector of \(\hat\Sigma\) with the \(k\)-th largest eigenvalue.

Turn the direction and watch the variance of the projection. The red segments are the distance from each point to its projection. The first principal component is the angle where the variance is largest, which is also where those segments are shortest.

In matrix notation

With all components at once,

\[ Z = X\Phi, \]

where \(\Phi\) is \(p \times p\) and its columns are the eigenvectors of \(\hat\Sigma\). The eigenvalue \(\lambda_k\) is the variance along component \(k\).

This is closely related to the singular value decomposition \(X = U S V^\top\), where \(U\) and \(V\) are orthogonal and \(S\) is diagonal. One can show that

\[ \Phi = V \qquad\text{and}\qquad Z = US. \]

In practice we compute the SVD of \(X\) rather than the eigenvectors of \(\hat\Sigma\). It is more accurate, and it avoids forming a \(p \times p\) matrix when \(p\) is large.

A change of basis

The clearest way to think about PCA is as a rotation of the coordinate system.

		size
data	\(X\)	\(n \times p\)	rows are coordinates in the standard basis
loadings	\(\Phi\)	\(p \times p\)	columns are the new basis vectors
scores	\(Z = X\Phi\)	\(n \times p\)	rows are coordinates in the new basis
reconstruction	\(X = Z\Phi^\top\)	\(n \times p\)	back to the standard basis

So far nothing is lost. This is a lossless change of basis.

The reduction comes from keeping only the first \(L\) columns,

\[ Z_L = X\Phi_L \in \mathbb{R}^{n \times L}, \qquad X_L = Z_L \Phi_L^\top \in \mathbb{R}^{n \times p}. \]

Now something is lost. One can show that \(\Phi_L\) is the choice that minimizes \(\|X - X\Phi_L\Phi_L^\top\|_2^2\), so PCA gives the \(L\) dimensional linear subspace that is closest to the data.

These are two different descriptions of the same thing. The directions of largest variance are also the subspace closest to the data.

An example

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

rng = np.random.default_rng(12)

cov = [[3.0, 2.2], [2.2, 2.0]]
X = rng.multivariate_normal([0, 0], cov, 300)
X = X - X.mean(axis=0)

pca = PCA().fit(X)
print("loadings (columns are components):")
print(pca.components_.T.round(3))
print("variance along each component:", pca.explained_variance_.round(3))

1

PCA needs centred data. sklearn centres internally, but it is worth being explicit.

2

sklearn stores the loadings as rows of components_, so we transpose to get our \(\Phi\).

loadings (columns are components):
[[ 0.781 -0.625]
 [ 0.625  0.781]]
variance along each component: [4.009 0.267]

Z = pca.transform(X)

fig, axes = plt.subplots(1, 2, figsize=(6.6, 3.2))
axes[0].plot(X[:, 0], X[:, 1], "o", ms=3, alpha=0.4)
for k in range(2):
    v = pca.components_[k] * np.sqrt(pca.explained_variance_[k]) * 2
    axes[0].annotate("", xy=v, xytext=(0, 0),
                     arrowprops=dict(arrowstyle="->", lw=2, color="#a8452f"))
axes[0].set_title("original coordinates", fontsize=9)
axes[1].plot(Z[:, 0], Z[:, 1], "o", ms=3, alpha=0.4)
axes[1].set_title("scores", fontsize=9)
for ax in axes:
    ax.set_aspect("equal")
    ax.set(xticks=[], yticks=[])
plt.show()

Figure 60.1: The two principal directions, scaled by their standard deviation.

In the new coordinates the data is uncorrelated. PCA decorrelates.

print("correlation before:", round(float(np.corrcoef(X.T)[0, 1]), 4))
print("correlation after: ", round(float(np.corrcoef(Z.T)[0, 1]), 4))

correlation before: 0.87
correlation after:  0.0

How many components

p = 20
A = rng.normal(0, 1, (p, 4))
D = rng.normal(0, 1, (400, 4)) @ A.T + rng.normal(0, 0.4, (400, p))

full = PCA().fit(D)
ratio = full.explained_variance_ratio_

fig, axes = plt.subplots(1, 2, figsize=(6.6, 2.8))
axes[0].plot(np.arange(1, p + 1), ratio, "o-")
axes[0].set(xlabel="component", ylabel="proportion of variance")
axes[1].plot(np.arange(1, p + 1), np.cumsum(ratio), "o-")
axes[1].axhline(0.95, ls="--", lw=1.2, color="#a8452f")
axes[1].set(xlabel="component", ylabel="cumulative proportion")
plt.show()

Figure 60.2: A scree plot and the cumulative proportion of variance explained.

The proportion of variance explained by component \(k\) is

\[ \frac{\lambda_k}{\sum_{j=1}^{p}\lambda_j} = \frac{s_k^2}{\sum_j s_j^2}. \]

There is an elbow after four components, which is how the data was generated. The elbow method is to look for the point where the curve flattens. It is a judgement, not a rule.

Scaling matters

PCA maximizes variance, and variance depends on units. A predictor measured in millimetres has a variance a million times larger than the same predictor in metres, and it will dominate the first component.

from sklearn.preprocessing import StandardScaler

W = rng.normal(0, 1, (200, 3))
W[:, 2] *= 1000

raw = PCA().fit(W)
scaled = PCA().fit(StandardScaler().fit_transform(W))

print("without scaling:", raw.explained_variance_ratio_.round(4))
print("with scaling:   ", scaled.explained_variance_ratio_.round(4))

1

The third column is measured in different units.

without scaling: [1. 0. 0.]
with scaling:    [0.3838 0.3383 0.278 ]

Without scaling, the first component is almost entirely the third column.

Standardize before PCA, unless all columns are already in the same units and the differences in variance are meaningful.

Biplots

A biplot shows the scores and the loadings in the same figure. The points are the observations in the new coordinates, and the arrows show how each original variable relates to the components.

from sklearn.datasets import load_wine

wine = load_wine()
Xw = StandardScaler().fit_transform(wine.data)
pw = PCA(n_components=2).fit(Xw)
Zw = pw.transform(Xw)

fig, ax = plt.subplots(figsize=(5.4, 4.6))
ax.plot(Zw[:, 0], Zw[:, 1], "o", ms=3, alpha=0.35, color="#7a838b")
for i, name in enumerate(wine.feature_names):
    v = pw.components_[:, i] * 3.2
    ax.annotate("", xy=v, xytext=(0, 0),
                arrowprops=dict(arrowstyle="->", lw=1, color="#a8452f"))
    ax.text(v[0] * 1.12, v[1] * 1.12, name, fontsize=6, ha="center", color="#a8452f")
ax.set(xlabel="PC 1 score", ylabel="PC 2 score")
plt.show()

Figure 60.3: A biplot. Arrows that point the same way belong to variables that are correlated.

Two arrows that point in the same direction belong to variables that are correlated. An arrow at right angles to another means the two are uncorrelated. A long arrow means the variable is well represented in these two components.