# Sketching with Kerdock’s crayons:
Fast sparsifying transforms for arbitrary linear maps

Tim Fuchs111Department of Mathematics, Technical University of Munich,
München, Germany David Gross222Institute for Theoretical Physics, University
of Cologne, Cologne, Germany Felix Krahmer11footnotemark: 1
Richard Kueng333Institute for Integrated Circuits, Johannes Kepler University
Linz, Linz, Austria Dustin G. Mixon444Department of Mathematics, The Ohio
State University, Columbus, Ohio, USA 555Translational Data Analytics
Institute, The Ohio State University, Columbus, Ohio, USA

###### Abstract

Given an arbitrary matrix $A\in\mathbb{R}^{n\times n}$, we consider the
fundamental problem of computing $Ax$ for any $x\in\mathbb{R}^{n}$ such that
$Ax$ is $s$-sparse. While fast algorithms exist for particular choices of $A$,
such as the discrete Fourier transform, there is currently no $o(n^{2})$
algorithm that treats the unstructured case. In this paper, we devise a
randomized approach to tackle the unstructured case. Our method relies on a
representation of $A$ in terms of certain real-valued mutually unbiased bases
derived from Kerdock sets. In the preprocessing phase of our algorithm, we
compute this representation of $A$ in $O(n^{3}\log n)$ operations. Next, given
any unit vector $x\in\mathbb{R}^{n}$ such that $Ax$ is $s$-sparse, our
randomized fast transform uses this representation of $A$ to compute the
entrywise $\epsilon$-hard threshold of $Ax$ with high probability in only
$O(sn+\epsilon^{-2}\|A\|_{2\to\infty}^{2}n\log n)$ operations. In addition to
a performance guarantee, we provide numerical results that demonstrate the
plausibility of real-world implementation of our algorithm.

##  1 Introduction

Computing matrix–vector products is a fundamental part of numerical linear
algebra. The naive algorithm takes $O(mn)$ operations to multiply an $m\times
n$ matrix by a vector. Many structured matrices admit a more efficient
implementation of this computation, the most well-known example being the fast
Fourier transform, which takes only $O(n\log n)$ operations. In some
applications, the desired Fourier transform of the given signal is nearly
$s$-sparse, and as we discuss below, a number of works have proposed methods
for such cases that are sublinear in the dimension $n$.

For the one-dimensional discrete Fourier transform, a randomized algorithm
with a runtime scaling quadratically in $s$ up to logarithmic factors in the
dimension $n$ has been provided in [15], while a deterministic approach with
similar complexity was found in [21, 22]. In later works, this could be
reduced to linear scaling in $s$ for both random [16, 21, 18, 19, 22] and
deterministic [27, 9] algorithms. For the $d$-dimensional Fourier transform
applied to signals in $n=N^{d}$ dimensions, the exponential scaling in $d$
presents an instance of the curse of dimensionality. Despite this, for random
signals, one may obtain runtimes that are linear in $sd$ up to logarithmic
factors in $N$ [5, 6]. For deterministic signals, various deterministic [32,
22] and random [20, 24, 7, 8, 23] sampling strategies have been proposed with
a computational complexity which scales polynomially in $d$, $s$ and $N$ up to
logarithmic factors.

Naturally, research on fast transforms is not restricted to Fourier structure.
For example, [35] proposes a multidimensional Chebyshev transform with reduced
runtime. In [7], a more general approach has been established that yields fast
sparse transforms for arbitrary bounded orthonormal product basis with a
runtime scaling polynomially in $s$ up to logarithmic factors. These results
have been generalized in [8] to signals with only an approximately $s$-sparse
representation while maintaining a computational complexity that is sublinear
in the dimension $n$. While covering a significantly larger class of
transforms than just the Fourier transform, all these approaches remain
restricted to a specific structure or class of structures of the
transformation matrix.

At the same time, data-driven sparsifying transforms, which have been
demonstrated to outperform predefined structured representation systems in a
variety of contexts [13, 2, 40, 31, 36, 34], typically do not have structural
properties that allow for the application of any of the above fast transform
methods. This issue was addressed in [37, 38] by imposing structure amenable
to fast transforms on the learned representation system $A$ so as to
facilitate the computation of $Ax$. At the same time, this imposed structure
significantly limits the space of admissible transforms, and the question
remains whether a fast transform can also be constructed for learned
representation systems beyond these restrictions.

In this paper, we consider cases where the desired product $Ax$ is
approximately sparse for a matrix $A$ that does not follow any preset
structural constraints, e.g., because it is learned from data. In particular,
we assume $A$ is arbitrary. Note that to compute the mapping $(A,x)\mapsto Ax$
for an arbitrary matrix $A$ and vector $x$, one must first read the input
$(A,x)\in\mathbb{R}^{m\times n}\times\mathbb{R}^{n}$. Since this already
requires $\Theta(mn)$ operations, naive matrix–vector multiplication is
optimally efficient when computing an individual matrix–vector product. To
obtain a speedup, we instead apply the same transform $x\mapsto Ax$ to a
stream of vectors $x\in\mathbb{R}^{n}$, which models the setting of many
applications. Our approach will require some upfront preprocessing given $A$
in exchange for a much faster per-vector computation. While little work has
been done in this vein, there has been quite a bit of work on related
problems, which we discuss below.

The first result following this strategy [42] concerns a matrix–vector
multiplication algorithm over finite semirings (for general matrices and
vectors, hence not assuming any kind of sparsity), which performs
$O(n^{2+\epsilon})$ operations of preprocessing on an $n\times n$ matrix
before multiplying with an arbitrary vector in $O(n^{2}/(\epsilon\log n)^{2})$
operations. The first and (to our knowledge) only algorithm that achieves a
comparable result for real $m\times n$ matrices is the mailman algorithm
introduced in [29]. Provided the matrix contains only a constant number of
distinct values, the algorithm takes $O(mn)$ operations of preprocessing and
then takes $O(mn/\log(m+n))$ operations to multiply the preprocessed matrix
with an arbitrary vector.

As an alternative, one might batch the stream of vectors into matrices and
then perform matrix multiplication. (Granted, such a batched computation is
unacceptable for many applications.) Research on matrix multiplication was
initiated by the seminal work of Strassen [39], which multiplies two arbitrary
$n\times n$ matrices in only $O(n^{2.808})$ operations (i.e., much faster than
the naive $O(n^{3})$ algorithm). Later algorithms [10, 11, 43, 28] improved
this computational complexity to its currently best known scaling of
$O(n^{2.373})$. After dividing by the batch size, this gives a per-vector cost
of $O(n^{1.373})$ operations. However, we note that such algorithms are
infeasible in practice.

A more feasible approach to matrix multiplication was proposed by Drineas,
Kanan, and Mahoney [12]. They compute a random approximation of the desired
product by multiplying two smaller matrices: one consisting of $k$ randomly
selected columns of the first matrix $A$, and the other consisting of
corresponding rows from the second matrix $B$. With high probability, the
Frobenius norm of the estimate error is $O(\|A\|_{F}\|B\|_{F}/\sqrt{k})$.
Unfortunately, if $B$ represents a batch of column vectors, then this
guarantee offers little control of the error in each vector. On the other
hand, if $B$ represents a single column vector $b$, then $\|Ab\|_{2}$ is
typically much smaller than $\|A\|_{F}\|b\|_{2}$, so the resulting relative
error is quite large even for relatively large values of $k$.

###  1.1 Our approach

Given $A\in\mathbb{R}^{m\times n}$, let $\Sigma(A,s,\delta)$ denote the set of
all unit vectors $x\in\mathbb{R}^{n}$ for which

$\inf\Big{\\{}\|Ax-v\|_{\infty}:v\in\mathbb{R}^{m},|\operatorname{supp}(v)|\leq
s\Big{\\}}\leq\delta,$

i.e., $Ax$ is $\delta$-close to being $s$-sparse. Given $\epsilon>0$, let
$h_{\epsilon}\colon\mathbb{R}\to\mathbb{R}$ denote the $\epsilon$-hard
thresholding function defined by
$h_{\epsilon}(t):=t\cdot\mathbf{1}\\{|t|\geq\epsilon\\}$. By abuse of
notation, we apply $h_{\epsilon}$ to the entries of a vector $v$ by writing
$h_{\epsilon}(v)$. We seek to solve the following:

######  Problem 1.

Given an arbitrary $A\in\mathbb{R}^{m\times n}$, $s\in\mathbb{N}$, and
$\epsilon>\delta\geq 0$, preprocess $A$ so that one may quickly compute
$h_{\epsilon}(Ax)$ for any $x\in\Sigma(A,s,\delta)$.

Our approach uses a specially designed random vector $z\in\mathbb{R}^{n}$ such
that $\mathbb{E}zz^{\top}=I$:

$Ax=A(\mathbb{E}zz^{\top})x=\mathbb{E}Azz^{\top}x.$ (1)

Denote the random vector $y:=Azz^{\top}x\in\mathbb{R}^{m}$. The fact that
$\mathbb{E}y=Ax$ suggests a Monte Carlo approach to estimate $Ax$. That is, we
will approximate the true average $Ax$ with an estimator determined by $N$
independent samples. To obtain a fast algorithm in this vein, we will select a
distribution for $z$ and an estimator $\hat{\mu}$ for $\mathbb{E}y$ that
together satisfy three properties:

  * (i)

the distribution of $z$ is discrete with small support,

  * (ii)

$\hat{\mu}$ can be computed in linear time from independent realizations
$\\{y_{j}\\}_{j\in[N]}$ of $y$, and

  * (iii)

for each $x\in\Sigma(A,s,\delta)$,
$\|\hat{\mu}-Ax\|_{\infty}<\frac{\epsilon-\delta}{2}$ with high probability,
even for small $N$.

Indeed, if (i)–(iii) hold, then one may compute $h_{\epsilon}(Ax)$ using the
following (fast) algorithm:

Let $\\{s_{\ell}:\ell\in[L]\\}\subseteq\mathbb{R}^{n}$ denote the support of
the distribution of $z$. Given $A\in\mathbb{R}^{m\times n}$, we run the
preprocessing step of computing $\\{As_{\ell}\\}_{\ell\in[L]}$ in $O(Lmn)$
operations. Granted, this is more expensive than computing $Ax$, but we only
need to compute $\\{As_{\ell}\\}_{\ell\in[L]}$ once, while we expect to
compute $Ax$ for a stream of $x$’s. Next, given $x\in\Sigma(A,s,\delta)$, we
draw independent realizations $\\{z_{j}\\}_{j\in[N]}$ of $z$. Since we already
computed $\\{As_{\ell}\\}_{\ell\in[L]}$, we may then compute the corresponding
realizations $\\{y_{j}\\}_{j\in[N]}$ of $y$ in $O(N(m+n))$ operations. Next,
by (ii), we may compute $\hat{\mu}$ from $\\{y_{j}\\}_{j\in[N]}$ in $O(mN)$
operations. Finally, let $S\subseteq[m]$ denote the indices of the $s$ entries
of $\hat{\mu}$ of largest magnitude. Then by (iii), it holds with high
probability that $|\hat{\mu}_{i}|>\frac{\epsilon+\delta}{2}$ for every
$i\in\operatorname{supp}(h_{\epsilon}(Ax))$ while
$|\hat{\mu}_{i}|<\frac{\epsilon+\delta}{2}$ for every $i\in[m]$ such that
$|(Ax)_{i}|\leq\delta$. Since $x\in\Sigma(A,s,\delta)$ by assumption, it
follows that $S\supseteq\operatorname{supp}(h_{\epsilon}(Ax))$, and so
$A_{S}x$ determines $h_{\epsilon}(Ax)$, which we compute in $O(sn)$ additional
operations. (Of course, $A_{S}x$ might determine other entries in the support
of $Ax$, and we would not discard this information in practice.)

To obtain (i)–(iii), we take $z$ to be uniformly distributed over an
appropriately scaled $n$-dimensional projection of a projective $2$-design,
and for $\hat{\mu}$, we partition $[N]$ into batches and compute the entrywise
median of means of $\\{y_{j}\\}_{j\in[N]}$ over these batches. The projective
$2$-design allows us to control the variance of each entry of the random
vector $y$; see Lemma 4. Next, the median-of-means estimator improves over the
sample mean by being less sensitive to outliers in the small random sample
$\\{y_{j}\\}_{j\in[N]}$. Thanks to this behavior, we can get away with drawing
only $N=O((\epsilon-\delta)^{-2}\|A\|_{2\to\infty}^{2}\log(m/\eta))$ samples,
where the induced norm $\|A\|_{2\to\infty}$ equals the largest $\ell_{2}$-norm
of the rows of $A$, and $\eta$ denotes the failure probability of the
randomized algorithm; see Theorem 5. As a bonus, the Kerdock set–based
projective $2$-design we use enjoys a fast matrix–vector multiplication
algorithm, yielding a preprocessing step of only $O(mn^{2}\log
n+n^{2}\log^{2}n)$ operations despite having $L=\Theta(n^{2})$; see Lemma 8.
See Algorithm 1 for a summary of our approach.

Data: Parameters $\epsilon>\delta\geq 0$, $s,J,K\in\mathbb{N}$, matrix
$A\in\mathbb{R}^{m\times n}$, stream of $x\in\Sigma(A,s,\delta)$

Result: Entrywise $\epsilon$-hard threshold of matrix–vector products
$h_{\epsilon}(Ax)$

Preprocessing step

Let $d$ denote the smallest power of $2$ that is at least $n$, and put
$L:=d(d/2+1)$

Let $\\{u_{\ell}\\}_{\ell\in[L]}$ denote a projective $2$-design for
$\mathbb{R}^{d}$ arising from a Kerdock set

Put $s_{\ell}:=\sqrt{d}\Pi u_{\ell}$, where $\Pi\in\mathbb{R}^{n\times d}$
denotes projection onto the first $n$ coordinates

Use Lemma 8 to compute $\\{As_{\ell}\\}_{\ell\in[L]}$

Streaming step

Draw $N:=JK$ indices $\\{\ell_{j}\\}_{j\in[N]}$ uniformly from $[L]$

Compute $y_{j}:=(As_{\ell_{j}})(s_{\ell_{j}}^{\top}x)$ for each $j\in[N]$

Compute the entrywise median of means $\hat{\mu}$ of $\\{y_{j}\\}_{j\in[N]}$
over $K$ batches of size $J$

Let $S\subseteq[m]$ denote the indices of the $s$ entries of $\hat{\mu}$ with
largest magnitude

Compute $A_{S}x$ and output the indices and values of entries with magnitude
at least $\epsilon$

Algorithm 1 Fast sparsifying transform for an arbitrary linear map

To quickly evaluate the utility of this algorithm, consider the following
model: $A$ is an arbitrary $n\times n$ orthogonal matrix, and $x$ is a random
vector such that the entries of $Ax$ are drawn independently from the
following mixture:

$e_{i}^{\top}Ax\sim\left\\{\begin{array}[]{cl}\mathcal{N}(0,1)&\text{with
probability }p\\\ 0&\text{with probability }1-p.\end{array}\right.$

Then the expected size of the support of $Ax$ is $pn$. In this model, our
algorithm provides a speedup over naive matrix–vector multiplication in the
regime

$1\prec pn\prec\frac{n}{\log n}.$

To see this, first put $s:=10pn$ (say). Then the multiplicative Chernoff bound
implies that $Ax$ is $s$-sparse with high probability, and so we take
$\delta=0$. Before selecting $\epsilon>0$, we normalize our vector so that
$\hat{x}:=x/\|x\|_{2}\in\Sigma(A,s,\delta)$ with high probability. Next,
standard tail bounds imply that $\|x\|_{2}^{2}\in\Theta(pn)$ with high
probability. This suggests the scaling $\epsilon:=\alpha/\sqrt{pn}$ with
$\alpha\in(0,1)$. Considering a fraction $\Theta(\alpha)$ of the support of
$Ax$ has magnitude $\Theta(\alpha)$, the hard threshold
$h_{\epsilon}(A\hat{x})$ serves as a decent estimate for the product
$A\hat{x}$:

$\|h_{\epsilon}(A\hat{x})-A\hat{x}\|_{2}^{2}=\Theta(\alpha^{3}).$

Furthermore, we obtain this quality of estimate with relatively little
computation: Since $m=n$, we have $O(n^{3}\log n)$ operations of
preprocessing, and then for each $x$, we compute $h_{\epsilon}(A\hat{x})$ in
$O(\alpha^{-2}sn\log n)$ operations since $\|A\|_{2\to\infty}=1$. By
comparison, if an oracle were to reveal the support of
$h_{\epsilon}(A\hat{x})$, then naive matrix–vector multiplication with the
appropriate submatrix of $A$ would cost $O(sn)$ operations.

###  1.2 Outline

In the next section, we review the necessary theory of projective designs, and
we show how they can be used in conjunction with a median-of-means estimator
to obtain a high-quality random estimate of a sparse matrix–vector product.
Next, Section 3 provides the details of a specific choice of projective
design, namely, one that arises from a Kerdock set described by Calderbank,
Cameron, Kantor, and Seidel [3]. This particular choice of projective design
allows us to leverage the fast Walsh–Hadamard transform to substantially speed
up the preprocessing step of our algorithm. We conclude in Section 4 with some
numerical results that demonstrate the plausibility of a real-world
implementation of our algorithm.

##  2 Projective designs and the median of means

Let $\sigma$ denote the uniform probability measure on the unit sphere
$S^{d-1}$ in $\mathbb{R}^{d}$, let $\operatorname{Hom}_{j}(\mathbb{R}^{d})$
denote the set of homogeneous polynomials of total degree $j$ in $d$ real
variables, and put

$c_{d,k}:=\frac{1\cdot 3\cdot 5\cdots(2k-1)}{d(d+2)\cdots(d+2(k-1))}.$

A projective $t$-design for $\mathbb{R}^{d}$ is defined to be any
$\\{u_{\ell}\\}_{\ell\in[L]}$ in $S^{d-1}$ that satisfies the following
equivalent properties:

######  Proposition 2.

Given $\\{u_{\ell}\\}_{\ell\in[L]}$ in $S^{d-1}$ and $t\in\mathbb{N}$, the
following are equivalent:

  * (a)

$\frac{1}{L}\sum_{\ell\in[L]}p(u_{\ell})=\int_{S^{d-1}}p(u)d\sigma(u)$ for
every $p\in\operatorname{Hom}_{2k}(\mathbb{R}^{d})$ and every
$k\in\\{0,\ldots,t\\}$.

  * (b)

$\frac{1}{L}\sum_{\ell\in[L]}\langle x,u_{\ell}\rangle^{2k}=c_{d,k}\|x\|^{2k}$
for every $x\in\mathbb{R}^{d}$ and every $k\in\\{1,\ldots,t\\}$.

  * (c)

$\frac{1}{L^{2}}\sum_{\ell,\ell^{\prime}\in[L]}\langle
u_{\ell},u_{\ell^{\prime}}\rangle^{2k}=c_{d,k}$ for every
$k\in\\{1,\ldots,t\\}$.

The proof of Proposition 2 is contained in Section 6.4 of [41] and references
therein, given the observation that $\\{u_{\ell}\\}_{\ell\in[L]}$ satisfies
Proposition 2(a) precisely when
$\\{u_{\ell}\\}_{\ell\in[L]}\cup\\{-u_{\ell}\\}_{\ell\in[L]}$ forms a so-
called spherical $2t$-design. We will see that the cubature rule in
Proposition 2(a) is what makes projective $t$-designs useful, while
Proposition 2(c) makes them easy to identify. We note that an analog of
Proposition 2(c) is used to define projective $t$-designs in a variety of
settings, such as complex projective space and the Cayley plane [33].

Our application of projective $2$-designs encourages us to take the size $L$
to be as small as possible. To this end, there is a general lower bound [1] of
$L\geq\binom{d+1}{2}$, but to date, equality is only known to be achieved for
$d\in\\{2,3,7,23\\}$; see [17] and references therein. Despite this scarcity,
there are infinite families of projective $2$-designs that take $L$ to be
slightly larger, specifically, $L=d(d/2+1)$ whenever $d$ is a power of $4$;
see [4, 3]. For these constructions, $\\{u_{\ell}\\}_{\ell\in[L]}$ takes the
form of a union of orthonormal bases. Orthonormal bases
$\\{x_{i}\\}_{i\in[d]}$ and $\\{y_{i}\\}_{i\in[d]}$ are said to be unbiased if
$|\langle x_{i},y_{j}\rangle|^{2}=1/d$ for every $i,j\in[d]$.

######  Proposition 3.

Suppose $\\{u_{b,i}\\}_{i\in[d]}$ in $\mathbb{R}^{d}$ is orthonormal for every
$b\in[d/2+1]$, and suppose further that $\\{u_{b,i}\\}_{i\in[d]}$ and
$\\{u_{b^{\prime},i}\\}_{i\in[d]}$ are unbiased for every
$b,b^{\prime}\in[d/2+1]$ with $b\neq b^{\prime}$. Then
$\\{u_{b,i}\\}_{b\in[d/2+1],i\in[d]}$ forms a projective $2$-design for
$\mathbb{R}^{d}$.

The proof of Propostion 3 follows from Propostion 2(c) and the definition of
unbiased. In Section 3, we will provide an explicit construction of this form.
In the meantime, we show how projective $2$-designs are useful in our
application. The following result defines the random vector $z$ in terms of a
projective $2$-design, and then uses this structure to control the variance of
each coordinate of the random vector $y:=Azz^{\top}x$.

######  Lemma 4.

Given $A\in\mathbb{R}^{m\times n}$, fix a projective $2$-design
$\\{u_{\ell}\\}_{\ell\in[L]}$ for $\mathbb{R}^{d}$ with $d\geq n$, let
$\Pi\colon\mathbb{R}^{d}\to\mathbb{R}^{n}$ denote projection onto the first
$n$ coordinates, and let $z$ denote a random vector with uniform distribution
over $\\{\sqrt{d}\Pi u_{\ell}:\ell\in[L]\\}$. Given a unit vector
$x\in\mathbb{R}^{n}$, define the random vector
$y:=Azz^{\top}x\in\mathbb{R}^{m}$. Then for each $i\in[m]$, it holds that

$\mathbb{E}[e_{i}^{\top}y]=e_{i}^{\top}Ax,\qquad\operatorname{Var}(e_{i}^{\top}y)\leq
2\|A\|_{2\to\infty}^{2}.$

###### Proof.

Let $\Pi^{*}$ denote the adjoint of $\Pi$, namely, the map that embeds
$\mathbb{R}^{n}$ into the first $n$ coordinates of $\mathbb{R}^{d}$. Fix
$i\in[m]$, let $a_{i}^{\top}$ denote the $i$th row of $A$, and put
$\tilde{a}_{i}:=\Pi^{*}a_{i}$ and $\tilde{x}:=\Pi^{*}x$. Then

$e_{i}^{\top}y=e_{i}^{\top}Azz^{\top}x=a_{i}^{\top}(\sqrt{d}\Pi
u_{\ell(z)})(\sqrt{d}\Pi
u_{\ell(z)})^{\top}x=d\tilde{a}_{i}^{\top}u_{\ell(z)}u_{\ell(z)}^{\top}\tilde{x}.$

To compute $\mathbb{E}[e_{i}^{\top}y]$ and $\mathbb{E}[(e_{i}^{\top}y)^{2}]$,
we will consider a random vector $u$ that is uniformly distributed on the unit
sphere $S^{d-1}$, as well as some rotation $Q\in\operatorname{O}(d)$ such that

$Q\tilde{x}=e_{1},\qquad
Q\tilde{a}_{i}=\alpha_{1}e_{1}+\alpha_{2}e_{2},\qquad\alpha_{1},\alpha_{2}\in\mathbb{R}.$

Since $v\mapsto d\tilde{a}_{i}^{\top}vv^{\top}\tilde{x}$ resides in
$\operatorname{Hom}_{2}(\mathbb{R}^{d})$, we may apply Proposition 2(a) to get

$\mathbb{E}[e_{i}^{\top}y]=\frac{1}{L}\sum_{\ell\in[L]}d\tilde{a}_{i}^{\top}u_{\ell}u_{\ell}^{\top}\tilde{x}=\mathbb{E}[d\tilde{a}_{i}^{\top}uu^{\top}\tilde{x}]=\mathbb{E}[d\tilde{a}_{i}^{\top}Q^{\top}Quu^{\top}Q^{\top}Q\tilde{x}].$

A change of variables $Qu\mapsto u$ then gives

$\mathbb{E}[e_{i}^{\top}y]=\mathbb{E}[d(Q\tilde{a}_{i})^{\top}uu^{\top}(Q\tilde{x})]=\mathbb{E}[d(\alpha_{1}e_{1}+\alpha_{2}e_{2})^{\top}uu^{\top}e_{1}]=d\alpha_{1}\mathbb{E}u_{1}^{2}.$

Considering $1=\mathbb{E}[\|u\|^{2}]=d\mathbb{E}u_{1}^{2}$ by symmetry and
linearity of expectation, it follows that

$\mathbb{E}[e_{i}^{\top}y]=\alpha_{1}=(\alpha_{1}e_{1}+\alpha_{2}e_{2})^{\top}e_{1}=\tilde{a}_{i}^{\top}\tilde{x}=a_{i}^{\top}x=e_{i}^{\top}Ax.$

Indeed, this behavior was the original motivation (1) for our approach. We
will apply the same technique to compute $\mathbb{E}[(e_{i}^{\top}y)^{2}]$.
Since $v\mapsto(d\tilde{a}_{i}^{\top}vv^{\top}\tilde{x})^{2}$ resides in
$\operatorname{Hom}_{4}(\mathbb{R}^{d})$, we may apply Proposition 2(a) to get

$\mathbb{E}[(e_{i}^{\top}y)^{2}]=\frac{1}{L}\sum_{\ell\in[L]}(d\tilde{a}_{i}^{\top}u_{\ell}u_{\ell}^{\top}\tilde{x})^{2}=\mathbb{E}[(d\tilde{a}_{i}^{\top}uu^{\top}\tilde{x})^{2}]=\mathbb{E}[(d\tilde{a}_{i}^{\top}Q^{\top}Quu^{\top}Q^{\top}Q\tilde{x})^{2}].$

A change of variables $Qu\mapsto u$ then gives

$\mathbb{E}[(e_{i}^{\top}y)^{2}]=\mathbb{E}[(d(\alpha_{1}e_{1}+\alpha_{2}e_{2})^{\top}uu^{\top}e_{1})^{2}]=d^{2}(\alpha_{1}^{2}\mathbb{E}u_{1}^{4}+\alpha_{2}^{2}\mathbb{E}u_{1}^{2}u_{2}^{2}).$

Next, the theorem in [14] implies $\mathbb{E}u_{1}^{4}=\frac{3}{d(d+2)}$ and
$\mathbb{E}u_{1}^{2}u_{2}^{2}=\frac{1}{d(d+2)}$, thereby implying

$\mathbb{E}[(e_{i}^{\top}y)^{2}]=\alpha_{1}^{2}\cdot\frac{3d}{d+2}+\alpha_{2}^{2}\cdot\frac{d}{d+2}.$

Finally, we recall $\mathbb{E}[e_{i}^{\top}y]=\alpha_{1}$ to compute the
desired variance:

$\operatorname{Var}(e_{i}^{\top}y)=\mathbb{E}[(e_{i}^{\top}y)^{2}]-(\mathbb{E}[e_{i}^{\top}y])^{2}=2\alpha_{1}^{2}\cdot\frac{d-1}{d+2}+\alpha_{2}^{2}\cdot\frac{d}{d+2}\leq
2\alpha_{1}^{2}+\alpha_{2}^{2}\leq 2\|a_{i}\|^{2}.$

The result then follows from the fact that
$\|A\|_{2\to\infty}^{2}=\max_{i\in[m]}\|a_{i}\|^{2}$. ∎

Now that we have control of the variance, we can obtain strong deviation
bounds on our median-of-means estimator $\hat{\mu}$ of $\mathbb{E}y=Ax$.

######  Theorem 5.

Given $A\in\mathbb{R}^{m\times n}$, fix a projective $2$-design
$\\{u_{\ell}\\}_{\ell\in[L]}$ for $\mathbb{R}^{d}$ with $d\geq n$, let
$\Pi\colon\mathbb{R}^{d}\to\mathbb{R}^{n}$ denote projection onto the first
$n$ coordinates, and let $z$ denote a random vector with uniform distribution
over $\\{\sqrt{d}\Pi u_{\ell}:\ell\in[L]\\}$. Given a unit vector
$x\in\mathbb{R}^{n}$, define the random vector
$y:=Azz^{\top}x\in\mathbb{R}^{m}$, select

$J\geq\frac{4e^{2}\|A\|_{2\to\infty}^{2}}{\gamma^{2}},\qquad K\geq
2\log\Big{(}\frac{m}{\eta}\Big{)},$

draw independent copies $\\{y_{jk}\\}_{j\in[J],k\in[K]}$ of $y$, and compute
the entrywise median of means:

$\hat{\mu}:=\operatorname{median}\\{\overline{y}_{k}\\}_{k\in[K]},\qquad\overline{y}_{k}:=\frac{1}{J}\sum_{j\in[J]}y_{jk}.$

Then $\|\hat{\mu}-Ax\|_{\infty}<\gamma$ with probability at least $1-\eta$.

###### Proof.

Fix $i\in[m]$. For notational convenience, we denote the random variable
$Y:=e_{i}^{\top}y$. Given independent copies $\\{Y_{j}\\}_{j\in[J]}$ of $Y$,
let $\overline{Y}$ denote their sample average. Then Chebyshev’s inequality
and Lemma 4 together imply the deviation inequality

$p:=\mathbb{P}\\{|\overline{Y}-(Ax)_{i}|\geq\gamma\\}=\mathbb{P}\\{|\overline{Y}-\mathbb{E}\overline{Y}|\geq\gamma\\}\leq\frac{\operatorname{Var}(\overline{Y})}{\gamma^{2}}\leq\frac{2\|A\|_{2\to\infty}^{2}\|x\|^{2}}{J\gamma^{2}}.$
(2)

Now take independent copies $\\{\overline{Y}_{k}\\}_{k\in[K]}$ of
$\overline{Y}$ and put
$\hat{\mu}_{i}:=\operatorname{median}\\{\overline{Y}_{k}\\}_{k\in[K]}$. Notice
that $\hat{\mu}_{i}\geq(Ax)_{i}+\gamma$ only if half of the
$\overline{Y}_{k}$’s satisfy $\overline{Y}_{k}\geq(Ax)_{i}+\gamma$. Similarly,
$\hat{\mu}_{i}\leq(Ax)_{i}-\gamma$ only if half of the $\overline{Y}_{k}$’s
satisfy $\overline{Y}_{k}\leq(Ax)_{i}-\gamma$. Thus, defining
$I_{k}:=\mathbf{1}\\{|\overline{Y}_{k}-(Ax)_{i}|\geq\gamma\\}$, we have

$\mathbb{P}\\{|\hat{\mu}_{i}-(Ax)_{i}|\geq\gamma\\}\leq\mathbb{P}\bigg{\\{}\sum_{k\in[K]}I_{k}\geq\frac{K}{2}\bigg{\\}}.$

Since $\\{I_{k}\\}_{k\in[K]}$ are independent Bernoulli random variables with
success probability $p$, we may continue with the help of the multiplicative
Chernoff bound:

$\mathbb{P}\bigg{\\{}\sum_{k\in[K]}I_{k}\geq(1+\lambda)Kp\bigg{\\}}\leq\Big{(}\frac{e^{\lambda}}{(1+\lambda)^{1+\lambda}}\Big{)}^{Kp}=e^{-Kp}\Big{(}\frac{e}{1+\lambda}\Big{)}^{(1+\lambda)Kp},\qquad\lambda>0.$

By our choice of $J$, equation (2) implies that $p\leq 1/(2e^{2})<1/2$. As
such, there exists $\lambda>0$ such that $(1+\lambda)Kp=K/2$. Combining the
above bounds then gives

$\mathbb{P}\\{|\hat{\mu}_{i}-(Ax)_{i}|\geq\gamma\\}\leq
e^{-Kp}(2ep)^{K/2}\leq(2ep)^{K/2}\leq\Big{(}\frac{4e\|A\|_{2\to\infty}^{2}\|x\|^{2}}{J\gamma^{2}}\Big{)}^{K/2}\leq
e^{-K/2}\leq\frac{\eta}{m},$

where the last steps follow from our choices for $J$ and $K$. Finally, since
our choice for $i\in[m]$ was arbitrary, the result follows from a union bound.
∎

##  3 Fast preprocessing with Kerdock sets

Select $k\in 2\mathbb{N}$, and consider the real vector space
$\mathbb{R}[\mathbb{F}_{2}^{k}]$ of functions
$f\colon\mathbb{F}_{2}^{k}\to\mathbb{R}$. Calderbank, Cameron, Kantor, and
Seidel [3] describe a projective $2$-design in this space that takes the form
of mutually unbiased orthonormal bases (à la Proposition 3). We explicitly
construct this projective $2$-design with some help from the underlying finite
field, and then we leverage its structure to speed up the preprocessing step
of our algorithm.

We say $M\in\mathbb{F}_{2}^{k\times k}$ is skew-symmetric if $M_{ii}=0$ and
$M_{ij}=M_{ji}$ for every $i,j\in[k]$. Given a skew-symmetric $M$, consider
the corresponding upper-triangular matrix $\tilde{M}\in\mathbb{F}_{2}^{k\times
k}$ defined by $\tilde{M}_{ij}:=M_{ij}\cdot\mathbf{1}\\{i<j\\}$ and the
quadratic form
$Q_{M}\colon\mathbb{F}_{2}^{k}\times\mathbb{F}_{2}^{k}\to\mathbb{F}_{2}$
defined by $Q_{M}(x):=x^{\top}\tilde{M}x$. Given a skew-symmetric
$M\in\mathbb{F}_{2}^{k\times k}$ and $w\in\mathbb{F}_{2}^{k}$, define
$u_{M,w}\in\mathbb{R}[\mathbb{F}_{2}^{k}]$ by

$u_{M,w}(x):=2^{-k/2}(-1)^{Q_{M}(x)+w^{\top}x}.$ (3)

The following result uses the vectors in (3) to form unbiased orthonormal
bases for $\mathbb{R}[\mathbb{F}_{2}^{k}]$; this result is contained in [25,
26, 3], but the proof is distributed over dozens of dense pages from multiple
papers, so we provide a direct and illustrative proof at the end of this
section.

######  Proposition 6.

Consider any skew-symmetric $M,M^{\prime}\in\mathbb{F}_{2}^{k\times k}$.

  * (a)

$\\{u_{M,w}\\}_{w\in\mathbb{F}_{2}^{k}}$ and
$\\{u_{M^{\prime},w}\\}_{w\in\mathbb{F}_{2}^{k}}$ are orthonormal bases.

  * (b)

If $M+M^{\prime}$ has full rank over $\mathbb{F}_{2}$, then
$\\{u_{M,w}\\}_{w\in\mathbb{F}_{2}^{k}}$ and
$\\{u_{M^{\prime},w}\\}_{w\in\mathbb{F}_{2}^{k}}$ are unbiased.

A Kerdock set is a collection $\mathcal{K}\subseteq\mathbb{F}_{2}^{k\times k}$
of $2^{k-1}$ skew-symmetric matrices such that $M+M^{\prime}$ has full rank
for every $M,M^{\prime}\in\mathcal{K}$ with $M\neq M^{\prime}$. (Note that one
cannot hope for a larger set with this property since the first row of each
matrix in $\mathcal{K}$ must be distinct, and the first entry of these rows
must equal zero.) By Proposition 6 and equation (3), a Kerdock set
$\mathcal{K}$ determines $2^{k-1}+1$ mutually unbiased bases in
$\mathbb{R}[\mathbb{F}_{2}^{k}]$, namely,
$\\{u_{M,w}\\}_{w\in\mathbb{F}_{2}^{k}}$ for each $M\in\mathcal{K}$ and the
identity basis $\\{e_{w}\\}_{w\in\mathbb{F}_{2}^{k}}$. By Proposition 3, these
bases combine to form a projective $2$-design. Below, we give the “standard”
Kerdock set described in Example 9.2 in [3]. (Note that [3] contains a typo
and omits the proof of this result; specifically, when they write $\alpha x$,
it should be $ax$; we will use $s$ instead of $a$ so as to clearly distinguish
from $\alpha$, and for completeness, we supply a proof at the end of this
section.) In what follows,
$\operatorname{tr}\colon\mathbb{F}_{2^{k-1}}\to\mathbb{F}_{2}$ denotes the
field trace, while for any finite set $B$, we let $\mathbb{F}_{2}^{B\times B}$
denote the set of $|B|\times|B|$ matrices with entries in $\mathbb{F}_{2}$
whose rows and columns are indexed by $B$.

######  Proposition 7.

Consider the $k$-dimensional vector space
$V:=\mathbb{F}_{2^{k-1}}\times\mathbb{F}_{2}$ over the scalar field
$\mathbb{F}_{2}$, and for each $s\in\mathbb{F}_{2^{k-1}}$, define the linear
map $L_{s}\colon V\to V$ by

$L_{s}(x,\alpha):=(s^{2}x+s\operatorname{tr}(sx)+\alpha
s,\operatorname{tr}(sx)).$

Next, select a basis $B$ for $V$, and consider the bilinear form

$(x,\alpha)\cdot(y,\beta):=\operatorname{tr}(xy)+\alpha\beta.$

Then $\\{M_{s}\in\mathbb{F}_{2}^{B\times B}:s\in\mathbb{F}_{2^{k-1}}\\}$
defined by $(M_{s})_{b,b^{\prime}}:=b\cdot L_{s}(b^{\prime})$ is a Kerdock
set.

Importantly, Kerdock sets provide speedups for the preprocessing step of our
algorithm:

######  Lemma 8.

Given $A\in\mathbb{R}^{m\times n}$, select the smallest $k\in 2\mathbb{N}$
satisfying $d:=2^{k}\geq n$. Let $\mathcal{K}\subseteq\mathbb{F}_{2}^{k\times
k}$ denote the Kerdock set described in Proposition 7, consider
$u_{M,w}\in\mathbb{R}[\mathbb{F}_{2}^{k}]$ defined by (3), let
$\\{e_{w}\\}_{w\in\mathbb{F}_{2}^{k}}$ denote the identity basis in
$\mathbb{R}[\mathbb{F}_{2}^{k}]$, and let
$\Pi\colon\mathbb{R}[\mathbb{F}_{2}^{k}]\to\mathbb{R}^{n}$ denote any
coordinate projection. Then

$\\{A\sqrt{d}\Pi
u_{M,w}\\}_{M\in\mathcal{K},w\in\mathbb{F}_{2}^{k}}\cup\\{A\sqrt{d}\Pi
e_{w}\\}_{w\in\mathbb{F}_{2}^{k}}$

can be computed in $O(mn^{2}\log n+n\log^{2}n)$ operations.

###### Proof.

We identify each member of $\mathbb{R}[\mathbb{F}_{2}^{k}]$ as a vector in
$\mathbb{R}^{d}$ with entries indexed by $\mathbb{F}_{2}^{k}$. By this
identification, the vectors $\\{u_{M,w}\\}_{w\in\mathbb{F}_{2}^{k}}$ appear as
the columns of the matrix product $D_{M}H^{\otimes k}$, where $\cdot^{\otimes
k}$ denotes the Kronecker power and $D_{M},H\in\mathbb{R}^{d\times d}$ are
defined by

$D_{M}:=\operatorname{diag}\\{(-1)^{Q_{M}(x)}\\}_{x\in\mathbb{F}_{2}^{k}},\qquad
H_{ab}:=2^{-1/2}(-1)^{ab},\qquad a,b\in\mathbb{F}_{2}.$

Let $a_{i}^{\top}$ denote the $i$th row of $A$. Then the following algorithm
runs in the claimed number of operations: For each $M\in\mathcal{K}$, (i)
compute $\\{(-1)^{Q_{M}(x)}\\}_{x\in\mathbb{F}_{2}^{k}}$ in $O(n\log^{2}n)$
operations, (ii) for each $i\in[m]$, use the fast Walsh–Hadamard transform to
compute $((\Pi^{*}a_{i})^{\top}D_{M})H^{\otimes k}$ in $O(n\log n)$
operations, and (iii) compute the columns of

$\sqrt{d}\sum_{i\in[m]}e_{i}(a_{i}^{\top}\Pi D_{M}H^{\otimes k})=A\sqrt{d}\Pi
D_{M}H^{\otimes k}$

in $O(mn)$ operations; finally, compute $\\{A\sqrt{d}\Pi
e_{w}\\}_{w\in\mathbb{F}_{2}^{k}}$ in $O(mn)$ operations. ∎

###### Proof of Proposition 6.

Define the Heisenberg group $H_{2k+1}(\mathbb{F}_{2})$ to be the set
$\mathbb{F}_{2}^{k}\times\mathbb{F}_{2}^{k}\times\mathbb{F}_{2}$ with
multiplication

$(p,q,\epsilon)\cdot(p^{\prime},q^{\prime},\epsilon^{\prime}):=(p+p^{\prime},q+q^{\prime},\epsilon+\epsilon^{\prime}+q^{\top}p^{\prime}).$

Fix any skew-symmetric $M\in\mathbb{F}_{2}^{k\times k}$. For each
$w\in\mathbb{F}_{2}^{k}$, we define $\psi_{M,w}\colon\mathbb{F}_{2}^{k}\to
H_{2k+1}(\mathbb{F}_{2})$ by

$\psi_{M,w}(p):=(p,Mp,Q_{M}(p)+w^{\top}p).$

One may apply the polarization identity

$Q_{M}(x+y)=Q_{M}(x)+x^{\top}\tilde{M}y+y^{\top}\tilde{M}x+Q_{M}(y)=Q_{M}(x)+x^{\top}My+Q_{M}(y).$
(4)

to verify that $\psi_{M,w}$ is a group homomorphism.

Next, define the Schrödinger representation $\rho\colon
H_{2k+1}(\mathbb{F}_{2})\to\operatorname{GL}(\mathbb{R}[\mathbb{F}_{2}^{k}])$
by

$\rho(p,q,\epsilon):=(-1)^{\epsilon}X_{p}Z_{q},$

where $X_{p}$ and $Z_{q}$ denote translation and modulation operators,
respectively:

$(X_{p}f)(x):=f(x+p),\qquad(Z_{q}f)(x):=(-1)^{q^{\top}x}f(x).$ (5)

Indeed, $\rho$ is a representation of $H_{2k+1}(\mathbb{F}_{2})$ as a
consequence of the easily verified relation

$Z_{q}X_{p}=(-1)^{q^{\top}p}X_{p}Z_{q}.$ (6)

It follows that $\rho\circ\psi_{M,w}$ is a representation of the additive
group $\mathbb{F}_{2}^{k}$. Explicitly, we have

$(\rho\circ\psi_{M,w})(p)=(-1)^{Q_{M}(p)+w^{\top}p}X_{p}Z_{Mp}.$ (7)

Next, we put

$P_{M,w}:=2^{-k}\sum_{p\in\mathbb{F}_{2}^{k}}(\rho\circ\psi_{M,w})(p).$ (8)

Then $P_{M,w}^{2}=P_{M,w}$ and $P_{M,w}^{*}=P_{M,w}$, and so $P_{M,w}$ is an
orthogonal projection. Decomposing into irreducible representations reveals
that $\operatorname{im}P_{M,w}$ equals the intersection of the eigenspaces of
$\\{(\rho\circ\psi_{M,w})(p)\\}_{p\in\mathbb{F}_{2}^{k}}$ with eigenvalue $1$.

We claim that $P_{M,w}=u_{M,w}u_{M,w}^{*}$. To see this, one may first apply
the definitions (7), (5), and (3) along with (4) and its consequence

$x^{\top}Mx=Q_{M}(x+x)+Q_{M}(x)+Q_{M}(x)=0$ (9)

to verify the eigenvector equation

$[(\rho\circ\psi_{M,w})(p)]u_{M,w}=u_{M,w}$

for each $p\in\mathbb{F}_{2}^{k}$. Next, one may apply the easily verified
fact that

$\operatorname{tr}(X_{p}Z_{q})=\left\\{\begin{array}[]{cl}2^{k}&\text{if
}p=q=0\\\ 0&\text{otherwise}\end{array}\right.$ (10)

to the definitions (8) and (7) to compute $\operatorname{tr}P_{M,w}=1$. This
proves our intermediate claim.

We are now ready to compute $|\langle
u_{M,w},u_{M^{\prime},w^{\prime}}\rangle|^{2}$ in various cases. First, we
cycle the trace:

$|\langle
u_{M,w},u_{M^{\prime},w^{\prime}}\rangle|^{2}=\operatorname{tr}(u_{M,w}^{*}u_{M^{\prime},w^{\prime}}u_{M^{\prime},w^{\prime}}^{*}u_{M,w})=\operatorname{tr}(u_{M,w}u_{M,w}^{*}u_{M^{\prime},w^{\prime}}u_{M^{\prime},w^{\prime}}^{*})=\operatorname{tr}(P_{M,w}P_{M^{\prime},w^{\prime}}).$

Next, we apply the definitions (8) and (7) along with identities (6), (10),
and (9) to obtain

$|\langle
u_{M,w},u_{M^{\prime},w^{\prime}}\rangle|^{2}=\operatorname{tr}(P_{M,w}P_{M^{\prime},w^{\prime}})=2^{-k}\sum_{p\in\operatorname{ker}(M+M^{\prime})}(-1)^{Q_{M}(p)+Q_{M^{\prime}}(p)+(w+w^{\prime})^{\top}p}.$

From here, we proceed in cases: If $M=M^{\prime}$, then
$\operatorname{ker}(M+M^{\prime})=\mathbb{F}_{2}^{k}$ and
$Q_{M}(p)=Q_{M^{\prime}}(p)$, and so $|\langle
u_{M,w},u_{M,w^{\prime}}\rangle|^{2}=\delta_{w,w^{\prime}}$. This gives (a).
For (b), if $M+M^{\prime}$ has full rank, then
$\operatorname{ker}(M+M^{\prime})=\\{0\\}$, and so $|\langle
u_{M,w},u_{M,w^{\prime}}\rangle|^{2}=2^{-k}$, as claimed. ∎

###### Proof of Proposition 7.

In what follows, we demonstrate three things:

  * (i)

For every $s,x\in\mathbb{F}_{2^{k-1}}$ and $\alpha\in\mathbb{F}_{2}$, it holds
that $(x,\alpha)\cdot L_{s}(x,\alpha)=0$.

  * (ii)

For every $s,x,y\in\mathbb{F}_{2^{k-1}}$ and $\alpha,\beta\in\mathbb{F}_{2}$,
it holds that $(x,\alpha)\cdot L_{s}(y,\beta)=L_{s}(x,\alpha)\cdot(y,\beta)$.

  * (iii)

For every $r,s\in\mathbb{F}_{2^{k-1}}$ with $r\neq s$,
$L_{r}(x,\alpha)=L_{s}(x,\alpha)$ implies $(x,\alpha)=(0,0)$.

Thanks to the non-degeneracy of the bilinear form, (i)–(iii) together imply
the result.

First, (i) is easily verified by applying three properties of the trace:
$\operatorname{tr}$ is $\mathbb{F}_{2}$-linear,
$\operatorname{tr}(z^{2})=\operatorname{tr}(z)$, and
$\operatorname{tr}(z)^{2}=\operatorname{tr}(z)$ since
$\operatorname{tr}(z)\in\mathbb{F}_{2}$. Also, (ii) quickly follows from the
linearity of the trace. For (iii), take $r,s\in\mathbb{F}_{2^{k-1}}$ with
$r\neq s$ and suppose $L_{r}(x,\alpha)=L_{s}(x,\alpha)$. The second argument
of this identity is $\operatorname{tr}(rx)=\operatorname{tr}(sx)$. Rearrange
the first argument to get

$0=(r^{2}+s^{2})x+r\operatorname{tr}(rx)+s\operatorname{tr}(sx)+\alpha(r+s)=(r+s)\Big{(}(r+s)x+\operatorname{tr}(rx)+\alpha\Big{)},$

where the last step applies the fact that
$\operatorname{tr}(rx)=\operatorname{tr}(sx)$. Since $r+s\neq 0$ by
assumption, it follows that

$(r+s)x+\operatorname{tr}(rx)=\alpha\in\mathbb{F}_{2}.$ (11)

Since $\alpha^{2}=\alpha$, the left-hand side of (11) satisfies the same
quadratic, which in turn implies $x\in\\{0,(r+s)^{-1}\\}$. Since $x$ also
satisfies $\operatorname{tr}((r+s)x)=0$, it follows that $x=0$. Plugging into
(11) then gives $\alpha=0$, as desired. ∎

##  4 Numerical results

$\|\hat{\mu}-Ax\|_{\infty}$ | $\|\hat{h}-Ax\|_{2}$ | ratio of runtimes
---|---|---
|  |
Figure 1:  Performance of fast sparsifying transform for a random orthogonal
matrix $A\in\mathbb{R}^{n\times n}$ with $n=4096$ and vectors
$x\in\Sigma(A,s,\delta)$ with $s=20$ and $\delta=0$. (left) We compute the
median-of-means estimator $\hat{\mu}\in\mathbb{R}^{n}$ and plot the worst-case
behavior over $1000$ trials. (middle) We compute $\hat{h}\in\mathbb{R}^{n}$ by
multiplying $x$ by the rows of $A$ corresponding to the $10s$ largest-
magnitude entries of $\hat{\mu}$, and we plot the worst-case behavior over
$1000$ trials. (right) We plot the quotient of our algorithm’s runtime with
the runtime of naive matrix–vector multiplication. (Here, we ignore the
runtime of randomly sampling $N$ vectors from the precomputed sketch
$\\{As_{\ell}\\}_{\ell\in[L]}$ since optimizing data structures is beyond the
scope of this paper.) Throughout, gray denotes choices of parameters for which
our algorithm is no faster than the naive algorithm.

In this section, we report the real-world performance of our fast sparsifying
transform. In our experiments, we take $A\in\mathbb{R}^{n\times n}$ to be a
random orthogonal matrix, and then we select a unit vector $x$ such that $Ax$
has exactly $s$ nonzero entries of the same size in random positions. (This is
straightforward to implement since $A^{-1}=A^{\top}$.) Due to limitations in
computing power and storage capacity, we restrict our experiments to dimension
$n=4096$, and we select sparsity level $s=20$.

What follows are some details about our implementation. Since $n$ is a power
of $2$, we have $d=n$, and so our projective $2$-design has size
$L=n(n/2+1)=8392704$. For such a large value of $L$, it turns out that the
runtime of selecting $N$ random members of the precomputed sketch
$\\{As_{\ell}\\}_{\ell\in[L]}$ is sensitive to the design of the underlying
data structure. In order to provide a useful runtime comparison, we therefore
assume that this random selection is performed by an oracle before we start
the runtime clock in our algorithm. To compute medians, we apply the
quickselect algorithm [30]. Finally, to boost performance, we take
$S\subseteq[n]$ to index the $10s$ largest-magnitude entries of $\hat{\mu}$
instead of the top $s$ entries.

The results of our experiments are summarized in Figure 1. Figure 1(left)
illustrates that the median-of-means estimator $\hat{\mu}$ performs better in
practice than predicted by Theorem 5. In particular, we can take $J$ and $K$
to be smaller than suggested by the bounds in our guarantee, which is good for
runtime considerations. (In fact, taking $\eta=1$ in Theorem 5 delivers a
lower bound on $K$ that is greater than $16$, meaning our theoretical
guarantees are off the scale in this plot.) Figure 1(middle) illustrates the
performance of our entire algorithm for different choices of $J$ and $K$.
Notably, we perfectly computed $Ax$ in all of our $1000$ random trials when
$K=2$ and $J=375$. Figure 1(right) illustrates the runtime of our algorithm
relative to naive matrix–vector multiplication. For this plot, we divide the
runtime of our algorithm by the runtime of the naive algorithm. Throughout,
gray denotes choices of $(K,L)$ for which our algorithm provides no speedup
over the naive algorithm. In particular, when $K=2$ and $J=375$, our method is
about twice as fast as the naive approach. Interestingly, the median-of-means
estimator reduces to the empirical mean when $K=2$; this might suggest an
opportunity to improve our theory, and perhaps even speed up our algorithm.

## Acknowledgments

We are grateful to Peter Jung for interesting discussions. TF and FK were
partially supported by the German Research Council (DFG) via contract KR
4512/2-2, DG was partially supported by the German Research Council (DFG) via
contract GR4334/1-1, and DGM was partially supported by AFOSR FA9550-18-1-0107
and NSF DMS 1829955.

## References

  * [1] E. Bannai, S. G. Hoggar, On tight $t$-designs in compact symmetric spaces of rank one, Proc. Japan Acad., 61, Ser. A (1985) 78–82.
  * [2] A. M. Bruckstein, D. L. Donoho, M. Elad, From sparse solutions of systems of equations to sparse modeling of signals and images, SIAM Rev. 51 (2009) 34–81.
  * [3] A. R. Calderbank, P. J. Cameron, W. M. Kantor, J. J. Seidel, $\mathbb{Z}_{4}$-Kerdock codes, orthogonal spreads, and extremal Euclidean line-sets, Proc. London Math. Soc. 75 (1997) 436–480.
  * [4] P. J. Cameron, J. J. Seidel, Quadratic forms over GF(2), Geometry and Combinatorics, Academic Press, 1991, 290–297.
  * [5] B. Choi, A. Christlieb, Y. Wang, High-dimensional sparse Fourier algorithms, Numer. Algor. 87 (2021) 161–186.
  * [6] B. Choi, A. Christlieb, Y. Wang, Multiscale High-Dimensional Sparse Fourier Algorithms for Noisy Data, arXiv:1907.03692
  * [7] B. Choi, M. Iwen, F. Krahmer. Sparse harmonic transforms: A new class of sublinear-time algorithms for learning functions of many variables, arXiv:1808.04932
  * [8] B. Choi, M. Iwen, T. Volkmer. Sparse Harmonic Transforms II: Best $s$-Term Approximation Guarantees for Bounded Orthonormal Product Bases in Sublinear-Time, arXiv:1909.09564
  * [9] A. Christlieb, D. Lawlor, Y. Wang, A multiscale sub-linear time Fourier algorithm for noisy data, Appl. Comput. Harmon. Anal. 40 (2016) 553–574.
  * [10] D. Coppersmith, S. Winograd, Matrix multiplication via arithmetic progressions, J. Symbolic Comput. 9 (1990) 251–280.
  * [11] A. M. Davie, A. J. Strothers, Improved bound for complexity of matrix multiplication, Proc. Roy. Soc. Edinburgh Sect. A 143 (2013) 351–369.
  * [12] R. Drineas, R. Kannan, M. W. Mahoney, Fast Monte Carlo algorithms for matrices I: Approximating matrix multiplication, SIAM J. Comput. 36 (2006) 132–157.
  * [13] M. Elad, M. Aharon, Image denoising via sparse and redundant representations over learned dictionaries, IEEE Trans. Image Process. 15 (2006) 3736–3745.
  * [14] G. B. Folland, How to integrate a polynomial over a sphere, Amer. Math. Monthly 108 (2001) 446–448.
  * [15] A. C. Gilbert, S. Guha, P. Indyk, S. Muthukrishnan, M. Strauss, Near-Optimal Sparse FourierRepresentations via Sampling, STOC 2002, 152–161.
  * [16] A. Gilbert, S. Muthukrishnan, and M. Strauss, Improved time bounds for near-optimal sparse Fourier representations, Proc. SPIE (2005).
  * [17] N. I. Gillespie, Equiangular lines, incoherent sets and quasi-symmetric designs, arXiv:1809.05739.
  * [18] H. Hassanieh, P. Indyk, D. Katabi, E. Price, Nearly optimal sparse Fourier transform, STOC 2012, 563–577.
  * [19] H. Hassanieh, P. Indyk, D. Katabi, E. Price, Simple and practical algorithm for sparse Fourier transform, SODA 2012, 1183–1194.
  * [20] P. Indyk, M. Kapralov, Sample-optimal Fourier sampling in any constant dimension, FOCS 2014, 514–523.
  * [21] M. Iwen, Combinatorial sublinear-time Fourier algorithms, Found. Comp. Math. 10 (2010) 303–338.
  * [22] M. Iwen, Improved approximation guarantees for sublinear-time Fourier algorithms, Appl. Comp. Harmon. Anal. 34 (2013) 57–82.
  * [23] L. Kämmerer, F. Krahmer, T. Volkmer, A sample efficient sparse FFT for arbitrary frequency candidate sets in high dimensions, arXiv:2006.13053
  * [24] L. Kämmerer, D. Potts, T. Volkmer, High-dimensional sparse FFT based on sampling along multiple rank-1 lattices, Appl. Comput. Harmon. Anal. 51 (2021) 225–257.
  * [25] W. M. Kantor, Spreads, translation planes and Kerdock sets. I, SIAM J. Alg. Discr. Meth. 3 (1982) 151–165.
  * [26] W. M. Kantor, Spreads, translation planes and Kerdock sets. II, SIAM J. Alg. Discr. Meth. 3 (1982) 308–318.
  * [27] D. Lawlor, Y. Wang, A. Christlieb, Adaptive sub-linear time Fourier algorithms, Adv. Adapt. Data Anal. 5 (2013) 1350003.
  * [28] F. Le Gall, Powers of tensors and fast matrix multiplication, ISSAC 2014, 296–303.
  * [29] E. Liberty, S. W. Zucker, The mailman algorithm: A note on matrix-vector multiplication, Inform. Process. Lett. 109 (2009) 179–182.
  * [30] H. M. Mahmoud, R. Modarres, R. T. Smythe, Analysis of quickselect: an algorithm for order statistics, RAIRO-Theor. Inf. Appl. 29 (1995) 255–276.
  * [31] J. Mairal, F. Bach, J. Ponce, Task-driven dictionary learning, IEEE Trans. Pattern Anal. Mach. Intell. 34 (2012) 791–804.
  * [32] L. Morotti, Explicit universal sampling sets in finite vector spaces, Appl. Comput. Harmon. Anal. 43 (2017) 354–369.
  * [33] A. Munemasa, Spherical designs, Handbook of Combinatorial Designs (2007) 637–643.
  * [34] S. Nam, M. E. Davies, M. Elad, R. Gribonval, The cosparse analysis model and algorithms, Appl. Comput. Harmon. Anal. 34 (2013) 30–56.
  * [35] D. Potts, T. Volkmer, Fast, exact and stable reconstruction of multivariate algebraic polynomials in Chebyshev form, 2015\.
  * [36] S. Ravishankar, Y. Bresler, Learning sparsifying transforms, IEEE Trans. Signal Process. 61 (2012) 1072–1086.
  * [37] C. Rusu, N. González-Prelcic, R. W. Heath, Fast orthonormal sparsifying transforms based on Householder reflectors, IEEE Trans. Signal Process. 64 (2016) 6589–6599.
  * [38] C. Rusu, J. Thompson, Learning fast sparsifying transforms, IEEE Trans. Signal Process. 65 (2017) 4367–4378.
  * [39] V. Strassen, Gaussian elimination is not optimal, Number. Math. 13 (1969) 354–346.
  * [40] I. Tošić, P. Frossard, Dictionary learning, IEEE Signal Process. Mag. 28 (2011) 27–38.
  * [41] S. Waldron, An Introduction to Finite Tight Frames, Applied and Numerical Harmonic Analysis, Springer, 2018.
  * [42] R. Williams, Matrix-vector multiplication in sub-quadratic time (some preprocessing required), SODA 2007, 995–1001.
  * [43] V. V. Williams, Multiplying matrices faster than Coppersmith-Winograd, STOC 2012, 887–898.