A randomized algorithm for principal component analysis

Principal component analysis (PCA) requires the computation of a low-rank approximation to a matrix containing the data being analyzed. In many applications of PCA, the best possible accuracy of any rank-deficient approximation is at most a few digits (measured in the spectral norm, relative to the spectral norm of the matrix being approximated). In such circumstances, efficient algorithms have not come with guarantees of good accuracy, unless one or both dimensions of the matrix being approximated are small. We describe an efficient algorithm for the low-rank approximation of matrices that produces accuracy very close to the best possible, for matrices of arbitrary sizes. We illustrate our theoretical results via several numerical examples.

💡 Research Summary

The paper addresses a fundamental challenge in large‑scale principal component analysis (PCA): obtaining a low‑rank approximation of a data matrix with provable accuracy while keeping computational costs modest. Classical deterministic algorithms (e.g., full singular value decomposition or Lanczos‑based methods) guarantee optimal spectral‑norm error but become prohibitive when either dimension of the matrix is large. Conversely, many fast heuristics lack rigorous error bounds, especially when the desired rank k is much smaller than the matrix dimensions and the best possible approximation is only accurate to a few digits.

To bridge this gap, the authors propose a randomized algorithm that combines oversampling, subspace iteration, and a final deterministic SVD on a small projected matrix. The method proceeds as follows. Given an m × n matrix A and a target rank k, one draws a random Gaussian matrix Ω of size n × ℓ, where ℓ = k + p and p (typically 5–10) is an oversampling parameter. The product Y = AΩ captures the column space of A with high probability. A QR factorization of Y yields an orthonormal basis Q ∈ ℝ^{m×ℓ}. The algorithm then forms the reduced matrix B = QᵀA (ℓ × n) and computes its exact SVD, B = ÛΣ̂V̂ᵀ. The final rank‑k approximation is Ã_k = Q Û_k Σ̂_k V̂_kᵀ.

The key theoretical contribution is a set of spectral‑norm error bounds that are essentially optimal. Without subspace iteration (q = 0), the authors prove
‖A − Ã_k‖₂ ≤ (1 + 4√{k + p}·√{min(m,n)}·σ_{k+1}/σ_k)·σ_{k+1},
where σ_{k+1} is the (k+1)‑st singular value of A. By applying q ≥ 1 steps of subspace iteration—replacing Y by Y ← (AAᵀ)^q AΩ—the bound improves exponentially:
‖A − Ã_k‖₂ ≤