Randomized numerical linear algebra
Randomized numerical linear algebra, commonly abbreviated RandNLA, studies algorithms that use random sampling or random projections to compute approximate solutions to problems in numerical linear algebra. Its principal objects include low-rank approximations of matrices, approximate solutions of overdetermined least-squares problems, and estimators for matrix-dependent quantities such as the trace or spectral density. Randomization changes how information is extracted from a matrix, while the accuracy of the resulting computation remains governed by deterministic properties such as singular-value decay and numerical conditioning.
Most RandNLA methods compress a large matrix into a smaller representation that preserves the subspace or quadratic form relevant to the calculation. This compression reduces the number of matrix entries that must be accessed and can replace expensive matrix factorizations with matrix–matrix multiplication. The resulting algorithms are particularly associated with data sets for which the matrix is too large to store in fast memory, is available only through products with vectors, or changes as a stream of observations.
Mathematical framework
Let (A\in\mathbb{R}^{m\times n}) have the singular value decomposition
[ A=U\Sigma V^{\mathsf T}, ]
where the singular values satisfy
[ \sigma_1\geq \sigma_2\geq\cdots\geq 0. ]
The best rank-(k) approximation in the spectral or Frobenius norm is the truncated matrix
[ A_k=U_k\Sigma_kV_k^{\mathsf T}. ]
The Eckart–Young–Mirsky theorem gives
[ \lVert A-A_k\rVert_2=\sigma_{k+1} ]
and
[ \lVert A-A_k\rVert_F^2=\sum_{j>k}\sigma_j^2. ]
Randomized low-rank approximation seeks a matrix (Q\in\mathbb{R}^{m\times \ell}), with (\ell) modestly larger than (k), whose columns approximately span the dominant left singular subspace of (A). The approximation then takes the form
[ A\approx QQ^{\mathsf T}A. ]
A standard randomized range finder uses a random test matrix (\Omega\in\mathbb{R}^{n\times \ell}) and forms
[ Y=A\Omega. ]
An orthonormal basis for the range of (Y) supplies (Q). The compressed matrix
[ B=Q^{\mathsf T}A ]
contains the information needed for a small singular-value decomposition. If
[ B=\widetilde U\Sigma V^{\mathsf T}, ]
then (Q\widetilde U\Sigma V^{\mathsf T}) is an approximate factorization of (A).
The test matrix may have independent Gaussian entries, or it may be drawn from a structured family that accelerates multiplication. A subsampled randomized Hadamard transform combines random sign changes with an orthogonal transform and coordinate sampling. A sparse embedding instead maps each input coordinate into a small number of output coordinates. These constructions differ in their computational costs and in the dimension required to obtain a specified embedding guarantee.
Subspace embeddings
A central concept is the subspace embedding. For a fixed (d)-dimensional subspace represented by the columns of (U\in\mathbb{R}^{m\times d}), a random matrix (S\in\mathbb{R}^{s\times m}) is an (\varepsilon)-embedding when
[ (1-\varepsilon)\lVert Ux\rVert_2^2 \leq \lVert SUx\rVert_2^2 \leq (1+\varepsilon)\lVert Ux\rVert_2^2 ]
holds for every (x\in\mathbb{R}^d). This statement preserves all Euclidean lengths in the subspace, rather than preserving only a finite set of designated vectors.
The embedding property connects RandNLA with the Johnson–Lindenstrauss lemma. A finite collection of points can be projected into a lower-dimensional Euclidean space while approximately preserving pairwise distances. In numerical linear algebra, the relevant set is often an entire subspace, and the embedding dimension therefore depends primarily on its dimension and the required distortion.
For an overdetermined system (Ax\approx b), sketching produces the smaller problem
[ \min_x\lVert S(Ax-b)\rVert_2. ]
When (S) embeds the column space of the augmented matrix ([A\ b]), the sketched residual approximates the original residual uniformly over all candidate vectors. The same geometric principle underlies randomized preconditioners, although those methods use the sketch to construct a better-conditioned transformed system rather than treating the sketched solution as the final approximation.
Error analysis
RandNLA error bounds separate deterministic matrix structure from the probability distribution of the random map. If (V_k) denotes the dominant right singular vectors of (A), then the success of a range finder depends on how well (\Omega) interacts with the subspace spanned by (V_k). For Gaussian test matrices, rotational invariance makes this interaction independent of the orientation of the singular vectors.
A typical expected Frobenius-norm estimate has the form
[ \mathbb{E}\lVert A-QQ^{\mathsf T}A\rVert_F \leq \left(1+\frac{k}{p-1}\right)^{1/2} \lVert A-A_k\rVert_F, ]
where (p=\ell-k) is the oversampling parameter and (p>1). Related high-probability bounds include additional terms controlled by (\sigma_{k+1}) and by the tail energy of the singular spectrum.
When the singular values decay slowly, a powered sample matrix can sharpen the separation between dominant and trailing subspaces:
[ Y=(AA^{\mathsf T})^qA\Omega. ]
Its singular values are proportional to (\sigma_j^{2q+1}), so the ratio between dominant and trailing components is amplified. In finite-precision arithmetic, repeated orthogonalization is associated with this construction because direct formation of high matrix powers can suppress components below the rounding threshold. The arithmetic cost also increases because each power iteration requires further applications of (A) and (A^{\mathsf T}).
Error criteria differ according to the problem being compressed. Relative Frobenius-norm error compares the approximation with the total energy outside the best rank-(k) model. Spectral-norm error instead controls the largest unresolved direction. For least-squares calculations, residual error and solution error are distinct because a small residual perturbation can correspond to a larger change in the solution when (A) is ill-conditioned.
Sampling by importance
Random projections mix coordinates before compression, whereas sampling methods select actual rows or columns. Uniform sampling is effective only when the relevant subspace is sufficiently dispersed across the coordinates. The concentration of that subspace is measured by statistical leverage scores.
For a rank-(k) matrix with right singular-vector matrix (V_k), the column leverage scores are
[ \ell_j=\lVert e_j^{\mathsf T}V_k\rVert_2^2. ]
Their sum is (k), and a large value identifies a coordinate that has substantial influence on the dominant row space. Sampling probabilities proportional to leverage scores prevent such coordinates from being systematically underrepresented. Approximate leverage scores can themselves be computed through randomized embeddings, making the sampling distribution part of the numerical approximation rather than an externally supplied quantity.
The related CUR decomposition expresses a matrix using selected columns and selected rows:
[ A\approx CUR. ]
Unlike factors obtained from a singular-value decomposition, the matrices (C) and (R) consist of entries from the original data matrix. Their interpretation is therefore tied to the original coordinate system, while the middle factor (U) accounts for interactions between the selected coordinates.
Computational model
The significance of a randomized algorithm depends on how matrix access is counted. For a dense matrix stored in memory, multiplication by a Gaussian matrix usually has the same leading-order dependence on (mn) as a deterministic traversal, although it reorganizes much of the work into matrix multiplication. For a sparse matrix, the cost can scale with the number of nonzero entries. For an implicit matrix, the relevant operation is the application of the associated linear map and its adjoint.
Communication costs can dominate arithmetic on hierarchical memory systems and distributed computers. Blocked randomized factorizations aggregate operations into matrix–matrix products, reducing repeated transfers between memory levels. A single-pass method obtains its compressed representation during one traversal of the input, although it generally retains less information than a method that accesses both (A\Omega) and (Q^{\mathsf T}A) in separate passes.
Randomness does not remove finite-precision effects. Orthogonalization can lose accuracy when sampled vectors become nearly dependent, and iterative refinement can be limited by conditioning. The numerical analysis therefore combines probabilistic events with backward-error concepts from conventional matrix computation.
Historical development
Random sampling entered matrix computation through several lines of research, including Monte Carlo estimators and randomized constructions in theoretical computer science. The 1990s development of fast low-distortion embeddings connected these ideas to approximate linear-algebraic computation. Work by Sanjoy Dasgupta and Anupam Gupta clarified geometric forms of the Johnson–Lindenstrauss lemma, while Dimitris Achlioptas established that sparse and discretized random projections could retain comparable distance-preservation properties.
During the 2000s, Petros Drineas and Michael W. Mahoney developed sampling-based algorithms for matrix multiplication, least-squares approximation, and CUR decompositions. Their analysis connected sampling probabilities with leverage and with the structure of singular subspaces. Nariankadu D. Shyamalkumar and David P. Woodruff contributed related developments in sketching and streaming models, where algorithms retain compact summaries rather than the full input.
Randomized range-finding methods were developed into numerically stable matrix algorithms by Vladimir Rokhlin, Arthur Szlam, and Mark Tygert. Nathan Halko, Per-Gunnar Martinsson, and Joel A. Tropp subsequently organized a broad framework for probabilistic low-rank approximation, including deterministic error identities and distribution-specific bounds.
In the early 2010s, You Watanabe and Ming Gu analyzed blocked range approximation under finite-precision and limited-pass computation. Their formulation related the propagation of local projection errors to the global residual and identified the dependence of the resulting approximation on block size, oversampling, and reorthogonalization. This work belonged to the same period in which communication-aware randomized factorizations were being adapted to multicore and distributed-memory architectures.
Later work integrated sketching with optimization, kernel approximation, and iterative linear solvers. These developments placed RandNLA between numerical analysis and the theory of randomized algorithms: the former supplies conditioning and stability models, while the latter supplies concentration inequalities and embedding dimensions.
Statistical interpretation
A matrix sketch can also be interpreted as a randomized summary of a statistical model. In linear regression, leverage scores measure the sensitivity of fitted values to individual observations. A sketch that respects leverage approximately preserves the geometry of the design matrix and therefore the objective function of the regression problem.
Randomized low-rank approximation is related to principal component analysis, since the leading right singular vectors of a centered data matrix define the principal directions. The randomized calculation changes the computational route to those directions but does not change the variance criterion defining them. Statistical uncertainty from finite data and algorithmic uncertainty from random sketching are separate quantities, even when both are expressed probabilistically.
Kernel methods produce matrices whose dimensions scale with the number of observations. Random feature maps and Nyström approximations reduce this cost by replacing the full kernel matrix with a lower-dimensional representation. These methods share the compression viewpoint of RandNLA, although their error analysis also depends on the sampling distribution and the spectral properties of the kernel operator.
Limitations
Randomized compression cannot recover information that is absent from the sampled subspace. Matrices with slowly decaying singular values require a larger sketch, additional power iterations, or a correspondingly larger approximation error. Highly nonuniform leverage can also make uniform row or column sampling unreliable.
A nominal reduction in arithmetic does not necessarily imply a reduction in execution time. Random-number generation, data movement, orthogonalization, and synchronization contribute to the total cost. The balance among these operations depends on the matrix representation and the computing architecture.
Probabilistic guarantees describe a distribution over algorithmic outcomes. Repetition can reduce failure probability, but it does not convert an inappropriate error criterion into an appropriate one. Applications that require deterministic certification can pair a randomized construction with residual estimation or with deterministic verification of the resulting subspace.