Principal Component Analysis

Take 5,000 houses with two features: size in square meters and number of rooms. The two columns are highly correlated — big houses have more rooms — so a scatter plot is a fat ellipse tilted at 45 degrees. PCA asks a single question: which direction through that cloud has the most variance, and what is the next direction orthogonal to it? Those two directions are the principal components, and in 2D the answer is visibly the ellipse's long axis and its short axis.

Here is the full mechanics, because interviewers routinely ask you to derive it. You have a data matrix X, n samples by p features, and you center every column first (subtract the mean — PCA measures variance, and uncentered data would make it measure distance from the origin instead). Define the sample covariance matrix:

S = X_c^T X_c / (n - 1)          # p x p, symmetric, positive semi-definite

You want the unit vector w that maximizes the variance of the projection w^T x, which is w^T S w subject to ||w|| = 1. Lagrange:

L(w, lambda) = w^T S w - lambda (w^T w - 1)
dL/dw = 2 S w - 2 lambda w = 0   =>   S w = lambda w

Substituting back, the variance along w is lambda. So the first principal component is the eigenvector of S with the largest eigenvalue, the second is the eigenvector with the next largest, and so on. PCA is an eigendecomposition of the covariance matrix, with eigenvectors sorted by eigenvalue. That is the whole algorithm.

Three structural facts fall out of this and are worth memorizing because they generate follow-up questions:

• Eigenvectors of a symmetric matrix are orthogonal, so the components are mutually uncorrelated: Cov(PC_i, PC_j) = w_i^T S w_j = lambda_j w_i^T w_j = 0 for i != j. Feeding PCs to a linear model removes multicollinearity by construction.
• The variance of PC_i is lambda_i, and the sum of all eigenvalues equals the trace of S — the total variance. Every component claims an exact share of the variance budget.
• Each component is a linear combination of the original features: PC1 = w_11x_1 + ... + w_1px_p. The weights are called loadings, they are only interpretable if your features are on comparable scales, and their signs are arbitrary (an eigenvector and its negation are both valid).

The numpy version, on correlated data:

import numpy as np

rng = np.random.default_rng(0)
n = 5000
size = rng.normal(120, 30, n)                 # m^2
rooms = 1.2 + size / 40 + rng.normal(0, 0.8, n)
X = np.column_stack([size, rooms])

Xc = X - X.mean(axis=0)
S = Xc.T @ Xc / (n - 1)

vals, vecs = np.linalg.eigh(S)                # symmetric solver -> real, orthogonal
order = np.argsort(vals)[::-1]
vals, vecs = vals[order], vecs[:, order]
print(vals)            # e.g. [901.3, 0.64]  — one axis dominates
print(vecs[:, 0])      # e.g. [0.9987, 0.0513] — nearly pure "size"

X_pca = Xc @ vecs                           # project onto both components

Note eigh, not eig — it exploits symmetry, guaranteeing real eigenvalues and an orthogonal basis. The first eigenvalue dwarfs the second because the data lies almost along one line; PC1 is essentially the size column rescaled, and PC2 is the residual jitter.

Under the hood, scikit-learn never forms S explicitly. It runs a singular value decomposition on the centered matrix, X_c = U Sigma V^T, where the columns of V are exactly these eigenvectors and lambda_i = sigma_i^2 / (n - 1). That is numerically stabler and cheaper when n < p — a case where the p x p covariance matrix would be rank-deficient but the thin SVD still works. So when someone asks "what is the relationship between PCA and SVD?", the precise answer is: PCA on centered data is the SVD of the centered data matrix; the components live in V and the explained variances are the squared singular values scaled by 1/(n-1).

The price of all this is interpretability: a component mixes your original features, so "PC1 rose by 2.3" is not a sentence a product manager can act on. Keep that cost in mind — it is the main reason Feature Selection Methods exist as a parallel path, and the reason the next chapter is about when the variance PCA preserves is actually the variance you care about.

You ran PCA on a 64-feature dataset and now face the practical question every pipeline eventually asks: how many components do you keep? The default answer in scikit-learn is a one-liner — PCA(n_components=0.95) — but what it actually computes, and when that criterion misleads, is where interviews separate candidates.

The quantity is the explained variance ratio: EVR_i = lambda_i / sum(lambda). Because the eigenvalues are sorted, the cumulative sum of EVR tells you what fraction of total variance the first k components capture:

import numpy as np
from sklearn.decomposition import PCA
from sklearn.datasets import load_digits
from sklearn.preprocessing import StandardScaler

X, y = load_digits(return_X_y=True)
Xs = StandardScaler().fit_transform(X)        # 64 pixel features, different scales

pca = PCA().fit(Xs)                           # keep all components for inspection
evr = pca.explained_variance_ratio_
cum = np.cumsum(evr)
for k in [1, 5, 10, 20, 30]:
    print(f"k={k:2d}  cumulative variance = {cum[k-1]:.3f}")
k95 = np.searchsorted(cum, 0.95) + 1
print("components for 95% variance:", k95)    # ~29 of 64

On the digits dataset, roughly 29 of 64 components hold 95% of the variance — a 2.2x compression with little visual loss. Setting n_components=0.95 automates exactly this: it picks the smallest k whose cumulative EVR reaches the threshold.

The same numbers govern reconstruction error, which is a nice way to prove you understand what PCA discarded. Reconstruct with x_hat = x_mean + W_k W_k^T (x - x_mean). The mean squared reconstruction error equals the variance thrown away:

pca20 = PCA(n_components=20).fit(Xs)
X20 = pca20.transform(Xs)
X_hat = pca20.inverse_transform(X20)
mse = ((Xs - X_hat) ** 2).sum(axis=1).mean()
discarded = pca20.explained_variance_[20:].sum()   # eigenvalues of dropped components
print(mse, discarded / (Xs.shape[1] - 1))           # match up to the n/(n-1) normalization

Dropping component i costs exactly lambda_i of reconstruction error, per feature, on average. That identity is why the cumulative-EVR rule is principled for compression tasks: you discard the directions that contribute least to describing your data.

Now the trap. EVR measures variance, not usefulness for the target, because PCA never looks at y. Construct two classes that differ by a small shift along one axis while a noisy axis carries ten times the variance, and PCA's first component will faithfully preserve the noise and may push the class-separating direction to component three or four. Keeping 95% variance can then discard the label signal. This is not hypothetical — it is the standard critique of PCA-before-classifier, and the correct rebuttal is always empirical: select k by cross-validated downstream performance, not by EVR alone.

from sklearn.pipeline import Pipeline
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score

pipe = Pipeline([
    ("scaler", StandardScaler()),
    ("pca", PCA(n_components=0.95)),
    ("clf", LogisticRegression(max_iter=2000)),
])
print(cross_val_score(pipe, X, y, cv=5).mean())
# grid over PCA(n_components=k) for k in [5, 10, 20, 40] to compare against EVR's choice

The other prerequisite is scaling, and it is the single most common PCA bug in industry code. PCA maximizes variance, so the feature with the largest units wins: express size in millimeters instead of meters and its variance multiplies by 10^6, making PC1 that column with nothing learned. Two equivalent fixes:

Default to standardizing unless all features share physical units and you deliberately want high-variance features to dominate. The one scenario where you skip both this worry and PCA entirely is when features are already decorrelated and regularized — that decision belongs to Feature Selection Methods, which follows this topic. The mechanism behind all of it, though, is the covariance eigendecomposition covered in this topic — if you cannot re-derive S w = lambda w from memory, do that before the next interview, and revisit why any of this matters in Curse of Dimensionality.

A candidate writes the classic leaky pipeline in a coding interview:

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score
import numpy as np

Xs = StandardScaler().fit_transform(X)          # fit on ALL data
Z = PCA(n_components=50).fit_transform(Xs)      # fit on ALL data
print(cross_val_score(LogisticRegression(max_iter=2000), Z, y, cv=5).mean())

The score is inflated, sometimes by several points, and the reason is subtle enough that even experienced engineers miss it. PCA and the scaler are unsupervised, so people assume they cannot leak. They can: the components are directions chosen to maximize variance over the entire dataset, which includes every test fold. Your model then receives features that were optimized with knowledge of the test data's spread. The fix is to make preprocessing a first-class part of the model — fitted inside each fold:

from sklearn.pipeline import Pipeline
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score

pipe = Pipeline([
    ("scaler", StandardScaler()),
    ("pca", PCA(n_components=50)),
    ("clf", LogisticRegression(max_iter=2000)),
])
print(cross_val_score(pipe, X, y, cv=5).mean())   # honest estimate

The same rule extends to every hyperparameter of PCA: if you tune n_components by cross-validated score, and then report that same score, you have selection bias — the honest protocol is tuning inside the CV loop (or nested CV) and reporting on a held-out test set.

A knob worth knowing is whiten=True. By default, PC1 has variance lambda_1, PC2 has lambda_2, and so on — the first components dominate any distance- or scale-sensitive downstream model. Whitening divides each component by sqrt(lambda_i), giving every retained component unit variance. It helps kernel SVMs and kNN on embeddings, and it changes regularization geometry: with whitening, a ridge penalty shrinks all directions equally instead of penalizing the already-small ones hardest. The cost is real, though — whitening amplifies the low-variance directions you kept, including residual noise in them.

When does PCA earn its place in a pipeline? Use it when p is large and collinear: pixel features, correlated sensor arrays, text or sentence embeddings, genomics — anywhere the covariance matrix is near-singular and a regularized linear model is fighting redundant columns. PCA decorrelates features, cuts training time (50 components instead of 784 features), and often denoises, because isotropic noise spreads its energy evenly across all directions, so the tiny eigenvalues are disproportionately noise.

Skip it when you must explain predictions in original feature units (components mix features — loadings are a partial patch), when features are already few and regularized models handle them fine, or when a low-variance direction might carry the label — EVR cannot see the target.

Cost-wise, scikit-learn's PCA runs an SVD of the centered matrix at O(min(n^2 p, n p^2)). Above ~10,000 features, svd_solver="randomized" approximates the top components with a truncated SVD (randomized range finder) for a large speedup at slight accuracy cost; for datasets too big for memory, IncrementalPCA fits in mini-batches.

The honest summary of PCA's ceiling: it is a linear projection. Unfold a swiss roll — a 2D sheet rolled into 3D — and no choice of linear components can flatten it; the intrinsic coordinates are a nonlinear function of the inputs. Nonlinear structure is the domain of t-SNE and UMAP later in this roadmap, and if your actual goal is keeping interpretable original features rather than compressed combinations, the right tool is the next topic, Feature Selection Methods. Revisit the failure modes of skipping all of this in Curse of Dimensionality if the motivation ever feels abstract.