09 - Principal Component Analysis (PCA) and Dimensionality Reduction

An 8x8 image of a handwritten digit has 64 numbers, a wine sample has 13, a gene expression profile has thousands, and almost always many of those numbers move together: two neighboring pixels are lit together, two chemical measurements grow together. Principal Component Analysis (PCA) looks for new axes, linear combinations of the original variables, ordered by how much of the data's variance (the dispersion around the mean) they capture and uncorrelated with one another.

Keeping the first ones and discarding the others compresses the data while losing, by construction, the least possible variance. Two mistakes make it useless or misleading: using it on variables with different units of measurement without standardizing, and mistaking the explained variance for the information useful to your problem. This article, the ninth in the ML Fundamentals series, computes it by hand on five points, reproduces it with the SVD and with scikit-learn, and measures the two mistakes on real data.

You need the concepts of clustering (the scikit-learn documentation suggests PCA before K-Means precisely when the dimensions are many, because Euclidean distances in high dimension get inflated) and a little linear algebra: eigenvalues and eigenvectors of a symmetric matrix. All the data are synthetic with a fixed seed or are datasets included in scikit-learn; every reported number comes from an executed listing.

PCA by hand on five points

Take five points in the plane already centered on the mean, $(-2,-1)$, $(-1,-1)$, $(0,0)$, $(1,1)$, $(2,1)$. The covariance matrix puts the variances of the two variables on the diagonal and their covariance off the diagonal: $S = X^\top X/(n-1)$. Here $S_{11} = 10/4 = 2.5$, $S_{22} = 4/4 = 1$ and $S_{12} = 6/4 = 1.5$.

The eigenvectors of $S$ are the directions of the principal components and the corresponding eigenvalues are the variances along those directions. For a 2x2 matrix they follow from the equation $\lambda^2 - \mathrm{tr}(S)\,\lambda + \det(S) = 0$, that is $\lambda^2 - 3.5\,\lambda + 0.25 = 0$, with solutions $\lambda = (3.5 \pm \sqrt{11.25})/2 = 3.4271$ and $0.0729$.

The listing redoes the computation in three ways: with the eigenvalues of the covariance, with the SVD (singular value decomposition) of the centered data, where the eigenvalues are $\sigma_i^2/(n-1)$, and with scikit-learn's PCA. Then it checks the three properties that define PCA: the variance along the first component is the first eigenvalue, the coordinates in the new axes are uncorrelated, and the total variance does not change, because a rotation neither creates nor destroys dispersion.

# PCA on five points: covariance, eigenvalues, SVD and sklearn give the same answer
import numpy as np
from sklearn.decomposition import PCA

X = np.array([[-2, -1], [-1, -1], [0, 0], [1, 1], [2, 1]], dtype=float)
n = len(X)
print("means:", X.mean(axis=0).tolist(), "(the data are already centered)")

S = X.T @ X / (n - 1)                                   # covariance matrix, divided by n-1
print("covariance:\n", S)
print("trace =", np.trace(S), " determinant =", round(float(np.linalg.det(S)), 4))

lam, V = np.linalg.eigh(S)                              # eigenvalues in ascending order
lam, V = lam[::-1], V[:, ::-1]
print("eigenvalues:", np.round(lam, 4).tolist())
print("eigenvector 1:", np.round(V[:, 0], 4).tolist(), " eigenvector 2:", np.round(V[:, 1], 4).tolist())

# same thing with the SVD of the centered data: eigenvalues = sigma^2 / (n-1)
U, s, Vt = np.linalg.svd(X, full_matrices=False)
print("sigma^2/(n-1):", np.round(s ** 2 / (n - 1), 4).tolist())

pca = PCA(n_components=2).fit(X)
print("sklearn explained_variance_:", np.round(pca.explained_variance_, 4).tolist())
print("sklearn explained_variance_ratio_:", np.round(pca.explained_variance_ratio_, 4).tolist())
print("sklearn components:", np.round(pca.components_, 4).tolist(), " (the sign may change)")

Z = X @ V                                               # coordinates in the new axes
print("variance along the 1st component:", round(float(Z[:, 0].var(ddof=1)), 4))
print("covariance between the two coordinates:", f"{float(np.cov(Z.T)[0, 1]):.1e}")
print("sum of the variances before and after:", float(np.trace(S)), round(float(Z.var(axis=0, ddof=1).sum()), 4))

# Output:
# means: [0.0, 0.0] (the data are already centered)
# covariance:
#  [[2.5 1.5]
#  [1.5 1. ]]
# trace = 3.5  determinant = 0.25
# eigenvalues: [3.4271, 0.0729]
# eigenvector 1: [-0.8507, -0.5257]  eigenvector 2: [0.5257, -0.8507]
# sigma^2/(n-1): [3.4271, 0.0729]
# sklearn explained_variance_: [3.4271, 0.0729]
# sklearn explained_variance_ratio_: [0.9792, 0.0208]
# sklearn components: [[0.8507, 0.5257], [-0.5257, 0.8507]]  (the sign may change)
# variance along the 1st component: 3.4271
# covariance between the two coordinates: 2.6e-16
# sum of the variances before and after: 3.5 3.5

The three methods agree: eigenvalues 3.4271 and 0.0729, equal to the hand computation. The first component explains 0.9792 of the total variance (3.4271 out of 3.5, the trace of the covariance); the second only 0.0208. The covariance between the new coordinates is 2.6e-16, that is zero up to rounding, and the sum of the variances is 3.5 before and after.

Note the sign of the eigenvector: eigh returns $(-0.8507, -0.5257)$ and scikit-learn $(0.8507, 0.5257)$. They are the same direction with the opposite orientation, and Jolliffe and Cadima warn that the coordinates can change sign with different software: never compare signs between two runs.

Why the SVD and not the covariance? The covariance is formed by multiplying $X^\top X$, a step that can lose precision, while the SVD works directly on the data. The scikit-learn documentation says that PCA centers the data but does not scale them, and that it chooses the solver according to the shape: covariance_eigh if there are fewer than 1000 features and more than ten times as many samples, randomized (Halko et al.'s truncated SVD) if the matrix exceeds 500x500 and you ask for less than 80% of the components, otherwise the full SVD. That «does not scale them» is the point of the next section.

Units of measurement decide the result

PCA looks for maximum variance, and variance depends on the unit: the same length in meters or in millimeters has variances that differ by a factor of a million. The wine dataset included in scikit-learn (178 samples, 13 chemical measurements, 3 wineries) shows it well: the variance of the 13 features goes from 0.015 to 98,610.

The listing computes PCA on the raw data and on the standardized data, and also checks that PCA on standardized data is PCA on the correlation matrix, the covariance of the variables with mean 0 and variance 1. Jolliffe and Cadima devote a section precisely to this choice between covariance and correlation.

# PCA on the wine dataset (178 samples, 13 chemical measurements): without and with standardization
import numpy as np
from sklearn.datasets import load_wine
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

wine = load_wine()
X, names = wine.data, np.array(wine.feature_names)
print("variance of the 13 features: min = %.3f, max = %.1f" % (X.var(axis=0).min(), X.var(axis=0).max()))

for label, M in [("raw", X), ("standardized", StandardScaler().fit_transform(X))]:
    p = PCA().fit(M)
    r = p.explained_variance_ratio_
    top = names[np.argsort(-np.abs(p.components_[0]))[:2]]
    print(f"{label:15s} explained variance PC1 = {r[0]:.3f}, PC2 = {r[1]:.3f}, "
          f"cumulative at 2 = {r[:2].sum():.3f}; PC1 weighs mostly: {top.tolist()}")

# PCA on standardized data = PCA on the correlation matrix
Z = StandardScaler().fit_transform(X)
corr = np.corrcoef(X, rowvar=False)
print("eigenvalues of the correlation, first 3:", np.round(np.linalg.eigvalsh(corr)[::-1][:3], 3).tolist())
# StandardScaler divides by n, PCA by n-1: the eigenvalues differ by the factor (n-1)/n = 177/178
print("explained_variance_ of PCA on standardized data (corrected), first 3:", np.round(PCA().fit(Z).explained_variance_[:3] * 177 / 178, 3).tolist())

# Output:
# variance of the 13 features: min = 0.015, max = 98609.6
# raw             explained variance PC1 = 0.998, PC2 = 0.002, cumulative at 2 = 1.000; PC1 weighs mostly: ['proline', 'magnesium']
# standardized    explained variance PC1 = 0.362, PC2 = 0.192, cumulative at 2 = 0.554; PC1 weighs mostly: ['flavanoids', 'total_phenols']
# eigenvalues of the correlation, first 3: [4.706, 2.497, 1.446]
# explained_variance_ of PCA on standardized data (corrected), first 3: [4.706, 2.497, 1.446]

On the raw data the first component explains 99.8% of the variance, but the largest weights are those of proline and magnesium: PCA discovered that one feature has large numbers, not that an important chemical direction exists. On the standardized data the first component drops to 36.2%, the second to 19.2%, and the weight goes to flavanoids and total_phenols, two measurements that are related to each other.

The eigenvalues of the correlation matrix (4.706; 2.497; 1.446) coincide with those of PCA on standardized data, after correcting the factor 177/178 that arises because StandardScaler divides by $n$ while PCA divides by $n-1$. The operational rule: if the variables have different units, standardize inside the Pipeline; if they all have the same unit (the pixels of an image), you can skip it.

How many components to keep

The ordered eigenvalues tell you how much variance each component carries. The explained variance ratio is $\lambda_i / \sum_j \lambda_j$, and its cumulative sum answers the question «keeping $k$ components, how much variance do I preserve?». The popular 90% or 95% rule is a convention, not a theorem, and its only justification is the variance.

If the data feed a classifier the right question is another: how much worse does the classifier get? The listing uses the 1,797 8x8 images of load_digits, with the pixels rescaled from 0..16 to 0..1 (they all have the same unit), and compares the cumulative variance with the accuracy of a logistic regression in 5-fold cross-validation, with PCA inside the pipeline and therefore recomputed on each training fold.

# How many components to keep: cumulative variance and downstream accuracy, on digits (1797 8x8 images = 64 pixels)
import numpy as np
from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import StratifiedKFold, cross_val_score
from sklearn.pipeline import make_pipeline

X, y = load_digits(return_X_y=True)
X = X / 16.0                       # pixels from 0..16 to 0..1: same unit for all, faster convergence
pca = PCA().fit(X)
cum = np.cumsum(pca.explained_variance_ratio_)
print("pixels with zero variance:", int((X.var(axis=0) == 0).sum()), " -> effective rank at most", 64 - int((X.var(axis=0) == 0).sum()))
for q in (0.80, 0.90, 0.95, 0.99):
    print(f"components to explain {q:.0%} of the variance: {int(np.searchsorted(cum, q) + 1)}")
print("PCA(n_components=0.95) chooses:", PCA(n_components=0.95).fit(X).n_components_)
print("first 5 eigenvalues:", np.round(pca.explained_variance_[:5], 3).tolist(), " numerically null eigenvalues (< 1e-10):", int((pca.explained_variance_ < 1e-10).sum()))

cv = StratifiedKFold(5, shuffle=True, random_state=0)
print("\n   k  cumulative variance  accuracy (5-fold CV, logistic regression)")
for k in (2, 5, 10, 20, 29, 40, 64):
    mod = make_pipeline(PCA(n_components=k), LogisticRegression(max_iter=1000))
    acc = cross_val_score(mod, X, y, cv=cv)
    print(f"{k:4d}  {cum[k - 1]:19.3f}  {acc.mean():.3f} (std {acc.std():.3f})")
without = cross_val_score(LogisticRegression(max_iter=1000), X, y, cv=cv)
print(f"without PCA (64 pixels): {without.mean():.3f} (std {without.std():.3f})")

# Output:
# pixels with zero variance: 3  -> effective rank at most 61
# components to explain 80% of the variance: 13
# components to explain 90% of the variance: 21
# components to explain 95% of the variance: 29
# components to explain 99% of the variance: 41
# PCA(n_components=0.95) chooses: 29
# first 5 eigenvalues: [0.699, 0.64, 0.554, 0.395, 0.272]  numerically null eigenvalues (< 1e-10): 3
#
#    k  cumulative variance  accuracy (5-fold CV, logistic regression)
#    2                0.285  0.605 (std 0.020)
#    5                0.545  0.850 (std 0.009)
#   10                0.738  0.928 (std 0.010)
#   20                0.894  0.953 (std 0.007)
#   29                0.955  0.962 (std 0.007)
#   40                0.988  0.968 (std 0.008)
#   64                1.000  0.969 (std 0.009)
# without PCA (64 pixels): 0.969 (std 0.009)

You need 21 components for 90% of the variance, 29 for 95%, 41 for 99%; three pixels have zero variance (they are always off), so the last three eigenvalues are null and the rank is at most 61. PCA(n_components=0.95) chooses 29 on its own. The accuracy does not follow the variance: with 2 components (28.5% of the variance) it is 0.605, with 10 (73.8%) it is 0.928, with 29 (95.5%) it is 0.962, against 0.969 without PCA.

The gain falls between 10 and 29 components and after 40 it stops. Here 95% of the variance costs less than one point of accuracy, but that is a property of these data and of this classifier, not of all: it is measured case by case, treating $k$ as a hyperparameter to be chosen in cross-validation.

Reconstruction and anomalies

Projecting onto $k$ components and going back with inverse_transform returns an approximate image, and the error has an exact form: the sum of the squared differences, divided by the total variance, equals $1 - \text{(explained variance)}$. This justifies reading PCA as compression: the unexplained variance is exactly the reconstruction error.

The same error, computed per image, can serve as an anomaly score: the images that the subspace of the digits cannot describe have a high error. The listing checks the identity for three values of $k$ and compares the error of real digits, not seen in training, with that of random noise images.

# Reconstruction: the error is exactly the variance you discarded; and it works as an anomaly score
import numpy as np
from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
from sklearn.metrics import roc_auc_score
from sklearn.model_selection import train_test_split

X, y = load_digits(return_X_y=True)
Xtr, Xte = train_test_split(X, test_size=0.3, random_state=0)

print("   k  relative error = 1 - explained variance   (measured on the training set)")
for k in (2, 10, 30):
    p = PCA(n_components=k).fit(Xtr)
    R = p.inverse_transform(p.transform(Xtr))
    err = ((Xtr - R) ** 2).sum() / ((Xtr - Xtr.mean(axis=0)) ** 2).sum()
    print(f"{k:4d}  error {err:.4f}   1 - ratio {1 - p.explained_variance_ratio_.sum():.4f}")

# Anomaly score: reconstruction error of each image with k = 20
p = PCA(n_components=20).fit(Xtr)
def error(M):
    return ((M - p.inverse_transform(p.transform(M))) ** 2).sum(axis=1)

rng = np.random.default_rng(0)
noise = rng.integers(0, 17, size=(len(Xte), 64)).astype(float)    # random pixels 0..16, like the digits
e_ok, e_noise = error(Xte), error(noise)
auc = roc_auc_score(np.r_[np.zeros(len(e_ok)), np.ones(len(e_noise))], np.r_[e_ok, e_noise])
print(f"median error of real digits = {np.median(e_ok):.1f}, noise images = {np.median(e_noise):.1f}")
print(f"AUC of the reconstruction error as an anomaly score: {auc:.3f}")

# Output:
#    k  relative error = 1 - explained variance   (measured on the training set)
#    2  error 0.7073   1 - ratio 0.7073
#   10  error 0.2573   1 - ratio 0.2573
#   30  error 0.0403   1 - ratio 0.0403
# median error of real digits = 124.3, noise images = 2322.9
# AUC of the reconstruction error as an anomaly score: 1.000

The identity holds to the fourth decimal digit: 0.7073 with two components, 0.2573 with ten, 0.0403 with thirty, equal to $1 - \text{variance ratio}$. With 20 components the median error of the real digits is 124 and that of the noise 2,323, and the AUC is 1.000.

Do not be fooled by the perfection: uniform noise is a very easy anomaly, because it has nothing of the structure of the strokes of a digit. For subtler anomalies (a badly written digit, a digit from another handwriting) the AUC would be lower, and it must be measured on real examples of your problem.

When linear PCA is not enough

PCA finds only straight axes. For two concentric circles, which no line separates, the two eigenvalues are almost equal (0.503 and 0.497): there is no privileged direction. Kernel PCA applies the same method to an implicit transformation of the variables, with the kernel trick seen in the SVM, and according to the scikit-learn documentation it obtains a non-linear reduction.

The listing uses the RBF kernel with three values of gamma and for each one measures how well a line (logistic regression in cross-validation) separates the two circles using the first $m$ components, from 1 to 4.

# When linear PCA is not enough: two concentric circles (inner circle at 0.4 of the outer radius)
import numpy as np
from sklearn.datasets import make_circles
from sklearn.decomposition import PCA, KernelPCA
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

X, y = make_circles(n_samples=400, noise=0.05, factor=0.4, random_state=0)
print("linear PCA, variance explained by the two components:", np.round(PCA(2).fit(X).explained_variance_ratio_, 3).tolist())

def acc(transform):
    """Accuracy of a logistic regression (that is, a LINE) on the components, in 5-fold CV."""
    return cross_val_score(make_pipeline(transform, StandardScaler(), LogisticRegression()), X, y, cv=5).mean()

print(f"line on the original data (PCA with 2 components): {acc(PCA(2)):.3f}")
print("\nKernelPCA RBF: accuracy of the line on the first m components")
print("gamma      m=1    m=2    m=3    m=4")
for g in (1.0, 3.0, 10.0):
    r = [acc(KernelPCA(m, kernel="rbf", gamma=g)) for m in (1, 2, 3, 4)]
    print(f"{g:5.1f}  " + "  ".join(f"{v:6.3f}" for v in r))

# Output:
# linear PCA, variance explained by the two components: [0.503, 0.497]
# line on the original data (PCA with 2 components): 0.435
#
# KernelPCA RBF: accuracy of the line on the first m components
# gamma      m=1    m=2    m=3    m=4
#   1.0   0.458   0.450   1.000   1.000
#   3.0   0.492   0.505   1.000   1.000
#  10.0   0.665   0.700   1.000   1.000

The line on the original data gives 0.435, below the 50% of chance: there is nothing linear to learn. With Kernel PCA the first two components are not enough (0.450-0.700 depending on gamma), and with three components the accuracy is 1.000 for all three values of gamma. The lesson is twofold: in kernel methods the number of components to keep cannot be deduced from the idea that the first one separates everything, and gamma changes the quality of the first two. You try it, as with any hyperparameter.

t-SNE: for drawing, not for measuring

To visualize high-dimensional data in two dimensions t-SNE (van der Maaten and Hinton, 2008) is often used, which builds a map in which points that are close in the original space stay close and which, according to the authors, reduces the tendency to crowd the points at the center. Unlike PCA it is not a linear projection, its cost function is not convex (with different initializations you get different maps, the documentation says) and it has no transform method for new points.

To measure the fidelity of the drawing I use scikit-learn's trustworthiness, which is 1 if the ten neighbors of every point in the drawing are also neighbors in the original space, and I compute it for the 2-component PCA and for three t-SNE runs with different seeds on load_digits.

# t-SNE against PCA for drawing digits in 2D: trustworthiness (1 = the neighbors in the drawing are truly neighbors)
import time
import numpy as np
from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE, trustworthiness

X, y = load_digits(return_X_y=True)
print("PCA with 2 components: trustworthiness =", round(trustworthiness(X, PCA(2).fit_transform(X), n_neighbors=10), 3))
T = []
for seed in (0, 1, 2):
    Z = TSNE(n_components=2, perplexity=30, init="pca", random_state=seed).fit_transform(X)
    T.append(trustworthiness(X, Z, n_neighbors=10))
print("t-SNE in 2 dimensions, three seeds: trustworthiness =", np.round(T, 3).tolist())
print("does TSNE have a transform method for new points?", hasattr(TSNE(), "transform"))
print("does PCA have a transform method for new points?", hasattr(PCA(), "transform"))

# Output:
# PCA with 2 components: trustworthiness = 0.83
# t-SNE in 2 dimensions, three seeds: trustworthiness = [0.993, 0.993, 0.993]
# does TSNE have a transform method for new points? False
# does PCA have a transform method for new points? True

The two-component PCA has trustworthiness 0.830, t-SNE 0.993 with all three seeds: for local proximity alone, t-SNE draws much better. The price is that distances between groups that are far apart on the map have no quantitative meaning, and that new points cannot be projected, so t-SNE is not a step to put in a prediction pipeline: it serves to look at the data. PCA, with transform, remains the tool for decorrelating, compressing and feeding a model.

Worked exercise

Three points are aligned on a line: $(-2,-1)$, $(0,0)$, $(2,1)$. Before running anything, compute the covariance matrix, the eigenvalues and the direction of the first component. What do you expect from the second component? Then compare with the listing.

# Exercise: three collinear points
import numpy as np
from sklearn.decomposition import PCA

X = np.array([[-2.0, -1.0], [0.0, 0.0], [2.0, 1.0]])
S = X.T @ X / (len(X) - 1)
print("covariance:", S.tolist())
lam, V = np.linalg.eigh(S)
print("eigenvalues:", np.round(lam[::-1], 6).tolist())
print("principal direction:", np.round(np.abs(V[:, -1]), 4).tolist(), " 2/sqrt(5), 1/sqrt(5) =", round(2 / 5 ** 0.5, 4), round(1 / 5 ** 0.5, 4))
print("variance explained by the first component:", round(float(PCA(1).fit(X).explained_variance_ratio_[0]), 6))

# Output:
# covariance: [[4.0, 2.0], [2.0, 1.0]]
# eigenvalues: [5.0, 0.0]
# principal direction: [0.8944, 0.4472]  2/sqrt(5), 1/sqrt(5) = 0.8944 0.4472
# variance explained by the first component: 1.0

What you can do now

After this article you can do five things, each with a numerical check ready:

• compute PCA by hand from covariance, eigenvalues and eigenvectors, and reproduce it with the SVD;
• decide whether to standardize, looking at the units of measurement and not at habit;
• choose $k$ by comparing cumulative variance and downstream accuracy in cross-validation;
• read the reconstruction error as discarded variance and use it as an anomaly score, with caution;
• recognize when a non-linear method is needed, and when a t-SNE plot is only a drawing.

The five abilities share a common thread: PCA is a rotation and a choice of which axes to throw away, and every doubt about it is settled by computing how much variance, and of which variable, you are keeping. After this article three stages of the series remain: feedforward neural networks, evaluation metrics and deployment in production. For PCA written from scratch and t-SNE, the Machine Learning Engineer series has a dedicated article.