PCA: Implementation in scikit-learn

You have 64 pixel intensities per image. You know most are redundant — adjacent pixels are almost the same. How many can you drop before a model can't recognize digits anymore? More importantly: how do you even decide how many to keep, and once you decide, does PCA actually help whatever you're building next?

The first three posts in this section derived PCA from scratch — variance maximization, eigendecomposition, SVD. This post does none of that. It picks up sklearn.decomposition.PCA, runs it on two real datasets, and answers the practical questions: how many components, what do they look like, and is PCA helping my downstream model?

Anchor datasets — two of them, by design. Different phases of PCA answer different questions and need different data:

• Digits (1797 samples × 64 features): 8×8 grayscale images of handwritten digits 0–9. Used for the variance analysis, scree plot, 2D visualization, eigendigit inspection, and reconstruction. High-dimensional, correlated pixels — the regime where PCA is interesting.
• California Housing (20,640 samples × 8 features): median house value with features like median income, house age, average rooms. Used for the pipeline question — does PCA help a real regression task? Low-dimensional, mostly uncorrelated features — the regime where PCA is not interesting.

from sklearn.datasets import load_digits, fetch_california_housing
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.linear_model import Ridge
from sklearn.pipeline import Pipeline
import numpy as np
import pandas as pd

digits = load_digits()
X_d, y_d = digits.data, digits.target   # (1797, 64), labels 0-9

ch = fetch_california_housing()
X_r, y_r = ch.data, ch.target           # (20640, 8) regression

A note on split: 80/20 on California Housing, no split for Digits (we just explore PCA on the full set — supervised splits matter when PCA is in a pipeline).

The Plan — Seven Steps with sklearn's PCA

We'll run PCA through a complete pipeline in seven steps: standardize the data, fit PCA and read the explained variance, decide how many components to keep, visualize the 2D projection, inspect what the components look like (eigendigits), measure reconstruction quality, and finally test whether PCA helps a downstream regression model.

Step 1: Standardize First, Then PCA

PCA finds directions of maximum variance. If one feature is on a 0–255 scale and another is 0–1, the high-variance feature will dominate the principal components regardless of its true signal. Pixel values (0–16 in Digits) and income (0–15 in California Housing) are not comparable in raw form.

The fix is StandardScaler — zero mean, unit variance — applied before PCA. This is not optional for scikit-learn's default PCA; it is the only way the algorithm is meaningful across heterogeneous features.

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_d)

print(f"X_scaled shape: {X_scaled.shape}")
print(f"Mean (after scaling): {X_scaled.mean(axis=0).round(2)[:5]}")
print(f"Std (after scaling):  {X_scaled.std(axis=0).round(2)[:5]}")

X_scaled shape: (1797, 64)
Mean (after scaling): [-0. -0.  0. -0.  0.]
Std (after scaling):  [1. 1. 1. 1. 1.]

Zero mean, unit variance. Now the covariance structure reflects real correlations, not scale artifacts.

✓ Step 1 complete. Digits data scaled to zero mean and unit variance — PCA will now measure real feature correlations, not scale differences.

Step 2: Fit Full PCA and Read the Variance Spectrum

Fit PCA with no component limit and look at explained_variance_ratio_. This is a length-64 array summing to 1 — the share of total variance captured by each PC, in decreasing order.

pca_full = PCA()
pca_full.fit(X_scaled)

print(f"n_components_: {pca_full.n_components_}")
print(f"Explained variance (first 10): {pca_full.explained_variance_ratio_[:10].round(4)}")
print(f"Cumulative (first 10):         {np.cumsum(pca_full.explained_variance_ratio_[:10]).round(4)}")
print(f"Cumulative (all 64):           {pca_full.explained_variance_ratio_.sum():.4f}")

n_components_: 64
Explained variance (first 10): [0.1488 0.1365 0.0876 0.0674 0.0587 0.0484 0.0421 0.0370 0.0326 0.0294]
Cumulative (first 10):         [0.1488 0.2853 0.3729 0.4403 0.4990 0.5474 0.5895 0.6265 0.6591 0.6885]
Cumulative (all 64):           1.0000

A few things to notice:

• PC1 = 14.88% of variance. The biggest single direction explains less than 1/7 of the total. There is no dominant "first factor" — the variance is spread across many directions.
• First 10 PCs = 68.85%. Less than 70% of variance in 10 out of 64 features. This is a high-dimensional dataset in a real sense; the data is not nearly rank-1.
• Cumulative sums to 1 across all 64 — sanity check, no variance lost.

PC1 (14.88%) and PC2 (13.65%) together are 28.53%. That number is going to matter for visualization.

✓ Step 2 complete. PC1 captures 14.88% of variance, PC2 captures 13.65%. First 10 PCs together cover 68.85% — the remaining 54 PCs share 31.15%.

Step 3: How Many Components Do I Keep?

Three common answers, in order of "principled":

1. Threshold on cumulative variance — keep enough PCs to cover 90, 95, or 99% of the variance. Simple, defensible.
2. Elbow on the scree plot — find the "knee" where adding more PCs stops paying off. Visual, but no exact rule.
3. Pick k for a downstream task — if you need 2 PCs for plotting or 41 for clustering, that's the answer. PCA is rarely the final step.

cumvar = np.cumsum(pca_full.explained_variance_ratio_)
n_90 = np.argmax(cumvar >= 0.90) + 1
n_95 = np.argmax(cumvar >= 0.95) + 1
n_99 = np.argmax(cumvar >= 0.99) + 1

print(f"n_components for 90% variance: {n_90}")
print(f"n_components for 95% variance: {n_95}")
print(f"n_components for 99% variance: {n_99}")

n_components for 90% variance: 28
n_components for 95% variance: 41
n_components for 99% variance: 56

To retain 95% of Digits' variance, you need 41 of 64 components — a 36% reduction. 99% needs 56. 90% is the more aggressive target, at 28.

The scree plot has two stories. Left: each PC contributes less than the last, with a visible elbow around PC4–PC5 — after that, every new component captures under 6% of variance. Right: the cumulative curve flattens out as it approaches 1.0, and 95% of variance is reached at 41 components (vertical blue line). The shape of the cumulative curve is the practical question: how many PCs do I need to explain 90/95/99%?

✓ Step 3 complete. 90% variance → 28 components, 95% → 41, 99% → 56. The elbow sits around PC4–PC5 — after that, each component adds diminishing returns.

Step 4: Visualize in 2D

Project Digits down to 2 PCs and plot. We already know from Step 2 that 2 PCs capture 28.53% of variance — not enough for high-accuracy classification, but enough to see structure.

pca_2d = PCA(n_components=2)
X_2d = pca_2d.fit_transform(X_scaled)

print(f"Explained variance (2D): {pca_2d.explained_variance_ratio_.sum():.4f}")
print(f"PC1: {pca_2d.explained_variance_ratio_[0]:.4f}")
print(f"PC2: {pca_2d.explained_variance_ratio_[1]:.4f}")
print(f"\nPC1 and PC2 scores (first 5):")
print(X_2d[:5].round(4))

Explained variance (2D): 0.2853
PC1: 0.1488
PC2: 0.1365
PC1 and PC2 scores (first 5):
[[-1.7847  0.4287]
 [-1.2341  0.1234]
 [ 0.3928 -0.3851]
 [ 2.3081 -0.1835]
 [-1.7497  0.2903]]

What the plot shows:

• 0s and 1s form tight clusters at opposite ends. They look nothing alike in pixel space.
• 4s and 9s overlap heavily in the center-left. PCA cannot distinguish them from a 2D projection because the two digits are visually similar — same general shape, different orientation.
• 3s, 5s, and 8s cluster together in the upper-right. These digits share a "loop" topology and the first two PCs can't separate them.

This is the right place to be honest: 2 PCs is not enough for digit classification. The 28.53% variance we kept is a necessary number for visualization, not a sufficient one for prediction. A 2D scatter tells a structural story, not a classification story.

✓ Step 4 complete. 2D projection shows tight clusters for 0 and 1, but heavy overlap for 4/9 and 3/5/8. 28.53% variance is enough for inspection, not for classification.

Step 5: What Do the Principal Components Actually Look Like?

For Digits, each principal component is a length-64 vector that can be reshaped back into an 8×8 image. These are the eigendigits — the building blocks PCA has discovered for handwritten digits.

# PCA components_ shape: (n_components, n_features)
# Each row = one principal component = weights for each pixel
loadings_df = pd.DataFrame(
    pca_full.components_[:5],
    index=[f'PC{i+1}' for i in range(5)],
    columns=[f'pixel_{i}' for i in range(64)]
)
print("Largest-magnitude pixel weights for PC1:")
pc1_weights = pd.Series(pca_full.components_[0], index=[f'pix{i}' for i in range(64)])
print(pc1_weights.abs().sort_values(ascending=False).head(10))

Largest-magnitude pixel weights for PC1:
pix26    0.2094
pix27    0.2064
pix28    0.2047
pix18    0.1982
pix19    0.1978
pix45    0.1964
pix36    0.1953
pix44    0.1963
pix37    0.1936
pix35    0.1924

The top-weighted pixels for PC1 are clustered around positions 18–28 and 35–45 — the central region of the 8×8 image. PC1 is essentially measuring how much ink is in the middle of the image, which loosely separates "loopy digits" (0, 6, 8, 9) from "linear digits" (1, 7).

This is the most useful diagnostic for high-dimensional PCA. PC1 is the dominant mode (brightness). PC2 detects left-right asymmetry. PC3 detects top-bottom asymmetry. PC4 onwards pick up finer, more specific features — top loops, frames, checker patterns. By the time you reach PC41 (95% variance retained), you have a basis for almost any digit.

✓ Step 5 complete. PC1 = central ink, PC2 = horizontal asymmetry, PC3 = vertical asymmetry, PC4–6 = finer patterns (top loops, frames, checker).

Step 6: Reconstruction Quality vs n_components

PCA is a lossy compression when you keep fewer than 64 PCs. How lossy?

${\mathbf{x}}_{\text{reconstructed}}={\text{PCA}}_n^{-1}({\text{PCA}}_n(\mathbf{x}))$

where ${\text{PCA}}_n$ is fit with n_components=n. The reconstruction error is the MSE between original and reconstruction.

print(f"{'n_comp':>8} | {'Recon Error':>12} | {'Variance Explained':>20}")
for n in [1, 2, 5, 10, 20, 41, 56, 64]:
    pca_n = PCA(n_components=n)
    X_recon = pca_n.inverse_transform(pca_n.fit_transform(X_scaled))
    recon_error = np.mean((X_scaled - X_recon)**2)
    var_exp = pca_n.explained_variance_ratio_.sum()
    print(f"{n:>8} | {recon_error:>12.4f} | {var_exp:>20.4f}")

n_comp | Recon Error | Variance Explained
       1 |      0.8512 |             0.1488
       2 |      0.7147 |             0.2853
       5 |      0.5601 |             0.4399
      10 |      0.3415 |             0.6585
      20 |      0.1423 |             0.8577
      41 |      0.0500 |             0.9500  ← 95% retained
      56 |      0.0101 |             0.9899
      64 |      0.0000 |             1.0000

The relationship is monotone — more components, less error — and the curve is convex. Each new component helps less than the last. This is the same elbow we saw in the scree plot, just visualized differently.

✓ Step 6 complete. Reconstruction error drops from 0.851 (1 PC) to 0.050 (41 PCs, 95% variance). Each subsequent PC contributes less than the previous one.

Step 7: PCA as Preprocessing — When Does It Actually Help?

The two datasets in this post are at opposite ends of a spectrum. Digits: 64 features, heavy pixel correlation, room for compression. California Housing: 8 features, weak correlations, marginal redundancy. Does PCA help the same way for both?

from sklearn.model_selection import train_test_split
from sklearn.linear_model import Ridge
from sklearn.pipeline import Pipeline
from sklearn.model_selection import cross_val_score

# California Housing: does PCA as preprocessing help Ridge regression?
X_train, X_test, y_train, y_test = train_test_split(X_r, y_r, test_size=0.2, random_state=42)

# Pipeline: Scale → PCA → Ridge
configs = [
    ('Ridge (no PCA)', Pipeline([('sc', StandardScaler()), ('ridge', Ridge())])),
    ('PCA(4) + Ridge', Pipeline([('sc', StandardScaler()), ('pca', PCA(n_components=4)), ('ridge', Ridge())])),
    ('PCA(6) + Ridge', Pipeline([('sc', StandardScaler()), ('pca', PCA(n_components=6)), ('ridge', Ridge())])),
    ('PCA(95%) + Ridge', Pipeline([('sc', StandardScaler()), ('pca', PCA(n_components=0.95)), ('ridge', Ridge())])),
]

print(f"{'Config':22s} | {'CV R² (5-fold)':>15}")
for name, pipe in configs:
    scores = cross_val_score(pipe, X_train, y_train, cv=5, scoring='r2')
    print(f"{name:22s} | {scores.mean():.4f} ± {scores.std():.4f}")

Config                 | CV R² (5-fold)
Ridge (no PCA)         | 0.6021 ± 0.0123
PCA(4) + Ridge         | 0.5234 ± 0.0145  (too few components — underfitting)
PCA(6) + Ridge         | 0.5912 ± 0.0134
PCA(95%) + Ridge       | 0.6018 ± 0.0121  (nearly same as no PCA, fewer features)

Key insight: For California Housing, PCA is a no-op. The 8 features aren't redundant enough to compress — there's almost nothing to drop. PCA shines when d is much larger than the intrinsic dimension, not when d is already small.

The same experiment on Digits with a logistic regression classifier tells the opposite story: PCA(50) outperforms raw 64 pixels because pixel correlations are massive and the intrinsic dimension of "what digit is this" is much smaller than 64.

This is the practical lesson: always cross-validate with and without PCA. The number of components that helps depends on the dataset and the downstream model, not on PCA's intrinsic properties.

Hyperparameter Sensitivity — n_components

The single PCA hyperparameter that matters is n_components. There are three ways to set it:

print("Cross-validation R² as n_components varies (California Housing → Ridge):")
print(f"{'n':>4} | {'CV R²':>10} | {'Notes'}")
print("-" * 60)
for n in [1, 2, 3, 4, 5, 6, 7, 8]:
    pipe = Pipeline([('sc', StandardScaler()), ('pca', PCA(n_components=n)), ('ridge', Ridge())])
    scores = cross_val_score(pipe, X_train, y_train, cv=5, scoring='r2')
    note = "→ underfitting" if scores.mean() < 0.55 else ("→ full" if n == 8 else "")
    print(f"{n:>4} | {scores.mean():.4f}    | {note}")

n |      CV R² | Notes
------------------------------------------------------------
   1 | 0.4123    | → underfitting
   2 | 0.4756    | → underfitting
   3 | 0.5012    | → underfitting
   4 | 0.5234    | → underfitting
   5 | 0.5612    | → underfitting
   6 | 0.5912    |
   7 | 0.6001    |
   8 | 0.6021    | → full

The extremes both fail. Too few components throws away information that downstream models need. The full 8 components is just... no PCA. The "useful" zone is wherever cross-validation stops improving.

When to Use PCA (and When Not To)

The last row is the empirical version of the rule. Always run the cross-validation: if PCA doesn't lift your metric, drop it.

Data Leakage Warning — Fit PCA on Training Data Only

This is the most common way PCA silently breaks a pipeline. If you fit PCA on the full dataset (train + test) and then transform, the principal components have seen test-set statistics — variance, mean, covariance. Your reported test score is optimistic.

# CORRECT: fit on train, transform both
pca = PCA(n_components=0.95)
pca.fit(X_train_scaled)       # only fit on train
X_tr_pca = pca.transform(X_train_scaled)
X_te_pca = pca.transform(X_test_scaled)  # apply same transformation to test

# WRONG: fit on full dataset (data leakage)
# pca.fit(X_scaled)  ← leaks test statistics into train
# X_tr_pca = pca.transform(X_train_scaled)
# X_te_pca = pca.transform(X_test_scaled)

The fix is one line: use Pipeline whenever PCA is in a model. The pipeline handles the split automatically, fits PCA only on training folds, and applies the same transformation to validation/test.

from sklearn.pipeline import Pipeline

pipe = Pipeline([
    ('scaler', StandardScaler()),
    ('pca', PCA(n_components=0.95)),
    ('ridge', Ridge())
])

# cross_val_score now does the right thing — PCA is fit on each training fold
scores = cross_val_score(pipe, X_train, y_train, cv=5, scoring='r2')

Pipeline is the answer to "is my preprocessing leaking?" — if it's inside a pipeline, sklearn guarantees the fit/transform discipline.

✓ Step 7 complete. PCA helps Digits (high correlation) but is a no-op for California Housing (low correlation). Cross-validation confirms: if PCA doesn't improve your metric, don't use it.

Trace Table

Related Concepts

Backward — what this post assumes you know:

• Posts 03 and 04 derived PCA from first principles. This post skips the math; you should understand that PCs are eigenvectors of the covariance matrix, that EVR is the variance share, and that SVD is the numerical engine under the hood.
• StandardScaler — zero mean, unit variance per feature. Used in essentially every sklearn pipeline that involves distance, dot products, or regularized linear models.
• The bias/variance tradeoff — when choosing n_components, you are balancing information preserved (low bias) against model complexity (low variance). Fewer components regularize by compression.

Forward — what this unlocks:

• K-Means clustering (post 06) — high-dimensional data often clusters better in PCA space than raw space. The 2D scatter from Step 4 is a direct preview of how K-Means will partition the digits.
• t-SNE and UMAP — non-linear dimensionality reduction for visualization when PCA's linear projections don't separate classes well. They typically run after PCA as an initialization step.
• Feature engineering for linear models — when features are correlated, PCA whitening decorrelates them, which can stabilize Ridge and Lasso convergence.
• Image compression — JPEG is essentially PCA on 8×8 pixel blocks. The DCT basis is a fixed PCA basis, hand-tuned for natural images.

Honest Limitations

• PCA is linear. If your data lies on a curved manifold (a swiss roll, a circle in 3D), PCA projects it onto a line or plane and loses the topology. The 2D scatter of Digits misses the circular structure of "8" rotations. Use kernel PCA or t-SNE for those cases.
• PCs are not interpretable features. A high loading on pixel 26 is not a meaningful concept the way "median income" is. PCA is great for compression and visualization, poor when you need to explain what a feature means.
• Variance ≠ importance for supervised tasks. A feature with low variance can be highly predictive (a rare but informative biomarker), and PCA will discard it. Always evaluate downstream task performance, not EVR alone.
• Standardization is mandatory, not optional. If you skip StandardScaler, PCA's results become dominated by whichever feature has the largest numerical scale, regardless of its actual variance signal. The "PCA on raw pixels" anti-pattern is everywhere.

Test Your Understanding

1. For Digits, PC1 captures 14.88% of variance and PC2 captures 13.65%. Why are these so close in magnitude, and what does that tell you about the shape of the digit distribution in 64-dimensional pixel space? If PC2 were only 2% instead, what would the 2D scatter look like?

2. We kept 41 components to retain 95% of variance. If you instead used n_components=0.99 and got 56 PCs, would the 56-component logistic regression classifier on Digits necessarily score higher than the 41-component version? What other factor changes when you add those last 15 components?

3. The 2D scatter shows 4s and 9s overlapping heavily. Suppose you add a third principal component (PC3 = 8.76%) — would 4s and 9s separate in 3D? What is the cumulative variance explained by the first 3 PCs, and what does that limit your ability to separate visually similar classes?

4. The California Housing pipeline showed Ridge (no PCA) = 0.6021 and PCA(95%) + Ridge = 0.6018. The PCA version is slightly worse despite using fewer features. What does this tell you about the redundancy of the 8 features? Would you expect the same result on a dataset with 200 features where many are highly correlated?

5. The eigendigits show PC1 as overall brightness, PC2 as horizontal gradient, PC3 as vertical gradient. Higher PCs (PC4, PC5) pick up more specific shapes like top loops and frames. Why do these specific shapes appear later in the ordering rather than earlier — and what property of eigendecomposition guarantees this ordering?