PaperlayerView exact revision on arXiv

Source-linked AI summary

Randomized Numerical Linear Algebra: Foundations & Algorithms

Per-Gunnar Martinsson, Joel Tropp

arXiv:2002.01387v3math.NA

TL;DR

The survey addresses how probabilistic methods can extend numerical linear algebra to large-scale and constrained data settings while retaining rigorous analysis and practical relevance. It synthesizes randomized algorithms and foundations across approximation, embeddings, iterative methods, factorizations, streaming, and linear-system solvers. The survey reports mainstream adoption and practical speedups, while noting that some randomized estimators cannot deliver high-precision approximations and that algorithm choice depends on the computational setting.

  • Problem

    Large datasets and constrained data access motivate alternatives to classical numerical linear algebra, while randomized methods historically raised concerns about accuracy and reproducibility.

  • Method

    The survey develops and synthesizes probabilistic linear algebra techniques, their mathematical analyses, and practical algorithms including embeddings, rangefinders, factorizations, and sparse solvers.

  • Results

    Randomized algorithms have joined mainstream NLA, with SPARSECHOLESKY achieving near-optimal runtime and storage guarantees for graph Poisson problems.

  • Takeaways & Limitations

    Randomized preconditioning and iterative sketching are preferred for high-precision sketching, while sketch-and-solve is appropriate when data access is severely constrained.

  • Takeaways & Limitations

    Matrix-sampling estimators require sample counts increasing with ε^-2, making them unsuitable for high-precision approximations.

Abstract

from arXiv · show

This survey describes probabilistic algorithms for linear algebra computations, such as factorizing matrices and solving linear systems. It focuses on techniques that have a proven track record for real-world problem instances. The paper treats both the theoretical foundations of the subject and the practical computational issues. Topics covered include norm estimation; matrix approximation by sampling; structured and unstructured random embeddings; linear regression problems; low-rank approximation; subspace iteration and Krylov methods; error estimation and adaptivity; interpolatory and CUR factorizations; Nyström approximation of positive-semidefinite matrices; single view ("streaming") algorithms; full rank-revealing factorizations; solvers for linear systems; and approximation of kernel matrices that arise in machine learning and in scientific computing.

1. INTRODUCTION

Randomized numerical linear algebra addresses scale, data-access, and hardware constraints that limit classical methods. This survey presents probabilistic algorithms, their mathematical foundations, and evidence of practical effectiveness.

  • The survey covers probabilistic techniques with mathematical analyses and a proven track record in practical computations.
  • Classical numerical linear algebra handles many standard problems, including linear systems, least squares, eigenvalue computations, and SVD.
  • Datasets at terabyte scale can exceed classical algorithms, while distributed, streaming, and out-of-core data require algorithms designed to minimize data movement.
  • Randomized algorithms moved from niche research into mainstream NLA, with methods appearing in textbooks and standard software libraries.
  • Randomized methods can reduce costs for selected problems, including least squares and Poisson systems, relative to classical methods.The survey cites O(mn + n3) versus O(mn2) for dense least squares and roughly quadratic versus cubic time for Poisson problems.
  • Randomization also supports single-pass truncated SVD with O(r(m + n)) working storage and reorganizes computation around parallelizable matrix–matrix multiplications.

4. TRACE ESTIMATION BY SAMPLING

Trace estimation uses randomized matrix–vector products to estimate a scalar trace, with variance reduction, probability guarantees, and extensions to matrix norms. The methods support streaming and parallel computation while offering empirical and a priori error assessment.

  • Core estimator: Randomized trace estimators construct unbiased samples and average independent copies to reduce variance.The sample mean remains unbiased, while variance decreases through independence.
  • Access model: Matrix–vector-product access makes randomized trace estimation useful when reading all diagonal entries or applying the primitive n times is expensive.The data-collection formulation uses Y = AΩ and can operate in one-pass or streaming settings.
  • Test-vector distributions: Rademacher sampling achieves the minimum variance among isotropic independent-coordinate estimators and avoids simulating normal variables.The comparison concerns the estimator generated by isotropic random vectors with statistically independent coordinates.
  • Error assessment: Trace estimators provide a priori probability bounds and a posteriori confidence intervals, with a practical rule of thumb k ≥30 and α ≥0.025.The interval is interpreted as containing trace(A) with probability roughly 1−2α over estimator randomness.
  • Practical design: Empirically, large-sample performance depends on the test-vector distribution through sample variance, a universality that permits convenient computational choices.The survey uses this phenomenon for heuristic error estimates and insights across distributions.
  • Norm extensions: Randomized trace estimation also yields unbiased estimators for the squared Frobenius norm and fourth power of the Schatten 4-norm.The resulting norm guarantees and error estimates can be expressed using stable rank, while structured test vectors can improve efficiency.

5. SCHATTEN p-NORM ESTIMATION BY SAMPLING

Randomized sampling can estimate Schatten 2p-norms from matrix–vector information, but higher-order estimators become costly or require large samples. For general matrices, the sample complexity grows polynomially with dimension when p > 2.

  • Scope and construction: Randomized sampling constructs unbiased estimators for the Schatten 2p-norm for every natural number p.For p ≥3, the estimators are combinatorial and may require many samples to control variance.
  • Sample complexity: k ≳ min{m, n}1−2/p samples are required for a 75%-probability constant-factor estimate of a general Schatten 2p-norm.Thus, for p > 2, the sample size grows polynomially with matrix dimension.
  • Scope boundary: Matrices with decaying singular values may be easier to handle than worst-case inputs, although the general theoretical bound can be pessimistic for such matrices.The survey identifies these spectral-decay cases as an area where smaller samples may suffice and questions remain open.
  • Computational cost: The direct U-statistic approach is computationally expensive, requiring O(min{m, n}p) operations in the worst case.Its statistic involves almost kp summands, although small fixed p can remain manageable.
  • Alternative estimator: The Kong–Valiant estimator sacrifices variance for computational efficiency, allowing larger p and larger sample sizes.Its increasing-index restriction permits a linear-algebraic reformulation, with cost often dominated by O(k2n) arithmetic.
  • Connection to spectral norm: Schatten 2p-norm estimates can approximate the spectral norm for sufficiently large p, but resampling is recommended because estimator reliability requires caution.The norms become equivalent when p is sufficiently large, with smaller p potentially adequate under spectral decay.

6. MAXIMUM EIGENVALUES AND TRACE FUNCTIONS

Randomized power and Krylov methods estimate maximum eigenvalues from matrix–vector products, while Lanczos-based techniques extend the framework to trace functions. Krylov methods converge faster but can require storing more data for eigenvectors.

  • Computational model: Sequential matrix–vector products make these iterative methods resistant to parallelization and unsuitable for one-pass or streaming environments.This contrasts with randomized methods based on linear information from a sample matrix.
  • Randomized power method: The randomized power method estimates the maximum eigenvalue by repeatedly applying a matrix to a random starting vector.It uses q sequential matrix–vector products and O(n) storage.
  • Power-method convergence: The power method needs a burn-in of q ≈log n iterations, after which error decreases with iteration count and faster when the spectral gap is bounded away from zero.The burn-in is necessary for any algorithm using q matrix–vector products in this setting.
  • Small-gap behavior: A reliable maximum-eigenvalue estimate remains possible when the relative spectral gap γ ≈0, although a reliable maximum-eigenvector estimate is generally not.The eigenvalue result therefore extends beyond analyses that require a separated leading eigenvalue.
  • Krylov acceleration: Krylov methods improve the post-burn-in rate to 1/q2, compared with 1/q for the power method, and reduce the per-iteration gap factor from e−γ to e−2√γ.Lanczos is an efficient Krylov implementation for self-adjoint matrices.
  • Resource tradeoff: Krylov eigenvalue estimation costs almost the same as the randomized power method when no eigenvector is needed, but storing a Krylov basis creates a time–data tradeoff.Eigenvector estimation commonly requires retaining a basis for Kq.
  • Trace functions: Stochastic Lanczos quadrature combines randomized trace estimation with Lanczos iteration to address traces of spectral functions.The survey presents trace-of-spectral-function computation as a broadly applicable problem.

7. MATRIX APPROXIMATION BY SAMPLING

Matrix approximation by sampling constructs structured unbiased estimators from simple summands and improves them by averaging independent samples. Matrix Bernstein analysis quantifies sample requirements, while applications include approximate multiplication and spectral graph sparsification.

  • Empirical approximation: Empirical approximation expresses a target matrix as a sum of simple, often sparse or low-rank, matrices and samples rescaled summands to form an unbiased structured estimator.A single sample preserves summand structure but is rarely accurate, motivating averages of independent copies.
  • Empirical approximation: Averaging independent estimators preserves unbiasedness and improves approximation quality.The sample average remains unbiased by linearity of expectation.
  • Matrix Monte Carlo: The matrix Bernstein bound requires a sample count proportional to the larger of the per-sample second moment v(X) and norm bound R.The required count also grows with ε^-2 and often with log(m + n).
  • Approximate multiplication: Approximate matrix multiplication samples rank-one column-row products, reducing explicit estimator cost to O(mnk) when k ≪I.The estimator is unbiased for BC, whereas direct multiplication costs O(mnI).
  • Approximate multiplication: Importance sampling needs a number of rank-one factors proportional to the maximum coherence, or to log(m + n) times total stable rank under the corresponding bound.Its sample-complexity bound always improves over uniform sampling under the stated normalization.
  • Graph sparsification: Random graph sampling produces an ε-spectral approximation with at most k nonzero weights; every n-vertex graph has a (1/2)-approximation with at most 12n log n weights.Related sparsification ideas also support practical fast Poisson solvers.

8. RANDOMIZED EMBEDDINGS

Randomized embeddings reduce dimension while approximately preserving the geometry of a target set. Gaussian-width analysis links the required embedding dimension to set geometry, while practical bounds and computational costs vary across embeddings and sets.

  • Embedding framework: An ℓ2 embedding maps a set in a high-dimensional space to a lower-dimensional space while preserving norms within distortion ε.The embedding should preserve the set’s geometry, but knowing the set, constructing the map, and applying it efficiently can conflict.
  • Gaussian embeddings: Gaussian embeddings provide a tractable theoretical model whose restricted singular values are controlled by the Gaussian width of the embedded set.The squared Gaussian width can be viewed as a measure of the set’s effective dimension.
  • Gaussian embeddings: Choosing d > (w(E) + 1)^2 makes a Gaussian embedding unlikely to annihilate points in E and unlikely to expand them by more than a factor of two.This reduces embedding-dimension selection to computing Gaussian widths.
  • Finite point sets: For a finite point set of cardinality N, the embedding dimension needed for distortion ε is logarithmic in N.The analysis also yields precise failure probabilities with additional calculation.
  • Limitations: The Johnson–Lindenstrauss bound can be impractical for small distortion because its required dimension scales with ε^-2 and may be prohibitively large.The uniform bound may remain large even when ε ≈1.
  • Subspace embeddings: For a k-dimensional subspace, d = 2k gives distortion 0.8 with failure probability e^-ck, and some applications use d = k + 5 or d = k + 10.These examples illustrate dimensions close to the subspace dimension.

9. STRUCTURED RANDOM EMBEDDINGS

Structured random embeddings seek Gaussian-like practical behavior with lower construction, storage, and application costs. Sparse maps and SRTTs offer theoretical subspace-embedding guarantees, while uniform and leverage-score sampling involve distinct coherence, cost, and adaptivity trade-offs.

  • Motivation: Gaussian embeddings and partial isometries perform well but can be expensive to construct, store, and apply, motivating structured alternatives.Structured constructions are emphasized for applications and often have similar practical performance.
  • Mixing and preprocessing: Structured embeddings commonly mix coordinates before sampling so that energy is homogenized and the sampled norm has low variance.Random sign flips, permutations, and restrictions support preconditioning and dimension reduction.
  • Sparse embeddings: Sparse sign matrices use about O(ζn log d) storage and O(ζn) application work, while retaining performance similar to Gaussian embeddings.The construction places ζ random signs in uniformly selected coordinates of each column.
  • Sparse embeddings: A sparse sign matrix is an oblivious constant-distortion embedding for any k-dimensional subspace when d = O(k log k) and ζ = O(log k).The result is attributed to Cohen (2016), with further improvements conjectured.
  • SRTTs: An SRTT combines a trigonometric transform, random sign flip, permutation, and restriction, and is a partial isometry.Its construction is designed to mimic a random partial isometry.
  • Sampling trade-offs: SRTTs provide oblivious constant-distortion embeddings for k-dimensional subspaces with d = O(k log k), although d = O(k) often suffices in practice without rigorous justification.Uniform sampling is attractive when coherence is small, but coherence is usually unavailable for choosing d in advance.
  • Sampling trade-offs: Leverage-score sampling can be costly and is rarely competitive for subspace embeddings, whereas coordinate sampling remains tractable for some kernel computations.Uniform sampling is a workhorse for large-scale kernel computation and may outperform expectations in practice.

10. HOW TO USE RANDOM EMBEDDINGS

Random embeddings support several least-squares algorithm templates that trade speed, precision, matrix views, and repeated work. Sketch-and-solve is fast but inaccurate, while iterative sketching and preconditioning recover high-precision behavior in suitable regimes.

  • Sketch and solve: Sketch-and-solve compresses an overdetermined least-squares problem and uses the reduced solution as a proxy, requiring one matrix view but often producing inaccurate results.For dense unstructured A, forming SA costs O(mn log d) and solving the reduced problem costs O(dn^2).
  • Sketch and solve: Sketch-and-solve improves computational cost when log n ≪n ≪m/log n and ε is constant, but it must accept large errors.The approach cannot make ε small in this regime.
  • Iterative sketching: Iterative sketching repeatedly sketches residual least-squares problems to reduce the error left by sketch-and-solve.Each iteration uses a fresh constant-distortion embedding and feeds its solution into the next subproblem.
  • Iterative sketching: Iterative sketching achieves relative error ε at O((mn + n^3) log(n) log(1/ε)) cost and can reach errors as small as traditional algorithms.It costs slightly more than sketch-and-solve for constant error and requires repeated sketches of A.
  • Sketch and precondition: Sketch-and-precondition uses a random proxy for A to precondition a classical iterative method and reduce its iteration count.For least-squares, it can be faster than both other sketching paradigms while maintaining comparable errors.
  • Sketch and precondition: In the regime log n ≪n ≪m/log n, sketch-and-precondition may be a factor of n/log(n) faster than classical methods while supporting very small ε.It requires repeated applications of A.
  • Practical guidance: For high-precision dense least-squares, randomized preconditioning is identified as the most efficient existing approach, while sketch-and-solve suits severely constrained data access.Iterative sketching remains useful for constrained least-squares problems.

11. THE RANDOMIZED RANGEFINDER

The randomized rangefinder seeks a subspace aligned with a matrix’s dominant left singular vectors, using random embeddings and orthogonalization. Its accuracy can match the best lower-dimensional approximation under suitable spectral conditions, while powering and Krylov variants improve performance.

  • Core algorithm: The randomized rangefinder multiplies B by a random embedding and orthogonalizes the resulting sample matrix.The survey analyzes Gaussian, structured, and other mixing embeddings, whose performance depends on more than basic subspace embedding properties.
  • Core problem: The rangefinder produces an orthonormal Q whose range aligns with the dominant left singular vectors of B.The approximation QQ∗B can then serve as a rank-ℓ proxy for B when the spectral-norm error is small.
  • RSVD consequences: The approximate SVD’s accuracy is determined entirely by the rangefinder error because subsequent computations are exact apart from floating-point errors.The smallest computed singular values and vectors may not estimate the true ones well, so truncation can be valuable.
  • Guarantees: The randomized rangefinder computes an ℓ-dimensional subspace whose captured action can approach that of the best k-dimensional subspace when ℓ is only slightly larger than k.The error is comparable with σk+1 when the tail singular values have small ℓ2 norm, as with rapidly decaying spectra.
  • Practical guidance: A posteriori error estimates are recommended to validate rangefinder performance, especially when structured embeddings have weaker a priori guarantees.Random coordinate sampling is not recommended unless the described random embeddings cannot be used.
  • Powered methods: Powering drives the rangefinder error toward σk+1 exponentially as q increases, with q = 2 or q = 3 often sufficient when the spectrum decays.For a flat spectral tail, the survey says q ≈log min{m, n} may be needed for a constant-factor error.
  • Krylov methods: Krylov rangefinders reach error O(ε) with depth q ≈log(n/ε)/√ε, versus about q ≈(log n)/ε iterations for the power method.The difference can be substantial when ε is small, and Krylov methods also support accurate singular-value and eigenvalue computation.

12. ERROR ESTIMATION AND ADAPTIVITY

This section addresses the practical problem of estimating approximation error and adapting rank when singular values or numerical rank are unknown. It uses inexpensive auxiliary random samples to certify accuracy and guide refinement.

  • Problem setting: Without advance knowledge of singular values, fixed-error approximation requires finding a low-rank factorization within a requested tolerance and determining the ε-rank.The rangefinder therefore needs an a posteriori estimator for ∥A − QQ∗A∥.
  • Adaptivity: If the estimate exceeds the requested tolerance, additional samples can enrich the existing basis or the computation can restart with a different strategy.Incremental basis construction is designed for cases where the rank is unknown in advance.
  • Error estimation: Error estimation draws a small independent Gaussian auxiliary sample of the residual and applies norm-estimation techniques to that sample.The framework supports spectral, Schatten 4-norm, and higher-order Schatten-norm estimates.
  • Error estimation: The auxiliary sample usually costs only s additional matrix–vector multiplications, with s suggested as a small fixed number such as s = 10.This cost is almost always much smaller than running the rangefinder itself.
  • Certification: Pairing fast structured random maps with a small Gaussian auxiliary sample can retain computational efficiency while providing strong confidence in the computed output.The estimator can also certify heuristic sampling strategies when their assumptions about the matrix are uncertain.
  • Adaptivity: Structured random matrices may make recycling existing samples difficult, so rank doubling from scratch is often sensible when a tolerance is missed.The recommendation is presented as a practical response to the structure of these test matrices.
  • Gaussian implementation: For Gaussian rangefinders, the test matrix can be split so that a thin block acts as an independent error-estimation sample.The resulting Q2R22 block samples the residual associated with the first test block, and a small estimate permits trusting the factorization.

13. FINDING NATURAL BASES: QR, ID, AND CUR

This section develops interpolative and CUR factorizations that preserve selected rows or columns of a matrix, trading orthogonality for interpretability and structural properties. It also describes randomized algorithms that compute near-optimal IDs and reduce computational cost.

  • Natural-basis factorizations use subsets of matrix columns or rows, preserving properties such as sparsity, nonnegativity, and interpretability.
  • Interpolative decompositions: The column ID A = CZ uses selected columns C and a coefficient matrix Z containing the identity as a submatrix; Z can be chosen with entries bounded by one in modulus.
  • Interpolative decompositions: The row ID analogously writes A = XR using selected rows R and an identity-containing matrix X; selected row and column indices are called skeleton index vectors.
  • CUR decomposition: The CUR decomposition A = CUR combines selected columns, selected rows, and a small linking matrix, often requiring less storage than dense ID factors, especially for sparse matrices.
  • Approximate rank: Rapid singular-value decay makes approximate ID ranks close to SVD ranks, whereas slow decay can create a substantial rank gap and further algorithmic suboptimality.
  • Randomized algorithms: Randomized rangefinding supports ID construction in O(mnk) time with Gaussian tests, or O(mn log k) time with structured random matrices.

14. NYSTR ¨OM APPROXIMATION

The randomized Nyström method constructs low-rank positive-semidefinite approximations from random samples, adapting the rangefinder framework. Its error approaches the optimal rank-k error and improves with more samples or powering.

  • The goal is a rank-k PSD approximation of A whose error is nearly that of the best rank-k approximation.
  • Construction: Nyström approximations are formed from a test matrix X, depend only on its sampled range, and have error equal to A's Schur complement with respect to X.
  • Construction: Randomized Nyström draws Ω, forms Y = AΩ, computes a stable factored approximation, and truncates its initial rank ℓ approximation to rank k.
  • Streaming: Because Algorithm 16 uses only linear information about A, it supports one-pass or streaming implementations and applications such as kernel PCA.
  • Guarantees: The Gaussian analysis shows that the rank-k approximation achieves nearly the optimal error λ_k+1, with error decreasing as the sample count ℓ increases and the eigenvalue tail shrinks.
  • Powering: Powering reduces Nyström error exponentially as q increases, but requires careful numerical implementation.

15. SINGLE-VIEW ALGORITHMS

Single-view algorithms approximate matrix decompositions while each matrix entry is accessed only once, addressing matrices that are too large to store or revisit. They use simultaneous randomized sketches, with accuracy and memory trade-offs compared with multi-pass methods.

  • Single-view SVD targets matrices that cannot be stored and permits each entry of A to be viewed only once, unlike basic randomized SVD.
  • Single-view SVD: The algorithm samples row and column information simultaneously, orthonormalizes the resulting sketches, and estimates the core matrix through least squares.
  • Single-view SVD: The resulting core estimate uses pseudoinverse-based least squares and can be postprocessed into an approximate SVD.
  • Streaming algorithms: Streaming algorithms retain an evolving low-dimensional linear sketch of matrix updates, discarding each innovation after processing.
  • Trade-offs: Single-view algorithms can reduce storage and data passes, but they inherently increase uncertainty and decrease accuracy; revisiting the matrix once is advisable when feasible.
  • Guarantees and implementation: For Gaussian tests with s ≥ 2ℓ, the theory provides an error guarantee relative to the best rank-k approximation, while structured tests can greatly reduce storage with similar observed errors.

16. FACTORING MATRICES OF FULL OR NEARLY FULL RANK

This section develops randomized full- and near-full-rank factorizations that preserve rank-revealing quality while improving practical speed through blocking and matrix–matrix multiplication.

  • 16. FACTORING MATRICES OF FULL OR NEARLY FULL RANK: Randomization reorganizes most arithmetic into matrix–matrix multiplications, reducing communication costs even when the asymptotic cost remains O(n3).Stable fast multiplication can further reduce the overall cost to O(nω).
  • 16. FACTORING MATRICES OF FULL OR NEARLY FULL RANK: Rank-revealing factorizations approximate the best rank-k truncation of A for every k, generalizing both the SVD and CPQR.They also support ill-conditioned solves, least squares, singular-spectrum estimation, and bases for fundamental subspaces.
  • 16.3. The powerURV algorithm: PowerURV is faster than CPQR in practice despite requiring approximately 5n3 flops versus 0.5 n3 flops when m = n and q = 2.The same passage reports that powerURV with q = 2 reveals numerical rank substantially better than CPQR.
  • 16.6. A randomized algorithm for computing a CPQR decomposition: Randomized blocked CPQR selects pivot columns whose quality is almost indistinguishable from traditional pivoting, with relation (16.2) holding to about the same accuracy.The diagonal entries of R need not strictly decay across block boundaries.
  • 16.7. A randomized algorithm for computing a URV decomposition: Randomized URV is nearly as fast as randomized CPQR while approximating numerical rank almost as well as the SVD and providing bases for both column and row spaces.Its R factor has small off-diagonal entries and diagonal entries that closely approximate the singular values of A.
  • 16.7. A randomized algorithm for computing a URV decomposition: The randomized URV process is incremental and can stop once a requested accuracy is reached, retaining practical advantages of blocked randomized CPQR.This makes it suitable when computation should terminate at a specified tolerance.

17. GENERAL LINEAR SOLVERS

This section surveys probabilistic approaches to linear systems, emphasizing their dependence on matrix structure and noting that randomized solvers for square systems remain comparatively preliminary.

  • 17. GENERAL LINEAR SOLVERS: Randomized linear-solver research is more preliminary than other parts of the survey, with attention mainly to square consistent systems and some regression problems.The section frames solver choice around matrix properties such as density, storage, and available matrix–vector access.
  • 17. GENERAL LINEAR SOLVERS: Well-conditioned or spectrally clustered systems are already solved rapidly by Krylov methods such as CG and GMRES, often at the cost of only a handful of matrix–vector multiplications.This makes randomized methods difficult to improve upon in such settings.
  • 17. GENERAL LINEAR SOLVERS: Random unitary left and right preconditioning enables blocked Gaussian elimination without pivoting, which is almost surely free of degenerate diagonal blocks.Blocking also reduces communication and improves the efficiency of matrix operations on modern computers.
  • 17. GENERAL LINEAR SOLVERS: Randomized rank-revealing factorizations can stabilize ill-conditioned solves by avoiding subspaces associated with small singular values, while randomized URV offers faster stabilization than an SVD.Truncating modes with σj ≤ ε permits a nonzero residual and suppresses large solution components that contribute little to satisfying the equation.
  • 17.3. Sketch and precondition: For square systems, randomized preconditioning has been less fruitful; useful techniques exist in special cases, while some spectral configurations remain open research problems.The survey specifically identifies difficult positive-definite spectra for which finding randomized preconditioners remains open.
  • 17.3. Sketch and precondition: For positive-definite systems, CG converges rapidly when eigenvalues are clustered, so preconditioning seeks a rapidly applicable M^-1 making M^-1/2AM^-1/2 tightly clustered.Randomized low-rank approximation provides excellent preconditioners when A has a few eigenvalues larger than the others.

18. LINEAR SOLVERS FOR GRAPH LAPLACIANS

SPARSECHOLESKY uses randomized incomplete Cholesky ideas to precondition graph-Laplacian systems, combining sparse clique sampling with PCG. The method offers near-optimal complexity while retaining a practical limitation in its initialization cost.

  • Method: SPARSECHOLESKY constructs a sparse approximate Cholesky factor for a graph Laplacian and uses it as a preconditioner for PCG.The approach is a randomized variant of incomplete Cholesky factorization.
  • Method: κ ≤3 when the approximate factor satisfies the required spectral relation, giving the preconditioned system a bounded condition number.This conditioning supports rapid convergence of PCG.
  • Performance: 1 + log3(1/ε) PCG iterations suffice for relative error ε, with O(m log2 n) arithmetic operations per iteration.Each iteration applies the Laplacian and solves a system involving CC∗.
  • Performance: O(m log2 n) nonzeros and O(m log3 n) expected operations characterize the produced factor and algorithm.The theorem states these guarantees with high probability for connected weighted graphs.
  • Method: Random clique sampling is unbiased, controls variance through bounded effective resistance, and avoids constructing the full quadratic-size clique.The sampled clique has no more multiedges than the eliminated star.
  • Limitations: The initialization splits each multiedge into O(log2 n) pieces, making the method uncompetitive for some instances; whether this step can be relaxed remains unclear.This is the stated main shortcoming of SPARSECHOLESKY.

19. KERNEL MATRICES IN MACHINE LEARNING

Kernel methods represent data-point similarities with positive-definite kernel matrices, but kernel evaluations can be expensive at scale. The survey presents Nyström coordinate sampling and random features as complementary approximation tools, with guarantees and costs tied to effective dimensions and feature counts.

  • Scope: The survey does not comprehensively cover kernel methods, learning theory, or broader computational considerations.Its goal is a brief introduction and a small selection of key references.
  • Kernel foundations: Positive-definite kernels encode pairwise data similarity in kernel matrices that generalize Gram matrices to feature spaces.A feature map represents k(x,y) as an inner product of mapped data points.
  • Kernel foundations: Kernel methods can transfer algorithms formulated for Euclidean inner products to other positive-definite kernels through kernel-matrix substitution.This supports applications across domains including images, text, and DNA sequences.
  • Approximation methods: Nyström coordinate sampling and random features provide two independent randomized approaches for approximating kernel matrices while controlling costly kernel evaluations.The survey treats coordinate sampling and random features in separate subsections.
  • Nyström approximation: Uniform Nyström sampling achieves optimal KRR learning guarantees when the coordinate count is proportional to ν_m log n.A corresponding result is also stated for KPCA.
  • Nyström approximation: Ridge leverage score sampling achieves optimal KRR learning guarantees with r proportional to ν_eff log n, improving over the uniform-sampling bound.Related theoretical results are stated for KPCA.
  • Random features: Random features form an unbiased kernel estimator, with O(rnd) construction cost versus O(n2d) for explicit kernel formation when r ≪n.The approximation can achieve relative spectral-norm error when r scales with the energetic dimensions of K.

20. HIGH-ACCURACY APPROXIMATION OF KERNEL MATRICES

The section develops randomized, rank-structured approximations for kernel matrices when fast matrix-vector products are available, targeting high relative accuracy without forming the full matrix. It presents hierarchical compression schemes ranging from black-box HODLR methods to linear-complexity HBS methods and butterfly-related formats.

  • High-accuracy kernel approximation: 10^-3 or 10^-6 relative accuracy is targeted for kernel-matrix approximations when fast evaluations of Kx and K∗x are available.The high-accuracy regime seeks small relative error in approximating K.
  • Rank-structured representations: Slowly decaying singular values preclude low-rank approximation of the full kernel matrix, motivating O(n)-block hierarchical representations with low-rank off-diagonal blocks.These representations are called rank-structured hierarchical matrices.
  • Downstream operations: The resulting sparse representations can support matrix inversion, LU factorization, and, in some cases, full spectral decompositions.These operations motivate constructing a rank-structured approximation when fast matrix-vector products are already available.
  • Low-rank block structure: Smooth kernels on separated subdomains yield approximate separations of variables, which produce low-rank factorizations for off-diagonal blocks K(It, Is).The factors B and C are assembled from kernel-dependent functions evaluated at target and source points.
  • Compression challenge: O(n^2) cost arises from forming the full matrix before compressing its off-diagonal blocks, whereas randomized compression avoids that necessarily quadratic preprocessing.The section frames randomized compression as a response to the expense of explicitly forming all compressible blocks.
  • Randomized compression algorithms: The HODLR black-box technique uses O(log n) applications of K to random n×(r+p) matrices with O(n log n) storage, while an HBS method reaches O(n) complexity with one application each of K and K∗.The linear-complexity method additionally evaluates O(rn) individual entries and applies only to a smaller matrix class.
Loading 2002.01387v3…