Faster Linear Algebra Algorithms with Structured Random Matrices

TL;DR

This paper introduces the OSI property to analyze structured random matrices, enabling faster randomized linear algebra algorithms with theoretical guarantees and practical efficiency for low-rank approximation and regression tasks.

Contribution

It develops a new OSI-based framework for analyzing structured random matrices, proving their effectiveness in randomized algorithms and identifying many practical matrix models.

Findings

Structured random matrices satisfy the OSI property.

Theoretical analysis shows near-optimal runtimes for key linear algebra tasks.

Empirical results demonstrate high performance of structured matrices in applications.

Abstract

To achieve the greatest possible speed, practitioners regularly implement randomized algorithms for low-rank approximation and least-squares regression with structured dimension reduction maps. Despite significant research effort, basic questions remain about the design and analysis of randomized linear algebra algorithms that employ structured random matrices. This paper develops a new perspective on structured dimension reduction, based on the oblivious subspace injection (OSI) property. The OSI property is a relatively weak assumption on a random matrix that holds when the matrix preserves the length of vectors on average and, with high probability, does not annihilate any vector in a low-dimensional subspace. With the OSI abstraction, the analysis of a randomized linear algebra algorithm factors into two parts: (i) proving that the algorithm works when implemented with an OSI; and (ii) proving that a given random matrix model has the OSI property. This paper develops both parts of the program. First, it analyzes standard randomized algorithms for low-rank approximation and least-squares regression under the OSI assumption. Second, it identifies many examples of OSIs, including random sparse matrices, randomized trigonometric transforms, and random matrices with tensor product structure. These theoretical results imply faster, near-optimal runtimes for several fundamental linear algebra tasks. The paper also provides guidance on implementation, along with empirical evidence that structured random matrices offer exemplary performance for a range of synthetic problems and contemporary scientific applications.

Tables

Table 1. Table 1: Sketching algorithms for linear algebra: Our runtimes. Sketching time and processing time for the algorithms in Theorems ˜ 1.3 and  1.4 applied to input matrix 𝑨 ∈ \mathbb ​ F n × d \bm{A}\in\mathbb{F}^{n\times d} and, for linear regression, response matrix 𝑩 ∈ \mathbb ​ F n × m \bm{B}\in\mathbb{F}^{n\times m} . The test matrices 𝛀 ∈ \mathbb ​ F d × k \bm{\Omega}\in\mathbb{F}^{d\times k} and 𝚿 ∈ \mathbb ​ F n × p \bm{\Psi}\in\mathbb{F}^{n\times p} . The function 𝒯 sk ​ ( 𝑴 , q ) {\mathcal{T}}_{\rm sk}(\bm{M},q) returns the cost of forming the sketch 𝑴 ​ 𝛀 \bm{M}\bm{\Omega} for a test matrix 𝛀 \bm{\Omega} with q q columns, and 𝒯 mv ​ ( 𝑴 ) {\mathcal{T}}_{\rm mv}(\bm{M}) is the time to compute the matrix–vector product 𝒙 ↦ 𝑴 ​ 𝒙 \bm{x}\mapsto\bm{M}\bm{x} . The last column lists the runtimes implied by our analysis when 𝑨 \bm{A} and 𝑩 \bm{B} are dense arrays and the test matrices are SparseStacks ( Definition ˜ 1.7 ) with row sparsity ζ = 𝒪 ​ ( log ⁡ r ) \zeta={\mathcal{O}}(\log r) and with sketch dimensions k , p = 𝒪 ​ ( r ) k,p={\mathcal{O}}(r) for low-rank approximation and with p = 𝒪 ​ ( d ) p={\mathcal{O}}(d) for linear regression. For legibility, the order notation 𝒪 {\mathcal{O}} is suppressed.

Task	Algorithm	Sketching time	Processing time	Total (dense data, SparseStack)
Low-rank	RSVD	${𝒯}_{\mathrm{sk}}(𝑨,k)$	$n{k}^2+k{𝒯}_{\mathrm{mv}}\left({𝑨}^{\ast }\right)$	$ndr$
approximation	Nyström (psd $𝑨$)	${𝒯}_{\mathrm{sk}}(𝑨,k)$	$n{k}^2$	$n^2\log \left(r\right)+n{r}^2$
	Gen. Nyström	${𝒯}_{\mathrm{sk}}(𝑨,k)+{𝒯}_{\mathrm{sk}}({𝑨}^{\ast },p)$	$n{k}^2+dkp$	$nd\log \left(r\right)+\left(n+d\right)r^2$
Least squares	Sketch-and-solve	${𝒯}_{\mathrm{sk}}({𝑨}^{\ast },p)$ + ${𝒯}_{\mathrm{sk}}({𝑩}^{\ast },p)$	$p{d}^2+pdm$	$nd\log \left(d\right)+d^2\left(d+m\right)$

Table 2. Table 2: Sketching algorithms for sparse linear algebra: Runtimes. For sketching algorithms with the SparseStack test matrix ( Definition ˜ 1.7 ), we compare our results with the existing analyses of Cohen  [ 36 ] and of Chenakkod et al.  [ 32 ] . Imported Theorem ˜ 1.8 justifies a row sparsity ζ = 𝒪 ​ ( log ⁡ d ) \zeta={\mathcal{O}}(\log d) at an embedding dimension k = 𝒪 ​ ( r ) k={\mathcal{O}}(r) , so we obtain both the faster sketching time of Cohen and the faster postprocessing time of Chenakkod et al. Matrices 𝑨 ∈ \mathbb ​ F n × d \bm{A}\in\mathbb{F}^{n\times d} and 𝑩 ∈ \mathbb ​ F n × d \bm{B}\in\mathbb{F}^{n\times d} can be dense or sparse, nnz ⁡ ( ⋅ ) \operatorname{nnz}(\cdot) returns the number of nonzero entries, and we employ the standard matrix multiplication algorithm. Blue terms arise from sketching costs, while magenta terms arise from postprocessing. We omit 𝒪 {\mathcal{O}} symbols for legibility. See Table ˜ 3 for another comparison.

	Best available runtime
Algorithm	Chenakkod et al. [32]	Cohen [36]	Ours
Sketch & solve (Alg. 3)	$\mathrm{nnz}\left(\left[𝑨𝑩\right]\right){\log }^3\left(d\right)+d^3$	$\mathrm{nnz}\left(\left[𝑨𝑩\right]\right)\log \left(d\right)+d^3\log d$	$\mathrm{nnz}\left(\left[𝑨𝑩\right]\right)\log \left(d\right)+d^3$
RSVD (Alg. 1)	$\mathrm{nnz}\left(𝑨\right)r+n{r}^2$	$\mathrm{nnz}\left(𝑨\right)r+n{r}^2{\log }^{2}r$	$\mathrm{nnz}\left(𝑨\right)r+n{r}^2$
Nyström (psd $𝑨$, Alg. 2)	$\mathrm{nnz}\left(𝑨\right){\log }^3\left(r\right)+n{r}^2$	$\mathrm{nnz}\left(𝑨\right)\log \left(r\right)+n{r}^2{\log }^{2}r$	$\mathrm{nnz}\left(𝑨\right)\log \left(r\right)+n{r}^2$
Gen. Nyström (Alg. 4)	$\mathrm{nnz}\left(𝑨\right){\log }^3\left(r\right)+\left(n+d\right)r^2$	$\mathrm{nnz}\left(𝑨\right)\log \left(r\right)+\left(n+d\right)r^2{\log }^{3}r$	$\mathrm{nnz}\left(𝑨\right)\log \left(r\right)+\left(n+d\right)r^2$

Table 3. Table 3: Sparse test matrices: Comparison. Theoretical results for sparse test matrices, including the type of guarantee (OSE or OSI), the embedding dimension k k and row sparsity ζ \zeta for sketching an r r -dimensional subspace, and the runtime for generalized Nyström ( Algorithm ˜ 4 ). OSI guarantees for the SparseStack test matrix ( Definition ˜ 1.7 ) permit both the smaller row sparsity ζ = 𝒪 ​ ( log ⁡ r ) \zeta={\mathcal{O}}(\log r) of Cohen’s result [ 36 ] and the smaller sketch dimension k = 𝒪 ​ ( r ) k={\mathcal{O}}(r) of Chenakkod et al. [ 32 ] .

Guarantee	Test matrix	Embedding dim. $k$	Row sparsity $\zeta$	Gen. Nyström runtime
OSE	CountSketch [78, 90]	$𝒪\left(r^2\right)$	$1$	$𝒪\left(\mathrm{nnz}\left(𝑨\right)+n{r}^4\right)$
	SparseStack [36]	$𝒪\left(r\log r\right)$	$𝒪\left(\log r\right)$	$𝒪\left(\mathrm{nnz}\left(𝑨\right)\log \left(r\right)+n{r}^2{\log }^{3}r\right)$
	SparseStack [33]	$𝒪\left(r\right)$	$𝒪\left({\log }^{3}r\right)$	$𝒪\left(\mathrm{nnz}\left(𝑨\right){\log }^3\left(r\right)+n{r}^2\right)$
OSI	SparseIID [115]	$𝒪\left(r\right)$	$𝒪\left(\log r\right)$	$𝒪\left(\mathrm{nnz}\left(𝑨\right)\log \left(r\right)+n{r}^2\right)$
	SparseStack (Thm. 1.8)	$𝒪\left(r\right)$	$𝒪\left(\log r\right)$	$𝒪\left(\mathrm{nnz}\left(𝑨\right)\log \left(r\right)+n{r}^2\right)$
None	SparseStack (Practice)	$2r$	$4$	$𝒪\left(\mathrm{nnz}\left(𝑨\right)+n{r}^2\right)$

Table 4. Table 4: Fast trigonometric transform test matrices: Comparison. Choices for the sparse matrix 𝑺 \bm{S} used within the trigonometric transform test matrix 𝛀 ∈ \mathbb ​ F d × k \bm{\Omega}\in\mathbb{F}^{d\times k} , described in  ˜ 4.1 . We report the embedding dimension k k and sparsity nnz ⁡ ( 𝑺 ) \operatorname{nnz}(\bm{S}) that suffice to justify the ( r , α ) (r,\alpha) -OSI property. The last column lists the runtimes for the sketch-and-solve method ( Algorithm ˜ 3 ) for an input matrix 𝑨 ∈ \mathbb ​ F n × d \bm{A}\in\mathbb{F}^{n\times d} and response matrix 𝑩 ∈ \mathbb ​ F n × d \bm{B}\in\mathbb{F}^{n\times d} . The SparseCol test matrix ( Definition ˜ 4.1 ), used in our SparseRTT construction ( Definition ˜ 1.10 ), simultaneously achieves the minimum embedding dimension k k and sparsity nnz ⁡ ( 𝑺 ) \operatorname{nnz}(\bm{S}) , resulting in improved runtimes for sketch-and-solve.

Sparse matrix $𝑺$	Embedding dim. $k$	Sparsity $\mathrm{nnz}\left(𝑺\right)$	Sketch-and-solve time
Subsampling [4]	$𝒪\left(r\log r\right)$	$𝒪\left(r\log r\right)$	$𝒪\left(nd\log \left(d\right)+d^3\log d\right)$
CountSketch [26]	$𝒪\left(r\right)$	$d$	$𝒪\left(nd\log \left(n\right)+d^3\right)$
LESS [33, 32]	$𝒪\left(r\right)$	$𝒪\left(r{\log }^{3}r\right)$	$𝒪\left(nd{\log }^3\left(d\right)+d^3\right)$
SparseCol (Def. 4.1)	$𝒪\left(r\right)$	$𝒪\left(r\log r\right)$	$𝒪\left(nd\log \left(d\right)+d^3\right)$

Table 5. Table 5 : Base distributions: Fourth-moment bounds, [ 79 , Tab. 1] . The exact fourth-moment constant C 𝝂 \mathrm{C}_{\boldsymbol{\nu}} for several base distributions 𝝂 ∈ \mathbb ​ F d 0 \boldsymbol{\nu}\in\mathbb{F}^{d_{0}} , including standard normal distributions, Rademacher distributions, the Steinhaus distribution, and the uniform distribution on the sphere of radius d 0 \sqrt{d_{0}} . The last two columns list the embedding dimension k k that suffices to obtain an ( r , 1 / 2 ) (r,\nicefrac{{1}}{{2}}) -OSI, where the 𝒪 {\mathcal{O}} is omitted for legibility. See  Theorem ˜ 5.1 for the definition of C 𝝂 \mathrm{C}_{\boldsymbol{\nu}} and Section ˜ 5.1 for the definitions of the distributions.

			Embedding dimension $k=𝒪\left(\cdots \right)$
Field ($\text{mathbb}F$)	Base distribution ($𝝂$)	$C_{𝝂}$	General $d_0$	$d_0=2$
$\text{mathbb}R$	$𝒩(\mathrm{𝟎},𝐈)$	$3$	$3^{\ell }r$	$3^{\ell }r$
	Real Rademacher	$3-\frac{2}{d_0}$	${\left(3-\frac{2}{d_0}\right)}^{\ell }r$	$2^{\ell }r$
	Real spherical	$3-\frac{6}{d_0+2}$	${\left(3-\frac{6}{d_0+2}\right)}^{\ell }r$	${1.5}^{\ell }r$
$\text{mathbb}C$	${𝒩}_{\text{mathbb}C}(\mathrm{𝟎},𝐈)$	$2$	$2^{\ell }r$	$2^{\ell }r$
	Complex Rademacher	$2-\frac{1}{d_0}$	${\left(2-\frac{1}{d_0}\right)}^{\ell }r$	${1.5}^{\ell }r$
	Steinhaus	$2-\frac{1}{d_0}$	${\left(2-\frac{1}{d_0}\right)}^{\ell }r$	${1.5}^{\ell }r$
	Complex spherical	$2-\frac{2}{d_0+1}$	${\left(2-\frac{2}{d_0+1}\right)}^{\ell }r$	$1.{\overline{33}}^{\ell }r$

Table 6. Table 6 : Gaussian matrix models, [ 115 , Secs. 3 and 6.3.2] . Entries with a leading ≤ \leq indicate inequalities; entries without a leading ≤ \leq are exact computations. The symbol μ ​ ( 𝑸 ) \mu(\bm{Q}) denotes the coherence of the matrix 𝑸 ∈ \mathbb ​ F d × r \bm{Q}\in{\mathbb F}^{d\times r} ( Section ˜ 1.8 ). See  Section ˜ 7.3.1 for other definitions.

Description	Notation $\left(𝒁\right)$	$\text{mathbb}E\left[𝒁\right]$	$Var\left[𝒁\right]\left(𝑴\right)$	$\text{mathbb}E\left[{\lambda }_{\max }\left(𝒁\right)\right]$	${\sigma }_{\ast }^2\left(𝒁\right)$
Scalar Gaussian	$g{𝐈}_d$	$\mathrm{𝟎}$	${\left(\mathrm{Tr}\left(𝑴\right)\right)}^2$	$0$	$1$
IID real diagonal	${𝑫}_d$	$\mathrm{𝟎}$	${\sum }_{i=1}^d{\left({𝑴}_{ii}\right)}^2$	$\le \sqrt{2\log d}$	$1$
Compressed diagonal	${𝑸}^{\ast }{𝑫}_d𝑸$	$\mathrm{𝟎}$	${\sum }_{i=1}^r{\left({\left(𝑸𝑴{𝑸}^{\ast }\right)}_{ii}\right)}^2$	$\le \sqrt{2\mu \left(𝑸\right)\log r}$	$\le \mu \left(𝑸\right)$
GUE	${𝑮}_{\mathrm{gue}}^{\left(d\right)}$	$\mathrm{𝟎}$	${\Vert 𝑴\Vert }_F^2$	$2\sqrt{d}$	$1$

Equations

$$\alpha \cdot {\Vert \mathbf{x}\Vert }_2^2\le {\Vert {\mathbf{\Omega }}^{\ast }\mathbf{x}\Vert }_2^2\le \beta \cdot {\Vert \mathbf{x}\Vert }_2^2\text{for all }\mathbf{x}\in \mathcal{V}.$$

$$E{\Vert {\mathbf{\Omega }}^{\ast }\mathbf{x}\Vert }_2^2={\Vert \mathbf{x}\Vert }_2^2\text{for all x∈Fd.}$$

$$\alpha \cdot {\Vert \mathbf{x}\Vert }_2^2\le {\Vert {\mathbf{\Omega }}^{\ast }\mathbf{x}\Vert }_2^2\text{for all }\mathbf{x}\in \mathcal{V}.$$

$${\Vert \mathbf{A}-\hat{\mathbf{A}}\Vert }_{\mathrm{F}}^2\le (\mathrm{C}/\alpha )\cdot \underset{\mathrm{rank}\mathbf{B}\le r}{\min }{\Vert \mathbf{A}-\mathbf{B}\Vert }_{\mathrm{F}}^2.$$

$${\Vert \mathbf{A}-\hat{\mathbf{A}}\Vert }_{\ast }\le (\mathrm{C}/\alpha )\cdot \underset{\mathrm{rank}\mathbf{B}\le r}{\min }{\Vert \mathbf{A}-\mathbf{B}\Vert }_{\ast }.$$

$${\Vert \mathbf{A}-\hat{\mathbf{A}}\Vert }_{\mathrm{F}}^2\le (\mathrm{C}/\alpha {)}^2\cdot \underset{\mathrm{rank}\mathbf{B}\le r}{\min }{\Vert \mathbf{A}-\mathbf{B}\Vert }_{\mathrm{F}}^2.$$

$${\min }_{\mathbf{X}\in {\mathbb{F}}^{d\times m}}{\Vert \mathbf{A}\mathbf{X}-\mathbf{B}\Vert }_{\mathrm{F}}^2.$$

$${\Vert \mathbf{A}\tilde{\mathbf{X}}-\mathbf{B}\Vert }_{\mathrm{F}}^2\le (\mathrm{C}/\alpha )\cdot {\min }_{\mathbf{X}}{\Vert \mathbf{A}\mathbf{X}-\mathbf{B}\Vert }_{\mathrm{F}}^2,$$

$$\mathbf{\Omega }:=\frac{1}{\sqrt{\zeta }}\left[\begin{array}{ccc}{\varrho }_{11}{\mathbf{e}}_{s_{11}}^{\ast }& \cdots & {\varrho }_{1\zeta }{\mathbf{e}}_{s_{1\zeta }}^{\ast }\\ {\varrho }_{21}{\mathbf{e}}_{s_{21}}^{\ast }& \cdots & {\varrho }_{2\zeta }{\mathbf{e}}_{s_{2\zeta }}^{\ast }\\ ⋮& \ddots & ⋮\\ {\varrho }_{d1}{\mathbf{e}}_{s_{d1}}^{\ast }& \cdots & {\varrho }_{d\zeta }{\mathbf{e}}_{s_{d\zeta }}^{\ast }\end{array}\right]\text{where}\begin{array}{rl}& \text{ϱij∼textscrademacher iid;}\\ & \text{sij∼textscuniform{1,…,b} iid.}\end{array}$$

$$\mathbf{\Omega }:=\mathbf{D}\mathbf{F}\mathbf{S}\in {\mathbb{F}}^{d\times k}\text{where}\begin{array}{rl}& \text{D∈Fd×d is a random diagonal matrix;}\\ & \text{F∈Fd×d is a trigonometric transform;}\\ & \text{S∈Fd×k is a sparse random matrix.}\end{array}$$

$$\mathbf{\omega }={\mathbf{\omega }}^{(1)}\otimes \cdots \otimes {\mathbf{\omega }}^{(\ell )}\text{where}{\mathbf{\omega }}^{(1)},\dots ,{\mathbf{\omega }}^{(\ell )}\in {\mathbb{F}}^{d_0}.$$

$$\mathbf{\Omega }:=\frac{1}{\sqrt{k}}\left[\begin{array}{ccc}& & \\ {\mathbf{\omega }}_1& \cdots & {\mathbf{\omega }}_k\\ & & \end{array}\right]\in {\mathbb{F}}^{d_0^{\ell }\times k}\text{where}{\mathbf{\omega }}_j:={\mathbf{\omega }}_j^{(1)}\otimes \cdots \otimes {\mathbf{\omega }}_j^{(\ell )}\in {\mathbb{F}}^{d_0^{\ell }}\text{and}\text{ωj(i)∼ν iid.}$$

$${\Vert \mathbf{B}-\tilde{\mathbf{B}}\Vert }_{\mathrm{F}}^2\le \mathrm{C}\cdot \underset{\mathbf{M}\in \mathcal{F}}{\min }{\Vert \mathbf{B}-\mathbf{M}\Vert }_{\mathrm{F}}^2.$$

$$\mathbf{\Omega }:=\frac{1}{\sqrt{k}}\left[\begin{array}{ccc}{\mathbf{\omega }}_1& \dots & {\mathbf{\omega }}_k\end{array}\right]\in {\mathbb{F}}^{d\times k}\text{where ωj∼ω iid.}$$

$${\Vert {\mathbf{\Omega }}^{\ast }\mathbf{x}\Vert }_2^2=\frac{1}{k}\sum _{i=1}^k|⟨{\mathbf{\omega }}_i,\mathbf{x}⟩{|}^2\text{for each fixed vector x∈Fd.}$$

$${\sigma }_{\min }^2({\mathbf{\Omega }}^{\ast }\mathbf{Q})={\lambda }_{\min }({\mathbf{Q}}^{\ast }\mathbf{\Omega }{\mathbf{\Omega }}^{\ast }\mathbf{Q})\ge \alpha >0\text{with probability at least nicefrac1920.}$$

$${\sigma }_{\min }(\mathbf{B}):=\underset{{\Vert \mathbf{u}\Vert }_2=1}{\min }{\Vert \mathbf{B}\mathbf{u}\Vert }_2\text{and}{\lambda }_{\min }(\mathbf{A}):=\underset{{\Vert \mathbf{u}\Vert }_2=1}{\min }{\mathbf{u}}^{\ast }\mathbf{A}\mathbf{u}.$$

$$\operatorname*{\mathbb{E}}\big{[}|\langle\boldsymbol{\omega},\bm{u}\rangle|^{4}\big{]}\leq h\quad\text{for each unit vector $\bm{u}\in\mathbb{F}^{d}$.}$$

$$\mathbb{P}\big{\{}|\langle\boldsymbol{\omega},\bm{u}\rangle|^{2}\leq\tau\big{\}}\leq 0.01\quad\text{for each unit vector $\bm{u}\in\mathbb{F}^{d}$.}$$

$$\mathbf{Y}:={\mathbf{Q}}^{\ast }\mathbf{\Omega }{\mathbf{\Omega }}^{\ast }\mathbf{Q}=\frac{1}{k}\sum _{j=1}^k{\mathbf{W}}_j\text{where the summands Wj=Q∗ωjωj∗Q are iid.}$$

$$\operatorname*{\mathbb{E}}[\bm{Z}]=\operatorname*{\mathbb{E}}[\bm{W}_{1}]\quad\text{and}\quad\operatorname*{Var}[\operatorname{Tr}(\bm{M}\bm{Z})]=k^{-1}\cdot\operatorname*{\mathbb{E}}\big{[}\operatorname{Tr}(\bm{M}\bm{W}_{1})^{2}\big{]}\quad\text{for all self-adjoint $\bm{M}\in\mathbb{F}^{r\times r}$.}$$

$${\Vert \mathbf{B}({\mathbf{Q}}_{\perp }^{\ast }\mathbf{\Omega })({\mathbf{Q}}^{\ast }\mathbf{\Omega }{)}^{†}\Vert }_{\mathrm{F}}^2\le (\mathrm{C}/\alpha )\cdot {\Vert \mathbf{B}\Vert }_{\mathrm{F}}^2.$$

$${\min }_{\mathbf{X}}{\Vert {\mathbf{\Omega }}^{\ast }(\mathbf{Q}\mathbf{X}-{\mathbf{Q}}_{\perp }{\mathbf{B}}^{\ast })\Vert }_{\mathrm{F}}^2.$$

$${\Vert \mathbf{B}({\mathbf{Q}}_{\perp }^{\ast }\mathbf{\Omega })({\mathbf{Q}}^{\ast }\mathbf{\Omega }{)}^{†}\Vert }_{\mathrm{F}}^2\le {\Vert \mathbf{B}{\mathbf{Q}}_{\perp }^{\ast }\mathbf{\Omega }\Vert }_{\mathrm{F}}^2\cdot {\Vert ({\mathbf{Q}}^{\ast }\mathbf{\Omega }{)}^{†}\Vert }_2^2\le \frac{{\Vert \mathbf{B}{\mathbf{Q}}_{\perp }^{\ast }\mathbf{\Omega }\Vert }_{\mathrm{F}}^2}{{\sigma }_{\min }^2({\mathbf{\Omega }}^{\ast }\mathbf{Q})}.$$

$$\operatorname*{\mathbb{E}}{}\bigl{[}\mathinner{\!\left\lVert\bm{B}\bm{Q}_{\perp}^{\smash{*}}\bm{\Omega}\right\rVert}_{\rm F}^{2}\bigr{]}=\operatorname*{\mathbb{E}}\big{[}\operatorname{Tr}(\bm{B}\bm{Q}_{\perp}^{\smash{*}}\bm{\Omega}\bm{\Omega}^{\smash{*}}\bm{Q}_{\perp}\bm{B}^{\smash{*}})\big{]}=\operatorname{Tr}(\bm{B}\bm{Q}_{\perp}^{\smash{*}}\bm{Q}_{\perp}\bm{B}^{\smash{*}})=\mathinner{\!\left\lVert\bm{B}\right\rVert}_{\rm F}^{2}.$$

$${\Vert \mathbf{B}{\mathbf{Q}}_{\perp }^{\ast }\mathbf{\Omega }\Vert }_{\mathrm{F}}^2\le 20\cdot {\Vert \mathbf{B}\Vert }_{\mathrm{F}}^2\text{with probability at least 1−201}.$$

$$\mathbf{Q}:=\mathrm{orth}(\mathbf{A}\mathbf{\Omega })\in {\mathbb{F}}^{n\times k^{\prime }}\text{where k′≤k.}$$

$$\hat{\mathbf{A}}:=\mathbf{Q}\mathbf{B}\in {\mathbb{F}}^{n\times d}\text{for }\mathbf{B}:={\mathbf{Q}}^{\ast }\mathbf{A}.$$

$${\Vert \mathbf{A}-\hat{\mathbf{A}}\Vert }_{\mathrm{F}}^2\le (\mathrm{C}/\alpha )\cdot {\Vert \mathbf{A}-[[\mathbf{A}{]]}_r\Vert }_{\mathrm{F}}^2.$$

$$\mathbf{A}=\mathbf{U}\left[\begin{array}{cc}{\mathbf{\Sigma }}_1& \\ & {\mathbf{\Sigma }}_2\end{array}\right]\left[\begin{array}{c}{\mathbf{V}}_1^{\ast }\\ {\mathbf{V}}_2^{\ast }\end{array}\right],$$

Full text

This paper develops a new perspective on structured dimension reduction, based on the oblivious subspace injection (OSI) property. The OSI property is a relatively weak assumption on a random matrix that holds when the matrix preserves the length of vectors on average and, with high probability, does not annihilate any vector in a low-dimensional subspace. With the OSI abstraction, the analysis of a randomized linear algebra algorithm factors into two parts: (i) proving that the algorithm works when implemented with an OSI; and (ii) proving that a given random matrix model has the OSI property.

This paper develops both parts of the program. First, it analyzes standard randomized algorithms for low-rank approximation and least-squares regression under the OSI assumption. Second, it identifies many examples of OSIs, including random sparse matrices, randomized trigonometric transforms, and random matrices with tensor product structure. These theoretical results imply faster, near-optimal runtimes for several fundamental linear algebra tasks. The paper also provides guidance on implementation, along with empirical evidence that structured random matrices offer exemplary performance for a range of synthetic problems and contemporary scientific applications.

1. 1 Introduction

2. 1.1 Empirical context: The power of structured random matrices

3. 1.2 Subspace embeddings (OSEs) and subspace injections (OSIs)

4. 1.3 Randomized linear algebra with OSIs

5. 1.4 OSIs suffice for (some) matrix computations

6. 1.5 A smörgåsbord of OSI constructions

7. 1.6 How to establish the OSI property

8. 1.7 Outline

9. 1.8 Notation & terminology

10. 2 Randomized matrix computations via OSIs

11. 2.1 Basic technical observation: OSIs respect orthogonality

12. 2.2 The randomized SVD

13. 2.3 Nyström approximation

14. 2.4 Sketch-and-solve for least-squares regression

15. 2.5 Generalized Nyström approximation

16. 3 Sparse test matrices

17. 3.1 Randomized linear algebra with SparseStacks

18. 3.2 Alternative sparse test matrices

19. 3.3 Numerical evidence

20. 4 SparseRTT: An optimized fast trigonometric transform test matrix

21. 4.1 The SparseRTT construction and its analysis

22. 4.2 Arithmetic costs

23. 4.3 Implementation guidance

24. 4.4 Alternative constructions

25. 4.5 Proof of Theorem˜4.2

26. 5 Khatri–Rao test matrices and their applications

27. 5.1 Khatri–Rao test matrices

28. 5.2 Applications

29. 5.3 Overwhelming orthogonality: When Khatri–Rao goes bad

30. 5.4 First result: Khatri–Rao embeddings with constant injectivity

31. 5.5 Second result: Khatri–Rao embeddings with proportional embedding dimension

32. 5.6 Application: Approximate matrix recovery from bilinear queries

33. 6 Scientific applications

34. 6.1 Application 1: Scientific data compression

35. 6.2 Application 2: Approximating the partition function of a quantum system

36. 7 The OSI property: Proof strategies

37. 7.1 Methodology 1: Fourth-moment bounds

38. 7.2 Methodology 2: Small-ball probabilities

39. 7.3 Methodology 3: Comparison with a Gaussian matrix

40. 7.4 Methodology 4: Matrix concentration

41. A Analysis of the SparseCol test matrix

42. B Analysis of the SparseStack test matrix

43. B.1 Moments of the CountSketch matrix

44. C Proof of Theorem˜7.3

45. C.1 Refined VC dimension bounds

46. C.2 Analysis in the complex field

47. D Implementation of scientific data compression experiment

48. E Implementation of partition function estimation experiment

1 Introduction

Randomized algorithms allow us to solve certain linear algebra problems faster and more reliably than ever before, with real implications for the practice of machine learning and scientific computing.

In the field of randomized linear algebra, many fundamental algorithms are based on randomized linear dimension reduction, often called sketching [105, 58, 123, 44, 73, 84, 65]. In this context, sketching algorithms compress a large input matrix $\bm{A}\in{\mathbb R}^{n\times d}$ by forming the product $\bm{Y}=\bm{A}\bm{\Omega}\in{\mathbb R}^{n\times k}$ with a smaller random test matrix $\bm{\Omega}\in{\mathbb R}^{d\times k}$, where the embedding dimension $k\ll d$. We can quickly extract salient information about the large input matrix $\bm{A}$ by manipulating the small sketch matrix $\bm{Y}$. This strategy leads to powerful algorithms for several important linear algebra problems, including low-rank approximation [51, 121, 42, 34, 58, 35, 87] and least-squares regression [105, 102, 11].

For many sketching algorithms, the cost of computing the sketch $\bm{Y}=\bm{A}\bm{\Omega}$ dominates the overall runtime. If we use an unstructured test matrix $\bm{\Omega}$, such as the standard Gaussian test matrix in [58], the sketching step requires an expensive dense matrix–matrix product. We can reduce the runtime by using a structured test matrix $\bm{\Omega}$ that supplies a faster algorithm for the matrix–matrix product. The most popular examples are sparse random matrices [2, 34, 90, 63] and randomized trigonometric transforms [4, 124, 45, 111]. For linear algebra problems with tensor structure, we can employ test matrices with a compatible tensor structure, such as random Khatri–Rao products [19, 16, 99, 24, 25, 3].

According to the cognoscenti [64, 73], sketching with a Gaussian test matrix is the gold standard for quality of output. Moreover, we can derive precise theoretical results that predict the empirical performance of Gaussian sketching methods [58, 118, 117, 119, 73, 40, 15]. In practice, algorithms that sketch with a structured test matrix produce outputs whose quality matches the Gaussian standard—at a dramatically lower computational cost [119, 46, 41, 11, 31, 84, 73]. Yet we cannot fully explain the fabulous performance of structured sketching methods through existing theoretical frameworks.

Theory.

This paper develops a new theoretical perspective on sketching with structured random matrices. We prove that several randomized linear algebra algorithms succeed when implemented with any test matrix that satisfies the oblivious subspace injection (OSI) property. The OSI condition is genuinely weaker than the oblivious subspace embedding (OSE) property that is standard in the literature. We explain how to establish the OSI property using modern tools from high-dimensional probability [92, 67, 114, 115, 120]. See Sections˜1.2, 1.4, and 1.6 for an overview of this work.

Design.

We make recommendations on how to design and implement structured sketching methods to attain the best empirical performance, with provable guarantees. For sparse test matrices, our analysis justifies sparsity levels that were previously unsupportable, and we can contemplate levels that had been considered impossible (Section˜1.5.1). Second, we construct a randomized trigonometric transform with the minimal embedding dimension and the fastest possible sketching time (Section˜1.5.2). Last, we produce the first qualitatively accurate bounds on the worst-case behavior of Khatri–Rao sketching methods, and we identify variations that offer exponential improvements over the most popular constructions (Section˜1.5.3).

Why it matters.

For practical computation, we must insist on algorithms that produce outputs with the highest quality, as quickly as possible, and with reliability guarantees. Sketching methods implemented with our recommended test matrices can meet these stringent requirements. We include two stylized applications that illustrate the benefits for scientific computing (Section˜6). First, we use sparse test matrices to compress a large PDE simulation quickly and accurately. Second, we use tensor product test matrices to estimate the partition function of a quantum mechanical system to high precision.

In sum, we hope that our work provides theoreticians with new insights on sketching methods for randomized linear algebra and that it equips practitioners with faster and more reliable algorithms.

1.1 Empirical context: The power of structured random matrices

Let us motivate our investigation by means of an empirical study. Our experiments align with the conclusions of two decades of empirical work on sketching algorithms, including [119, 46, 41, 11, 31, 84]. Structured test matrices provide the same accuracy as Gaussian test matrices, while the structure allows us to design far more efficient algorithms.

To confirm the first point, Figure˜1 charts the approximation performance of the randomized SVD algorithm (RSVD, Algorithm˜1) [58, p. 227], implemented with several structured test matrices. The testbed consists of each matrix from the SuiteSparse collection [38] whose dimensions both fall between $300$ and $500{,}000$; these $2{,}314$ matrices arise from a range of applications and exhibit a variety of singular value distributions. For each structured test matrix and for each input matrix, we compare the rank $k=200$ approximation error to the baseline error achieved with a (scaled) standard Gaussian matrix, and we perform three trials. In most cases, the approximation quality is almost indistinguishable from the Gaussian baseline. Even for the hardest examples, the structured test matrix yields at most $4\times$ higher error.

Even though we suffer no loss in output quality by implementing sketching algorithms with structured test matrices, we gain dramatic improvements in efficiency. For example, Figure˜3 illustrates that sparse test matrices can reduce the sketching time by two orders of magnitude, as compared with Gaussian test matrices. The benefits only accrue as the problem instances become larger.

1.2 Subspace embeddings (OSEs) and subspace injections (OSIs)

To understand the behavior of structured sketching methods for linear algebra, we adopt a perspective that departs from the existing literature. Our assumptions are genuinely weaker, and they open a portal to a more expansive universe of sketching techniques.

In his 2006 paper on sketch-and-solve algorithms for least-squares regression [105], Sarlós introduced the concept of an oblivious subspace embedding (OSE), where “oblivious” means “subspace independent.” He proved that the sketch-and-solve algorithm produces a near-optimal least-squares solution whenever it is implemented with a random test matrix that satisfies the OSE property. Since then, researchers have employed OSEs to analyze other randomized linear algebra algorithms, including the RSVD and the generalized Nyström approximation [123, 44, 73]. In our paper, the scalar field $\mathbb{F}\in\{{\mathbb R},{\mathbb C}\}$.

Definition 1.1** (Oblivious subspace embedding).**

A random matrix $\bm{\Omega}\in\mathbb{F}^{d\times k}$ is called an $(r,\alpha,\beta)$-OSE with subspace dimension $r$, embedding dimension $k\geq r$, injectivity $\alpha\in(0,1]$, and dilation $\beta\geq 1$ if the following condition holds for each fixed $r$-dimensional subspace ${\mathcal{V}}\subseteq\mathbb{F}^{d}$. With probability at least $\frac{19}{20}$,

[TABLE]

As it happens, the two inequalities in the OSE property ˜1.1 play asymmetric roles. The truly essential condition is the lower bound, involving the injection parameter $\alpha$, because it ensures that the random matrix $\bm{\Omega}$ does not annihilate any vector $\bm{x}$ in the fixed subspace ${\mathcal{V}}$. Upper bounds on the dilation $\beta$ are far less important. This insight motivates us to focus on the injectivity property.

Definition 1.2** (Oblivious subspace injection).**

A random matrix $\bm{\Omega}\in\mathbb{F}^{d\times k}$ is called an $(r,\alpha)$-OSI with subspace dimension $r$, embedding dimension $k\geq r$, and injectivity $\alpha\in(0,1]$ when it meets two conditions:

Isotropy. On average, the matrix preserves the squared length of each vector:

[TABLE] 2. 2.

Injectivity. For each fixed $r$-dimensional subspace ${\mathcal{V}}\subseteq\mathbb{F}^{d}$, with probability at least $\frac{19}{20}$,

[TABLE]

Morally, Definition˜1.2 of an OSI is weaker than Definition˜1.1 of an OSE. Indeed, the isotropy property of an OSI ˜1.2 holds on average, while the dilation property of an OSE in ˜1.1 must hold uniformly. We have also obtained numerical evidence that the OSI property is more permissive. In particular, Figure˜2 suggests that an ultra-sparse test matrix can achieve constant injectivity, independent of the subspace dimension $r$, while its dilation grows with $r$. See ˜1.9 for a more precise statement.

In this paper, all test matrices satisfy the isotropy condition ˜1.2. We sometimes refer to an OSI as an embedding, even though its dilation factor $\beta$ is uncontrolled. Section˜1.6 outlines mathematical methods for establishing the OSI property.

1.3 Randomized linear algebra with OSIs

With the OSI abstraction (Definition˜1.2), we can factor the analysis of a randomized sketching algorithm into two independent parts. First, we argue that the algorithm succeeds when it is implemented using any test matrix $\bm{\Omega}\in\mathbb{F}^{d\times k}$ that is an OSI with specific parameters $(r,\alpha)$. Second, we prove that a particular construction of the random test matrix $\bm{\Omega}\in\mathbb{F}^{d\times k}$ enjoys the $(r,\alpha)$-OSI property.

To take full advantage of this theory, we must develop constructions of OSI test matrices that have a favorable computational profile. Here is a summary of the desiderata for an OSI:

OSI parameters.

The embedding dimension is proportional to the subspace dimension: $k={\mathcal{O}}(r)$. The injectivity parameter is constant: $\alpha\geq\mathrm{c}$, where $\mathrm{c}\in(0,1)$ is universal. 2. 2.

Arithmetic.

For each vector $\bm{x}\in\mathbb{F}^{d}$, we can form the sketch $\bm{\Omega}^{{\smash{*}}}\bm{x}$ with $\widetilde{{\mathcal{O}}}(d)$ arithmetic operations. 3. 3.

Implementation.

The OSI construction must be practical and implementable with available program libraries (e.g., sparse arithmetic, fast trig transforms). It should perform well in applications. 4. 4.

*Other resources.

We also hope to control the amount of storage, randomness, and communication required to build and apply the test matrix $\bm{\Omega}$. Our work does not address these questions.

This paper carries out our program to study OSIs. First, we analyze randomized linear algebra algorithms that use OSIs (Section˜1.4). Then we explore the universe of OSI constructions (Section˜1.5).

1.4 OSIs suffice for (some) matrix computations

Several widely used sketching algorithms for randomized linear algebra can be implemented with OSIs. This section summarizes the main conclusions; we postpone full descriptions of the algorithms and the theoretical analysis to Section˜2. Our first result concerns three standard randomized algorithms for computing a low-rank approximation of an input matrix.

Theorem 1.3** **(OSIs for low-rank approximation;

summary of Theorems˜2.2, 2.5, and 2.8).

Fix an input matrix $\bm{A}\in\mathbb{F}^{n\times d}$ and a target rank $r\leq\min\{n,d\}$. Suppose that the random matrix $\bm{\Omega}\in\mathbb{F}^{d\times k}$ is an ($r,\alpha$)-OSI with subspace dimension $r$ and embedding dimension $k\geq r$. Then each of the following statements holds with probability at least $\>\nicefrac{{4}}{{5}}$.

Randomized SVD.* RSVD (Algorithm˜1) with test matrix $\bm{\Omega}$ yields a rank-$k$ approximation $\smash{{\widehat{\smash{\boldsymbol{\bm{A}}}\mathstrut\mkern 4.0mu}\mkern-4.0mu}{}}$ where*

[TABLE] 2. 2.

Nyström approximation.* If $\bm{A}$ is square and positive semidefinite, then the Nyström approximation method (Algorithm˜2) with test matrix $\bm{\Omega}$ returns a rank-$k$ approximation $\smash{{\widehat{\smash{\boldsymbol{\bm{A}}}\mathstrut\mkern 4.0mu}\mkern-4.0mu}{}}$ that satisfies*

[TABLE] 3. 3.

Generalized Nyström approximation.* Let $\bm{\Psi}\in\mathbb{F}^{n\times p}$ be a $(k,\alpha)$-OSI with subspace dimension $k$ and embedding dimension $p\geq k$. With test matrices $\bm{\Omega}$ and $\bm{\Psi}$, the generalized Nyström approximation method (Algorithm˜4) constructs a rank-$k$ approximation $\smash{{\widehat{\smash{\boldsymbol{\bm{A}}}\mathstrut\mkern 4.0mu}\mkern-4.0mu}{}}$ that satisfies*

[TABLE]

Here and elsewhere, $\mathrm{C}\geq 1$ is a universal constant. See Table˜1 for a summary of the arithmetic costs.

We write $\mathinner{\!\left\lVert\cdot\right\rVert}_{\rm F}$ for the Frobenius norm, while $\mathinner{\!\left\lVert\cdot\right\rVert}_{*}$ denotes the Schatten 1-norm. As is standard for randomized low-rank approximation algorithms [58, 73, 123, 55], Theorem˜1.3 presents a bicriteria guarantee where the approximation rank $k$ exceeds the target rank $r$. If desired, we can promote the rank-$k$ approximations to rank-$r$ approximations using a standard truncation method [58, Thm. 9.3].

Our second result treats the sketch-and-solve algorithm [105] for solving a linear regression problem.

Theorem 1.4** (OSIs for least-squares; summary of Theorem˜2.7).**

Fix a design matrix $\bm{A}\in\mathbb{F}^{n\times d}$ and a response matrix $\bm{B}\in\mathbb{F}^{n\times m}$, and consider the linear least-squares problem

[TABLE]

Let $\bm{\Psi}\in\mathbb{F}^{n\times p}$ be a $(d,\alpha)$-OSI with subspace dimension $d$ and embedding dimension $p\geq d$. With probability at least $\nicefrac{{9}}{{10}}$, sketch-and-solve (Algorithm˜3) with test matrix $\bm{\Psi}$ returns an approximate solution $\smash{{\tilde{\smash{\boldsymbol{\bm{X}}}\mathstrut\mkern 4.0mu}\mkern-4.0mu}{}}\in\mathbb{F}^{d\times m}$ that satisfies

[TABLE]

where $\mathrm{C}\geq 1$ is a universal constant. See Table˜1 for a summary of the arithmetic costs.

Together, Theorems˜1.3 and 1.4 prove that we can use OSIs to implement several fundamental sketching algorithms from randomized linear algebra. The runtime analysis (Table˜1) indicates that we can control the total sketching cost by constructing OSIs that support fast matrix–vector products (matvecs), while we can limit postprocessing costs by minimizing the embedding dimension of the OSI. This approach ultimately allows us to obtain faster algorithms because the OSI property is weaker than the familiar OSE property that drives most existing analyses.

A number of prior works have studied the performance of sketching algorithms under conditions less stringent than the OSE property. In particular, the tutorial [44] establishes parts of Theorem˜1.3 under conditions weaker than the OSE. The results [58, Thm. 11.2] and [72, Lem. 2] implicitly justify that RSVD works under an injectivity hypothesis, without control on the dilation. Saibaba & Międlar have undertaken a direct analysis of test matrices with iid columns that satisfy moment conditions [103]. Nevertheless, the OSI condition has not been isolated or systematically exploited.

Remark 1.5* (Constant-factor approximations).*

Theorems˜1.3 and 1.4 demonstrate that several sketching algorithms produce errors that lie within a constant factor of the optimal error. Many papers in the randomized linear algebra literature seek stronger approximation guarantees, within a factor of $1+\varepsilon$ of the optimal error, while the algorithm costs exhibit a $\mathrm{poly}(\varepsilon^{-1})$ dependence in the parameter $\varepsilon>0$. This type of high-accuracy guarantee is presently out of reach of our proof techniques.

Regardless, a constant-factor error bound suffices for typical applications. Indeed, it is most fruitful to deploy sketching algorithms in settings where the optimal error is tiny [88, Sec. 1.1]. To reach high accuracy (i.e., small $\varepsilon$), practitioners do not use simple sketching algorithms. For low-rank approximation, block Krylov methods [86, 116] are preferred. For linear regression, sketch-and-precondition methods [102, 11, 48] are the standard.

1.5 A smörgåsbord of OSI constructions

Section˜1.4 advances the first part of the OSI program, establishing that we can implement several randomized sketching algorithms using OSIs. To address the second part of the program, this section introduces three examples of structured OSIs, based on sparse matrices (Section˜1.5.1), fast trigonometric transforms (Section˜1.5.2), and tensor products (Section˜1.5.3). These constructions meet the desiderata from Section˜1.3, allowing us to perform fast linear algebra computations in several environments.

1.5.0 Baseline: Gaussian test matrices

As the gold standard for output quality, Gaussian sketching methods provide a point of comparison in both the theoretical and the numerical work.

Definition 1.6** (Gaussian test matrix).**

Fix ambient dimension $d\in{\mathbb N}$ and embedding dimension $k\in{\mathbb N}$. A real Gaussian test matrix $\bm{\Omega}\in{\mathbb R}^{d\times k}$ has iid ${\mathcal{N}}(0,\nicefrac{{1}}{{k}})$ entries.

The OSI properties of a Gaussian test matrix reflect the excellent quality of output. According to [73, Sec. 8.7.2], Gaussian test matrices attain constant injectivity $\alpha\geq\nicefrac{{1}}{{2}}$ for subspaces of dimension $r$ at an embedding dimension $k={\mathcal{O}}(r)$, which motivates the first desideratum (Item˜1). On the other hand, for a vector $\bm{x}\in{\mathbb R}^{d}$, forming the Gaussian sketch $\bm{\Omega}^{\smash{*}}\bm{x}\in{\mathbb R}^{k}$ is expensive, requiring ${\mathcal{O}}(dk)$ arithmetic operations (cf. Item˜2). We aim to design test matrices that match the OSI properties and output quality associated with a Gaussian test matrix, while permitting more economical computations.

1.5.1 Sparse test matrices

For our first in-depth demonstration of the OSI framework, we study sparse test matrices. Sparse test matrices offer compelling theoretical benefits for both sparse and dense linear algebra, and they exhibit spectacular empirical performance [119, 31, 46, 74, 84, 85, 41]. Yet there remain core questions about their design and analysis. Our results justify the application of sparse test matrices at levels of sparsity that were previously unsupportable, and they hint at sparsity levels that had been considered impossible.

The basic case for sparse sketching is simple. If the test matrix $\bm{\Omega}\in\mathbb{F}^{d\times k}$ has only $\zeta$ nonzero entries per row, then we can form the sketch $\bm{A}\bm{\Omega}$ with just ${\mathcal{O}}(\zeta\cdot\operatorname{nnz}(\bm{A}))$ arithmetic operations. For contrast, sketching with an unstructured dense test matrix typically costs ${\mathcal{O}}(k\cdot\operatorname{nnz}(\bm{A}))$ operations. If the sparsity is much smaller than the embedding dimension ($\zeta\ll k$), then sparse test matrices can offer significant acceleration (when implemented with a high-performance sparse arithmetic library). Figure˜3 provides empirical evidence of this speedup.

Among several types of sparse test matrices, we recommend using the SparseStack, an existing construction due to Kane & Nelson [63, Fig. 1(c)], because of its excellent empirical performance. To analyze this sparse test matrix, we want to establish an OSI property at the minimal row sparsity $\zeta$ (to control sketching costs) and with minimal embedding dimension $k$ (to control postprocessing costs).

Definition 1.7** (SparseStack).**

Fix the ambient dimension $d\in{\mathbb N}$, the row sparsity $\zeta\in{\mathbb N}$, and the block size $b\in{\mathbb N}$. The embedding dimension $k\coloneqq b\>\zeta$. A SparseStack is a random matrix $\bm{\Omega}\in\mathbb{F}^{d\times k}$ of the form

[TABLE]

In this display, $\mathbf{e}_{i}\in\mathbb{F}^{b}$ is the $i$th standard basis vector; iid means independent and identically distributed.

The SparseStack test matrix consists of $\zeta$ independent copies of a (scaled) CountSketch test matrix [35], stacked side by side, like books on a shelf. Section˜3.1 discusses implementation matters, while Section˜3.2 outlines reasons to prefer the SparseStack over other sparse matrix constructions.

In a recent paper, Tropp studied the injectivity properties of random sparse matrices [115, Sec. 6]. Without proof, he stated a result that translates into an OSI guarantee for the SparseStack matrix. For reference, we present the full analysis in Appendix˜B.

Imported Theorem 1.8** (SparseStacks are OSIs, [115, Rem. 6.5]).**

A SparseStack test matrix $\bm{\Omega}\in\mathbb{F}^{d\times k}$ serves as an $(r,\nicefrac{{1}}{{2}})$-OSI for some embedding dimension $k={\mathcal{O}}(r)$ and row sparsity $\zeta={\mathcal{O}}(\log r)$. With these parameters, building the sketch $\bm{A}\bm{\Omega}$ has arithmetic cost ${\mathcal{O}}(\operatorname{nnz}(\bm{A})\cdot\log r)$.

Imported Theorem˜1.8 compares favorably with results for the OSE property. Cohen [36] proved that a SparseStack is an $(r,\nicefrac{{1}}{{2}},\nicefrac{{3}}{{2}})$-OSE for some $k={\mathcal{O}}(r\log r)$ and $\zeta={\mathcal{O}}(\log r)$. Chenakkod et al. [32] proved a similar OSE bound for some $k={\mathcal{O}}(r)$ and $\zeta={\mathcal{O}}(\log^{3}r)$. By shifting the goal from OSE to OSI, Imported Theorem˜1.8 achieves the better embedding dimension and row sparsity.

We can now establish end-to-end guarantees for sketching algorithms based on sparse test matrices by combining Tropp’s result (Imported Theorem˜1.8) with our new analysis (Theorems˜1.3 and 1.4). We learn that randomized linear algebra algorithms implemented with the SparseStack test matrix can achieve near-optimal error bounds with fast runtimes. For sparse matrices, our approach yields sketching times within a single logarithmic factor of the best possible. See Table˜2 for a summary of our theoretical results and a comparison with existing research.

Empirical results.

This paper also contains new evidence that SparseStack matrices have exemplary empirical performance. Figures˜2 and 3 illustrate that a SparseStack test matrix with constant row sparsity $\zeta=4$ reliably achieves an OSI guarantee, while the sketching time is up to 88$\times$ faster than the Gaussian benchmark. Section˜3.3 documents additional experiments on difficult synthetic test problems, where SparseStack matrices achieve near-Gaussian output quality on low-rank approximation tasks.

These benefits persist in practical computations. Section˜6.1 showcases an application that employs the generalized Nyström method to compress and analyze scientific simulation data. For this problem SparseStack matrices offer up to a 12$\times$ speedup over the Gaussian baseline, while attaining the same accuracy. Together, these results demonstrate that SparseStack test matrices are highly effective for linear algebra computations, while they are significantly faster than alternatives.

Constant sparsity?

While Imported Theorem˜1.8 allows for row sparsity $\zeta={\mathcal{O}}(\log r)$, practitioners [48, 74, 41, 31, 119] have long observed that sparse test matrices are reliable even with constant row sparsity. For the SparseStack model, our numerical work (Figures˜5 and 2) suggests that the parameters $k=2r$ and $\zeta=4$ already suffice to achieve the OSI property. Thus, we frame a conjecture:

Conjecture 1.9** (SparseStack: OSI with constant row sparsity).**

A SparseStack test matrix is an $(r,\nicefrac{{1}}{{2}})$-OSI for some embedding dimension $k={\mathcal{O}}(r)$ and some sparsity level $\zeta={\mathcal{O}}(1)$.

In contrast, a lower bound of Nelson & Nguyễn [91, Thm. 7] implies that a sparse test matrix with proportional embedding dimension $k={\mathcal{O}}(r)$ cannot serve as an $(r,\alpha,\beta)$-OSE with constant injection and dilation $\alpha,\beta=\Theta(1)$ unless the row sparsity $\zeta=\Omega(\log r)$. Some papers interpret this lower bound as a total prohibition on sparse sketching methods with constant row sparsity. Yet the evidence in Figure˜2 is consistent with both the impossibility result for OSEs and with ˜1.9 on OSIs. Since the OSI property already suffices to implement many sketching algorithms, a proof of the conjecture would justify the deployment of sparse test matrices with constant row sparsity, leading to optimal runtimes for several linear algebra tasks.

1.5.2 Test matrices based on fast trigonometric transforms

Next, we study test matrices built from trigonometric transforms with fast algorithms, such as the discrete Fourier transform (DFT), the discrete cosine transform (DCT), or the Walsh–Hadamard transform (WHT). Historically, these constructions were among the first examples of structured test matrices [124, 4, 111, 58]. Sketching methods based on trigonometric transforms remain popular because they show strong empirical performance [118, 87, 88]. They also benefit from the wide availability of software [52] and hardware [108] for fast transforms. Within this ecosystem, we will design an effective test matrix that achieves superior theoretical guarantees by exploiting the OSI framework.

The standard template for a transform-based test matrix consists of three matrix factors:

[TABLE]

The diagonal matrix $\bm{D}$ is usually populated with iid Rademacher variables. Typical choices for the trigonometric transform $\bm{F}$ include the WHT or the DCT (when $\mathbb{F}\in\{{\mathbb R},{\mathbb C}\}$) or the DFT (when $\mathbb{F}={\mathbb C}$). The sparse matrix $\bm{S}$ enacts a fast sampling operation that reduces the dimension and rescales the output so that $\bm{\Omega}$ is isotropic. For an input matrix $\bm{A}\in\mathbb{F}^{n\times d}$, computing the sketch $\bm{A}\bm{\Omega}$ usually costs ${\mathcal{O}}(nd\log d)$ operations for the trigonometric transform, plus the cost of applying the sparse sampling matrix.

To complete the construction ˜1.4, we must design a sparse sampling matrix $\bm{S}\in\mathbb{F}^{d\times k}$ that delivers fast sketching times along with injectivity guarantees. The most common choice for $\bm{S}$ is a restriction matrix that selects $k$ coordinates from $\{1,\ldots,d\}$ uniformly at random; the resulting test matrix is called a subsampled randomized trigonometric transform (SRTT)  [73, Sec. 9.3]. To achieve an OSI with subspace dimension $r$, standard SRTTs require embedding dimension $k=\Omega(r\log r)$ to guard against the worst-case input matrix [111, Sec. 3.3]. This fact has prompted researchers to explore alternative forms of the sparse matrix [26, 33, 32] that allow for the minimal embedding dimension $k={\mathcal{O}}(r)$.

We recommend using the SparseRTT, a transform-based test matrix that achieves the optimal embedding dimension and the minimal sketching cost within the class ˜1.4.

Definition 1.10** (SparseRTT).**

Fix a column sparsity parameter $\xi\in{\mathbb N}$. A SparseRTT is a random matrix $\bm{\Omega}\coloneqq\bm{D}\bm{F}\bm{S}\in\mathbb{F}^{d\times k}$, where the diagonal matrix $\bm{D}\in\mathbb{F}^{d\times d}$ has iid Rademacher entries, and the trigonometric transform $\bm{F}\in\mathbb{F}^{d\times d}$ is a unitary WHT, DCT, or DFT. The sparse matrix $\bm{S}\in\mathbb{F}^{d\times k}$ is a SparseCol random matrix (Definition˜4.1) with exactly $\xi$ nonzero entries per column.

Theorem 1.11** (SparseRTTs are OSIs, variant of Theorem˜4.2).**

Let the subspace dimension $r\geq\log d$. A SparseRTT test matrix $\bm{\Omega}\in\mathbb{F}^{d\times k}$ serves as an $(r,\nicefrac{{1}}{{2}})$-OSI for some embedding dimension $k={\mathcal{O}}(r)$ and some column sparsity $\xi={\mathcal{O}}(\log r)$. For a matrix $\bm{A}\in\mathbb{F}^{n\times d}$, the sketch $\bm{A}\bm{\Omega}$ has arithmetic cost ${\mathcal{O}}(nd\log r)$.

In the setting of Theorem˜1.11, the random sparse matrix $\bm{S}$ has $\Theta(r\log r)$ entries. This is the minimum sparsity level possible because of the SRTT counterexample [111, Sec. 3.3]. In particular, the sketching cost cannot be reduced further with an (oblivious) test matrix of the form ˜1.4. For a comparison with prior theoretical work [4, 33, 32], Table˜4 lists the costs of the sketch-and-solve method (Algorithm˜3) with various trigonometric test matrices of the form ˜1.4. Our numerical work (Figure˜6) confirms that the SparseRTT outperforms other trigonometric transform test matrices, with accuracy matching a Gaussian test matrix. We discuss additional implementation details in Section˜4.3.

For a dense unstructured input matrix $\bm{A}$, the SparseRTT and the SparseStack (Definition˜1.7) yield OSIs with the same arithmetic costs, up to constant factors. Sparse test matrices are typically faster than those based on randomized trigonometric transforms because they exploit hardware-level parallelism more effectively and involve less data movement. Therefore, for most applications, we recommend the SparseStack over the SparseRTT, but SparseRTT test matrices remain useful in environments that lack sparse linear algebra libraries.

1.5.3 Test matrices with tensor structure

Many linear algebra problems in machine learning and scientific computing exhibit tensor structure. Because of the towering scale of these problems, randomized dimension reduction serves as a powerful ally. Nevertheless, sketching methods for tensors remain poorly understood. This section employs the OSI framework to develop and analyze sketching methods for tensor data.

Khatri–Rao test matrices.

The Kronecker matvec access model [79, 81] offers a useful abstraction for operations on tensor data. Here is the simplest version of the model. Fix the base dimension $d_{0}\in{\mathbb N}$ and the tensor order $\ell\in{\mathbb N}$. We interact with an input matrix $\bm{A}\in\mathbb{F}^{n\times d_{0}^{\ell}}$ by forming the matvec $\bm{A}\boldsymbol{\omega}\in\mathbb{F}^{n}$ with a vector $\boldsymbol{\omega}\in\mathbb{F}^{d_{0}^{\ell}}$ that has Kronecker product structure:

[TABLE]

Section˜5.2 lists several applications that are well described by the Kronecker matvec access model.

Let us outline a natural construction of a random test matrix that we can implement in the Kronecker matvec access model. First, choose an isotropic random vector $\boldsymbol{\nu}\in\mathbb{F}^{d_{0}}$, called the base distribution. The columns of the test matrix are tensor products of independent copies of the base distribution. Common examples of base distributions include the real and complex Rademacher distributions, as well as the real and complex Gaussian distributions. We also treat the real and complex spherical distributions; that is, \boldsymbol{\nu}\sim\textsc{uniform}\big{\{}\bm{x}\in\mathbb{F}^{d_{0}}\mathrel{\mathop{\ordinarycolon}}\mathinner{\!\left\lVert\bm{x}\right\rVert}_{2}=\sqrt{d_{0}}\big{\}}. For additional details, see Sections˜1.8 and 5.1.

Definition 1.12** (Khatri–Rao test matrix).**

Fix base dimension $d_{0}\in{\mathbb N}$, tensor order $\ell\in{\mathbb N}$, and embedding dimension $k\in{\mathbb N}$. Select an isotropic base distribution $\boldsymbol{\nu}$ on $\mathbb{F}^{d_{0}}$. Define the Khatri–Rao test matrix

[TABLE]

For an input matrix $\bm{A}\in\mathbb{F}^{n\times d_{0}^{\ell}}$, we can form the sketch $\bm{A}\bm{\Omega}\in\mathbb{F}^{n\times k}$ by performing $k$ queries in the Kronecker matvec access model, one for each column of the test matrix. Moreover, the test matrix $\bm{\Omega}$ inherits the isotropy property ˜1.2 from the base distribution $\boldsymbol{\nu}$. The matrix $\bm{\Omega}$ can also be written as a Khatri–Rao product [66, p. 462] of independent random matrices, but we avoid this formalism.

OSI properties.

There is an extensive literature on randomized linear algebra computations with Khatri–Rao test matrices; e.g., see [12, 3, 99, 28, 24]. To advance our understanding, we establish that Khatri–Rao test matrices are OSIs. Here is an essential special case, extracted from Theorems˜5.1 and 5.2.

Theorem 1.13** (Gaussian Khatri–Rao test matrices are OSIs).**

Assume that the base distribution $\boldsymbol{\nu}\in{\mathbb R}^{d_{0}}$ is the real standard normal distribution, and fix the tensor order $\ell$. For an embedding dimension $k$, construct the Khatri–Rao test matrix $\bm{\Omega}\in{\mathbb R}^{d_{0}^{\ell}\times k}$ described by Definition˜1.12.

1. (a)

Exponential embedding dimension, constant injectivity.* The test matrix $\bm{\Omega}$ serves as an $(r,\nicefrac{{1}}{{2}})$-OSI for some embedding dimension $k={\mathcal{O}}(3^{\ell}r)$. .* 2. (b)

Linear embedding dimension, exponentially-small injectivity.* The test matrix $\bm{\Omega}$ serves as an $(r,\mathrm{c}^{\ell})$-OSI for some embedding dimension $k={\mathcal{O}}(r)$, where $\mathrm{c}\in(0,1)$ is a universal constant.*

Item˜(a)* extends to many other distributions (Theorem˜5.1), while Item˜(b) is more delicate (Theorem˜5.2).*

Theorem˜1.13 suggests that the quality of a Khatri–Rao test matrix degrades quickly as the tensor order $\ell$ increases. Item˜(a) states that we can achieve constant injectivity $\alpha$ when the embedding dimension $k$ is exponential in the tensor order $\ell$. Item˜(b) states that we can reach the optimal embedding dimension $k={\mathcal{O}}(r)$ at the cost of exponentially small injectivity $\alpha$. The exponential scaling is not an artifact of the analysis—it is necessary for worst-case examples [81]. Khatri–Rao sketching also has the surprising feature [79, 81] that changing the base distribution can result in an exponentially large difference in the embedding dimension $k$ required for constant injectivity $\alpha\geq\mathrm{c}$. See Section˜5 for details.

We are not aware of any result similar to Theorem˜1.13, aside from the concurrent research [104], although the literature contains OSE guarantees in special cases. For example, Bujanović et al. [24] prove that a Khatri–Rao test matrix with a real Gaussian base distribution and tensor order $\ell=2$ is an OSE when the embedding dimension $k={\mathcal{O}}(r^{3/2})$.

Practitioner recommendations.

In the literature, the most common base distributions for Khatri–Rao test matrices are the real Gaussian and real Rademacher distributions. We advise against using these base distributions. Instead, we recommend the complex spherical distribution. If complex arithmetic is not available, then we recommend using the real spherical distribution.

Why? To guarantee the OSI property with constant injectivity $\alpha\geq\nicefrac{{1}}{{2}}$, the real Gaussian distribution requires embedding dimension $k={\mathcal{O}}(\smash{3^{\ell}}r)$. By using a spherical distribution, the base of the exponent shrinks by a factor depending on the base dimension $d_{0}$. In particular, in the important case $d_{0}=2$, the real spherical distribution allows for $k={\mathcal{O}}(\smash{1.5^{\ell}}r)$ and the complex spherical distribution permits $k={\mathcal{O}}(\smash{1.\overline{33}}^{\ell}r)$. The argument against using a Rademacher base distribution is even more brutal, as Khatri–Rao test matrices over this distribution can fail catastrophically. Indeed, for a worst-case subspace, the Rademacher Khatri–Rao test matrix has injectivity $\alpha=0$ unless the embedding dimension $k=\Omega(2^{\ell})$. Section˜5.3.1 presents numerical experiments that document these effects.

Computational consequences.

While the exponential factors in Theorem˜1.13 may seem deadly, Khatri–Rao test matrices remain powerful tools for scientific computing.

First, many applications of tensor linear algebra concern tensors with small order, such as $\ell=2$ or $\ell=3$. In these cases, the exponential factors in $\ell$ become constants, and we recover near-optimal approximations for some embedding dimension $k={\mathcal{O}}(r)$. As a concrete example, the next theorem treats the problem of structured matrix approximation under a bilinear query model.

Theorem 1.14** (Matrix recovery from bilinear queries, variant of Theorem˜5.3.).**

Choose a $d$-dimensional subspace of matrices ${\mathcal{F}}\coloneqq\operatorname{span}\{\bm{M}_{1},\ldots,\bm{M}_{d}\}\subseteq{\mathbb R}^{n\times n}$. Fix a target matrix $\bm{B}\in{\mathbb R}^{n\times n}$, accessible with bilinear queries: $(\bm{x},\bm{y})\mapsto\bm{y}^{\sf T}\bm{B}\bm{x}$. Algorithm˜6 performs ${\mathcal{O}}(d)$ bilinear queries on $\bm{B}$, and it returns a near-optimal orthogonal projection $\smash{{\tilde{\smash{\boldsymbol{\bm{B}}}\mathstrut\mkern 4.0mu}\mkern-4.0mu}{}}\in{\mathbb R}^{n\times n}$ of the target matrix $\bm{B}$ onto the subspace ${\mathcal{F}}$ with probability at least $\>\nicefrac{{9}}{{10}}$. That is,

[TABLE]

As usual, $\mathrm{C}\geq 1$ is a universal constant.

The subspace ${\mathcal{F}}$ of structured matrices is sometimes called a linearly parameterized family. Typical examples include the set of Toeplitz matrices or the set of matrices with fixed bandwidth. The bilinear query model arises in many settings [110, 29, 10, 70, 110], and it has been the focus of several theoretical studies [122, 101, 89]. To prove Theorem˜1.14, we observe that a family of random bilinear queries amounts to a Khatri–Rao sketch of (the vectorization of) the target matrix $\bm{B}$. Theorem˜1.13 provides an OSI for this test matrix, which implies the optimal query complexity ${\mathcal{O}}(d)$. In contrast, existing OSE results for Khatri–Rao test matrices [24] yield a weaker query complexity of ${\mathcal{O}}(d^{3/2})$. The concurrent research [104] offers an improvement, but it also falls short of the optimal bound.

Second, although Khatri–Rao test matrices behave exponentially badly for worst-case problem instances, their empirical performance is often substantially better than the worst-case theory suggests. As evidence, Figure˜1 applies (complex spherical) Khatri–Rao test matrices with large tensor order $\ell$ to accurately approximate thousands of matrices from SparseSuite. In Section˜6.2, we successfully perform variance-reduced trace estimation [80, 95, 49] to compute the partition function of a quantum system. These experiments require Khatri–Rao products with tensor order $\ell=16$, and yet the methods are powerful enough to achieve 12 digits of precision in regimes where baseline methods fail to achieve even 1 digit of precision.

1.6 How to establish the OSI property

So far, we have argued that many random test matrices exhibit the OSI property. Moreover, when we use these test matrices to implement randomized linear algebra algorithms, we can achieve near-optimal error bounds and state-of-the-art runtimes. Still, some conceptual questions nag us. How is the OSI property fundamentally different from the OSE property? Why is the OSI property more likely to hold? How can we prove that a test matrix is an OSI?

Throughout this section, we assume that the test matrix has iid columns. More precisely, let $\boldsymbol{\omega}\in\mathbb{F}^{d}$ be an isotropic random vector: $\operatorname*{\mathbb{E}}{}[\boldsymbol{\omega}\boldsymbol{\omega}^{\smash{*}}]=\mathbf{I}$. Consider the random test matrix

[TABLE]

Intuition: Dilation versus injection.

On closer inspection of the embedding property ˜1.1, we discover a fundamental asymmetry between the dilation and injection conditions. To see why, write

[TABLE]

Note that the right-hand side is the average of $k$ iid nonnegative summands, each with expectation $\mathinner{\!\left\lVert\bm{x}\right\rVert}_{2}^{2}$.

We realize that a sufficient condition for the test matrix $\bm{\Omega}$ to stretch out the vector $\bm{x}$ by a large factor is for one summand in ˜1.6 to be unusually large. When linear marginals of $\boldsymbol{\omega}$ have heavy tails, the dilation parameter $\beta$ tends to be large. On the other hand, a necessary condition for the test matrix $\bm{\Omega}$ to annihilate the vector $\bm{x}$ is for all summands in ˜1.6 to equal zero. Therefore, the injectivity parameter $\alpha$ hates to be zero. In summary: It is almost a law of nature that a test matrix is an OSI (Definition˜1.2) when the embedding dimension $k$ exceeds the subspace dimension $r$ by a constant factor; cf. [67, p. 12994].

Spectral formulations.

We still need to identify strategies for establishing subspace injectivity. To that end, observe that the OSI property admits an equivalent formulation as a spectral inequality.

Proposition 1.15** (OSI: Spectral formulation).**

Suppose that $\bm{\Omega}\in\mathbb{F}^{d\times k}$ is an isotropic random matrix, as in ˜1.2. The following conditions are equivalent.

The test matrix $\bm{\Omega}$ is an $(r,\alpha)$-OSI. 2. 2.

For each fixed matrix $\bm{Q}\in\mathbb{F}^{d\times r}$ with orthonormal columns,

[TABLE]

We define the minimum singular value and the minimum eigenvalue as

[TABLE]

In particular, the minimum singular value of a matrix is zero if and only the matrix has a trivial null space.

Therefore, we can check the OSI property by obtaining lower bounds for the minimum singular value of certain random matrices. Fortuitously, the literature on high-dimensional probability contains several powerful tools for studying the minimum singular value of a random matrix.

Fourth-moment bounds (Section˜7.1).

Oliveira [92, Thm. 1.1] established that fourth-moment bounds suffice to control the minimum eigenvalue of a sample covariance matrix. In detail, let $\boldsymbol{\omega}\in\mathbb{F}^{d}$ be an isotropic random vector where the fourth moments of the linear marginals are uniformly bounded. That is, for a fixed number $h\geq 1$, we assume that

[TABLE]

Consider a random test matrix $\bm{\Omega}\in\mathbb{F}^{d\times k}$ with iid columns, as in ˜1.5. Then the random matrix $\bm{\Omega}$ is an $(r,\nicefrac{{1}}{{2}})$-OSI for some embedding dimension $k={\mathcal{O}}(hr)$. We use this strategy to establish Item˜(a) of Theorem˜1.13 on Khatri–Rao test matrices with constant injectivity.

Small-ball probabilities (Section˜7.2).

Koltchinskii & Mendelson [67, Thm. 1.3] have demonstrated that small-ball probability bounds allow us to control the minimum singular value of a random matrix. We employ a variant of their approach, adapted from [56, 114, 120]. Consider an isotropic random vector $\boldsymbol{\omega}\in\mathbb{F}^{d}$ whose linear marginals admit a uniform bound on the small-ball probability. That is, for a fixed number $\tau\in(0,1)$, we assume that

[TABLE]

Consider a random test matrix $\bm{\Omega}\in\mathbb{F}^{d\times k}$ with iid columns, as in ˜1.5. Then the test matrix is an $(r,\nicefrac{{\tau}}{{10}})$-OSI for some embedding dimension $k={\mathcal{O}}(r)$. We apply this strategy to prove Item˜(b) of Theorem˜1.13 on Khatri–Rao test matrices with proportional embedding dimension.

Gaussian comparison (Section˜7.3).

Tropp [115, Thm. 2.3] recently devised a method for bounding the minimum eigenvalue of a random psd matrix by comparison with the minimum eigenvalue of a suitable Gaussian random matrix. This approach is the most elaborate of the three techniques we consider. Here is the simplest setting. Consider a random test matrix $\bm{\Omega}\in\mathbb{F}^{d\times k}$ with iid columns, as in ˜1.5. For a fixed matrix $\bm{Q}\in\mathbb{F}^{d\times r}$ with orthonormal columns, we introduce the random psd matrix

[TABLE]

The idea is to construct a Gaussian matrix $\bm{Z}\in\mathbb{F}^{r\times r}$ whose statistics are determined by the summands:

[TABLE]

In an appropriate sense, the minimum eigenvalue $\lambda_{\min}(\bm{Y})$ is stochastically bounded below by the minimum eigenvalue $\lambda_{\min}(\bm{Z})$. Since $\bm{Z}$ is Gaussian, we can employ specialized methods to study its minimum eigenvalue. This strategy is effective for analyzing sparse test matrices. We use it to establish Theorem˜1.11 about SparseRTTs and to confirm Imported Theorem˜1.8 about SparseStacks.

1.7 Outline

This paper carries out the OSI program for studying randomized linear algebra algorithms with structured test matrices. Section˜2 shows that the OSI property suffices to justify several fundamental algorithms. Afterward, we establish OSI properties and discuss constructions of sparse test matrices (Section˜3), fast trigonometric transforms (Section˜4), and Khatri–Rao products (Section˜5). Section˜6 highlights two scientific applications of randomized linear algebra that benefit from structured test matrices. Last, Section˜7 elaborates on the proof strategies that we use to check the OSI property. ¡Vámonos!

1.8 Notation & terminology

Our results are valid in either the real field ${\mathbb R}$ or the complex field ${\mathbb C}$; the symbol $\mathbb{F}$ refers to either field. The bar $\overline{z}$ denotes the complex conjugate of $z\in{\mathbb C}$. For vectors and matrices, the asterisk ∗ denotes the conjugate transpose, while the superscript T denotes the plain transpose. The vector space ${\mathbb F}^{d}$ is equipped with the canonical inner product $\langle\bm{x},\bm{y}\rangle\coloneqq\bm{x}^{\smash{*}}\bm{y}$, which induces the $\ell_{2}$ norm $\mathinner{\!\left\lVert\cdot\right\rVert}_{2}$. We define the standard Kronecker product $\otimes$ and the vectorization operation $\operatorname{vec}(\cdot)$.

A positive semidefinite (psd) matrix is self-adjoint and has nonnegative eigenvalues. An orthonormal matrix $\bm{Q}\in\mathbb{F}^{n\times r}$ satisfies the condition $\bm{Q}^{\smash{*}}\bm{Q}=\mathbf{I}$. The coherence of an orthonormal matrix is the numerical quantity $\mu(\bm{Q})\coloneqq\max_{1\leq i\leq n}\ \mathinner{\!\left\lVert\smash{\mathbf{e}_{i}^{*}\bm{Q}}\right\rVert}_{2}^{2}$, where $\mathbf{e}_{i}$ denotes the $i$th canonical basis vector. We work with the spectral norm $\mathinner{\!\left\lVert\cdot\right\rVert}_{2}$, the Frobenius norm $\mathinner{\!\left\lVert\cdot\right\rVert}_{\rm F}$, and the nuclear norm $\mathinner{\!\left\lVert\cdot\right\rVert}_{*}$. For a matrix $\bm{A}$, the symbol $\llbracket\bm{A}\rrbracket_{r}$ denotes a best rank-$r$ approximation with respect to the Frobenius norm. The adjective “rank-$r$” means “with rank at most $r$.”

We employ two sketching matrices, which are always called $\bm{\Omega}$ and $\bm{\Psi}$, each of which is tall and skinny. Consequently, $\bm{\Omega}$ always acts on the right ($\bm{A}\mapsto\bm{A}\bm{\Omega}$), while $\bm{\Psi}^{\smash{*}}$ acts on the left ($\bm{A}\mapsto\bm{\Psi}^{\smash{*}}\bm{A}$). The symbol $\zeta$ represents the (expected) number of nonzeros per row of a sparse random matrix, while $\xi$ is the (expected) number of nonzeros per column.

The operator $\mathbb{P}$ returns the probability of an event, while $\operatorname*{\mathbb{E}}$ and $\operatorname*{Var}$ compute the expectation and variance of a random variable. The symbol $\sim$ means “has the distribution.” As always, iid abbreviates “independent and identically distributed.” The real standard normal distribution ${\mathcal{N}}(0,1)$ describes a real Gaussian random variable with expectation zero and variance one. The complex standard normal distribution ${\mathcal{N}}_{{\mathbb C}}(0,1)$ induces a complex random variable of the form $(Z_{1}+\mathrm{i}Z_{2})/\sqrt{2}$, where $Z_{1},Z_{2}\sim{\mathcal{N}}(0,1)$ iid. The real Rademacher distribution $\textsc{rademacher}=\textsc{uniform}\{\pm 1\}$. The complex Rademacher distribution $\textsc{rademacher}_{{\mathbb C}}$ induces a random variable of the form $(\varrho_{1}+\mathrm{i}\varrho_{2})/\sqrt{2}$ where $\varrho_{1},\varrho_{2}\sim\textsc{rademacher}$ iid. We also define another complex distribution, $\textsc{steinhaus}=\textsc{uniform}\{z\in{\mathbb C}\mathrel{\mathop{\ordinarycolon}}|z|=1\}$. The spherical distribution on $\mathbb{F}^{d}$ is uniform over the set of vectors with $\ell_{2}$ norm $\sqrt{d_{0}}$.

We employ the standard computer-science definitions of the order symbols ${\mathcal{O}},\Omega,\Theta$. The notation $\widetilde{}{\mathcal{O}}(\cdot)$ suppresses polylogarithmic factors. As usual, the upright letters $\mathrm{c},\mathrm{C}$ denote positive universal constants that may change from appearance to appearance.

2 Randomized matrix computations via OSIs

In this section, we review standard algorithms for least-squares regression and low-rank approximation that employ randomized dimension reduction. When these algorithms are implemented with OSIs, we can prove that they achieve errors within a constant factor of the optimal value. The runtime of these algorithms also depends on the choice of the OSI, so it is productive to isolate the sketching costs from the auxiliary arithmetic performed by each algorithm; see Table˜1.

Section˜2.1 contains our main technical observation, which states that OSIs “respect orthogonality.” In Section˜2.2, we discuss the randomized SVD algorithm, which is a two-pass algorithm for computing a low-rank approximation. Sections˜2.3 and 2.5 discuss the Nyström and generalized Nyström approximations, which produce low-rank approximations after a single pass over the matrix. Section˜2.4 treats the sketch-and-solve algorithm for least-squares, a fundamental method that also figures in the analysis of the generalized Nyström approximation.

2.1 Basic technical observation: OSIs respect orthogonality

Before treating specific algorithms, we highlight a core theoretical insight:

Lemma 2.1** (OSIs respect orthogonality).**

Fix two orthonormal matrices $\bm{Q}\in\mathbb{F}^{n\times r}$ and $\bm{Q}_{\perp}\in\mathbb{F}^{n\times s}$ whose ranges are mutually orthogonal: $\bm{Q}^{\smash{*}}\bm{Q}_{\perp}=\bm{0}$. Choose an arbitrary matrix $\bm{B}\in{\mathbb F}^{t\times s}$. Draw a random $(r,\alpha)$-OSI test matrix $\bm{\Omega}\in{\mathbb F}^{n\times k}$. With probability at least $\frac{9}{10}$, the matrix $\bm{Q}^{\smash{*}}\bm{\Omega}$ has full row rank, and

[TABLE]

We can interpret the mysterious quantity $F\coloneqq\mathinner{\!\left\lVert\smash{\bm{B}(\bm{Q}_{\perp}^{\smash{*}}\bm{\Omega})(\bm{Q}^{\smash{*}}\bm{\Omega})^{\smash{\dagger}}}\right\rVert}_{\rm F}$ as the minimum Frobenius norm of a solution to the sketched least-squares problem

[TABLE]

In words, $F$ describes how well we can approximate the point $\bm{Q}_{\perp}\bm{B}^{\smash{*}}$ by a point in the orthogonal subspace $\operatorname{range}(\bm{Q})$ after sketching both with $\bm{\Omega}$. The lemma states that an OSI $\bm{\Omega}$ roughly preserves the orthogonality between the point and the subspace. As we will discover, the result is valuable because $F$ provides an exact expression for the residual in the sketch-and-solve method for least-squares regression. It also controls the error in randomized algorithms for low-rank approximation of matrices.

Proof.

Using the operator-norm bound for the Frobenius norm, we have

[TABLE]

Since $\bm{Q}$ has $r$ columns and $\bm{\Omega}$ is an $(r,\alpha)$-OSI, Proposition˜1.15 ensures that $\sigma_{\rm min}^{2}(\bm{\Omega}^{\smash{*}}\bm{Q})\geq\alpha$ with probability at least $1-\frac{1}{20}$. After we expand the numerator, the isotropy property ˜1.2 of the OSI yields

[TABLE]

Indeed, $\operatorname*{\mathbb{E}}{}[\bm{\Omega}\bm{\Omega}^{{\smash{*}}}]=\mathbf{I}$. By Markov’s inequality,

[TABLE]

With a union bound, we can control the numerator and denominator in ˜2.1 simultaneously with probability at least $1-\frac{1}{10}$. On this event, we can combine the inequalities to reach the stated result. ∎

The rest of this section establishes that we can use OSIs to implement RSVD, Nyström, generalized Nyström, and sketch-and-solve regression. The proofs all hinge on Lemma˜2.1.

2.2 The randomized SVD

The randomized SVD [58], or RSVD, is a basic and widely used randomized algorithm for computing a low-rank approximation of a matrix.

Fix an input matrix $\bm{A}\in\mathbb{F}^{n\times d}$ and an approximation rank $k\leq\min\{n,d\}$. We draw a random test matrix $\bm{\Omega}\in\mathbb{F}^{d\times k}$, and we acquire information about the input matrix via the sketch $\bm{Y}\coloneqq\bm{A}\bm{\Omega}\in\mathbb{F}^{n\times k}$. The range of the sketch $\bm{Y}$ captures important directions in the column space of $\bm{A}$. The orthogonal projection of the input matrix $\bm{A}$ onto the range of the sketch $\bm{Y}$ induces a rank-$k$ approximation of the form $\widehat{\bm{A}}\coloneqq\bm{Y}\bm{Y}^{{\smash{\dagger}}}\bm{A}\in\mathbb{F}^{n\times d}.$ The RSVD algorithm reports a compact SVD of this approximation: $\widehat{\bm{A}}=\bm{U}\bm{\Sigma}\bm{V}^{\smash{*}}$.

For numerical stability, the RSVD algorithm accomplishes this task by forming an orthonormal basis for the range of the test matrix:

[TABLE]

The rank-$k$ approximation $\widehat{\bm{A}}$ of the input matrix can be expressed with the QB factorization:

[TABLE]

To obtain the approximation, we need to compute the matrix–matrix product $\bm{Q}^{*}\bm{A}=(\bm{A}^{\smash{*}}\bm{Q})^{\smash{*}}\in\mathbb{F}^{k^{\prime}\times d}$, which is often the dominant cost of the RSVD algorithm. Afterward, the RSVD algorithm manipulates the QB factorization to obtain a compact SVD.

See Algorithm˜1 for reliable RSVD pseudocode. The RSVD algorithm performs two passes over the input matrix $\bm{A}$, including one sketching operation ($\bm{A}\bm{\Omega}$) and one matrix–matrix product ($\bm{Q}^{*}\bm{A}$). It also requires ${\mathcal{O}}(k^{2}n)$ operations for orthogonalization and SVD computations. While we can exploit properties of the input matrix $\bm{A}$ in the matrix–matrix product $\bm{Q}^{*}\bm{A}$, this step does not benefit from the structure of the test matrix $\bm{\Omega}$. See Table˜1 for a summary of the arithmetic costs. Sections˜2.3 and 2.5 describe low-rank approximation algorithms that can achieve arithmetic costs strictly smaller than the RSVD.

Here is our main theoretical result for the RSVD. Recall that $\llbracket\bm{A}\rrbracket_{r}$ denotes a best rank-$r$ approximation of a matrix $\bm{A}$ with respect to the Frobenius norm.

Theorem 2.2** (Randomized RSVD with an OSI).**

Fix an input matrix $\bm{A}\in\mathbb{F}^{n\times d}$ and a target rank $r\leq\min\{n,d\}$. Draw a random $(r,\alpha)$-OSI test matrix $\bm{\Omega}\in\mathbb{F}^{d\times k}$ where $k\geq r$. Construct the rank-$k$ RSVD approximation $\widehat{\bm{A}}\in\mathbb{F}^{n\times d}$ determined by ˜2.2 and ˜2.3. Then, with probability at least $\frac{9}{10}$,

[TABLE]

Theorem˜2.2 provides a bicriteria guarantee for the RSVD approximation. The approximation error for the rank-$k$ output of the RSVD lies within a constant factor of the best rank-$r$ approximation error. This type of bound is typical for the RSVD and related algorithms [58, 116, 55, 123]. We can upgrade the rank-$k$ approximation to a rank-$r$ approximation by a standard truncation technique [58, Thm. 9.3].

The proof of Theorem˜2.2 combines the standard RSVD error bound with the core technical observation about OSIs from Section˜2.1. Halko et al. employed a similar approach in their analysis of the RSVD with SRFT test matrices [58, Thm. 11.2], as do Saibaba & Międlar [103]. A related idea appears in the tutorial of Drineas & Mahoney [44, Rem. 98].

Imported Theorem 2.3** (RSVD error decomposition [58, Thm. 9.1], [72, Lem. 2]).**

Fix an input matrix $\bm{A}\in\mathbb{F}^{n\times d}$. Choose an arbitrary test matrix $\bm{\Omega}\in\mathbb{F}^{d\times k}$, let $\bm{Q}\in\mathbb{F}^{n\times k}$ be an orthonormal matrix whose range contains the range of $\bm{A}\bm{\Omega}$, and construct the approximation $\smash{{\widehat{\smash{\boldsymbol{\bm{A}}}\mathstrut\mkern 4.0mu}\mkern-4.0mu}{}}=\bm{Q}\bm{Q}^{\smash{*}}\bm{A}\in\mathbb{F}^{n\times n}$.

To analyze the RSVD, fix $r\leq k$, and introduce the partitioned SVD of the input matrix:

[TABLE]

where $\bm{\Sigma}_{1}\in\mathbb{F}^{r\times r}$ and $\bm{V}_{1}\in\mathbb{F}^{d\times r}$ contain the top $r$ singular values and right singular vectors of $\bm{A}$, while $\bm{\Sigma}_{2}$ and $\bm{V}_{2}$ contain the other singular values and right singular vectors. Provided that $\bm{V}_{1}^{\smash{*}}\bm{\Omega}$ has full row-rank,

[TABLE]

Proof of Theorem˜2.2.

To proceed, we simply invoke Lemma˜2.1 to bound the excess error term from Imported Theorem˜2.3. Thus, with probability at least $\frac{9}{10}$, we can compute

[TABLE]

In the last line, we applied the Eckart–Young theorem to rewrite $\mathinner{\!\left\lVert\bm{\Sigma}_{2}\right\rVert}_{\rm F}^{2}=\mathinner{\!\left\lVert\bm{A}-\llbracket\bm{A}\rrbracket_{r}\right\rVert}_{\rm F}^{2}$. Finally, note that $1+\mathrm{C}/\alpha\leq\mathrm{C}^{\prime}/\alpha$, where $\mathrm{C}^{\prime}$ is a larger constant. ∎

2.3 Nyström approximation

When the input matrix is psd, we can obtain a low-rank approximation in a single pass using the Nyström approximation [117, 55, 43, 58].

Fix a psd input matrix $\bm{A}\in\mathbb{F}^{n\times n}$ and an approximation rank $k\leq n$. Draw a random test matrix $\bm{\Omega}\in\mathbb{F}^{n\times k}$, and form the sketch $\bm{Y}\coloneqq\bm{A}\bm{\Omega}\in\mathbb{F}^{n\times k}$. The Nyström approximation is the rank-$k$ psd matrix

[TABLE]

The Nyström approximation algorithm manipulates the approximation ˜2.4 to obtain its eigenvalue decomposition $\widehat{\bm{A}}=\bm{U}\bm{\Lambda}\bm{U}^{\smash{*}}$. A numerically stable implementation of this procedure demands care. Algorithm˜2 presents reliable pseudocode, adapted from [69, 117]. The algorithm performs a single pass over the input matrix $\bm{A}$. It requires one sketch of the input matrix, a sketch of intermediate data, and ${\mathcal{O}}(nk^{2})$ additional arithmetic operations. See Table˜1.

To analyze the Nyström approximation, we exploit an algebraic relationship with the RSVD approximation [54, 55]. This equivalence, recently termed the Gram correspondence [47], allows us to transfer Frobenius-norm error bounds for the RSVD to nuclear-norm guarantees for the Nyström approximation.

Imported Theorem 2.4** (Gram correspondence [55, Lem. 1]).**

Let $\bm{A}\in\mathbb{F}^{n\times n}$ be a psd matrix. Choose an arbitrary test matrix $\bm{\Omega}\in\mathbb{F}^{n\times k}$, and form the orthonormal matrix $\bm{Q}=\operatorname{orth}(\bm{A}^{1/2}\bm{\Omega})$. Then the Nyström approximation $\widehat{\bm{A}}\in\mathbb{F}^{n\times n}$, defined in ˜2.4, satisfies

[TABLE]

Put simply, the Nyström approximation applied to $\bm{A}$ is equivalent to the RSVD applied to $\bm{A}^{1/2}$. Combining the Gram correspondence (Imported Theorem˜2.4) and the RSVD analysis (Theorem˜2.2), we immediately obtain an error bound for the single-pass Nyström approximation. This argument is entirely standard [54, 55, 117, 73, 116, 47, 96].

Corollary 2.5** (Randomized Nyström approximation with an OSI).**

Fix a psd input matrix $\bm{A}\in\mathbb{F}^{n\times n}$ and a target rank $r\leq n$. Draw a random $(r,\alpha)$-OSI test matrix $\bm{\Omega}\in\mathbb{F}^{n\times k}$ where $k\geq r$. Construct the rank-$k$ psd Nyström approximation $\widehat{\bm{A}}\in\mathbb{F}^{n\times n}$ determined by ˜2.4. Then, with probability at least $\frac{9}{10}$,

[TABLE]

2.4 Sketch-and-solve for least-squares regression

Sketch-and-solve [105] is a standard randomized algorithm for least-squares problems. If we use an OSI to implement the sketch-and-solve method, the residual error lies within a constant factor of optimal.

Suppose that we are given a (tall) design matrix $\bm{A}\in\mathbb{F}^{n\times d}$ with $n\geq d$ and a response matrix $\bm{B}\in\mathbb{F}^{n\times m}$. We wish to solve an overdetermined least-squares problem:

[TABLE]

To quickly obtain an approximate solution $\smash{{\tilde{\smash{\boldsymbol{\bm{X}}}\mathstrut\mkern 4.0mu}\mkern-4.0mu}{}}\in\mathbb{F}^{d\times m}$, we can employ the sketch-and-solve method. Draw a random test matrix $\bm{\Psi}\in\mathbb{F}^{n\times p}$, and extract the two sketches $\bm{\Psi}^{\smash{*}}\bm{A}\in\mathbb{F}^{p\times d}$ and $\bm{\Psi}^{\smash{*}}\bm{B}\in\mathbb{F}^{p\times m}$. We report the minimum-Frobenius-norm solution to the sketched least-squares problem:

[TABLE]

See Algorithm˜3 for numerically stable pseudocode. This algorithm requires sketching both $\bm{A}$ and $\bm{B}$, plus $\mathcal{O}(d^{2}(p+m))$ arithmetic operations. See Table˜1.

To analyze the sketch-and-solve procedure, we make a change of variables. Construct a reduced SVD $\bm{A}=\bm{Q}\bm{\Sigma}\bm{V}^{\smash{*}}$, where the square diagonal matrix $\bm{\Sigma}$ is invertible. After the transformation $\bm{\Sigma}\bm{V}^{{\smash{*}}}\bm{X}\mapsto\bm{Z}$, the original least-squares problem ˜2.5 and the sketched-least-squares problem ˜2.6 become

[TABLE]

The (unique) solution $\smash{{\tilde{\smash{\boldsymbol{\bm{Z}}}\mathstrut\mkern 4.0mu}\mkern-4.0mu}{}}$ to ˜2.8 induces the minimum-norm solution $\smash{{\tilde{\smash{\boldsymbol{\bm{X}}}\mathstrut\mkern 4.0mu}\mkern-4.0mu}{}}=\bm{V}\bm{\Sigma}^{-1}\smash{{\tilde{\smash{\boldsymbol{\bm{Z}}}\mathstrut\mkern 4.0mu}\mkern-4.0mu}{}}$ to the sketched least-squares problem ˜2.6. The next lemma provides an expression for the difference between the solutions of ˜2.7 and ˜2.8. The argument is adapted from [118, Lemma A.4], where it appears in the analysis of the generalized Nyström approximation.

Lemma 2.6** (Sketch-and-solve error formula).**

Instate the prevailing notation. Let $\bm{Q}_{\perp}$ be an orthonormal matrix that satisfies $\bm{Q}^{\vphantom{{\smash{*}}}}_{\perp}\bm{Q}_{\perp}^{\smash{*}}=\mathbf{I}-\bm{Q}\bm{Q}^{\smash{*}}$. Define the matrices $\bm{\Psi}_{1}\;{\vcentcolon=}\;\bm{\Psi}^{\smash{*}}\bm{Q}$ and $\bm{\Psi}_{2}\;{\vcentcolon=}\;\bm{\Psi}^{\smash{*}}\bm{Q}_{\perp}$. Provided that $\bm{\Psi}_{1}$ has full column rank,

[TABLE]

Proof.

Starting from the definition ˜2.8 of $\smash{{\tilde{\smash{\boldsymbol{\bm{Z}}}\mathstrut\mkern 4.0mu}\mkern-4.0mu}{}}$, we calculate that

[TABLE]

This is the advertised statement. ∎

With this lemma in place, we can now prove the central result on sketch-and-solve. Several prior works [123, 45, 94, 35] have obtained similar guarantees under alternative hypotheses on $\bm{\Psi}$. The standard analysis assumes that $\bm{\Psi}$ is an OSE for the range of the concatenated matrix $[\bm{A}~\bm{B}]$. Parulekar et al. [94] assume that $\bm{\Psi}$ is an isotropic matrix that is also an OSE for the range of $\bm{A}$. Drineas et al. [45] assume that $\bm{\Psi}$ is a subspace injection for the range of $\bm{A}$ and that $\bm{\Psi}$ has an approximate matrix multiplication property. We are not aware of any sources which recognize that the OSI property is enough.

Theorem 2.7** (Sketch-and-solve with an OSI).**

Fix matrices $\bm{A}\in\mathbb{F}^{n\times d}$ and $\bm{B}\in\mathbb{F}^{n\times m}$. Draw a random $(d,\alpha)$-OSI test matrix $\bm{\Psi}\in\mathbb{F}^{n\times p}$ where $p\geq d$, and let $\smash{{\tilde{\smash{\boldsymbol{\bm{X}}}\mathstrut\mkern 4.0mu}\mkern-4.0mu}{}}\in\mathbb{F}^{d\times m}$ be the output ˜2.6 of the sketch-and-solve procedure. Then, with probability at least $\tfrac{9}{10}$,

[TABLE]

Proof.

We adopt the notation and the reparameterization outlined above. The goal is to show that

[TABLE]

The optimal least-squares residual $\bm{Q}\bm{Z}_{\star}-\bm{B}$ is orthogonal to the range of $\bm{Q}$, so

[TABLE]

The last identity follows from Lemma˜2.6. Finally, apply Lemma˜2.1 to deduce that

[TABLE]

where the bound holds with probability at least $\frac{9}{10}$. ∎

2.5 Generalized Nyström approximation

Last, we treat low-rank approximation of a matrix via the generalized Nyström method [124, 34, 87, 118].

Fix an (arbitrary) input matrix $\bm{A}\in\mathbb{F}^{n\times d}$ and an approximation rank parameter $k\leq\min\{n,d\}$. We draw two random test matrices $\bm{\Omega}\in\mathbb{F}^{d\times k}$ and $\bm{\Psi}\in\mathbb{F}^{n\times p}$ where the embedding dimensions $k\leq p\leq\min\{n,d\}$. Collect two sketches of the input matrix: $\bm{Y}\coloneqq\bm{A}\bm{\Omega}\in\mathbb{F}^{n\times k}$ and $\bm{X}\coloneqq\bm{\Psi}^{\smash{*}}\bm{A}\in\mathbb{F}^{p\times d}$. The generalized Nyström approximation takes the form

[TABLE]

We present two numerically stable implementations of ˜2.9 that prepare the approximation in outer-product form (Algorithm˜4) or in compact SVD form (Algorithm˜5). Each algorithm employs two sketches of the input matrix, as well as a sketch of intermediate data. The outer-product algorithm requires ${\mathcal{O}}(nk^{2}+dpk)$ extra arithmetic, while the SVD algorithm requires ${\mathcal{O}}(nk^{2}+dp^{2})$ extra arithmetic. Since $p\geq k$, both algorithms are more efficient when $\bm{A}$ is wide ($n\geq d$); otherwise, it is faster to apply them to the adjoint $\bm{A}^{\smash{*}}$. See Table˜1 for a summary of arithmetic costs.

To motivate the construction of the generalized Nyström approximation ˜2.9 and our analysis, we adopt the perspective from [123, 118, 87] that the RSVD approximation ˜2.3 can be written as

[TABLE]

We recognize the matrix $(\bm{A}\bm{\Omega})^{\smash{\dagger}}\bm{A}$ as a solution to the least-squares problem

[TABLE]

The generalized Nyström method can be derived by approximating the solution of this least-squares problem with sketch-and-solve, resulting in a further approximation

[TABLE]

Hence, the generalized Nyström approximation ˜2.9 is a proxy for the RSVD approximation ˜2.10.

With this observation, the error bound for the generalized Nyström approximation is an immediate corollary of our analyses of RSVD (Theorem˜2.2) and sketch-and-solve (Theorem˜2.7). Our argument parallels [118, Thm. 4.3], although the cited work only treats Gaussian test matrices.

Corollary 2.8** (Generalized Nyström with OSIs).**

Fix an input matrix $\bm{A}\in{\mathbb R}^{n\times d}$ and a target rank $r\leq\min\{n,d\}$. Draw a random $(r,\alpha)$-OSI test matrix $\bm{\Omega}\in{\mathbb R}^{d\times k}$ and a random $(k,\alpha)$-OSI test matrix $\bm{\Psi}\in{\mathbb R}^{n\times p}$ where $r\leq k\leq p$. Construct the generalized Nyström approximation $\widehat{\bm{A}}\in\mathbb{F}^{n\times d}$ determined by ˜2.9. Then, with probability at least $\frac{4}{5}$,

[TABLE]

Proof.

Using the interpretation ˜2.11 of generalized Nyström as the solution to a sketched least-squares problem, the bound follows from our results for sketch-and-solve and RSVD with OSIs:

[TABLE]

The sketch-and-solve bound and the RSVD bound each holds with probability $\tfrac{9}{10}$, so they hold simultaneously with probability $\tfrac{4}{5}$. ∎

Remark 2.9* (Truncated generalized Nyström).*

For a class of $(r,\alpha)$-OSI matrices, the ratio $\gamma\coloneqq k/r$ of the embedding dimension to the subspace dimension is called the oversampling factor. For the generalized Nyström algorithm, Corollary˜2.8 requires the test matrices $\bm{\Omega}$ and $\bm{\Psi}$ to have embedding dimensions $k=\gamma r$ and $p=\gamma k=\gamma^{2}r$, respectively. In case the oversampling parameter $\gamma={\mathcal{O}}(1)$, as with the SparseStack (Definition˜1.7) and SparseRTT (Definition˜1.10) test matrices, the larger embedding dimension still satisfies $p={\mathcal{O}}(r)$. However, for test matrices with oversampling $\gamma\gg 1$, such as an SRTT [111] or a Khatri–Rao construction (Definition˜1.12), the larger embedding dimension $p=\gamma^{2}r$ may be prohibitively expensive. In theory, we can avoid this cost by using the truncated generalized Nyström approximation [123, Thm. 4.3], which we can analyze by adapting the proof of Corollary˜2.8. This algorithm allows the test matrix $\bm{\Psi}$ to have embedding dimension $p={\mathcal{O}}(\gamma r)$. In practice, we cannot endorse this algorithm because the truncation step severely impairs the quality of the resulting low-rank approximation. Designing an algorithm with the theoretical guarantees of truncated generalized Nyström estimator and the accuracy of the standard generalized Nyström estimator is an open problem.

3 Sparse test matrices

Section˜2 shows that we can implement several sketching algorithms with OSIs. As a first example, we study the implications for randomized linear algebra with sparse test matrices. Sparse test matrices are fundamental because they allow us to design sketching methods whose runtime is (almost) linear in the sparsity of the input matrix [35, 90, 36]. In practice, sparse embeddings exhibit superb empirical performance for both sparse and dense linear algebra [46, 48, 41, 31, 74, 119]. See Figure˜3.

This section supplements the discussion of sparse test matrices from Section˜1.5.1 with additional information about theoretical guarantees, implementation, and empirical performance. Section˜3.1 reviews the SparseStack construction; Section˜3.2 presents alternative constructions; and Section˜3.3 provides an empirical study.

3.1 Randomized linear algebra with SparseStacks

We highly recommend the SparseStack test matrix, introduced in Definition˜1.7, for randomized linear algebra. This section summarizes the construction, the content of our analysis, and the overall benefits.

The parameters of a SparseStack test matrix are the ambient dimension $d$, the row sparsity $\zeta$, and the block size $b$; the embedding dimension $k\coloneqq b\zeta$. The matrix takes the form

[TABLE]

In this formula, $\mathbf{e}_{i}\in\mathbb{F}^{b}$ refers to the $i$th standard basis vector. For an input matrix $\bm{A}\in\mathbb{F}^{n\times d}$, the sketch $\bm{A}\bm{\Omega}$ requires ${\mathcal{O}}(\zeta\cdot\operatorname{nnz}(\bm{A}))$ operations because $\bm{\Omega}$ has $\zeta$ nonzero entries per row.

In the complex field ($\mathbb{F}={\mathbb C}$), we recommend switching to the distributions $\varrho_{ij}\sim\textsc{steinhaus}$ or $\varrho_{ij}\sim\textsc{rademacher}_{\mathbb C}$ for greater reliability. In all these cases, our theoretical guarantees remain valid.

Imported Theorem˜1.8, adapted from [115, Rem. 6.5], states that the SparseStack serves as an $(r,\nicefrac{{1}}{{2}})$-OSI for some embedding dimension $k={\mathcal{O}}(r)$ and row sparsity $\zeta={\mathcal{O}}(\log r)$. When we implement Algorithms˜1, 2, 3, 4, and 5 using SparseStacks with these parameter choices, we can certify that the algorithms run faster than previously known. Tables˜2 and 3 provide a comparison with existing work.

In addition to the OSI guarantee and the superior runtime bounds, the SparseStack test matrix has several other merits (and one demerit):

($+$) Construction.

We form the SparseStack by drawing and storing $2d\zeta$ simple random variables. The storage cost is fixed in advance by the parameters, which aids with memory management. In addition, the SparseStack is extensible, which means that we can increase the ambient dimension $d$ on the fly, just by tacking on additional independent rows of the same form.

($+$) Reliability.

The SparseStack is extremely reliable in practice. Numerical experiments in Figure˜2 indicate that constant row sparsity $\zeta=4$ and proportional embedding dimension $k=2r$ suffice to obtain an OSI with subspace dimension $r$. See Section˜3.3 for evidence, which gives credence to ˜1.9.

($-$) Library support.

SparseStack matrices are practical only when they are implemented with high-quality sparse linear algebra libraries. Sadly, these software libraries are less widely available than one would hope [1].

When fast sparse linear algebra primitives are not available, we recommend the SparseRTT test matrix (Definition˜1.10) as an alternative to the SparseStack.

3.2 Alternative sparse test matrices

The literature contains many constructions for sparse dimension reduction maps; for instance, see [33, 63, 37, 115, 35, 78, 90]. This section reviews two natural alternatives: SparseUniform (Section˜3.2.1) and SparseIID (Section˜3.2.2). We state reasons for (and against) preferring the SparseStack construction.

3.2.1 The SparseUniform test matrix

By and large, practitioners [84, 48, 31, 41, 119] have worked with the SparseUniform test matrix, also proposed by Kane & Nelson [63, Fig. 1(b)].

Definition 3.1** (SparseUniform).**

Fix the ambient dimension $d$, embedding dimension $k$, and row sparsity $\zeta$. The SparseUniform test matrix $\bm{\Omega}\in\mathbb{F}^{d\times k}$ is a random sparse matrix with exactly $\zeta$ independent copies of the random variable $\zeta^{-1/2}\varrho$ in each row, placed in uniformly random positions, sampled without replacement, where $\varrho\sim\textsc{rademacher}$.

Figure˜4 illustrates how the sparsity patterns of the SparseUniform and SparseStack constructions differ. In our experience, the SparseUniform and SparseStack test matrices have similar injectivity $\alpha$ and dilation $\beta$, and they attain similar accuracy on downstream tasks, such as low-rank approximation.

Although the performance is comparable, we prefer the SparseStack test matrix for two reasons. First, the SparseUniform construction involves repeated sampling of $\zeta$ indices without replacement, which is hard to implement efficiently [31, 46, 84]. Second, we presently have stronger theoretical guarantees for the SparseStack construction. Indeed, the best OSI result for a SparseStack (Imported Theorem˜1.8) allows for an embedding dimension $k={\mathcal{O}}(r)$ and a row sparsity $\zeta={\mathcal{O}}(\log r)$, while SparseUniform admits an OSE guarantee [36] for some embedding dimension $k={\mathcal{O}}(r\log r)$ and row sparsity $\zeta={\mathcal{O}}(\log r)$.

3.2.2 The SparseIID test matrix

Another alternative is the SparseIID model, a construction with roots in the work of Achlioptas [2].

Definition 3.2** (SparseIID).**

Fix the ambient dimension $d$, embedding dimension $k$, and the expected row sparsity $\zeta\in(0,k]$. Define a random variable $\varphi\coloneqq\zeta^{-1/2}\varrho\delta$, where $\varrho\sim\textsc{rademacher}$ and $\delta\sim\textsc{bernoulli}(\zeta/k)$. The entries of the SparseIID test matrix $\bm{\Omega}\in\mathbb{F}^{d\times k}$ are iid copies of $\varphi$.

Figure˜4 illustrates the sparsity pattern of the SparseIID test matrix. Tropp [115, Thm. 6.3] proved that the SparseIID model admits an OSI guarantee when the expected row sparsity $\zeta$ is adapted to the subspace coherence, defined in Section˜1.8.

Imported Theorem 3.3** (SparseIID matrices are OSIs, [115, Thm. 6.3]).**

Fix an orthonormal matrix $\bm{Q}\in{\mathbb F}^{d\times r}$, with coherence $\mu(\bm{Q})\in[\nicefrac{{r}}{{d}},1]$. A SparseIID test matrix $\bm{\Omega}\in\mathbb{F}^{d\times k}$ serves as an $(r,\nicefrac{{1}}{{2}})$-OSI for some embedding dimension $k={\mathcal{O}}(r)$ and some expected row sparsity $\zeta={\mathcal{O}}(\mu(\bm{Q})\log r)$.

Because of the role of the subspace coherence, Imported Theorem˜3.3 for SparseIID provides a stronger guarantee than Imported Theorem˜1.8 offers for SparseStack.

Nevertheless, we deprecate the use of the SparseIID construction because of three significant flaws. First, it is hard to develop efficient implementations of the iid sampling procedure that have expected runtime ${\mathcal{O}}(d\zeta)$. Second, the storage costs are variable. Third, the SparseIID model cannot serve as an OSI with constant expected row sparsity $\zeta={\mathcal{O}}(1)$ because there is a high probability that the random matrix has rows that equal zero; this failure mode is visible in Figure˜5. As a consequence, we believe that the SparseStack construction remains preferable.

3.3 Numerical evidence

Figure˜5 charts the error in the RSVD approximation (Algorithm˜1) using the SparseStack test matrix with a range of row sparsity parameters $\zeta$ and embedding dimensions $k$. The SparseStack with $\zeta=1$ agrees with the CountSketch test matrix [35, 78, 90, 27]. The plot also includes comparisons with SparseIID test matrices and the Gaussian baseline. The input matrices $\bm{A}$ are diagonal matrices from the testbed in [117, 118], as these instances are adversarial for sparse test matrices [91]. The experiments suggest that constant row sparsity $\zeta=4$ is sufficient for the SparseStack test matrix to achieve the same low-rank approximation error as a Gaussian test matrix. This evidence supports our surmise about the reliability of the SparseStack with constant row sparsity (˜1.9).

For more targeted evidence, Figure˜2 tracks numerical approximations of the injectivity and dilation parameters of a SparseStack matrix as the subspace dimension $r$ varies and the embedding dimension $k=2r$. We estimate both parameters by applying the test matrix to the adversarial orthonormal matrix $\bm{Q}=[\mathbf{e}_{1}~\cdots~\mathbf{e}_{r}]\in{\mathbb R}^{d\times r}$, which is likely the worst-case instance. When $\zeta\geq 4$, the injectivity of SparseStack appears to plateau as $r$ increases. When $\zeta=8$, the SparseStack injectivity almost matches the asymptotic injectivity $\alpha=(1-1/\sqrt{2})^{2}$ of a Gaussian test matrix (dashed line). However, for all values of $\zeta$, the SparseStack dilation deviates from the asymptotic dilation $\beta=(1+1/\sqrt{2})^{2}$ of a Gaussian. This disparity strongly suggests that the SparseStack matrix with constant row sparsity behaves as an OSI—but not as an OSE.

In Section˜6.1, we showcase a scientific application of SparseStack matrices for the compression and analysis of a numerical simulation of a Bose–Einstein condensate. In this context, the generalized Nyström method runs up to 12$\times$ faster with SparseStack test matrices, as compared with the baseline cost using Gaussian matrices.

4 SparseRTT: An optimized fast trigonometric transform test matrix

This section elaborates on the SparseRTT test matrix (Definition˜1.10). Trigonometric transform-based test matrices, such as the SparseRTT, are widely used in practice [118, 87, 88], and they are valuable for computing environments that lack high-quality sparse arithmetic libraries.

Section˜4.1 summarizes the design and analysis of SparseRTTs; Section˜4.2 contains an accounting of their computational cost; and Section˜4.3 offers implementation guidance. Section˜4.4 outlines several related constructions, along with numerical comparisons. Last, Section˜4.5 contains a proof of the main theorem on SparseRTTs.

4.1 The SparseRTT construction and its analysis

In Section˜1.5.2, we introduced random test matrices based on fast trigonometric transforms. For ambient dimension $d$ and embedding dimension $k$, the standard construction takes the form

[TABLE]

Among constructions of this type, we recommend using the SparseRTT test matrix. The SparseRTT has the fastest possible sketching time in the class ˜4.1, and it enjoys the $(r,\nicefrac{{1}}{{2}})$-OSI property at a proportional embedding dimension $k={\mathcal{O}}(r)$.

In the SparseRTT construction (Definition˜1.10), we populate the random diagonal matrix $\bm{D}$ with iid rademacher random variables; the trigonometric transform $\bm{F}$ is a WHT or DCT (when $\bm{F}\in\{{\mathbb R},{\mathbb C}\}$) or a DFT (when ${\mathbb F}={\mathbb C}$); the sparse matrix $\bm{S}$ follows the SparseCol distribution, defined below.

Definition 4.1** (SparseCol).**

Fix the dimension $d$, the embedding dimension $k$, and the column sparsity $\xi\in{\mathbb N}$. Construct a random sparse vector

[TABLE]

The random variables $\varrho_{1},...,\varrho_{\xi}\sim\textsc{rademacher}$ iid, and the selectors $s_{1},...,s_{\xi}$ are sampled uniformly from $\{1,\dots,d\}$ without replacement. As usual, $\mathbf{e}_{i}\in\mathbb{F}^{d}$ is the $i$th standard basis vector. Up to scaling, the columns of the SparseCol matrix $\bm{S}\in\mathbb{F}^{d\times k}$ are iid copies of the random sparse vector $\boldsymbol{\omega}$:

[TABLE]

Our main result contains an OSI guarantee for the SparseRTT test matrix.

Theorem 4.2** (SparseRTTs are OSIs).**

The SparseRTT test matrix $\bm{\Omega}\in\mathbb{F}^{d\times k}$ is an $(r,\nicefrac{{1}}{{2}})$-OSI for some embedding dimension $k={\mathcal{O}}(r)$ and some column sparsity \smash{\xi={\mathcal{O}}\big{(}\big{(}1+\tfrac{\log d}{r}\big{)}\cdot\log r\big{)}}. With these parameters, when the input matrix $\bm{A}\in\mathbb{F}^{n\times d}$, the sketch $\bm{A}\bm{\Omega}$ has arithmetic cost ${\mathcal{O}}(nd\log r)$.

We justify the bound on the sketching cost in Section˜4.2; the proof of the OSI guarantee appears in Section˜4.5. An example of Tropp [111, Sec. 3.3] implies that a column sparsity $\xi={\mathcal{O}}(\log r)$ is optimal for the SparseRTT construction at embedding dimension $k={\mathcal{O}}(r)$. For trigonometric transform embeddings of the form ˜4.1 with embedding dimension $k={\mathcal{O}}(r)$, the best available results [33, 32] establish the OSE property with a weaker column sparsity guarantee $\xi={\mathcal{O}}(\log^{3}r)$ using a slightly different sparse matrix $\bm{S}$.

4.2 Arithmetic costs

This section justifies the sketching cost for the SparseRTT, stated in Theorem˜4.2. Consider a SparseRTT matrix $\bm{\Omega}\in\mathbb{F}^{d\times k}$ with embedding dimension $k\leq d$. It suffices to work out the cost of forming $\bm{y}=\bm{\Omega}^{\smash{*}}\bm{x}$ for a single vector $\bm{x}\in\mathbb{F}^{d}$. A direct implementation using a standard fast trigonometric transform requires ${\mathcal{O}}(d\log d)$ operations. Here is an accounting of the cost:

Form $\bm{y}_{1}\leftarrow\bm{D}^{\smash{*}}\bm{x}$. Since $\bm{D}$ is diagonal, this step expends ${\mathcal{O}}(d)$ operations. 2. 2.

Compute $\bm{y}_{2}\leftarrow\bm{F}^{{\smash{*}}}\bm{y}_{1}$ using a fast trig transform algorithm, at a cost of ${\mathcal{O}}(d\log d)$ operations. 3. 3.

Evaluate $\bm{y}\leftarrow\bm{S}^{\smash{*}}\bm{y}_{2}$. This step requires ${\mathcal{O}}(k\xi)$ operations.

Assuming the optimal column sparsity $\xi={\mathcal{O}}(\log k)$, we obtain a final operation count of ${\mathcal{O}}(d\log d)$ Practitioners almost always use the implementation just described.

To optimize the asymptotic cost, we can employ a fast subsampled trigonometric transform. See [109] for the subsampled DFT and DCT; see [4] for the subsampled WHT. We adjust steps 2 and 3:

1. 2’.

Let ${\mathcal{R}}$ be the set of nonzero rows of $\bm{S}$. Compute $\bm{y}_{2}^{\prime}\leftarrow(\bm{F}(\mathrel{\mathop{\ordinarycolon}},{\mathcal{R}}))^{\smash{*}}\bm{y}_{1}$ using the fast subsampled trigonometric transform. The cost is ${\mathcal{O}}(d\log k)$ operations. 2. 3’.

Evaluate $\bm{y}\leftarrow(\bm{S}({\mathcal{R}},\mathrel{\mathop{\ordinarycolon}}))^{{\smash{*}}}\bm{y}_{2}^{\prime}$ at a cost of ${\mathcal{O}}(k\xi)$ operations.

Indeed, the set ${\mathcal{R}}$ contains at most $k\xi$ coordinates. Assuming the optimal column sparsity $\xi={\mathcal{O}}(\log k)$, the total cost is ${\mathcal{O}}(d\log k)$ operations, consistent with Theorem˜4.2.

4.3 Implementation guidance

This section offers some advice on the proper implementation of SparseRTT test matrices. First, let us discuss the distribution of the iid entries of the random diagonal matrix $\bm{D}$. We remark that real Rademacher variables are effective, regardless of the field $\mathbb{F}$. Unfortunately, real and complex Rademacher random variables both have failure modes for discrete data; cf. [84, Sec. 2.4]. It is often preferable to employ a continuous distribution, such as $\textsc{uniform}[-\sqrt{3},+\sqrt{3}]$ when $\mathbb{F}={\mathbb R}$ or the Steinhaus distribution when $\mathbb{F}={\mathbb C}$. Theoretical results for SparseRTTs with these distributions are the same with the results for real Rademachers, modulo constant factors.

For the distribution of entries in the sparse matrix $\bm{S}$, we can make similar recommendations. When ${\mathbb F}={\mathbb C}$, we endorse the use of $\textsc{rademacher}_{\mathbb C}$ or steinhaus variables in place of real Rademachers. The complex distributions increase the reliability and yield equivalent theoretical guarantees.

For the trigonometric transform $\bm{F}$, we generally recommend the DCT when $\mathbb{F}={\mathbb R}$ and the DFT when $\mathbb{F}={\mathbb C}$. Although the WHT is simpler arithmetically than the DCT or DFT, it is only defined when the dimension $d=2^{m}$. In practice, the cost of applying the SparseRTT is usually dominated by the cost of the standard fast trigonometric transform. While subsampled trigonometric transforms are potentially faster, high-quality implementations of these algorithms are not available at the time of writing.

To construct the sparse matrix $\bm{S}$, the user must select an explicit value for the column sparsity $\xi$. We recommend $\xi=\lceil 1.5\log k\rceil$ to avoid worst-case examples identified by Tropp [111, Rem. 1.5]. When the subspace dimension $r$ is known, we generally set the embedding dimension $k=2r$. The numerical evidence in Figure˜6 suggests that more aggressive parameter choices can be effective.

In the available runtime comparisons [41, 46, 31], sparse test matrices are significantly faster than test matrices based on trigonometric transforms. On the other hand, to achieve competitive performance, sparse test matrices require a high-quality sparse linear algebra library, and it may be necessary to implement the sketching method in a low-level programming language, such as C. Absent a sparse linear algebra library, the SparseRTT construction can be the better alternative. Although the SparseRTT involves a sparse matrix multiplication, the matrix $\bm{S}$ is so sparse that special libraries are unnecessary.

4.4 Alternative constructions

The literature describes many trigonometric transform test matrices that fit the template ˜4.1. These constructions make different choices for the sparse matrix $\bm{S}$. In this section, we summarize these alternatives, their theoretical guarantees (Table˜4), and their empirical performance (Figure˜6).

Here are several options for the sparse matrix $\bm{S}$ in a fast trigonometric transform ˜4.1. The tag describes the style of the line and data marker in Figure˜6.

Subsampling (solid line, open square). The classic subsampled trigonometric transform test matrix [4, 124] arises when $\bm{S}$ is a (scaled) restriction matrix. This matrix has minimal sparsity $\operatorname{nnz}(\bm{S})=k$, but it requires a larger embedding dimension $k=\Omega(r\log r)$ to serve as an OSI. 2. 2.

CountSketch (solid line, solid circle). Cartis et al. [26] choose $\bm{S}$ to be a CountSketch embedding. This test matrix is reliable, but it has higher sparsity ($\operatorname{nnz}(\bm{S})=d$) than the SparseCol. It is also more expensive, requiring ${\mathcal{O}}(nd\log d)$ operations to form the sketch $\bm{A}\bm{\Omega}$ of the matrix $\bm{A}\in\mathbb{F}^{n\times d}$. 3. 3.

SparseCol (dashed line, solid circle). The SparseRTT construction (Definition˜1.10) relies on the SparseCol matrix (Definition˜4.1). 4. 4.

SparseIID (solid line, solid square). If we replace SparseCol with the SparseIID map defined in Section˜3.2.2, our main result for the SparseRTT (Theorem˜4.2) holds with the same embedding dimension $k$ and (expected) column sparsity $\xi$. Empirically, the SparseIID map requires higher column sparsity $\xi$ than SparseCol to perform reliably.

Table˜4 summarizes these constructions and the influence on the arithmetic cost of the sketch-and-solve method (Algorithm˜3). We see that the SparseRTT test matrix simultaneously allows the minimum embedding dimension $k={\mathcal{O}}(r)$, while the SparseCol matrix $\bm{S}$ attains the minimum possible sparsity $\operatorname{nnz}(\bm{S})={\mathcal{O}}(r\log r)$. In addition, among the trigonometric transform sketches, the SparseRTT offers the fastest runtime for sketch-and-solve.

Figure˜6 provides an empirical comparison of the RSVD (Algorithm˜1) with trigonometric transform sketches ˜4.1 as we vary the sparse matrix $\bm{S}$. The left panels show the performance of when $\bm{S}$ is a SparseIID matrix, while the right panels show the recommended SparseCol matrix. For reference, all panels display the results for the existing constructions (subsampling $\bm{S}$, CountSketch $\bm{S}$, Gaussian $\bm{\Omega}$). For these examples, the recommended SparseRTT test matrix with column sparsity $\xi=4$ achieves accuracy comparable with a Gaussian test matrix. The SparseIID sampling method is less effective than the SparseCol sampling method at every value of $\xi$.

4.5 Proof of Theorem˜4.2

The proof of Theorem˜4.2 parallels the analysis of the SRTT construction with the WHT matrix. We rely on the following incoherence result, adapted from [111, Lem. 3.3].

Imported Theorem 4.3** (Randomized trig transforms create incoherence, [111, Lem. 3.3]).**

Fix an orthonormal matrix $\bm{Q}\in\mathbb{F}^{d\times r}$. Select $\bm{D}\in\mathbb{F}^{d\times d}$ and $\bm{F}\in\mathbb{F}^{d\times d}$ as in Definition˜1.10. With probability at least $1-\delta$, the coherence of the orthonormal matrix $\bm{Q}_{\rm F}\;{\vcentcolon=}\;(\bm{D}\bm{F})^{\smash{*}}\bm{Q}\in{\mathbb R}^{d\times k}$ satisfies

[TABLE]

Imported Theorem˜4.3 is based on the result [111, Lem. 3.3], which is stated under the assumptions that $\bm{D}\in\mathbb{F}^{d\times d}$ consists of real Rademacher variables and that $\bm{F}\in\mathbb{F}^{d\times d}$ is a WHT. Both conditions can be relaxed without any additional insight. To handle complex Rademacher variables, we split them into real and imaginary parts. The DCT and the DFT (like the WHT) are unitary matrices whose entries are bounded in magnitude by ${\mathcal{O}}(d^{-1/2})$. The proof only depends on these two features of the WHT.

Imported Theorem˜4.3 states that the (unsampled) random trigonometric transform $(\bm{D}\bm{F})^{\smash{*}}$ converts an arbitrary orthonormal matrix into an incoherent orthonormal matrix. Thus, to prove that a SparseRTT is an OSI, it suffices to confirm that a SparseCol matrix $\bm{S}$ acts as an injection on an incoherent subspace. This is the content of the next result.

Theorem 4.4** (SparseCol: Injectivity for incoherent matrices).**

Fix an orthonormal matrix $\bm{Q}\in\mathbb{F}^{d\times r}$ with coherence $\mu(\bm{Q})$. For some choice of embedding dimension $k={\mathcal{O}}(r)$ and column sparsity $\xi={\mathcal{O}}(\frac{d}{r}\mu(\bm{Q})\log r)$, a random SparseCol matrix $\bm{S}\in\mathbb{F}^{d\times k}$ satisfies

[TABLE]

Section˜7.3.2 contains a proof of Theorem˜4.4 based on the Gaussian comparison theorem (Section˜7.3). With this preparation, we can establish Theorem˜4.2.

Proof of Theorem˜4.2..

Introduce the orthonormal matrix $\bm{Q}_{\rm F}\;{\vcentcolon=}\;(\bm{D}\bm{F})^{\smash{*}}\bm{Q}\in\mathbb{F}^{d\times k}$. Imported Theorem˜4.3 states that this matrix has coherence $\mu(\bm{Q}_{\rm F})={\mathcal{O}}((r+\log d)/d)$ with probability at least $\nicefrac{{39}}{{40}}$. We construct a SparseCol matrix $\bm{S}\in\mathbb{F}^{d\times k}$ with column sparsity \xi={\mathcal{O}}\big{(}\big{(}1+\frac{\log d}{r}\big{)}\cdot\log r\big{)}={\mathcal{O}}\big{(}\frac{d}{r}\mu(\bm{Q}_{\rm F})\cdot\log r\big{)}. According to Theorem˜4.4, this choice of $\xi$ is sufficient to ensure that the SparseRTT acts injectively on $\operatorname{range}(\bm{Q})$:

[TABLE]

with probability at least $\nicefrac{{39}}{{40}}$. Take a union bound over the two events to complete the argument. ∎

5 Khatri–Rao test matrices and their applications

In recent years, we have witnessed an explosion of research on linear algebra problems with tensor structure [25, 5, 18, 3, 24, 62, 13, 61, 79, 81, 98, 50, 71, 99, 100]. For these problems, the simplest sketching techniques are grossly inefficient or entirely inapplicable. Instead, we must resort to special test matrices that are compatible with the tensor structure.

Section˜5.1 summarizes the construction of the Khatri–Rao test matrix, and Section˜5.2 describes some applications where tensor-structured sketching is essential. As we explain in Section˜5.3, Khatri–Rao test matrices can degrade exponentially as the tensor order increases. Sections˜5.4 and 5.5 present two theorems that quantify the exponential loss in complementary ways. Throughout, we will argue in favor of constructing Khatri–Rao test matrices with a spherical base distribution. We conclude with an application of this theory to a matrix recovery problem in Section˜5.6.

5.1 Khatri–Rao test matrices

Recall that the Kronecker matvec access model [79, 81] provides a useful abstraction for a variety of tensor applications. In this model, we can only access the input matrix $\bm{A}$ via the matvec $\bm{A}\boldsymbol{\omega}$ with a vector $\boldsymbol{\omega}$ that has Kronecker product structure:

[TABLE]

We call $d_{0}$ the base dimension, and we call $\ell$ the tensor order.

In the Kronecker matvec access model, the natural choice of embedding is a Khatri–Rao test matrix (Definition˜1.12). For an isotropic random vector $\boldsymbol{\nu}\in\mathbb{F}^{d_{0}}$, called the base distribution, the columns of the random matrix $\bm{\Omega}$ are Kronecker products of iid copies of $\boldsymbol{\nu}$. That is,

[TABLE]

Since the random vector $\boldsymbol{\nu}$ is isotropic, each column $\boldsymbol{\omega}_{i}$ of the test matrix $\bm{\Omega}$ is isotropic, and the test matrix $\bm{\Omega}$ satisfies the isotropy property (Definition˜1.2). The technical challenge is to understand when the test matrix ˜5.1 is a subspace injection.

We treat seven isotropic base distributions, three in ${\mathbb R}^{d_{0}}$ and four in ${\mathbb C}^{d_{0}}$. In the real case ${\mathbb F}={\mathbb R}$, we consider vectors with iid standard normal and Rademacher entries, as well as spherical random vectors. In the complex case ${\mathbb F}={\mathbb C}$, we consider the complex analogues of these three distributions, as well as vectors with iid Steinhaus entries. See Section˜1.8 for some background.

These base distributions yield test matrices that exhibit disparate performance, both in theory and in practice. For real arithmetic (${\mathbb F}={\mathbb R}$), we recommend using the real spherical distribution. For complex arithmetic (${\mathbb F}={\mathbb C}$), we recommend the complex spherical distributions. In support of these choices, we offer numerical evidence in Section˜5.3.1 and theoretical evidence in Sections˜5.4 and 5.5.

5.2 Applications

Khatri–Rao test matrices fill a different niche in the computational ecosystem, as compared with sparse test matrices and trigonometric transform test matrices. Indeed, for certain applications, our best option is to use Khatri–Rao sketching. Here are several examples.

Sylvester equations. For fixed matrices $\bm{B},\bm{C},\bm{Y}$, the linear equation $\bm{B}\bm{X}+\bm{X}\bm{C}=\bm{Y}$ in the matrix variable $\bm{X}$ is called a Sylvester equation. Equivalently,

[TABLE]

The solution operator $\bm{A}$ is often numerically low rank, so we can approximate it by means of the generalized Nyström method (Algorithm˜4). To sketch the solution operator, we can solve the Sylvester equation with (random) rank-one matrices $\bm{Y}$, for example, via the fADI algorithm [17]. The vectorization of a rank-one matrix has tensor-product structure (with order $\ell=2$), so this is an instance of the Kronecker matvec access model. For other applications of Khatri–Rao sketching to matrix equations, see [24]. 2. 2.

Tensor networks. When simulating entangled quantum systems, scientists often employ tensor networks to represent exponentially large matrices [93]. The tensor network format is an instance of the Kronecker matvec access model, so Khatri–Rao sketching offers an efficient method for compressing tensor networks. For examples, see the recent papers [25, 50]. 3. 3.

CP decompositions. The canonical polyadic decomposition (CPD) is a powerful tool for analyzing tensor-structured data [66, Sec. 3.5]. Algorithms for computing CPDs repeatedly solve overdetermined least-squares problems with tensor structure [66, Sec. 3.4]. Existing randomized algorithms depend on leverage-score sampling [68, 18]. Khatri–Rao sketching offers a natural alternative.

To complement these fancy examples, Section˜5.6 describes an elementary matrix recovery problem that involves Khatri–Rao test matrices. Section˜6.2 highlights a scientific application of Khatri–Rao test matrices for estimating the partition function of a quantum-mechanical system. Both applications benefit from a deeper theoretical understanding of tensor-structured sketches.

5.3 Overwhelming orthogonality: When Khatri–Rao goes bad

Figure˜1, in the introduction, provides empirical evidence that Khatri–Rao test matrices serve almost as well as Gaussian test matrices for computing low-rank approximations of sparse matrices that arise in applications. Nevertheless, tensor-structured test matrices exhibit behaviors that are fundamentally different from the constructions we studied before. To explain why, let us describe a phenomenon called overwhelming orthogonality [81]. Here is the most vivid example.

Consider the $2$-dimensional Rademacher base distribution $\boldsymbol{\nu}\sim\textsc{uniform}\{(\pm 1,\pm 1)\}\in{\mathbb R}^{2}$. Note that

[TABLE]

Form the tensor product $\boldsymbol{\omega}=\boldsymbol{\omega}^{(1)}\otimes\cdots\otimes\boldsymbol{\omega}^{(\ell)}$ where the factors $\boldsymbol{\omega}^{(i)}\sim\boldsymbol{\nu}$ iid. For the vector of ones $\bm{1}\in{\mathbb R}^{2^{\ell}}$,

[TABLE]

Note that the inner product equals zero when any one of the random vectors $\boldsymbol{\omega}^{(i)}$ is orthogonal to $\bm{1}$. It follows that the orthogonality probability is almost one:

[TABLE]

Thus, to certify that a vector in ${\mathbb R}^{2^{\ell}}$ is nonzero, it requires $\Omega(2^{\ell})$ tensor Rademacher measurements!

The same phenomenon persists for any isotropic base distribution $\boldsymbol{\nu}\in\mathbb{F}^{d_{0}}$ that is sometimes orthogonal to a fixed vector:

[TABLE]

Define the tensor product vector $\bm{q}\coloneqq\bm{a}\otimes\cdots\otimes\bm{a}$ and the tensor product distribution $\boldsymbol{\omega}=\boldsymbol{\omega}^{(1)}\otimes\cdots\otimes\boldsymbol{\omega}^{(\ell)}$ where the factors $\boldsymbol{\omega}^{(i)}\sim\boldsymbol{\nu}$ iid. Then the measurement $\langle\boldsymbol{\omega},\bm{q}\rangle$ is overwhelmingly likely to equal zero unless we repeat the measurement process an exponential number of times (with respect to the tensor order $\ell$).

Similar behavior occurs for any base distribution, even when the orthogonality probability ˜5.2 equals zero (rather than $\tau>0$). For example, suppose that the base distribution $\boldsymbol{\nu}\in{\mathbb R}^{d_{0}}$ is the real spherical distribution and that $\bm{q}=\bm{a}\otimes\cdots\otimes\bm{a}$ for any fixed unit vector $\bm{a}\in{\mathbb R}^{d_{0}}$. Meyer et al. [81, Lem. 8] proved that $|\langle\boldsymbol{\omega},\bm{q}\rangle|\leq\tau^{\ell}$ with probability at least $1-\tau^{\ell}$ for some $\tau\in(0,1)$. As a consequence, Khatri–Rao embeddings have trouble picking up information about vectors with Kronecker structure. For these worst-case instances, we must anticipate an exponential dependence on the tensor order $\ell$ in either the embedding dimension $k$ or in the injectivity parameter $\alpha$.

5.3.1 Numerical demonstration

Figure˜7 illustrates the problem of overwhelming orthogonality and compares several base distributions. In this experiment, we fix an orthonormal matrix $\bm{Q}\in\mathbb{F}^{d_{0}^{\ell}\times r}$ whose columns have Kronecker product structure. The base dimension $d_{0}=2$, the tensor order $\ell=10$, and the subspace dimension $r=50$. We construct a Khatri–Rao embedding $\bm{\Omega}\in\mathbb{F}^{d_{0}^{\ell}\times k}$ with each of the seven base distributions $\boldsymbol{\nu}$. The left panel charts the injectivity $\sigma_{\min}^{2}(\bm{\Omega}^{\smash{*}}\bm{Q})$ as a function of the embedding dimension $k$; the right panel shows the relative error in the RSVD approximation of the matrix $\bm{A}=\bm{Q}\bm{Q}^{\smash{*}}+10^{-4}\mathbf{I}$.

Since the columns of the input matrix $\bm{Q}$ have Kronecker structure, overwhelming orthogonality manifests fiercely. For all base distributions, observe that embedding dimension $k\gg r$ is required to achieve injectivity similar with that of a Gaussian test matrix. Nevertheless, we can achieve competitive RSVD performance using a Khatri–Rao embedding with the complex spherical base distribution. We also observe that the injectivity is very sensitive to the base distribution. In particular, the complex Rademacher distribution fails to achieve injectivity $\alpha>0$ until the embedding dimension $k\geq 200$, and the real Rademacher distribution fails catastrophically even when $k=1000$. Notably, the Steinhaus base distribution, another complex generalization of the Rademacher, avoids this failure mode entirely.

5.4 First result: Khatri–Rao embeddings with constant injectivity

In this section, we take the first step toward understanding the OSI properties of a Khatri–Rao test matrix. We obtain a sufficient condition on the embedding dimension $k$ to obtain a constant bound for the injectivity, say $\alpha\geq\tfrac{1}{2}$. The result reflects the choice of base distribution, and—because of overwhelming orthogonality—the embedding dimension depends exponentially on the tensor order $\ell$.

Theorem 5.1** (Khatri–Rao: Constant injectivity).**

Consider an isotropic base distribution $\boldsymbol{\nu}\in\mathbb{F}^{d_{0}}$ that satisfies the fourth-moment bound

[TABLE]

For tensor order $\ell$, the Khatri–Rao test matrix $\bm{\Omega}\in\mathbb{F}^{d_{0}^{\ell}\times k}$ enjoys the $(r,\nicefrac{{1}}{{2}})$-OSI property for some embedding dimension $k={\mathcal{O}}({\mathrm{C}}_{\boldsymbol{\nu}}^{\ell}r)$.

The proof of Theorem˜5.1 appears in Section˜7.1. The argument depends on the fourth-moment theorem of Oliveira [92, Thm. 1.1] and variance bounds for Kronecker quadratic forms due to Meyer & Avron [79]. For the embedding dimension columns, we omit ${\mathcal{O}}(\cdot)$ symbols for visual clarity.

Meyer & Avron [79] have calculated the exact value of the fourth-moment constant $\mathrm{C}_{\boldsymbol{\nu}}$ for each of several base distributions; Table˜5 summarizes their results. The relative ordering of the fourth-moment constants is wholly consistent with the numerical experiments on injectivity in Figure˜7.

In contrast, the RSVD error induces a different ordering of the base distributions in Figure˜7. The reason for the discrepancy is that RSVD (Algorithm˜1) and psd Nyström (Algorithm˜2) only depend on the range of the test matrix $\bm{\Omega}$, so Gaussian and spherical base distributions result in the same performance in exact arithmetic. On the other hand, sketch-and-solve (Algorithm˜3) and generalized Nyström (Algorithm˜4) are sensitive to the scaling of columns of the test matrix $\bm{\Omega}$, so spherical distributions are the best choice for general-purpose use.

In quantum science, we frequently work in the regime where the base dimension $d_{0}=2$ and the tensor order $\ell$ is very large [93, 50]. In this case, real standard normal vectors have the worst predicted embedding dimension bound $k={\mathcal{O}}(3^{\ell}r)$, and the complex spherical distribution has the best embedding dimension bound $k={\mathcal{O}}(\smash{1.\overline{33}{}^{\ell}r})$. It is remarkable that changes in the scalar field and the base distribution can lead to exponential changes in the embedding dimension of a Khatri–Rao test matrix (and the downstream cost of linear algebra algorithms).

5.5 Second result: Khatri–Rao embeddings with proportional embedding dimension

Theorem˜5.1 requires the embedding dimension $k$ of a Khatri–Rao test matrix to scale exponentially with the tensor order $\ell$. This cost is acceptable for small $\ell$, but we can rarely afford the exponential scaling when $\ell$ is large. Instead, we would prefer to employ a Khatri–Rao test matrix whose embedding dimension is proportional to the subspace dimension: $k={\mathcal{O}}(r)$, without dependence on $\ell$.

To achieve this goal, we must avoid base distributions that have a zero orthogonality probability, as in ˜5.2. Therefore, we restrict our attention to continuous base distributions that have an anticoncentration property. The resulting Khatri–Rao embedding achieves the linear scaling of the embedding dimension with the subspace dimension. Nevertheless, as predicted by overwhelming orthogonality, we must accept an exponentially small injectivity parameter.

Theorem 5.2** (Khatri–Rao: Proportional embedding dimension).**

Consider a real, isotropic base distribution $\boldsymbol{\nu}\in{\mathbb R}^{d_{0}}$ whose entries are independent, continuous random variables, with densities uniformly bounded by $M$. The associated Khatri–Rao embedding $\bm{\Omega}\in{\mathbb R}^{d_{0}^{\ell}\times k}$ serves as an $(r,\alpha)$-OSI with injectivity $\alpha=(\mathrm{C}M)^{-2\ell}$ at an embedding dimension $k={\mathcal{O}}(r)$. The same result holds with $M=1$ when the base distribution $\boldsymbol{\nu}$ is isotropic and log-concave, or when $\boldsymbol{\nu}$ is the real spherical distribution.

We prove Theorem˜5.2 using a version of the small-ball method [120, Exer. 8.27]. The argument requires small-ball probability bounds for tensor products due to Hu & Paouris [59]. See Section˜7.2 for details.

In particular, a Khatri–Rao embedding with a real standard normal base distribution is an $(r,\alpha)$-OSI with injectivity $\alpha=\mathrm{C}^{-\ell}$ at an embedding dimension $k={\mathcal{O}}(r)$. Recall that the randomized linear algebra algorithms (Section˜2) produce errors that are within a $\mathrm{poly}(\alpha^{-1})$ factor of optimal. As a consequence, Gaussian Khatri–Rao embeddings with $k=\mathcal{O}(r)$ can yield errors that are suboptimal by an exponential factor $\mathrm{C}^{\ell}$. This debility is visible in Figure˜7 (left) when the embedding dimension $k\approx r=50$.

Observe that discrete base distributions, such as Rademacher distributions, do not meet the hypotheses of Theorem˜5.2. This limitation reflects a genuine failure mode of the Khatri–Rao embedding, arising from the overwhelming orthogonality phenomenon outlined in Section˜5.3. For a discrete base distribution supported on a finite number of points, the associated Khatri–Rao test matrix can never achieve a proportional embedding dimension $k={\mathcal{O}}(r)$ that is independent of the tensor order $\ell$. See [81, Thm. 19] for confirmation of this point.

In practice, we can often avoid the exponential scaling of the Khatri–Rao embedding dimension or the injectivity parameter. First, some applications of Khatri–Rao sketching involve small tensor order $\ell$. For example, we have $\ell=2$ for Sylvester equations (Section˜5.2) and for the matrix recovery problem discussed below (Section˜5.6). Second, even when $\ell$ is large, Khatri–Rao embeddings only exhibit exponential scaling for worst-case subspaces. Indeed, Figure˜1 already provides clear evidence that the RSVD with a Khatri–Rao embedding can reliably approximate matrices arising from applications.

5.6 Application: Approximate matrix recovery from bilinear queries

In recent years, matrix recovery problems have captured the imagination of many researchers in computational mathematics [8, 57, 30, 21, 20, 7]. Matrix recovery problems ask us to find a near-optimal approximation of an input matrix $\bm{B}\in\mathbb{F}^{n\times n}$ from within a class ${\mathcal{F}}\subseteq\mathbb{F}^{n\times n}$ of structured matrices:

[TABLE]

The parameter $\Gamma\geq 1$ is a factor of suboptimality. Examples of structured classes ${\mathcal{F}}$ include patterned matrices (e.g., circulant, Toeplitz, Hankel) or matrices with a fixed bandwidth.

In many setting, we can only access the target matrix $\bm{B}\in\mathbb{F}^{n\times n}$ via a restricted class of queries. For example, in the matvec access model, the algorithm can interact with the matrix by extracting the product $\bm{x}\mapsto\bm{B}\bm{x}$ with a test vector $\bm{x}\in\mathbb{F}^{n}$. The goal is to produce the approximation $\tilde{\bm{B}}$ after performing the minimum number of queries. This challenge arises in operator learning problems [21, 20].

In this paper, we consider the problem of approximating a matrix from a linearly parameterized matrix family [57, Def. 2.1]. More precisely, the class ${\mathcal{F}}$ is a $d$-dimensional linear subspace spanned by a linearly independent family of matrices $\bm{M}_{1},\dots,\bm{M}_{d}\in\mathbb{F}^{n\times n}$:

[TABLE]

We study the bilinear access model, where we interact with the input matrix $\bm{B}\in\mathbb{F}^{n\times n}$ by evaluating bilinear forms $(\bm{x},\bm{y})\mapsto\bm{y}^{\sf T}\bm{B}\bm{x}$ for input vectors $\bm{x},\bm{y}\in\mathbb{F}^{n}$ of our choice. Let us emphasize that the bilinear form involves the plain transpose, without a conjugate; this choice reflects the fact that tensor products are multilinear. This access model provides an abstraction for resource-constrained operator learning problems where we can only test the operator by applying it to elements of the computational basis, measuring the output in another computational basis.

A simple dimension-counting argument shows that $d$ bilinear queries are necessary and sufficient to perfectly reconstruct a target matrix $\bm{B}\in{\mathcal{F}}$ that belongs to the class. For an arbitrary target matrix $\bm{B}\in\mathbb{F}^{n\times n}$, we can prove that ${\mathcal{O}}(d)$ bilinear queries suffice to find a near-optimal approximation $\tilde{\bm{B}}\in{\mathcal{F}}$ from the structured class. This result is a consequence of our OSI guarantees for Khatri–Rao embeddings.

Theorem 5.3** (Approximate matrix recovery from bilinear queries).**

Let ${\mathcal{F}}$ be a $d$-dimensional linearly parametrized matrix family, as in ˜5.5, and let $\bm{B}\in{\mathbb F}^{n\times n}$ be any matrix. With $p={\mathcal{O}}(d)$ bilinear queries, Algorithm˜6 returns a matrix $\smash{{\tilde{\smash{\boldsymbol{\bm{B}}}\mathstrut\mkern 4.0mu}\mkern-4.0mu}{}}\in\mathbb{F}^{n\times n}$ that satisfies

[TABLE]

To the best of our knowledge, no one has investigated the approximate matrix recovery problem in the bilinear access model. Thus, we believe Theorem˜5.3 is the first result to achieve a nontrivial query complexity for this problem. Using the best existing OSE result for Khatri–Rao embeddings [24, Thm. 5], we can only attain the weaker query complexity bound $p={\mathcal{O}}(d^{3/2})$. The concurrent research [104] yields an improvement, but it still falls short of the optimal query complexity.

Proof.

The matrix recovery problem ˜5.4 is as an overdetermined linear least-squares problem. Define

[TABLE]

Given an input matrix $\bm{B}$, the best Frobenius-norm approximation $\bm{B}_{\star}$ from the subspace ${\mathcal{F}}$ admits a representation as a linear combination of the basis matrices: $\bm{B}_{\star}=\sum_{j=1}^{d}x_{j}^{\star}\bm{M}_{j}$. The vector $\bm{x}^{\star}\coloneqq(x_{1}^{\star},\dots,x_{d}^{\star})$ of optimal coefficients is obtained as the (unique) solution to the least-squares problem

[TABLE]

Construct the Khatri–Rao test matrix $\bm{\Omega}\in\mathbb{F}^{n^{2}\times k}$ whose columns $(\bm{x}_{i}\otimes\bm{y}_{i}\mathrel{\mathop{\ordinarycolon}}i=1,\dots,k)$ are tensor products of independent standard normal random vectors. According to Theorem˜5.1 and Table˜5, (the entrywise complex conjugate of) this test matrix satisfies the $(d,\nicefrac{{1}}{{2}})$-OSI property for some $k=\mathcal{O}(d)$.

Algorithm˜6 attacks the vectorized least-squares problem by sketching $\bm{F}=\bm{\Omega}^{\sf T}\bm{A}$ and $\bm{g}=\bm{\Omega}^{\sf T}\bm{b}$. It returns the approximate solution $\tilde{}\bm{x}=(\bm{\Omega}^{\sf T}\bm{A})^{\smash{\dagger}}(\bm{\Omega}^{\sf T}\bm{g})$. The resulting matrix solution $\tilde{\bm{B}}=\sum_{j=1}^{d}\tilde{x}_{j}\bm{M}_{j}$ satisfies

[TABLE]

The inequality follows from our analysis of sketch-and-solve with an OSI (Theorem˜2.7). ∎

6 Scientific applications

We present two scientific applications of the structured test matrices studied in this paper. In Section˜6.1, we employ the SparseStack test matrix for compression and analysis of scientific simulation data. This approach produces an accurate approximation of the input data faster than baseline algorithms that use Gaussian test matrices. We apply this methodology to data from a simulation of Bose–Einstein condensates in condensed matter physics.

In Section˜6.2, Khatri–Rao test matrices allow us to estimate the partition function of a large quantum-mechanical system, a key quantity that describes the system’s thermodynamic behavior. The tensor structure of the quantum Hamiltonian motivates us to deploy algorithms based on tensor-structured sketches. We can use Khatri–Rao test matrices to implement several variance-reduced trace estimation algorithms [80, 95, 49] that include a low-rank approximation stage. For the quantum systems that we study, these approaches achieve much better accuracy than simpler trace estimation methods.

Our experiments were run on the UC Berkeley SCF cluster, using 64 cores of a dual-socket AMD EPYC 7543 (2.8 GHz) with 500 GB of memory. All experiments were run in MATLAB v24.1 Our code is publicly available at https://github.com/chriscamano/OSI-and-Structured-Random-Matrices.

6.1 Application 1: Scientific data compression

In our first application, we use the SparseStack test matrix (Definition˜1.7) with the generalized Nyström algorithm (Algorithm˜5) to compress and analyze matrices of scientific simulation data. Our experiments document up to a 12$\times$ speedup relative to algorithms with Gaussian test matrices (Definition˜1.6).

Here is the basic setup. Consider a dynamical system $\bm{x}(t)\in{\mathbb F}^{n}$, and assemble a data matrix

[TABLE]

by collating the states of the system at $d$ equally spaced times $t_{1},\ldots,t_{d}$. The snapshots $\bm{x}(t_{i})$ are typically obtained as outputs of a numerical differential equation solver, although they could also be measurements from experiments. We are interested in two tasks:

Compression. To reduce the storage cost for the data matrix, we can compute and store a rank-$k$ approximation to $\bm{A}$, where $k\ll\min\{n,d\}$. 2. 2.

POD modes. To understand the dynamics, we can compute the $r$ dominant left singular vectors of $\bm{A}$, which are known as the proper orthogonal decomposition (POD) modes in this setting. The POD modes help us visualize the “important structures” in the dynamical system, and they provide a convenient basis for reduced-order modeling [39].

For both tasks, we can invoke the generalized Nyström approximation (Section˜2.5). For compression, we use Algorithm˜4, which returns an approximation in outer-product form $\smash{{\widehat{\smash{\boldsymbol{\bm{A}}}\mathstrut\mkern 4.0mu}\mkern-4.0mu}{}}=\bm{F}\bm{G}^{\smash{*}}$. For the POD modes, we need the output in SVD form, which is provided by Algorithm˜5. Algorithms˜4 and 5 have essentially the same asymptotic runtime, but the former is significantly faster because it avoids costly orthogonalization steps. As an aside, the same methodology applies in the streaming setting [119], where the snapshots $\bm{x}(t_{i})$ are delivered sequentially, processed immediately, and then discarded. Our paper does not include experiments with streaming data.

To demonstrate this approach, we apply it to data from the simulation of a rotating Bose–Einstein condensate (BEC), a quantum superfluid that exhibits quantum phenomena at macroscopic scales [9]. To compute the ground state of this system, condensed matter physicists often solve the imaginary-time Gross–Pitaevskii equation, a nonlinear partial differential equation whose solutions converge to a ground-state in the large-time limit. These ground states, visualized in the left panel of Figure˜8, are known to display complex dynamics, such as quantized vortices whose cores are spatial regions where the probability density for the BEC vanishes. These vortices are the “holes” visible in this plot.

Here is a summary of the experimental setup; see Appendix˜D for details. We simulate the imaginary-time Gross–Pitaevskii equation for $d=40{,}000$ timesteps on a discretized $256\times 256$ grid, resulting in a data matrix $\bm{A}\in{\mathbb C}^{256^{2}\times 40{,}000}$. Following standard practice, we center each snapshot so that each column of $\bm{A}$ has zero mean. We perform both the data compression and POD-mode tasks using generalized Nyström approximation (Algorithm˜4 or 5) with SparseStack test matrices with row sparsity $\zeta=4$.

The right-hand panels in Figure˜8 contain visualizations of the top $6$ POD modes (left singular vectors). These POD modes reveal symmetric density patterns in the condensate, isolating the dominant dynamics that emerge during imaginary-time evolution. They reflect the formation of quantized vortices in the BEC.

Figure˜9 charts the errors and runtime of this method. As compared to a baseline experiment using Gaussian test matrices (open circles), the right panel shows that SparseStack test matrices (closed squares) improve runtimes by up to $\boldsymbol{12\times}$ for compression (orange) and $\boldsymbol{4\times}$ for computation of POD modes (blue). Absolute runtimes appear in the top left panel. The bottom left panel confirms that sparse test matrices result in no loss of accuracy.

6.2 Application 2: Approximating the partition function

of a quantum system

In our second application, we exploit Khatri–Rao test matrices to approximate the partition function of a quantum-mechanical system. For this computation, we apply variance-reduced trace estimators to a tensor-structured matrix. As compared with earlier Khatri–Rao sketching methods, our approach improves the accuracy by six to twelve orders of magnitude—with minimal increase in runtime.

Let $\bm{H}\in{\mathbb R}^{2^{\ell}\times 2^{\ell}}$ denote the (symmetric) Hamiltonian matrix for a system of $\ell$ interacting quantum particles. We are interested in computing the partition function $Z(\beta)\;{\vcentcolon=}\;\operatorname{Tr}(\exp(-\beta\bm{H}))$ of the system, where $\beta>0$ is a parameter called the inverse temperature. The partition function is a central quantity in statistical mechanics that reflects the properties of a quantum system in thermal equilibrium. Because the dimension of $\bm{H}$ is exponential in the number of particles $\ell$, direct computation with $\bm{H}$ rapidly becomes infeasible.

Fortunately, most quantum Hamiltonians natively possess tensor structure, and $\bm{H}$ can often be represented using a compact tensor or tensor network format [93]. In particular, these tensor (network) formats allow us to access the matrix $\bm{H}$ via Kronecker matvecs, enabling us to efficiently sketch $\bm{H}$ with Khatri–Rao test matrices of order $\ell$. Moreover, we can apply iterative methods to approximate tensor-structured matvecs with $\bm{A}\coloneqq\exp(-\beta\bm{H})$ [6, 29]. Thus, computing the partition function $Z(\beta)$ amounts to estimating $\operatorname{Tr}(\bm{A})$ using a small number of Kronecker matvecs.

The most basic trace estimation algorithm is the Girard–Hutchinson trace estimator [53, 60]:

[TABLE]

This formula yields an unbiased estimate of $\operatorname{Tr}(\bm{A})$, and it requires $t$ matvecs with $\bm{A}$. This estimator is unbiased for any isotropic test matrix $\bm{\Omega}$ and requires $t$ matrix–vector products with $\bm{A}$ to compute. This estimator has been applied to partition function estimation in the past with ordinary (i.e., non-Kronecker) matvecs [83, 23, 107]. These earlier works do not fully exploit the tensor structure, so the methods are limited to small tensor order $\ell$. In order to reach larger values of $\ell$, we can take $\bm{\Omega}$ to be a random Khatri–Rao test matrix. The Girard–Hutchinson estimator $H_{t}$ with a Khatri–Rao test matrix provably converges as the parameter $t$ increases, but the convergence is very slow for large $\ell$ [50, 79], requiring an egregious number of matvecs to attain even one digit of accuracy; see Figure˜10.

A recent line of work studies variance-reduced algorithms for trace estimation, which augment the Girard–Hutchinson trace estimator with a low-rank approximation stage. This strategy is exemplified by the NA-Hutch++estimator [80], which returns

[TABLE]

In this expression, $\smash{{\widehat{\smash{\boldsymbol{\bm{A}}}\mathstrut\mkern 4.0mu}\mkern-4.0mu}{}}$ is a rank-$t/6$ generalized Nyström approximation of $\bm{A}$, computed with test matrices $\bm{\Omega}\in{\mathbb F}^{n\times(t/6)}$ and with $\bm{\Psi}\in{\mathbb F}^{n\times(t/3)}$. The algorithm expends $t$ matvecs to apply $\bm{A}$ and $\bm{A}^{*}$.

The variance of the NA-Hutch++ estimator equals $\operatorname*{Var}[H_{t/2}(\bm{A}-\smash{{\widehat{\smash{\boldsymbol{\bm{A}}}\mathstrut\mkern 4.0mu}\mkern-4.0mu}{}})]$, which tends to be much smaller than the variance of the Girard–Hutchinson trace estimator when the input matrix $\bm{A}$ is close to a low-rank matrix. Fortunately, this is the case for quantum Hamiltonians, where the eigenvalues of the matrix $\bm{A}=\exp(-\beta\bm{H})$ often decay rapidly.

As we showed in Section˜5.5, Khatri–Rao test matrices can serve as OSIs, so they offer guarantees for low-rank approximation, and they can be used to implement variance-reduced trace estimators, such as NA-Hutch++. In the worst case, Khatri–Rao test matrices increase errors by a factor exponential in the tensor order $\ell$. Even so, we can still hope that variance-reduced trace estimators with Khatri–Rao test matrices perform well in practice.

Figure˜10 puts this hope to the test. We deploy several variance-reduced trace estimators using real spherical Khatri–Rao test matrices with tensor order $\ell$. We estimate the partition function for the transverse field Ising model [97], a 1D quantum system whose Hamiltonian is given as

[TABLE]

Here, $\bm{X}_{i}$ and $\bm{Z}_{i}$ denote Pauli operators acting on the $i$th site and $h$ is a coupling coefficient. For simplicity, this experiment is restricted to the moderate value $\ell=16$, and our code works with an explicit representation of $\bm{H}$ as a sparse matrix. We test three variance-reduced trace estimators that are applicable to the Kronecker matvec access model, including NA-Hutch++, Nyström++ [95], and XNysTrace [49]. Appendix˜E contains further implementation details.

Figure˜10 shows the results. We see that variance-reduced estimators significantly outperform the Girard–Hutchinson estimator when using the same number of matvecs with the input matrix. On the most extreme example, variance-reduced trace estimators attain accuracy $\boldsymbol{10^{-12}}$, while the Kronecker Girard–Hutchinson estimator fails to achieve even a single digit of accuracy. Additionally, we see that NA-Hutch++, which uses generalized Nyström, has consistently worse error than the other two variance-reduced estimators, which rely on psd Nyström. In the right subplot, we even see the error in NA-Hutch++ increase with the number of matvecs. For partition function estimation, the overall evidence suggests that variance-reduced trace estimators with Khatri–Rao test matrices can be extremely effective.

7 The OSI property: Proof strategies

In Section˜2, we saw that several randomized algorithms for low-rank approximation and linear regression are effective when implemented with OSI test matrices. To complete the design of these algorithms, we must develop strategies for establishing the OSI property.

Proposition˜1.15 states that the OSI property admits an equivalent spectral formulation. As a consequence, we can establish the OSI by producing lower bounds on the minimum singular value of the test matrix, compressed to a subspace. The literature on high-dimensional probability contains several powerful techniques for bounding minimum singular values. This section elaborates on these technical methods, and it derives OSI guarantees for each of the major classes of test matrices discussed in this paper (SparseStack, SparseRTT, Khatri–Rao).

Section˜7.1 describes a method based on fourth-moment bounds, which we use to identify Khatri–Rao test matrices with constant injectivity. Section˜7.2 introduces the small-ball method, which allows us to prove that certain Khatri–Rao test matrices are subspace injections at an embedding dimension $k={\mathcal{O}}(r)$. Section˜7.3 outlines the Gaussian comparison method, which lets us analyze sparse test matrices. Last, Section˜7.4 comments on some matrix concentration tools that may be relevant for other problems.

7.1 Methodology 1: Fourth-moment bounds

The first strategy deduces the injectivity property from a moment comparison assumption. This theorem was established by Oliveira [92, Thm. 1.1] using the PAC-Bayesian method. Tropp [115, Thm. 5.2] recently obtained an alternative proof based on Gaussian comparison.

Imported Theorem 7.1** (OSI via fourth moment bound, [115, Thm. 5.2 and Rem. 5.1]).**

Let $\boldsymbol{\omega}\in\mathbb{F}^{d}$ be an isotropic random vector. Suppose that the second and fourth moments of all linear marginals of $\boldsymbol{\omega}$ are comparable up to a constant factor $h_{4}\geq 1$. That is,

[TABLE]

Consider the random test matrix

[TABLE]

Then $\bm{\Omega}$ satisfies the $(r,\alpha)$-OSI property when the embedding dimension $k\geq\mathrm{C}\cdot(1-\alpha)^{-2}\cdot h_{4}r$.

The equivalence ˜7.1 of all second and fourth moments is a rather strong assumption. It requires the random vector $\boldsymbol{\omega}$ to have controlled variance in all directions. For example, standard Gaussian and standard Rademacher vectors satisfy the comparison ˜7.1 with $h_{4}=3$. On the other hand, for a sparse Rademacher random vector in $\mathbb{F}^{d}$ with $\xi$ nonzero entries, $h_{4}\geq d/\xi$, so Imported Theorem˜7.1 gives pessimistic bounds for the embedding dimension of sparse test matrices.

7.1.1 Application: Khatri–Rao embeddings with constant injectivity

As an example of the fourth-moment theorem, we can establish Theorem˜5.1, which states that Khatri–Rao test matrices are OSIs when the embedding dimension is exponential in the order $\ell$ of the tensor product. In fact, Theorem˜5.1 follows immediately when we combine Imported Theorem˜7.1 with the following result of Meyer & Avron [79].

Imported Theorem 7.2** (Random Kronecker products: Fourth-moment bounds, [79, Thm. 4]).**

Consider an isotropic random vector $\boldsymbol{\nu}\in\mathbb{F}^{d_{0}}$ that satisfies the moment bound

[TABLE]

Construct the Kronecker product $\boldsymbol{\omega}\coloneqq\boldsymbol{\omega}^{(1)}\otimes\dots\otimes\boldsymbol{\omega}^{(\ell)}\in\mathbb{F}^{d_{0}^{\ell}}$ where the factors $\boldsymbol{\omega}^{(i)}\sim\boldsymbol{\nu}$ are iid. Then

[TABLE]

Equivalently,

[TABLE]

To recover the original formulation [79, Thm. 4] of Imported Theorem˜7.2, we confirm the equivalence of ˜7.2 and ˜7.3. By choosing a rank-one matrix $\bm{M}=\bm{u}\bm{u}^{\smash{*}}$, it is clear that the matrix formulation ˜7.3 implies the vector formulation ˜7.2. For the converse, rewrite ˜7.3 as

[TABLE]

A convex function on a compact set attains its supremum at an extreme point. In this case, each extreme point is a rank-one matrix $\bm{M}=\bm{u}\bm{u}^{{\smash{*}}}$ where $\bm{u}$ is a unit-norm vector.

For a random Kronecker product vector, Imported Theorem˜7.2 indicates that the moment comparison factor $h_{4}=\mathrm{C}_{\boldsymbol{\nu}}^{\ell}$ scales exponentially with the tensor order $\ell$. As a consequence, our result Theorem˜5.1 provides an $(r,\nicefrac{{1}}{{2}})$-OSI guarantee for a Khatri–Rao test matrix at some embedding dimension $k=\Theta(\mathrm{C}_{\boldsymbol{\nu}}^{\ell}r)$, while it offers no information about the test matrix when $k=o(\mathrm{C}_{\boldsymbol{\nu}}^{\ell}r)$.

This limitation reflects a genuine phenomenon. Consider a Khatri–Rao test matrix $\bm{\Omega}\in{\mathbb R}^{2^{\ell}\times k}$, generated from a $2$-dimensional Rademacher base distribution $\boldsymbol{\nu}$ on ${\mathbb R}^{2}$. For this base distribution, the comparison factor $\mathrm{C}_{\boldsymbol{\nu}}=2$, and Theorem˜5.1 demands an embedding dimension $k={\mathcal{O}}(2^{\ell}r)$ to achieve an $(r,\nicefrac{{1}}{{2}})$-OSI. The discussion in Section˜5.3 implies that we must choose the embedding dimension $k=\Omega(2^{\ell})$ to achieve a nontrivial OSI guarantee; see also [81, Thm. 4]. Confirming these predictions, our empirical work (Section˜5.3.1) provides clear evidence that the Rademacher Khatri–Rao test matrix fails catastrophically when $k=\Theta(r)$ and the tensor order $\ell$ is large.

7.2 Methodology 2: Small-ball probabilities

We must exploit an alternative structure to obtain guarantees for Khatri–Rao test matrices with proportional embedding dimension $k={\mathcal{O}}(r)$. To that end, we introduce an injectivity theorem based on the small-ball method. This strategy was originally developed by Mendelson and coauthors [75, 76, 67, 77]. The version here is adapted from several other sources: a result [56, Thm. 1.6] of Guédon et al., Tropp’s lecture notes [114, Thm. 19.1], and Vershynin’s textbook [120, Exer. 8.27]. Our result improves the constants and extends the bound to the complex setting.

Theorem 7.3** (OSI via small-ball method).**

Fix a constant $\tau>0$. Assume that the isotropic random vector $\boldsymbol{\omega}\in\mathbb{F}^{d}$ satisfies the small-ball probability bound

[TABLE]

Consider the random test matrix

[TABLE]

Then $\bm{\Omega}$ enjoys the $(r,\nicefrac{{\tau}}{{10}})$-OSI property when the embedding dimension satisfies

[TABLE]

The small-ball method captures a fundamental insight about injectivity. Indeed, the main obstacle to injectivity is the existence of a fixed unit-norm vector that is (almost) orthogonal to the random vector $\boldsymbol{\omega}$. The statistic $\tau$ in the small-ball condition ˜7.4 quantifies the severity of this issue.

Proof sketch.

The argument relies on two ideas. First, we can bound the minimum singular value below by introducing indicator random variables:

[TABLE]

The first identity follows from the variational definition of the minimum singular value. The lower bound depends on the numerical inequality $t\geq\tau\cdot\mathbbm{1}_{\{t\geq\tau\}}$, valid for $t\geq 0$.

Second, we regard the indicators as Boolean classifiers, which suggests an argument based on the VC dimension. The VC law of large numbers [120, Thm. 8.3.15] states that the empirical classification probability concentrates uniformly around its expectation:

[TABLE]

The family of linear classifiers in $\mathbb{F}^{d}$ has VC dimension $\Theta(d)$. With some additional arguments to handle the absolute values, we determine that $k=\Theta(d)$ samples suffice to obtain concentration simultaneously for all $\bm{u}\in\mathbb{F}^{d}$. See Appendix˜C for the details. ∎

The small-ball method yields OSI guarantees for a wide range of random test matrices. It is well-suited for problems where the random vector $\boldsymbol{\omega}$ is continuous, with a bounded density. As with Imported Theorem˜7.1, when the random vector $\boldsymbol{\omega}$ has controlled variance in every direction, the small-ball condition ˜7.4 holds. In fact, we can sometimes activate the small-ball condition in settings where the random vector $\boldsymbol{\omega}$ lacks finite polynomial moments of any order.

On the other hand, the small-ball method is ill-suited for problems where the distribution of the random vector $\boldsymbol{\omega}$ is spiky. The condition ˜7.4 never holds for a (very) sparse random vector $\boldsymbol{\omega}$. The small-ball condition also fails when $\boldsymbol{\omega}$ is a Kronecker product of iid Rademacher vectors because of the overwhelming orthogonality property (Section˜5.3).

7.2.1 Application: Khatri–Rao embeddings with proportional embedding dimension

We can deploy the small-ball technology to study Khatri–Rao test matrices. To that end, we need a small-ball probability bound for the Kronecker product of random vectors. This property is somewhat delicate, especially when the base distribution is discrete. As a consequence, it is natural to assume that the base distribution is continuous, with an anticoncentration property. Recently, Hu & Paouris [59] established a result of this type for real random vectors.

Imported Theorem 7.4** (Random Kronecker products: Small-ball probability, [59, Thms. 1.1 and 1.4]).**

Consider a real random vector $\boldsymbol{\nu}\in{\mathbb R}^{d_{0}}$ with independent entries, each drawn from an absolutely continuous distribution whose density is uniformly bounded by $M$. Construct the Kronecker product $\boldsymbol{\omega}\coloneqq\boldsymbol{\omega}^{(1)}\otimes\cdots\otimes\boldsymbol{\omega}^{(\ell)}\in{\mathbb R}^{d_{0}^{\ell}}$ where the factors $\boldsymbol{\omega}^{(i)}\sim\boldsymbol{\nu}$ iid. Then we have the small-ball probability bound

[TABLE]

Alternatively, if the distribution of $\boldsymbol{\nu}\in{\mathbb R}^{d_{0}}$ is log-concave, then the same bound holds with $M=1$.

To establish our result on Khatri–Rao products at the proportional embedding dimension (Theorem˜5.2), we simply invoke the small-ball method (Theorem˜7.3) with the small-bound probability bound from Imported Theorem˜7.4.

We can extend Imported Theorem˜7.4 to handle the real spherical base distribution. Here is a sketch of the argument. Consider the base distribution $\boldsymbol{\nu}$ that is uniform on the ball with radius $\sqrt{d_{0}}$. This base distribution is log-concave, so Imported Theorem˜7.4 provides a small-ball probability bound with $M=1$. To pass to the spherical distribution, we dilate each tensor factor $\boldsymbol{\omega}^{(i)}$ to increase its norm to $\sqrt{d_{0}}$. This operation only decreases the small-ball probability.

7.3 Methodology 3: Comparison with a Gaussian matrix

The first two methodologies work well for some types of test matrices, but they are unable to capture the behavior of sparse test matrices. To treat sparse random matrix models, we rely on the Gaussian comparison method, recently devised by Tropp [115]. This technique allows us to obtain a lower bound for the minimum eigenvalue of a random psd matrix by comparison with a suitable Gaussian matrix.

To state Tropp’s result, we introduce several functions that capture the variability of a random matrix. For a random self-adjoint matrix $\bm{X}\in\mathbb{F}^{d\times d}$, we define the variance and the second moment in the direction of the self-adjoint matrix $\bm{M}\in\mathbb{F}^{d\times d}$:

[TABLE]

The weak variance is defined as $\sigma_{*}^{2}(\bm{X})\coloneqq\max_{\ \mathinner{\!\left\lVert\bm{u}\right\rVert}_{2}=1}\ \operatorname*{Var}[\bm{u}^{\smash{*}}\bm{X}\bm{u}]$.

Imported Theorem 7.5** (Gaussian comparison for random psd matrices, [115, Thm. 2.3]).**

Consider a random psd matrix $\bm{W}\in\mathbb{F}^{d\times d}$ with two finite moments: $\operatorname*{\mathbb{E}}{}\mathinner{\!\left\lVert\bm{W}\right\rVert}_{\rm F}^{2}<+\infty$. Construct a self-adjoint Gaussian matrix $\bm{X}\in\mathbb{F}^{d\times d}$ that satisfies

[TABLE]

Define the random matrices

[TABLE]

Then, with probability at least $1-\delta$, the minimum eigenvalues admit the comparison

[TABLE]

Imported Theorem˜7.5 requires more effort than the other methodologies, but it has a wide scope of application. (For instance, Tropp invokes the result to establish Imported Theorem˜7.1.) In our setting, the Gaussian comparison method can be used to verify the $(r,\alpha)$-OSI property for a test matrix $\bm{\Omega}$ whose outer product decomposes as a sum of iid terms:

[TABLE]

As an example of this strategy, Section˜7.3.2 provides a detailed analysis of the SparseCol test matrix. Appendix˜B presents the analysis of the SparseStack test matrix.

7.3.1 Additional background

To employ Imported Theorem˜7.5, we collect a few additional facts from Tropp’s paper [115]. Let $\bm{W},\bm{X}\in\mathbb{F}^{d\times d}$ be random self-adjoint matrices, and let $\bm{Q}\in\mathbb{F}^{d\times r}$ be a fixed matrix. Then

[TABLE]

The inequalities between the variance function and the second-moment function hold pointwise for each self-adjoint argument $\bm{M}\in\mathbb{F}^{d\times d}$. In words, first and second moments are equivariant under conjugation [115, Prop. 3.1].

Next, we summarize the statistics of several common Gaussian matrix models. Table˜6 lists the mean, variance, expected maximum eigenvalue, and weak variance for each instance. The first example is a scalar Gaussian matrix, defined as $g\mathbf{I}_{d}$ where $g\sim{\mathcal{N}}(0,1)$. The second is an iid diagonal Gaussian matrix, namely the diagonal matrix $\bm{D}_{d}\in\mathbb{F}^{d\times d}$ with iid real diagonal entries drawn from ${\mathcal{N}}(0,1)$. The third is the compressed diagonal Gaussian matrix $\bm{Q}^{\smash{*}}\bm{D}_{d}\bm{Q}$, where $\bm{Q}\in\mathbb{F}^{d\times r}$ is an orthonormal matrix that depends on the context. Last, we consider a Gaussian unitary ensemble (GUE) matrix:

[TABLE]

The GUE matrix is rotationally invariant.

Imported Lemma 7.6** (Rotational invariance, [115, Sec. 3.7.4]).**

Suppose that $\bm{Q}\in{\mathbb C}^{d\times r}$ is orthonormal. Then the matrix $\bm{Q}^{\smash{*}}\bm{G}_{\rm gue}^{(d)}\bm{Q}$ follows the same distribution as $\bm{G}_{\rm gue}^{(r)}$.

7.3.2 Application: SparseCol embeddings

This section deploys the Gaussian comparison method to establish Theorem˜4.4, which studies the injectivity properties of the SparseCol test matrix (Definition˜4.1). The argument depends on a moment calculation, encapsulated in the next lemma. Throughout this section, we work in the complex field, as the real field is a special case.

Lemma 7.7** (SparseCol: Moment computations).**

Suppose that $d\geq 2$. Let $\boldsymbol{\omega}\in{\mathbb C}^{d}$ be the sparse random vector with $\xi$ nonzero entries, appearing in Definition˜4.1 of the SparseCol matrix. Define the random matrix $\bm{W}=k^{-1}\boldsymbol{\omega}\boldsymbol{\omega}^{\smash{*}}\in{\mathbb C}^{d\times d}$. The first moment satisfies  $\operatorname*{\mathbb{E}}{}[\bm{W}]=k^{-1}\mathbf{I}$. The second moment in the direction of the self-adjoint matrix $\bm{M}\in{\mathbb C}^{d\times d}$ satisfies

[TABLE]

The proof of Lemma 7.7 is similar with moment computations from [115], so we postpone it to Appendix˜A. Let us continue with the proof of Theorem˜4.4.

Proof of Theorem˜4.4.

To apply Imported Theorem˜7.5, we must develop a Gaussian model that compares with the SparseCol random matrix. Define the random matrix $\bm{W}$ as in Lemma˜7.7. Consider a Gaussian random matrix $\bm{X}\in{\mathbb C}^{d\times d}$ of the form

[TABLE]

As in Table˜6, the independent Gaussian random matrices $\bm{D}_{d}$ and $g\mathbf{I}_{d}$ and $\bm{G}_{\rm gue}^{(d)}$ are diagonal, scalar, and GUE. Using this dictionary, we quickly determine that $\operatorname*{\mathbb{E}}[\bm{X}]=k^{-1}\mathbf{I}=\operatorname*{\mathbb{E}}[\bm{W}]$ and that

[TABLE]

The inequality follows from Lemma˜7.7. In other words, the moments of $\bm{X}$ compare with the moments of $\bm{W}$. Now, fix any matrix $\bm{Q}\in{\mathbb F}^{d\times r}$ with orthonormal columns. Since moments are equivariant under conjugation ˜7.6, the random matrices $\bm{Q}^{{\smash{*}}}\bm{X}\bm{Q}$ and $\bm{Q}^{{\smash{*}}}\bm{W}\bm{Q}$ admit the same moment comparison.

Introduce a SparseCol random matrix $\bm{S}\in\mathbb{F}^{d\times k}$, as in the statement of Theorem˜4.4. Observe that

[TABLE]

We wish to control the quantity $\lambda_{\rm min}(\bm{Q}^{\smash{*}}\bm{S}\bm{S}^{\smash{*}}\bm{Q})=\lambda_{\rm min}(\sum_{i=1}^{k}\bm{Q}^{\smash{*}}\bm{W}_{i}\bm{Q})$. To that end, select

[TABLE]

Both $\bm{Y}$ and $\bm{Z}$ are the sums of iid matrices in $\mathbb{F}^{r\times r}$ with the appropriate moment relationships, so we can invoke Imported Theorem˜7.5. With probability at least $1-\delta$,

[TABLE]

It remains to bound the expected minimum eigenvalue below and to bound the weak variance above.

To achieve this goal, we must examine the distribution of $\bm{Z}$. Owing to the definition ˜7.7 of $\bm{X}$,

[TABLE]

This statement follows from the stability of Gaussian distributions. Imported Lemma˜7.6 states that $\bm{Q}^{\smash{*}}\bm{G}_{\rm gue}^{(d)}\bm{Q}\sim\bm{G}_{\rm gue}^{(r)}$ because $\bm{Q}\in\mathbb{F}^{d\times r}$ is an orthonormal matrix. Thus,

[TABLE]

According to Table˜6, the expected minimum eigenvalue of $\bm{Z}$ admits the lower bound

[TABLE]

The coherence $\mu(\bm{Q})$ is defined in Section˜1.8. An upper bound for the weak variance $\sigma_{*}^{2}(\bm{Z})$ follows from the subadditivity property: $\sigma_{*}^{2}(\bm{A}+\bm{B})\leq\sigma_{*}^{2}(\bm{A})+\sigma_{*}^{2}(\bm{B})$ for statistically independent $\bm{A}$ and $\bm{B}$ [115, Sec. 3.6]. Using the calculations from Table˜6,

[TABLE]

Combine these results to reach a bound for the minimum eigenvalue of $\bm{Y}$. With probability at least $1-\delta$,

[TABLE]

In particular, we can choose the column sparsity parameter \xi={\mathcal{O}}\bigl{(}\frac{k}{d}\mu(\bm{Q})\log(r/\delta)\bigr{)} and the embedding dimension $k={\mathcal{O}}(\max\{r,\log(2r/\delta)\})$. We arrive at the bound

[TABLE]

with probability at least $1-\delta$. ∎

7.4 Methodology 4: Matrix concentration

We would be remiss not to mention another established technique for obtaining injectivity properties of random matrices, based on matrix concentration tools.

Several of the standard matrix concentration inequalities from Tropp’s monograph [113] yield lower bounds on the minimum eigenvalue of a self-adjoint random matrix, expressed as an independent sum. In particular, the lower matrix Chernoff inequality [113, Thm. 5.1.1] provides a good bound on the minimum eigenvalue of a psd random matrix under minimal conditions. For example, Tropp used this strategy to analyze SRTT test matrices in his paper [111]. Unfortunately, classic matrix concentration tools are usually not powerful enough to achieve the $(r,\alpha)$-OSI property unless the embedding dimension $k=\Omega(r\log r)$. Since we aim to reach the optimal scaling $k={\mathcal{O}}(r)$, we must employ alternative methods.

Over the last several years, van Handel and coauthors have developed more advanced matrix concentration tools, including universality methods [22] and intrinsic freeness results [14]. These techniques can also be used to prove injectivity results. For example, Chenakkod et al. [33, 32] employed the universality method to study the SparseStack test matrix. They achieve an $(r,1-\varepsilon,1+\varepsilon)$-OSE with optimal dependence on $\varepsilon$, with an optimal embedding dimension $k={\mathcal{O}}(r)$, but with a suboptimal row sparsity $\zeta={\mathcal{O}}(\log^{3}r)$. As compared with the analysis via Gaussian comparison (Imported Theorem˜1.8), this approach results in excess logarithmic factors. It is likely that contemporary matrix concentration techniques can be used to establish injectivity properties for other types of random test matrices.

Appendix A Analysis of the SparseCol test matrix

This section establishes Lemma˜7.7, which controls the first two moments of the SparseCol test matrix.

Proof of Lemma˜7.7.

Recall that $\bm{W}=k^{-1}\boldsymbol{\omega}\boldsymbol{\omega}^{\smash{*}}$, where $\boldsymbol{\omega}\in\mathbb{F}^{d}$ is the sparse random vector introduced in the construction of the SparseCol random matrix (Definition˜4.1). We can easily obtain the first moment of the matrix $\bm{W}$ from the linearity of expectation, so we focus on the second-moment bound. The quantity of interest can be expressed as

[TABLE]

Here, $m_{ij}$ and $\omega_{i}$ refer to the entries of $\bm{M}$ and $\boldsymbol{\omega}$. In view of Definition˜4.1, the entries of the sparse random vector take the form $\omega_{i}=\sqrt{d/\xi}\sum_{q=1}^{\xi}\varrho_{q}\mathbbm{1}_{\{s_{q}=i\}}$. Therefore, the product of four such entries expands as a sum:

[TABLE]

Each random variable $\varrho_{q}\sim\textsc{rademacher}$, so

[TABLE]

Therefore, in ˜A.2, it suffices to examine the sum over tuples $(q_{1},q_{2},q_{3},q_{4})$ where the indicators are nonzero. We proceed through each of the four expectations. In the first case, we have $q_{1}=q_{2}=q_{3}=q_{4}$, whence

[TABLE]

In the second case, we have $q_{1}=q_{2}\neq q_{3}=q_{4}$, whence

[TABLE]

The third and fourth cases are symmetric to the second, but with the $i,j,u,t$ symbols permuted appropriately. Substitute the last three displays into ˜A.2 and then into ˜A.1, and simplify the result. ∎

Appendix B Analysis of the SparseStack test matrix

In this section, we prove Imported Theorem˜1.8, which establishes an OSI guarantee for the SparseStack test matrix (Definition˜1.7). We work in the complex field, as the real field is just a special case. The SparseStack random matrix $\bm{\Omega}\in{\mathbb C}^{d\times k}$ has a sparsity parameter $\zeta\in{\mathbb N}$ and a block size parameter $b\in{\mathbb N}$. The embedding dimension $k\coloneqq b\zeta$.

We can express the SparseStack matrix in terms of a CountSketch matrix [27, 35, 78, 90]:

[TABLE]

As usual, $\mathbf{e}_{i}\in{\mathbb C}^{b}$ is a standard basis vector. Then the SparseStack matrix takes the form

[TABLE]

This characterization is valuable because it allow us to write $\bm{\Omega}\bm{\Omega}^{{\smash{*}}}=\zeta^{-1}\sum_{i=1}^{\zeta}\bm{\Phi}_{i}\bm{\Phi}_{i}^{\smash{*}}$, which is an average of iid random psd matrices.

To prove Imported Theorem˜1.8, we employ the Gaussian comparison approach (Imported Theorem˜7.5). As with the SparseCol matrix, the first step requires us to compute first and second moments of our test matrix.

Lemma B.1** (CountSketch moments).**

Let $\bm{\Phi}\in{\mathbb C}^{d\times b}$ be the CountSketch matrix ˜B.1. Then  $\operatorname*{\mathbb{E}}[\bm{\Phi}\bm{\Phi}^{\smash{*}}]=\mathbf{I}$. In addition, for each self-adjoint matrix $\bm{M}\in{\mathbb C}^{d\times d}$,

[TABLE]

We defer proof of Lemma˜B.1 to Section˜B.1. With the result at hand, we turn to the proof of Imported Theorem˜1.8.

Proof of Imported Theorem˜1.8.

Let $\bm{\Omega}\in{\mathbb C}^{k\times d}$ be the SparseStack test matrix constructed in ˜B.2. It is easy to verify that $\bm{\Omega}$ is isotropic ˜1.2, so we focus on establishing the injectivity property ˜1.3.

To begin, write the outer product $\bm{\Omega}\bm{\Omega}^{\smash{*}}$ as an iid sum. Define $\bm{W}\coloneqq\zeta^{-1}\bm{\Phi}\bm{\Phi}^{\smash{*}}\in{\mathbb C}^{d\times d}$, and note that

[TABLE]

Lemma˜B.1 implies that the first moment $\operatorname*{\mathbb{E}}[\bm{W}]=\zeta^{-1}\mathbf{I}$ and that the second moment

[TABLE]

Now, we build a Gaussian comparison model. Let $\bm{X}\in{\mathbb C}^{d\times d}$ be the Gaussian matrix

[TABLE]

where $g\mathbf{I}_{d}$ and $\bm{G}_{\rm gue}^{(d)}$ are the scaled identity and GUE matrices defined in Table˜6. We can compare the moments of $\bm{X}$ and $\bm{W}$. The first moments are equal: $\operatorname*{\mathbb{E}}[\bm{X}]=\frac{1}{\zeta}\mathbf{I}_{d}=\operatorname*{\mathbb{E}}[\bm{W}]$. The second moments satisfy

[TABLE]

In other words, $\bm{X}$ is a Gaussian comparison model for $\bm{W}$. For any orthonormal matrix $\bm{Q}\in{\mathbb C}^{d\times r}$, we can use equivariance ˜7.6 to confirm that $\bm{Q}^{\smash{*}}\bm{X}\bm{Q}$ is a Gaussian comparison model for $\bm{Q}^{\smash{*}}\bm{W}\bm{Q}$.

We seek to establish control on the minimum eigenvalue $\lambda_{\rm min}(\bm{Q}^{\smash{*}}\bm{\Omega}\bm{\Omega}^{\smash{*}}\bm{Q})=\lambda_{\rm min}(\sum_{i=1}^{\zeta}\bm{Q}^{\smash{*}}\bm{W}_{i}\bm{Q})$. To do so, we introduce iid copies $\bm{X}_{1},\ldots,\bm{X}_{\zeta}$ of the Gaussian matrix $\bm{X}$, and we define

[TABLE]

Since $\bm{Y}$ and $\bm{Z}$ are sums of iid matrices in ${\mathbb C}^{d\times d}$, Imported Theorem˜7.5 implies that

[TABLE]

To prove the theorem, we must establish a lower bound on the expected minimum eigenvalue $\operatorname*{\mathbb{E}}[\lambda_{\rm min}(\bm{Z})]$ and an upper bound on the weak variance $\sigma_{*}^{2}(\bm{Z})$. To do so, we first examine the distribution of $\bm{Z}$ more closely. By linearity, Gaussianity, and the definition $k=b\zeta$,

[TABLE]

Now, we conjugate by $\bm{Q}\in\mathbb{F}^{d\times r}$. By rotational invariance (Imported Lemma˜7.6), $\bm{Q}^{\smash{*}}\bm{G}_{\rm gue}^{(d)}\bm{Q}\sim\bm{G}_{\rm gue}^{(r)}$. Thus,

[TABLE]

By symmetry of the Gaussian distribution with respect to negation, we can bound the expected minimum eigenvalue of $\bm{Z}$ below:

[TABLE]

We have introduced results from Table˜6. Next, we bound the weak variance above via the subadditivity property [115, Sec. 3.6] and the results from Table˜6:

[TABLE]

where the last inequality relies on the fact $\zeta\leq k$. We conclude that, with probability at least $1-\delta$,

[TABLE]

Set $\delta=\nicefrac{{1}}{{20}}$, and note that the right-hand side exceeds $\nicefrac{{1}}{{2}}$ at an embedding dimension $k={\mathcal{O}}(r)$ and a sparsity level $\zeta={\mathcal{O}}(\log r)$. This argument confirms the injectivity property ˜1.3. ∎

B.1 Moments of the CountSketch matrix

In this section, we calculate the first two moments of the CountSketch matrix.

Proof of Lemma˜B.1.

Introduce the CountSketch matrix, $\bm{\Phi}\in{\mathbb C}^{d\times b}$ described in ˜B.1, whose entries are denoted $\phi_{i\alpha}$. The first moment $\operatorname*{\mathbb{E}}{}[\bm{\Phi}\bm{\Phi}^{\smash{*}}]=\mathbf{I}$ by linearity of expectation, so we focus on the second moment. To begin, we use the definition of the second moment function to reach the expansion

[TABLE]

In the last line, we defined the quantities $\theta_{ijtu}\;{\vcentcolon=}\;\sum_{\alpha,\beta=1}^{b}\operatorname*{\mathbb{E}}[\phi_{i\alpha}\phi_{j\alpha}\phi_{t\beta}\phi_{u\beta}]$. We will compute each of these terms before returning to ˜B.3.

By construction of the CountSketch matrix, its entries take the form $\phi_{i\alpha}=\varrho_{i}\mathbbm{1}_{\{s_{i}=\alpha\}}$. Thus,

[TABLE]

Each random variable $\varrho_{q}\sim\textsc{rademacher}$, so

[TABLE]

This equation follows from the fact that the rows of $\bm{\Phi}$ are independent. Therefore, in ˜B.4, it suffices to compute $\theta_{ijtu}$ for tuples $(i,j,t,u)$ that make the indicators nonzero. We proceed through each of the four indicators. In the first case, we have $i=j=t=u$. Then

[TABLE]

Thus, $\theta_{iiii}=\sum_{\alpha,\beta=1}^{b}b^{-1}\mathbbm{1}_{\{\alpha=\beta\}}=1$. In the second case, $i=j\neq t=u$, so

[TABLE]

Thus, $\theta_{iitt}=\sum_{\alpha,\beta=1}^{b}b^{-2}=1$. In the third case, $i=t\neq j=u$, so

[TABLE]

Thus, $\theta_{ijij}=\sum_{\alpha,\beta=1}^{b}b^{-2}\mathbbm{1}_{\{\alpha=\beta\}}=b^{-1}$. In the fourth case, $i=u\neq j=t$. This case is the same as the third, and we determine that $\theta_{ijji}=b^{-1}$. Substitute these identities for $\theta_{ijtu}$ back into ˜B.3 to reach

[TABLE]

In the last step, since $\bm{M}$ is self-adjoint, we may bound $\sum_{i\neq j}m_{ij}^{2}=\sum_{i\neq j}\mathrm{Re}(m_{ij}^{2})\leq\mathinner{\!\left\lVert\bm{M}\right\rVert}_{\rm F}^{2}$. ∎

Appendix C Proof of Theorem˜7.3

To prove Theorem˜7.3, we use the following result from VC theory:

Theorem C.1** (Symmetric classifiers).**

Let ${\mathcal{F}}$ be the class of symmetric two-sided linear classifiers

[TABLE]

and let $\bm{x}_{1},\ldots,\bm{x}_{k}$ be iid copies of a random vector $\bm{x}\in{\mathbb R}^{r}$. With probability at least $1-\delta$, all $f\in{\mathcal{F}}$ have

[TABLE]

Asymptotically equivalent results are known [114, 82]; we provide sharper constants. Our version of Theorem˜C.1 is proven in Section˜C.1. This result simplifies the results in [112] and obtains tighter constants. The proof is by the small-ball method [56, 114, 120]. With this preparation in place, we frame a bound on the minimum singular value of a matrix with iid columns:

Theorem C.2** (Minimum singular value by small-ball method, real case).**

Let $\bm{X}=[\bm{x}_{1}~\cdots~\bm{x}_{k}]\in{\mathbb R}^{r\times k}$, where $\bm{x}_{i}\in{\mathbb R}^{r}$ are iid copies of random vector $\bm{x}$. Then, for all $\tau>0$ and $\delta\in(0,1)$, we have that

[TABLE]

holds with probability at least $1-\delta$.

Proof.

We expand the minimum singular value of $\bm{X}$:

[TABLE]

Then, since $t\geq\tau\mathbbm{1}_{\{t\geq\tau\}}$ for positive $t$, we have that

[TABLE]

Notice that the term inside the minimum is an empirical mean of indicator variables. Thus, the stated result follows immediately by Theorem˜C.1. ∎

Notice that Theorem˜C.2 is locked into the real field ${\mathbb F}={\mathbb R}$. The real case of Theorem˜7.3 follows immediately from this result and a few short numerical calculations. We discuss the complex case below in Section˜C.2. Of potential independent interest, we also provide a high-probability version of this result with small explicit constants:

Theorem C.3** (High probability subspace injection by small-ball method).**

Let $\bm{\Omega}=\frac{1}{\sqrt{k}}[\boldsymbol{\omega}_{1}~\cdots~\boldsymbol{\omega}_{k}]\in{\mathbb F}^{d\times k}$ be generated from iid copies $\boldsymbol{\omega}_{i}$ of an isotropic random vector $\boldsymbol{\omega}$. Suppose that \mathbb{P}\bigl{\{}|\langle\bm{u},\boldsymbol{\omega}\rangle|^{2}\geq\tau\bigr{\}}\geq 0.99 for all unit vectors $\bm{u}\in{\mathbb F}^{d}$ and some value $\tau>0$. Then for all orthonormal matrices $\bm{Q}\in{\mathbb F}^{d\times r}$, it holds that $\sigma_{\rm min}(\bm{Q}^{\smash{*}}\bm{\Omega})\geq\frac{\tau}{10}$ with probability at least $1-\delta$ provided $k\geq\max\{4r,10\log(1/\delta)\}$ if ${\mathbb F}={\mathbb R}$ or $k\geq\max\{18r,57\log(1/\delta)\}$ if $\>{\mathbb F}={\mathbb C}$.

C.1 Refined VC dimension bounds

In this section, we develop Theorem˜C.1 using standard tools from VC theory carefully. We first prove the real case, and extend our results to the complex case in Section˜C.2. We start by importing useful definitions and concentration results. We start with the growth function:

Definition C.4** (Growth function, [82, Def. 3.6]).**

Let ${\mathcal{X}}$ be a set, and let ${\mathcal{F}}$ be a set of functions mapping from ${\mathcal{X}}$ to $\{0,1\}$. Then, the growth function of ${\mathcal{F}}$ is defined as

[TABLE]

The growth function $\Pi_{{\mathcal{F}}}(k)$ counts the number of ways functions in ${\mathcal{F}}$ can assign labels to $k$ points in ${\mathcal{X}}$.

We list two standard properties of the growth function:

Imported Theorem C.5** (Concentration of empirical mean, adapted from [82, Cor. 3.9]).**

Let ${\mathcal{X}}$ be a set, and let ${\mathcal{F}}$ be a set of functions that map from ${\mathcal{X}}$ to $\{0,1\}$. Let $\bm{x}_{1},\ldots,\bm{x}_{i}$ be iid copies of random variable $\bm{x}\in{\mathcal{X}}$. Then, with probability at least $1-\delta$, it holds simultaneously for all $f\in{\mathcal{F}}$ that

[TABLE]

Imported Lemma C.6** (Self-intersection/union, [82, Exercise 3.23]).**

Let ${\mathcal{F}}$ be a set of functions that maps from ${\mathcal{X}}$ to $\{0,1\}$, and define the self-intersection of ${\mathcal{F}}$ as ${\mathcal{F}}_{\cap}\;{\vcentcolon=}\;\left\{\bm{x}\mapsto f_{1}(\bm{x})f_{2}(\bm{x})~\mathrel{\mathop{\ordinarycolon}}~f_{1},f_{2}\in{\mathcal{F}}\right\}.$ Then, it holds that $\Pi_{{\mathcal{F}}_{\cap}}\leq\Pi_{{\mathcal{F}}}^{2}.$ The same bound holds for the self-union ${\mathcal{F}}_{\cup}\;{\vcentcolon=}\;\left\{\bm{x}\mapsto\max\{f_{1}(\bm{x}),f_{2}(\bm{x})\}~\mathrel{\mathop{\ordinarycolon}}~f_{1},f_{2}\in{\mathcal{F}}\right\}$.

The product and maximum in the above result represent the boolean “and” and “or” operations on $\{0,1\}$. The bound for the self-union follows from the bound for the self-intersection by De Morgan’s law. Next, we import the definition of the VC dimension and recall its relation to the growth function:

Definition C.7** (VC dimension, [82, Def. 3.10]).**

Let ${\mathcal{X}}$ be a set, and let ${\mathcal{F}}$ be a set of functions that map from ${\mathcal{X}}$ to $\{0,1\}$. The VC dimension of ${\mathcal{F}}$ is $\text{VCdim}({\mathcal{F}})\;{\vcentcolon=}\;\max\{k~\mathrel{\mathop{\ordinarycolon}}~\Pi_{{\mathcal{F}}}(k)=2^{k}\}.$

Imported Lemma C.8** (VC dimensions of linear classifiers, [106, Thm. 9.2]).**

The VC dimension of the set of linear classifiers ${\mathcal{F}}\;{\vcentcolon=}\;\left\{\bm{x}\mapsto\mathbbm{1}_{\{\langle\bm{x},\bm{u}\rangle>\tau\}}~\mathrel{\mathop{\ordinarycolon}}~\bm{u}\in{\mathbb R}^{r},\tau\in{\mathbb R}\right\}$ is $r$.

Imported Lemma C.9** (Sauer’s Lemma, [82, Cor. 3.18]).**

Let ${\mathcal{F}}$ be a set of function with VC dimension $r$. Then for $k\geq r$, the growth function is bounded $\Pi_{{\mathcal{F}}}(k)\leq(\mathrm{e}k/r)^{r}.$

Proof of Theorem˜C.1.

Let ${\mathcal{H}}$ be the set of linear classifiers in Imported Lemma˜C.8. Notice that for any two-sided linear classifier $f\in{\mathcal{F}}$, we can write

[TABLE]

We conclude ${\mathcal{F}}\subseteq{\mathcal{H}}_{\cap}$. Thus, by Imported Lemmas˜C.6, C.8, and C.9, we find that

[TABLE]

Therefore, by Imported Theorem˜C.5, we know that simultaneously for all $f\in{\mathcal{F}}\subseteq{\mathcal{H}}_{\cap}$,

[TABLE]

Taking a supremum of bound sides yields the stated result. ∎

C.2 Analysis in the complex field

To adapt our results to the complex setting, we first develop a version of Theorem˜C.1 for unions of symmetric two-sided linear classifiers and a complex-valued version of Theorem˜C.2.

Theorem C.10** (Union of symmetric classifiers).**

Let ${\mathcal{F}}$ be the class of unions of symmetric classifiers

[TABLE]

and let $\bm{x}_{1},\ldots,\bm{x}_{k}$ be iid copies of a random vector $\bm{x}\in{\mathbb R}^{r}$. Then, with probability at least $1-\delta$, we have that

[TABLE]

Proof.

Notice that ${\mathcal{F}}$ is exactly the self-union of the class of two-sided symmetric linear classifiers as defined in Theorem˜C.1. The proof of the advertise result then follows the proof of Theorem˜C.1 but instead using Imported Lemma˜C.6 to bound $\Pi_{\mathcal{F}}(k)\leq(\mathrm{e}k/r)^{4r}$. ∎

Theorem C.11** (Minimum singular value by small-ball method, complex case).**

Let $\bm{X}=[\bm{x}_{1}~\cdots~\bm{x}_{k}]\in{\mathbb C}^{r\times k}$, where $\bm{x}_{i}\in{\mathbb C}^{r}$ are iid copies of random vector $\bm{x}$. Then, for all $\tau>0$ and $\delta\in(0,1)$, we have that

[TABLE]

holds with probability at least $1-\delta$.

Proof.

As in the proof of Theorem˜C.2, we lower bound the minimum singular value of $\bm{X}$ as

[TABLE]

Note that $|\langle\bm{u},\bm{x}_{i}\rangle|^{2}\geq\tau$ only if $|\mathrm{Re}\>\langle\bm{u},\bm{x}_{i}\rangle|^{2}\geq\tau/2$ or $|\mathrm{Im}\>\langle\bm{u},\bm{x}_{i}\rangle|^{2}\geq\tau/2$. Now let $\varphi\mathrel{\mathop{\ordinarycolon}}{\mathbb C}^{r}\rightarrow{\mathbb R}^{2r}$ be the map $\varphi(\bm{x})\;{\vcentcolon=}\;\left[\begin{smallmatrix}\mathrm{Re}\bm{x}\\ \mathrm{Im}\bm{x}\end{smallmatrix}\right]$. We have that $\mathrm{Re}\langle\bm{u},\bm{x}_{i}\rangle=\langle\varphi(\bm{u}),\varphi(\bm{x}_{i})\rangle$ and $\mathrm{Im}\langle\bm{u},\bm{x}_{i}\rangle=\langle\bm{P}\varphi(\bm{u}),\varphi(\bm{x}_{i})\rangle$ for permutation matrix $\bm{P}=\left[\begin{smallmatrix}\bm{0}&\mathbf{I}\\ \mathbf{I}&\bm{0}\end{smallmatrix}\right]$. By substituting in this map, we get

[TABLE]

Notice that the term inside the minimum is an empirical mean of indicator variables involving vectors in ${\mathbb R}^{2r}$. Thus, the stated result follows immediately by Theorem˜C.10. ∎

The complex case of Theorem˜7.3 now follows as a corollary of Theorem˜C.11.

Appendix D Implementation of scientific data compression experiment

In this appendix, we describe the experiment from Section˜6.1 in greater detail. The dynamics of a Bose–Einstein condensate trapped in an external potential $V_{\mathrm{trap}}$ rotating around the $z$-axis at an angular frequency $\Omega$ are governed by the Gross–Pitaevskii equation (GPE):

[TABLE]

where $\hbar$ is the reduced Planck constant, $m$ is the atomic mass, $g$ is the interparticle interaction strength, and L_{z}=-i\hbar\bigl{(}x\,\partial_{y}-y\,\partial_{x}\bigr{)} is the $z$-component of the angular-momentum operator.

To study ground states, we follow the standard procedure [125] and consider the evolution of this equation in imaginary time. Which choose the initial state $\psi_{0}(\bm{x},0)$ to be a Gaussian profile, carrying a phase defect in a harmonic potential $V_{h}$. In order to simulate the imaginary-time GPE evolution, we use a Fourier pseudospectral method in space and a semi-implicit backward/forward Euler scheme in time following [125]. The initial state $\psi_{0}(\bm{x},0)$ is evolved for $40{,}000$ time steps using a time step of $\Delta t=10^{-4}$ on a $256\times 256$ domain. The pde data $\bm{x}(t)\in{\mathbb C}^{65{,}536}$ is natively complex, but we convert to real data by concatenating the real and imaginary parts $(\mathrm{Re}\bm{x}(t),\mathrm{Im}\bm{x}(t))\in{\mathbb R}^{131{,}072}$.

Appendix E Implementation of partition function estimation experiment

In this appendix, we describe the experimental setup for Section˜6.2 in greater detail. Partition function estimation for the transverse-field Ising model was studied in [49]. Following their lead, we shift the Hamiltonian matrix by a constant $b=(1+h)\ell$ before applying stochastic trace estimation so that $\bm{H}+b\mathbf{I}$ is psd, allowing us to use variance-reduced trace estimators for psd matrices $\bm{A}$ like Nyström++ [95] and XNysTrace [49]. To compute $\exp(-\beta\bm{H})\bm{\Omega}$ we follow [6] evaluating the matrix exponential matvec implicitly via a truncated Taylor approximation. We observe that when paired with Khatri–Rao test matrices, the previously mentioned variance-reduced trace estimators strongly outperform the naive Girard–Hutchinson trace estimator for this problem.

Acknowledgements

CC was supported by the Caltech Kortschak Scholars Program and the NSF GRFP Award 2139433. RAM, ENE, and JAT were supported by ONR Award N00014-24-1-2223, a Caltech Center for Sensing to Intelligence grant, and the Caltech Carver Mead New Adventures Fund. ENE was also supported by the DOE CSGF under Award DE-SC0021110.