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:

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.

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:

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.