Dimensionality Reduction and Principal Component Analysis

Why reduce dimensionality?

High-dimensional data is expensive to store and process, hard to visualise, prone to overfitting, and suffers from the curse of dimensionality (see k-NN). Often the data actually lies near a lower-dimensional subspace because features are correlated.

PCA: the idea

Find new orthogonal axes — principal components — ordered so that the first captures as much variance as possible, the second as much of the remaining variance as possible, and so on. Keep the first $k$ and discard the rest.

Two equivalent views

1. Maximum variance: the first direction $\mathbf v_1$ maximises the variance of the projections, $\max_{\lVert\mathbf v\rVert = 1}\mathbf v^\top\Sigma\mathbf v$. By the Rayleigh quotient (see quadratic forms) the answer is the top eigenvector, with variance $\lambda_1$.
2. Minimum reconstruction error: the $k$-dimensional subspace minimising $\sum_i\lVert\mathbf x_i - \hat{\mathbf x}_i\rVert^2$ (squared perpendicular distances) is spanned by the top $k$ eigenvectors; the error equals the sum of the discarded eigenvalues.

The algorithm

1. Centre the data: $\mathbf x_i \leftarrow \mathbf x_i - \bar{\mathbf x}$ (and standardise if features have different units).
2. Compute the covariance matrix $\Sigma = \frac{1}{n-1}X^\top X$ ($p \times p$).
3. Eigendecompose $\Sigma = V\Lambda V^\top$ with $\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_p \ge 0$.
4. Keep the first $k$ eigenvectors $V_k = [\mathbf v_1 \cdots \mathbf v_k]$.
5. Scores (new coordinates): $Z = XV_k$ ($n\times k$). Reconstruction: $\hat X = ZV_k^\top$ (+ mean).

Equivalently, take the SVD $X = U\Sigma V^\top$ of the centred data: principal directions are the columns of $V$, variances $\sigma_j^2/(n - 1)$, scores $U\Sigma$. The SVD route is numerically preferred.

import numpy as np

def pca(X, k):
    Xc = X - X.mean(axis=0)
    U, s, Vt = np.linalg.svd(Xc, full_matrices=False)
    var = s ** 2 / (len(X) - 1)
    return Xc @ Vt[:k].T, Vt[:k], var / var.sum()     # scores, directions, explained ratios

Worked example

Take the points $(2, 0), (0, 2), (3, 3), (-1, -1)$. Their mean is $(1, 1)$, so the centred points are $(1, -1), (-1, 1), (2, 2), (-2, -2)$.

• Covariance ($n - 1 = 3$): $\Sigma = \frac13\begin{pmatrix}1+1+4+4 & -1-1+4+4\\ \cdot & 1+1+4+4\end{pmatrix} = \frac13\begin{pmatrix}10 & 6\\6 & 10\end{pmatrix}$.
• Eigenvalues: $\frac{16}{3}$ along $\frac{1}{\sqrt2}(1, 1)$ and $\frac43$ along $\frac{1}{\sqrt2}(1, -1)$.
• PC1 explains $\frac{16/3}{20/3} = 80\%$ of the variance.
• Scores on PC1: $(1, -1) \to 0$, $(-1, 1) \to 0$, $(2, 2) \to 2\sqrt2$, $(-2, -2) \to -2\sqrt2$.

How many components?

• Explained variance ratio $\frac{\lambda_j}{\sum_i\lambda_i}$; keep enough components for a cumulative threshold (e.g. 90–95%).
• Scree plot: plot eigenvalues against index and look for the elbow.
• Kaiser rule (standardised data): keep components with $\lambda > 1$.
• Or choose $k$ by downstream cross-validated performance.

Practical points

Limitations

• Linear only — curved manifolds need kernel PCA, Isomap, t-SNE, UMAP or autoencoders.
• Maximises variance, which need not align with class separation (use LDA when labels exist).
• Sensitive to outliers and to feature scaling; components can be hard to interpret.

Quick revision