# 🧭 PCA via SVD: Explained Variance and Component Selection

By Better Wise

You already know that a matrix has a singular value decomposition. This lesson is where that fact starts earning its keep: PCA is the first unsupervised model in this course, and it is the reducer your capstone will use to shrink features before the k-means and MLP stages. There is one discipline I want you to carry over from the gradient work in Module 1: no code until the quantity you are computing is written down on paper. Here the quantity is variance, and the whole lesson is about reading it off a singular value spectrum correctly.

**After this lesson you can:**
- Center a data matrix and factor it with `np.linalg.svd` instead of eigendecomposing a covariance matrix.
- Convert singular values into explained variance and explained variance ratios.
- Read a scree plot and a cumulative curve, and defend a component count with numbers.
- Project, reconstruct, and quantify exactly what you threw away.
- Decide when to standardize and when centering alone is the honest choice.

## Center First, Then Factor

PCA looks for the orthonormal directions that maximize projected variance. Run it on an uncentered matrix and the first direction points at the mean vector, because the squared norm of `X` includes `n * ||mean||^2`. That is a location, not a spread. One subtraction fixes it:

```python
import numpy as np

def pca_fit(X, n_components=None):
    X = np.asarray(X, dtype=float)
    n, d = X.shape
    mu = X.mean(axis=0)
    Xc = X - mu
    U, s, Vt = np.linalg.svd(Xc, full_matrices=False)  # thin SVD
    var = s**2 / (n - 1)
    k = len(s) if n_components is None else n_components
    return {
        "mean": mu,
        "components": Vt[:k],           # (k, d), rows are the axes
        "singular_values": s[:k],
        "explained_variance": var[:k],
        "explained_variance_ratio": var[:k] / var.sum(),
    }

def pca_transform(model, X):
    return (np.asarray(X, dtype=float) - model["mean"]) @ model["components"].T
```

`full_matrices=False` gives you the thin factors, which is all you need when `n` and `d` differ by orders of magnitude. Never form `Xc.T @ Xc` and call `eigh` on it: squaring the matrix squares the condition number, and the small components, the ones you actually want to discard deliberately, get destroyed by rounding first.

## From Singular Values to Explained Variance

If `Xc = U S Vt`, then the rows of `Vt` are the principal directions, `U @ diag(s)` are the scores (the coordinates in the new basis), and the sample covariance is `Xc.T @ Xc / (n - 1) = V diag(s**2 / (n - 1)) Vt`. That identity is the whole lesson:

$$
\text{variance}_i = \frac{s_i^2}{n-1}, \qquad \text{ratio}_i = \frac{s_i^2}{\sum_j s_j^2}
$$

Variance has units, ratios do not, and the ratios always sum to something at most 1. Grab the ratios for your project matrix, and note that `np.var(pca_transform(m, X)[:, 0], ddof=1)` reproduces `m["explained_variance"][0]` up to floating point. If it does not, you skipped the centering.

The variance along the i-th principal component equals _____.
The explained variance ratio of component i is _____.
The rows of _____ hold the principal directions, and the scores are _____.
The truncated SVD of rank k is the best rank-k approximation in the _____ by the Eckart-Young theorem.
Total squared reconstruction error of a rank-k truncation equals _____.

**Exercise 1.** Your project matrix has shape `(500, 30)` and you fit without subtracting the mean. What does the first component look like, and how would you catch the bug from the spectrum alone?
**Answer**

The first right singular vector tilts toward the mean vector, and `s[0]` jumps because the total squared norm now contains the between-mean energy. Refit centered and compare `s[:3]` side by side: an uncentered fit typically inflates `s[0]` and shrinks every subsequent ratio. A cheap assertion is `np.allclose(pca_transform(m, X).mean(axis=0), 0, atol=1e-10)`.

## Choosing k from the Spectrum

Two plots, thirty seconds, one decision. The scree plot shows `explained_variance` per component in order; the elbow is where the drop flattens. The cumulative plot shows `np.cumsum(ratio)` and is what you quote in a report.

| Rule | Signal it reads | Where it fails |
| --- | --- | --- |
| Cumulative threshold (95%, 99%) | total variance retained | ignores which directions the downstream model needs |
| Elbow on the scree plot | curvature of the drop-off | subjective on smooth spectra |
| Kaiser (eigenvalue > 1) | variance above one standardized unit | only valid on the correlation matrix, so it requires standardization |

**PCA is unsupervised, so variance is not the same as usefulness.** A low-variance direction can carry the class-separation signal your logistic regression needs, and a high-variance direction can be pure measurement noise. Use the 95% rule when you need a defensible default, then confirm k by downstream validation rather than trusting it.

**Exercise 2.** You observe `s = [10, 4, 2, 1, 1]` with `n = 101`. What are the explained variance of component 1 and the cumulative ratio at `k = 2`?
**Answer**

`sum(s**2) = 100 + 16 + 4 + 1 + 1 = 122`. Component 1 variance is `100 / 100 = 1.0`, and the cumulative ratio at `k = 2` is `(100 + 16) / 122 = 0.951`. Two components already clear a 95% threshold, but the spectrum says nothing about whether those two components are the ones your classifier wants.

## Projection, Reconstruction, and What You Lose

`Z = Xc @ Vt[:k].T` is the projection; `Z @ Vt[:k] + mu` is the reconstruction, and the Eckart-Young theorem guarantees no other rank-k matrix gets closer in Frobenius norm. The cost is explicit:

```python
s_discarded = s[k:]
rel_error = np.sqrt((s_discarded**2).sum() / (s**2).sum())
```

Report `rel_error` next to `k` every time, because it converts "we kept 95% of the variance" into "we accept 22% relative reconstruction error". Two practical footnotes. Standardize columns before the fit when features have different units, but know that this redefines variance as correlation and can amplify near-constant columns. And singular vectors are unique only up to sign: `U[:, i]` and `Vt[i]` can both flip with no change in the decomposition.

> **Important:** 
> Never compare two PCA fits by raw component values. Compare by absolute cosine similarity, or adopt scikit-learn's `svd_flip` convention (make the largest-magnitude entry of each `Vt` row positive) so your capstone table is reproducible across runs.

**Exercise 3.** Out of 30 components, your cumulative curve crosses 95% at `k = 12`, but a downstream logistic regression on 12 components scores worse than on 6. What do you ship?
**Answer**

Ship `k = 6`, and report both numbers. The 95% rule is a statement about retained variance, not about discriminative content; components 7 through 12 are low-variance directions that spend degrees of freedom on label noise. Then repeat the comparison against the untransformed baseline and against `sklearn.decomposition.PCA`, which is exactly the protocol your capstone needs.

## Numeric Habits That Save You a Silent Bug

- Center inside the fold. If you compute the mean on the full dataset before cross-validation, you leak information into the test split.
- Fit PCA on train only, then call `pca_transform` on test with the stored mean. Recomputing `mu` at test time is the second most common leak.
- Assert `ratio.sum() <= 1` and that ratios are non-increasing. Both fail loudly when you hand the function the wrong axis.
- Keep `singular_values` in the model object. You cannot recover `rel_error` later without them.
- For wide matrices (`d >> n`), prefer the thin SVD or a randomized solver over the full one; the full `U` costs `O(n * d)` memory for nothing.

Front: Why must you center before running PCA?
Back: Otherwise the first component points toward the mean vector, which measures location rather than spread.

Front: In Xc = U S Vt, what are the principal directions?
Back: The rows of Vt, the right singular vectors.

Front: What are the PCA scores?
Back: U scaled columnwise by s, equivalently Xc @ Vt.T.

Front: Variance along component i in terms of singular values?
Back: s_i squared divided by (n - 1).

Front: Explained variance ratio formula?
Back: s_i squared divided by the sum of all squared singular values.

Front: Why use the SVD of Xc instead of eigh of Xc.T @ Xc?
Back: Forming the covariance squares the condition number and loses small components to rounding.

Front: What does the Eckart-Young theorem guarantee?
Back: The rank-k truncated SVD is the best rank-k approximation in Frobenius and spectral norm.

Front: Squared Frobenius reconstruction error of a rank-k truncation?
Back: The sum of the squared discarded singular values.

Front: What does full_matrices=False return?
Back: The thin SVD: U is n by min(n,d) and Vt is min(n,d) by d, saving memory.

Front: When should you standardize rather than only center?
Back: When features have different units or scales, so raw variances are not comparable.

Front: The Kaiser rule (eigenvalue > 1) applies to which matrix?
Back: The correlation matrix, that is, only after standardizing to unit variance.

Front: Why do two PCA runs return components with flipped signs?
Back: Singular vectors are unique only up to sign; U column and Vt row flip together, which svd_flip resolves.

Front: What does a scree plot show?
Back: Explained variance per component in descending order, used to look for an elbow.

## What Carries Over to the Reduction Stage

Centering is not optional, singular values squared over `n - 1` is the variance, and the ratio is that quantity divided by the total. SVD on the centered matrix beats eigendecomposition of the covariance on both stability and speed, and the discarded singular values tell you exactly what the truncation costs. Choose k with a rule, then justify it with a downstream number.

Next up is k-means: the same reduced matrix you just built becomes the input space, and the distance computations there are the ones PCA was quietly optimizing for you.

---

## Contents

Part of [🧮 Mathematics for Machine Learning 4: Build Core ML Algorithms from Scratch](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch.md)

- **📈 Linear and Cross-Entropy Models from Scratch**
  - [📉 Linear Regression with Normal Equations and Gradient Descent](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/linear-regression-normal-equations-and-gradient-descent.md)
  - [🔍 Gradient Checking with a Finite-Difference Harness](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/gradient-checking-finite-difference-harness.md) (sign-in required)
  - [🎚️ Logistic Regression to Softmax: Cross-Entropy Done Right](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/logistic-regression-to-softmax-cross-entropy.md)
  - [✅ Checkpoint: Linear and Cross-Entropy Models](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/checkpoint-linear-and-cross-entropy-models.md)
- **🔵 Unsupervised Learning from Scratch**
  - [🧭 PCA via SVD: Explained Variance and Component Selection](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/pca-via-svd-explained-variance.md) — this lesson
  - [🌀 k-means: Lloyd's Algorithm and k-means++ Initialization](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/k-means-lloyd-and-k-means-plus-plus.md)
  - [🔔 Gaussian Mixture Models with Expectation Maximization](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/gaussian-mixture-models-expectation-maximization.md)
  - [✅ Checkpoint: PCA, k-means and GMM](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/checkpoint-pca-kmeans-and-gmm.md)
- **🧠 Backpropagation and Benchmark Capstone**
  - [🔗 Backpropagation for a One-Hidden-Layer MLP](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/backpropagation-for-one-hidden-layer-mlp.md)
  - [📊 Cross-Validation and a Fair Baseline Comparison Protocol](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/cross-validation-and-baseline-comparison-protocol.md)
  - [🏁 Capstone: NumPy Implementations vs Library Baselines](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/capstone-numpy-vs-library-baselines.md) (premium content)
  - [✅ Checkpoint: Backpropagation and Model Comparison](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/checkpoint-backpropagation-and-model-comparison.md) (premium content)
- [Course cheatsheet](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/cheatsheet.md) (premium content)