Randomized Numerical Linear Algebra : A Perspective on the Field With an Eye to Software

Release History

• The first draft circulated as an SC22 technical report on 11/13/2022; arXiv V1 on 02/22/2023 incorporated feedback across Sections 3.2.2, 4, 5.2, and Appendix B.2.Source: Release History · p. 5
• ArXiv V2 on 04/12/2023 revised Section 5.1, rewrote and expanded Section 5.3, added software references in Sections 4.1.1 and 4.5, and clarified Section 5.2.2.Source: Release History · p. 5

Introduction

• The introduction motivates the subject, defines the monograph’s scope and structure, summarizes sections, offers a survey of surveys, and recommends consulting notation and definitions.Its three principal goals are motivation and clarifying misconceptions, explaining scope and structure, and guiding readers through section summaries.Source: Introduction · p. 8

1.1 Our world

• NLA’s scaling faces hardware limits that make parallelization essential and a mature field where useful problem structures and algorithmic breakthroughs are increasingly difficult to find.Even purpose-built accelerators require programming, and high-performance implementations of classical NLA algorithms are difficult.Source: §1.1 · p. 9
• RandNLA offers efficient, reliable approximate solutions and fine-grained control over the accuracy–cost tradeoff, with implementations adaptable to new hardware.Elementary MATLAB or Python implementations can demonstrate practicality, while high-performance implementations are easier to build with suitable building blocks.Source: §1.1 · p. 9

Introduction

• RandNLA algorithms use randomness to uncover structures hidden from classical NLA, without assuming a distribution over inputs or intrinsic input uncertainty.For highly overdetermined least squares, randomization can quickly identify a basis in which the matrix is nearly column-orthonormal, enabling a classical iterative method to solve reliably to error ϵ.Source: Introduction · p. 10
• The monograph addresses gaps in transferring RandNLA research into software by providing a principled, practical foundation for high-quality RandBLAS and RandLAPACK libraries.Many described algorithms have been implemented and tested in MATLAB and Python.Source: Introduction · pp. 10–11

1.1.1 Four value propositions of randomization

• For many important linear algebra problems, RandNLA offers new ways to compute approximate solutions with linear or nearly-linear complexity.This is presented as the first value proposition of randomization.Source: §1.1.1 · p. 11

1.1. Our world

• For n × n positive definite matrices, the number of free parameters is about 50 million at n = 10,000 and 5 billion at n = 100,000, making superlinear complexity crippling.The authors urge users to think about problem size in terms of the matrix’s free parameters.Source: §1.1 · p. 12
• Although sparse approximate Cholesky for graph Laplacians is a RandNLA success, this monograph omits sparse-only algorithms and includes approximate eigendecompositions for regularized positive semidefinite matrices that can speed some linear-system solves.Source: §1.1 · p. 12
• Randomized algorithms are asymptotically fastest among accurate methods for several important linear algebra problems, including highly overdetermined least squares.Source: §1.1 · p. 12
• Randomization can reduce and redirect data movement, making randomized algorithms significantly faster by wallclock time than the best available deterministic methods.The RAM model misses the cost of moving data through memory and cache, so algorithms with similar operation counts can differ greatly in wallclock time.Source: §1.1 · p. 12
• RandNLA can use lower-precision arithmetic on matrix sketches while preserving final accuracy and nearly avoiding degeneracies that can cause common deterministic algorithms to fail.Source: §1.1 · p. 13

1.1.2 What is, and isn’t, subject to randomness

• RandNLA algorithms generally use randomness only to sample from a sketching distribution, and specifying the random-number-generator seed makes them deterministic.Sketching operators are random linear dimension-reduction maps, ranging from row or column selection to more complex operators.Source: §1.1.2 · p. 13
• Sketch processing can yield rough or high-accuracy approximations, sometimes exact under mild assumptions, and analyses often bound the probability that a sketch has a property needed for the desired output.A common analysis first characterizes the output using a sketch property, then uses random matrix theory to bound the probability that the property holds.Source: §1.1.2 · p. 13

1.2. This monograph, from an astronaut’s-eye view

• Randomized algorithms generally expose parameters that let users control either accuracy or cost, and large-scale applications tend to become more predictable.When cost is controllable, parameters can influence accuracy; when accuracy is controllable, they can influence cost.Source: §1.2 · p. 14

1.2 This monograph, from an astronaut’s-eye view

• The monograph aggregates over 300 references and adds unifying taxonomies, candidate application areas, theoretical results, and algorithms while informing plans for RandBLAS and RandLAPACK.It is designed to be modular and accessible, with background in each section and technical material and implementation notes in appendices.Source: §1.2 · pp. 14–15

A word on “drivers” and “computational routines”

• Drivers solve higher-level traditional linear algebra problems using computational routines, but this design makes driver-level performance guarantees difficult to state without prescribing routine choices.The authors consequently make relatively few claims about driver-level performance guarantees or computational complexity.Source: A word on “drivers” and “computational… · p. 15

Basic Sketching (Section 2)

• The RandBLAS discussion proposes data-oblivious sketching functionality and covers sketching concepts, operator distributions, and elementary operations.Appendix A.1 introduces effective distortion for characterizing randomized algorithms in least squares and related problems.Source: Basic Sketching (Section 2) · pp. 16–17

Least Squares and Optimization (Section 3)

• For least squares and optimization, the monograph covers low-accuracy sketch-and-solve, high-accuracy random preconditioning, and novel low-memory preconditioners for regularized saddle point problems.It also reframes two approximate kernel ridge regression algorithms as sketch-and-solve methods for regularized problems.Source: Least Squares and Optimization… · p. 17, §1.3 · p. 18

Low-rank Approximation (Section 4)

• The low-rank approximation section clarifies which parts of the broad randomized-algorithm design space suit different situations and surveys drivers and their computational building blocks.It covers SVD, Hermitian eigendecomposition, CUR and two-sided interpolative decomposition drivers, with QB decomposition and column subset selection among the basic building blocks.Source: Low-rank Approximation (Section 4) · pp. 18–19

Further Possibilities for Drivers (Section 5)

• Further driver topics include randomized QR and SVD-surrogate decompositions, direct and iterative methods for unstructured linear systems, and randomized trace estimation.The trace discussion includes methods for estimating the trace of f(B) when B is Hermitian.Source: Further Possibilities for Drivers… · p. 19

Advanced Sketching: Leverage Score Sampling (Section 6); Advanced Sketching: Tensor Product Structures (Section 7)

• Advanced sketching covers leverage-score sampling and efficient sketching of Kronecker- and Khatri–Rao-structured matrices, including exact leverage-score sampling without explicitly forming tall matrices.The tensor-structure material also discusses motivating tensor-decomposition applications.Source: Advanced Sketching: Leverage Score… · p. 20, Advanced Sketching: Tensor Product… · p. 20

1.4.1 Tutorials, light on prerequisites

• Drineas and Mahoney [DM18] present representative RandNLA results with proofs for least squares and low-rank approximation, emphasizing separation of deterministic and probabilistic analysis.The resource is notable for brevity and self-contained introductions to linear algebra and probability.Source: §1.4.1 · p. 21

1.4.2 Broad and proof-heavy resources

• Woodruff develops RandNLA problem by problem, including lower bounds; Tropp introduces matrix concentration, while Mahoney [Mah16] covers theory and practice.Tropp’s course notes build on with subsequent theory and applications; Mahoney’s notes cover results through 2013, including concentration, sketching, least squares, and low-rank approximation.Source: §1.4.2 · p. 21

1.4.3 Perspectives on theory, light on proofs

• Mahoney [Mah11] emphasizes concepts, interpretations, and qualitative proof strategies to explain RandNLA’s theoretical contributions to least squares and low-rank approximation.It also connects RandNLA theory to developments in other disciplines.Source: §1.4.3 · p. 21

1.4. Recommended reading

• Dereziński and Mahoney present RandNLA theory through determinantal point processes and statistical data analysis, taking a prospective rather than retrospective view.Source: §1.4 · p. 22

1.4.4 Deep investigations of specific topics

• Halko, Martinsson, and Tropp offer an influential, thorough treatment of low-rank approximation, while Kannan and Vempala survey sampling and tensor methods.Martinsson [Mar18] covers practical advances in low-rank approximation and algorithms for full-rank matrix decomposition; is described as partly superseded by but retains proved results and some distinct algorithms.Source: §1.4.4 · p. 22

1.4.5 Randomized numerical linear algebra: Foundations and Algorithms

• Martinsson and Tropp’s survey spans topics with technical and historical depth and informed the authors’ plans for RandBLAS and RandLAPACK.The authors distinguish their own similarly long monograph because they see value in clarifying how it differs from.Source: §1.4.5 · p. 22

Basic sketching.

• Compared with, the monograph emphasizes implementation over theory, reviews literature relevant to sketching implementation, and adds novel technical contributions in its appendices.The comparison refers to its coverage of basic sketching.Source: Basic sketching · p. 22

Least squares and optimization.

• The optimization coverage is comprehensive insofar as optimization reduces to linear algebra, and includes novel technical contributions and relevant software, unlike ’s limited treatment.The passage points to Section 3 and Appendix B as the monograph’s coverage.Source: Least squares and optimization · p. 23

Low-rank approximation.

• This section provides broad low-rank approximation scaffolding and an overview of implementing software, but fewer opportunities for mathematical explanation than its approach to might allow.The authors describe their approach as very different from.Source: Low-rank approximation · p. 23

Full-rank matrix decompositions

• The monograph covers a broader range of full-rank decompositions and algorithms than, including one novel algorithm with proofs, while provides more mathematical detail.Source: Low-rank approximation · p. 23

Kernel methods

• For kernel methods, the monograph covers positive definite machine-learning kernels but not indefinite scientific-computing kernels, emphasizing the resulting linear algebra problems.Both topics are addressed in.Source: Low-rank approximation · p. 23

Linear system solvers

• The monograph covers slightly more material on unstructured linear-system solvers than, but excludes sparse-specific methods, including approximate Cholesky methods for sparse graph Laplacians.Source: Low-rank approximation · p. 23

Trace estimation

• Two recently developed trace-estimation methods are described as the first major advances since the late 1980s; coverage is generally less deep than, except for stochastic Lanczos quadrature.Source: Low-rank approximation · p. 23

Advanced sketching

• The monograph gives leverage-score sampling and tensor-product sketching ten pages each, treating them as complementary because implicit tensor-structured matrices are practical candidates for leverage-score sampling.Kernel matrices from machine learning are also described as strong practical candidates.Source: §1.5 · p. 24

1.5 Notation and terminology; Matrices and vectors

• Notation is summarized in Table 1.1; conventions include adjoints and pseudoinverses, real or complex matrices, column-oriented vectors, QR decompositions, and compact SVDs.The text defines Hermitian and positive semidefinite matrices, and also describes QR with column pivoting.Source: §1.5 · p. 24, Matrices and vectors · p. 24, Table 1.1 · p. 26

Probability and our usage of the term “random.”

• “Random” matrix operations are deterministic after construction; unqualified “random” for finite-outcome actions means uniform randomness over the possible actions.For abstract linear operators, randomness applies to constructing the data that defines the operator.Source: Probability and our usage of the term… · pp. 24–25

Basic Sketching

• The proposed RandBLAS is to RandNLA as BLAS is to classical NLA, initially supporting sketching operators represented by dense arrays or structured sparse matrices.Row and column sampling are treated as particular forms of sparse sketching.Source: Basic Sketching · pp. 28–29

2.1 A high-level plan

• RandBLAS should facilitate high-level RandNLA algorithms while initially limiting its scope to sketching dense data matrices.Data-aware operators, including power-iteration sketches and leverage-score sampling, are important but should largely be handled by higher-level libraries.Source: §2.1 · p. 29
• The design groups data-oblivious sketching into sampling an operator, applying it to a data matrix, and sketching not expressed as a single linear operator.These categories offer increasing opportunities for machine-specific optimization, but the authors do not advocate formal RandBLAS levels.Source: §2.1 · p. 29

2.1.1 Random number generation

• Counter-based random number generators are recommended for reproducibility because keyed counters allow flexible, embarrassingly parallel generation of sketch entries.For a dense d×m sketch, an entry can be generated using counter c+(i+dj), where c is a user-provided offset.Source: §2.1.1 · p. 30

Shift-register RNGs

• Shift-register RNGs can make sparse Gaussian sketching faster than Intel MKL, but downstream processing may dominate; state-based RNGs matter when sketching is the bottleneck, a situation not yet observed in high-performance RandNLA.CBRNGs trade longer runtimes for reproducibility across compute environments with different levels of parallelism; future hardware-accelerated sketching could warrant reconsidering fast state-based RNGs for RandBLAS.Source: Shift-register RNGs · pp. 30–31

2.1.2 Portability, reproducibility and exception handling

• RandBLAS is planned as a portable procedural API using elementary structs, while RandLAPACK will expose it through an object-oriented operator interface that minimizes assumptions about sketch distributions.The interface also reduces coupling between the procedural API and higher-level prototyping.Source: §2.1.2 · p. 31
• RandBLAS aims for reproducibility through CBRNGs, while most routines propagate NaNs and Infs; extremely sparse sketches may miss them, and bitwise reproducibility depends partly on BLAS behavior.The text notes that summation order can vary with parallelism and vectorization, and leaves reproducible summation approaches to future work.Source: §2.1.2 · p. 31

Prototypical left-sketching and right-sketching

• Left-sketching requires a wide operator and reduces rows while preserving columns; right-sketching requires a tall operator and reduces columns while preserving rows, and the two views are adjoint-equivalent.Left sketches can compress a matrix range and support right-singular-vector estimation, whereas right sketches can sample its range and support left-singular-vector estimation.Source: §2.2 · p. 32, Prototypical left-sketching and… · p. 32, Equivalence of left-sketching and… · p. 33

The embedding and sampling regimes

• Embedding-regime sketches are larger than the data, while sampling-regime sketches are far smaller; the former is common in least squares and the latter underlies low-rank approximation, without fixed left/right assignments.Sketch quality should be assessed through how an operator’s distribution induces quality measures for a given application and embedding dimension.Source: The embedding and sampling regimes · p. 33, §2.2.2 · pp. 33–34

Subspace embeddings

• A subspace embedding preserves a subspace’s geometry by approximating the Gram matrix, but requires sketch dimension at least the subspace dimension and substantial reduction only when the ambient dimension is much larger.Standard distortion is scale-sensitive; effective distortion minimizes distortion over positive rescalings of the sketch and is more relevant to least-squares and optimization algorithms.Source: Subspace embeddings · p. 34, Effective distortion · p. 34

Oblivious subspace embeddings

• OSE theory bounds sketch dimension for given distortion, subspace dimension, and failure probability, but bounds for sophisticated sparse sketches can be pessimistic; effective distortion is more relevant for least squares and optimization.The authors note that practical sketch-distribution selection often works with simple rules of thumb rather than finding the best distribution.Source: Oblivious subspace embeddings · p. 35

Johnson–Lindenstrauss embeddings

• The JL lemma permits data-oblivious dimension reduction with dimension logarithmic in the point-set size; sketch distributions commonly use mean-zero entries, and covariance scale often matters less when algorithms depend on relative norms.For partial sketching, scale can matter: in lifted ridge regression, sampling only the data block treats its rows differently from the regularization block and induces nonuniform sampling.Source: Johnson–Lindenstrauss embeddings · p. 35, §2.2.3 · p. 36, Scale sensitivity: partial sketching · p. 36

2.3 Dense sketching operators

• The authors recommend that RandBLAS support uniform sampling of row- or column-orthonormal matrices; Gaussian sketches have stronger theory, but the three iid entrywise distributions perform similarly in practice.The text notes that implementations such as use uniform sketches and cites theoretical support from universality principles,, [MT20, §8.8],.Source: §2.3 · pp. 37–38

Sampling iid-dense sketching operators

• For Gaussian sampling, Box-Muller is easy to implement and parallelizes well, whereas the more efficient single-threaded Ziggurat method does not parallelize well; Gaussian sampling likely need not exceed single precision.Explicit dense sketches should be applied with an appropriate BLAS routine, most likely gemm; wrappers that generate and apply the sketch could be an initial standard feature.Source: Sampling iid-dense sketching operators · p. 38, Applying iid-dense sketching operators · p. 38

Sampling and applying Haar operators

• Haar operators are uniform over row-orthonormal matrices for left-sketching and column-orthonormal matrices for right-sketching; Householder construction costs O(dm), versus O(d2m) for Gaussian sampling and QR.The authors recommend Haar operators for a first RandBLAS version, both for sketching and for generating test data.Source: Sampling and applying Haar operators · pp. 38–39

Intended use-cases

• Dense sketches suit sampling-regime low-rank approximation, ridge regression, and some full-rank decompositions, but are less useful for embedding dense matrices because they cost more to apply.They may be useful for embedding sparse or otherwise structured data matrices.Source: Intended use-cases · p. 39

2.4 Sparse sketching operators

• The proposed sparse-operator taxonomy is agnostic to left or right sketching and distinguishes short-axis-sparse from long-axis-sparse operators by which vectors are independent.Short-axis vectors have a fixed small number of nonzeros; long-axis vectors sample nonzero indices with replacement from a prescribed distribution.Source: §2.4 · p. 39

2.4.2 Long-axis-sparse sketching operators

• For wide LASOs, sampling indices with replacement yields at most dk nonzeros, versus mk for a SASO with the same parameters, while sketch quality under uniform sampling depends on matrix leverage scores.Using the matrix’s leverage-score distribution defines a LESS embedding [DLD+21]; uniform sampling gives a LESS-uniform operator [DLP+21].Source: §2.4.2 · p. 41

Traditional SRFTs

• Traditional SRFTs can be applied in O(mn log d) time, but are harder to implement efficiently than SASOs when k ≪ log d; the authors do not consider them essential for a preliminary standard.Block SRFTs can use parallel hardware but change their distribution with the processor count, complicating reproducibility across platforms.Source: Traditional SRFTs · p. 42, Block SRFTs · p. 43

2.6 Multi-sketch and quadratic-sketch routines

• RandBLAS should provide multi-sketch primitives and optimized quadratic sketches, which compute linear sketches of A∗A or AA∗ and could reduce data movement.The proposed multi-sketch operations include AS and A∗AS; AS1 and S2A; and AS1, S2A, and S3AS4. Quadratic sketching is described as composing syrk and gemm, and should support dense operators.Source: §2.6 · p. 44

3.1.1 Minimizing regularized quadratics

• For regularized quadratics, higher-accuracy methods repeatedly multiply by G, while lower-accuracy methods may select columns or compute one product GS; suitability depends on how many eigenvalues of G exceed µ.The low-accuracy methods may be useful for kernel ridge regression, where an inaccurate solution can still help downstream tasks; practical KRR data usually have far fewer than n eigenvalues above µ.Source: §3.1.1 · p. 48, Amenable problem structures · p. 48

Pros and cons of this viewpoint

• The saddle point formulations generalize least squares to include linear terms and make primal and dual problems equivalent for most algorithmic purposes.Their optimality conditions can be expressed as a saddle point system over (x, y).Source: §3.1.2 · p. 48, Pros and cons of this viewpoint · p. 49

3.2 Drivers

• Two sketch-and-precondition driver families target accurate approximations, while two cheaper sketch-and-solve families target rough approximations and use sketch-and-precondition methods as core subroutines.The sketch-and-solve drivers in Section 3.2.4 are described as novel.Source: §3.2 · p. 50

3.2.1 Sketch-and-solve for overdetermined least squares

• Sketch-and-solve estimates least-squares solutions from a sketched problem, with error bounded when the sketch is a subspace embedding; it is most compelling for extremely large, row-accessible data.Error estimates can guide downstream processing or indicate whether a more accurate method is needed; one proposed estimator uses statistical bootstrapping.Source: §3.2.1 · p. 50, Application to tensor decomposition · p. 51, Error estimation · p. 51

3.2.2 Sketch-and-precondition for least squares and saddle point problems

• For very tall A, sketch-and-precondition uses a sketched QR factorization to build a preconditioner that makes AM nearly orthogonal, enabling rapid iterative solution with complexity independent of A’s condition number.The approach is intended for very tall matrices; on a personal laptop, MATLAB or Python implementations can compete with LAPACK when m ≥50n ≥105.Source: §3.2.2 · p. 51

Your attention, please!

• With regularization and rapid spectral decay, Nyström preconditioning can find a saddle-point preconditioner in less than O(n3) time for any m ≥n, relaxing the very-tall-matrix requirement.The section also notes applications to block-projection linear-system methods and random-features ridge regression, where Nyström preconditioning has been reported effective when n ≲m.Source: Your attention, please! · p. 52, Applications · p. 55

Algorithms

• A novel algorithm uses sketching data from standard sketch-and-precondition to transform primal or dual saddle point problems into an equivalent primal problem with c = 0.The conversion enables iterative solvers with numerical properties specific to least squares; the authors state that no such conversion routine had been described in the literature.Source: Algorithms · p. 53
• When m ≥n2 and an accurate solution is unnecessary, sketch-and-solve may be preferable because even one sketch-and-precondition iteration can cost as much as the full sketch-and-solve algorithm.The runtime discussion assumes SRFT sketching and compares sketch-and-precondition with direct QR methods costing O(mn2).Source: Algorithms · p. 55, When is sketch-and-precondition… · p. 55

When is sketch-and-precondition asymptotically faster than QR?

• For SRFT sketch-and-precondition, high-probability runtime bounds use d ∈ Ω(n log n), while practice often uses d = sn for small s > 1.The original theoretical analysis required d ∈ Ω(n^2); later analysis reduced this to Ω(n log n), whereas d = sn is preferred in practice because preconditioned matrices tend to be extremely well-conditioned.Source: When is sketch-and-precondition… · pp. 55–56

3.2.3 Nystr¨om PCG for minimizing regularized quadratics

• Nyström PCG is most effective when the eigenvalue threshold λℓ(G) < µ occurs at ℓ ≪ n, making the preconditioned condition number O(1).A rank-ℓ approximation yields a preconditioner stored and applied in O(ℓn) space and time; an adaptive method can select ℓ from (G, µ). A good preconditioner is possible when λℓ(G)/µ is O(1).Source: §3.2.3 · p. 56, Details on the preconditioner · p. 57

3.2.4 Sketch-and-solve for minimizing regularized quadratics

• For KRR, the authors present two novel sketch-and-solve drivers that use sketch-and-precondition as a core subroutine, a form of nested randomization.The section also describes a one-shot Nyström fallback whose approximation can be computed from a single sketch GS and whose dual saddle-point problem can be approached with sketch-and-precondition; a random-subspace constraint likewise reduces to an overdetermined least-squares problem.Source: §3.2.4 · p. 57, A one-shot fallback on Nystr¨om… · p. 57, Applying a random subspace constraint · p. 58

3.3 Computational routines

• The computational-routines discussion introduces preconditioner-generation families and deterministic iterative methods, while omitting routines for sketch-and-solve drivers because they are simple to implement.It also notes that Nyström approximations and data-aware kernel-matrix sketches are treated elsewhere.Source: §3.3 · p. 58, Routines not detailed here · p. 58

3.3.1 Technical background: optimality conditions for saddle point problems

• Saddle-point systems characterize primal- and dual-optimal solutions under the stated existence conditions, while augmented-data reformulations reduce reliance on products with A∗A.For primal problems, randomized preconditioner generation can compute the needed near-minimum-norm solution in O(mn) extra time; when µ = 0, system inconsistency is equivalent to the associated saddle-point problem having no optimal solution.Source: §3.3.1 · p. 59, Inconsistent saddle point systems · p. 59

Sketch and orthogonalize

• The general tall-matrix preconditioning framework applies to least-squares and saddle-point problems when m ≫ n, and a valid preconditioner makes the spectrum of AM equal that of (SU)†.The rank requirement is nearly universal in practice and is conjectured to hold with exponentially-high probability for suitable SASOs and SRFTs; the authors provide what they describe as the first general proof of the proposition, including its regularized application.Source: §3.3.2 · p. 60, Sketch and orthogonalize · p. 60

QR-based preconditioning in the full-rank case

• For full-rank QR preconditioning with µ > 0, the triangular factor can be computed from [Ro; √µI] using specialized QR algorithms with O(n) extra workspace.A Cholesky alternative based on the regularized Gram matrix is numerically unstable in sketch-and-solve preprocessing because it requires solving the normal equations [Bj¨o96].Source: QR-based preconditioning in the… · p. 61

QR-based preconditioning in the rank-deficient case

• In the rank-deficient case, orthogonalizing preconditioners allow a recovered least-squares solution, but a preconditioner independent of b cannot ensure that it is minimum-norm.When k < n, the least-squares minimizers form a nontrivial affine subspace, so the norm of a recovered minimizer may be far from the smallest possible.Source: QR-based preconditioning in the… · pp. 61–62

SVD-based preconditioners

• SVD-based preconditioning supports minimum-norm solutions in rank-deficient cases, regularized construction from singular values and right singular vectors, and recovery of the regularized problem’s SVD for preprocessing.The method also yields a uniquely optimal, well-posed preconditioned solution.Source: SVD-based preconditioners · pp. 62–63

3.3.3 Preconditioning least squares and saddle point problems: data matrices with fast spectral decay

• Nyström preconditioners are positive definite, condition the system on a subspace containing dominant eigenspaces, and act as identity on its orthogonal complement; effectiveness requires limited action outside that subspace.For the saddle-point-compatible construction, the corresponding condition is ∥A∥L⊥/√µ = O(1); the authors leave construction details to future work.Source: §3.3.3 · pp. 63–64

3.3.4 Deterministic preconditioned iterative solvers

• CG, CGLS, LSQR, Chebyshev semi-iteration, and LSMR serve different system classes and stability or parallelism needs; RandLAPACK should include CG and LSQR, with richer libraries adding CS, LSMR, and blocked solvers.Appropriate termination criteria matter, and preconditioning must account for any regularization term.Source: §3.3.4 · pp. 64–65

Facilitating second-order optimization algorithms

• Second-order optimization can amortize randomized-algorithm costs across changing saddle-point systems; Newton Sketch and subsampled Newton methods resketch each iteration, while traditional sketch-and-precondition methods converge faster for least squares.The reweighting-based approach requires applying W and its adjoint to m-vectors in O(m) time, a condition met in some logistic-regression, linear-programming, and iteratively reweighted least-squares algorithms; the least-squares Newton Sketch specialization is not planned for RandLAPACK.Source: Facilitating second-order optimization… · p. 65, Stochastic Newton and subsampled… · pp. 65–66

Random features preconditioning for KRR

• Random-features preconditioning can support iterative KRR solvers, but empirical results suggest Nyström PCG performs better on shared-memory machines, so RandLAPACK does not plan to support random features currently.Both approaches require matrix-vector access to G.Source: Random features preconditioning for KRR · p. 66

Utilities for iterative refinement

• Iterative refinement can reduce rounding error through higher-precision residuals and correction solves, but its best role in randomized algorithms remains unresolved; RandNLA libraries may need methods beyond LAPACK to exploit low precision and accelerators.Sketch-and-precondition may introduce distinct numerical issues.Source: Utilities for iterative refinement · pp. 66–67

3.5 Existing libraries

• Existing sketch-and-precondition libraries cover least squares but not the considered saddle-point problems or Nyström preconditioning; Ski-LLS supports several sketching and factorization options and reports speed advantages on large overdetermined problems.The supplied passages report Ski-LLS as 10 times faster than Sparse QR and incomplete-Cholesky-preconditioned LSQR for specified sparse cases.Source: §3.5 · pp. 67–68

Low-rank Approximation

• Low-rank approximation represents A by factored matrices of rank k ≪ min{m, n}, with user-controlled rank trading approximation accuracy against compression and interpretation needs.For numerically low-rank matrices, the approximation can be built iteratively until the desired error is reached; the section covers efficient and reliable algorithms for both use cases.Source: Low-rank Approximation · pp. 70–71

Singular value decomposition

• Truncating the SVD yields an optimal rank-k approximation in any unitarily invariant norm; small trailing singular values are often treated as noise for denoising.The truncated SVD also supplies canonical approximations at every lower rank without further computation.Source: Singular value decomposition · p. 73

Hermitian eigendecomposition

• For Hermitian matrices, eigendecomposition is preferable to SVD because rank-k forms use almost half the storage and can exploit structure, especially for positive semidefinite matrices.A Hermitian eigendecomposition can also be converted to an SVD with matching spectral- and Frobenius-norm errors.Source: Hermitian eigendecomposition · pp. 73–74

Connections to principal component analysis

• PCA retains the top covariance eigenvectors, computed directly when covariance is given or via the right singular vectors when centered data is given.Randomized SVD and eigendecomposition can scale PCA to large problems; SciKit-Learn used a randomized default for d ≥ 500 and k < 0.8d from 2016 through at least the time of writing.Source: Connections to principal component… · p. 74

4.1.2 Submatrix-oriented decompositions

• Submatrix-oriented decompositions can reduce storage, particularly for sparse data, and make results more interpretable by retaining actual rows or columns.The overview groups CUR and interpolative decompositions in this class.Source: §4.1.2 · p. 75

CUR decomposition

• CUR links selected columns and rows of the target with a middle factor, preserving their data meaning, but high-accuracy approximation of low-numerical-rank matrices can make all factors ill-conditioned.For exact full-rank CUR of tall or wide matrices, the principal freedom is selecting the spanning rows or columns, respectively.Source: CUR decomposition · pp. 75–76

Interpolative decompositions

• Interpolative decompositions feature a target submatrix and interpolation matrices whose regularity conditions affect approximation guarantees; bounding their entries is shown to help control those guarantees.One-sided IDs use selected columns or rows, while a two-sided ID combines them as Z A[I, J] X.Source: Interpolative decompositions · pp. 76–78

On relationships between submatrix-oriented decompositions

• One-sided ID algorithms serve as computational routines underlying CUR and two-sided ID drivers, which are dual in placing submatrices in the outer versus inner factor.CUR leaves its linking factor unrestricted, whereas two-sided ID uses a target submatrix as the inner factor and interpolation matrices as the outer factors.Source: On relationships between… · p. 78
• A two-sided ID can be better numerically behaved than CUR, while CUR can offer greater storage reduction for sparse matrices.For randomized CSS-based CUR, the projection linking matrix is preferable to submatrix inversion; a column-ID-plus-CSS approach instead uses one pseudoinverse rather than the two needed for Uproj.Source: On relationships between… · p. 78, CUR by falling back on CSS · p. 86, CUR by a combination of column ID and… · p. 86

4.1.3 On accuracy metrics

• Low-rank approximation error is usually measured by spectral- or Frobenius-norm distance to the target, not distance to a nominally optimal approximation.Absolute target distance can be large even for a near-optimal result, such as in PCA, with slow spectral decay, or with highly constrained factors; relative distance to a reference is useful, though references may be unavailable or nonoptimal due to rank or structure constraints.Source: §4.1.3 · p. 79, Distance relative to that of a… · p. 80

Distance to optimal approximations

• Truncated SVD approximations can be non-unique, and near-equal singular values can make their singular vectors highly sensitive to small perturbations.For submatrix-oriented decompositions, defining an optimal ID also requires regularity constraints; finding a rank-k column ID with |Xij| ≤ M for M = 1 is NP-hard even when the approximation has exact rank k.Source: Distance to optimal approximations · p. 79

A flexible method

• The flexible SVD method randomizes to form QB, computes its compact SVD deterministically, and truncates to rank k, with the randomized phase capped at rank k + s.Positive oversampling can complicate interpreting the error tolerance; power iteration can improve trailing singular-vector estimates, or users can set k ← k + s and s ← 0.Source: §4.2 · p. 80, A flexible method · p. 81

Sacrificing accuracy for speed

• ID-based SVD conversion can access A with one matrix-matrix multiplication but is less accurate than typical QB methods; single-pass alternatives use multi-sketching.Single-pass methods may not yield highly accurate singular vectors or values, motivating error estimation such as the described bootstrap-based estimator.Source: Sacrificing accuracy for speed · p. 82

A flexible method for Hermitian indefinite matrices

• For Hermitian indefinite matrices, the QB-backed eigendecomposition scales the input tolerance by 1/2 so that with s = 0, ∥QB − A∥ ≤ ϵ implies ∥ˆA − A∥ ≤ ϵ.ID-based and multi-sketching approaches are also described for cases where moving A through fast memory once is important.Source: A flexible method for Hermitian… · pp. 82–83, Sacrificing accuracy for speed with… · p. 83

Nystr¨om approximations for positive semidefinite matrices

• For positive semidefinite matrices, Nyström approximations tend to be more accurate than comparable general Hermitian eigendecomposition methods using the same sketching operator type.The practical EVD2 algorithm truncates a rank-(k + s) approximation to rank k but has no approximation-error control.Source: Nystr¨om approximations for positive… · p. 84, Algorithms · p. 85

CUR by falling back on CSS

• For CUR built via randomized column subset selection, Uproj gives no greater error than submatrix inversion and is better conditioned to compute.The text therefore generally prefers Uproj as the linking matrix in this approach.Source: CUR by falling back on CSS · p. 86

Two-sided ID via one-sided ID

• Two-sided IDs reduce to a randomized column ID followed by a deterministic row ID, and can be converted to CUR using index sets and either a linking matrix or an interpolation matrix.The deterministic row ID is inexpensive under the stated assumption |J| ≪ min{m, n}.Source: Two-sided ID via one-sided ID · pp. 86–87

4.3.1 Power iteration

• Power iteration makes sketch ranges better aligned with dominant singular-vector subspaces, and practical implementations need stabilization between applications of A and A∗.The general formulation allows an arbitrary number of passes over A, including zero passes.Source: §4.3.1 · p. 88

4.3.2 Orthogonal projections: QB and rangefinders

• A QB decomposition approximates A by QQ∗A, while a rangefinder computes Q alone; full-rank QB outputs can be unstructured, though QB remains fundamental for low-rank approximation.In a QB decomposition, Q is column-orthonormal and B = Q∗A.Source: §4.3.2 · p. 88

Iterative QB algorithms

• Iterative QB methods add blocks to Q and B until an error criterion is met; YGL18 methods control Frobenius error or reduce data movement, with the latter risking wasted work at excessive rank limits.The error-controlled method avoids modifying A and updates the Frobenius error efficiently; the data-movement-saving variant can terminate early and, without power iteration, use a single pass.Source: Iterative QB algorithms · p. 89

Stopping criteria for iterative QB algorithms

• When Frobenius error is hard to compute for an abstract operator, randomized norm estimation can support QB stopping criteria; high-quality approximations often instead target spectral-norm error.The Frobenius norm is usually cheap for explicitly stored dense or sparse matrices.Source: Stopping criteria for iterative QB… · p. 89

Approximate single-pass QB via Type 2 multi-sketching

• A Type 2 multi-sketch enables approximate single-pass QB by setting Q = orth(Y1) and B = (S2Q)†Y2, without requiring B = Q∗A.The method relies on Y2 ≈ S2QQ∗A when QQ∗A approximates A.Source: Approximate single-pass QB via Type 2… · p. 90

More details on column-pivoted decomposition algorithms

• Strong RRQR [GE96] gives the best theoretical guarantees for low-rank approximation, but truncating GEQP3 is faster and more common; LUPP is faster yet can yield similar quality for randomized CSS and column ID.The reported similarity applies to some randomized algorithms seeking nearly full-rank decompositions of wide sketches.Source: More details on column-pivoted… · p. 92

4.3.4 One-sided ID and CSS

• Column ID and CSS are nearly equivalent, and the deterministic template derives an ID from a column-pivoted decomposition; the section lists five randomized algorithms for these tasks.A CSS method yields column ID by forming X = (A[:, J])†A, while a column ID method can be adapted to CSS by omitting calculations used only to form X.Source: §4.3.4 · p. 92, Template deterministic algorithms · p. 93, Randomized algorithms · p. 93

Randomized algorithms

• Power-iteration sketching makes a simple, practical low-rank column ID approach toward truncated-SVD error, with ∥A − Ã∥2 ≈ σ_k+1(A).The method computes a nearly full-rank column ID of a sketch, then uses the original data to form the low-rank ID; the authors recommend a data-aware sketch based on power iteration.Source: Randomized algorithms · p. 94

On fixed-accuracy one-sided ID

• Randomized low-rank one-sided ID lacks the specified-accuracy capability of deterministic QRCP-based methods, though an iterative QB approach offers a workaround.An ID of B can certify the target error once the rank-k QB decomposition meets the stated condition; the approach has a moving termination target and requires spectral-norm bounds for implicit operators.Source: On fixed-accuracy one-sided ID · p. 94

4.3.5 Estimating matrix norms

• Random-probe methods estimate spectral and Frobenius norms using matrix-vector products, while the Kong-Valiant estimator uses a smaller data-oblivious sketch for Schatten norms.The Frobenius estimator is unbiased with controllable variance; the cheap spectral bound holds with probability at least 1 − β^-r.Source: §4.3.5 · p. 95

4.3.6 Oblique projections

• Oblique projections form each approximation column by projecting onto the range of AS1 under an inner product induced by S2, connecting the construction to CUR and Nyström methods.Independent sketch operators can cause ill-conditioned computations; the described stable approach uses oversampling for S2 and splits the linking matrix into two factors.Source: §4.3.6 · p. 96

4.4 Other low-rank approximations

• RandLAPACK’s prospective low-rank methods include domain-specific factorizations and selected Cholesky, UTV, and LU approaches, but exclude unpivoted low-rank QR.DMD is expected to receive highest priority among the named specialized methods; low-rank LU plans emphasize Gaussian-elimination-related algorithms, likely suited to fixed-accuracy requirements.Source: §4.4 · p. 97, §4.5 · p. 98

4.5 Existing libraries

• Existing RandNLA libraries offer varied low-rank factorizations, while RandLAPACK plans broader QB and ID/CUR methods, Nyström eigendecomposition, and portability through LAPACK++.RandLAPACK’s near-term priority is algorithms for the factorizations in Section 4.2; longer-term plans include making its C++ implementation callable from Python.Source: §4.5 · pp. 98–99

Further Possibilities for Drivers

• The monograph’s remaining driver topics contribute a Cholesky-QR-based pivoted QR algorithm for ill-conditioned or rank-deficient matrices, introductions to iterative linear-system methods, and current trace-estimation algorithms.It does not cover general randomized eigenvalue problems or randomized matrix-function actions.Source: Further Possibilities for Drivers · p. 100

5.1.1 QR decomposition of tall matrices

• The decomposition section seeks shorter time-to-solution on modern hardware rather than asymptotic speedups; for very tall matrices, unpivoted QR followed by GEQP3 on R yields an equivalent pivoted decomposition and is much faster.The section introduces two types of randomized algorithms for unpivoted QR.Source: §5.1 · p. 101, §5.1.1 · p. 101

Orthogonality in the standard inner product

• For very tall matrices, Cholesky QR can substantially outperform Householder QR and TSQR, while random preconditioning can stabilize it for ill-conditioned inputs [FGL21].The cited analysis of RCholeskyQR2 is in.Source: Orthogonality in the standard inner… · pp. 101–102

Orthogonality in a sketched inner product

• RGS produces sketch-orthogonal columns and can make Q nearly orthonormal in the standard inner product; it is faster than classical Gram–Schmidt and as stable as modified Gram–Schmidt.For A = QR, the sketch-orthogonality condition is (SQ)∗(SQ) = In.Source: Orthogonality in a sketched inner… · p. 102

5.1.2 QR decomposition with column pivoting

• QRCP selects pivots so R’s diagonal approximates A’s singular values and Q’s columns approximate its left singular vectors, supporting low-rank approximation and ill-conditioned least squares.The widely used column-norm pivoting method has finite-precision failure cases that can be resolved by restructuring norm calculations.Source: §5.1.2 · pp. 102–103

An established randomized algorithm for general matrices

• HQRRP can run much faster than standard QRCP by processing matrices in column blocks, allowing the overwhelming majority of operations to use BLAS 3 rather than about half BLAS 2.Typical user-provided parameters are b = 64 and s = 10; HQRRP begins by forming a thin (b + s) × n sketch Y = SA.Source: An established randomized algorithm… · p. 103

A novel randomized algorithm for very tall matrices

• A QRCP-based randomized preconditioner removes the full-rank limitation of preconditioned Cholesky QR [FGL21], and Algorithm 7 returns A[:, J] = QR when rank(SA) = rank(A).The singular values of Apre are reciprocals of those of SU; applying Cholesky QR to blocks in HQRRP would require restoring Q’s Householder representation.Source: A novel randomized algorithm for very… · pp. 104–105
• Cholesky QR may lose orthogonality or fail when cond(Apre) ≥ ϵ^-1/2; practical implementations also need a tolerance-based rank estimate for Rsk.In exact arithmetic, Apre’s condition number depends on neither A nor Ask, so a suitable sketch distribution should yield an accurate decomposition with extremely high probability.Source: A novel randomized algorithm for very… · p. 104

5.1.3 UTV, URV, and QLP decompositions

• UTV, URV, and QLP decompose A as UTV with column-orthogonal U and V and square triangular T, encompassing the SVD and QRCP forms; Stewart studied these decompositions.Source: §5.1.3 · p. 105

Deterministic algorithms

• Stewart’s deterministic UTV applies QRCP twice; the second call gives T a decreasing diagonal, which can track A’s singular values better than R1’s diagonal and permits early stopping [Ste99, §3; Ste99, §5].The complete orthogonal decomposition instead uses one QRCP followed by unpivoted QR, chiefly to apply A† when A is rank-deficient.Source: Deterministic algorithms · pp. 105–106

Randomized algorithms

• Randomized UTV methods include PowerURV, which improves SVD approximations over Stewart’s UTV without pivoted QR [GM18, §3], and a sketching-based incremental full-rank method [MQH19, Figure 4].Using HQRRP in Stewart’s QRCP calls sacrifices guaranteed diagonal decrease across block boundaries; a standard second QRCP can restore it, with speedups expected for matrices of a few thousand dimensions.Source: Randomized algorithms · p. 106

Further Possibilities for Drivers / 5.2 Solving unstructured linear systems

• The discussion turns to solving unstructured linear systems, where randomized solver strategies either reduce safeguards’ cost for fast direct methods or restructure general-purpose iterative methods.Source: Further Possibilities for Drivers · p. 107, §5.2 · p. 107

5.2.1 Direct methods

• LU and LDL factorizations can solve systems in O(n^2) operations and suit parallel processing, but some matrices make them nonexistent or unstable in finite precision, motivating reliability safeguards.Source: §5.2.1 · p. 107

Stability through randomized pivoting

• Randomized partial-pivoted LU achieves GEPP-comparable efficiency with GECP-like element-growth bounds with high probability; randomized pivoted LDL is as stable as GECP and only slightly slower than Bunch–Kaufman variants [FXG18].Source: Stability through randomized pivoting · pp. 107–108

Stability through randomized rotations

• Random orthogonal recursive butterfly transformations precondition nonsingular matrices so unpivoted LU succeeds with high probability, and the approach can handle sparse matrices without catastrophic fill-in.Their FFT-like application and unpivoted decompositions’ data locality have drawn HPC interest; the same preconditioning idea also applies to LDL.Source: Stability through randomized rotations · pp. 108–109

Randomized GMRES: Arnoldi decompositions in a sketch-orthogonal basis

• Sketched Arnoldi builds an S-orthogonal basis and Arnoldi-like decomposition, yielding residual norm at most (1 + δ) times GMRES’s when S has effective distortion δ on Kp+1.Its asymptotic complexity matches classic Arnoldi, but it can use up to half the flops, better exploit BLAS 2, and reduce synchronization points while retaining classic Arnoldi’s reliability.Source: Randomized GMRES: Arnoldi… · pp. 109–110

Randomized GMRES: handling general non-orthogonal bases

• For general non-orthogonal Krylov bases, sketch-and-solve can undercut sketched Arnoldi’s O(np²) cost, but lax basis construction can make AVp ill-conditioned as p grows.In exact arithmetic, factoring SAVp gives the same solutions as sketched Arnoldi; solving the sketched problem to high accuracy via sketch-and-precondition can also produce a solution close to GMRES.Source: Randomized GMRES: handling general… · p. 110

Nested randomization in block-projection and block-descent methods

• Sketch-and-project and block-projection or block-descent solvers have convergence rates that worsen on larger problems, so they may suit cases where an m × m matrix cannot fit in memory [GR15] [PJM22].Their subproblems can be structured to use Algorithm 2, and the general analysis framework also covers methods applicable to underdetermined problems.Source: Nested randomization in… · p. 111

5.3.1 Trace estimation by sampling; Compress and trace

• Girard–Hutchinson trace estimation uses Ω(1/ε²) samples for ε error at constant failure probability, while compress-and-trace can improve relative error when spectral decay is sufficiently fast and a power-iteration rangefinder is well aligned [Epp23; SAI17].For traces of matrix functions such as log(B + I), low-rank compression also yields bounds under suitable spectral-decay assumptions, with some techniques extending to operator-monotone functions such as the matrix square root.Source: §5.3.1 · p. 112, Compress and trace · p. 112

Split, trace, and approximate

• Hutch++ estimates the trace of a positive semidefinite matrix to relative error ε with O(1/ε) matrix-vector products, improving on Girard–Hutchinson’s O(1/ε²) requirement [MMM+21, Theorems 4.1 and 4.2].Its sample complexity cannot be improved for a large class of algorithms. An adaptive extension handles symmetric indefinite matrices but gives additive error guarantees, since relative guarantees are essentially impossible when tr(A) = 0 [PCK21].Source: Split, trace, and approximate · p. 113

Leveraging the exchangability principle

• XTrace symmetrizes Hutch++ by averaging m runs constructed from shared samples, and a careful implementation achieves Hutch++’s asymptotic complexity [ETW23].XTrace also provides adaptive stopping and variance estimation; the motivation is that Hutch++ uses separate randomness stages and does not satisfy the exchangability principle.Source: Leveraging the exchangability principle · pp. 113–114

5.3.3 Estimating the trace of f(B) via integral quadrature

• Stochastic Lanczos quadrature estimates tr(f(B)) for sufficiently smooth f and Hermitian B; guarantees usually require stronger conditions such as analyticity over B’s spectral interval or positive semidefiniteness [CH22] [PK22] [CTU22].Gaussian quadrature is exact for polynomials through degree 2s−1 using s samples, and its nodes and weights can be recovered from the eigendecomposition of an s × s Jacobi matrix [GM10, Theorem 6.2].Source: §5.3.3 · p. 114, Technical background · p. 115

Stochastic Lanczos quadrature : approximating Girard–Hutchinson

• SLQ estimates traces by approximating Girard–Hutchinson sample integrals with Gaussian quadrature, using Lanczos polynomials and a Jacobi matrix to obtain quadrature nodes and weights [BFG96; UCS17; GM10].The Lanczos basis spans a Krylov subspace; standard tridiagonal eigensolvers then provide the nodes and weights.Source: Stochastic Lanczos quadrature :… · pp. 115–116
• SLQ sample accuracy depends on quadrature nodes, while s Lanczos steps require s − 1 matrix-vector products and trade full reorthogonalization’s O(ns) storage and O(ns^2) arithmetic for O(n) storage and O(ns) arithmetic without reorthogonalization.The stated costs exclude the matrix-vector products with B.Source: Stochastic Lanczos quadrature :… · p. 116

Beyond stochastic Lanczos quadrature

• For positive definite A, SLQ inherits the Girard–Hutchinson estimator’s O(1/ϵ^2) sample requirement, while Hutch++-style variance reduction can reduce queries to O(1/ϵ); SLQ nodes and weights also transfer across functions for spectral-density estimation [CH22; PK22; Lin16; CTU22].The same quadrature nodes and weights used for tr(f(B)) can estimate tr(g(B)); the associated spectral distribution encodes B’s spectrum.Source: Beyond stochastic Lanczos quadrature · p. 117

Advanced Sketching: Leverage Score Sampling

• Leverage scores measure subspace alignment with coordinates and induce data-aware sampling distributions that are relatively insensitive to matrix numerical properties, unlike uniform sampling, which can perform poorly [DMM+12] [BMD09].Computing scores can be expensive, however, and their approximation costs may offset the efficiency of row or column selection; they are most useful when selection is required.Source: Advanced Sketching: Leverage Score… · pp. 118–119

Probabilistic guarantees of sketching via row sampling

• When q approaches the leverage-score distribution, row-sampling sketches are more likely to be good; uniform sampling can require more samples and be much worse, while q_j ≥ βp_j(A) permits d = d_lev/β for the same guarantees.Exact leverage-score sampling has r = 1 in the stated bound; approximate scores are one practical route to a data-aware distribution.Source: Probabilistic guarantees of sketching… · pp. 120–121

Preconditioned leverage score sampling, hidden in plain sight

• Rotations act as sampling preconditioners: an SRFT rotates the input before coordinate sampling, and embedding dimension O(n log n) suffices for constant distortion with high probability.The passage also describes Haar sketching as an orthogonal transform followed by coordinate sampling.Source: Preconditioned leverage score… · p. 121

6.1.2 Subspace leverage scores

• Subspace leverage scores address limits of standard scores for low-rank approximation, but their rank-k spaces may be discontinuous or nonunique; a spectral gap or approximation to a nearby rank-k space can mitigate these issues [DMM+12, §5; DMM08].They can be nonuniform for both range(A) and range(A∗) when k < rank(A), and row-sampling projections may require d ≫ k.Source: §6.1.2 · pp. 122–123

6.1.3 Ridge leverage scores; 6.2 Approximation schemes; 6.2.1 Standard leverage scores

• Ridge leverage scores [AM15] support data-aware column sampling for Nyström approximations of regularized PSD matrices; the effective rank tr(K(K + λI)−1) determines approximation requirements, and low-rank approximations make inverse application cheaper.Source: §6.1.3 · p. 123

6.2.1 Standard leverage scores

• Approximating both the pseudoinverse and the product AA† reduces leverage-score estimation to O(d1n^2 + d2mn), with sketch dimensions chosen for useful accuracy guarantees.[DMM+12] gives the dimension requirements; and describe faster alternatives for the initial sketch.Source: §6.2.1 · pp. 124–125

6.2.2 Subspace leverage scores

• A QB decomposition yields approximate rank-k leverage scores by taking row norms of Q times the leading k left singular vectors of B; the original method targets Frobenius-near-optimal rank-k approximation.The user specifies oversampling s ∈ O(k), and the decomposition has rank k + s.Source: §6.2.2 · p. 125

6.2. Approximation schemes

• [DMM+12] estimates rank-k leverage scores for a spectrally near-optimal rank-k approximation using power iteration, then estimates leverage scores of the resulting sketch.A prescribed q gives a factor guarantee for the estimated scores relative to those of a rank-k matrix satisfying the stated approximation condition.Source: §6.2 · p. 126

6.2.3 Ridge leverage scores

• Rank-n ridge-leverage-score sampling can produce, with high probability, a rank-n approximation with spectral-norm error at most λ in O(mn^2) time using n column samples.Here n ∈ O(d log d), where d = tr(K(K + λI)^−1) is the effective rank of K. Iterative methods can improve the accuracy of ridge-leverage-score estimates, but are not summarized in detail.Source: §6.2.3 · p. 127

6.3.1 Leverage score sparsified embeddings

• When k ≈ n = rank(A), LESS sketches are nearly indistinguishable from dense sub-Gaussian operators with high probability, while reducing application time from O(dmn) to O(dn^2).Approximate leverage scores suffice, but estimating them typically dominates LESS embedding cost; using them in the sparsification pattern is essential for the theoretical guarantee.Source: §6.3.1 · p. 127

6.3.2 Determinantal point processes

• DPP sampling jointly selects diverse indices and yields Nyström approximations with nuclear-norm error at most (1 + ϵ) times the best rank-k error when the sample size is at least k.Efficient DPP sampling remains challenging; ridge-leverage-score-based approaches incur overhead because some intermediate samples must be discarded.Source: §6.3.2 · p. 128

6.3.3 Further variations on leverage scores

• Leverage-score sampling has extensions including Lewis weights for preserving p-norms and probabilities for estimating variability in sketch-and-solve least-squares solutions.The passage also points readers to practical guidance and explicit-constant results in [LK20, §4 and Appendix A], and to the pylspack library [SG22].Source: §6.3 · p. 127, §6.3.3 · p. 129

Advanced Sketching: Tensor Product Structures

• Kronecker and Khatri–Rao products can grow exponentially with the number of factors, motivating efficient sketches for implicit tensor-structured matrices; row-structured operators apply efficiently to these products.The discussion treats tensor routines as RandLAPACK utility functions rather than general tensor-computation algorithms, and covers several sketch families for structured matrices.Source: Advanced Sketching: Tensor Product… · p. 130, §7.1 · p. 131, A note on scope · p. 131

Khatri–Rao products of elementary sketching operators

• Sun et al. give guarantees for elementary sketches at L = 2, while vectorized CP and tensor-train sketches enable efficient application to Kronecker and Khatri–Rao products.Iwen et al. [INR+20] propose a two-stage operator that first applies a row-structured Kronecker sketch and then maps the intermediate representation to dimension d.Source: Khatri–Rao products of elementary… · p. 132, Row-wise vectorized tensors · p. 133, Two-stage operators · p. 133

7.2.2 The Kronecker SRFT

• Kronecker SRFTs apply efficiently without explicitly forming product vectors, while TensorSketch accelerates Kronecker-product sketching at the cost of a larger embedding dimension than CountSketch.For Kronecker SRFTs, only the sampled entries of the factor transforms need to be computed, using subsampled FFT methods.Source: §7.2.2 · pp. 133–134, §7.2.3 · p. 134

7.2.4 Recursive sketching

• Recursive sketching reduces the required embedding dimension from exponential to polynomial in L, while a graph-based generalization achieves linear scaling and supports tensor-network vectors.The cited recursive results consider equal-length vectors; the text notes they should extend to unequal lengths.Source: §7.2.4 · pp. 134–135

7.2.5 Leverage score sampling for implicit matrices with tensor product structures

• For implicit tensor-structured least squares, leverage-score sampling can avoid accessing every entry of the unstructured right-hand side, and Kronecker-product leverage scores permit exact sampling without forming the matrix.This is useful when accessing all right-hand-side entries is prohibitively expensive; sampled rows can be computed and rescaled directly.Source: §7.2.5 · pp. 135–136, Kronecker product structure · p. 136

Khatri–Rao product structure

• Malik et al. [MBM22] give a method for sampling Khatri–Rao product matrices by their exact leverage scores without forming the matrix, alongside earlier coarse-score approaches.The earlier methods use Kronecker-product leverage scores as approximations and provide sketch-and-solve guarantees; Malik improves sampling and computational complexities, and [MBM22] generalizes the method.Source: Khatri–Rao product structure · p. 137

7.3 Partial updates to Kronecker product sketches

• Structured sketches can reuse smaller component sketches across multiple calls, a functionality motivated by tensor-decomposition algorithms.The discussion aims to identify functionality useful for structured sketches in the tensor community.Source: §7.3 · pp. 137–138

7.3.1 Background on the CP decomposition

• CP decomposition represents a tensor as a sum of rank-1 tensors, and alternating least squares updates one factor matrix at a time as a practical heuristic.The general tensor-decomposition problem is computationally intractable; ALS cycles through factor matrices until a termination criterion is met.Source: §7.3.1 · p. 138

Formulating and solving the least squares problem

• In CP-ALS, computing X(n)A̸=n is the most expensive part of the update, motivating row-sampling sketching for least-squares subproblems.The Gram matrix can be computed efficiently, so this tensor setting makes the matrix product—not the normal-equations Gram computation—the bottleneck.Source: Formulating and solving the least… · p. 139

7.3.2 Sketching for the CP decomposition

• For CP updates, Battaglino et al. [BBK18] reuse fixed Fourier factors and redraw the sampling operator, avoiding repeated expensive Kronecker SRFT applications.The fixed factors are F1, . . . , FL; the sampling operator R is drawn anew for each least-squares problem.Source: §7.3.2 · p. 140
• Larsen and Kolda [LK20] use leverage-score sampling for CP updates, reusing each sampling distribution across L −1 least-squares problems before recomputing it.The sampling probabilities are formed from the leverage-score distributions of the factor matrices.Source: §7.3.2 · p. 140

7.3.4 Sketching for the Tucker decomposition

• Reusing two CountSketch sets for Tucker updates improves runtime and can enable a single pass over the tensor.The two sets allow different embedding dimensions for the two least-squares problems, with d2 > d1; relevant sketches of the unfoldings can be computed at the start.Source: §7.3.4 · p. 141

A.1 Subspace embeddings and effective distortion

• Effective distortion addresses the scale dependence of ordinary subspace-embedding distortion by optimizing over a positive rescaling of the sketch.The definition uses an infimum over t > 0 because D(S; L) = D(−S; L).Source: §A.1 · p. 142, Details on Basic Sketching · p. 143

A.1.1 Effective distortion of Gaussian operators

• For any fixed n-dimensional subspace, Gaussian sketching results justify aggressively small embedding dimensions in randomized least-squares algorithms.The cited concentration results are attributed to Silverstein and Geman; Meng, Saunders, and Mahoney reach the same conclusion for LSRN [MSM14, Theorem 4.4].Source: §A.1.1 · p. 144

Constructing SASOs column-wise

• Wide SASOs are cheap to construct and store, can be generated columnwise in parallel, and favor CSC storage; CSR may be constructed columnwise with extra data for conversion.The authors recommend considering both CSR and CSC, with CSC the more natural choice for columnwise construction.Source: Constructing SASOs column-wise · p. 144, Remarks on storage formats · p. 145

Applying a wide SASO

• Wide SASO application has mediocre single-thread performance but excellent many-core scaling for one format pairing, while CSR sketches with column-major data may require conversion to CSC.The reported scaling observation is preliminary.Source: Applying a wide SASO · p. 146

Selecting parameters for SASOs

• For SASOs, sketch quality rises with sparsity k before plateauing; in one tested embedding, k above eight brought no real benefit and restricted condition numbers were O(1).The example embeds the range of a 100,000 × 2,000 matrix into ambient dimension d = 6,000.Source: Selecting parameters for SASOs · p. 147

A.3 Theory for sketching by row selection

• Theorem A.3.1 gives eigenvalue-tail bounds for sums of independent positive semidefinite random matrices, parameterized by summand and mean-matrix eigenvalues.The setup assumes an almost-sure upper bound on each summand’s largest eigenvalue and defines the mean matrix’s largest and smallest eigenvalues.Source: §A.3 · p. 147
• The row-selection proposition bounds the distortion of sketched norms for every vector in the range of a full-rank matrix.Its proof normalizes sampled row outer products by the Gram matrix and uses the relation ℓ_j(A) = n p_j(A).Source: §A.3 · pp. 147–148

B.2.2 Basic sensitivity analysis for least squares problems

• Appendix B.2 provides classical numerical-linear-algebra background on accuracy for overdetermined and underdetermined least-squares solutions.Source: §B.2.2 · p. 150

B.1 Quality of preconditioners

• When rank(Ask) = rank(A), the singular values of AM are reciprocals of those of SU, given the stated orthonormal-basis conditions.The proposition is a linear algebraic result; randomness enters randomized algorithms only through construction of the sketch.Source: §B.1 · p. 151
• For ridge-regularized problems, the preconditioner result applies to augmented matrices, with S replaced by the augmented sketching operator that maps  to Âsk.Source: §B.1 · p. 151
• If AskM has orthonormal columns spanning range(Ask) and Ask = AskMB, then the lemma establishes M = B†.This lemma supplies an explicit expression for the orthogonalizing map under its stated full-rank assumptions.Source: §B.1 · p. 151

B.1.1 Effective distortion in sketch-and-precondition

• In saddle-point sketch-and-precondition for tall A, the preconditioned matrix’s condition number equals the restricted condition number of S on range(A).This is the stated algorithmic consequence when the sketch-and-precondition methodology is used.Source: §B.1.1 · p. 153
• The authors connect effective distortion to the convergence rate of LSQR or CGLS for least squares with data matrix AM.They identify the right-hand side of their effective-distortion expression as that convergence rate.Source: §B.1.1 · p. 153

• Finite-precision arithmetic makes computed solutions differ from exact solutions, even when reliable algorithms are used.This is an unavoidable consequence of finite-precision computation.Source: §B.2 · p. 153
• Large-scale computations may trade computational complexity for accuracy, yielding approximate solutions even in exact arithmetic.Source: §B.2 · p. 153
• The appendix summarizes error-analysis literature because some randomized algorithms claim to solve problems to any desired accuracy, with basic least squares as its focus.Appendices B.2.1–B.2.4 cover material typically taught in introductory numerical analysis; B.2.5 mentions additional topics.Source: §B.2 · p. 153

B.2.1 Concepts: forward and backward error

• Forward error measures an approximate solution’s distance from the exact solution, but is often unavailable and algorithm-specific analyses may not transfer.These limitations motivate backward error and sensitivity analysis.Source: §B.2.1 · p. 153
• Backward error asks how much the data must change for the computed solution to be exact, while sensitivity asks how small data changes affect the exact solution.Source: §B.2.1 · p. 154
• Combining backward-error bounds with sensitivity analysis estimates forward error for perturbations from rounding, data uncertainty, or deliberate approximation.Sensitivity analysis is agnostic to the source of the perturbation.Source: §B.2.1 · p. 154
• Much algorithm-specific error analysis can focus on backward error, and many problems permit cheap runtime upper bounds on a solution’s backward error.Source: §B.2.1 · p. 154
• Together, backward error and sensitivity analysis yield a priori behavior guarantees and a posteriori solution-quality guarantees, but runtime forward-error bounds may require unavailable data such as extreme singular values.Source: §B.2.1 · p. 154

B.2.2 Basic sensitivity analysis for unregularized least squares problems

• For overdetermined least squares, the forward-error bound for δx depends on the square of A’s condition number, a fundamental obstacle to high-accuracy solutions.Theorem B.2.2 assumes b is neither in the range nor the kernel of A∗ and perturbations satisfy ∥δA∥2 < σn(A); the theorem statements require full-rank A.Source: §B.2.2 · pp. 154–155
• Theorem B.2.3 gives a sensitivity result for A(δx), indicating more hope for high forward-error accuracy in dual than primal saddle point problems.Under full rank, the dual solution is y = b − Ax, and δy = −A(δx) transfers the perturbed primal solution to the dual.Source: §B.2.2 · p. 155
• For underdetermined least squares, the adapted [GV13] Theorem 5.6.1 provides an even more favorable sensitivity analysis result.Theorem B.2.4 considers y = (A∗)†c for nonzero c and perturbations to A and c.Source: §B.2.2 · pp. 155–156

B.2.3 Sharper sensitivity analysis for overdetermined problems

• For overdetermined least squares, sharper bounds normalize relative error by the perturbed solution; when the matrix ranges match, solution sensitivity depends linearly on κ(Ao).The discussion also distinguishes sensitivity to matrix and right-hand-side perturbations and points to componentwise analysis using |δAij| ≤ α|Aij| [Hig02, §20.1].Source: §B.2.3 · p. 156

B.2.4 Simple constructions to bound backward error

• Backward-error constructions make approximate solutions exact for perturbed data, providing normwise bounds useful as iterative-solver stopping criteria.For inconsistent problems, the construction uses a rank-1 perturbation whose spectral norm is readily computed; with LSQR in exact arithmetic, δbo is zero [Ste77, Theorem 3.2; PS82, §6.2].Source: §B.2.4 · p. 156, Inconsistent overdetermined problems · p. 157

Consistent overdetermined problems: a word of warning

• For nearly consistent problems, the inconsistent-problem perturbations are unsuitable, and setting ϵA = ϵb can produce unreliable behavior.The LSQR construction uses tolerances ϵA and ϵb tied to allowable data perturbations; the authors recommend ϵA = 0, or ϵA = (ϵb)^2 if that is too extreme [PS82, §6.1].Source: Consistent overdetermined problems: a… · pp. 157–158

B.2.5 More advanced concepts

• Backward error plus sensitivity can overestimate forward error, while runtime sensitivity bounds require accurate extreme-singular-value estimates.Krylov methods can estimate forward error or condition numbers; LSQR estimates ∥L∥F∥L†∥F, and sketch singular values can approximate reciprocals of those of Ao, though cheap smallest-singular-value estimates may be poor.Source: §B.2.5 · p. 158

B.3 Ill-posed saddle point problems

• When A is rank-deficient and µ = 0, saddle point problems may be infeasible or unbounded, so the appendix defines canonical solutions by limiting analysis.The canonical dual solution is y(0) = (A*)†c + (I − AA†)b; the primal solution is defined by a limit even when c is not orthogonal to ker(A), in which case ∥x(µ)∥ diverges as µ tends to zero.Source: §B.3 · pp. 159–160

B.4.2 Efficient sketch-and-solve for regularized quadratics

• The appendix presents a sketch-and-precondition procedure for regularized quadratics using a Nyström approximation, including preconditioner generation when A is represented implicitly.The implementation avoids applying R−1 in the iterative solver.Source: §B.4 · p. 160, §B.4.2 · pp. 162–163

B.4.1 A primer on kernel ridge regression

• Kernel ridge regression is formulated as a regularized linear algebra problem, and ordinary ridge regression is a special case with kernel matrix X*X.In the Hilbert-space formulation, the solution lies in the span of the kernel functions {kxi}.Source: §B.4.1 · p. 160, A more technical description · p. 161, Why is ridge regression a special case… · p. 162

C.1.1 Approximation quality in low-rank ID

• For a rank-k approximation Atilde with spectral-norm error at most ϵ, replacing it by a column ID formed from its selected columns yields a low-rank ID with a corresponding error bound.The supplied passages establish that the selected-column construction preserves a controlled spectral-norm error, but do not provide the bound’s full expression.Source: §C.1.1 · pp. 164–165
• Bound (4.17) can be sharpened: M need only bound entries of X outside its identity submatrix.Source: §C.1.1 · p. 165

C.1.2 Truncation in column-pivoted matrix decompositions

• For diagonal basis changes, g_k(D) equals the standard truncation (F[:, :k])(T[:k, :])P∗.Source: §C.1.2 · p. 165
• The identity-basis truncation error is ∥F2T2∥, and the identity is optimal exactly when range(F2) and range(F1) are orthogonal.Source: §C.1.2 · p. 165

C.2 Computational routine interfaces and implementations

• The routines use core interfaces for orthogonalization, sketch generation, and stabilization, with stabilization requiring only improved conditioning while preserving range.An economic QR factor can implement Orth; Stabilizer may use a cheaper pivoted LU factor. SketchOpGen returns an oblivious sketching operator.Source: §C.2 · p. 167

C.2.1 Power iteration for data-aware sketching

• Power-iteration TallSketchOpGen produces a sketch aligned with A’s leading right-singular subspace; p = 0 gives an oblivious sketch, while p > 0 is recommended for slowly decaying singular values.The parameter q controls how many multiplications by A or A∗ occur before stabilization; examples for p include 2 or 3.Source: §C.2.1 · pp. 168–169

C.2.2 RangeFinders and QB decompositions

• RF1 orthogonalizes a sketch to approximate A’s leading left-singular subspace, while QB decomposition aims for ∥A − QB∥ ≤ ϵ with rank(QB) ≤ min{k, rank(A)}.QB1 delegates accuracy control to its RangeFinder; QB2 is fully adaptive with fine-grained Frobenius-error control, and QB3 is pass-efficient and partially adaptive, using tolerance for early stopping and disallowing k = min{m, n}.Source: §C.2.2 · pp. 170–171

C.2.3 ID and subset selection

• The Python prototype implements randomized one-sided ID and subset-selection algorithms that apply QRCP to sketches; the appendix presents OSID1 and ROCS1, with two other randomized ID algorithms omitted.The omitted prototype algorithms are for one-sided and two-sided ID.Source: §C.2.3 · pp. 173–175

Correctness of Preconditioned Cholesky QRCP

• The appendix begins proving Proposition 5.1.2 by restating the algorithm under assumptions including A ∈ R^{m×n}, S ∈ R^{d×m}, n ≤ d ≪ m, and a column-orthonormal Qsk.The setup also specifies that J is a permutation of the n indices and k = rank(Ask).Source: Correctness of Preconditioned Cholesky… · p. 176
• When rank(SA)=rank(A), the proof establishes A[:,J]=QR, with Q an orthonormal basis for range(A) when rank(A)=k.The argument uses the assumption that sketching does not drop rank to show that dependencies among sketched columns also hold among the original columns.Source: Correctness of Preconditioned Cholesky… · pp. 176–177

Bootstrap Methods for Error Estimation

• A posteriori error estimates use available problem data to assess solution accuracy and support adaptive computation, such as deciding whether to iterate further or tune parameters.They can help refine an initial solution until the desired accuracy is reached with just enough work.Source: Bootstrap Methods for Error Estimation · p. 178
• Bootstrap error-estimation methods have been developed for only a handful of sketch-and-solve algorithms, leaving broader coverage an open research direction.The appendix records the authors’ considerations and offers non-experts a starting point for surveying this evolving literature.Source: Bootstrap Methods for Error Estimation · p. 178

E.1 Bootstrap methods in a nutshell

• An error estimate seeks a bound on the random, unknown error that holds with probability near 1 − α, using the empirical (1 − α)-quantile of bootstrap samples.The bootstrap generates approximate error samples when independent samples of the actual error cannot be obtained.Source: §E.1 · p. 179
• Bootstrap methods estimate error quantiles using approximate samples to mimic the distribution of an unknown random error, with asymptotic validity and quantitative guarantees in broad settings.The framework treats an exact linear algebra solution as an unknown parameter and a randomized estimate as an estimate of that parameter.Source: §E.1 · pp. 179–180

E.2 Sketch-and-solve least squares

• For sketch-and-solve least squares, the exact solution is the target, the sketched matrix and right-hand side are the dataset, and the sketch-and-solve solution is the estimate.Error can be measured with a preferred metric such as one induced by the ℓ2 or ℓ∞ norms; extensions to iterative randomized least-squares algorithms are discussed in.Source: §E.2 · p. 180
• The least-squares bootstrap resamples d rows with replacement from the sketched matrix and right-hand side, then returns a bootstrap quantile estimate for the solution error.The procedure uses B samples to estimate the (1 − α)-quantile of ∥ˆx−x⋆∥.Source: §E.2 · p. 180
• The least-squares bootstrap iterations are embarrassingly parallel, use only sketched quantities, and have cost independent of the number of rows of A.Source: §E.2 · p. 180

E.3 Sketch-and-solve one-sided SVD

• For sketch-and-solve one-sided SVD, the bootstrap method estimates (1 − α)-quantiles of errors in the top k singular values and right singular vectors [LEM20].The errors are measured uniformly over the top k singular values and right singular vectors.Source: §E.3 · p. 181
• A variant estimates errors for approximate left singular vectors [LEM20], while bootstrapped spectral-norm error bounds for ˆA*ˆA − A*A imply bounds for all sketched singular values and right singular vectors [LEM23].The latter implication uses the Weyl and Davis–Kahan inequalities.Source: §E.3 · p. 181