Randomized SVD

The primary advantage of Randomized SVD is its sub-quadratic computational complexity for approximating the top-k singular values/vectors. For a large m×n matrix, classical deterministic SVD requires O(m n min(m, n)) operations. Randomized SVD reduces this to approximately O(m n log(k) + (m+n) k²) through the use of random projections. This makes it feasible for massive datasets where full SVD is computationally prohibitive. The algorithm is particularly effective when the numerical rank k is much smaller than the matrix dimensions.

The core technique involves constructing a low-dimensional subspace that captures the action of the input matrix. The algorithm proceeds in three key stages:

• Random Sampling: A random test matrix Ω is used to form a sketch Y = AΩ, projecting the high-dimensional data into a lower-dimensional space.
• Subspace Iteration (Power Method): To improve accuracy for matrices with slowly decaying singular values, the sketch is refined via repeated multiplication with A and its transpose: Y = (AAᵀ)^q AΩ. This power iteration amplifies the dominant singular directions.
• Orthogonalization: The columns of Y are orthonormalized (e.g., via QR decomposition) to form basis Q for the approximate range of A.

Randomized SVD is explicitly designed for the fixed-rank approximation problem, where the target rank k is specified in advance. The algorithm directly computes an approximation A ≈ UΣVᵀ, where U and V have k orthonormal columns and Σ is a k×k diagonal matrix. This is in contrast to computing a full SVD and then truncating. The approximation quality is controlled by two parameters: the target rank k and the number of oversampling parameters p (e.g., using k+p random vectors to provide a safety margin). The oversampling ensures the probabilistic capture of the dominant subspace.

The algorithm requires only the ability to compute the matrix-vector products Ax and Aᵀy. It never needs to explicitly form or store the full dense matrix A in fast memory. This matrix-free property enables:

• Streaming data applications where the matrix is generated on-the-fly.
• Out-of-core computation for matrices larger than available RAM, as data can be loaded in blocks from disk.
• Application to implicitly defined linear operators (e.g., fast Fourier transforms, integral operators). This makes it a cornerstone of modern numerical linear algebra for big data.

Randomized SVD provides rigorous probabilistic guarantees on approximation error. For a given target rank k and oversampling parameter p, the algorithm produces an approximation satisfying ‖A - UΣVᵀ‖ ≤ (1 + ε) σ_{k+1} with high probability, where σ_{k+1} is the (k+1)-th true singular value. The error decays exponentially with the number of power iterations q. These bounds are derived from the Johnson-Lindenstrauss lemma and matrix concentration inequalities, ensuring reliability for scientific computing. The failure probability is user-controllable and typically extremely small (e.g., 10^-15).

In practice, Randomized SVD is implemented in major numerical libraries (e.g., sklearn.utils.extmath.randomized_svd, SciPy). Key practical considerations include:

• Choosing oversampling: Typically p=10 or p=0.1*k provides a good balance.
• Power iterations: 1-4 iterations (q=1 to 4) suffice for most matrices with reasonable spectral gaps.
• Blocked operations: Using a block of random vectors (instead of a single vector) improves performance via optimized BLAS-3 matrix-matrix operations.
• Adaptive rank determination: Variants can estimate the numerical rank on-the-fly by monitoring the decay of singular values in the small projected problem. It is the default method for PCA on large datasets in many machine learning pipelines.

Singular Value Decomposition (SVD) is the fundamental matrix factorization that Randomized SVD approximates. For any matrix (A \in \mathbb{R}^{m \times n}), the full SVD is (A = U \Sigma V^T), where (U) and (V) are orthogonal matrices containing the left and right singular vectors, and (\Sigma) is a diagonal matrix of non-negative singular values. This decomposition is exact and reveals the matrix's rank and fundamental subspaces, but computing it for large matrices is computationally expensive, motivating randomized approximations.

Truncated SVD is the target output of the Randomized SVD algorithm. Instead of computing the full decomposition, truncated SVD retains only the top (k) largest singular values and their corresponding singular vectors, producing a rank-(k) approximation (A_k = U_k \Sigma_k V_k^T). This is the optimal rank-(k) approximation in both the spectral and Frobenius norms, as guaranteed by the Eckart–Young theorem. Randomized SVD provides a fast, approximate method to compute this truncated decomposition for matrices where traditional methods are prohibitive.

Random projection is the core probabilistic technique that enables Randomized SVD's efficiency. It uses a random matrix (\Omega) to project the high-dimensional input matrix (A) onto a lower-dimensional subspace, capturing its essential range with high probability (per the Johnson-Lindenstrauss lemma).

• In Randomized SVD, this creates a sketch matrix (Y = A\Omega).
• The columns of (Y) form a basis for an approximate subspace spanning the dominant singular vectors of (A).
• This step reduces the problem size, allowing standard dense SVD to be performed on a much smaller matrix.

Power iteration (generalized to subspace iteration) is an optional but crucial refinement step in Randomized SVD to improve accuracy, especially for matrices with slowly decaying singular values. The basic randomized sketch (Y = A\Omega) may poorly capture the dominant subspace if singular values are close. Applying (q) steps of power iteration, by computing (Y = (AA^T)^q A\Omega), rapidly amplifies the dominant singular directions, yielding a much higher quality basis for the approximation. This comes at the cost of (2q) additional matrix-matrix multiplications.

The Nyström method is a closely related randomized algorithm for approximating large symmetric positive semidefinite (PSD) matrices, such as kernel matrices in machine learning. While Randomized SVD approximates the matrix from its column space, the Nyström method samples a subset of columns (C) and uses them to construct a low-rank approximation of the form (\tilde{A} = C W^+ C^T), where (W) is the intersection of the sampled columns and rows. For PSD matrices, it can be viewed as a variant of Randomized SVD that enforces a specific structural constraint on the approximation.

Krylov subspace methods, like the Lanczos algorithm, are classical deterministic iterative alternatives for computing a few extreme singular values/vectors of large, sparse matrices. They build an orthogonal basis for the Krylov subspace (\mathcal{K}_k(A, v)). Compared to Randomized SVD:

• Lanczos: Excellent for sparse matrices, highly accurate for extreme singular values, but suffers from numerical instability (loss of orthogonality) and requires matrix-vector products which can be irregular.
• Randomized SVD: More robust, inherently parallelizable (uses matrix-matrix products), and often faster for dense matrices or when a moderate-precision approximation suffices. It is less sensitive to the matrix spectrum's gap.