Singular Value Decomposition › Randomized SVD for Scale

Even with sparse algorithms, exact SVD becomes slow for matrices with millions of rows or columns. Modern production NLP systems may want to compute embeddings for vocabularies of 500,000 to 1 million words. Randomized SVD, introduced by Halko, Martinsson, and Tropp in their influential 2011 paper "Finding structure with randomness," provides a principled approximation that dramatically reduces computation while maintaining high accuracy.

The algorithm works in two stages that together are much cheaper than direct SVD.

Stage 1: Range finding. The goal of this stage is to find a small orthonormal matrix $\mathbf{Q}$ whose columns approximately span the range (column space) of $\mathbf{M}$, restricted to the $k$ dominant directions. Generate a random matrix $\mathbf{G}$ of shape $n \times (k + p)$, where $p$ is a small oversampling parameter, typically 5 to 10. The oversampling provides insurance that the random projection does not miss any important directions. Form the sample matrix:

Compute an orthonormal basis $\mathbf{Q}$ for the column space of $\mathbf{Y}$ using QR decomposition. The key insight is that multiplying $\mathbf{M}$ by a random matrix produces a "sketch" of the column space of $\mathbf{M}$. If $\mathbf{M}$ has $k$ dominant singular directions with much larger singular values than the rest, then $\mathbf{Y}$ will be approximately a linear combination of those $k$ directions plus small random perturbations. QR decomposition extracts an orthonormal basis for this sketch.

Stage 2: SVD of the projected matrix. Form the smaller matrix:

$\mathbf{B}$ has dimensions $(k + p) \times n$. Rather than being enormous, $\mathbf{B}$ is small in its first dimension: only $k + p$ rows. Compute the exact SVD of $\mathbf{B}$:

Recover the approximate left singular vectors as $\mathbf{U} = \mathbf{Q} \hat{\mathbf{U}}$.

The computational saving is substantial. The expensive operations in the randomized algorithm are the matrix-vector products in forming $\mathbf{Y}$ and $\mathbf{B}$, each costing $O(\text{nnz})$ for a sparse matrix. The full SVD of $\mathbf{B}$ costs $O((k+p)^2 n)$, which is small when $k \ll m$. Compared to exact SVD at $O(k \cdot \text{nnz})$, the asymptotic cost is similar, but the constant factors are much smaller because the algorithm avoids the iterative Lanczos refinement. For very large matrices, the difference is an order of magnitude in wall time.

The theoretical guarantee is that the error of randomized SVD is, with high probability, only slightly larger than the error of the best rank-$k$ approximation. Specifically:

where the expectation is over the random choice of $\mathbf{G}$ and $\sigma_{k+1}$ is the $(k+1)$-th singular value. When $p$ is even modestly large (say, $p = 10$) and the singular values decay rapidly (as they do in NLP matrices), the bound is tight and the approximation is excellent.

An optional but often used power iteration step can further improve accuracy: instead of computing $\mathbf{Y} = \mathbf{M} \mathbf{G}$, compute $\mathbf{Y} = (\mathbf{M} \mathbf{M}^T)^q \mathbf{M} \mathbf{G}$ for small $q$ (typically 1 or 2). Each power iteration squares the singular values, making the dominant ones even more prominent relative to the small ones and producing a better sketch of the dominant subspace. The cost is $2q$ additional matrix-vector products.

Scikit-learn's TruncatedSVD uses randomized SVD internally for large matrices, making it the practical tool of choice for LSA at scale. The algorithm='randomized' option (the default) invokes this randomized approach.