Home › Machine learning › Classical ML
algorithm low

PCA and SVD

Problem Someone with high-dimensional data, such as embeddings, wants to shrink it to fewer dimensions, denoise it or plot it, while keeping as much of the variation as possible.

Find orthogonal directions of maximum variance by an SVD of the centred data; randomized SVD gets the top k components in a few passes.

full SVD O(n·d·min(n,d)); randomized top-k O(n·d·(k+p)·(q+1)) with oversampling p and q power iterations
pcasvddimensionality-reductionrandomized-svdlow-ranklinear-algebra
Overview

Any matrix factors as X = U Σ Vᵀ with orthonormal U, V and decreasing singular values. For centred data, V's columns are the principal axes, σᵢ²/(n−1) the variance along each, and the top k give the best rank-k approximation (Eckart-Young). For large matrices, multiply X by a thin random matrix, orthonormalize, and a tiny SVD of the projection finishes the job.

Key facts

Use when

  • Compressing embeddings (768 → 64 before k-means or IVF-PQ), denoising, 2-D plots, whitening; population structure in genetics (Patterson et al.), latent factors for recommender systems (Koren et al.).

Avoid when

  • Non-linear structure: UMAP or t-SNE to plot, autoencoders to compress.
  • Mixed units without standardizing: PCA finds the largest-unit feature.

Deep dive

How it works

In: matrix X (n × d), target rank k · Out: k principal axes, singular values and scores · State: the column means (to centre new data) and V_k · Cost: exact O(n·d·min(n,d)); randomized O(n·d·(k+p)·(q+1))

Exact PCA. Centre X (n×d). Compute the thin SVD (LAPACK gesdd). Components = first k rows of Vᵀ. Scores = U_k Σ_k = X V_k. Explained variance ratio = σᵢ² / Σσ².

Randomized SVD (rank k, l = k + p columns):

  1. Draw Gaussian Ω (d×l). Y = XΩ (n×l).
  2. Power iterations, q times: Y ← X(XᵀY), re-orthonormalizing between multiplies. Each pass multiplies the spectrum by σ², separating the top k from the tail.
  3. Q = orth(Y) (n×l, orthonormal columns).
  4. B = QᵀX (l×d), small.
  5. SVD of B, or eigen-decomposition of the l×l matrix BBᵀ = U_B Σ² U_Bᵀ. Then U = Q·U_B, V = BᵀU_B Σ⁻¹.
  6. Keep the first k.
X (n×d) ──·Ω──► Y (n×l) ──orth──► Q ──Qᵀ·X──► B (l×d) ──small SVD──► U_B, Σ, V
                 ▲      power its: Y ← X·(Xᵀ·Y)       error ≈ σ_{k+1} · (factor shrinking with q)

Code

/// Row-major dense helpers. a: r×c.
fn matmul(a: &[f64], b: &[f64], r: usize, m: usize, c: usize) -> Vec<f64> {   // (r×m)(m×c)
    let mut o = vec![0.0; r * c];
    for i in 0..r { for k in 0..m { let v = a[i * m + k]; for j in 0..c { o[i * c + j] += v * b[k * c + j]; } } }
    o
}
fn transpose(a: &[f64], r: usize, c: usize) -> Vec<f64> {
    let mut t = vec![0.0; r * c];
    for i in 0..r { for j in 0..c { t[j * r + i] = a[i * c + j]; } }
    t
}
/// Orthonormalize the l columns of y (n×l) in place: modified Gram-Schmidt, run twice for stability.
fn orth(y: &mut [f64], n: usize, l: usize) {
    for _ in 0..2 { for j in 0..l {
        for p in 0..j {
            let d: f64 = (0..n).map(|i| y[i * l + j] * y[i * l + p]).sum();
            for i in 0..n { y[i * l + j] -= d * y[i * l + p]; }
        }
        let nrm = (0..n).map(|i| y[i * l + j].powi(2)).sum::<f64>().sqrt();
        for i in 0..n { y[i * l + j] /= nrm; }
    } }
}
/// Cyclic Jacobi eigen-decomposition of a symmetric m×m matrix: returns (eigenvalues, eigenvectors as columns).
fn jacobi_eig(mut a: Vec<f64>, m: usize) -> (Vec<f64>, Vec<f64>) {
    let mut v = vec![0.0; m * m];
    for i in 0..m { v[i * m + i] = 1.0; }
    for _sweep in 0..100 {
        let off: f64 = (0..m).flat_map(|i| (0..m).filter(move |&j| j != i).map(move |j| (i, j))).map(|(i, j)| a[i * m + j].powi(2)).sum();
        if off < 1e-22 { break; }
        for p in 0..m { for q in p + 1..m {
            if a[p * m + q].abs() < 1e-300 { continue; }
            let theta = (a[q * m + q] - a[p * m + p]) / (2.0 * a[p * m + q]);
            let t = theta.signum() / (theta.abs() + (theta * theta + 1.0).sqrt()); // smaller rotation angle
            let (c, s) = (1.0 / (t * t + 1.0).sqrt(), t / (t * t + 1.0).sqrt());
            for k in 0..m {                                         // A ← A·J then Jᵀ·A, V ← V·J
                let (akp, akq) = (a[k * m + p], a[k * m + q]);
                a[k * m + p] = c * akp - s * akq; a[k * m + q] = s * akp + c * akq;
            }
            for k in 0..m {
                let (apk, aqk) = (a[p * m + k], a[q * m + k]);
                a[p * m + k] = c * apk - s * aqk; a[q * m + k] = s * apk + c * aqk;
            }
            for k in 0..m {
                let (vkp, vkq) = (v[k * m + p], v[k * m + q]);
                v[k * m + p] = c * vkp - s * vkq; v[k * m + q] = s * vkp + c * vkq;
            }
        } }
    }
    ((0..m).map(|i| a[i * m + i]).collect(), v)
}

/// Randomized PCA of x (n×d). Returns (top-k components as k×d rows, explained variances).
pub fn randomized_pca(x: &mut [f64], n: usize, d: usize, k: usize, p: usize, q: usize, seed: u64) -> (Vec<f64>, Vec<f64>) {
    for j in 0..d {                                                  // centre columns
        let mu = (0..n).map(|i| x[i * d + j]).sum::<f64>() / n as f64;
        for i in 0..n { x[i * d + j] -= mu; }
    }
    let l = k + p;
    let mut s = seed | 1;
    let mut gauss = || { s ^= s << 13; s ^= s >> 7; s ^= s << 17;
        let u1 = ((s >> 11) as f64 + 1.0) / (1u64 << 53) as f64; s ^= s << 13; s ^= s >> 7; s ^= s << 17;
        (-2.0 * u1.ln()).sqrt() * (6.283185307179586 * (s >> 11) as f64 / (1u64 << 53) as f64).cos() };
    let omega: Vec<f64> = (0..d * l).map(|_| gauss()).collect();
    let xt = transpose(x, n, d);
    let mut y = matmul(x, &omega, n, d, l);
    for _ in 0..q {                                                  // power iterations
        orth(&mut y, n, l);
        let mut z = matmul(&xt, &y, d, n, l);
        orth(&mut z, d, l);
        y = matmul(x, &z, n, d, l);
    }
    orth(&mut y, n, l);
    let b = matmul(&transpose(&y, n, l), x, l, n, d);               // B = Qᵀ X, l×d
    let (evals, ub) = jacobi_eig(matmul(&b, &transpose(&b, l, d), l, d, l), l); // B Bᵀ = U Σ² Uᵀ
    let mut order: Vec<usize> = (0..l).collect();
    order.sort_by(|&i, &j| evals[j].partial_cmp(&evals[i]).unwrap());
    let mut comps = vec![0.0; k * d];
    for (r, &i) in order.iter().take(k).enumerate() {
        let sigma = evals[i].max(0.0).sqrt();
        for c in 0..d { comps[r * d + c] = (0..l).map(|t| b[t * d + c] * ub[t * l + i]).sum::<f64>() / sigma; } // v = Bᵀu/σ
        let lead = comps[r * d..(r + 1) * d].iter().cloned().fold(0.0f64, |m, v| if v.abs() > m.abs() { v } else { m });
        if lead < 0.0 { comps[r * d..(r + 1) * d].iter_mut().for_each(|v| *v = -*v); } // sign convention
    }
    (comps, order.iter().take(k).map(|&i| evals[i] / (n - 1) as f64).collect())
}

fn main() {
    let (n, d) = (2000, 20);
    let mut s = 12345u64;
    let mut u = || { s ^= s << 13; s ^= s >> 7; s ^= s << 17; (s >> 11) as f64 / (1u64 << 53) as f64 - 0.5 };
    let scales: Vec<f64> = (0..d).map(|j| if j < 3 { [10.0, 5.0, 2.0][j] } else { 0.3 }).collect();
    let mut x: Vec<f64> = (0..n * d).map(|i| u() * 12f64.sqrt() * scales[i % d] + 3.0).collect(); // var = scale²
    let (comps, var) = randomized_pca(&mut x, n, d, 3, 10, 2, 99);
    println!("explained variance {:.2?} (truth about [100, 25, 4])", var);
    println!("PC1 leading coords {:.3?}", &comps[..4]);             // ≈ e0
}

Complexity and constants

  • Full thin SVD: O(n·d²) for n ≥ d, with a small constant (Golub and Van Loan list the exact counts). n = 10⁶, d = 1000 is order of 10¹² flops: tens of seconds to minutes on a multicore CPU.
  • Randomized: (2q + 2) passes over X, each an n×d by d×l GEMM: 2·n·d·l·(2q+2) flops. For k = 50, p = 10, q = 2 it is 720·n·d flops; far below the full SVD when l ≪ d.
  • Passes over the data matter more than flops for out-of-core matrices: q + 2 reads (one for Y, q for power iterations, one for B; each power iteration reads X twice if Xᵀ is not stored).
  • Memory: the n×l sketch Y plus the l×d matrix B.

Efficiency

  • Cost: full thin SVD is O(n·d²); n = 10⁶, d = 1000 is order of 10¹² flops, tens of seconds to minutes on a multicore CPU. Randomized SVD with k = 50, p = 10, q = 2 costs 720·n·d flops (above) and q + 2 reads of X. Bound: flops for the full SVD, memory bandwidth and passes for the randomized one.
  • Levers:
    • Randomized range finder — cost scales with l = k + p instead of d (Halko et al. 2011) — k ≪ d.
    • Power iterations q = 1-2 — sharper spectrum for slowly decaying singular values, at 2 extra passes each (above) — noisy data.
    • Covariance route (eigendecompose XᵀX) — O(n·d²) once, then O(d³) — only for small d; it squares the condition number.
    • Incremental or streaming PCA in batches — memory O(batch·d) — data that does not fit in RAM.
    • fp32 instead of fp64 — half the memory and bandwidth — when the spectrum is well separated.
  • Naive — full eigendecomposition of the d×d covariance in fp64 · O(n·d² + d³), unstable for small singular values
  • Standard — LAPACK thin SVD of centered X · O(n·d²)
  • State of the art — randomized SVD with a few power iterations on GPU or out-of-core · O(n·d·l·q)

Implementation notes

  • Centre (and usually scale) features, and store the means: new data must be transformed with the training means.
  • Re-orthonormalize (QR) between every multiply in the power iteration. Without it, all columns converge to the top singular vector and the subspace collapses numerically.
  • With n ≫ d (tall data), the d×d covariance in f64 plus a symmetric eigensolver (syevd) is often fastest. It squares the condition number, so do not trust components with tiny variance.
  • Sign: flip each component so its largest-magnitude coordinate is positive (as scikit-learn's svd_flip does with U). Downstream code that caches projections relies on it.
  • Streaming or larger-than-memory data: incremental PCA (batch SVD updates) or the randomized method with a few out-of-core passes.
  • For embedding compression, keep components until the explained-variance curve flattens, then measure retrieval recall, not just variance.

Pitfalls

  • Skipping centring: the first component points at the mean, not along the spread.
  • PCA on raw features with different units: the result reflects units, not structure.
  • Fitting PCA on train+test: information leaks into evaluation.
  • Using σᵢ²/n vs σᵢ²/(n−1) inconsistently when comparing explained variance with a library.
  • q = 0 on a slowly decaying spectrum (text, embeddings): the randomized estimate is noticeably worse than optimal. Use q = 2 or more.
  • Whitening (dividing by σᵢ) components with tiny σᵢ amplifies noise. Add ε or drop them.

Variants and alternatives

Truncated SVD on uncentred sparse data (LSA)keeps sparsity, no mean subtraction.
Lanczos / ARPACK (scipy.sparse.linalg.svds)iterative top-k for sparse matrices, high accuracy, more passes.
Kernel PCAPCA in a kernel feature space. O(n²) memory.
Probabilistic PCA / factor analysisGaussian latent model, handles missing data via EM.
Autoencoders, UMAP, t-SNEnon-linear reduction.

In production

LAPACK gesdd/gesvd/syevd (via MKL, OpenBLAS, Accelerate), scikit-learn PCA (full, randomized, ARPACK) and TruncatedSVD, torch.linalg.svd and torch.svd_lowrank, Faiss PCAMatrix before IVF-PQ, Spark MLlib RowMatrix.computeSVD; EIGENSOFT/smartpca for genotype data; SVD-based matrix factorization in recommenders.

Implementations

namelangrepoinstallverdictlicencehealthas_ofnotes
linfa-reductionRustrust-ml/linfacargo add linfa linfa-reduction————PCA, diffusion maps, random projection
faerRustsarah-quinones/faer-rscargo add faer————dense SVD and eigendecomposition, pure Rust
nalgebraRustdimforge/nalgebracargo add nalgebra————small and medium SVD, no BLAS needed
ndarray-linalgRustrust-ndarray/ndarray-linalgcargo add ndarray-linalg————LAPACK-backed SVD on ndarray; needs a BLAS/LAPACK backend feature
scikit-learnPythonscikit-learn/scikit-learnpip install scikit-learn————PCA with svd_solver="randomized", IncrementalPCA

Sources

  • Patterson, Price, Reich, Population Structure and Eigenanalysis, PLoS Genetics 2006. paper
  • Koren, Bell, Volinsky, Matrix Factorization Techniques for Recommender Systems, IEEE Computer 2009. paper
  • Halko, Martinsson, Tropp, Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions, SIAM Review 2011. paper · arXiv
  • Eckart and Young, The Approximation of One Matrix by Another of Lower Rank, Psychometrika 1936. paper
  • Golub and Van Loan, Matrix Computations, 4th ed., 2013.
  • Jolliffe, Principal Component Analysis, 2nd ed., Springer 2002.
  • Tipping and Bishop, Probabilistic Principal Component Analysis, JRSS B 1999. paper

Used by

attitude determination · plannedEstimate orientation from vector observations (TRIAD, QUEST/Wahba) and gyros
camera calibration · plannedIntrinsic/extrinsic calibration (Zhang's method, distortion models, bundle adjustment)
cma es · plannedCovariance matrix adaptation evolution strategy for derivative-free continuous optimization
fiducial marker localization · plannedLocalize from known floor/wall markers (AprilTag, 2D codes) fused with odometry
inverse kinematics · plannedCCD/FABRIK/Jacobian solvers posing joint chains to reach targets
quantum gate synthesis · plannedDecompose unitaries into native gates (KAK, Solovay-Kitaev, Clifford+T synthesis)
tensor network contraction · plannedContract tensor networks with good contraction orders; MPS and slicing for circuit simulation

Systems (12):

model compressionmodel monitoringeigensolverarray computingstatistics environmentrisk engineportfolio managementsingle cell analysishts screeningpredictive maintenancefab mesinterpretability platform

Alternatives and related

k-means clusteringAn engineer needs to split a large set of vectors or points into groups of similar items, for example to build a compact codebook, a colour palette or coarse partitions for search.IVF-PQVector search over hundreds of millions to billions of embeddings, where full vectors do not fit in memory, needs compact storage and fast approximate nearest-neighbour lookup, on CPU or GPU.EmbeddingsA neural network has to take categorical inputs, such as words, users, items or board squares, and represent each as a useful numeric vector, without paying for a huge one-hot input.Linear and logistic regressionA data team wants a simple, fast and interpretable model that predicts a number or a probability from features, and a baseline that any fancier model must beat to justify its complexity.LU and Cholesky factorizationSolving dense linear systems Ax = b, often repeatedly with the same matrix and different right-hand sides, or testing and sampling with covariance matrices, without inverting the matrix.Conjugate gradientLarge sparse linear systems with a symmetric positive definite matrix (finite elements, pressure solves, graph Laplacians) are too big for direct factorization. This solves them using only matrix-vector products.Weight quantizationRunning a large language model for one user at a time is slow and needs a lot of GPU memory, because every weight must be read for each generated token. Engineers want faster, smaller deployments.