PCA, in depth: variance maximisation and reconstruction error, a worked ten-point example, SVD done right, choosing k, scaling to big data and the failure modes

One idea, two derivations

Start with centred data: subtract each column's mean so the cloud sits at the origin. Projecting a centred row x onto a unit vector w gives the scalar w·x. Ask for the w that makes the variance of those scalars as large as possible. That variance is w^TCw, where C = X^TX/(n-1) is the sample covariance matrix, and maximising it under the constraint |w| = 1 with a Lagrange multiplier gives Cw = λw. The best direction is the eigenvector of C with the largest eigenvalue, and the eigenvalue is the variance along it. Repeat under the constraint of being orthogonal to the earlier directions and you get the remaining eigenvectors in order.

There is a second derivation that explains why PCA is a compression method. Project every row onto a k-dimensional subspace and measure the total squared distance between each row and its projection. The subspace that minimises that reconstruction error is spanned by the same top-k eigenvectors, and the error left over equals (n-1) times the sum of the discarded eigenvalues. Maximum variance kept and minimum error lost are one problem seen from two sides, which is why the explained-variance ratio λ₁+…+λ_k over the sum of all eigenvalues is the natural dial for choosing k.

Two consequences follow straight from the algebra. PCA is linear: every component is a weighted sum of the original columns, so it can only find flat structure. And it is driven entirely by variance, so a column measured in millimetres will dominate one measured in metres for no reason other than its units.

A worked example by hand

Take the ten points plotted in the figure below. The column means are 1.81 and 1.91. After centring, the sample covariance matrix is [[0.6166, 0.6154], [0.6154, 0.7166]]: both columns vary by similar amounts and move strongly together.

For a 2×2 symmetric matrix the eigenvalues come from the trace and determinant: λ = tr/2 ± sqrt(tr²/4 - det). That gives λ₁ = 1.2840 and λ₂ = 0.0491. The first component carries 96.3% of the total variance. Its eigenvector is (0.678, 0.735), a line tilted slightly steeper than 45 degrees, and the second is perpendicular to it.

Keeping k = 1 turns each point into one number, its score along PC1. The first point (2.5, 2.4) scores 0.828 and the second (0.5, 0.7) scores -1.778. Mapping the scores back (mean plus score times PC1) lands every point on the orange line; the dashed segments are what was thrown away. Their squared lengths sum to 0.4418, which is exactly (n-1)λ₂ = 9 × 0.0491, as the reconstruction argument promised.

Computing PCA correctly

The textbook recipe builds C and calls an eigensolver. In code, prefer the singular value decomposition of the centred data matrix itself: X_c = UΣV^T. The rows of V^T are the principal directions, the squared singular values divided by n-1 are the eigenvalues, and UΣ holds the scores. Forming X^TX squares the condition number, so directions with small variance lose about twice as many significant digits as they need to; the SVD route avoids that. The eigen route is still fine, and faster, when n is huge and d is small, because the d×d covariance is cheap to build in one pass.

import numpy as np

def pca_fit(X, k):
    """X: (n, d) training rows. Returns mean, components (k, d), explained variance (k,)."""
    X = np.asarray(X, dtype=np.float64)
    mean = X.mean(axis=0)
    Xc = X - mean                                   # centring is not optional
    U, S, Vt = np.linalg.svd(Xc, full_matrices=False)
    # Deterministic signs: make the largest-magnitude loading of each component positive.
    idx = np.argmax(np.abs(Vt), axis=1)
    signs = np.sign(Vt[np.arange(Vt.shape[0]), idx])
    Vt *= signs[:, None]
    var = S**2 / (X.shape[0] - 1)                   # eigenvalues of the covariance
    return mean, Vt[:k], var[:k], var.sum()

def pca_transform(X, mean, comps):
    return (np.asarray(X) - mean) @ comps.T         # (n, k) scores

def pca_inverse(Z, mean, comps):
    return Z @ comps + mean                         # back to (n, d), rank k

X = np.array([[2.5,2.4],[0.5,0.7],[2.2,2.9],[1.9,2.2],[3.1,3.0],
              [2.3,2.7],[2.0,1.6],[1.0,1.1],[1.5,1.6],[1.1,0.9]])
mean, comps, var, total = pca_fit(X, k=1)
print(var / total)          # [0.9632]
print(comps)                # [[0.6779 0.7352]]

Note the sign fix-up. An eigenvector is only defined up to sign, and different libraries, versions or even BLAS builds can return either. On this very data, np.linalg.eigh on the covariance returns PC1 pointing one way and np.linalg.svd returns it pointing the other. Any code that stores components, compares them across refits, or labels an axis "higher means more X" must pin the sign explicitly, as above. scikit-learn applies its own deterministic sign convention, but only within scikit-learn.

Choosing k

There is no universal k. Three rules work in practice, and they answer different questions.

• Variance threshold. Keep the smallest k whose cumulative explained-variance ratio passes a target such as 0.95. Good for compression; in scikit-learn, passing a float such as PCA(n_components=0.95) does this for you.
• Scree elbow. Plot eigenvalues in order and look for where the curve flattens into a noise floor. Good for exploration, subjective for automation.
• Downstream metric. Sweep k and measure what you actually care about: classifier accuracy, retrieval recall@10, reconstruction error on held-out rows. This is the only rule that knows whether the low-variance directions you drop carry the signal. In classification the separating direction can have small variance, and PCA will happily discard it.

For embedding compression the downstream rule dominates. Truncating 1,024-dimensional sentence embeddings to 256 principal components can keep most of the variance and still cost measurable recall on rare queries, so measure recall against the uncompressed index before shipping. Whitening (dividing each score by the square root of its eigenvalue) makes every kept direction unit-variance; it helps some distance-based methods and hurts others because it amplifies the noisiest kept directions, so treat it as a hyperparameter.

Scaling to large data

Full SVD of an n×d matrix costs O(n d min(n, d)) time and needs the whole matrix in memory. That is fine to tens of thousands of rows and a few thousand columns. Beyond that, pick the method by shape:

scikit-learn 1.8 exposes these as svd_solver values full, covariance_eigh, arpack and randomized, with auto choosing among them from the data shape and n_components; IncrementalPCA handles the out-of-core case through partial_fit. On GPUs, the same algebra runs in torch.linalg.svd or torch.pca_lowrank; keep the accumulation in float32 or better, because a covariance built in half precision loses the small eigenvalues entirely.

PCA inside a machine learning pipeline

PCA is a fitted transform, exactly like a scaler, and it leaks the same way. Fitting it on the full dataset before a train/test split lets the test rows shape the axes the model is trained on. Fit on the training split only, store the mean and components, and apply them unchanged to validation, test and live traffic.

Decide on scaling before fitting. If columns share units and their variances are meaningful (pixel intensities, embedding coordinates), centre only. If they mix units (age in years, income in dollars, latency in milliseconds), standardise each to unit variance first, which is the same as running PCA on the correlation matrix. Neither choice is wrong in general; choosing without thinking is.

Treat the fitted artefact as a versioned model. A vector index built on PCA-128 projections is only valid for the exact components that produced it; refitting PCA on new data changes the axes (and possibly their signs and order), so every stored vector must be re-projected together with the new components, never piecemeal.

Failure modes

• Forgot to centre. The first component then points at the mean rather than along the spread. Truncated SVD on raw data does exactly this, which is right for sparse text matrices and wrong almost everywhere else.
• One column dominates. A feature with a large numeric range becomes PC1 by itself. Check each component's loadings; a single loading near ±1 is the symptom.
• Outliers steer the axes. Squared error lets a handful of extreme rows rotate a component towards them. Clip, winsorise, or use a robust variant, and inspect the rows with the largest scores.
• Curved structure. Data on a curve or a manifold needs many linear components to describe it; the explained-variance curve decays slowly and the 2-D plot folds distinct regions on top of each other.
• Reading components as causes. A component is a direction of variance, not a factor of the world. Loadings are not unique when eigenvalues are nearly tied, and small data changes can rotate tied components arbitrarily.
• Missing values. One NaN makes the column mean NaN, which poisons that whole column after centring; NumPy's SVD then fails to converge and scikit-learn rejects the input. Impute inside the same training-only pipeline, or use a method designed for missing entries.
• Drift. Components fitted last quarter may not describe this quarter. Monitor the reconstruction error of live data; a rising error means the subspace no longer fits and a refit, with full re-projection, is due.

Trade-offs against other reducers

Related reading on this site: optimized product quantization (PCA-style rotations before quantizing), the maths inside Faiss, quantizing embeddings, k-means clustering (often run on PCA scores) and power iteration in eigenvector centrality, the same tool that finds a top principal component without a full decomposition.

What to do next

1. Reproduce the ten-point example with the code above and confirm the 96.3% / 3.7% split.
2. On your own dataset, decide centre-only versus standardise, and write the reason down.
3. Wrap scaling and PCA in a single pipeline fitted on the training split only.
4. Plot cumulative explained variance, then sweep k against your real downstream metric.
5. Pin component signs and version the fitted mean and components with the model.
6. Log reconstruction error on live data and alert when it drifts from the training baseline.