PCA

Principal Component Analysis (PCA) finds a set of orthogonal directions (principal components) that capture the most variance in the data, giving a linear low-dimensional representation. It is the workhorse linear dimensionality reduction method.

Variance Maximization

Assume the data $X\in {\mathbb{R}}^{n\times d}$ is mean-centered. The first principal component is the unit vector $w$ that maximizes the variance of the projected data:

$$w_1=\arg \underset{\Vert w\Vert =1}{\max }{w}^{\top }\Sigma w,\Sigma =\frac{1}{n}{X}^{\top }X,$$

where $\Sigma$ is the $d\times d$ covariance matrix. Subsequent components maximize the same quantity subject to being orthogonal to the earlier ones.

Covariance Eigen-Decomposition

The Lagrangian for the constrained problem gives $\Sigma w=\lambda w$: the principal components are the eigenvectors of $\Sigma$, and the variance captured by each is its eigenvalue $\lambda$. Ranking eigenvectors by decreasing $\lambda$ gives the components in order of importance.

Relation to Singular Value Decomposition

In practice we do not form $\Sigma$; we take the SVD of the centered data directly:

$$X=US{V}^{\top }.$$

Then the principal components (directions) are the columns of $V$, and $\Sigma =\frac{1}{n}V{S}^2{V}^{\top }$, so the variance along component $j$ is ${\sigma }_j^2/n$ where ${\sigma }_j$ is the $j$-th singular value. SVD is more numerically stable than eigen-decomposing $X^{\top }X$, which squares the condition number.

Code (numpy and torch via SVD)

import numpy as np
import torch
 
def pca(X, k):
    Xc = X - X.mean(axis=0, keepdims=True)        # center
    U, S, Vt = np.linalg.svd(Xc, full_matrices=False)
    components = Vt[:k]                            # (k, d)
    scores = Xc @ components.T                     # (n, k) projection
    explained = (S[:k] ** 2) / (S ** 2).sum()
    return scores, components, explained
 
# torch version, projection written with einsum
def pca_torch(X, k):
    Xc = X - X.mean(dim=0, keepdim=True)
    U, S, Vh = torch.linalg.svd(Xc, full_matrices=False)
    V = Vh[:k]                                     # (k, d)
    scores = torch.einsum('nd,kd->nk', Xc, V)      # project onto components
    return scores, V, (S[:k] ** 2) / (S ** 2).sum()

Whitening

Whitening rescales each component to unit variance, decorrelating the features:

$$X_{\text{white}}=X_{c}V{S}^{-1}\sqrt{n}.$$

The result has identity covariance. This helps downstream methods that assume isotropic inputs, but it amplifies low-variance (noisy) directions, so it is often combined with dropping small components.

Choosing the Number of Components

Pick $k$ by the cumulative explained variance ratio:

$$r(k)=\frac{{\sum }_{j=1}^k{\sigma }_j^2}{{\sum }_{j=1}^d{\sigma }_j^2},$$

and choose the smallest $k$ with $r(k)$ above a threshold (for example 0.95). A scree plot (eigenvalue vs index) often shows an “elbow” marking a natural cutoff.

Pitfalls

• Always center the data first; forgetting to subtract the mean mixes the mean into the first component.
• PCA is scale-sensitive. Standardize features when they have different units, otherwise large-magnitude features dominate.
• PCA is linear: it cannot unfold curved manifolds. For nonlinear structure use kernel PCA or t-SNE / UMAP for visualization.

Related

• Singular Value Decomposition is the computational backbone.
• t-SNE is often applied to PCA-reduced data for visualization.