Principal Component Analysis

The previous chapter ended with an observation: four housing features — sqft, rooms, bedrooms, and baths — inter-correlate at r > 0.7 because they all measure roughly the same underlying concept, “house size”. The features are redundant but the concept is real. Ideally we would have one feature called size that combines all four, keeps most of what they collectively carry, and makes the downstream regression’s coefficients stable. Principal component analysis — PCA — is the tool that constructs that size feature from the correlation structure alone, with no domain knowledge required.

The geometric picture

PCA is most easily understood as a rotation of the coordinate axes. Take the scatter of sqft vs rooms. Because the two features correlate at r = 0.85, the cloud is stretched along a diagonal. The direction of maximum spread — what we would call the long axis of the cloud — is not aligned with either the sqft axis or the rooms axis. It sits at roughly 45° between them.

If we define two new features — PC1, the position along the diagonal, and PC2, the perpendicular distance to the diagonal — then PC1 carries almost all the spread in the cloud and PC2 carries almost none. Specifically: with r = 0.85, PC1’s variance (its eigenvalue) is 1.85 in standardised units and PC2’s is 0.15. Of the total variance in the dataset (sum of eigenvalues = 2.0), PC1 explains 92.5% and PC2 the remaining 7.5%.

This is the entire pedagogical core of PCA. In two dimensions it is just a 45° rotation. In ten dimensions it is a rotation in 10-dimensional space to a new basis whose first axis captures the most variance of any possible axis, the second captures the most variance among directions perpendicular to the first, and so on. Every axis in the new basis is uncorrelated with every other — the correlation matrix of the transformed features is the identity plus zero everywhere else.

The algebra

The directions of maximum variance are exactly the eigenvectors of the data’s covariance matrix. For a d-feature dataset, the covariance matrix is d × d: the entry at row i, column j is the covariance between feature i and feature j. The diagonal entries are the variances of each feature alone. If the features have been standardised (mean 0, variance 1), the covariance matrix is identical to the correlation matrix we looked at in the previous chapter.

For the sqft/rooms pair with correlation 0.85, the covariance matrix is:

The eigenvectors of this 2×2 matrix turn out to be exactly (1, 1)/√2 and (1, −1)/√2 — the 45° and −45° directions — and the eigenvalues are 1 + 0.85 = 1.85 and 1 − 0.85 = 0.15. This is not a coincidence: for any symmetric matrix with equal diagonal entries, the eigenvectors are always the ±45° directions with eigenvalues equal to the diagonal plus and minus the off-diagonal term.

More generally, the matrices appearing in PCA are always symmetric (covariance matrices are symmetric by construction — Cov(X, Y) = Cov(Y, X)) and positive-semi-definite (variances are never negative). Symmetric positive-semi-definite matrices have real, non-negative eigenvalues and orthogonal eigenvectors, so PCA always produces real perpendicular principal components. This is why PCA is computationally stable and why it exists as a well-defined procedure at all: not every matrix has such nice structure.

Compression: keeping only the top k components

The real operational value of PCA comes from keeping only the first k principal components — the ones with the largest eigenvalues. Projecting the data onto those k directions gives a lower-dimensional representation that retains as much variance as possible. The variance retained is the sum of the first k eigenvalues divided by the sum of all eigenvalues.

Consider a real housing dataset with 10 features. After standardising each feature and running PCA, the eigenvalues might come out as 5.70, 2.20, 0.80, 0.55, 0.35, 0.18, 0.12, 0.06, 0.03, 0.01 — a total of 10.0 (as expected when every feature contributes variance 1 after standardisation).

The scree plot (left panel) visualises the rapid decay of eigenvalues. The cumulative curve (right panel) reads off how much variance you retain for any choice of k. For this dataset:

• k = 3: 87% of variance retained. A good default choice — past the elbow, below the 95% threshold, uses only 30% of the original feature count.
• k = 5: 96% retained. More conservative; appropriate if 4% variance loss is too much for the downstream task.
• k = 7: 99.1% retained. Nearly lossless but wastes most of PCA’s compression benefit.

A common rule: choose k by the variance threshold (80% or 95% depending on task sensitivity), then verify by checking downstream performance. If the model’s validation accuracy is equal at k=3 and k=5, pick the smaller one. If k=3 costs meaningful accuracy, use k=5. The scree plot is the reason to have a short candidate list of k values, not the final word on which one to use.

A worked example: five houses

Take five houses from a dataset, standardised and described only by sqft and rooms:

PC1 score is the projection onto PC1: x · PC1 = (sqft_z + rooms_z)/√2. PC2 is (sqft_z − rooms_z)/√2. The PC1 scores neatly range from −1.84 to +2.12, separating houses by size. The PC2 scores range only from −0.21 to +0.14 — a tenth the spread of PC1, reflecting that 92.5% of the variance lives along PC1.

The price class visibly tracks PC1 position in both panels. In the original space, a regression on sqft and rooms would have had the correlation problem we saw in the last chapter. After projection onto PC1 alone, the regression becomes a simple univariate fit — price as a function of the size index, with no correlation to worry about.

The preprocessing order matters

PCA fits on training data to learn the rotation, then applies that rotation to validation and test data. The same pipeline discipline from the Data Splits chapter applies: fitting PCA on the full dataset before splitting is a form of data leakage — the rotation has been informed by validation and test rows.

import numpy as np
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split

# X: 10,000 house feature rows with 10 columns (sqft, rooms, bedrooms, baths,
# lot_size, garages, floors, age, year_built, crime_rate)
# y: house price in GBP

X_trainval, X_test, y_trainval, y_test = train_test_split(
  X, y, test_size=0.15, random_state=42,
)

# The pipeline: standardise -> PCA -> regression
# PCA is fitted INSIDE the pipeline, so it only sees X_trainval when calling .fit
pipe = Pipeline([
  ('scale', StandardScaler()),
  ('pca',   PCA(n_components=3)),
  ('reg',   LinearRegression()),
])

pipe.fit(X_trainval, y_trainval)

# Inspect what PCA learned
pca_step = pipe.named_steps['pca']
print("Variance ratio per PC:", pca_step.explained_variance_ratio_.round(3))
# Output (example): [0.57 0.22 0.08]
print(f"Cumulative at k=3:     {pca_step.explained_variance_ratio_.sum():.1%}")
# Output: Cumulative at k=3:     87.0%

# Evaluate on test
print("R2 on test:", pipe.score(X_test, y_test))

# View the PC1 loadings — what the 'size' component weights each original feature
pc1_loadings = pca_step.components_[0]
feature_names = ['sqft', 'rooms', 'bedrooms', 'baths', 'lot_size',
               'garages', 'floors', 'age', 'year_built', 'crime_rate']
for f, w in sorted(zip(feature_names, pc1_loadings), key=lambda x: -abs(x[1])):
  print(f"  {f:14s} {w:+.3f}")
# Output (example):
#   sqft           +0.420
#   rooms          +0.405
#   bedrooms       +0.385
#   baths          +0.350
#   lot_size       +0.220
#   garages        +0.190
#   floors         +0.150
#   age            -0.130
#   year_built     +0.125
#   crime_rate     -0.050

The PC1 loadings (last block) confirm the interpretation: PC1 is a positively-weighted combination of all the size-related features, with smaller coefficients on the less-correlated features. It is a ‘size index’ that summarises what the four redundant size features collectively contained — in one column instead of four.

Common pitfalls

Three failure modes show up repeatedly in practice:

Unscaled features. If one column has variance 10,000 (dollars) and another has variance 1 (z-scored feature), PCA will essentially recover the dollar feature as PC1 and everything else as small corrections. Always standardise before PCA, or use a Pipeline that does it automatically.

Leakage from fitting on the full dataset. Fitting PCA on train + test together means the rotation was informed by test rows. Subtle, but test performance becomes optimistic. The Pipeline above is leakage-free because pipe.fit(X_trainval, ...) never sees X_test.

Over-interpreting loadings. A PC1 with positive weights on every feature is usually a ‘general intensity’ component. A PC2 that contrasts old-construction features against new-construction features might look like an ‘era’ component but is really just a mathematical artifact of the covariance structure. Interpretations are hypotheses to verify with domain knowledge, not causal truths.

PCA is a linear method and only captures linear correlation structure. For data that lies on a curved manifold — a non-linear lower-dimensional structure — PCA will miss the structure. The next chapter compares PCA against non-linear alternatives (t-SNE, UMAP) and the supervised alternative (LDA), developing the judgement for when each is the right tool.