Randomized Numerical Linear Algebra (RNLA)

• Randomized Numerical Linear Algebra (RNLA) leverages randomness for approximating, factorizing, solving, and analyzing large matrices, targeting computational regimes constrained by memory, data access, communication, and work.
• Applications range from least-squares regression, low-rank approximation, and kernel methods to scientific computing, graph algorithms, and quantum information.
• Key benefits include reduced computational complexity and memory requirements, making RNLA suitable for large-scale and distributed systems, though it may not replace deterministic methods for exact solutions.

Randomized Numerical Linear Algebra (RNLA), also called Randomized Linear Algebra (RLA) or RandNLA, is a family of numerical methods that uses randomness as an algorithmic resource for approximating, factoring, solving, and analyzing large matrices and linear operators. Its central mechanism is to replace an expensive matrix computation by a smaller randomized representation—a sketch, random projection, sampled subset, stochastic estimator, or randomized iterate—that preserves the geometry or observables relevant to the downstream task. The resulting reduced problem is then solved by deterministic numerical linear algebra, iterative refinement, or further randomized computation. RNLA addresses computational regimes constrained by memory capacity, data-access latency, communication, synchronization, and floating-point work, with applications ranging from least-squares regression and low-rank approximation to kernel methods, scientific computing, graph algorithms, inverse problems, quantum information, and semidefinite optimization (Mahoney, 2016, 2002.01387, Murray et al., 2023).

1. Scope, computational model, and historical development

RNLA differs from conventional randomized numerical analysis in that randomness is deliberately injected into the algorithm rather than assumed to originate in the input data. Once the random seed or random operator is fixed, the computation is deterministic. The objective is typically to reduce the dimensions of a problem before applying QR, SVD, Cholesky, eigendecomposition, or an iterative solver.

For a matrix $A$, a sketch may be formed by left multiplication,

$SA,$

which reduces the number of rows, or by right multiplication,

$AS,$

which reduces the number of columns. The sketching operator may be a sampling-and-rescaling matrix, a Gaussian or sub-Gaussian projection, a sparse embedding, a randomized Hadamard or Fourier transform, a hashing operator, or a data-aware leverage-score sampler. The relevant approximation target may be a norm, inner product, subspace, least-squares objective, dominant singular subspace, matrix product, trace, or matrix function.

The contemporary field developed from randomized algorithms, streaming computation, dimensionality reduction, and approximate matrix multiplication. The Johnson–Lindenstrauss principle established that a finite set of high-dimensional points can be embedded into a lower-dimensional space while approximately preserving pairwise distances. RNLA extended this principle from finite point sets to matrix ranges and subspaces through subspace embeddings, which preserve the geometry of every vector in a specified low-dimensional subspace (Murray et al., 2023).

RNLA subsequently integrated with classical numerical linear algebra. Randomization is commonly used to construct:

• a smaller least-squares problem solved by QR or SVD;
• a right preconditioner for LSQR, CG, CGLS, LSMR, or Chebyshev iteration;
• a randomized range approximation followed by a deterministic small SVD;
• a pivot structure for QR, LU, or related factorizations;
• a low-rank approximation or Nyström factorization;
• a trace or matrix-function estimator;
• a sparse or distributed representation suitable for streaming and parallel computation.

This integration is significant because the randomized component need not replace numerical analysis. Instead, it can produce a coordinate system, basis, or preconditioner in which classical stable methods operate more efficiently. Recent work has extended this perspective to backward-stable randomized least-squares solvers: FOSSILS is proved backward stable under explicit assumptions, while SPIR is empirically backward stable and is supported by a proposed finite-precision framework (Epperly et al., 2024).

The access model is central. In the pass-efficient model, an algorithm accesses the input through a small constant number of sequential passes, uses space sublinear in the input size, and performs additional computation sublinear in the input length (Mahoney, 2016). Other RNLA algorithms assume matrix–vector products, sampled entries, sampled rows or columns, dynamic length-squared sampling, or an oracle capable of producing nonzero entries of a matrix acting on a sparse representation.

2. Sketches, embeddings, and probabilistic guarantees

A sketching matrix $S$ is useful when it approximately preserves the geometry relevant to a computation. If $U\in\mathbb{R}^{n\times d}$ has orthonormal columns, then $S$ is an $\epsilon$-subspace embedding for $\operatorname{range}(U)$ when

$\left\|I-(SU)^T(SU)\right\|_2\leq \epsilon.$

Equivalently,

$$(1-\epsilon)\|x\|_2^2 \leq \|SUx\|_2^2 \leq (1+\epsilon)\|x\|_2^2$$

for every $SA,$0. For a matrix $SA,$1, this implies approximate preservation of the Gram matrix and singular geometry. An oblivious subspace embedding is drawn without knowledge of the particular subspace and succeeds with probability at least $SA,$2 for every subspace of the prescribed dimension (Murray et al., 2023).

Sketching distributions

Dense sketches include Gaussian, Rademacher, uniform sub-Gaussian, and Haar-distributed orthogonal or unitary matrices. They have strong theoretical guarantees and are compatible with matrix–matrix multiplication, but explicit dense application can cost $SA,$3 for a matrix with $SA,$4 rows and $SA,$5 columns.

Sparse sketches reduce storage, arithmetic, and data movement. They include CountSketch, sparse Johnson–Lindenstrauss transforms, Clarkson–Woodruff transforms, OSNAPs, row sampling, and column sampling. A Clarkson–Woodruff sketch maps each input row to a randomly selected output row and multiplies it by an independent random sign. Its application can be performed through signed accumulation rather than dense matrix multiplication, with cost proportional to $SA,$6 for a dense input and $SA,$7 for a sparse input in the implementation described by “Sketch ’n Solve” (Lavaee, 2024).

Subsampled randomized trigonometric transforms combine a random sign or phase modulation, a fast Fourier, cosine, or Walsh–Hadamard transform, and coordinate restriction. Randomized Hadamard transforms spread the mass of vectors and approximately flatten leverage scores. Their fast application is compatible with streaming, distributed, and memory-sensitive computation (Mahoney, 2016, Drineas et al., 2017).

The main accuracy tools are concentration inequalities and perturbation theory. Markov and Chebyshev inequalities convert expectation bounds into probability bounds but often yield weak dependence on the failure probability. Chernoff, Hoeffding, Bernstein, Azuma–Hoeffding, and matrix concentration inequalities provide sharper high-probability guarantees. Matrix concentration is needed when the target is a spectral norm or subspace embedding rather than an entrywise or Frobenius-norm quantity (Mahoney, 2016).

The principal RNLA trade-off is

$SA,$8

The relevant parameters include sketch dimension, oversampling, coherence, leverage-score domination, spectral decay, condition number, effective dimension, stable rank, precision, and the norm used to measure error. The same sketch can be adequate for a Frobenius-norm approximation but inadequate for a spectral-norm guarantee.

3. Sampling, leverage scores, and structured random subsets

Sampling-based RNLA represents a large matrix or matrix product using a small number of rescaled rows, columns, entries, or rank-one terms. If

$SA,$9

then one may sample $AS,$0 terms using probabilities $AS,$1 and form an unbiased estimator. For column–row sampling, the optimal probabilities for minimizing the expected Frobenius error are proportional to

$AS,$2

With nearly optimal probabilities, the expected error decreases as $AS,$3, while high-probability bounds introduce logarithmic dependence on the failure probability (Mahoney, 2016, Drineas et al., 2017).

Leverage scores

For an orthonormal basis $AS,$4 of a $AS,$5-dimensional subspace, the leverage score of row $AS,$6 is

$AS,$7

The scores satisfy

$AS,$8

Leverage scores measure the contribution of rows to the geometry of a subspace. Sampling proportionally to leverage scores provides near-optimal subspace embeddings and removes the dependence on coherence that affects uniform sampling. Uniform sampling can fail when the subspace is concentrated on a small number of coordinates. The randomized Hadamard transform offers a data-oblivious alternative by flattening leverage scores before uniform sampling (Drineas et al., 2017, Murray et al., 2023).

Ridge leverage scores regularize this construction. For a positive regularization parameter $AS,$9, the ridge leverage score of row $S$0 is

$S$1

Their sum is the effective dimension,

$S$2

Regularization therefore replaces exact rank by a statistical or effective dimension. Ridge leverage sampling is useful for ridge regression, kernel approximation, low-rank approximation, and preconditioning (2002.01387, Chepurko et al., 2020).

Determinantal sampling

Determinantal point processes (DPPs) replace independent sampling by negatively dependent sampling. A DPP assigns subset probabilities proportional to determinants:

$S$3

Because determinants measure squared volume, DPPs favor subsets that are informative and mutually nonredundant. Projection DPPs associated with a rank-$S$4 subspace select exactly $S$5 rows, whereas independent leverage-score sampling uses the same marginal probabilities but typically requires $S$6 samples. The DPP’s negative dependence reduces redundancy and rank deficiency (Dereziński et al., 2020).

Volume sampling is the projection-DPP specialization for least squares. If $S$7 has full column rank and $S$8 is a size-$S$9 volume sample, then

$U\in\mathbb{R}^{n\times d}$0

Under the stated general-position assumptions, the expected loss satisfies

$U\in\mathbb{R}^{n\times d}$1

For ridge regression, an $U\in\mathbb{R}^{n\times d}$2-ensemble DPP produces an unbiased minimum-norm estimator whose expectation equals the ridge solution. DPPs also yield expected-error guarantees for low-rank approximation and Nyström approximation with sample sizes that are worst-case optimal for the stated results (Dereziński et al., 2020).

Sampling design can also be optimized for statistical inference rather than only geometric preservation. Asymptotic analyses of randomized least-squares estimators distinguish inference for the population parameter $U\in\mathbb{R}^{n\times d}$3 from approximation of the full-data OLS estimator. The resulting optimal distributions include inverse-covariance, root-leverage, predictor-length, inverse-covariance negative-leverage, root-leverage negative-leverage, and predictor-length negative-leverage sampling (Ma et al., 2020).

4. Core computational tasks

Matrix multiplication and Gram approximation

Randomized matrix multiplication estimates a product by sampling rank-one outer products or applying a random projection. The estimator is unbiased when sampling-and-rescaling matrices satisfy an expectation identity such as

$U\in\mathbb{R}^{n\times d}$4

Approximate matrix multiplication is useful when approximate products suffice, matrices arrive in streams, or communication is expensive. Its basic sampling error typically scales as $U\in\mathbb{R}^{n\times d}$5, making it more suitable for coarse approximation than for high-precision spectral computation (2002.01387).

Graph sparsification is a related application. A graph Laplacian can be written as a sum of weighted rank-one edge Laplacians. Sampling edges according to effective resistance produces an unbiased sparse Laplacian, and $U\in\mathbb{R}^{n\times d}$6 sampled edges yield a spectral approximation with high probability in the stated construction (2002.01387).

Least-squares regression

For

$U\in\mathbb{R}^{n\times d}$7

classical QR costs $U\in\mathbb{R}^{n\times d}$8. RNLA offers three principal paradigms.

Sketch-and-solve forms $U\in\mathbb{R}^{n\times d}$9 and $S$0, then solves

$S$1

If $S$2 is a suitable subspace embedding, the sketched residual is close to the optimal residual. This approach reduces the number of constraints but does not generally provide high-accuracy solution-vector approximation for ill-conditioned $S$3.

Sketch-and-precondition factors $S$4, typically by QR or SVD, and uses the resulting factor to precondition the original least-squares problem. The original matrix is retained, while the sketch supplies a nearly orthogonal coordinate system. The method is particularly effective for highly overdetermined problems and is implemented in methods such as Blendenpik, LSRN, and related solvers (Murray et al., 2023).

Iterative sketching repeatedly corrects a current solution using sketched subproblems, connecting RNLA to stochastic optimization, randomized Kaczmarz, sketch-and-project, and stochastic gradient methods.

Recent stability work distinguishes forward accuracy, residual accuracy, strong forward stability, and backward stability. FOSSILS combines a randomized subspace embedding, preconditioned iterations on implicitly applied normal equations, a sketch-and-solve initial guess, and iterative refinement. Under explicit assumptions it is backward stable with complexity $S$5 up to logarithmic factors. SPIR uses LSQR or a related Krylov method with iterative refinement and is empirically backward stable, although the complete LSQR-based finite-precision theorem remains open (Epperly et al., 2024).

Statistical RNLA analyses additionally derive asymptotic normality for sampling estimators. Under stated regularity conditions, randomized least-squares estimators are asymptotically unbiased, and their asymptotic covariance includes both the ordinary OLS variance and a sampling-induced variance term. This permits approximate confidence intervals and hypothesis tests, not only pointwise computational error bounds (Ma et al., 2020).

Low-rank approximation and randomized SVD

The randomized range finder draws a random test matrix $S$6, forms

$S$7

and computes an orthonormal basis $S$8 for $S$9. The approximation

$\epsilon$0

is refined by computing a small SVD of $\epsilon$1. With target rank $\epsilon$2 and oversampling $\epsilon$3, $\epsilon$4 usually has $\epsilon$5 columns. Oversampling improves rank preservation and conditioning; power iteration improves accuracy when singular values decay slowly.

The structural error bound separates the optimal rank-$\epsilon$6 residual from a sketch-dependent term:

$\epsilon$7

Randomized SVD is well suited to BLAS-3 operations, GPUs, distributed systems, and matrices whose numerical rank is substantially smaller than their dimensions (2002.01387, Drineas et al., 2017).

Power iteration replaces the test matrix by a powered sample involving $\epsilon$8 or $\epsilon$9. It amplifies spectral gaps but requires additional matrix passes and intermediate orthogonalization. Block Krylov methods can converge faster than power iteration for clustered spectra, at the cost of additional storage and orthogonalization (2002.01387).

Matrix functions and trace estimation

For a positive semidefinite matrix $\operatorname{range}(U)$0, randomized trace estimation uses random vectors $\operatorname{range}(U)$1 satisfying

$\operatorname{range}(U)$2

and estimates

$\operatorname{range}(U)$3

Gaussian and Rademacher probes provide unbiased estimators with variance controlled by the Frobenius norm and, for positive semidefinite matrices, intrinsic dimension. Hutch++ combines low-rank approximation with residual sampling and improves the matrix-vector-query complexity for positive semidefinite trace estimation (2002.01387, Murray et al., 2023).

Stochastic Lanczos quadrature combines random trace estimation with Lanczos iterations to approximate $\operatorname{range}(U)$4. It is used for log determinants, traces of inverses, Schatten norms, and matrix-function quantities.

Randomized polynomial approximations provide a specific application to von Neumann entropy,

$\operatorname{range}(U)$5

for large density matrices. Taylor and Chebyshev approximations use randomized trace estimation and randomized power iteration, while a random-projection method approximates the nonzero eigenvalues of low-rank density matrices. The polynomial methods require a positive lower spectral bound for their stated relative-error theorems; the low-rank projection method avoids that assumption but requires low rank (Kontopoulou et al., 2018).

Generalized and nonstandard tasks

RNLA has also been applied to positive-semidefinite projection. A symmetric matrix $\operatorname{range}(U)$6 can be approximately projected onto the PSD cone by first computing a randomized low-rank approximation, then eigendecomposing a small compressed matrix and truncating its negative eigenvalues. A scaled variant shifts and rescales $\operatorname{range}(U)$7 so that positive eigenvalues are preferentially captured by the randomized range finder. The resulting projection can reduce the dominant cost from $\operatorname{range}(U)$8 to approximately $\operatorname{range}(U)$9 for target rank $\left\|I-(SU)^T(SU)\right\|_2\leq \epsilon.$0, supporting first-order semidefinite and sum-of-squares optimization (Jones et al., 2024).

Fast Randomized Iteration (FRI) extends RNLA in a different direction by repeatedly compressing the iterate itself. Given

$\left\|I-(SU)^T(SU)\right\|_2\leq \epsilon.$1

FRI uses

$\left\|I-(SU)^T(SU)\right\|_2\leq \epsilon.$2

where $\left\|I-(SU)^T(SU)\right\|_2\leq \epsilon.$3 has only $\left\|I-(SU)^T(SU)\right\|_2\leq \epsilon.$4 nonzero entries. This targets problems in which the matrix, matrix–vector product, or even the solution vector cannot be stored. FRI has been applied to eigenvalue problems, linear systems, matrix exponentiation, Fokker–Planck equations, diffusion problems, and formal dimensions exceeding $\left\|I-(SU)^T(SU)\right\|_2\leq \epsilon.$5 (Lim et al., 2015).

5. Applications and specialized extensions

RNLA is used in applications where large matrices are sparse, implicitly represented, streamed, distributed, or numerically low rank.

In machine learning, RNLA supports randomized PCA, partial SVD, kernel ridge regression, random features, Hessian approximation, Newton Sketch, gradient compression, and low-rank model representations. Randomized sketches can reduce the cost of Hessian formation and trace estimation when only Hessian–vector products are available. Modern analyses use random matrix theory, algorithmic Gaussianization, inversion-bias analysis, and statistical modeling to explain behavior in proportional-dimensional regimes and iterative optimization (Dereziński et al., 2024).

Kernel matrices benefit from coordinate Nyström approximation, ridge-leverage sampling, random features, and hierarchical randomized compression. Nyström methods approximate a positive-semidefinite kernel matrix using a small number of sampled columns, while ridge leverage scores adapt the sampling to regularization and effective dimension (2002.01387).

In graphs, RNLA provides spectral sparsification, approximate Laplacian solvers, graph embeddings, clustering, PageRank-like computations, and randomized trace estimators. Effective-resistance sampling produces sparse spectral approximations, while sparse approximate Cholesky methods provide preconditioners for graph Laplacian systems (2002.01387).

In scientific computing, applications include computed tomography, hyperspectral unmixing, genome-wide association studies, linear mixed-effect models, SNP principal-component analysis, operator inference, and data assimilation. Randomized Kaczmarz changes the access pattern of CT reconstruction; randomized SVD and randomized block Krylov methods extract dominant genetic components; randomized trace estimators approximate coefficients involving genetic covariance matrices; randomized low-rank approximations construct reduced-order models; and randomized preconditioners accelerate data-assimilation systems (Patel et al., 19 Jun 2025).

PDE-based inverse problems introduce an infinite-dimensional setting in which both the parameter space and experimental design space may be continuous. After linearization, the Jacobian becomes a design matrix or operator, and its Gram matrix is the Fisher information matrix. Leverage-score sampling, sketching, tensor-structured sketching, randomized Hessian or Fisher-information approximation, column subset selection, randomized SVD, compressed sensing, and matrix completion can therefore be used for data selection. PDE-specific source–detector factorizations can produce tensor or Khatri–Rao structure, motivating structure-aware RNLA methods (Hellmuth et al., 2 Oct 2025).

Specialized hardware has also been investigated. Photonic co-processors such as LightOn Optical Processing Units implement dense random projections through optical scattering. The OPU is used to accelerate the randomization stage of approximate matrix multiplication, trace estimation, triangle estimation, and randomized SVD, while conventional processors perform orthogonalization, reduced matrix products, and small factorizations. The reported OPU execution time is approximately $\left\|I-(SU)^T(SU)\right\|_2\leq \epsilon.$6 ms for the cited technology, with support up to roughly one million input dimensions and two million output dimensions; the paper reports approximately $\left\|I-(SU)^T(SU)\right\|_2\leq \epsilon.$7 tera-operations per second at $\left\|I-(SU)^T(SU)\right\|_2\leq \epsilon.$8 W (Hesslow et al., 2021).

RNLA also operates in extreme memory-avoidance regimes. Quantum-inspired algorithms for recommendation and least-squares regression can be interpreted as leverage-score or ridge-leverage-score sampling implemented through dynamic length-squared access structures. These structures support row sampling, conditional entry sampling, and entry queries under dynamic updates. The resulting algorithms use subspace embeddings, approximate matrix multiplication, projection-cost-preserving sketches, and sketch-and-solve methods (Chepurko et al., 2020).

6. Software, implementation, limitations, and future directions

RNLA software is increasingly organized into reusable sketching primitives and higher-level numerical drivers. The RandBLAS proposal targets data-oblivious dense sketching, sparse sketch application, multi-sketch operations, and quadratic sketches. RandLAPACK is intended to provide drivers for randomized least squares, preconditioning, low-rank SVD and eigendecomposition, CUR and interpolative decompositions, randomized QR, linear systems, trace estimation, leverage-score sampling, and tensor sketches (Murray et al., 2023).

Existing implementations and integrations include MATLAB and Python randomized low-rank routines, Blendenpik, LSRN, LibSkylark, Ski-LLS, RSVDPACK, Ristretto, Falkon, IMATE, MATLAB randomized methods, NAG Library routines, NVIDIA cuSOLVER randomized SVD functionality, and SciKit-Learn randomized SVD. The proposed software design emphasizes substitutable sketch distributions, factorization methods, iterative solvers, stopping criteria, and random-number generators.

Sketch ’n Solve provides a Python implementation of sketch-and-solve methods for large overdetermined least squares. It supports dense and sparse matrices, Gaussian and Hadamard-type sketches, Clarkson–Woodruff sketches, and solver families described as sketch-and-precondition and sketch-and-apply. In reported benchmarks on ill-conditioned synthetic matrices with $\left\|I-(SU)^T(SU)\right\|_2\leq \epsilon.$9, the package achieved speedups of up to $$(1-\epsilon)\|x\|_2^2 \leq \|SUx\|_2^2 \leq (1+\epsilon)\|x\|_2^2$$0 over SciPy’s LSQR while maintaining high reported accuracy. These are benchmark outcomes rather than universal performance or stability guarantees (Lavaee, 2024).

Randomness introduces reproducibility requirements. Counter-based random-number generators permit reproducible parallel sketch construction without storing all random coefficients. Nevertheless, bitwise reproducibility also depends on floating-point reduction order and the underlying BLAS implementation. Practical software should record random seeds or counter keys, sketch dimensions, sampling probabilities, precision, repetitions, stopping criteria, and fallback behavior.

RNLA is not universally advantageous. Important limitations include:

• Approximation requirements: exact factorizations or exact solutions remain necessary in some applications.
• Input-access costs: sketching cannot remove the cost of reading data when input movement dominates.
• Conditioning: residual accuracy may not imply solution accuracy for ill-conditioned systems.
• Spectral decay: low-rank methods lose their advantage for slowly decaying or flat spectra.
• Failure probability: many guarantees are probabilistic or expectation-based rather than deterministic.
• Sketch quality: undersized or poorly matched sketches can lose rank or produce weak preconditioners.
• Data-aware preprocessing: exact leverage scores may be nearly as expensive as deterministic factorizations.
• Single-pass constraints: one-pass methods usually have weaker guarantees and limited opportunities for correction.
• Finite precision: normal equations, Cholesky QR, low precision, and loss of orthogonality require explicit safeguards.
• Application structure: generic sketches may not preserve nonnegativity, sparsity, physical constraints, conservation laws, or tensor structure.
• DPP complexity: DPP sampling provides diversity and exact expected-error identities but is more expensive than independent sampling.
• FRI limitations: repeated compression can amplify instability, general signed or complex matrices are more difficult, and full-vector accuracy is generally unavailable.
• Hardware constraints: optical and accelerator approaches remain sensitive to encoding, readout, synchronization, and post-processing costs.

The field’s future directions emphasize structure-aware algorithms, hardware and mixed-precision co-design, and modular software infrastructure. Application-driven work calls for sketches that preserve sparsity, nonnegativity, Hermitian or positive-semidefinite structure, physical constraints, temporal or spatial structure, and tensor-product geometry (Patel et al., 19 Jun 2025). Modern theory seeks sharper explanations of practical sketch dimensions, inversion bias, proportional-dimensional behavior, and algorithmic Gaussianization (Dereziński et al., 2024).

Adaptive algorithms are an important direction. Desired capabilities include choosing the sketch dimension, rank, number of trace probes, power-iteration count, precision, and stopping threshold from an application-level tolerance and a posteriori error estimate. Independent residual sketches, randomized norm estimators, bootstrap methods, and posterior backward-error estimators provide mechanisms for certification. The resulting design pattern is:

$$(1-\epsilon)\|x\|_2^2 \leq \|SUx\|_2^2 \leq (1+\epsilon)\|x\|_2^2$$1

RNLA therefore constitutes a complementary layer of numerical linear algebra rather than a universal replacement for deterministic methods. Its principal contribution is to reorganize computation around limited information access, reduced communication, memory avoidance, parallelism, and controlled approximation. At moderate scale, classical QR, SVD, Krylov methods, multigrid, and deterministic eigensolvers may remain preferable. At extreme scale—when the matrix cannot be stored, the solution vector cannot be represented, or only selected observables are needed—randomized sketches, stochastic iterates, structured samples, and randomized preconditioners can make otherwise inaccessible computations feasible.