2.1. Basics

We will work in the real field ( ${\mathbb {R}}$ ) or the complex field ( $\mathbb {C}$ ). The symbol $\mathbb {F}$ refers to either the real or complex field, in cases where the precise choice is unimportant. As usual, scalars are denoted by lowercase italic Roman ( $a, b, c$ ) or Greek ( $\alpha , \beta $ ) letters.

Vectors are elements of $\mathbb {F}^n$ , where n is a natural number. We always denote vectors with lowercase bold Roman ( $\mathbf {\mathsf {a}}, \mathbf {\mathsf {b}}, \mathbf {\mathsf {u}}, \mathbf {\mathsf {v}}$ ) or Greek ( $\mathbf {{\alpha }}, \mathbf {{\beta }}$ ) letters. We write $\mathbf {\mathsf {0}}$ for the zero vector and $\mathbf {\mathsf {1}}$ for the vector of ones. The standard basis vectors are denoted as $\mathbf {{\delta }}_1, \ldots , \mathbf {{\delta }}_n$ . The dimensions of these special vectors are determined by context.

A general matrix is an element of $\mathbb {F}^{m \times n}$ , where $m, n$ are natural numbers. We always denote matrices with uppercase bold Roman ( $\mathbf {\mathsf {A}}, \mathbf {\mathsf {B}}, \mathbf {\mathsf {C}}$ ) or Greek ( $\mathbf {\mathsf{\Delta }}, \mathbf {\mathsf{\Lambda }}$ ) letters. We write $\mathbf {\mathsf {0}}$ for the zero matrix and ${\mathbf {\mathsf {I}}}$ for the identity matrix; their dimensions are determined by a subscript or by context.

The parenthesis notation is used for indexing into vectors and matrices: $(\mathbf {\mathsf {a}})_i$ is the ith coordinate of vector $\mathbf {\mathsf {a}}$ , while $(\mathbf {\mathsf {A}})_{ij}$ is the $(i, j)$ th coordinate of matrix $\mathbf {\mathsf {A}}$ . In some cases it is more convenient to invoke the functional form of indexing. For example, $\mathbf {\mathsf {A}}(i, j)$ also refers to the $(i, j)$ th coordinate of the matrix $\mathbf {\mathsf {A}}$ .

The colon notation is used to specify ranges of coordinates. For example, $(\mathbf {\mathsf {a}})_{1:i}$ and $\mathbf {\mathsf {a}}(1:i)$ refer to the vector comprising the first i coordinates of $\mathbf {\mathsf {a}}$ . The colon by itself refers to the entire range of coordinates. For instance, $(\mathbf {\mathsf {A}})_{i:}$ denotes the ith row of $\mathbf {\mathsf {A}}$ , while $(\mathbf {\mathsf {A}})_{:j}$ denotes the jth column.

The symbol ${}^{*}$ is reserved for the (conjugate) transpose of a matrix of vector. A matrix that satisfies $\mathbf {\mathsf {A}} = \mathbf {\mathsf {A}}^{*}$ is said to be self-adjoint. It is convenient to distinguish the space $\mathbb {H}_n$ of self-adjoint $n \times n$ matrices over the scalar field. We may write $\mathbb {H}_n(\mathbb {F})$ if it is necessary to specify the field.

The operator ${}^{\dagger }$ extracts the Moore–Penrose pseudoinverse of a matrix. More precisely, for $\mathbf {\mathsf {A}} \in \mathbb {F}^{m \times n}$ , the pseudoinverse $\mathbf {\mathsf {A}}^{\dagger } \in \mathbb {F}^{n \times m}$ is the unique matrix that satisfies the following:

1. (1) $\mathbf {\mathsf {AA}}^{\dagger }$ is self-adjoint;

2. (2) $\mathbf {\mathsf {A}}^{\dagger } \mathbf {\mathsf {A}}$ is self-adjoint;

3. (3) $\mathbf {\mathsf {AA}}^{\dagger } \mathbf {\mathsf {A}} = \mathbf {\mathsf {A}}$ ;

4. (4) $\mathbf {\mathsf {A}}^{\dagger } \mathbf {\mathsf {A}} {\mathbf {\mathsf {A}}}^{\dagger } = \mathbf {\mathsf {A}}^{\dagger }$ .

If $\mathbf {\mathsf {A}}$ has full column rank, then $\mathbf {\mathsf {A}}^{\dagger } = ({\mathbf {\mathsf {A}}^{*}} \mathbf {\mathsf {A}})^{-1} {\mathbf {\mathsf {A}}^{*}}$ , where $(\cdot )^{-1}$ denotes the ordinary matrix inverse.

2.2. Eigenvalues and singular values

A positive semidefinite matrix is a self-adjoint matrix with nonnegative eigenvalues. We will generally abbreviate positive semidefinite to PSD. Likewise, a positive definite (PD) matrix is a self-adjoint matrix with positive eigenvalues.

The symbol denotes the semidefinite order on self-adjoint matrices. The relation means that $\mathbf {\mathsf {B}} - \mathbf {\mathsf {A}}$ is PSD.

We write $\lambda _1 \geq \lambda _2 \geq \cdots $ for the eigenvalues of a self-adjoint matrix. We write $\sigma _1 \geq \sigma _2 \geq \cdots $ for the singular values of a general matrix. If the matrix is not clear from context, we may include it in the notation, so that $\sigma _{j}(\mathbf {\mathsf {A}})$ is the jth singular value of $\mathbf {\mathsf {A}}$ .

Let $f \colon {\mathbb {R}} \to {\mathbb {R}}$ be a function on the real line. We can extend f to a spectral function $f \colon \mathbb {H}_n \to \mathbb {H}_n$ on self-adjoint matrices. Indeed, for a matrix $\mathbf {\mathsf {A}} \in \mathbb {H}_n$ with eigenvalue decomposition

$$ \begin{align*} \mathbf{\mathsf{A}} = \sum\nolimits_{i=1}^n \lambda_i \, \mathbf{\mathsf{u}}_i {\mathbf{\mathsf{u}}}_i^{*}, \quad\text{we define } f(\mathbf{\mathsf{A}}) := \sum\nolimits_{i=1}^n f(\lambda_i) \, \mathbf{\mathsf{u}}_i {\mathbf{\mathsf{u}}}_i^{*}. \end{align*} $$

The Pascal notation for definitions ( $:=$ and $=:$ ) is used sparingly, when we need to emphasize that the definition is taking place.

2.3. Inner product geometry

We equip $\mathbb {F}^n$ with the standard inner product and the associated $\ell _2$ norm. For all vectors $\mathbf {\mathsf {a}}, \mathbf {\mathsf {b}} \in \mathbb {F}^{n}$ ,

$$ \begin{align*} {\langle { \mathbf{\mathsf{a}} }, \, { \mathbf{\mathsf{b}} } \rangle} := \mathbf{\mathsf{a}} \cdot \mathbf{\mathsf{b}} := \sum\nolimits_{i=1}^n ({\mathbf{\mathsf{a}}})_i^{*} (\mathbf{\mathsf{b}})_i \quad\text{and}\quad {\Vert \mathbf{\mathsf{a}} \Vert}^2 := {\langle { \mathbf{\mathsf{a}} }, \, { \mathbf{\mathsf{a}} } \rangle}. \end{align*} $$

We write $\mathbb {S}^{n-1}$ for the set of vectors in $\mathbb {F}^n$ with unit $\ell _2$ norm. If needed, we may specify the field: $\mathbb {S}^{n-1}(\mathbb {F})$ .

The trace of a square matrix is the sum of its diagonal entries:

$$ \begin{align*} \operatorname{trace}(\mathbf{\mathsf{A}}) := \sum\nolimits_{i=1}^n (\mathbf{\mathsf{A}})_{ii} \quad\text{for }\mathbf{\mathsf{A}} \in \mathbb{F}^{n \times n}. \end{align*} $$

Nonlinear functions bind before the trace. We equip $\mathbb {F}^{m \times n}$ with the standard trace inner product and the Frobenius norm. For all matrices $\mathbf {\mathsf {A}}, \mathbf {\mathsf {B}} \in \mathbb {F}^{m \times n}$ ,

$$ \begin{align*} {\langle { \mathbf{\mathsf{A}} }, \, { \mathbf{\mathsf{B}} } \rangle} := \operatorname{trace}(\mathbf{\mathsf{A}}^{*} \mathbf{\mathsf{B}}) \quad\text{and}\quad {{\Vert {{ \mathbf{\mathsf{A}} }} \Vert}^2_{\mathrm{F}}} := {\langle { \mathbf{\mathsf{A}} }, \, { \mathbf{\mathsf{A}} } \rangle}. \end{align*} $$

For vectors, these definitions coincide with the ones in the last paragraph.

We say that a matrix $\mathbf {\mathsf {U}}$ is orthonormal when its columns are orthonormal with respect to the standard inner product. That is, $\mathbf {\mathsf {U}}^{*} \mathbf {\mathsf {U}} = {\mathbf {\mathsf {I}}}$ . If $\mathbf {\mathsf {U}}$ is also square, we say instead that $\mathbf {\mathsf {U}}$ is orthogonal ( $\mathbb {F} = \mathbb {R}$ ) or unitary ( $\mathbb {F} = \mathbb {C}$ ).

2.4. Norms on matrices

Several different norms on matrices arise during this survey. We use consistent notation for these norms.

• • The unadorned norm ${\Vert \cdot \Vert }$ refers to the spectral norm of a matrix, also known as the $\ell _2$ operator norm. It reports the maximum singular value of its argument. For vectors, it coincides with the $\ell _2$ norm.

• • The norm ${\Vert \cdot \Vert }_{*}$ is the nuclear norm of a matrix, which is the dual of the spectral norm. It reports the sum of the singular values of its argument.

• • The symbol ${{\Vert \cdot \Vert }_{\mathrm {F}}}$ refers to the Frobenius norm, defined in the last subsection. The Frobenius norm coincides with the $\ell _2$ norm of the singular values of its argument.

• • The notation ${\Vert \cdot \Vert }_{p}$ denotes the Schatten p-norm for each $p \in [1, \infty ]$ . The Schatten p-norm is the $\ell _p$ norm of the singular values of its argument. Special cases with their own notation include the nuclear norm (Schatten $1$ ), the Frobenius norm (Schatten $2$ ) and the spectral norm (Schatten ${\infty }$ ).

Occasionally other norms may arise, and we will define them explicitly when they do.

2.5. Approximation in the spectral norm

Throughout this survey we will almost exclusively use the spectral norm to measure the error in matrix computations. Let us recall some of the implications that follow from spectral norm bounds.

Suppose that $\mathbf {\mathsf {A}} \in \mathbb {F}^{m \times n}$ is a matrix and $\widehat {\mathbf {\mathsf {A}}} \in \mathbb {F}^{m \times n}$ is an approximation. If the approximation satisfies the spectral norm error bound

$$ \begin{align*} {\Vert { \mathbf{\mathsf{A}} - \widehat{\mathbf{\mathsf{A}}} } \Vert} \leq \varepsilon, \end{align*} $$

then we can transfer the following information.

• • Linear functionals: ${\vert {{\langle {\mathbf {\mathsf {F}}}, \, {\mathbf {\mathsf {A}}} \rangle } - {\langle {\mathbf {\mathsf {F}}}, \, {\widehat {\mathbf {\mathsf {A}}}} \rangle } } \vert } \leq \varepsilon {\Vert \mathbf {\mathsf {F}} \Vert }_{*}$ for every matrix $\mathbf {\mathsf {F}} \in \mathbb {F}^{m \times n}$ .

• • Singular values: ${\vert { \sigma _j(\mathbf {\mathsf {A}}) - \sigma _j(\widehat {\mathbf {\mathsf {A}}}) } \vert } \leq \varepsilon $ for each index j.

• • Singular vectors: if the jth singular value $\sigma _j(\mathbf {\mathsf {A}})$ is well separated from the other singular values, then the jth right singular vector of $\mathbf {\mathsf {A}}$ is well approximated by the jth right singular vector of $\widehat {\mathbf {\mathsf {A}}}$ ; a similar statement holds for the left singular vectors.

Detailed statements about the singular vectors are complicated, so we refer the reader to Bhatia (Reference Bhatia1997, Chapters VII, X) for his treatment of perturbation of spectral subspaces.

For one-pass and streaming data models, it may not be possible to obtain good error bounds in the spectral norm. In this case we may retrench to Frobenius norm or nuclear norm error bounds. These estimates give weaker information about linear functionals, singular values and singular vectors.

Remark 2.1 (Frobenius norm approximation)

In the literature on randomized NLA, some authors prefer to bound errors with respect to the Frobenius norm because the arguments are technically simpler. In many instances, these bounds are less valuable because the error can have the same scale as the matrix that we wish to approximate.

For example, let us consider a variant of the spiked covariance model that is common in statistics applications (Johnstone Reference Johnstone2001). Suppose we need to approximate a rank-one matrix contaminated with additive noise: $\mathbf {\mathsf {A}} = \mathbf {\mathsf {uu}}^{*} + \varepsilon \mathbf {\mathsf {G}} \in {\mathbb {R}}^{n \times n}$ , where ${\Vert \mathbf {\mathsf {u}} \Vert } = 1$ and $\mathbf {\mathsf {G}} \in {\mathbb {R}}^{n \times n}$ has independent ${\small NORMAL}(0, n^{-1})$ entries. It is well known that ${\Vert \mathbf {\mathsf {G}} \Vert } \approx 2$ , while ${{\Vert {{\mathbf {\mathsf {G}}}} \Vert }_{\mathrm {F}}} \approx \sqrt {n}$ . With respect to the Frobenius norm, the zero matrix is almost as good an approximation of $\mathbf {\mathsf {A}}$ as the rank-one matrix $\mathbf {\mathsf {uu}}^{*}$ :

$$ \begin{align*} \operatorname{\mathbb{E}} \Vert {{{ \mathbf{\mathsf{A}} - \mathbf{\mathsf{uu}}^{*} }} \Vert}_{\mathrm{F}}^2 = \varepsilon^2 n \quad\text{and}\quad \operatorname{\mathbb{E}} \Vert {{ \mathbf{\mathsf{A}} - \mathbf{\mathsf{0}} }} \Vert_{\mathrm{F}}^2 = \varepsilon^2 n + 1. \end{align*} $$

The difference is visible only when the size of the perturbation $\varepsilon \approx n^{-1/2}$ . In contrast, the spectral norm error can easily distinguish between the good approximation $\mathbf {\mathsf {uu}}^{*}$ and the vacuous approximation $\mathbf {\mathsf {0}}$ , even when $\varepsilon = O(1)$ .

For additional discussion, see Tropp (Reference Tropp2015, Section 6.2.3) and Li et al. (Reference Li, Linderman, Szlam, Stanton, Kluger and Tygert2017, Appendix).

2.6. Intrinsic dimension and stable rank

Let $\mathbf {\mathsf {A}} \in \mathbb {H}_n$ be a PSD matrix. We define its intrinsic dimension:

(2.1) $$ \begin{align} \operatorname{intdim}(\mathbf{\mathsf{A}}) := \displaystyle \frac{\operatorname{trace}(\mathbf{\mathsf{A}})}{{\Vert \mathbf{\mathsf{A}} \Vert}}. \end{align} $$

The intrinsic dimension of a nonzero matrix satisfies the inequalities $1 \leq \operatorname {intdim}(\mathbf {\mathsf {A}}) \leq \operatorname {rank}(\mathbf {\mathsf {A}})$ ; the upper bound is saturated when $\mathbf {\mathsf {A}}$ is an orthogonal projector. We can interpret the intrinsic dimension as a continuous measure of the rank, or the number of energetic dimensions in the matrix.

Let $\mathbf {\mathsf {B}} \in \mathbb {F}^{m \times n}$ be a rectangular matrix. Its stable rank is

(2.2) $$ \begin{align} \operatorname{srank}(\mathbf{\mathsf{B}}) := \operatorname{intdim}(\mathbf{\mathsf{B}}^{*}\mathbf{\mathsf{B}}) = \displaystyle \frac{{{\Vert {{\mathbf{\mathsf{B}}}} \Vert}^2_{\mathrm{F}}}}{{\Vert \mathbf{\mathsf{B}} \Vert}^2}. \end{align} $$

Similar to the intrinsic dimension, the stable rank provides a continuous measure of the rank of $\mathbf {\mathsf {B}}$ .

2.7. Schur complements

Schur complements arise from partial Gaussian elimination and partial least-squares. They also play a key role in several parts of randomized NLA. We give the basic definitions here, referring to Zhang (Reference Zhang2005) for a more complete treatment.

Let $\mathbf {\mathsf {A}} \in \mathbb {F}^{n \times n}$ be a PSD matrix, and let $\mathbf {\mathsf {X}} \in \mathbb {F}^{n \times k}$ be a fixed matrix. First, define the PSD matrix

(2.3) $$ \begin{align} \mathbf{\mathsf{A}}{\langle {\mathbf{\mathsf{X}}} \rangle} := (\mathbf{\mathsf{AX}})(\mathbf{\mathsf{X}}^{*} \mathbf{\mathsf{A}} \mathbf{\mathsf{X}})^{\dagger} (\mathbf{\mathsf{AX}})^{*}. \end{align} $$

The Schur complement of $\mathbf {\mathsf {A}}$ with respect to $\mathbf {\mathsf {X}}$ is the PSD matrix

(2.4) $$ \begin{align} \mathbf{\mathsf{A}} / \mathbf{\mathsf{X}} := \mathbf{\mathsf{A}} - \mathbf{\mathsf{A}}\langle \mathbf{\mathsf{X}} \rangle. \end{align} $$

The matrices $\mathbf {\mathsf {A}}{\langle {\mathbf {\mathsf {X}}} \rangle }$ and $\mathbf {\mathsf {A}}/\mathbf {\mathsf {X}}$ depend on $\mathbf {\mathsf {X}}$ only through its range. They also enjoy geometric interpretations in terms of orthogonal projections with respect to the $\mathbf {\mathsf {A}}$ semi-inner-product.

2.8. Miscellaneous

We use big-O notation following standard computer science convention. For instance, we say that a method has (arithmetic) complexity $O (n^{\alpha })$ if there is a finite C for which the number of floating-point operations (flops) expended is bounded by $Cn^{\alpha }$ as the problem size $n\rightarrow \infty $ .

We use MATLAB-inspired syntax in summarizing algorithms. For instance, the task of computing an SVD $\mathbf {\mathsf {A}} = \mathbf {\mathsf {U}}\mathbf {\mathsf{\Sigma }}\mathbf {\mathsf {V}}^{*}$ of a given matrix $\mathbf {\mathsf {A}}$ is written as $[\mathbf {\mathsf {U}},\mathbf {\mathsf{\Sigma }},\mathbf {\mathsf {V}}]=\texttt {svd}(\mathbf {\mathsf {A}})$ . We have taken the liberty to modify the syntax when we believe that this improves clarity. For instance, we write $[\mathbf {\mathsf {Q}},\mathbf {\mathsf {R}}] = \texttt {qr_econ}(\mathbf {\mathsf {A}})$ to denote the economy-size QR factorization where the matrix $\mathbf {\mathsf {Q}}$ has size $m\times \min (m,n)$ for an input matrix $\mathbf {\mathsf {A}} \in \mathbb {F}^{m\times n}$ . Arguments that are not needed are replaced by ‘ $\sim $ ’, so that, for example, $[\mathbf {\mathsf {Q}},\sim ] = \texttt {qr_econ}(\mathbf {\mathsf {A}})$ returns only the matrix $\mathbf {\mathsf {Q}}$ whose columns form an orthonormal basis for the range of $\mathbf {\mathsf {A}}$ .

3.1. Basics

We work in a master probability space that is rich enough to support all of the random variables that are defined. We will not comment further about the underlying model.

In this paper the unqualified term random variable encompasses random scalars, vectors and matrices. Scalar-valued random variables are usually (but not always) denoted by uppercase italic Roman letters ( $X, Y, Z$ ). A random vector is denoted by a lowercase bold letter ( $\mathbf {\mathsf {x}}, \mathbf {{\omega }}$ ). A random matrix is denoted by an uppercase bold letter ( $\mathbf {\mathsf {X}}, \mathbf {\mathsf {Y}}, \mathbf {\mathsf{\Gamma }}, \mathbf {\mathsf{\Omega }}$ ). This notation works in concert with the notation for deterministic vectors and matrices.

The map $\mathbb {P}(E)$ returns the probability of an event E. We usually specify the event using the compact set builder notation that is standard in probability. For example, ${\mathbb {P}\{ { X> t } \}}$ is the probability that the scalar random variable X exceeds a level t.

The operator $\operatorname {\mathbb {E}}$ returns the expectation of a random variable. For vectors and matrices, the expectation can be computed coordinate by coordinate. The expectation is linear, which justifies relations like

$$ \begin{align*} \operatorname{\mathbb{E}}[ \mathbf{\mathsf{AX}} ] = \mathbf{\mathsf{A}} \operatorname{\mathbb{E}}[ \mathbf{\mathsf{X}} ] \quad\text{when }\mathbf{\mathsf{A}}\text{ is deterministic and }\mathbf{\mathsf{X}}\text{ is random.} \end{align*} $$

We use the convention that nonlinear functions bind before the expectation; for instance, $\operatorname {\mathbb {E}} X^2 = \operatorname {\mathbb {E}}[ X^2 ]$ . The operator $\operatorname {Var}[\cdot ]$ returns the variance of a scalar random variable.

We say that a random variable is centred when its expectation equals zero. A random vector $\mathbf {\mathsf {x}}$ is isotropic when $\operatorname {\mathbb {E}}[ \mathbf {\mathsf {xx}}^{*} ] = {\mathbf {\mathsf {I}}}$ . A random vector is standardized when it is both centred and isotropic. In particular, a scalar random variable is standardized when it has expectation zero and variance one.

When referring to independent random variables, we often include the qualification ‘statistically independent’ to make a distinction with ‘linearly independent’. We abbreviate the term (statistically) ‘independent and identically distributed’ as i.i.d.

3.2. Distributions

To refer to a named distribution, we use small capitals. In this context, the symbol $\sim $ means ‘has the same distribution as’.

We write ${\small\rm UNIF}$ for the uniform distribution over a finite set (with counting measure). In particular, a scalar Rademacher random variable has the distribution ${\small\rm UNIF}\{ \pm 1 \}$ . A Rademacher random vector has i.i.d. coordinates, each distributed as a scalar Rademacher random variable. We sometimes require the uniform distribution over a Borel subset of $\mathbb {F}^n$ , equipped with Lebesgue measure.

We write ${\small\rm NORMAL}(\mathbf {{\mu }}, \mathbf {\mathsf {C}})$ for the normal distribution on $\mathbb {F}^n$ with expectation $\mathbf {{\mu }} \in \mathbb {F}^n$ and PSD covariance matrix $\mathbf {\mathsf {C}} \in \mathbb {H}_n(\mathbb {F})$ . A standard normal random variable or random vector has expectation zero and covariance matrix equal to the identity. We often use the term Gaussian to refer to normal distributions.

3.3. Concentration inequalities

Concentration inequalities provide bounds on the probability that a random variable is close to its expectation. A good reference for the scalar case is the book by Boucheron, Lugosi and Massart (Reference Boucheron, Lugosi and Massart2013). In the matrix setting, closeness is measured in the spectral norm. For an introduction to matrix concentration, see Tropp (Reference Tropp2015, Reference Tropp, Yurtsever, Udell and Cevher2019). These results play an important role in randomized linear algebra.

3.4. Gaussian random matrix theory

On several occasions we use comparison principles to study the action of a random matrix with i.i.d. Gaussian entries. In particular, these methods can be used to control the largest and smallest singular values. The main classical comparison theorems are associated with the names Slepian, Chevet and Gordon. For accounts of this, see Ledoux and Talagrand (Reference Ledoux and Talagrand1991, Section 3.3), Davidson and Szarek (Reference Davidson and Szarek2001, Section 2.3) or Vershynin (Reference Vershynin2018, Sections 7.2–7.3). More recently, it has been observed that Gordon’s inequality can be reversed in certain settings (Thrampoulidis, Oymak and Hassibi Reference Thrampoulidis, Oymak and Hassibi2014).

In several instances we require more detailed information about Gaussian random matrices. General resources include Muirhead (Reference Muirhead1982) and Bai and Silverstein (Reference Bai and Silverstein2010). Most of the specific results we need are presented in Halko et al. (Reference Halko, Martinsson and Tropp2011a ).

4.1. Overview

We will focus on the problem of estimating the trace of a nonzero PSD matrix $\mathbf {\mathsf {A}} \in \mathbb {H}_n$ . Our goal is to to produce an approximation of $\operatorname {trace}(\mathbf {\mathsf {A}})$ , along with a measure of quality.

Trace estimation is easy in the case where we have inexpensive access to the entries of the matrix $\mathbf {\mathsf {A}}$ because we can simply read off the n diagonal entries. But there are many environments where the primitive operation is the matrix–vector product $\mathbf {\mathsf {u}} \mapsto \mathbf {\mathsf {A}}\mathbf {\mathsf {u}}$ . For example, $\mathbf {\mathsf {u}} \mapsto \mathbf {\mathsf {A}}\mathbf {\mathsf {u}}$ might be the solution of a (discretized) linear differential equation with initial condition $\mathbf {\mathsf {u}}$ , implemented by some computer program. In this case we would really prefer to avoid n applications of the primitive. (We can obviously compute the trace by applying the primitive to each standard basis vector $\mathbf {{\delta }}_i$ .)

The methods in this section all use linear information about the matrix $\mathbf {\mathsf {A}}$ . In other words, we will extract data from the input matrix by computing the product $\mathbf {\mathsf {Y}} = \mathbf {\mathsf {A}}\mathbf {\mathsf{\Omega }}$ , where $\mathbf {\mathsf{\Omega }} \in \mathbb {F}^{n \times k}$ is a (random) test matrix. All subsequent operations involve only the sample matrix $\mathbf {\mathsf {Y}}$ and the test matrix $\mathbf {\mathsf{\Omega }}$ . Since the data collection process is linear, we can apply randomized trace estimators in the one-pass or streaming environments. Moreover, parts of these algorithms are trivially parallelizable.

The original application of randomized trace estimation was to perform a posteriori error estimation for large least-squares computations. More specifically, it was used to accelerate cross-validation procedures for estimating the optimal regularization parameter in a smoothing spline (Girard Reference Girard1989, Hutchinson Reference Hutchinson1990). See Fitzsimons, Osborne, Roberts and Fitzsimons (Reference Fitzsimons, Osborne, Roberts, Fitzsimons, Globerson and Silva2018) for a list of contemporary applications in machine learning, uncertainty quantification and other fields.

4.2. Trace estimation by randomized sampling

Randomized trace estimation is based on the insight that it is easy to construct a random variable whose expectation equals the trace of the input matrix.

Consider a random test vector $\mathbf {{\omega }} \in \mathbb {F}^n$ that is isotropic: $\operatorname {\mathbb {E}}[ \mathbf {{\omega }} {\mathbf {{\omega }}}^{*} ] = {\mathbf {\mathsf {I}}}$ . By the cyclicity of the trace and by linearity,

(4.1) $$ \begin{align} X = {\mathbf{{\omega}}}^{*} (\mathbf{\mathsf{A}} \mathbf{{\omega}}) \quad\text{satisfies}\quad \operatorname{\mathbb{E}} X = \operatorname{trace}(\mathbf{\mathsf{A}}). \end{align} $$

In other words, the random variable X is an unbiased estimator of the trace. Note that the distribution of X depends on the unknown matrix $\mathbf {\mathsf {A}}$ .

A single sample of X is rarely adequate because its variance, $\operatorname {Var}[X]$ , will be large. The most common mechanism for reducing the variance is to average k independent copies of X. For $k \in \mathbb {N}$ , define

(4.2) $$ \begin{align} \bar{X}_k = \displaystyle \frac{1}{k} \sum\nolimits_{i=1}^k X_i \quad\text{where }X_i \sim X\text{ are i.i.d.} \end{align} $$

By linearity, $\bar {X}_k$ is also an unbiased estimator of the trace. The individual samples are statistically independent, so the variance decreases. Indeed,

$$ \begin{align*} \operatorname{\mathbb{E}}[ \bar{X}_k ] = \operatorname{trace}(\mathbf{\mathsf{A}}) \quad\text{and}\quad \operatorname{Var}[ \bar{X}_k ] = \displaystyle \frac{1}{k} \operatorname{Var}[X]. \end{align*} $$

The estimator (4.2) can be regarded as the most elementary method in randomized linear algebra. See Algorithm 1.

To compute $\bar {X}_k$ , we must simulate k independent copies of the random vector $\mathbf {{\omega }} \in \mathbb {F}^n$ and perform k matrix–vector products with $\mathbf {\mathsf {A}}$ , plus $O (kn)$ additional arithmetic.

Example 4.1 (Girard Reference Girard1989)

Consider a standard normal random vector $\mathbf {{\omega }} \sim {\small\rm NORMAL}(\mathbf {\mathsf {0}}, {\mathbf {\mathsf {I}}})$ . The variance of the resulting trace estimator (4.1)–(4.2) satisfies

(4.3) $$ \begin{align} \operatorname{Var}[ \bar{X}_k ] = \displaystyle \frac{2}{k} \sum\nolimits_{i,j = 1}^n {\vert { (\mathbf{\mathsf{A}})_{ij} } \vert}^2 = \displaystyle \frac{2}{k} {{\Vert {{\mathbf{\mathsf{A}}}} \Vert}^2_{\mathrm{F}}} \leq \displaystyle \frac{2}{k} {\Vert \mathbf{\mathsf{A}} \Vert} \operatorname{trace}(\mathbf{\mathsf{A}}).\ \ \end{align} $$

The rotational invariance of the standard normal distribution allows us to characterize the behaviour of this estimator in full detail.

Example 4.2 (Hutchinson Reference Hutchinson1990)

Consider a Rademacher random vector $\mathbf {{\omega }} \sim {\small\rm UNIF}\{\pm 1\}^n$ . The variance of the resulting trace estimator (4.1)–(4.2) satisfies

$$ \begin{align*} {\operatorname{Var}[ \bar{X}_k ] = \displaystyle \frac{4}{k} \sum\nolimits_{1 \leq i < j \leq n} {\vert { (\mathbf{\mathsf{A}})_{ij} } \vert}^2 < \displaystyle \frac{2}{k} {{\Vert {{\mathbf{\mathsf{A}}}} \Vert}^2_{\mathrm{F}}} \leq \displaystyle \frac{2}{k} {\Vert \mathbf{\mathsf{A}} \Vert} \operatorname{trace}(\mathbf{\mathsf{A}}).} \end{align*} $$

This is the minimum variance trace estimator generated by an isotropic random vector $\mathbf {{\omega }}$ with statistically independent coordinates. It also avoids the simulation of normal variables.

Girard (Reference Girard1989) also studied the estimator obtained by drawing $\mathbf {{\omega }} \in \mathbb {F}^n$ uniformly at random from the sphere $\sqrt {n} \, \mathbb {S}^{n-1}(\mathbb {F})$ for $\mathbb {F} = {\mathbb {R}}$ . When $\mathbb {F} = \mathbb {C}$ , this approach has the minimax variance among all trace estimators of the form (4.1)–(4.2). We return to this example in Section 4.7.1.

4.3. A priori error estimates

We can use theoretical analysis to obtain prior guarantees on the performance of the trace estimator. These results illuminate what features of the input matrix affect the quality of the trace estimate, and they tell us how many samples k suffice to achieve a given error tolerance. Note, however, that these bounds depend on properties of the input matrix that are often unknown to the user of the trace estimator.

Regardless of the distribution of the isotropic test vector $\mathbf {{\omega }}$ , Chebyshev’s inequality delivers a simple probability bound for the trace estimator:

(4.4) $$ \begin{align} {\mathbb{P}\{ { {\vert { \bar{X}_k - \operatorname{trace}(\mathbf{\mathsf{A}}) } \vert} \geq t } \}} \leq \displaystyle \frac{ \operatorname{Var}[X] }{k t^2} \quad\text{for }t> 0. \end{align} $$

We can specialize this result to specific trace estimators by inserting the variance.

Example 4.4 (Girard trace estimator)

If the test vector $\mathbf {{\omega }}$ is standard normal, the trace estimator $\bar {X}_k$ satisfies

$$ \begin{align*} {\mathbb{P}\{ { {\vert { \bar{X}_k - \operatorname{trace}(\mathbf{\mathsf{A}}) } \vert} \geq t \cdot \operatorname{trace}(\mathbf{\mathsf{A}}) } \}} \leq \displaystyle \frac{ 2 }{k \operatorname{intdim}(\mathbf{\mathsf{A}}) \, t^2}. \end{align*} $$

The bound follows from (4.3), (4.4) and (2.1). In words, the trace estimator achieves a relative error bound that is sharpest when the intrinsic dimension (2.1) of $\mathbf {\mathsf {A}}$ is large.

For specific distributions of the random test vector $\mathbf {{\omega }}$ , we can obtain much stronger probability bounds for the resulting trace estimator using exponential concentration inequalities. Here is a recent analysis for Girard’s estimator based on fine properties of the standard normal distribution.

Theorem 4.5 (Gratton and Titley-Peloquin Reference Gratton and Titley-Peloquin2018)

Let $\mathbf {\mathsf {A}} \in \mathbb {H}_n({\mathbb {R}})$ be a nonzero PSD matrix. Consider the trace estimator (4.1)–(4.2) obtained from a standard normal test vector $\mathbf {{\omega }} \in {\mathbb {R}}^n$ . For $\tau> 1$ and $k \leq n$ ,

$$ \begin{align*} {\mathbb{P}\{ { \bar{X}_k \geq \tau \operatorname{trace}(\mathbf{\mathsf{A}}) } \}} &\leq \exp\bigl( -\tfrac{1}{2} k \operatorname{intdim}(\mathbf{\mathsf{A}}) (\sqrt{\tau} - 1 )^2 \bigr){,} \\ {\mathbb{P}\{ { \bar{X}_k \leq \tau^{-1} \operatorname{trace}(\mathbf{\mathsf{A}}) } \}} &\leq \exp\bigl( -\tfrac{1}{4} k \operatorname{intdim}(\mathbf{\mathsf{A}}) (\tau^{-1} - 1 )^2 \bigr). \end{align*} $$

When $\mathbf {\mathsf {A}} \in \mathbb {H}_n(\mathbb {C})$ is PSD and $\mathbf {{\omega }} \in \mathbb {C}^n$ is complex standard normal, the same bounds hold with an extra factor two in the exponent. (So the estimator works better in the complex setting.)

4.4. Universality

Empirically, for a large sample, the performance of the trace estimator $\bar {X}_k$ only depends on the distribution of the test vector $\mathbf {{\omega }}$ through the variance of the resulting sample X. In a word, the estimator exhibits universality. As a consequence, we can select the distribution that is most convenient for computational purposes.

Classical probability theory furnishes justification for these claims. The strong law of large numbers tells us that

$$ \begin{align*} \bar{X}_k \to \operatorname{trace}(\mathbf{\mathsf{A}}) \quad\text{almost surely as }k \to \infty. \end{align*} $$

Concentration inequalities (Boucheron et al. Reference Boucheron, Lugosi and Massart2013) allow us to derive rates of convergence akin to Theorem 4.5.

To understand the sampling distribution of the estimator $\bar {X}_k$ , we can invoke the central limit theorem:

$$ \begin{align*} \sqrt{k} ( \bar{X}_k - \operatorname{trace}(\mathbf{\mathsf{A}}) ) \to {\small\rm NORMAL}(0, \operatorname{Var}[X]) \quad\text{in distribution as }k \to \infty. \end{align*} $$

We can obtain estimates for the rate of convergence to normality using the Berry–Esséen theorem and its variants (Ross Reference Ross2011, Chen, Goldstein and Shao Reference Chen, Goldstein and Shao2011).

Owing to the universality phenomenon, we can formally use the normal limit to obtain heuristic error estimates and insights for trace estimators constructed with test vectors from any distribution. This strategy becomes even more valuable when the linear algebra problem is more complicated.

4.5. A posteriori error estimates

In practice, we rarely have access to all the information required to activate a priori error bounds. It is wiser to assess the quality of the estimate from the information that we actually collect. Since we have full knowledge of the random process that generates the trace estimate, we can confidently use approaches from classical statistics.

At the most basic level, the sample variance is an unbiased estimator for the variance of the individual samples:

$$ \begin{align*} S_k = \displaystyle \frac{1}{k-1} \sum\nolimits_{i=1}^k (X_i - \bar{X}_k)^2 \quad\text{satisfies}\quad \operatorname{\mathbb{E}}[ S_k ] = \operatorname{Var}[X]. \end{align*} $$

The variance of $S_k$ depends on the fourth moment of the random variable X. A standard estimate is

$$ \begin{align*} \operatorname{Var}[S_k] \leq \displaystyle \frac{1}{k} \operatorname{\mathbb{E}}[ (X - \operatorname{\mathbb{E}} X)^4 ]. \end{align*} $$

Bounds and empirical estimates for the variance of $S_k$ can also be obtained using the Efron–Stein inequality (Boucheron et al. Reference Boucheron, Lugosi and Massart2013, Section 3.1).

For $\alpha \in (0, 1/2)$ , we can construct the (symmetric, Student’s t) confidence interval at level $1 - 2\alpha $ :

$$ \begin{align*} \operatorname{trace}(\mathbf{\mathsf{A}}) \in \bar{X}_k \pm t_{\alpha, k-1} \sqrt{S_k}{,} \end{align*} $$

where $t_{\alpha , k-1}$ is the $\alpha $ quantile of the Student’s t-distribution with $k-1$ degrees of freedom. We interpret this result as saying that $\operatorname {trace}(\mathbf {\mathsf {A}})$ lies in the specified interval with probability roughly $1 - 2\alpha $ (over the randomness in the trace estimator). The usual rule of thumb is that the sample size should be moderate (say, $k \geq 30$ ), while $\alpha $ cannot be too small (say, $\alpha \geq 0.025$ ).

4.6. Bootstrapping the sampling distribution

Miles Lopes has proposed a sweeping program that uses the bootstrap to construct data-driven confidence sets for randomized NLA algorithms (Lopes Reference Lopes2019). For trace estimation, this approach is straightforward to describe and implement; see Algorithm 2.

Let $\mathcal {X} = ( X_1, \ldots , X_k )$ be the empirical sample from (4.2). The bootstrap draws further random samples from $\mathcal {X}$ to elicit more information about the sampling distribution of the trace estimator X, such as confidence sets.

1. (1) For each $b = 1, \ldots , B$ :

   • – Draw a bootstrap replicate $( X_1^{*}, \ldots , X_k^{*} )$ uniformly from $\mathcal {X}$ with replacement.

   • – Compute the error estimate $e_b^{*} = \bar {X}_k^{*} - \bar {X}_k$ , where $\bar {X}_k^{*}$ is the sample average of the bootstrap replicate.

2. (2) Compute quantiles $q_\alpha $ and $q_{1-\alpha }$ of the error distribution $( e_1^{*}, \ldots , e_B^{*} )$ .

3. (3) Report the $1-2\alpha $ confidence set $[\bar {X}_k + q_{\alpha }, \bar {X}_k + q_{1 - \alpha }]$ .

Typical values are $k \geq 30$ samples and $B \geq 1000$ bootstrap replicates when $\alpha \geq 0.025$ . This method is effective for a wide range of distributions on the test vector, and it extends to other problems. See Efron (Reference Efron1982) for an introduction to resampling methods.

4.7. Structured distributions for test vectors

As discussed, there is a lot of flexibility in designing the distribution of the test vector. We can exploit this freedom to achieve additional computational goals. For example, we might:

• • minimize the variance $\operatorname {Var}[X]$ of each sample;

• • minimize the number of random bits required to construct $\mathbf {{\omega }}$ ;

• • design a test vector $\mathbf {{\omega }}$ that is ‘compatible’ with the input matrix $\mathbf {\mathsf {A}}$ to facilitate the matrix–vector product. For example, if $\mathbf {\mathsf {A}}$ has a tensor product structure, we might require $\mathbf {{\omega }}$ to share the tensor structure.

Let us describe a general construction for test vectors that can help achieve these desiderata. The ideas come from frame theory and quantum information theory. The approach here extends the work in Fitzsimons et al. (Reference Fitzsimons, Osborne, Roberts, Fitzsimons, Globerson and Silva2018).

4.7.1. Optimal measurement systems

In this section we work in the complex field. Consider a discrete set $\mathcal {U} := \{ \mathbf {\mathsf {u}}_1, \ldots , \mathbf {\mathsf {u}}_m \} \subset \mathbb {C}^n$ of vectors, each with unit $\ell _2$ norm. We say that the $\mathcal {U}$ is an optimal measurement system when

(4.5) $$ \begin{align} \displaystyle \frac{1}{m} \sum\nolimits_{i=1}^m ({\mathbf{\mathsf{u}}}_i^{*}\mathbf{\mathsf{M}} {\mathbf{\mathsf{u}}}_i) \, {\mathbf{\mathsf{u}}}_i {\mathbf{\mathsf{u}}}_i^{*} = \displaystyle \frac{1}{(n+1) n} [ \mathbf{\mathsf{M}} + \operatorname{trace}(\mathbf{\mathsf{M}}) \, {\mathbf{\mathsf{I}}} ] \end{align} $$

for all $\mathbf {\mathsf {M}} \in \mathbb {H}_n(\mathbb {C})$ . The reproducing property (4.5) shows that the system of vectors acquires enough information to reconstruct an arbitrary self-adjoint matrix. A similar definition is valid for an infinite system of unit vectors, provided we replace the sum in (4.5) with an integral.

Now, suppose that we draw a random test vector $\mathbf {{\omega }} = \sqrt {n} \, \mathbf {\mathsf {u}}$ , where $\mathbf {\mathsf {u}}$ is drawn uniformly at random from an optimal measurement system $\mathcal {U}$ . Then the resulting trace estimator is unbiased:

(4.6) $$ \begin{align} X = {\mathbf{{\omega}}}^{*} (\mathbf{\mathsf{A}} \mathbf{{\omega}}) \quad\text{satisfies}\quad \operatorname{\mathbb{E}} X = \operatorname{trace}(\mathbf{\mathsf{A}}). \end{align} $$

The variance of this trace estimator satisfies

(4.7) $$ \begin{align} \operatorname{Var}[X] = \displaystyle \frac{n}{n+1} \biggl[ {{\Vert {{\mathbf{\mathsf{A}}}} \Vert}^2_{\mathrm{F}}} - \displaystyle \frac{1}{n} \operatorname{trace}(\mathbf{\mathsf{A}})^2 \biggr]. \end{align} $$

The identities (4.6) and (4.7) follow quickly from (4.5). As it happens, this is the minimax variance achievable for a best isotropic distribution on test vectors (Kueng Reference Kueng2019)

4.8. Extension: The Frobenius norm and the Schatten 4-norm

The randomized trace estimators developed in this section can also be deployed to estimate a couple of matrix norms.

Consider a rectangular matrix $\mathbf {\mathsf {B}} \in \mathbb {F}^{m \times n}$ , accessed via the matrix–vector product $\mathbf {\mathsf {u}} \mapsto \mathbf {\mathsf {B}} \mathbf {\mathsf {u}}$ . Let us demonstrate how to estimate the Frobenius norm (i.e. Schatten 2-norm) and the Schatten 4-norm of $\mathbf {\mathsf {B}}$ .

For concreteness, suppose that we extract test vectors from the standard normal distribution. Draw a standard normal matrix $\mathbf {\mathsf{\Omega }} \in \mathbb {F}^{n \times k}$ with columns $\mathbf {{\omega }}_i$ . Construct the random variable

$$ \begin{align*} \bar{X}_k := \displaystyle \frac{1}{k} {{\Vert {{ \mathbf{\mathsf{B}} \mathbf{\mathsf{\Omega}} }} \Vert}^2_{\mathrm{F}}} = \displaystyle \frac{1}{k} \sum\nolimits_{i=1}^k {\mathbf{{\omega}}}_i^{*} (\mathbf{\mathsf{B}}^{*} \mathbf{\mathsf{B}}) \mathbf{{\omega}}_i =: \displaystyle \frac{1}{k} \sum\nolimits_{i=1}^k X_i. \end{align*} $$

We compute $\bar {X}_k$ by simulating $nk$ standard normal variables, taking k matrix–vector products with $\mathbf {\mathsf {B}}$ , and performing ${{O}}(km)$ additional arithmetic.

To analyse $\bar {X}_k$ , note that it is an instance of the randomized trace estimator, where $\mathbf {\mathsf {A}} = \mathbf {\mathsf {B}}^{*} \mathbf {\mathsf {B}}$ . In particular, its statistics are

$$ \begin{align*} \operatorname{\mathbb{E}}[ \bar{X}_k ] = {{\Vert {{\mathbf{\mathsf{B}}}} \Vert}^2_{\mathrm{F}}} \quad\text{and}\quad \operatorname{Var}[ \bar{X}_k ] = \displaystyle \frac{2}{k} {\Vert \mathbf{\mathsf{B}} \Vert}_{4}^4. \end{align*} $$

We see that $\bar {X}_k$ provides an unbiased estimate for the squared Frobenius norm of the matrix $\mathbf {\mathsf {B}}$ . Meanwhile, by rescaling the sample variance $S_k^2$ of the data $( X_1, \ldots , X_k)$ , we obtain an unbiased estimate for the fourth power of the Schatten 4-norm of $\mathbf {\mathsf {B}}$ .

Our discussion shows how we can obtain a priori guarantees and a posteriori error estimates for these norm computations; the results can be expressed in terms of the stable rank (2.2) of $\mathbf {\mathsf {B}}$ . We can also obtain norm estimators that are more computationally efficient using structured distributions for the test vectors, such as elements from an optimal measurement system.

5.1. Overview

Consider a general matrix $\mathbf {\mathsf {B}} \in \mathbb {F}^{m \times n}\!$ , accessed via the matrix–vector product $\mathbf {\mathsf {u}} \mapsto \mathbf {\mathsf {B}}\mathbf {\mathsf {u}}$ . For a sample size k, let $\mathbf {\mathsf{\Omega }} \in \mathbb {F}^{n \times k}$ be a (random) test matrix that does not depend on $\mathbf {\mathsf {B}}$ . For a natural number $p \geq 3$ , we consider the problem of estimating the Schatten $2p$ -norm ${\Vert \mathbf {\mathsf {B}} \Vert }_{2p}$ from the sample matrix $\mathbf {\mathsf {Y}} = \mathbf {\mathsf {B}}\mathbf {\mathsf{\Omega }}$ . The key idea is to cleverly process the sample matrix $\mathbf {\mathsf {Y}}$ to form an unbiased estimator for the $2p$ th power of the norm. Since the methods in this section use linear information, they can be parallelized, and they are applicable in the one-pass and streaming environments.

Of course, if we are given a singular value decomposition (SVD) of $\mathbf {\mathsf {B}}$ , it is straightforward to extract the Schatten $2p$ -norm. We are interested in methods that are much less expensive than the $O (\min \{mn^2, nm^2\})$ cost of computing an SVD with a classical direct algorithm. Randomized SVD or URV algorithms (Sections 11.2, 15 and 16) can also be used for Schatten norm estimation, but this approach is typically overkill.

5.2. Interlude: Lower bounds

How large a sample size k is required to estimate ${\Vert \mathbf {\mathsf {B}} \Vert }_{2p}$ up to a fixed constant factor with 75% probability over the randomness in $\mathbf {\mathsf{\Omega }}$ ? The answer, unfortunately, turns out to be . In other words, for a general matrix, we cannot approximate the Schatten $2p$ -norm for any $p> 2$ unless the sample size k grows polynomially with the dimension.

A more detailed version of this statement appears in Li, Nguyen and Woodruff (Reference Woodruff2014a , Theorem 3.2). The authors exhibit a particularly difficult type of input matrix (a standard normal matrix with a rank-one spike, as in Remark 2.1) to arrive at the negative conclusion.

If you have a more optimistic nature, you can also take the inspiration that other types of input matrices might be easier to handle. For example, it is possible to use a small random sample to compute the Schatten norm of a matrix that enjoys some decay in the singular value spectrum.

5.3. Estimating Schatten norms the hard way

First, let us describe a technique from classical statistics that leads to an unbiased estimator of ${\Vert \mathbf {\mathsf {B}} \Vert }_{2p}^{2p}$ . This estimator is both highly variable and computationally expensive, so we must proceed with caution.

For the rest of this section, we assume that the random test matrix $\mathbf {\mathsf{\Omega }} \in \mathbb {F}^{n \times k}$ has isotropic columns $\mathbf {{\omega }}_i$ that are i.i.d. Form the sample matrix $\mathbf {\mathsf {Y}} = \mathbf {\mathsf {B}}\mathbf {\mathsf{\Omega }}$ . Abbreviate $\mathbf {\mathsf {A}} = \mathbf {\mathsf {B}}^{*}\mathbf {\mathsf {B}}$ and $\mathbf {\mathsf {X}} = \mathbf {\mathsf {Y}}*\mathbf {\mathsf {Y}}$ . Observe that

$$ \begin{align*} (\mathbf{\mathsf{X}})_{ij} = ({\mathbf{\mathsf{Y}}^{*}} \mathbf{\mathsf{Y}})_{ij} = {\mathbf{{\omega}}}_i^{*} \mathbf{\mathsf{A}} {\mathbf{{\omega}}}_j. \end{align*} $$

Therefore, for any natural numbers that satisfy $1 \leq i_1, \ldots , i_p \leq k$ ,

$$ \begin{align*} (\mathbf{\mathsf{X}})_{i_1 i_2} ({\mathbf{\mathsf{X}})}_{i_2 i_3} \cdots {(\mathbf{\mathsf{X}})}_{i_p i_1} = \operatorname{trace}[ {\mathbf{{\omega}}}_{i_1} {\mathbf{{\omega}}}_{i_1}^{*} \mathbf{\mathsf{A}} \cdots {\mathbf{{\omega}}}_{i_p} {\mathbf{{\omega}}}_{i_p}^{*} \mathbf{\mathsf{A}} ]. \end{align*} $$

If we assume that $i_1, \ldots , i_p$ are distinct, we can use independence and isotropy to compute the expectation:

$$ \begin{align*} \operatorname{\mathbb{E}}[ {(\mathbf{\mathsf{X}})}_{i_1 i_2} {(\mathbf{\mathsf{X}})}_{i_2 i_3} \cdots {(\mathbf{\mathsf{X}})}_{i_p i_1} ] = \operatorname{trace}[ {\mathbf{\mathsf{A}}}^p ] = {\Vert \mathbf{\mathsf{B}} \Vert}_{2p}^{2p} . \end{align*} $$

By averaging over all sequences of distinct indices, we obtain an unbiased estimator:

$$ \begin{align*} U_p = \displaystyle \frac{(k-p)!}{k!} \sum\nolimits_{1 \leq i_1, \ldots, i_p \leq k}^{\circ} {(\mathbf{\mathsf{X}})}_{i_1 i_2} {(\mathbf{\mathsf{X}})}_{i_2 i_3} \cdots {(\mathbf{\mathsf{X}})}_{i_p i_1}. \end{align*} $$

The circle over the sum indicates that the indices must be distinct. Since $\operatorname {\mathbb {E}} U_p = {\Vert \mathbf {\mathsf {B}} \Vert }_{2p}^{2p}$ , our hope is that

$$ \begin{align*} U_p^{1/(2p)} \approx {\Vert \mathbf{\mathsf{B}} \Vert}_{2p}. \end{align*} $$

To ensure that the approximation is precise, the standard deviation of $U_p$ should be somewhat smaller than the mean of $U_p$ .

To understand the statistic $U_p$ , we can use tools from the theory of U-statistics (Koroljuk and Borovskich Reference Koroljuk and Borovskich1994). For instance, when p is fixed, we have the limit

$$ \begin{align*} k \operatorname{Var}[ U_p ] \to p^2 \operatorname{Var}[ {\mathbf{{\omega}}}^{*} {\mathbf{\mathsf{A}}}^p \mathbf{{\omega}} ] \quad\text{as }k \to \infty. \end{align*} $$

In particular, if the test vector $\mathbf {{\omega }} \sim {\small\rm NORMAL}(\mathbf {\mathsf {0}}, {\mathbf {\mathsf {I}}})$ , then

$$ \begin{align*} k \operatorname{Var}[ U_p ] \to 2p^2 {\Vert \mathbf{\mathsf{B}} \Vert}_{4p}^{4p} \quad\text{as }k \to \infty. \end{align*} $$

We can reduce the variance further by using test vectors from an optimal measurement system.

Unfortunately, it is quite expensive to compute the statistic $U_p$ because it involves almost $k^p$ summands. When p is a small constant (say, $p = 4$ or $p = 5$ ), it is not too onerous to form $U_p$ . On the other hand, in the worst case, we cannot beat the lower bound , where computation of $U_p$ requires ${{O}}(\min \{m,n\}^p)$ operations. Good (Reference Good1977) proposes some small economies in this computation.

5.4. Estimating Schatten norms the easy way

Next, we describe a more recent approach that was proposed by Kong and Valiant (Reference Kong and Valiant2017). This method suffers from even higher variance, but it is computationally efficient. As a consequence, we can target larger values of p and exploit a larger sample size k.

Superficially, the Kong and Valiant estimator appears similar to the statistic $U_p$ . With the same notation, they restrict their attention to increasing sequences $i_1 < i_2 < \cdots < i_p$ of indices. In other words,

$$ \begin{align*} V_p = \binom{k}{p}^{-1} \sum\nolimits_{1 \leq i_1 < \cdots < i_p \leq k} {(\mathbf{\mathsf{X}})}_{i_1 i_2} {(\mathbf{\mathsf{X}})}_{i_2 i_3} \cdots {(\mathbf{\mathsf{X}})}_{i_p i_1}. \end{align*} $$

Much as before, $V_p$ gives an unbiased estimator for ${\Vert \mathbf {\mathsf {B}} \Vert }_{2p}^{2p}$ .

Although $V_p$ still appears to be combinatorial, the restriction to increasing sequences allows for a linear algebraic reformulation of the statistic. Let $\mathcal {T} \colon \mathbb {H}_k \to \mathbb {F}^{k \times k}$ be the linear map that reports the strict upper triangle of a self-adjoint matrix. Then

$$ \begin{align*} V_p = \binom{k}{p}^{-1} \operatorname{trace}[ \mathcal{T}(\mathbf{\mathsf{X}})^{p-1} \mathbf{\mathsf{X}} ]. \end{align*} $$

The cost of computing $V_p$ is usually dominated by the ${{O}}(k^2 n)$ arithmetic required to form $\mathbf {\mathsf {X}}$ given $\mathbf {\mathsf {Y}}$ . See Algorithm 3 for the procedure.

Kong and Valiant (Reference Kong and Valiant2017) obtain bounds for the variance of the estimator $V_p$ to justify its employment when the number k of samples satisfies . This bound is probably substantially pessimistic for matrices that exhibit spectral decay, but these theoretical and computational questions remain open.

5.5. Bootstrapping the sampling distribution

Given the lack of prior guarantees for the estimators $U_p$ and $V_p$ described in this section, we recommend that users apply resampling methods to obtain empirical information about the sampling distribution. See Arcones and Gine (Reference Arcones and Gine1992) for reliable bootstrap procedures for U-statistics.

5.6. Extension: Estimating the spectral norm by random sampling

Recall that the spectral norm of $\mathbf {\mathsf {B}}$ is comparable with the Schatten $2p$ -norm of $\mathbf {\mathsf {B}}$ for an appropriate choice of p. Indeed,

$$ \begin{align*} {\Vert \mathbf{\mathsf{B}} \Vert}_{2p} \geq {\Vert \mathbf{\mathsf{B}} \Vert} \geq \min\{m,n\}^{-1/(2p)} {\Vert \mathbf{\mathsf{B}} \Vert}_{2p} \quad\text{for }p \geq 1/2. \end{align*} $$

Thus the Schatten $2p$ -norm is equivalent to the spectral norm when . In fact, when the matrix $\mathbf {\mathsf {B}}$ has a decaying singular value spectrum, the Schatten $2p$ -norm may already be comparable with the spectral norm for much smaller values of p.

Thus we can try to approximate the spectral norm by estimating the Schatten $2p$ -norm for a sufficiently large value of p. Resampling techniques can help ensure that the estimate is reliable. Nevertheless, this method should be used with caution.

6.1. Overview

Consider a PSD matrix $\mathbf {\mathsf {A}} \in \mathbb {H}_n$ with decreasingly ordered eigenvalues $\lambda _1 \geq \lambda _2 \geq \cdots \geq \lambda _n \geq 0$ . We are interested in the problem of estimating the maximum eigenvalue $\lambda _1$ . As in the last two sections, we assume access to $\mathbf {\mathsf {A}}$ via the matrix–vector product $\mathbf {\mathsf {u}} \mapsto \mathbf {\mathsf {A}} \mathbf {\mathsf {u}}$ .

In contrast to the methods in Sections 4 and 5, the algorithms in this section require sequential applications of the matrix–vector product. In other words, we now demand nonlinear information about the input matrix $\mathbf {\mathsf {A}}$ . As a consequence, these algorithms resist parallelization, and they cannot be used in the one-pass or streaming environments.

The theoretical treatment in this section also covers the case of estimating the spectral norm ${\Vert \mathbf {\mathsf {B}} \Vert }$ of a rectangular matrix $\mathbf {\mathsf {B}} \in \mathbb {F}^{m \times n}$ . Indeed, we can simply pass to the PSD matrix $\mathbf {\mathsf {A}} = \mathbf {\mathsf {B}}^{*} \mathbf {\mathsf {B}}$ . From an applied point of view, however, it is important to develop separate algorithms that avoid squaring the matrix $\mathbf {\mathsf {B}}$ . For brevity, we omit a discussion about estimating spectral norms; see Golub and Van Loan (Reference Golub and Van Loan2013, Section 10.4).

6.2. The randomized power method

The randomized power method is a simple iterative algorithm for estimating the maximum eigenvalue.

6.2.1. Procedure

First, draw a random test vector $\mathbf {{\omega }} \in \mathbb {F}^n$ . Since the maximum eigenvalue is unitarily invariant, it is most natural to draw the test vector $\mathbf {{\omega }}$ from a rotationally invariant distribution, such as $\mathbf {{\omega }} \sim {\small\rm NORMAL}(\mathbf {\mathsf {0}}, {\mathbf {\mathsf {I}}})$ . The power method iteratively constructs the sequence

$$ \begin{align*} {\mathbf{\mathsf{y}}}_0 = \displaystyle \frac{\mathbf{{\omega}}}{{\Vert \mathbf{{\omega}} \Vert}} \quad\text{and}\quad {\mathbf{\mathsf{y}}}_q = \displaystyle \frac{\mathbf{\mathsf{A}} {\mathbf{\mathsf{y}}}_{q-1}}{{\Vert {\mathbf{\mathsf{A}} {\mathbf{\mathsf{y}}}_{q-1}} \Vert}} \quad\text{{for }}q \geq 1. \end{align*} $$

At each step, we obtain an eigenvalue estimate

$$ \begin{align*} \xi_q = {\mathbf{\mathsf{y}}}_q^{*} \mathbf{\mathsf{A}} {\mathbf{\mathsf{y}}}_q = \displaystyle \frac{{\mathbf{{\omega}}}^{*} {\mathbf{\mathsf{A}}^{2q+1}} \mathbf{{\omega}}}{{\mathbf{{\omega}}}^{*} {\mathbf{\mathsf{A}}^{2q}} {\mathbf{{\omega}}}} \quad\text{for }q \geq 0. \end{align*} $$

The randomized power method requires simulation of a single random test vector $\mathbf {{\omega }}$ . To perform q iterations, it takes q sequential matrix–vector products with $\mathbf {\mathsf {A}}$ and lower-order arithmetic. It operates with storage ${{O}}(n)$ . See Algorithm 4 for pseudocode.

6.2.2. Analysis

The question is how many iterations q suffice to make $\xi _q$ close to the maximum eigenvalue $\lambda _1$ . More precisely, we aim to control the relative error $e_q$ in the eigenvalue estimate $\xi _q$ :

$$ \begin{align*} e_q = \displaystyle \frac{\lambda_1 - \xi_q}{\lambda_1}. \end{align*} $$

The error $e_q$ is always nonnegative because of the Rayleigh theorem. Note that the computed vector ${\mathbf {\mathsf {y}}}_q$ always has a substantial component in the invariant subspace associated with eigenvalues larger than $\xi _q$ , but it may not be close to any maximum eigenvector, even when $\xi _q \approx \lambda _1$ .

Kuczyński and Woźniakowski (Reference Kuczyński and Woźniakowski1992) have established several remarkable results about the evolution of the error.

Theorem 6.1 (randomized power method)

Let $\mathbf {\mathsf {A}} \in \mathbb {H}_n({\mathbb {R}})$ be a real PSD matrix. Draw a real test vector $\mathbf {{\omega }} \sim {\small\rm NORMAL}(\mathbf {\mathsf {0}}, {\mathbf {\mathsf {I}}})$ . After q iterations of the randomized power method, the error $e_q$ in the maximum eigenvalue estimate $\xi _q$ satisfies

(6.1) $$ \begin{align} \operatorname{\mathbb{E}} e_q \leq 0.871 \cdot \displaystyle \frac{\log n}{q - 1} \quad\text{for }q \geq 2. \end{align} $$

Furthermore, if $\gamma = (\lambda _1 - \lambda _2) / \lambda _1$ is the relative spectral gap, then

(6.2) $$ \begin{align} \operatorname{\mathbb{E}} e_q \leq 1.254 \cdot \sqrt{n} \, \gamma \mathrm{e}^{-q \gamma} \quad\text{for }q \geq 1. \end{align} $$

The second result (6.2) does not appear explicitly in Kuczyński and Woźniakowski (Reference Kuczyński and Woźniakowski1992), but it follows from related ideas (Tropp Reference Tropp2018).

6.3. Randomized Krylov methods

The power method uses the qth power of the matrix to estimate the maximum eigenvalue. A more sophisticated approach allows any polynomial with degree q. Algorithms based on this general technique are often referred to as Krylov subspace methods. The most famous instantiation is the Lanczos method, which is an efficient implementation of a Krylov subspace method for estimating the eigenvalues of a self-adjoint matrix.

6.3.1. Abstract procedure

Draw a random test vector $\mathbf {{\omega }} \in \mathbb {F}^n$ . It is natural to use a rotationally invariant distribution, such as $\mathbf {{\omega }} \sim {\small\rm NORMAL}(\mathbf {\mathsf {0}}, {\mathbf {\mathsf {I}}})$ . For a depth parameter $q \in \mathbb {N}$ , a randomized Krylov subspace method implicitly constructs the subspace

$$ \begin{align*} K_{q+1} := \operatorname{span} \{ \mathbf{{\omega}}, \mathbf{\mathsf{A}}\mathbf{{\omega}}, \ldots, \mathbf{\mathsf{A}}^q \mathbf{{\omega}} \}. \end{align*} $$

We can estimate the maximum eigenvalue of $\mathbf {\mathsf {A}}$ as

$$ \begin{align*} \xi_q = \max_{\mathbf{\mathsf{u}} \in K_{q+1}} \displaystyle \frac{\mathbf{\mathsf{u}}^{*} \mathbf{\mathsf{A}} \mathbf{\mathsf{u}}}{\mathbf{\mathsf{u}}^{*}\mathbf{\mathsf{u}}} = \max_{\deg p \leq q} \displaystyle \frac{{\mathbf{{\omega}}}^{*} \mathbf{\mathsf{A}} p^2(\mathbf{\mathsf{A}}) \mathbf{{\omega}}}{{\Vert { p(\mathbf{\mathsf{A}}) \mathbf{{\omega}} } \Vert}^2}. \end{align*} $$

The maximum occurs over polynomials p with coefficients in $\mathbb {F}$ and with degree at most q. The notation $p(\mathbf {\mathsf {A}})$ refers to the spectral function induced by the polynomial p. We will discuss implementations of Krylov methods below in Section 6.3.4.

6.3.2. Analysis

Since the Krylov subspace is invariant to shifts of the spectrum of $\mathbf {\mathsf {A}}$ , it is more natural to compute the error relative to the spectral range of the matrix:

$$ \begin{align*} f_q = \displaystyle \frac{\lambda_1 - \xi_q}{\lambda_1 - \lambda_n}. \end{align*} $$

The error $f_q$ is always nonnegative because it is a Rayleigh quotient.

Kuczyński and Woźniakowski (Reference Kuczyński and Woźniakowski1992) have established striking results for the maximum eigenvalue estimate obtained via a randomized Krylov subspace method.

Theorem 6.2 (randomized Krylov method)

Let $\mathbf {\mathsf {A}} \in \mathbb {H}_n({\mathbb {R}})$ be a real PSD matrix. Draw a real test vector $\mathbf {{\omega }} \sim {\small\rm NORMAL}(\mathbf {\mathsf {0}}, {\mathbf {\mathsf {I}}})$ . After q iterations of the randomized Krylov method, the error $f_q$ in the maximum eigenvalue estimate $\xi _q$ satisfies

(6.3) $$ \begin{align} \operatorname{\mathbb{E}} f_q \leq 2.575 \cdot \biggl(\displaystyle \frac{\log n}{q - 1} \biggr)^2 \quad\text{for }q \geq 4. \end{align} $$

Furthermore, if $\gamma = (\lambda _1 - \lambda _2) / \lambda _1$ is the relative spectral gap, then

(6.4) $$ \begin{align} \operatorname{\mathbb{E}} f_q \leq 2.589 \cdot \sqrt{n} \, \mathrm{e}^{-2(q-1) \sqrt{\gamma}} \quad\text{for }q \geq 1. \end{align} $$

The second result (6.4) is a direct consequence of a more detailed formula reported in Kuczyński and Woźniakowski (Reference Kuczyński and Woźniakowski1992).

6.3.4. Implementing Krylov methods

What do we have to pay for the superior performance of the randomized Krylov method? If we only need an estimate of the maximum eigenvalue, without an associated eigenvector estimate, the cost is almost the same as for the randomized power method! On the other hand, if we desire the eigenvector estimate, it is common practice to store a basis for the Krylov subspace $K_q$ . This is a classic example of a time–data tradeoff in computation.

We present pseudocode for the randomized Lanczos method, which is an efficient formulation of the Krylov method. Algorithm 5 is a direct implementation of the Lanczos recursion, but it exhibits complicated performance in floating-point arithmetic. Algorithm 5 includes the option to add full reorthogonalization; this step removes the numerical shortcomings at a substantial price in arithmetic (and storage).

If we use the Lanczos method without orthogonalization, then q iterations require q matrix–vector multiplies with $\mathbf {\mathsf {A}}$ plus ${{O}}(qn)$ additional arithmetic. The orthogonalization step adds a total of ${{O}}(q^2 n)$ additional arithmetic. Computing the maximum eigenvalue (and eigenvector) of the tridiagonal matrix $\mathbf {\mathsf {T}}$ can be performed with ${{O}}(q)$ arithmetic (Golub and Van Loan Reference Golub and Van Loan2013, Section 8.4).

If we do not require the maximum eigenvector, the Lanczos method without orthogonalization operates with storage ${{O}}(n)$ . If we need the maximum eigenvector or we add orthogonalization, the storage cost grows to ${{O}}(qn)$ . It is possible to avoid the extra storage by recomputing the Lanczos vectors, but this approach requires great care. One of the main thrusts in the literature on Krylov methods is to reduce these storage costs while maintaining rapid convergence and numerical stability.

Warning 6.3 (Lanczos method)

Algorithm 5 should be used with caution. For a proper discussion about designing Krylov methods, we recommend the books by Parlett (Reference Parlett1998), Bai et al. (Reference Bai, Demmel, Dongarra, Ruhe and van der Vorst1987) and Golub and Van Loan (Reference Golub and Van Loan2013). There is also some recent theoretical work on the numerical stability of Lanczos methods by Musco, Musco and Sidford (Reference Musco, Musco and Sidford2018) and Carmon, Duchi, Sidford and Tian (Reference Carmon, Duchi, Sidford, Tian, Beygelzimer and Hsu2019).

6.4. The minimum eigenvalue

The randomized power method and the randomized Krylov method can be used to estimate the minimum eigenvalue $\lambda _{n}$ of the PSD matrix $\mathbf {\mathsf {A}} \in \mathbb {H}_n$ .

The first approach is to apply the randomized power method to the shifted matrix $\nu {\mathbf {\mathsf {I}}} - \mathbf {\mathsf {A}}$ , where the shift is chosen so that $\nu \geq \lambda _1$ . In this case the algorithm produces an approximation for $\nu - \lambda _n$ . Note that the error in Theorem 6.1 is relative to $\nu - \lambda _n$ , rather than $\lambda _n$ .

The second approach begins with the computation of the Krylov subspace $K_{q+1}$ . Instead of maximizing the Rayleigh quotient over the Krylov subspace, we minimize it:

$$ \begin{align*} \zeta_q = \min_{\mathbf{\mathsf{u}} \in K_{q+1}} \displaystyle \frac{ \mathbf{\mathsf{u}}^{*}\mathbf{\mathsf{A}}\mathbf{\mathsf{u}} }{\mathbf{\mathsf{u}}^{*}\mathbf{\mathsf{u}}}. \end{align*} $$

This approach directly produces an estimate $\zeta _q$ for $\lambda _n$ . The Krylov subspace is invariant under affine transformations of the spectrum of $\mathbf {\mathsf {A}}$ , so we can obtain an error bound for $\zeta _q$ by applying Theorem 6.2 formally to $\lambda _n {\mathbf {\mathsf {I}}} - \mathbf {\mathsf {A}}$ .

6.5. Block methods

The basic power method and Krylov method can be extended by applying the iteration simultaneously to a larger number of (random) test vectors. The resulting algorithms are called subspace iteration and block Krylov methods, respectively. Historically, the reason for developing block methods was to resolve repeated or clustered eigenvalues.

In randomized linear algebra, we discover additional motivations for developing block methods. When the test vectors are drawn at random, block methods may converge slightly faster, and they succeed with much higher probability. On modern computer architectures, the cost of a block method may be comparable with the cost of a simple vector iteration, which makes this modification appealing. We will treat this class of algorithm more thoroughly in Section 11, so we postpone a full discussion. See also Tropp (Reference Tropp2018).

6.6. Estimating trace functions

Finally, we turn to the problem of estimating the trace of a spectral function of a PSD matrix.

6.6.1. Overview

Consider a PSD matrix $\mathbf {\mathsf {A}} \in \mathbb {H}_n$ with eigenpairs $(\lambda _j, \mathbf {\mathsf {u}}_j)$ for $j = 1, \ldots , n$ . Let $f \colon {\mathbb {R}}_+ \to {\mathbb {R}}$ be a function, and suppose that we wish to approximate

$$ \begin{align*} \operatorname{trace} f(\mathbf{\mathsf{A}}) = \sum\nolimits_{j=1}^n f(\lambda_j). \end{align*} $$

We outline an incredible approach to this problem, called stochastic Lanczos quadrature (SLQ), that marries the randomized trace estimator (Section 4) to the Lanczos iteration (Algorithm 5).

This algorithm was devised by Golub and Meurant (Reference Golub and Meurant1994). Our presentation is based on Golub and Van Loan (Reference Golub and Van Loan2013, Section 10.2) and Ubaru, Chen and Saad (Reference Ubaru, Chen and Saad2017). Golub and Meurant (Reference Golub and Meurant2010) provide a more complete treatment, and Musco et al. (Reference Musco, Musco and Sidford2018) give a theoretical discussion about stability.

Related ideas can be used to estimate the trace of a spectral function of a rectangular matrix; that is, the sum of a function of each singular value of the matrix. For brevity, we omit all details on the rectangular case.

6.6.3. Procedure

Let us summarize the mathematical ideas that lead to SLQ. As usual, we draw an isotropic random vector $\mathbf {{\omega }} \in \mathbb {F}^n$ . Then the random variable

$$ \begin{align*} X = {\mathbf{{\omega}}}^{*} f(\mathbf{\mathsf{A}}) \mathbf{{\omega}} \quad\text{satisfies}\quad \operatorname{\mathbb{E}} X = \operatorname{trace} f(\mathbf{\mathsf{A}}). \end{align*} $$

Using the spectral resolution of $\mathbf {\mathsf {A}}$ , we can rewrite X in the form

$$ \begin{align*} X = \sum\nolimits_{j=1}^n f(\lambda_j) \, {\vert {\mathbf{\mathsf{u}}_j^{*} \mathbf{{\omega}}} \vert}^2 = \int_{{\mathbb{R}}_+} f(\lambda) \, {\nu}({\mathrm{d}{\lambda}}) \end{align*} $$

for an appropriate measure $\nu $ on ${\mathbb {R}}_+$ that depends on $\mathbf {\mathsf {A}}$ and $\mathbf {{\omega }}$ . Although we cannot generally compute the integral directly, we can approximate it by using a numerical quadrature rule:

$$ \begin{align*} X \approx \sum\nolimits_{\ell=1}^{q+1} \tau_\ell^2 f(\theta_\ell) =: Z. \end{align*} $$

What is truly amazing is that the weights $\tau _\ell ^2$ and the nodes $\theta _\ell $ for the quadrature rule can be extracted from the tridiagonal matrix $\mathbf {\mathsf {T}} \in \mathbb {H}_{q+1}$ produced by q iterations of the Lanczos iteration with starting vector $\mathbf {{\omega }}$ . This point is not obvious, but a full explanation exceeds our scope.

SLQ approximates the trace function by averaging independent copies of the simple approximation:

$$ \begin{align*} \operatorname{trace} f(\mathbf{\mathsf{A}}) \approx \displaystyle \frac{1}{k} \sum\nolimits_{i=1}^k Z_i \quad\text{where } Z_i \sim Z\text{ are i.i.d.} \end{align*} $$

The analysis of the SLQ approximation requires heavy machinery from approximation theory. See Golub and Meurant (Reference Golub and Meurant2010) and Ubaru et al. (Reference Ubaru, Chen and Saad2017) for more details.

Algorithm 6 contains pseudocode for SLQ. The dominant cost is ${{O}}(kq)$ matrix–vector multiplies with $\mathbf {\mathsf {A}}$ , plus ${{O}}(k q^2)$ additional arithmetic. We recommend using structured random test vectors to reduce the variance of the resulting approximation. The storage cost is ${{O}}(qn)$ numbers.

7.1. Empirical approximation

We begin with a high-level discussion of the method of empirical approximation of a matrix. The examples in this section are all instances of this general idea.

Let $\mathbf {\mathsf {B}} \in \mathbb {F}^{m \times n}$ be a target matrix that we wish to approximate by a more ‘structured’ matrix. Suppose that we can express the matrix $\mathbf {\mathsf {B}}$ as a sum of ‘simple’ matrices:

(7.1) $$ \begin{align} \mathbf{\mathsf{B}} = \sum\nolimits_{i=1}^I {\mathbf{\mathsf{B}}}_i. \end{align} $$

In the cases we will study, the summands $\mathbf {\mathsf {B}}_i$ will be sparse or low-rank.

Next, consider a probability distribution $\{ p_i \colon i = 1, \ldots , I \}$ over the indices in the sum (7.1). For now, we treat this distribution as given. Construct a random matrix $\mathbf {\mathsf {X}} \in \mathbb {F}^{m \times n}$ that takes values

$$ \begin{align*} \mathbf{\mathsf{X}} = p_i^{-1} {\mathbf{\mathsf{B}}}_i \quad\text{with probability }p_i\text{ for each }i = 1, \ldots, I. \end{align*} $$

(Enforce the convention that $0/0 = 0$ .) It is clear that $\mathbf {\mathsf {X}}$ is an unbiased estimator for $\mathbf {\mathsf {B}}$ :

$$ \begin{align*} \operatorname{\mathbb{E}} \mathbf{\mathsf{X}} = \sum\nolimits_{i=1}^I (p_i^{-1} {\mathbf{\mathsf{B}}}_i) p_i = \mathbf{\mathsf{B}}. \end{align*} $$

The random matrix $\mathbf {\mathsf {X}}$ enjoys the same kind of structure as the summands $\mathbf {\mathsf {B}}_i$ . On the other hand, a single draw of the matrix $\mathbf {\mathsf {X}}$ is rarely a good approximation of the matrix $\mathbf {\mathsf {B}}$ .

To obtain a better estimator for $\mathbf {\mathsf {B}}$ , we average independent copies of the initial estimator:

$$ \begin{align*} \bar{\mathbf{\mathsf{X}}}_k = \displaystyle \frac{1}{k} \sum\nolimits_{i=1}^k {\mathbf{\mathsf{X}}}_i \quad\text{where }{\mathbf{\mathsf{X}}}_i \sim \mathbf{\mathsf{X}}\text{ are i.i.d.} \end{align*} $$

By linearity of expectation, $\bar {\mathbf {\mathsf {X}}}_k$ is also unbiased:

$$ \begin{align*} \operatorname{\mathbb{E}} \bar{\mathbf{\mathsf{X}}}_k = \mathbf{\mathsf{B}}. \end{align*} $$

If the parameter k remains small, then $\bar {\mathbf {\mathsf {X}}}_k$ inherits some of the structure from the $\mathbf {\mathsf {B}}_i$ . The question is how many samples k we need to ensure that $\bar {\mathbf {\mathsf {X}}}_k$ also approximates $\mathbf {\mathsf {B}}$ well.

7.2. The matrix Bernstein inequality

The main tool for analysing the sample average estimator $\bar {\mathbf {\mathsf {X}}}_k$ from the last section is a variant of the matrix Bernstein inequality. The following result is drawn from Tropp (Reference Tropp2015).

Theorem 7.2 (matrix Monte Carlo)

Let $\mathbf {\mathsf {B}} \in \mathbb {F}^{m \times n}$ be a fixed matrix. Construct a random matrix $\mathbf {\mathsf {X}} \in \mathbb {F}^{m \times n}$ that satisfies

$$ \begin{align*} \operatorname{\mathbb{E}} \mathbf{\mathsf{X}} = \mathbf{\mathsf{B}} \quad\text{and}\quad {\Vert \mathbf{\mathsf{X}} \Vert} \leq R. \end{align*} $$

Define the per-sample second moment:

$$ \begin{align*} v(\mathbf{\mathsf{X}}) := \max\{ {\Vert { \operatorname{\mathbb{E}}[\mathbf{\mathsf{X}} \mathbf{\mathsf{X}}^{*}] } \Vert}, {\Vert {\operatorname{\mathbb{E}}[\mathbf{\mathsf{X}}^{*}\mathbf{\mathsf{X}}]} \Vert} \}. \end{align*} $$

Form the matrix sampling estimator

$$ \begin{align*} \bar{\mathbf{\mathsf{X}}}_k = \displaystyle \frac{1}{k} \sum\nolimits_{i=1}^k {\mathbf{\mathsf{X}}}_i \quad\text{where }{\mathbf{\mathsf{X}}}_i \sim \mathbf{\mathsf{X}}\text{ are i.i.d.} \end{align*} $$

Then

$$ \begin{align*} \operatorname{\mathbb{E}} {\Vert \bar{\mathbf{\mathsf{X}}}_k - \mathbf{\mathsf{B}} \Vert} \leq \sqrt{\displaystyle \frac{2 v(\mathbf{\mathsf{X}}) \log(m+n)}{k}} + \displaystyle \frac{2R \log(m+n)}{3k}. \end{align*} $$

Furthermore, for all $t \geq 0$ ,

$$ \begin{align*} {{\mathbb{P}\{ { {\Vert \bar{\mathbf{\mathsf{X}}}_k - \mathbf{\mathsf{B}} \Vert} \geq t } \}} \leq (m + n) \exp\biggl( \displaystyle \frac{-kt^2/2}{v(\mathbf{\mathsf{X}}) + 2Rt/3} \biggr).} \end{align*} $$

To explain the meaning of this result, let us determine how large k should be to ensure that the expected approximation error lies below a positive threshold $\varepsilon $ . The bound

$$ \begin{align*} {k \geq \displaystyle \frac{2v(\mathbf{\mathsf{X}})\log(m+n)}{\varepsilon^2} \vee \displaystyle \frac{2R \log(m+n)}{3\varepsilon}} \end{align*} $$

implies that $\operatorname {\mathbb {E}} {\Vert {\bar {\mathbf {\mathsf {X}}}_k - \mathbf {\mathsf {B}} } \Vert } \leq \varepsilon + \varepsilon ^2$ . In other words, the number k of samples should be proportional to the larger of the second moment $v(\mathbf {\mathsf {X}})$ and the upper bound R.

This fact points toward a disappointing feature of empirical matrix approximation: to make $\varepsilon $ small, the number k of samples must increase with $\varepsilon ^{-2}$ and often with $\log (m + n)$ as well. This phenomenon is an unavoidable consequence of the central limit theorem and the geometry induced by the spectral norm. It means that matrix sampling estimators are not suitable for achieving high-precision approximations. See Tropp (Reference Tropp2015, Section 6.2.3) for further discussion. The logarithmic terms are also necessary in the worst case.

7.3. Approximate matrix multiplication

A first application of empirical matrix approximation is to approximate the product of two matrices: $\mathbf {\mathsf {M}} = \mathbf {\mathsf {BC}}$ where $\mathbf {\mathsf {B}} \in \mathbb {F}^{m \times I}$ and $\mathbf {\mathsf {C}} \in \mathbb {F}^{I \times n}$ . Computing this product by direct matrix–matrix multiplication requires $O (mnI)$ arithmetic operations. When the inner dimension I is very large as compared with m and n, we might try to reduce the cost by sampling.

In our own work, we have not encountered situations where approximate matrix multiplication is practical because the quality of the approximation is very low. Nevertheless, the theory serves a dual purpose as the foundation for subspace embedding by discrete sampling (Section 9.6).

7.3.1. Matrix multiplication by sampling

To simplify the analysis, let us pre-scale the matrices $\mathbf {\mathsf {B}}$ and $\mathbf {\mathsf {C}}$ so that each one has spectral norm equal to one:

$$ \begin{align*} {\Vert \mathbf{\mathsf{B}} \Vert} = {\Vert \mathbf{\mathsf{C}} \Vert} = 1. \end{align*} $$

This step can be performed efficiently using the spectral norm estimators outlined in Section 6. This normalization will remain in force for the rest of Section 7.3.

Observe that the matrix–matrix product satisfies

$$ \begin{align*} \mathbf{\mathsf{BC}} = \sum\nolimits_{i=1}^I (\mathbf{\mathsf{B}})_{:i} (\mathbf{\mathsf{C}})_{i:}, \end{align*} $$

where $(\mathbf {\mathsf {B}})_{:i}$ is the ith column of $\mathbf {\mathsf {B}}$ and $(\mathbf {\mathsf {C}})_{i:}$ is the ith row of $\mathbf {\mathsf {C}}$ . This expression motivates us to approximate the product by sampling terms at random from the sum.

Let $\{ p_i \colon i = 1, \ldots , I \}$ be a sampling distribution, to be specified later. Form a random rank-one unbiased estimator for the product:

$$ \begin{align*} \mathbf{\mathsf{X}} = p_i^{-1} \cdot (\mathbf{\mathsf{B}})_{:i} (\mathbf{\mathsf{C}})_{i:} \quad\text{with probability }p_i\text{ for each }i = 1, \ldots, I. \end{align*} $$

By construction, $\operatorname {\mathbb {E}} \mathbf {\mathsf {X}} = \mathbf {\mathsf {BC}}$ . We can average k independent copies of $\mathbf {\mathsf {X}}$ to obtain a better approximation $\bar {\mathbf {\mathsf {X}}}_k$ to the matrix product. The cost of computing the estimator $\bar {\mathbf {\mathsf {X}}}_k$ of the matrix product explicitly is only $O (mnk)$ operations, so it is more efficient than the full multiplication when $k \ll I$ .

Theorem 7.2 yields an easy analysis of this approach. The conclusion depends heavily on the choice of sampling distribution. Regardless, we can never expect to attain a high-accuracy approximation of the product by sampling because the number k of samples must scale proportionally with the inverse square $\varepsilon ^{-2}$ of the accuracy parameter $\varepsilon $ .

7.3.2. Uniform sampling

The easiest way to approximate matrix multiplication is to choose uniform sampling probabilities: $p_i = 1/I$ for each $i = 1, \ldots , I$ . To analyse this case, we introduce the coherence parameter:

$$ \begin{align*} \mu(\mathbf{\mathsf{B}}) := I \cdot \max_{i=1, \ldots, I} {\Vert {(\mathbf{\mathsf{B}})_{:i}} \Vert}^2. \end{align*} $$

Up to scaling, this is the maximum squared norm of a column of $\mathbf {\mathsf {B}}$ . Since $\mathbf {\mathsf {B}} \in \mathbb {F}^{m \times I}$ and ${\Vert \mathbf {\mathsf {B}} \Vert }=1$ , the coherence lies in the range $[m, I]$ . The difficulty of approximating matrix multiplication by uniform sampling increases with the coherence of $\mathbf {\mathsf {B}}$ and $\mathbf {\mathsf {C}}^{*}$ .

To invoke Theorem 7.2, observe that the per-sample second moment and the spectral norm of the estimator $\mathbf {\mathsf {X}}$ satisfy the bound

$$ \begin{align*} \max\{ v(\mathbf{\mathsf{X}}), {\Vert \mathbf{\mathsf{X}} \Vert} \} \leq \max\{\mu(\mathbf{\mathsf{B}}), \mu(\mathbf{\mathsf{C}}^{*})\}. \end{align*} $$

Let $\varepsilon \in (0, 1]$ be an accuracy parameter. If

$$ \begin{align*} k \geq 2 \varepsilon^{-2} \max\{ \mu(\mathbf{\mathsf{B}}), \mu(\mathbf{\mathsf{C}}^{*}) \} \log(m + n), \end{align*} $$

then

$$ \begin{align*} \operatorname{\mathbb{E}} {\Vert { \bar{\mathbf{\mathsf{X}}}_k - \mathbf{\mathsf{BC}} } \Vert} \leq 2\varepsilon. \end{align*} $$

In other words, the number k of rank-one factors we need to obtain a relative approximation of the matrix product is proportional to the maximum coherence of $\mathbf {\mathsf {B}}$ and $\mathbf {\mathsf {C}}^{*}$ . In the best scenario, the number of samples is proportional to $\max \{m,n\} \log (m+n)$ ; in the worst case, the sample complexity can be as large as I.

7.3.3. Importance sampling

If the norms of the columns of the factors $\mathbf {\mathsf {B}}$ and $\mathbf {\mathsf {C}}^{*}$ vary wildly, we may need to use importance sampling to obtain a nontrivial approximation bound. Define the sampling probabilities

$$ \begin{align*} p_i = \displaystyle \frac{{\Vert {(\mathbf{\mathsf{B}})_{:i}} \Vert}^2 + {\Vert {(\mathbf{\mathsf{C}})_{:i}} \Vert}^2}{{{\Vert {{\mathbf{\mathsf{B}}}} \Vert}^2_{\mathrm{F}}} + {{\Vert {{\mathbf{\mathsf{C}}}} \Vert}^2_{\mathrm{F}}}} \quad\text{for }i = 1, \ldots, I. \end{align*} $$

These probabilities are designed to balance terms arising from Theorem 7.2. In most cases we compute the sampling distribution by directly evaluating the formula, at the cost of $O ((m + n)I)$ operations.

With the importance sampling distribution, the per-sample second moment and the spectral norm of the estimator $\mathbf {\mathsf {X}}$ satisfy

$$ \begin{align*} \max\{ v(\mathbf{\mathsf{X}}), {\Vert \mathbf{\mathsf{X}} \Vert} \} \leq \displaystyle \frac{1}{2} (\operatorname{srank}(\mathbf{\mathsf{B}}) + \operatorname{srank}(\mathbf{\mathsf{C}})). \end{align*} $$

The stable rank is defined in (2.2). Let $\varepsilon \in (0, 1]$ be an accuracy parameter. If

$$ \begin{align*} k \geq \varepsilon^{-2} ( \operatorname{srank}(\mathbf{\mathsf{B}}) + \operatorname{srank}(\mathbf{\mathsf{C}}) ) \log(m + n), \end{align*} $$

then

$$ \begin{align*} \operatorname{\mathbb{E}} {\Vert { \bar{\mathbf{\mathsf{X}}}_k - \mathbf{\mathsf{BC}} } \Vert} \leq 2 \varepsilon. \end{align*} $$

In other words, the number of rank-one factors we need to approximate the matrix product by importance sampling is $\log (m+n)$ times the total stable rank of the matrices.

The sample complexity bound for importance sampling always improves over the bound for uniform sampling. Indeed, $\operatorname {srank}(\mathbf {\mathsf {B}}) \leq \mu (\mathbf {\mathsf {B}})$ under the normalization ${\Vert \mathbf {\mathsf {B}} \Vert }=1$ . There are also many situations where the stable rank is smaller than either dimension of the matrix. These are the cases where one might consider using approximate matrix multiplication by importance sampling.

7.4. Approximating a graph by a sparse graph

As a second application of empirical matrix approximation, we will show how to take a dense graph and find a sparse graph that serves as a proxy. This procedure operates by replacing the Laplacian matrix of the graph with a sparse Laplacian matrix. Beyond its intrinsic interest as a fact about graphs, this technique plays a central role in randomized solvers for Laplacian linear systems (Section 18).

7.4.1. Graphs and Laplacians

We will consider weighted, loop-free, undirected graphs on the vertex set $V = \{1, \ldots , n\}$ . We can specify a weighted graph G on V by means of a nonnegative weight function $w \colon V \times V \to {\mathbb {R}}_+$ . The graph is loop-free when $w_{ii} = 0$ for each vertex $i \in V$ . The graph is undirected if and only if $w_{ij} = w_{ji}$ for each pair $(i, j)$ . The sparsity of an undirected graph is the number of strictly positive weights $w_{ij}$ with $i \leq j$ .

Alternatively, we can work with graph Laplacians. The elementary Laplacian $\mathbf {\mathsf{\Delta }}_{ij}$ on a vertex pair $(i, j) \in V \times V$ is the PSD matrix

$$ \begin{align*} \mathbf{\mathsf{\Delta}}_{ij} = (\mathbf{{\delta}}_{i} - \mathbf{{\delta}}_j)(\mathbf{{\delta}}_i - \mathbf{{\delta}}_j)^{*} \in \mathbb{H}_n({\mathbb{R}}). \end{align*} $$

Let w be the weight function of a loop-free, undirected graph G. The Laplacian associated with the graph G is the matrix

(7.2) $$ \begin{align} {\mathbf{\mathsf{L}}}_{G} = \sum\nolimits_{1 \leq i < j \leq n} {w_{ij}} {\mathbf{\mathsf{\Delta}}_{ij}} \in \mathbb{H}_n({\mathbb{R}}). \end{align} $$

The Laplacian $\mathbf {\mathsf {L}}_G$ is PSD because it is a nonnegative linear combination of PSD matrices.

The Laplacian of a graph is analogous to the Laplacian differential operator. You can think about $\mathbf {\mathsf {L}}_G$ as an analogue of the heat kernel that models the diffusion of a particle on the graph. The Poisson problem $\mathbf {\mathsf {L}}_G \mathbf {\mathsf {x}} = \mathbf {\mathsf {f}}$ serves as a primitive for answering a wide range of questions involving undirected graphs (Teng Reference Teng2010). Applications include clustering and partitioning data, studying random walks on graphs, and solving finite-element discretizations of elliptic PDEs.

7.4.2. Spectral approximation

Fix a parameter $\varepsilon \in (0, 1)$ . We say that a graph H is an $\varepsilon $ -spectral approximation of a graph G if their Laplacians are comparable in the PSD order:

(7.3)

The relation (7.3) ensures that the graph H and the graph G are close cousins.

In particular, under (7.3), the matrix $\mathbf {\mathsf {L}}_H$ serves as an excellent preconditioner for the Laplacian $\mathbf {\mathsf {L}}_G$ . In other words, if we can easily solve (consistent) linear systems of the form $\mathbf {\mathsf {L}}_H \mathbf {\mathsf {y}} = \mathbf {\mathsf {b}}$ , then we can just as easily solve the Poisson problem $\mathbf {\mathsf {L}}_G \mathbf {\mathsf {x}} = \mathbf {\mathsf {f}}$ .

For an arbitrary input graph G, we will demonstrate that there is a sparse graph H that is a good spectral approximation of G. For several reasons, this construction does not immediately lead to effective methods for designing preconditioners. Nevertheless, related ideas have resulted in practical, fast solvers for Poisson problems on undirected graphs. We will make this connection in Section 18.

7.4.3. The normalizing map

It is convenient to present a few more concepts from spectral graph theory. Let us introduce the normalizing map

$$ \begin{align*} \mathbf{\mathsf{K}}_G(\mathbf{\mathsf{A}}) := (\mathbf{\mathsf{L}}_G^{\dagger})^{1/2} \mathbf{\mathsf{A}} (\mathbf{\mathsf{L}}_G^{\dagger})^{1/2} \quad\text{for }\mathbf{\mathsf{A}} \in \mathbb{H}_n({\mathbb{R}}). \end{align*} $$

As usual, $(\cdot )^{\dagger }$ is the pseudoinverse and $(\cdot )^{1/2}$ is the PSD square root of a PSD matrix. We use the normalizing map to compare graph Laplacians. Indeed,

(7.4) $$ \begin{align} {\Vert { \mathbf{\mathsf{K}}_G(\mathbf{\mathsf{L}}_H - \mathbf{\mathsf{L}}_G) } \Vert} \leq \varepsilon \end{align} $$

implies that the graph H is an $\varepsilon $ -spectral approximation of the graph G. This claim follows easily from (7.3).

7.4.4. Effective resistance

Next, define the effective resistance $\varrho _{ij}$ of a vertex pair $(i, j)$ in the graph G as

(7.5) $$ \begin{align} \varrho_{ij} := \operatorname{trace}[ \mathbf{\mathsf{K}}_G(\mathbf{\mathsf{\Delta}}_{ij}) ] \geq 0 \quad\text{for each }1 \leq i < j \leq n. \end{align} $$

We can compute the family of effective resistances in time $O (n^3)$ by means of a Cholesky factorization of $\mathbf {\mathsf {L}}_G$ , although faster algorithms are now available (Kyng Reference Kyng2017).

To understand the terminology, let us regard the graph G as an electrical network where $w_{ij}$ is the conductivity of the wire connecting the vertex pair $(i, j)$ . The effective resistance $\varrho _{ij}$ is the resistance of the entire electrical network G against passing a unit of current from vertex i to vertex j.

7.4.5. Sparsification by sampling

We are now prepared to construct a sparse approximation of a loop-free, undirected graph G specified by the weight function w. The representation (7.2) of the graph Laplacian immediately suggests that we can apply the empirical approximation paradigm. To do so, we must design an appropriate sampling distribution.

Define the sampling probabilities

$$ \begin{align*} p_{ij} = \displaystyle \frac{w_{ij} \, \varrho_{ij}}{\operatorname{rank}(\mathbf{\mathsf{L}}_G)} \quad\text{for each }1 \leq i < j \leq n. \end{align*} $$

It is clear that $p_{ij} \geq 0$ . With some basic matrix algebra, we can also confirm that $\sum \nolimits _{i < j} p_{ij} = 1$ .

Following our standard approach, we construct the random matrix

$$ \begin{align*} \mathbf{\mathsf{X}} = \displaystyle \frac{w_{ij}}{p_{ij}} \, \mathbf{\mathsf{\Delta}}_{ij} \quad\text{with probability }p_{ij}\text{ for each }i < j. \end{align*} $$

Next, we average k independent copies of the estimator:

$$ \begin{align*} \bar{\mathbf{\mathsf{X}}}_k = \displaystyle \frac{1}{k} \sum\nolimits_{i=1}^k {\mathbf{\mathsf{X}}}_i \quad\text{where }{\mathbf{\mathsf{X}}}_i \sim \mathbf{\mathsf{X}}\text{ are i.i.d.} \end{align*} $$

By construction, the estimator is unbiased: $\operatorname {\mathbb {E}} \bar {\mathbf {\mathsf {X}}}_k = \mathbf {\mathsf {L}}_G$ . Moreover, $\bar {\mathbf {\mathsf {X}}}_k$ is itself the graph Laplacian of a (random) graph H with sparsity at most k. We just need to determine the number k of samples that are sufficient to make H a good spectral approximation of G.

Given the effective resistances, computation of the sampling probabilities requires $O (s)$ time, where s is the sparsity of the graph. The cost of sampling k copies of $\mathbf {\mathsf {X}}$ is $\tilde {O}(s + k)$ . It is natural to represent $\bar {\mathbf {\mathsf {X}}}_k$ using a sparse matrix data structure, with storage cost $O (k \log n)$ .

7.4.6. Analysis

The analysis involves a small twist. In view of (7.4), we need to demonstrate that $\mathbf {\mathsf {K}}_G(\mathbf {\mathsf {L}}_H) \approx \mathbf {\mathsf {K}}_G(\mathbf {\mathsf {L}}_G)$ . Therefore, instead of considering the random matrix $\bar {\mathbf {\mathsf {X}}}_k$ , we pass to the random matrix $\mathbf {\mathsf {K}}_G(\bar {\mathbf {\mathsf {X}}}_k)$ , which is an unbiased estimator for $\mathbf {\mathsf {K}}_G(\mathbf {\mathsf {L}}_G)$ .

Note that the random matrix $\mathbf {\mathsf {K}}_G(\mathbf {\mathsf {X}})$ satisfies

$$ \begin{align*} \max\{ v( \mathbf{\mathsf{K}}_G(\mathbf{\mathsf{X}}) ), {\Vert { \mathbf{\mathsf{K}}_G(\mathbf{\mathsf{X}}) } \Vert}_2 \} \leq \operatorname{rank}(\mathbf{\mathsf{L}}_G) < n. \end{align*} $$

Suppose that we choose

$$ \begin{align*} k \geq 3 \varepsilon^{-2} n \log(2n) \quad\text{where }\varepsilon \in (0, 1). \end{align*} $$

Then Theorem 7.2 implies

$$ \begin{align*} \operatorname{\mathbb{E}} {\Vert { \mathbf{\mathsf{K}}_G( \bar{\mathbf{\mathsf{X}}}_k - \mathbf{\mathsf{L}}_G ) } \Vert}_2 \leq \varepsilon. \end{align*} $$

We conclude that the random graph H with Laplacian $\mathbf {\mathsf {L}}_H = \bar {\mathbf {\mathsf {X}}}_k$ has at most k nonzero weights, and it is an $\varepsilon $ -spectral approximation to the graph G.

In other words, every graph on n vertices has a $(1/2)$ -spectral approximation with at most $12 n \log n$ nonzero weights. Modulo the precise constant, this is the tightest result that can be obtained if we form the graph H via random sampling.

8.1. What is a random embedding?

Let $E \subseteq \mathbb {F}^n$ be a set, and let $\varepsilon \in (0, 1)$ be a distortion parameter. We say that a linear map $\mathbf {\mathsf {S}} : \mathbb {F}^{n} \to \mathbb {F}^d$ is an ( $\ell _2$ ) embedding of E with distortion $\varepsilon $ when

(8.1) $$ \begin{align} (1 - \varepsilon) \, {\Vert \mathbf{\mathsf{x}} \Vert} \leq {\Vert { \mathbf{\mathsf{S}} \mathbf{\mathsf{x}} } \Vert} \leq (1 + \varepsilon) \, {\Vert \mathbf{\mathsf{x}} \Vert} \quad\text{for all }\mathbf{\mathsf{x}} \in E. \end{align} $$

It is sometimes convenient to abbreviate this kind of two-sided inequality as ${\Vert \mathbf {\mathsf {S}}\mathbf {\mathsf {x}} \Vert } = (1 \pm \varepsilon ) {\Vert \mathbf {\mathsf {x}} \Vert }$ .

We usually think about the case where $d \ll n$ , so the map $\mathbf {\mathsf {S}}$ enacts a dimension reduction. In other words, $\mathbf {\mathsf {S}}$ transfers data from the high-dimensional space $\mathbb {F}^n$ to the low-dimensional space $\mathbb {F}^d$ . As we will discuss, the low-dimensional representation of the data can be used to obtain fast, approximate solutions to computational problems.

The relation (8.1) expresses the idea that the embedding $\mathbf {\mathsf {S}}$ should preserve the geometry of the set E. Unfortunately, we do not always know the set E in advance. Moreover, we would like the map $\mathbf {\mathsf {S}}$ to be easy to construct, and it should be computationally efficient to apply $\mathbf {\mathsf {S}}$ to the data. These goals may be in tension.

We can resolve this dilemma by drawing the embedding $\mathbf {\mathsf {S}}$ from a probability distribution. Many types of probability distributions serve. In particular, we can use highly structured random matrices (Section 9) that are easy to build, to store, and to apply to vectors. Section 10 presents a case study about how random embeddings can be applied to solve overdetermined least-squares problems.

8.2. Restricted singular values

Our initial goal is to understand something about the theoretical behaviour of randomized embeddings. To that end, let us introduce quantities that measure how much an embedding distorts a set. Let $\mathbf {\mathsf {S}} \in \mathbb {F}^{d \times n}$ be a linear map, and let $E \subseteq \mathbb {S}^{n-1}(\mathbb {F})$ be an arbitrary subset of the unit sphere in $\mathbb {F}^n$ . The minimum and maximum restricted singular value are, respectively, defined as

(8.2) $$ \begin{align} \sigma_{\min}(\mathbf{\mathsf{S}}; E) := \min_{\mathbf{\mathsf{x}} \in E}\ {\Vert { \mathbf{\mathsf{S}} \mathbf{\mathsf{x}} } \Vert} \quad\text{and}\quad \sigma_{\max}(\mathbf{\mathsf{S}}; E) := \max_{\mathbf{\mathsf{x}} \in E}\ {\Vert { \mathbf{\mathsf{S}} \mathbf{\mathsf{x}} } \Vert}. \end{align} $$

If E composes the entire unit sphere, then these quantities coincide with the ordinary minimum and maximum singular value of $\mathbf {\mathsf {S}}$ . More generally, the restricted singular values describe how much the linear map $\mathbf {\mathsf {S}}$ can contract or expand a point in E.

8.3. Gaussian embeddings

Our theoretical treatment of random embeddings focuses on the most highly structured case. A Gaussian embedding is a random matrix of the form

$$ \begin{align*} \mathbf{\mathsf{\Gamma}} \in \mathbb{F}^{d \times n} \quad\text{with i.i.d. entries }(\mathbf{\mathsf{\Gamma}})_{ij} \sim {\small\rm NORMAL}(0, d^{-1}). \end{align*} $$

The cost of explicitly storing a Gaussian embedding is ${{O}}(dn)$ , and the cost of applying it to a vector is ${{O}}(dn)$ .

The scaling of the matrix ensures that

$$ \begin{align*} \operatorname{\mathbb{E}} {\Vert { \mathbf{\mathsf{\Gamma}} \mathbf{\mathsf{x}} } \Vert}^2 = {\Vert \mathbf{\mathsf{x}} \Vert}^2 \quad\text{for each }\mathbf{\mathsf{x}} \in \mathbb{F}^n. \end{align*} $$

We wish to understand how large to choose the embedding dimension d so that the map $\mathbf {\mathsf{\Gamma }}$ approximately preserves the norms of all points in a given set E. We can do so by obtaining bounds for the restricted singular values. For a Gaussian embedding, we will see that $\sigma _{\min }(\mathbf {\mathsf{\Gamma }}; E)$ and $\sigma _{\max }(\mathbf {\mathsf{\Gamma }}; E)$ are controlled by the geometry of the set E.

8.4. The Gaussian width

For the remainder of Section 8 we will work in the real field ( $\mathbb {F} = {\mathbb {R}}$ ). Given a set $E \subseteq \mathbb {S}^{n-1}({\mathbb {R}})$ , define the Gaussian width $w(E)$ via

$$ \begin{align*} w(E) := \operatorname{\mathbb{E}} \sup_{\mathbf{\mathsf{x}} \in E}\ {\langle { \mathbf{\mathsf{g}} }, \, { \mathbf{\mathsf{x}} } \rangle} \quad\text{where }\mathbf{\mathsf{g}} \in {\mathbb{R}}^n\text{ is standard normal.} \end{align*} $$

The Gaussian width is a measure of the content of the set E. It plays a fundamental role in the performance of randomized embeddings.

Here are some basic properties of the Gaussian width.

• • The width is invariant under rotations: $w(\mathbf {\mathsf {Q}}E) = w(E)$ for each orthogonal matrix $\mathbf {\mathsf {Q}}$ .

• • The width is increasing with respect to set inclusion: $E \subseteq F$ implies that $w(E) \leq w(F)$ .

• • The width lies in the range $0 \leq w(E) \leq \operatorname {\mathbb {E}} {\Vert \mathbf {\mathsf {g}} \Vert } < \sqrt {n}$ .

The width can be calculated accurately for many sets of interest. In particular, if L is an arbitrary k-dimensional subspace of ${\mathbb {R}}^n$ , then

(8.3) $$ \begin{align} \sqrt{k-1} < w( L \cap \mathbb{S}^{n-1} ) < \sqrt{k}. \end{align} $$

Indeed, it is productive to think about the squared width $w^2(E)$ as a measure of the ‘dimension’ of the set E.

8.5. Restricted singular values of Gaussian matrices

In a classic work on Banach space geometry, Gordon (Reference Gordon1988) showed that the Gaussian width controls both the minimum and maximum restricted singular values of a subset of the sphere.

Theorem 8.4 (restricted singular values: Gaussian matrix)

Fix a subset $E \subseteq \mathbb {S}^{n-1}({\mathbb {R}})$ of the unit sphere. Draw a Gaussian matrix $\mathbf {\mathsf{\Gamma }} \in {\mathbb {R}}^{d \times n}$ whose entries are i.i.d. ${\small\rm NORMAL}(0, d^{-1})$ . For all $t> 0$ ,

$$ \begin{align*} {\mathbb{P}\biggl\{ { \sigma_{\min}(\mathbf{\mathsf{\Gamma}}; E ) \leq 1 - \displaystyle \frac{w(E) + 1}{\sqrt{d}} - t } \biggr\}} &\leq \mathrm{e}^{-dt^2/2}{,} \\* {\mathbb{P}\biggl\{ { \sigma_{\max}(\mathbf{\mathsf{\Gamma}}; E) \geq 1 + \displaystyle \frac{w(E)}{\sqrt{d}} + t } \biggr\}} &\leq \mathrm{e}^{-dt^2/2}. \end{align*} $$

Theorem 8.4 yields the relations

In other words, the embedding dimension should satisfy $d> (w(E) + 1)^2$ to ensure that the map $\mathbf {\mathsf{\Gamma }}$ is unlikely to annihilate any point in E. For this choice of d, the random embedding is unlikely to dilate any point in E by more than a factor of two.

As a consequence, we have reduced the problem of computing embedding dimensions for Gaussian maps to the problem of computing Gaussian widths. In the next two subsections, we work out two important examples.

Remark 8.5 (optimality)

The statements in Theorem 8.4 are nearly optimal. One way to see this is to consider the set $E = \mathbb {S}^{n-1}({\mathbb {R}})$ , for which the theorem implies that

The Bai–Yin law (Bai and Silverstein Reference Bai and Silverstein2010, Section 5.2) confirms that the first and last inequality are sharp as $n, d \to \infty $ with $n/d \to \mathrm {const.} \in [0,1]$ .

Moreover, if E is spherically convex (i.e. the intersection of a convex cone with the unit sphere), then the minimum restricted singular value satisfies the reverse inequality

$$ \begin{align*} {\mathbb{P}\biggl\{ { \sigma_{\min}(\mathbf{\mathsf{\Gamma}}; E) \geq 1 - \displaystyle \frac{w(E)}{\sqrt{d}} + t } \biggr\}} \leq 2\mathrm{e}^{-dt^2/2}. \end{align*} $$

This result is adapted from Thrampoulidis et al. (Reference Thrampoulidis, Oymak and Hassibi2014).

In addition, fifteen years of computational experiments have also shown that the predictions from Theorem 8.4 are frequently sharp. See Oymak and Tropp (Reference Oymak and Tropp2018) for some examples and references.

8.6. Example: Johnson–Lindenstrauss

As a first application of Theorem 8.4, let us explain how it implies the classic dimension reduction result of Johnson and Lindenstrauss (Reference Johnson and Lindenstrauss1984).

8.6.1. Overview

Let $\{ \mathbf {\mathsf {a}}_1, \ldots , \mathbf {\mathsf {a}}_N \} \subset {\mathbb {R}}^n$ be a discrete point set. We would like to know when a Gaussian embedding $\mathbf {\mathsf{\Gamma }} \in {\mathbb {R}}^{d \times n}$ approximately preserves all the pairwise distances between these points:

(8.4) $$ \begin{align} 1 - \varepsilon \leq \displaystyle \frac{{\Vert { \mathbf{\mathsf{\Gamma}}(\mathbf{\mathsf{a}}_i - \mathbf{\mathsf{a}}_j) } \Vert}}{ {\Vert { \mathbf{\mathsf{a}}_i - \mathbf{\mathsf{a}}_j } \Vert} } \leq 1 + \varepsilon \quad\text{for all }i \neq j. \end{align} $$

The question is how large we must set the embedding dimension d to achieve distortion $\varepsilon \in (0, 1)$ .

8.6.2. Analysis

We can solve this problem using the machinery described in Section 8.5. Consider the set E of normalized chords:

$$ \begin{align*} E = \biggl\{ \displaystyle \frac{\mathbf{\mathsf{a}}_i - \mathbf{\mathsf{a}}_j}{{\Vert {\mathbf{\mathsf{a}}_i - \mathbf{\mathsf{a}}_j} \Vert}} : 1 \leq i < j \leq N \biggr\}. \end{align*} $$

By the definition (8.2) of the restricted singular values,

$$ \begin{align*} \sigma_{\min}(\mathbf{\mathsf{\Gamma}}; E) \leq \displaystyle \frac{{\Vert { \mathbf{\mathsf{\Gamma}}(\mathbf{\mathsf{a}}_i - \mathbf{\mathsf{a}}_j) } \Vert}}{ {\Vert { \mathbf{\mathsf{a}}_i - \mathbf{\mathsf{a}}_j } \Vert} } \leq \sigma_{\max}(\mathbf{\mathsf{\Gamma}}; E). \end{align*} $$

Therefore we can invoke Theorem 8.4 to determine how the embedding dimension controls the distortion.

Let us summarize the argument. First, observe that the Gaussian width of the set E satisfies

$$ \begin{align*} w(E) = \operatorname{\mathbb{E}} \max_{\mathbf{\mathsf{x}} \in E}\ {\langle { \mathbf{\mathsf{g}} }, \, { \mathbf{\mathsf{x}} } \rangle} \leq \sqrt{2 \log \# E} < 2 \sqrt{\log (N/2)}. \end{align*} $$

As a consequence,

$$ \begin{align*} {\mathbb{P}\{ { \sigma_{\min}(\mathbf{\mathsf{\Gamma}}; E) \leq 1 - (1 + 2 \sqrt{\log(N/2)})/\sqrt{d} - t } \}} &\leq \mathrm{e}^{-d t^2 /2}{,} \\* {\mathbb{P}\{ { \sigma_{\max}(\mathbf{\mathsf{\Gamma}}; E) \geq 1 + 2 \sqrt{\log(N/2)}/\sqrt{d} + t } \}} &\leq \mathrm{e}^{-d t^2 /2}. \end{align*} $$

To achieve distortion $\varepsilon $ with high probability, it is sufficient to choose

$$ \begin{align*} d \geq 8 \varepsilon^{-2} \log N. \end{align*} $$

In other words, the embedding dimension only needs to be logarithmic in the cardinality N of the point set. With some additional calculation, we can also extract precise failure probabilities from this analysis.

8.7. Example: Subspace embedding

Next, we consider a question at the heart of randomized linear algebra. Can we embed an unknown subspace into a lower-dimensional space?

8.7.1. Overview

Suppose that L is a k-dimensional subspace in ${\mathbb {R}}^n$ . We say that a dimension reduction map $\mathbf {\mathsf {S}}$ is a subspace embedding for L with distortion $\varepsilon \in (0, 1)$ if

(8.5) $$ \begin{align} (1 - \varepsilon) \, {\Vert \mathbf{\mathsf{x}} \Vert} \leq {\Vert { \mathbf{\mathsf{S}}\mathbf{\mathsf{x}} } \Vert} \leq (1 + \varepsilon) \, {\Vert \mathbf{\mathsf{x}} \Vert} \quad\text{for every }\mathbf{\mathsf{x}} \in L. \end{align} $$

We say that $\mathbf {\mathsf {S}}$ is oblivious if it can be constructed without knowledge of the subspace L, except for its dimension.

Two questions arise. First, what types of dimension reduction maps yield (oblivious) subspace embeddings? Second, how large must we choose the embedding dimension to achieve this outcome?

8.7.2. Analysis

Gaussian dimension reduction maps yield very good oblivious subspace embeddings. Theorem 8.4 easily furnishes the justification. Consider the unit sphere in the subspace: $E = L \cap \mathbb {S}^{n-1}({\mathbb {R}})$ . Then construct the Gaussian dimension reduction map $\mathbf {\mathsf{\Gamma }} \in {\mathbb {R}}^{d \times n}$ . In view of (8.3), we have

$$ \begin{align*} {\mathbb{P}\{ { \sigma_{\min}( \mathbf{\mathsf{\Gamma}}; E ) \leq 1 - (1+\sqrt{k})/\sqrt{d} - t } \}} &\leq \mathrm{e}^{-d t^2/2}{,} \\* {\mathbb{P}\{ { \sigma_{\max}( \mathbf{\mathsf{\Gamma}}; E ) \geq 1 + \sqrt{k}/\sqrt{d} + t } \}} &\leq \mathrm{e}^{-d t^2/2}. \end{align*} $$

As a specific example, we can set the embedding dimension $d = 2k$ to ensure that ${\Vert \mathbf {\mathsf{\Gamma }}\mathbf {\mathsf {x}} \Vert } = (1 \pm 0.8) {\Vert \mathbf {\mathsf {x}} \Vert }$ simultaneously for all points $\mathbf {\mathsf {x}} \in L$ , except with probability $\mathrm {e}^{-\mathrm {c} k}$ . In some applications of subspace embeddings, we can even choose the dimension as small as $d = k + 5$ or $d = k + 10$ .

Many theoretical papers on randomized NLA use subspace embeddings as a primitive for designing algorithms for other linear algebra problems. For example, Section 10 describes several ways to use subspace embeddings to solve overdetermined least-squares problems.

8.8. Universality of the minimum restricted singular value

We have seen how to apply Gaussian dimension reduction for embedding discrete point sets and for embedding subspaces. Theorem 8.4 contains precise theoretical results on the behaviour of Gaussian maps in terms of the Gaussian width. To what extent can we transfer this analysis to other types of random embeddings?

The following theorem (Oymak and Tropp Reference Oymak and Tropp2018, Theorem 9.1) shows that the bound on the minimum restricted singular value in Theorem 8.4 is universal for a large class of random embeddings. In particular, this class includes sparse random matrices, whose nonzero entries compose a vanishing proportion of the total.

Theorem 8.7 (universality)

Fix a set $E \subseteq \mathbb {S}^{n-1}$ . Let $\mathbf {\mathsf {S}} \in {\mathbb {R}}^{d \times n}$ be a random matrix whose entries are independent random variables that satisfy

$$ \begin{align*} \operatorname{\mathbb{E}}[ (\mathbf{\mathsf{S}})_{ij} ] = 0, \quad \operatorname{\mathbb{E}}[ (\mathbf{\mathsf{S}})_{ij}^2 ] = d^{-1}, \quad \operatorname{\mathbb{E}}[ (\mathbf{\mathsf{S}})_{ij}^5 ] \leq R. \end{align*} $$

When $d \leq n$ , with high probability,

$$ \begin{align*} \sigma_{\min}(\mathbf{\mathsf{S}}; E) \geq 1 - \displaystyle \frac{w(E)}{\sqrt{d}} - o(\sqrt{n/d}). \end{align*} $$

The constant in $o(\sqrt {n/d})$ depends only on R. A matching lower bound for $\sigma _{\min }(\mathbf {\mathsf {S}}; E)$ holds when E is spherically convex.

In other words, if E is a moderately large set, the distribution of the entries of the random map $\mathbf {\mathsf {S}}$ does not have an impact on the embedding dimension d sufficient to ensure no point in E is annihilated.

Theorem 8.7 is confirmed by extensive numerical experiments (Oymak and Tropp Reference Oymak and Tropp2018), which demonstrate that dimension reduction maps with independent, standardized entries have identical performance for a wide range of examples.

It is perhaps surprising that the bound on the maximum restricted singular value from Theorem 8.4 is not universal. For some sets E, the quantity $\sigma _{\max }(\mathbf {\mathsf {S}}; E)$ depends heavily on the distribution of the entries of $\mathbf {\mathsf {S}}$ .

8.9. Random partial isometries

Last, we consider a variant of Gaussian embedding that is more suitable when the embedding dimension d is close to the ambient dimension n. In this section we allow the field $\mathbb {F}$ to be real or complex.

First, suppose that $d \leq n$ , and let $\mathbf {\mathsf{\Gamma }} \in \mathbb {F}^{d \times n}$ be a Gaussian embedding. Almost surely, the co-range of $\mathbf {\mathsf{\Gamma }}$ is a uniformly random d-dimensional subspace of $\mathbb {F}^n$ . Construct an embedding $\mathbf {\mathsf {S}} \in \mathbb {F}^{d \times n}$ with orthonormal rows that span the co-range of $\mathbf {\mathsf{\Gamma }}$ ; for example, by QR factorization.

Similarly, we can consider a Gaussian embedding $\mathbf {\mathsf{\Gamma }} \in \mathbb {F}^{d \times n}$ with $d \geq n$ . In this case the range of $\mathbf {\mathsf{\Gamma }}$ is almost surely a uniformly random n-dimensional subspace of $\mathbb {F}^d$ . Construct an embedding $\mathbf {\mathsf {S}} \in \mathbb {F}^{d \times n}$ with orthonormal columns that span the range of $\mathbf {\mathsf{\Gamma }}$ , for example by QR factorization.

In each case, we call $\mathbf {\mathsf {S}}$ a random partial isometry. The cost of storing a random partial isometry is ${{O}}(dn)$ , and the cost of applying it to a vector is ${{O}}(dn)$ . (We should warn the punctilious reader that QR factorization of $\mathbf {\mathsf{\Gamma }}$ may not produce a matrix $\mathbf {\mathsf {S}}$ that is Haar-distributed on the Stiefel manifold. To achieve this guarantee, use the algorithms from Mezzadri (Reference Mezzadri2007).)

When $d \approx n$ , random partial isometries are better embeddings than Gaussian maps (because the nonzero singular values of a partial isometry are all equal). When d and n are significantly different, the two models are quite similar to each other.

Thrampoulidis and Hassibi (Reference Thrampoulidis and Hassibi2015) have established some theoretical results on the embedding behaviour of real partial isometries $(\mathbb {F} = {\mathbb {R}})$ . Unfortunately, the situation is more complicated than in Theorem 8.4. More relations between Gaussian matrices and partial isometries follow from the Marcus–Pisier comparison theorem (Marcus and Pisier Reference Marcus and Pisier1981); see also Tropp (Reference Tropp2012a ).

9.1. General techniques

Many types of structured random embeddings operate on the same principle, articulated in Ailon and Chazelle (Reference Ailon and Chazelle2009). When we apply the random embedding to a fixed vector, it should homogenize (‘mix’) the coordinates so that each one carries about the same amount of energy. Then the embedding can sample coordinates at random to extract a lower-dimensional vector whose norm is proportional to the norm of the original vector and has low variance. Random embeddings differ in how they perform the initial mixing step. Regardless of how it is done, mixing is very important for obtaining embeddings that work well in practice.

With this intuition at hand, let us introduce a number of pre- and post-processing transforms that help us design effective random embeddings. These approaches are used in many of the constructions below.

We say that a random variable is a random sign if it is $ {\small\rm UNIFORM}\{ z \in \mathbb {F} \colon {\vert {z} \vert } = 1 \}$ . A random matrix $\mathbf {\mathsf {E}} \in \mathbb {F}^{n \times n}$ is called a random sign flip if it is diagonal, and the diagonal entries are i.i.d. random sign variables.

A random permutation $\mathbf {{\Pi }} \in \mathbb {F}^{n \times n}$ is a matrix drawn uniformly at random from the set of permutation matrices. That is, each row and column of $\mathbf {{\Pi }}$ has a single nonzero entry, which equals one, and all such matrices are equally likely.

For $d \leq n$ , a random matrix $\mathbf {\mathsf {R}} \in \mathbb {F}^{d \times n}$ is called a random restriction if it selects d uniformly random entries from its input. With an abuse of terminology, we extend the definition to the case $d \geq n$ by making $\mathbf {\mathsf {R}} \in \mathbb {F}^{d \times n}$ the matrix that embeds its input into the first n coordinates of the output. That is, $(\mathbf {\mathsf {R}})_{ij} = 1$ when $i = j$ and zero otherwise.

Random sign flips and permutations are useful for preconditioning the input to a random embedding. Random restrictions are useful for reducing the dimension of a vector that has already been homogenized.

9.2. Sparse sign matrices

Among the earliest proposals for non-Gaussian embedding is to use a sparse random matrix whose entries are random signs.

Here is an effective construction of a sparse sign matrix $\mathbf {\mathsf {S}} \in \mathbb {F}^{d \times n}$ . Fix a sparsity parameter $\zeta $ in the range $2 \leq \zeta \leq d$ . The random embedding takes the form

$$ \begin{align*} \mathbf{\mathsf{S}} = \sqrt{\displaystyle \frac{n}{\zeta}} \begin{bmatrix} \mathbf{\mathsf{s}}_1 & \cdots & \mathbf{\mathsf{s}}_n \end{bmatrix} \in \mathbb{F}^{d \times n}. \end{align*} $$

The columns $\mathbf {\mathsf {s}}_i \in \mathbb {F}^d$ are i.i.d. random vectors. To construct each column, we draw $\zeta $ i.i.d. random signs, and we situate them in $\zeta $ uniformly random coordinates. Tropp, Yurtsever, Udell and Cevher (Reference Tropp, Yurtsever, Udell and Cevher2019) recommend choosing $\zeta = \min \{d, 8\}$ in practice.

We can store a sparse embedding using about ${{O}}(\zeta n \log d)$ numbers. We can apply it to a vector in $\mathbb {F}^n$ with ${{O}}(\zeta n)$ arithmetic operations. The main disadvantage is that we must use sparse data structures and arithmetic to achieve these benefits. Sparse sign matrices have similar performance to Gaussian embeddings.

Cohen (Reference Cohen2016) has shown that a sparse sign matrix serves as an oblivious subspace embedding with constant distortion for an arbitrary k-dimensional subspace of ${\mathbb {R}}^n$ when the embedding dimension $d = {{O}}(k \log k)$ and the per-column sparsity $\zeta = {{O}}(\log k)$ . It is conjectured that improvements are still possible.

9.3. Subsampled trigonometric transforms

Another type of structured randomized embeddings is designed to mimic the performance of a random partial isometry. One important class of examples consists of the subsampled randomized trigonometric transforms (SRTTs).

To construct a random embedding $\mathbf {\mathsf {S}} \in \mathbb {F}^{d \times n}$ with $d \leq n$ , we select a unitary trigonometric transform $\mathbf {\mathsf {F}} \in \mathbb {F}^{n \times n}$ . Then we form

$$ \begin{align*} \mathbf{\mathsf{S}} = \sqrt{\displaystyle \frac{n}{d}} \mathbf{\mathsf{RFE}}\mathbf{{\Pi}}, \end{align*} $$

where $\mathbf {\mathsf {R}} \in \mathbb {F}^{d\times n}$ is a random restriction, $\mathbf {\mathsf {E}} \in \mathbb {F}^{n \times n}$ is a random sign flip, and $\mathbf {{\Pi }} \in \mathbb {F}^{n \times n}$ is a random permutation. Note that $\mathbf {\mathsf {S}}$ is a partial isometry.

The trigonometric transform $\mathbf {\mathsf {F}}$ can be any one of the usual suspects. In the complex case ( $\mathbb {F} = \mathbb {C}$ ), we often use a discrete Fourier transform. In the real case ( $\mathbb {F} = {\mathbb {R}}$ ), common choices are the discrete cosine transform (DCT2) or the discrete Hartley transform (DHT). When n is a power of two, we can consider a Walsh–Hadamard transform (WHT). The paper (Avron et al. Reference Avron, Maymounkov and Toledo2010) reports that the DHT is the best option in the real case.

The cost of storing a SRTT is ${{O}}(n \log n)$ , and it can be applied to a vector in ${{O}}(n \log d)$ operations using a fast subsampled trigonometric transform algorithm. The main disadvantage is that it requires a good implementation of the fast transform.

Tropp (Reference Tropp2011b ) has shown that an SRTT serves as an oblivious subspace embedding with constant distortion for an arbitrary k-dimensional subspace of $\mathbb {F}^n$ provided that $d = {{O}}(k \log k)$ . This paper focuses on the Walsh–Hadamard transform, but the analysis extends to other SRTTs. In practice, it often suffices to choose $d = {{O}}(k)$ , but no rigorous justification is available.

Remark 9.2 (rerandomization)

It is also common to repeat the randomization and trigonometric transformations:

$$ \begin{align*} \mathbf{\mathsf{S}} = \sqrt{\displaystyle \frac{n}{d}} \mathbf{\mathsf{R}}\mathbf{\mathsf{F}}{\mathbf{\mathsf{E}}^{\prime}}{\mathbf{{\Pi}}^{\prime}}\mathbf{\mathsf{FE}}\mathbf{{\Pi}}, \end{align*} $$

with an independent sign flip $\mathbf {\mathsf {E}}'$ and an independent permutation $\mathbf {{\Pi }}'$ . This enhancement can make the embedding more robust, although it is not always necessary.

9.4. Tensor random projections

Next, we describe a class of random embeddings that are useful for very large linear algebra and multilinear algebra problems. This approach invokes tensor products to form a random embedding for a high-dimensional space from a family of random embeddings for lower-dimensional spaces.

Let $\mathbf {\mathsf {S}}_1 \in \mathbb {F}^{d \times m_1}$ and $\mathbf {\mathsf {S}}_2 \in \mathbb {F}^{d \times m_2}$ be statistically independent random embeddings. We define the tensor random embedding

$$ \begin{align*} \mathbf{\mathsf{S}} := \mathbf{\mathsf{S}}_1 \odot \mathbf{\mathsf{S}}_2 \in \mathbb{F}^{d \times n} \quad\text{where } n = m_1 m_2 \end{align*} $$

to be the Khatri–Rao product of $\mathbf {\mathsf {S}}_1$ and $\mathbf {\mathsf {S}}_2$ . That is, the ith row of $\mathbf {\mathsf {S}}$ is

$$ \begin{align*} (\mathbf{\mathsf{S}})_{i:} = \begin{bmatrix} (\mathbf{\mathsf{S}}_1)_{i1} (\mathbf{\mathsf{S}}_2)_{i:} & \cdots & (\mathbf{\mathsf{S}}_1)_{im} (\mathbf{\mathsf{S}}_2)_{i:} \end{bmatrix} \quad\text{for }i=1,\ldots,d. \end{align*} $$

Under moderate assumptions on the component embeddings $\mathbf {\mathsf {S}}_1$ and $\mathbf {\mathsf {S}}_2$ , the tensor random embedding $\mathbf {\mathsf {S}}$ preserves the squared Euclidean norm of an arbitrary vector in $\mathbb {F}^{n}$ .

A natural extension of this idea is to draw many component embeddings $\mathbf {\mathsf {S}}_i \in \mathbb {F}^{d \times m_i}$ for $i = 1, \ldots , k$ and to form the tensor random embedding

$$ \begin{align*} \mathbf{\mathsf{S}} := \mathbf{\mathsf{S}}_1 \odot \mathbf{\mathsf{S}}_2 \odot \cdots \odot \mathbf{\mathsf{S}}_k \in \mathbb{F}^{d \times n} \quad\text{where } n = \prod\nolimits_{i=1}^k m_i. \end{align*} $$

This embedding also inherits nice properties from its components.

The striking thing about this construction is that the tensor product embedding operates on a much larger space than the component embeddings. The storage cost for the component embeddings is ${{O}}(d (\sum \nolimits _{i=1}^k m_i))$ , or less. We can apply the tensor random embedding to a vector directly with ${{O}}(dn)$ arithmetic. We can accelerate the process by using component embeddings that have fast transforms, and we can obtain improvements for vectors that have a compatible tensor product structure. Some theoretical analysis is available, but results are not yet complete.

9.5. Other types of structured random embeddings

We have described the random embeddings that have received the most attention in the NLA literature. Yet there are other types of random embeddings that may be useful in special circumstances. Some examples include random filters (Tropp et al. Reference Tropp, Wakin, Duarte, Baron and Baraniuk2006, Krahmer and Ward Reference Krahmer and Ward2011, Rauhut, Romberg and Tropp Reference Tropp2012, Mendelson, Rauhut and Ward Reference Mendelson, Rauhut and Ward2018), the Kac random walk (Kac Reference Kac1956, Rosenthal Reference Rosenthal1994, Oliveira Reference Oliveira2009 b, Pillai and Smith Reference Pillai and Smith2017) and sequences of random reflections (Sloane Reference Sloane1983, Porod Reference Porod1996). Liberty (Reference Liberty2009) discusses a number of other instances.

9.6. Coordinate sampling

So far we have discussed random embeddings that mix up the coordinates of a vector. It is sometimes possible to construct embeddings just by sampling coordinates at random. Coordinate sampling can be appealing in specialized situations (e.g. kernel computations), where we only have access to individual entries of the data. On the other hand, this approach requires strong assumptions, and it is far less reliable than random embeddings that mix coordinates. In this section we summarize the basic facts about subspace embedding via random coordinate sampling.

A note on terminology: we will use the term coordinate sampling to distinguish these maps from random embeddings that mix coordinates.

9.6.1. Coherence and leverage

Let $L \subset \mathbb {F}^{n}$ be a k-dimensional subspace. The coherence $\mu (L)$ of the subspace with respect to the standard basis is

$$ \begin{align*} \mu(L) := n \cdot \max_{i = 1, \ldots, n} {\Vert { {\mathbf{\mathsf{P}}}_L {\mathbf{{\delta}}}_i } \Vert}^2, \end{align*} $$

where $\mathbf {\mathsf {P}}_L \in \mathbb {H}_n$ is the orthogonal projector onto L and $\mathbf {{\delta }}_i$ is the ith standard basis vector. The coherence $\mu (L)$ lies in the range $[k, n]$ . The behaviour of coordinate sampling methods degrades as the coherence increases.

Next, define the (subspace) leverage score distribution with respect to the standard coordinate basis:

$$ \begin{align*} p_i = \displaystyle \frac{1}{k} {\Vert { {\mathbf{\mathsf{P}}}_L {\mathbf{{\delta}}}_i } \Vert}^2 \quad\text{for }i = 1, \ldots, n. \end{align*} $$

It is straightforward to verify that $(p_1, \ldots , p_n)$ is a probability distribution. In most applications, it is expensive to compute or estimate subspace leverage scores because we typically do not have a basis for the subspace L at hand.

9.6.2. Uniform sampling

We can construct an embedding $\mathbf {\mathsf {S}} \in \mathbb {F}^{d \times n}$ by sampling each output coordinate uniformly at random. That is, the rows of $\mathbf {\mathsf {S}}$ are i.i.d., and each row takes values $\mathbf {{\delta }}_i / \sqrt {d}$ , each with probability $1/n$ . (We can also sample coordinates uniformly without replacement; this approach performs slightly better but requires more work to analyse.)

The embedding dimension d must be chosen to ensure that

(9.1) $$ \begin{align} {\Vert {\mathbf{\mathsf{S}} \mathbf{\mathsf{x}} } \Vert}^2 = (1 \pm \varepsilon) {\Vert \mathbf{\mathsf{x}} \Vert}^2 \quad\text{for all }\mathbf{\mathsf{x}} \in L. \end{align} $$

To achieve this goal, it suffices that

$$ \begin{align*} d \geq 2 \varepsilon^{-2} \mu(L) \log(2k). \end{align*} $$

In other words, the embedding dimension is proportional to the coherence of the subspace, up to a logarithmic factor. We expect uniform sampling to work well precisely when the coherence is small ( $\mu (L) \approx k$ ).

To prove this result, let $\mathbf {\mathsf {U}} \in \mathbb {F}^{n \times k}$ be an orthonormal basis for the subspace L. We can approximate the product ${\mathbf {\mathsf {I}}}_k = \mathbf {\mathsf {U}}^{*}\mathbf {\mathsf {U}}$ by sampling columns of $\mathbf {\mathsf {U}}$ uniformly at random. The analysis in Section 7.3.2 furnishes the conclusion.

9.6.3. Leverage score sampling

When the coherence is large, it seems more natural to sample with respect to the leverage score distribution $(p_1, \ldots , p_n)$ described above. That is, the embedding $\mathbf {\mathsf {S}} \in \mathbb {F}^{d \times n}$ has i.i.d. rows, and each row takes value $\mathbf {{\delta }}_i / \sqrt {d}$ with probability $p_i$ .

To achieve the embedding guarantee (9.1) with this sampling distribution, we should choose the embedding dimension

$$ \begin{align*} d \geq 2 \varepsilon^{-2} k \log(2k). \end{align*} $$

In other words, it suffices that the embedding dimension is proportional to the dimension k of the subspace, up to the logarithmic factor. This result follows from the analysis of matrix multiplication by importance sampling (Section 7.3.3).

9.7. But how does it work in theory?

Structured random embeddings and random coordinate sampling lack the precise guarantees that we can attribute to Gaussian embedding matrices. So how can we apply them with confidence?

First, we advocate using a posteriori error estimators to assess the quality of the output of a randomized linear algebra computation. These error estimators are often quite cheap, yet they can give (statistical) evidence that the computation was performed correctly. We also recommend adaptive algorithms that can detect when the accuracy is insufficient and make refinements. With this approach, it is not pressing to produce theory that justifies all of the internal choices (e.g. the specific type of random embedding) in the NLA algorithm. See Section 12 for further discussion.

Even so, we would like to have a priori predictions about how our algorithms will behave. Beyond that, we need reliable methods for selecting algorithm parameters, especially in the streaming setting where we cannot review the data and repeat the computation.

Here is one answer to these concerns. As a practical matter, we can simply invoke the lessons from the Gaussian theory, even when we are using a different type of random embedding. The universality result, Theorem 8.7, gives a rationale for this approach in one special case. We also recommend undertaking computational experiments to verify that the Gaussian theory gives an adequate description of the observed behaviour of an algorithm.

10.1. Overdetermined least-squares

Overdetermined least-squares problems sometimes arise in statistics and data-analysis applications. We may imagine that some of the data in these problems is redundant. As such, it seems plausible that we could reduce the size of the problem to accelerate computation without too much loss in accuracy.

Consider a matrix $\mathbf {\mathsf {A}} \in \mathbb {F}^{m \times n}$ with $m \gg n$ and a vector $\mathbf {\mathsf {b}} \in \mathbb {F}^m$ . An overdetermined least-squares problem has the form

(10.1) $$ \begin{align} \underset{\mathbf{\mathsf{x}} \in \mathbb{F}^n}{\text{minimize}} \ \displaystyle \frac{1}{2} {\Vert { \mathbf{\mathsf{A}} \mathbf{\mathsf{x}} - \mathbf{\mathsf{b}} } \Vert}^2. \end{align} $$

Following Pilanci and Wainwright (Reference Pilanci and Wainwright2015), we can also rewrite the least-squares problem to emphasize the role of the matrix:

(10.2) $$ \begin{align} \underset{\mathbf{\mathsf{x}} \in \mathbb{F}^n}{\text{minimize}} \ \displaystyle \frac{1}{2} {\Vert { \mathbf{\mathsf{A}} \mathbf{\mathsf{x}} } \Vert}^2 - {\langle { \mathbf{\mathsf{x}} }, \, { \mathbf{\mathsf{A}}^{*} \mathbf{\mathsf{b}} } \rangle}. \end{align} $$

Write $\mathbf {\mathsf {x}}_{\star }$ for an arbitrary solution to the problem (10.1).

To make a clear comparison among algorithm design templates, we will assume that $\mathbf {\mathsf {A}}$ is dense and unstructured. In this case the classical approach to solving (10.1) is based on factorization of the coefficient matrix (such as QR or SVD) at a cost of ${{O}}(mn^2)$ arithmetic operations.

When $\mathbf {\mathsf {A}}$ is sparse, we would typically use iterative methods (such as CG), which have a different computational profile. For sparse matrices, we would also make different design choices in a sketch-based algorithm. Nevertheless, for simplicity, we will not discuss the sparse case.

10.2. Subspace embeddings for least-squares

To design a sketching algorithm for the overdetermined least-squares problem (10.1), we need to construct a subspace embedding $\mathbf {\mathsf {S}} \in \mathbb {F}^{d \times m}$ that preserves the geometry of the range of the matrix $\mathbf {\mathsf {A}} \in \mathbb {F}^{m \times n}$ . In some cases we may also need the embedding to preserve the range of the bordered matrix $\begin {bmatrix} \mathbf {\mathsf {A}} & \mathbf {\mathsf {b}} \end {bmatrix} \in \mathbb {F}^{m \times (n+1)}$ .

Since the matrix $\mathbf {\mathsf {A}}$ is dense and unstructured, we will work with a structured subspace embedding, such as an SRTT (Section 9.3). More precisely, we will assume that evaluating the product $\mathbf {\mathsf {SA}}$ costs only ${{O}}(mn \log d)$ arithmetic operations. The best theoretical results for these structured sketches require that the embedding dimension $d \sim n \log (n)/\varepsilon ^2$ to achieve distortion $\varepsilon $ , although the logarithmic factor seems to be unnecessary in practice.

Throughout this section we use the heuristic notation $\sim $ to indicate quantities that are proportional. We also write $\ll $ to mean ‘much smaller than’.

10.3. Sketch and solve

The sketch-and-solve paradigm maps the overdetermined least-squares problem (10.1) into a smaller space. Then it uses the solution to the reduced problem as a proxy for the solution to the original problem. This approach can be very fast, and we only need one view of the matrix $\mathbf {\mathsf {A}}$ . On the other hand, the results tend to be very inaccurate.

Let $\mathbf {\mathsf {S}} \in \mathbb {F}^{d \times m}$ be a subspace embedding for the range of $\begin {bmatrix} \mathbf {\mathsf {A}} & \mathbf {\mathsf {b}} \end {bmatrix}$ with distortion $\varepsilon $ . Consider the compressed least-squares problem

(10.3) $$ \begin{align} \underset{\mathbf{\mathsf{x}} \in {\mathbb{R}}^n}{\text{minimize}} \ \displaystyle \frac{1}{2} {\Vert { \mathbf{\mathsf{S}} (\mathbf{\mathsf{A}} \mathbf{\mathsf{x}} - \mathbf{\mathsf{b}}) } \Vert}^2. \end{align} $$

Since $\mathbf {\mathsf {S}}$ preserves geometry, we may hope that the solution $\widehat {\mathbf {\mathsf {x}}}$ to the sketched problem (10.3) can replace the solution $\mathbf {\mathsf {x}}_{\star }$ to the original problem (10.2). A typical theoretical bound is

$$ \begin{align*} {\Vert { \mathbf{\mathsf{A}} \widehat{\mathbf{\mathsf{x}}} - \mathbf{\mathsf{b}} } \Vert} \leq (1 + \varepsilon) {\Vert { \mathbf{\mathsf{A}} \mathbf{\mathsf{x}}_{\star} - \mathbf{\mathsf{b}} } \Vert} \quad\text{when }d \sim n \log(n) / \varepsilon^2. \end{align*} $$

See Sarlós (Reference Sarlós2006). Although the residuals are comparable, it need not be the case that $\widehat {\mathbf {\mathsf {x}}} \approx \mathbf {\mathsf {x}}_{\star }$ , even when the solution to (10.1) is unique.

The sketch-and-solve paradigm requires us to form the matrix $\mathbf {\mathsf {SA}}$ at a cost of ${{O}}(mn \log d)$ operations. We would typically solve the (dense) reduced problem with a direct method, using ${{O}}(dn^2)$ operations. Assuming $d \sim n \log (n)/\varepsilon ^2$ , the total arithmetic cost is ${{O}}(mn \log (n/\varepsilon ^2) + n^3 \log (n)/\varepsilon ^2)$ .

In summary, we witness an improvement in computational cost over classical methods if $\log n \ll n \ll m/\log n$ and $\varepsilon $ is constant. But we must also be willing to accept large errors, because we cannot make $\varepsilon $ small.

10.4. Iterative sketching

Iterative sketching attempts to remediate the poor accuracy of the sketch-and-solve paradigm by using it repeatedly to reduce the residual error.

First we construct an initial solution $\mathbf {\mathsf {x}}_0 \in \mathbb {F}^n$ using the sketch-and-solve paradigm with a constant distortion embedding. For each iteration i, draw a fresh random subspace embedding $\mathbf {\mathsf {S}}_i \in \mathbb {F}^{d \times m}$ for $\operatorname {range}(\mathbf {\mathsf {A}})$ , with constant distortion. We can solve a sequence of least-squares problems

(10.4) $$ \begin{align} \underset{\mathbf{\mathsf{x}} \in {\mathbb{R}}^n}{\text{minimize}} \ \displaystyle \frac{1}{2} {\Vert { \mathbf{\mathsf{S}}_i \mathbf{\mathsf{A}}(\mathbf{\mathsf{x}} - \mathbf{\mathsf{x}}_{i-1}) } \Vert}^2 + {\langle { \mathbf{\mathsf{x}} - \mathbf{\mathsf{x}}_{i-1} }, \, { \mathbf{\mathsf{A}}^{*} (\mathbf{\mathsf{b}} - \mathbf{\mathsf{A}}\mathbf{\mathsf{x}}_{i-1}) } \rangle}. \end{align} $$

The solution $\mathbf {\mathsf {x}}_{i}$ to this subproblem is fed into the next subproblem. Without the sketch, each subproblem is equivalent to solving (10.2) with $\mathbf {\mathsf {b}}$ replaced by the residual $\mathbf {\mathsf {r}}_{i-1} = \mathbf {\mathsf {b}} - \mathbf {\mathsf {A}} \mathbf {\mathsf {x}}_{i-1}$ . The sketch $\mathbf {\mathsf {S}}_i$ preserves the geometry while reducing the problem size. A typical theoretical error bound would be

$$ \begin{align*} {\Vert { \mathbf{\mathsf{A}} \mathbf{\mathsf{x}}_j - \mathbf{\mathsf{b}} } \Vert} \leq (1 + \varepsilon) {\Vert { \mathbf{\mathsf{A}} \mathbf{\mathsf{x}}_{\star} - \mathbf{\mathsf{b}} } \Vert} \quad\text{when }j \sim \log(1/\varepsilon)\text{ and }d \sim n \log n. \end{align*} $$

See Pilanci and Wainwright (Reference Pilanci and Wainwright2016) for related results.

In each iteration, the iterative sketching approach requires us to form $\mathbf {\mathsf {S}}_i \mathbf {\mathsf {A}}$ at a cost of ${{O}}(mn \log d)$ . We compute $\mathbf {\mathsf {A}}^{*} \mathbf {\mathsf {r}}_{i-1}$ at a cost of ${{O}}(mn)$ . Although it is unnecessary to solve each subproblem accurately, we cannot obtain reliable behaviour without using a dense method at a cost of ${{O}}(dn^2)$ per iteration. With the theoretical parameter choices, the total arithmetic is ${{O}}((mn + n^3) \log (n) \log (1/\varepsilon ))$ to achieve relative error $\varepsilon $ .

The interesting parameter regime is $\log n \ll n \ll m / \log n$ , but we can now allow $\varepsilon $ to be tiny. In this setting, iterative sketching costs slightly more than the sketch-and-solve paradigm to achieve constant relative error, while it is faster than the classical approach. At the same time, it can produce errors as small as traditional least-squares algorithms. A shortcoming is that this method requires repeated sketches of the matrix $\mathbf {\mathsf {A}}$ .

For overdetermined least-squares, we can short-circuit the iterative sketching approach. In this setting, we can sketch the input matrix just once and factorize it. We can use the same factorized sketch in each iteration to solve the subproblems faster. For problems more general than least-squares, it may be necessary to extract a fresh sketch at each iteration, as we have done here.

10.5. Sketch and precondition

The sketch-and-precondition paradigm uses random embedding to find a proxy for the input matrix. We can use this proxy to precondition a classical iterative algorithm so it converges in a minimal number of iterations.

Let $\mathbf {\mathsf {S}} \in \mathbb {F}^{d \times m}$ be a subspace embedding for $\operatorname {range}(\mathbf {\mathsf {A}})$ with constant distortion. Compress the input matrix $\mathbf {\mathsf {A}}$ , and then compute a (pivoted) QR factorization:

$$ \begin{align*} \mathbf{\mathsf{Y}} = \mathbf{\mathsf{SA}} \quad\text{and}\quad \mathbf{\mathsf{Y}} = \mathbf{\mathsf{QR}}. \end{align*} $$

Since $\mathbf {\mathsf {S}}$ preserves the range of $\mathbf {\mathsf {A}}$ when $d \sim n \log n$ , we anticipate that ${\mathbf {\mathsf {Y}}^{*}} \mathbf {\mathsf {Y}} \approx \mathbf {\mathsf {A}}^{*}\mathbf {\mathsf {A}}$ . As a consequence, $\mathbf {\mathsf {A}}\mathbf {\mathsf {R}}^{\dagger }$ should be close to an isometry. Thus we can pass to the preconditioned problem

(10.5) $$ \begin{align} \underset{\mathbf{\mathsf{x}} \in {\mathbb{R}}^n}{\text{minimize}} \ \displaystyle \frac{1}{2} {\Vert { (\mathbf{\mathsf{A}}\mathbf{\mathsf{R}}^{\dagger})(\mathbf{\mathsf{R}}\mathbf{\mathsf{x}}) - \mathbf{\mathsf{b}} } \Vert}^2. \end{align} $$

Construct an initial solution $\mathbf {\mathsf {x}}_0 \in \mathbb {F}^n$ using the sketch-and-solve paradigm with the embedding $\mathbf {\mathsf {S}}$ . From this starting point, we solve (10.5) using preconditioned LQSR. The jth iterate satisfies

$$ \begin{align*} {\Vert { \mathbf{\mathsf{A}} \mathbf{\mathsf{x}}_j - \mathbf{\mathsf{b}} } \Vert} \leq (1 + \varepsilon) {\Vert { \mathbf{\mathsf{A}} \mathbf{\mathsf{x}}_{\star} - \mathbf{\mathsf{b}} } \Vert} \quad\text{when }j \sim \log(1/\varepsilon). \end{align*} $$

This statement is a reinterpretation of the theory in Rokhlin and Tygert (Reference Rokhlin and Tygert2008).

The cost of sketching the input matrix and performing the QR decomposition is ${{O}}(mn \log d + dn^2)$ . Afterwards, we pay ${{O}}(mn)$ for each iteration of PCG. With the theoretical parameter settings, the total cost is ${{O}}(mn \log (n/\varepsilon ) + n^3 \log n)$ operations.

Once again, the interesting regime is $\log n \ll n \ll m / \log n$ , and the value of $\varepsilon $ can be very small. For overdetermined least-squares, this approach is faster than both the sketch-and-solve paradigm and the iterative sketching paradigm. The sketch-and-precondition approach leads to errors that are comparable with classical linear algebra algorithms, but it may be a factor of $n / \log (n)$ faster. On the other hand, it requires repeated applications of the matrix $\mathbf {\mathsf {A}}$ .

10.6. Comparisons

If we seek a high-precision solution to a dense, unstructured, overdetermined least-squares problem, randomized preconditioning leads to the most efficient existing algorithm. For the same problem, if we can only view the input matrix once, then the sketch-and-solve paradigm still allows us to obtain a low-accuracy solution. Although iterative sketching is less efficient than its competitors in this setting, it remains useful for solving constrained least-squares problems, and it has further connections with optimization.

10.7. Summary

From the perspective of a numerical analyst, randomized preconditioning and iterative sketching should be the preferred methods for designing sketching algorithms because they allow for high precision. The sketch-and-solve approach is appropriate only when data access is severely constrained.

In spite of this fact, a majority of the literature on randomized NLA develops algorithms based on the sketch-and-solve paradigm. There are far fewer works on randomized preconditioning or iterative sketching. This discrepancy points to an opportunity for further research.

11.1. The rangefinder: Problem statement

Let $\mathbf {\mathsf {B}} \in \mathbb {F}^{m \times n}$ be an input matrix, and let $\ell \leq \min \{m, n\}$ be the subspace dimension. The goal of the rangefinder problem is to produce an orthonormal matrix $\mathbf {\mathsf {Q}} \in \mathbb {F}^{m \times \ell }$ whose range aligns with the dominant left singular vectors of $\mathbf {\mathsf {B}}$ .

To measure the quality of $\mathbf {\mathsf {Q}}$ , we use the spectral norm error

(11.1) $$ \begin{align} {\Vert { \mathbf{\mathsf{B}} - \mathbf{\mathsf{QQ}}^{*} \mathbf{\mathsf{B}} } \Vert} = {\Vert { ({\mathbf{\mathsf{I}}} - \mathbf{\mathsf{QQ}}^{*}) \mathbf{\mathsf{B}} } \Vert}. \end{align} $$

If the error measure (11.1) is small, then the rank- $\ell $ matrix $\widehat {\mathbf {\mathsf {B}}} = \mathbf {\mathsf {QQ}}^{*} \mathbf {\mathsf {B}}$ can serve as a proxy for $\mathbf {\mathsf {B}}$ . See Section 11.2 for an important application.

11.1.1. The randomized rangefinder: A pseudoalgorithm

Using randomized methods, it is remarkably easy to find an initial solution to the rangefinder problem. We simply multiply the target matrix by a random embedding and then orthogonalize the resulting matrix.

More rigorously: consider a target matrix $\mathbf {\mathsf {B}} \in \mathbb {F}^{m \times n}$ and a subspace dimension $\ell $ . We draw a random test matrix $\mathbf {\mathsf{\Omega }} \in \mathbb {F}^{n \times \ell }$ , where $\mathbf {\mathsf{\Omega }}^{*}$ is a mixing random embedding. We form the product $\mathbf {\mathsf {Y}} = \mathbf {\mathsf {B}}\mathbf {\mathsf{\Omega }} \in \mathbb {F}^{m \times \ell }$ . Then we compute an orthonormal basis $\mathbf {\mathsf {Q}} \in \mathbb {F}^{m \times \ell }$ for the range of $\mathbf {\mathsf {Y}}$ using a QR factorization method. See Algorithm 7 for pseudocode.

In a general setting, the arithmetic cost of this procedure is dominated by $O (mn\ell )$ operations for the matrix–matrix multiplication. The QR factorization of $\mathbf {\mathsf {Y}}$ requires $O (m\ell ^2)$ arithmetic, and we also need to simulate the $n \times \ell $ random matrix $\mathbf {\mathsf{\Omega }}$ . Economies are possible when either $\mathbf {\mathsf {B}}$ or $\mathbf {\mathsf{\Omega }}$ admits fast multiplication.

11.2. The randomized singular value decomposition (RSVD)

Before we continue with our discussion of the rangefinder, let us summarize one of the key applications: the randomized SVD algorithm.

Low-rank approximation problems often arise when a user seeks an incomplete matrix factorization that exposes structure, such as a truncated eigenvalue decomposition or a partial QR factorization. The randomized rangefinder, described in Section 11.1.1, can be used to perform the heavy lifting in these computations. Afterwards, we perform some light post-processing to reach the desired factorization.

To illustrate how this works, suppose that we want to compute an approximate rank- $\ell $ truncated singular value decomposition of the input matrix $\mathbf {\mathsf {B}} \in \mathbb {F}^{m \times n}$ . That is,

$$ \begin{align*} \mathbf{\mathsf{B}} \approx \mathbf{\mathsf{U}}\mathbf{\mathsf{\Sigma}}\mathbf{\mathsf{V}}^{*}, \end{align*} $$

where $\mathbf {\mathsf {U}} \in \mathbb {F}^{m \times \ell }$ and $\mathbf {\mathsf {V}} \in \mathbb {F}^{n \times \ell }$ are orthonormal matrices and $\mathbf {\mathsf{\Sigma }} \in \mathbb {F}^{\ell \times \ell }$ is a diagonal matrix whose diagonal entries approximate the largest singular values of $\mathbf {\mathsf {A}}$ .

Choose a target rank $\ell $ , and suppose that $\mathbf {\mathsf {Q}} \in \mathbb {F}^{m \times \ell }$ is a computed solution to the rangefinder problem. The rangefinder furnishes an approximate rank- $\ell $ factorization of the input matrix: $\mathbf {\mathsf {B}} \approx \mathbf {\mathsf {Q}}(\mathbf {\mathsf {Q}}^{*}\mathbf {\mathsf {A}})$ . To convert this representation into a truncated SVD, we just compute an economy-size SVD of the matrix $\mathbf {\mathsf {C}} := \mathbf {\mathsf {Q}}^{*}\mathbf {\mathsf {A}} \in \mathbb {F}^{\ell \times n}$ and consolidate the factors.

In symbols, once $\mathbf {\mathsf {Q}}$ is available, the computation proceeds as follows:

$$\begin{equation*}\begin{array}{l@{}ll@{}l} \mathbf{\mathsf{B}} &\approx \mathbf{\mathsf{Q}}\mathbf{\mathsf{Q}}^{*}\mathbf{\mathsf{B}} &\quad &\{\mbox{matrix--matrix multiplication: }\mathbf{\mathsf{C}} = \mathbf{\mathsf{Q}}^{*}\mathbf{\mathsf{B}}\} \\[1pt] &= \mathbf{\mathsf{Q}}\mathbf{\mathsf{C}} &&\{\mbox{{economy-size} SVD: }\mathbf{\mathsf{C}} = \widehat{\mathbf{\mathsf{U}}}\mathbf{\mathsf{\Sigma}}\mathbf{\mathsf{V}}^{*}\} \\[1pt] &= \mathbf{\mathsf{Q}}\widehat{\mathbf{\mathsf{U}}}\mathbf{\mathsf{\Sigma}}\mathbf{\mathsf{V}}^{*} &&\{\mbox{matrix--matrix multiplication: }\mathbf{\mathsf{U}} = \mathbf{\mathsf{Q}}\widehat{\mathbf{\mathsf{U}}}\} \\[1pt] &= \mathbf{\mathsf{U}}\mathbf{\mathsf{\Sigma}}\mathbf{\mathsf{V}}^{*}. \end{array}\end{equation*}$$

After the rangefinder step, the remaining computations are all exact (modulo floating-point arithmetic errors). Therefore

$$ \begin{align*} \|\mathbf{\mathsf{B}} - \mathbf{\mathsf{U}}\mathbf{\mathsf{\Sigma}}\mathbf{\mathsf{V}}^{*}\| = \|\mathbf{\mathsf{B}} - \mathbf{\mathsf{Q}}\mathbf{\mathsf{Q}}^{*}\mathbf{\mathsf{B}}\|. \end{align*} $$

In words, the accuracy of the approximate SVD is determined entirely by the error in the rangefinder computation!

Empirically, the smallest singular values and singular vectors of the approximate SVD contribute to the accuracy of the approximation, but they are not good estimates for the true singular values and vectors of the matrix. Therefore it can be valuable to truncate the rank by zeroing out the smallest computed singular values. We omit the details. See Halko et al. (Reference Halko, Martinsson and Tropp2011a ), Gu (Reference Gu2015) and Tropp et al. (Reference Tropp, Yurtsever, Udell and Cevher2019) for further discussion.

Algorithm 8 contains pseudocode for the randomized SVD. In a general setting, the dominant cost after the rangefinder step is the matrix–matrix multiply, which requires $O (mn\ell )$ operations. The storage requirements are $O ((m+n)\ell )$ numbers.

The rangefinder primitive allows us to perform other matrix computations as well. For example, in Section 13 we explain how to use the rangefinder to construct matrix factorizations where a subset of the rows/columns are picked to form a basis for the row/column spaces. This approach gives far better results than the more obvious randomized algorithms based on coordinate sampling.

11.3. The rangefinder and Schur complements

Why does the randomized rangefinder work? We will demonstrate that the procedure has its most natural expression in the language of Schur complements. This point is implicit in the analysis in Halko et al. (Reference Halko, Martinsson and Tropp2011a ), and it occasionally appears more overtly in the literature (e.g. in Gittens Reference Gittens2013 and Tropp, Yurtsever, Udell and Cevher Reference Tropp, Yurtsever, Udell, Cevher and Guyon2017a ). Nevertheless, this connection has not been explored in a systematic way.

Proposition 11.1 (rangefinder: Schur complements)

Let $\mathbf {\mathsf {Y}} = \mathbf {\mathsf {BX}}$ for an arbitrary test matrix $\mathbf {\mathsf {X}} \in \mathbb {F}^{n \times \ell }$ , and let $\mathbf {\mathsf {P}}_{\mathbf {\mathsf {Y}}}$ be the orthogonal projector onto the range of $\mathbf {\mathsf {Y}}$ . Define the approximation error as

$$ \begin{align*} \mathbf{\mathsf{E}} := \mathbf{\mathsf{E}}(\mathbf{\mathsf{B}}, \mathbf{\mathsf{X}}) := ({\mathbf{\mathsf{I}}} - \mathbf{\mathsf{P}}_{\mathbf{\mathsf{Y}}}) \mathbf{\mathsf{B}}. \end{align*} $$

Then the squared error can be written as a Schur complement (2.4):

$$ \begin{align*} {\vert {\mathbf{\mathsf{E}}} \vert}^2 := \mathbf{\mathsf{E}}^{*} \mathbf{\mathsf{E}} = (\mathbf{\mathsf{B}}^{*} \mathbf{\mathsf{B}}) / \mathbf{\mathsf{X}}. \end{align*} $$

We emphasize that ${\vert {\mathbf {\mathsf {E}}} \vert }^2$ is a PSD matrix, not a scalar.

Proof. This result follows from a short calculation. We can write the orthogonal projector $\mathbf {\mathsf {P}}_{\mathbf {\mathsf {Y}}}$ in the form

$$ \begin{align*} \mathbf{\mathsf{P}}_{\mathbf{\mathsf{Y}}} = (\mathbf{\mathsf{BX}})((\mathbf{\mathsf{B}}\mathbf{\mathsf{X}})^{*} (\mathbf{\mathsf{BX}}))^{\dagger} (\mathbf{\mathsf{BX}})^{*}. \end{align*} $$

Abbreviating $\mathbf {\mathsf {A}} = \mathbf {\mathsf {B}}^{*} \mathbf {\mathsf {B}}$ , we have

$$ \begin{align*} \mathbf{\mathsf{E}}^{*} \mathbf{\mathsf{E}} = \mathbf{\mathsf{B}}^{*} ({\mathbf{\mathsf{I}}} - \mathbf{\mathsf{P}}_{\mathbf{\mathsf{Y}}}) \mathbf{\mathsf{B}} = \mathbf{\mathsf{A}} - (\mathbf{\mathsf{A}} \mathbf{\mathsf{X}})(\mathbf{\mathsf{X}}^{*} \mathbf{\mathsf{A}} \mathbf{\mathsf{X}})^{\dagger} (\mathbf{\mathsf{AX}})^{*}. \end{align*} $$

This is precisely the definition (2.4) of the Schur complement $\mathbf {\mathsf {A}} / \mathbf {\mathsf {X}}$ .□

Proposition 11.1 gives us access to the deep theory of Schur complements (Zhang Reference Zhang2005). In particular, we have a beautiful monotonicity property that follows instantly from Ando (Reference Ando and Zhang2005, Theorem 5.3).

Corollary 11.2 (monotonicity)

Suppose that

with respect to the semidefinite order

. For each fixed test matrix $\mathbf {\mathsf {X}}$ ,

In particular, the error increases if we increase any singular value of $\mathbf {\mathsf {B}}$ while retaining the same right singular vectors; the error decreases if we decrease any singular value of $\mathbf {\mathsf {B}}$ . The left singular vectors do not play a role here. This observation allows us to identify which target matrices are hardest to approximate.

Example 11.3 (extremals)

Consider the parametrized matrix $\mathbf {\mathsf {B}}(\mathbf {{\sigma }}) = \mathbf {\mathsf {U}} \operatorname {diag}(\mathbf {{\sigma }}) \mathbf {\mathsf {V}}^{*} \in \mathbb {F}^{n \times n}$ , where $\mathbf {\mathsf {U}}, \mathbf {\mathsf {V}}$ are unitary. Suppose that we fix $\sigma _1$ and $\sigma _{k+1}$ . For each test matrix $\mathbf {\mathsf {X}}$ , the error ${\vert {\mathbf {\mathsf {E}}(\mathbf {\mathsf {B}}(\mathbf {{\sigma }}), \mathbf {\mathsf {X}})} \vert }^2$ is maximal in the semidefinite order when

$$ \begin{align*} \mathbf{{\sigma}} = (\underbrace{\sigma_1, \ldots, \sigma_1}_k, \underbrace{\sigma_{k+1}, \ldots, \sigma_{k+1}}_{n-k}). \end{align*} $$

It has long been appreciated that Example 11.3 is the hardest matrix to approximate; see Martinsson et al. (Reference Martinsson, Rokhlin and Tygert2006a , Section 5, Examples 4, 5). The justification of this insight is new.

11.4. A priori error bounds

Proposition 11.1 shows that the error in the rangefinder procedure can be written as a Schur complement. Incredibly, the Schur complement of a PSD matrix with respect to a random subspace tends to be quite small. In this section we summarize a theoretical analysis, due to Halko et al. (Reference Halko, Martinsson and Tropp2011a ), that explains why this claim is true.

11.4.1. Master error bound

First, we present a deterministic upper bound on the error incurred by the rangefinder procedure. This requires some notation.

Without loss of generality, we may assume that $m = n$ by extending $\mathbf {\mathsf {B}}$ with zeros. For any $k \leq \ell $ , construct a partitioned SVD of the target matrix:

$$ \begin{align*} \mathbf{\mathsf{B}} = \mathbf{\mathsf{U}} \begin{bmatrix} \mathbf{\mathsf{\Sigma}}_1 & \\ & \mathbf{\mathsf{\Sigma}}_{2} \end{bmatrix} \begin{bmatrix} \mathbf{\mathsf{V}}_1 & \mathbf{\mathsf{V}}_{2} \end{bmatrix}^{*} \quad\text{with }\mathbf{\mathsf{\Sigma}}_1 \in {\mathbb{R}}^{k \times k}\text{ and }\mathbf{\mathsf{V}}_1 \in \mathbb{F}^{n \times k}. \end{align*} $$

The factors $\mathbf {\mathsf {U}}, \mathbf {\mathsf{\Sigma }}, \mathbf {\mathsf {V}}$ are all square matrices. As usual, the entries of $\mathbf {\mathsf{\Sigma }} = \operatorname {diag}(\sigma _1, \sigma _2, \ldots )$ are arranged in weakly decreasing order. So $\mathbf {\mathsf{\Sigma }}_1$ lists the first k singular values, and $\mathbf {\mathsf{\Sigma }}_2$ lists the remaining $n - k$ singular values. The matrix $\mathbf {\mathsf {V}}_1$ contains the first k right singular vectors; the matrix $\mathbf {\mathsf {V}}_2$ contains the remaining $n - k$ right singular vectors.

For any test matrix $\mathbf {\mathsf {X}} \in \mathbb {F}^{n \times \ell }$ , define

$$ \begin{align*} {\mathbf{\mathsf{X}}}_1 = \mathbf{\mathsf{V}}_1^{*} \mathbf{\mathsf{X}} \quad\text{and}\quad {\mathbf{\mathsf{X}}}_{2} = \mathbf{\mathsf{V}}_{2}^{*} \mathbf{\mathsf{X}}. \end{align*} $$

These matrices reflect the alignment of the test matrix $\mathbf {\mathsf {X}}$ with the matrix $\mathbf {\mathsf {V}}_1$ of dominant right singular vectors of $\mathbf {\mathsf {B}}$ . We assume ${\mathbf {\mathsf {X}}}_1$ has full row rank.

With this notation, we can present a strong deterministic bound on the error in Algorithm 7.

Theorem 11.4 (rangefinder: deterministic bound)

Let $\mathbf {\mathsf {Y}} = \mathbf {\mathsf {BX}}$ be the sample matrix obtained by testing $\mathbf {\mathsf {B}}$ with $\mathbf {\mathsf {X}}$ . With the notation and assumptions above, for all $k \leq \ell $ ,

(11.2) $$ \begin{align} {\Vert { ({\mathbf{\mathsf{I}}} - \mathbf{\mathsf{P}}_{\mathbf{\mathsf{Y}}}) \mathbf{\mathsf{B}} } \Vert} \leq \sigma_{k+1} + {\Vert { \mathbf{\mathsf{\Sigma}}_{2} {\mathbf{\mathsf{X}}}_{2} {\mathbf{\mathsf{X}}}_1^{\dagger} } \Vert}. \end{align} $$

A related inequality holds for every quadratic unitarily invariant norm.

The result and its proof are drawn from Halko et al. (Reference Halko, Martinsson and Tropp2011a , Theorem 9.1). The same bound was obtained independently in Boutsidis, Mahoney and Drineas (Reference Mahoney and Drineas2009) by means of a different technique.

Theorem 11.4 leads to sharp bounds on the performance of the rangefinder in most situations of practical interest. Let us present a sketch of the argument. Our approach can be modified to obtain matching lower and upper bounds, but they do not give any additional insight into the performance.

Proof. In view of Proposition 11.1, we want to bound the spectral norm of

$$ \begin{align*} {\vert {({\mathbf{\mathsf{I}}} - \mathbf{\mathsf{P}}_{\mathbf{\mathsf{Y}}}) \mathbf{\mathsf{B}}} \vert}^2 = (\mathbf{\mathsf{B}}^{*}\mathbf{\mathsf{B}}) / \mathbf{\mathsf{X}}. \end{align*} $$

First, change coordinates so that the right singular vectors of $\mathbf {\mathsf {B}}$ are the identity: $\mathbf {\mathsf {V}} = {\mathbf {\mathsf {I}}}$ . In particular, $\mathbf {\mathsf {B}}^{*} \mathbf {\mathsf {B}} = \mathbf {\mathsf{\Sigma }}^2$ is diagonal. By homogeneity of (11.2), we may assume that $\sigma _1 = 1$ . Next, using Corollary 11.2, we may also assume that $\sigma _1 = \cdots = \sigma _k = 1$ , which leads to the worst-case error. Thus it suffices to bound the spectral norm of the PSD matrix

$$ \begin{align*} \mathbf{\mathsf{S}} := \begin{bmatrix} {\mathbf{\mathsf{I}}}_k & \mathbf{\mathsf{0}} \\ \mathbf{\mathsf{0}} & \mathbf{\mathsf{\Sigma}}_{2}^2\end{bmatrix} / \mathbf{\mathsf{X}}. \end{align*} $$

To accomplish this task, we may as well take the Schur complement of the diagonal matrix with respect to a test matrix that has a smaller range than $\mathbf {\mathsf {X}}$ ; see Ando (Reference Ando and Zhang2005, Theorem 5.9(iv)). Define $\widetilde {\mathbf {\mathsf {X}}} := \mathbf {\mathsf {X}} {\mathbf {\mathsf {X}}}_{1}^{\dagger } = [{\mathbf {\mathsf {I}}}_k; {\mathbf {\mathsf {X}}}_{2}{\mathbf {\mathsf {X}}}_1^{\dagger }]$ . Since $\operatorname {range}(\widetilde {\mathbf {\mathsf {X}}}) \subseteq \operatorname {range}(\mathbf {\mathsf {X}})$ ,

Using the definition of the Schur complement, we can write out the matrix on the right-hand side in block form. With the abbreviation $\mathbf {\mathsf {F}} := \mathbf {\mathsf{\Sigma }}_{2} {\mathbf {\mathsf {X}}}_{2} {\mathbf {\mathsf {X}}}_1^{\dagger }$ ,

$$ \begin{align*} \widetilde{\mathbf{\mathsf{S}}} = \begin{bmatrix} {\mathbf{\mathsf{I}}} - {({\mathbf{\mathsf{I}}} + \mathbf{\mathsf{FF}}^{*})}^{-1} & \star \\ \star & {\mathbf{\mathsf{\Sigma}}}_{2}^2 - \mathbf{\mathsf{F}}{({\mathbf{\mathsf{I}}} + \mathbf{\mathsf{FF}}^{*})}^{-1} {\mathbf{\mathsf{F}}}^{*} \end{bmatrix}\!. \end{align*} $$

The $\star $ symbol denotes matrices that do not play a role in the rest of the argument. We can bound the block matrix above in the PSD order:

The inequality for the top-left block holds because $1 - (1+a)^{-1} \leq a$ for all numbers $a \geq 0$ . Last, take the spectral norm:

$$ \begin{align*} {\Vert { \mathbf{\mathsf{S}} } \Vert} \leq {\Vert { \widetilde{\mathbf{\mathsf{S}}} } \Vert} \leq {\biggl\Vert { \begin{bmatrix} \mathbf{\mathsf{FF}}^{*} & \star \\ \star & {\mathbf{\mathsf{\Sigma}}^2_{2}} \end{bmatrix} } \biggr\Vert} \leq {\Vert {{\mathbf{\mathsf{FF}}}^{*}} \Vert} + {\Vert {{\mathbf{\mathsf{\Sigma}}^2_{2}}} \Vert}. \end{align*} $$

This bound is stronger than the stated result. □

11.4.2. Gaussian test matrices

We can obtain precise results for the behaviour of the randomized rangefinder when the test matrix is (real) standard normal. Let us present a variant of Halko et al. (Reference Halko, Martinsson and Tropp2011a , Theorem 10.1).

Theorem 11.5 (rangefinder: Gaussian analysis)

Fix a matrix $\mathbf {\mathsf {B}} \in {\mathbb {R}}^{m \times n}$ with singular values $\sigma _1 \geq \sigma _2 \geq \cdots .$ Draw a standard normal test matrix $\mathbf {\mathsf{\Omega }} \in \mathbb {F}^{n \times \ell }$ , and construct the sample matrix $\mathbf {\mathsf {Y}} = \mathbf {\mathsf {B}}\mathbf {\mathsf{\Omega }}$ . Choose $k < \ell - 1$ , and introduce the random variable

$$ \begin{align*} Z = {\Vert { \mathbf{\mathsf{\Gamma}}^{\dagger} } \Vert} \quad\text{where }\mathbf{\mathsf{\Gamma}} \in {\mathbb{R}}^{k \times \ell}\text{ is standard normal.} \end{align*} $$

Then the expected error in the random rangefinder satisfies

$$ \begin{align*} \operatorname{\mathbb{E}} {\Vert {({\mathbf{\mathsf{I}}} - \mathbf{\mathsf{P}}_{\mathbf{\mathsf{Y}}}) \mathbf{\mathsf{B}}} \Vert} \leq \biggl( 1 + \sqrt{\displaystyle \frac{k}{\ell-k-1}} \biggr) \sigma_{k+1} + (\operatorname{\mathbb{E}} Z) \Bigl(\sum\nolimits_{j> k} \sigma_{j}^2\Bigr)^{1/2}. \end{align*} $$

In other words, the randomized rangefinder computes an $\ell $ -dimensional subspace that captures as much of the action of the matrix $\mathbf {\mathsf {B}}$ as the best k-dimensional subspace. If we think about k as fixed and $\ell $ as the variable, we only need to choose $\ell $ slightly larger than k to enjoy this outcome.

The error is comparable with $\sigma _{k+1}$ , the error in the best rank-k approximation, provided that the tail singular values $\sigma _{j}$ for $j> k$ have small $\ell _2$ norm. This situation occurs, for example, when $\mathbf {\mathsf {B}}$ has a rapidly decaying spectrum.

Proof. Here is a sketch of the argument. Since the test matrix $\mathbf {\mathsf{\Omega }}$ is standard normal, the matrices $\mathbf {\mathsf{\Omega }}_{1} := \mathbf {\mathsf {V}}_1^{*} \mathbf {\mathsf{\Omega }}$ and $\mathbf {\mathsf{\Omega }}_{2} := \mathbf {\mathsf {V}}_2^{*}\mathbf {\mathsf{\Omega }}$ are independent standard normal matrices because $\mathbf {\mathsf {V}}_1$ and $\mathbf {\mathsf {V}}_2$ are orthonormal and mutually orthogonal. Using Chevet’s theorem (Halko et al. Reference Halko, Martinsson and Tropp2011a , Proposition 10.1),

$$ \begin{align*} \operatorname{\mathbb{E}} \Vert \mathbf{\mathsf{\Sigma}}_{2} {\mathbf{\mathsf{\Omega}}_{2}} {\mathbf{\mathsf{\Omega}}^{\dagger}_1 \Vert} &= \operatorname{\mathbb{E}}_{{\mathbf{\mathsf{\Omega}}_{1}}} \operatorname{\mathbb{E}}_{{\mathbf{\mathsf{\Omega}}_{2}}} \big[ {\Vert { \mathbf{\mathsf{\Sigma}}_{2} {\mathbf{\mathsf{\Omega}}_{2}} {\mathbf{\mathsf{\Omega}}^{\dagger}_1} } \Vert} \big] \\* &\leq \operatorname{\mathbb{E}} \big[ {\Vert {\mathbf{\mathsf{\Sigma}}_{2}} \Vert} {{\Vert {{\mathbf{\mathsf{\Omega}}_1^{\dagger}}} \Vert}_{\mathrm{F}}} + {{\Vert {{\mathbf{\mathsf{\Sigma}}_1}} \Vert}_{\mathrm{F}}} {\Vert {\mathbf{\mathsf{\Omega}}_1^{\dagger}} \Vert} \big] \\* &\leq \sqrt{\displaystyle \frac{k}{\ell-k-1}} {\Vert { \mathbf{\mathsf{\Sigma}}_{2} } \Vert} + (\operatorname{\mathbb{E}} Z) {{\Vert {{\mathbf{\mathsf{\Sigma}}_{2}}} \Vert}_{\mathrm{F}}}. \end{align*} $$

The last inequality involves a well-known estimate for the trace of an inverted Wishart matrix (Halko et al. Reference Halko, Martinsson and Tropp2011a , Proposition 10.2). □

To make use of the result, we simply insert estimates for the expectation of the random variable Z. For instance,

$$ \begin{align*} \operatorname{\mathbb{E}} Z \leq \displaystyle \frac{\mathrm{e}\sqrt{\ell}}{\ell - k} \quad\text{when }2 \leq k < \ell\quad \text{and}\quad \operatorname{\mathbb{E}} Z \approx \displaystyle \frac{1}{\sqrt{\ell} - \sqrt{k}} \quad\text{for }k \ll \ell. \end{align*} $$

These estimates lead to very accurate performance bounds across a wide selection of matrices and parameters.

The rangefinder also operates in the regime $k \in \{ \ell - 1, \ell \}$ . In this case it attains significantly larger errors. A heuristic is

This point follows because ${{\Vert {{\mathbf {\mathsf{\Omega }}_1^{\dagger }}} \Vert }_{\mathrm {F}}} \approx k$ and ${\Vert {\mathbf {\mathsf{\Omega }}_1^{\dagger }} \Vert } \approx \sqrt {k}$ when $k \approx \ell $ .

For relevant results about Gaussian random matrices, we refer the reader to Edelman (Reference Edelman1989), Davidson and Szarek (Reference Davidson and Szarek2001), Chen and Dongarra (Reference Chen and Dongarra2005), Bai and Silverstein (Reference Bai and Silverstein2010) and Halko et al. (Reference Halko, Martinsson and Tropp2011a ).

11.5. Other test matrices

In many cases it is too expensive to use Gaussian test matrices to implement Algorithm 7. Instead, we may prefer to apply (the adjoint of) one of the structured random embeddings discussed in Section 9.