# Fast Sampling Based Sketches for Tensors††thanks: Research of William
Swartworth and David P. Woodruff was supported by a Simons Investigator Award,
a Google Faculty Award, and NSF Grant No. CCF-2335412.

William Swartworth (CMU<EMAIL_ADDRESS>David P. Woodruff (CMU)
<EMAIL_ADDRESS>
###### Abstract

We introduce a new approach for applying sampling-based sketches to two and
three mode tensors. We illustrate our technique to construct sketches for the
classical problems of $\ell_{0}$ sampling and producing $\ell_{1}$ embeddings.
In both settings we achieve sketches that can be applied to a rank one tensor
in $(\mathbb{R}^{d})^{\otimes q}$ (for $q=2,3$) in time scaling with $d$
rather than $d^{2}$ or $d^{3}$. Our main idea is a particular sampling
construction based on fast convolution which allows us to quickly compute sums
over sufficiently random subsets of tensor entries.

##  1 Introduction

In the modern area of enormous data sets, space is often at a premium, and one
would like algorithms that either avoid storing all available data, or that
compress existing data. A common and widely-applied strategy is sketching.
Given a vector $x\in\mathbb{R}^{n}$ consisting of the relevant data, a
(linear) sketch of $x$ is given by $Sx$ where $S$ is a linear map down to a
dimension much smaller than $n$. Typically the goal is to design $S$ so that
some useful statistic of $x$ can be computed from the sketched vector $Sx$,
even though most of the information about $x$ has been discarded. Linear
sketches are particularly useful in the context of streaming algorithms, since
linear updates to $x$ can be translated to the sketch, simply by sketching the
vector of updates, and adding it to the previous value of the sketch. Sketches
have also found important applications in speeding up linear-algebraic (and
related) computations [Woo+14]. Here, the idea is to first apply a sketch to
reduce the dimensionality of the problem, while approximately preserving a
quantity of interest (e.g., a regression solution). Then standard algorithms
can be applied to the smaller problem resulting in improved runtimes.

A lot of work on sketching has focused on efficiently applying sketches to
structured data. For example, if the underlying data is sparse, one might hope
for a sketch that can be applied in input-sparsity time. A different type of
structure that has been widely studied in this context is tensor structure. A
line of work has been devoted to developing fast tensor sketches [PP13,
Ahl+20, AK19, MSW19], which are sketches that can be applied quickly to low-
rank tensors. A rank-one tensor in $\mathbb{R}^{n}\otimes\mathbb{R}^{n}$ for
example (i.e., a rank-one matrix) only requires $O(n)$ parameters to specify.
However, a naive sketch might require $O(n^{2})$ time to apply if each entry
of the tensor must be calculated to form the sketch. The goal with tensor
sketches is to do better then expanding into a vector form prior to sketching.

Much of the work based on sketching tensors has focused on $\ell_{2}$-norm
settings. The most studied example is Johnson-Lindenstrauss (JL) embeddings,
which have been studied extensively for tensors. For JL sketches, the ideal
sketch is a dense Gaussian sketch. However, such sketches are slow to apply
and so the main interest is in approximating the properties of a Gaussian
sketch with a sketch that can be applied faster.

We study a completely different type of “sampling-based” sketch, which is
embodied by the $\ell_{0}$ sampling problem. The goal here is to observe a
linear sketch of $x$, and then to (nearly) uniformly sample an element from
$x$’s support. A dense Gaussian sketch is completely useless as we would like
a row of our sketch to single out an element of $\textnormal{supp}(x)$. To
allow for this, the standard idea is to take a sketch which performs sparse
samples at various scales. We ask: Are there sampling-based sketches which can
be applied quickly to rank-one tensors?

For tensors with two and three modes (which cover many tensors of practical
importance), we provide a new approach for constructing sampling-based
sketches in the tensor setting. To illustrate our approach, we focus on two
fundamental problems: $\ell_{0}$ sampling, and $\ell_{1}$ embeddings. We
recall the setup for these problems.

####  $\ell_{0}$ Sampling.

For the $\ell_{0}$ sampling problem one would like to construct a distribution
over sketching matrices $\mathbf{S}$ so that observing $\mathbf{S}x$ allows us
to return a uniformly random element from the support of $x$. We allow
constant distortion, and $\delta$ failure probability, so conditioned on not
failing (which must occur with probability at least $1-\delta$), we should
have that for all $i$ in the support of $x$, the sampling algorithm outputs
$i$ with probability in
$[c_{1}/|\textnormal{supp}(x)|,c_{2}/|\textnormal{supp}(X)|]$.

####  $\ell_{1}$ Embeddings.

For the $\ell_{1}$ embedding problem, the goal is to construct a sketch
$\mathbf{S}$, such that for any $x\in\mathbb{R}^{n}$, we have

$c_{1}\left\|x\right\|_{1}\leq\left\|\mathbf{S}x\right\|_{1}\leq
c_{2}\left\|x\right\|_{1}$

for absolute constants $0<c_{1}<c_{2}.$

###  1.1 Our Results

Our main idea is constructing a way to subsample a tensor at a given sampling
rate $p$, and then to sum over the resulting sampling. Specifically we
introduce a sampling primitive that we call a “$p$-sample” which samples each
element with probability approximately $p$, and does so in a nearly pairwise
independent manner. As we show, this will be sufficient for constructing both
$\ell_{0}$-samplers and $\ell_{1}$ embeddings. Constructing $p$-samples is
straightforward for tensors which are given entrywise – one can simply sample
each entry independently with probability $p.$ Our key novelty is in showing
that it is possible to construct $p$-samples for two and three mode tensors,
which for rank-one tensors can be summed over in nearly linear time. We
discuss the ideas for constructing samples below.

Given our sampling primitive, we show how to construct $\ell_{0}$ samplers and
$\ell_{1}$ embeddings.

For the $\ell_{0}$ sampling problem we show the following result:

######  Theorem 1.1.

For $q=2,3$ there there is a linear sketch of $X\in\mathbb{R}^{n^{q}}$ with
sketching dimension $m=O(\log\frac{1}{\delta}\log^{2}n(\log\log
n+\log\frac{1}{\delta}))$ space, and a sampling algorithm that succeeds with
probability $1-\delta$, which, conditioned on succeeding, outputs an index $i$
in $\textnormal{supp}(X).$ Conditioned on succeeding, the probability that the
algorithm outputs $i\in\textnormal{supp}(X)$ is in
$[\frac{0.5}{|\textnormal{supp}(X)|},\frac{2}{|\textnormal{supp}(X)|}].$ The
algorithm also returns the value of $X_{i}.$ Moveover, the entries of the
sketching matrix can be taken to be in $\\{0,+1,-1\\}.$

For $q=2$, our sketch can be applied to rank-one tensors in $O(mn)$ time, and
when $q=3$ can be applied in $O(mn\log^{2}n)$ time.

For $\ell_{1}$ embeddings, we show the following:

######  Theorem 1.2.

There is an $O(1)$-distortion $\ell_{1}$ embedding sketch $\mathbf{S}$
$\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ with sketching dimension
$m=O(\log^{4}n+\log^{2}(1/\delta)\log n)$ that satisfies
$\left\|\mathbf{S}x\right\|_{1}\leq c\left\|x\right\|_{1}$ with constant
probability, and satisfies $\left\|\mathbf{S}x\right\|_{1}\geq
c\left\|x\right\|_{1}$ with probability at least $1-\delta.$

Moreover our sketch can be applied to rank one tensors in
$\mathbb{R}^{k}\otimes\mathbb{R}^{k}$ in $O(mk)$ time, and to rank one tensors
in $\mathbb{R}^{k}\otimes\mathbb{R}^{k}\otimes\mathbb{R}^{k}$ in
$O(mk\log^{2}k)$ time.

Our main novelty here lies in the $p$-sample construction, which is applied in
a similar way as for our $\ell_{0}$-sample, although the details are more
complicated. We therefore defer the proof of this result to the appendix.

#### Regression

We note that this error guarantee of “no contraction” with high probability,
and “no dilation” with constant probability is know to be sufficient to
construct algorithms for solving $\ell_{1}$ regression. In particular, by
setting $\delta=\exp(-d)$, a standard net argument given in [CW14] for
example, shows that “no contraction” holds with high probability over a
$d$-dimensional subspace. As long as “no dilation” holds for the solution
vector, then it is well-known that this yields a dimension-reduction for
$\ell_{1}$ regression. Thus our $\ell_{1}$ sketch can be used to speed up
sketching for an $\ell_{1}$ regression problem of the form
$\min_{x}\left\|Ax-b\right\|_{1},$ where $A$ has $d$ columns, each of which
have low-rank structure.

Finally, while we choose to illustrate our approach on the problems of
$\ell_{0}$ sampling and $\ell_{1}$ embeddings, we believe that this technique
could be applied more generally whenever there is a known sketch that is built
from random subsampling and taking random linear combinations within the
sample.

###  1.2 Techniques

The main idea behind both $\ell_{0}$-sampling and $\ell_{1}$-embeddings is to
subsample coordinates at various scales in order to isolate some coordinates
of $x.$ As a warmup, we begin by describing $\ell_{0}$ sampling. Then we
describe how our techniques extend to constructing $\ell_{1}$ embeddings. We
note that our basic approach to constructing $\ell_{0}$ samplers and
$\ell_{1}$ embeddings are not novel. However our ability to quickly sketch
rank-one tensors is.

####  $\ell_{0}$ Sampling

In the general form of $\ell_{0}$ sampling, we are given a vector (or in our
case a tensor) $X$, and the goal is to design a sketch that allows us sample
an entry of $X$ nearly uniformly from the support of $X.$ The idea is to take
a random sample $S$ of $X$’s entries and then to store a random linear
combination of the values in $S$ as an entry of the sketch. If $S$ intersects
$\textnormal{supp}(X)$ in a single element then this will allow us to recover
a value in $\textnormal{supp}(X)$. If $\textnormal{supp}(X)$ has size $k$,
then we would like $S$ to sample each value with probability roughly $1/k$ in
order to have a good probability of isolating a single element. Since $k$ is
unknown to us initially, the idea is to perform the sampling procedure at each
of $\log n$ sampling levels $1,\frac{1}{2},\frac{1}{4},\ldots$.

The challenge for us is to design samples $S$ that allow for fast sketching on
rank-one tensors, e.g., tensors of the form $x\otimes y$. As we later show,
the step of computing a random linear combination of values in $S$ can be
achieved simply by randomly flipping the signs along each mode. The harder
part is designing $S$ so that we can sum the resulting tensor over $S$ in
roughly $O(n)$ time. In other words we need $S$ to be such that we can quickly
compute $\sum_{(i,j)\in S}x_{i}y_{j}.$ If $S$ just samples each entry of $X$
i.i.d. with probability $p$, then it is not clear how to do better than
computing the sum term-by-term, which would require $\Omega(pn^{2})$ time in
our example.

The next most natural thing to try is to take $S=S_{1}\times S_{2}$ to be a
random rectangle, since rectangles allow the above summation to be computed in
$O(n)$ time. To achieve sampling probability $p$ one could take a rectangle of
dimension $\sqrt{p}d\times\sqrt{p}d,$ where the subsets of indices $S_{1}$ and
$S_{2}$ corresponding to each dimension are chosen randomly. Thus,
construction succeeds in sampling each entry of $X$ with probability $p$;
however this is not sufficient to isolate a single entry of
$\textnormal{supp}(X)$ with good probability. For example
$\textnormal{supp}(X)$ could consist of a single row of $X$. If we set $p=1/d$
then $S$ has a good probability of sampling some entry from
$\textnormal{supp}(X)$. However, when this occurs, $S$ nearly always contains
more than one such element. Indeed, $S_{2}$ would typically contain around
$\sqrt{d}$ elements, and in this situation it is not possible to sample a
single entry from a row of $X.$ One could try to fix this by randomizing the
sampling scale on each side of the rectangle for a fixed sampling probability.
In the case of $\ell_{0}$ sampling this would at least result in worse
logarithmic factors in the sketching dimension than what we achieve. In the
case of $\ell_{1}$ embeddings, it is unclear how to obtain a constant-factor
embedding with this approach. For instance, in the two-mode case, a given
sampling probability $p$ can be realized by $O(\log d)$ choices of scales for
the two rectangle dimensions. If a given level of $X$ has size $1/p$ and is
distributed uniformly throughout $X$, then all rectangles with sampling
probability $p$ will pick out single elements of that level. On the other
hand, if the level is distributed along a single $1$-dimensional fiber, then
it might be the case that only one of the rectangles is good for that level.
This suggests that the distortion would likely need to be $\Omega(\log d)$ on
some inputs. It does not seem clear how one could arrange arrange for an
$O(1)$ distortion $\ell_{1}$ embedding that handles these situations
simultaneously.

Surprisingly, at least for $2$ and $3$ mode tensors, it is possible to
construct samples $S$ that have better sampling properties than rectangles,
but which can still be summed over in nearly linear time. For brevity, we
refer to the sampling that we need as a $p$-sample, when it samples each index
with probability (approximately) $p.$ The idea is to design $S$ in such a way
that we can employ algorithms for fast convolution using the Fast Fourier
Transform. For example, in the $2$-mode case, one can compute

$\sum_{i+j\in[0,T-1]}x_{i}y_{j}$

for a constant $T$ (chosen to give the appropriate sampling probability) in
$O(n\log n)$ time, by first calculating the convolution $x*y$ and then summing
over the indices in $[0,T-1].$ This only gives a sum over a fixed set, but by
randomly permuting the indices along each mode, and choosing $T$
appropriately, this allows us to achieve sampling probabilities down to
$p=1/n.$ Smaller sampling probabilities result in samples of size $O(n)$, and
so the sum can just be calculated explicitly. (As it turns out, in the
$2$-mode case, the runtime can be improved to $O(n)$ by a simple
optimization.)

The $3$-mode case is somewhat more complicated. A similar convolution trick
involving a convolution of three vectors allows us to construct $p$-samples
with fast summation down to $1/n.$ And sampling probabilities below $1/n^{2}$
allow for direct computation. What about sampling probabilities between
$1/n^{2}$ and $1/n$? Here, we show how to compute sums of the form

$\sum_{i+j+k=0,j-i\in[0,T-1]}x_{i}y_{j}z_{k}$

in $O(n\log^{2}T)$ time. The rough approach is to show that we can reduce
computing this quantity to $n/T$ instances of multiplication by the top-left
corner ($\\{i,j:i\geq j\\}$) of a roughly $T\times T$ Toeplitz matrix. Each
top corner matrix can be decomposed into a sum of $T\log T$ Toeplitz matrices
(up to some zero-padding) each of which has an $O(T\log T)$ multiplication
time.

####  $\ell_{1}$ embeddings

For $\ell_{1}$ embeddings our approach roughly follows the sketch introduced
by [VZ12] for producing constant distortion $\ell_{1}$ embeddings. The main
idea is to think of a vector $x\in\mathbb{R}^{n}$ as decomposed into
approximate level sets, with the size of entries in each level set decreasing
exponentially. For instance, when $\left\|x\right\|_{1}=1$, the level sets
$L_{i}$ could be taken to contain entries in $[q^{-(i+1)},q^{-i}]$ for some
$q<1.$ For each level set, our sketch has an associated level with sampling
probability larger than $q$, which is designed to capture most of the mass of
$L_{i}$. Since we would like our sampling level to capture the mass of $L_{i}$
with high probability, we need to oversample the entries of $L_{i}$
substantially. The usual approach here is to choose a fairly large $q$, and
then hash the sampled values into separate buckets so as to minimize
cancellations. This is typically done with CountSketch; however it is not
clear how to efficiently compose a CountSketch with our fast sample
constructions. Therefore we take a slightly different approach. Instead of
hashing each sampling level into $T$ buckets, we simply take $T$ independent
samples, each with sampling probability $q/T.$. The analysis is quite similar
to the analysis using CountSketch, but we retain the the ability to quickly
sketch rank-one tensors. However, we do pay a price – our sketching time
scales linearly with the size of the sketch. In contrast, for vectors,
CountSketch constructions allow for sketching in time that scales as
$O(n\cdot({\text{number of sampling levels}}))$. However these sketches are
not efficient for rank-one tensors, since $n$ would be replaced by $n^{2}$ or
$n^{3}.$

To avoid too much dilation, the rough idea is that each sampling level picks
up roughly the right amount of mass from its corresponding $L_{i}$ with good
probability. Also, large levels have a small numbers of elements and so they
are unlikely to be picked up by a smaller sampling level. The main issue is in
bounding the contribution of the $L_{i}$’s containing many small values. To
avoid, this [VZ12] introduced the idea of hashing each sampling level into
buckets with random signs in order to induce cancellation among the
contributions from such $L_{i}$’s. We employ the same technique here. However
in order to make applying the signs efficient, our approach is to first apply
random signs along the modes of the tensor, and then to sum over appropriate
$p$-samples in order to compute a “bucket”. For a constant number of modes,
this turns out to be enough randomness to induce cancellation.

This outline for constructing $\ell_{1}$ embeddings has appeared in the
literature several times. Our main novelty lies in constructing samples that
admit fast linear combinations for rank-one tensors.

Finally, one might wonder if our techniques are really necessary. We could
hope to design fast sketches using Kronecker-structured-sketches of the form
$S_{1}\otimes S_{2}$, which are easy to apply to a rank-one tensor $x\otimes
y$ – one simply computes $S_{1}x\otimes S_{2}y.$ It does not appear that a
Kronecker-structured sketch can match our bounds. Indeed by a lower bound
given in [LWY21] , an $O(1)$ $\ell_{1}$ embedding for a single vector
$x\in\mathbb{R}^{n}$ requires at least $\operatorname{poly}(n)$ sketching
dimension for Kronecker-structured sketches. On the other hand, our sketch can
still be applied to rank-one tensors in near-linear time, but requires only
$\operatorname{poly}\log n$ space to embed $x.$ One could also ask whether
Kahtri-Rao sketches (i.e., a sketch with each row of the sketch a rank one
tensor) could be applied to match our sketching dimension. While we do not
provide a lower bound against Kahtri-Rao measurements, we also do not see how
to match our bounds for either $\ell_{0}$ sampling or $\ell_{1}$ embeddings.
We leave the question of whether such a lower bound exists as an interesting
open question.

For larger $q$ we also note that our $\ell_{1}$ embedding procedure could be
applied to triples of modes at a time, giving a tree construction similar to
[Ahl+20]. Unfortunately however, the distortion grows like $2^{\log_{3}q}$,
and the success probability becomes $c^{-q}$ by a union bound. Whether these
parameters can be improved while allowing for fast sketching is an interesting
problem.

###  1.3 Additional Open Questions

An interesting open question is whether a similar set of ideas can be applied
to tensors with a larger number of modes. For four modes, the most natural
approach here would be to attempt to develop a fast sum that can sum over a
(subset of) a codimension $2$ subspaces of the tensor entries. For example, we
might wish to calculate sums of the form

$\sum_{i,j,k,\ell\in Q}w_{i}x_{j}y_{k}z_{\ell}$

where $Q$ is a set of entries satisfying two linear constraints.
Unfortunately, we are not aware of a way to calculate such a sum faster than
$O(n^{2})$ time. It would be interesting to either give a fast summation
algorithm, or to find a new technique that gets around this issue.

###  1.4 Related Work

A substantial literature has been devoted to obtaining fast sketches for
tensors. Tensor sketching was initiated in [PP13] where it was shown that
CountSketch can be quickly applied to rank-one tensors using the Fast Fourier
Transform. More recently, $\ell_{2}$ embeddings were constructed by [AK19] and
improved in [Ahl+20] with applications for sketching polynomial kernels.

In earlier work [IM08] gives a Kronecker-structured for the $\ell_{1}$ norm in
the context of independence testing, using a variation on the previously-known
Cauchy sketch [Ind06]. We note that this sketch requires taking a median in
order to recover the $\ell_{1}$ norm, and thus does not give an $\ell_{1}$
embedding, which may be more suitable for optimization applications.

[Ind06] studies the problem of $\ell_{1}$ estimation using Cauchy sketches.
Later, [VZ12] gave a construction of an oblivious $\ell_{1}$ embedding. Our
general approach for constructing $\ell_{1}$ embeddings is largely based on
theirs. This approach is also generalized in [CW14] to M-estimators, and in
particular is applied to construct $\ell_{1}$ (and more general) subspace
embeddings. These bounds for $\ell_{1}$ embeddings were recently improved in
[MOW21].

[LWY21] expanded on the work and considered independence testing for higher
order tensors. They also give a $\operatorname{poly}(d)$ lower bound for
constructing $\ell_{1}$ embeddings for a single vector in $\mathbb{R}^{d}$
using Kronecker-structured measurements.

###  1.5 Notation and Preliminaries

We use the notation $[N]$ to refer to the set $\\{1,\ldots,N\\}.$ For two
vectors $x$ and $y$, the notation $x\otimes y$ refers to the Kronecker product
of $x$ and $y.$ We use the notation $x*y$ for the circular convolution of $x$
and $y$. That is

$(x*y)_{k}=\sum_{i+j=k}x_{i}y_{j},$

where the sum $i+j$ is interpreted mod $N.$ In such situations it is sometimes
convenient to index vectors starting at $0$ so we will do this occasionally.
However unless otherwise stated, we index starting at $1.$

We typically identify a $q$-mode tensor with a vector in $\mathbb{R}^{n^{q}}$,
assuming that the dimension along each tensor mode is $n.$ We will typically
make this assumption for the sake of convenience. We say that a tensor in
$\mathbb{R}^{n^{q}}$ has rank one if it can be expressed as
$x_{1}\otimes\ldots\otimes x_{q}$ for some vectors $x_{1},\ldots,x_{q}.$

We will use $c$ throughout to refer to an absolute constant, which might be
different between uses (even within an equation).

The notation $\tilde{O}(f)$ means $O(f\log^{c}f)$ for some constant $c.$

####  $\ell_{0}$ sampling.

$\ell_{0}$ sampling has received extensive attention. [CF14] surveys the
standard recipe for constructing $\ell_{0}$-samplers, which we roughly follow
here.

[WZ18] considers solving the $\ell_{0}$-sampling problem in a distributed
setting for a product of two matrices which are held on different servers,
however this is different from our setting, as their communication scheme is
not a linear sketch.

##  2 $p$-sample constructions

We first define our main sampling primitive, which we call a $p$-sample. This
can be viewed as a slightly weaker version of a pairwise independent sample.

######  Definition 2.1.

Let $T$ be an arbitrary finite set. We say that a random subset of $S$ of $T$
is a $p$-sample for $T$ if

  1. 1.

For all entries $i\in T$,

$p/2\leq\operatorname{Pr\/}(i\in S)\leq p.$

  2. 2.

For all $j\neq i$ in $T$,

$\operatorname{Pr\/}(j\in S|i\in S)\leq 2p.$

######  Definition 2.2.

Fix an $n$ and $m$ and suppose that $S$ is a subset of
$[n]^{m}=[n]\times\ldots\times[n]$. Let
$x^{(1)},\ldots,x^{(m)}\in\mathbb{R}^{n}$ be arbitrary vectors. We say that
$S$ admits fast summation if there is an algorithm which computes
$\sum_{(i_{1},\ldots,i_{m})\in S}x_{i_{1}}^{(1)}x_{i_{2}}^{(2)}\ldots
x_{i_{m}}^{(m)}$ in $\tilde{O}(mn)$ time.

To perform one of our constructions we will need a particular type of
Toeplitz-like fast matrix multiplication.

######  Lemma 2.3.

Let $A\in\mathbb{R}^{n\times n}$ be a matrix where $A_{i,j}=A_{i+1,j+2}$ for
all $0\leq i\leq n-1$, $0\leq j\leq n-2.$ Let $B$ be defined by
$B_{i,j}=A_{i,j}$ if $j\geq i$ and $B_{i,j}=0$ otherwise.

Given $v\in\mathbb{R}^{n}$ the matrix product $Av$ can be computed in $O(n\log
n)$ time, and the product $Bv$ can be computed in $O(n\log^{2}n)$ time.

###### Proof.

To see that $A$ admits fast multiplication, note that the rows of $A$ coincide
with the even rows of a Toeplitz matrix of size $2n\times 2n.$ Since
multiplication by a Toeplitz matrix can be carried out in $O(n\log n)$ time
[KS99], the same holds for $A.$

To get a fast algorithm for $B$ we decompose it into a sum of matrices, each
with the same structure of $A$ up to some zero-padding. We call these matrices
$B_{1},B_{2}$ and $B_{3}.$ Take $B_{1}$ to be $B$ but with all indices outside
of $\\{(i,j):1\leq i\leq\lceil n/2\rceil,\lceil n/2\rceil\leq j\leq n\\}$
replaced with $0.$ Visually, the support of $A$ is a right triangle and the
support of $B_{1}$ corresponds to the largest square inscribed in $A.$ Let
$B_{2}$ and $B_{3}$ correspond to the two triangles that are left after
removing the square. That is, $B_{2}$ has support contained in
$\\{(i,j):i>\lfloor\frac{n}{2}\rfloor\\}$ and $B_{3}$ has support contained in
$\\{(i,j):j>\lfloor\frac{n}{2}\rfloor\\}.$

By construction, $B=B_{1}+B_{2}+B_{3}.$ Then $B_{1}$ has the structure of $A$
from the first part of the lemma, so admits $O(n\log n)$ time multiplication.
After removing the zero padding, both $B_{2}$ and $B_{3}$ have the same
structure as $B$ but with half the dimension. Thus we have a runtime
recurrence of the form $T(n)=2T(\lfloor n/2\rfloor)+O(n\log n)$, which gives
an overall runtime of $O(n\log^{2}n).$

∎

######  Theorem 2.4.

For $m=2,3$, and for all $n$ and $p$, there exists a $p$-sample for $[n]^{m}$
that admits fast summation. Additionally

  1. 1.

When $m=2$ the summation runs in $O(n)$ time

  2. 2.

When $m=3$ the summation runs in $O(n\log^{2}n)$ time

###### Proof.

We first show how to construct sets of indices that can be summed over
quickly. We will then show how to use these sets to construct $p$-samples with
fast summation.

m=2 Constructions. Let $T$ be an arbitrary positive integer, and consider the
set

$A_{T}=\\{(i,j)\in[n]\times[n]:i+j\in[T]\\}\subseteq(\mathbb{Z}/n)^{2}.$

(Note that we are treating all indices as values in $\mathbb{Z}/n.$) We will
first show that $A_{T}$ admits $O(n)$ summation for all $T.$ We would like to
compute a sum of the form

$\sum_{i+j\in[T]}x_{i}y_{j}=\sum_{i}\sum_{j:j-i\in
T}x_{i}y_{j}=\sum_{i}x_{i}\sum_{j\in[1+i,T+i]}y_{j}.$

Let $S(i)$ denote the value of the inner sum for a fixed $i$. Note that
$S(i+1)=S(i)+y_{T+j+1}-y_{i+1}.$ Then $S(0)$ can be computed in $O(n)$ time
and each of $S(1),S(2),\ldots,S(n-1)$ can be computed in turn, each in $O(1)$
time using the recurrence. This gives an $O(n)$ algorithm for computing the
original sum.

m=3 Constructions. Let $T$ be arbitrary and set
$B_{T}=\\{(i,j,k):i+j+k\in[0,T-1]\\}$.consider a sum of the form

$\sum_{(i,j,k)\in
B_{T}}x_{i}y_{j}z_{k}=\sum_{t\in[0,T-1]}\sum_{i+j+k=t}x_{i}y_{j}z_{k}.$

Each of the inner sums occurs as an entry in the circular convolution $x*y*z$.
A circular convolution can be computed in $O(n\log n)$ time using the Fast
Fourier Transform, and so $B_{T}$ can be computed in $O(n\log n)$ time.

For $m=3$, we also need a sparser construction. For this, define
$C_{T}=\\{(i,j,k):i+j+k=0,j-i\in[0,T-1]\\}.$ We are interested in the sum

$\displaystyle\sum_{(i,j,k)\in C_{T}}x_{i}y_{j}z_{k}$
$\displaystyle=\sum_{i+j+k=0,j-i\in[T]}x_{i}y_{j}z_{k}$
$\displaystyle=\sum_{k}z_{-k}\sum_{i+j=k,j-i\in[T]}x_{i}y_{j}.$

As the algorithm for fast multiplication is slightly more involved we defer to
Lemma 2.5.

#### Constructing $p$-samples.

In each case we will begin by applying a random function
$P_{i}:[n]\rightarrow[n]$ independently along each mode $i$. Then we will take
the indices that land in either $A_{T}$, $B_{T}$, or $C_{T}.$

For $A_{T}$ we take the set of indices

$A_{T}^{\prime}=\\{(i,j):(P_{1}(i),P_{2}(j))\in A_{T}\\}.$

The probability that $(i,j)\in A_{T}$ is $T/n$ since $(i,j)$ is uniformly
random. For $(i,j)\neq(k,\ell)$ the values of $P_{1}(i)+P_{2}(j)$ and
$P_{1}(k)+P_{2}(\ell)$ are independent. To see this, suppose without loss of
generality that $j\neq\ell.$ Then conditioned on $P_{1}(i),P_{2}(j)$ and
$P_{1}(k)$, we have that $P_{2}(\ell)$ is uniform over $[n].$ Therefore the
events $(P_{1}(i),P_{2}(j))\in A_{T}$ and $(P_{1}(k),P_{2}(\ell))\in A_{T}$
are independent, and so

$\operatorname{Pr\/}((k,\ell)\in A_{T}^{\prime}|(i,j)\in A_{T}^{\prime})=T/n.$

For $B_{T}$, precisely the same construction and argument applies.

$C_{T}$ requires a bit more work. This time we choose $P_{1},P_{2},P_{3}$ to
be random permutations. Define $P$ by
$P((i,j,k))=(P_{1}(i),P_{2}(j),P_{3}(k)).$ As before, for any $(i,j,k)$,
$(P_{1}(i),P_{2}(j),P_{3}(k))$ is uniformly random, so
$(P_{1}(i),P_{2}(j),P_{3}(k))\in C_{T}$ with $T/n^{2}$ probability. Now
consider two pairs $u_{1}=(i_{1},j_{1},k_{1})$ and $u_{2}=(i_{2},j_{2},k_{2})$
with $u_{1}\neq u_{2}.$ Suppose that $u_{1}$ and $u_{2}$ differ in only a
single coordinate. Then since $P_{1},P_{2},P_{3}$ are permutations,

$\operatorname{Pr\/}(P(u_{1})\in C_{T}|P(U_{2})\in C_{T})=0,$

since $C_{T}$ intersects each single-mode fiber in at most one coordinate.

Next we consider the case where $u_{1}$ and $u_{2}$ differ in precisely two
coordinates. We start the case where $u_{1}$ and $u_{2}$ agree in the third
coordinate. Without loss of generality, we assume that $u_{1}=(0,0,0)$ and
$u_{2}=(1,1,0).$ Then $\operatorname{Pr\/}(P(u_{1})\in
C_{T}\,\,\text{and}\,\,P(u_{2})\in C_{T})$ is the probability that the
following events occur:

  1. 1.

$P_{1}(0)+P_{2}(0)=-P_{3}(0)$

  2. 2.

$P_{2}(0)-P_{1}(0)\in[0,T-1]$

  3. 3.

$P_{1}(1)+P_{2}(1)=-P_{3}(0)$

  4. 4.

$P_{2}(1)-P_{1}(1)\in[0,T-1]$

For any fixed value of $P_{3}(0)$, there are exactly $T$ pairs
$(P_{1}(0),P_{2}(0))$ that satisfy the first two equations, and similarly
there are $T$ pairs $(P_{1}(1),P_{2}(1))$ that satisfy the last two equations.
There are $T(T-1)$ ways to choose the two pairs, so the probability that
$P_{1}$ and $P_{2}$ give two such pairs is
$T(T-1)(\frac{1}{n}\frac{1}{n-1}\frac{1}{n}\frac{1}{n-1}).$ This implies that

$\operatorname{Pr\/}(u_{2}\in C_{T}|u_{1}\in C_{T})=\frac{T-1}{(n-1)^{2}}.$

A similar calculation gives the same probability when $u_{1}$ and $u_{2}$
agree on precisely the first or second coordinates.

Finally consider the case where $u_{1}$ and $u_{2}$ differ on all three
coordinates. Then $P_{1}(u_{1})$ is uniform over all possible triples, and
conditioned on $P_{1}(u_{1})$, $P_{2}(u_{2})$ is uniform over all triples that
differ from $P_{1}(u_{1})$ in all coordinates. There are $(n-1)^{3}$ such
triples and at most $Tn$ of these are in $C_{T}.$ So

$\operatorname{Pr\/}(u_{2}\in C_{T}|u_{1}\in C_{T})\leq\frac{Tn}{(n-1)^{3}}.$

For two modes, we have constructed $p$-samples for $p=n/T$. This is sufficient
to give a $p$-sample for all sampling probabilities down to $1/n.$ For
sampling probabilities smaller than $1/n$, the size of the sample will be
$O(n)$ with high probability, and so it admits fast summation by direct
computation.

For three modes, our two constructions give $p$-samples down to sampling
probability $1/n^{2}.$ For smaller $p$, we can again compute the desired sum
explicitly in $O(n)$ time. ∎

We now verify that $C_{T}$ indeed admits fast summation.

######  Lemma 2.5.

Let $T$ be a positive integer, and let $x,y,z\in\mathbb{R}^{n}$. Let

$C_{T}=\\{(i,j,k):i+j+k=0\,\,\text{and}\,\,i-j\in\\{0,\ldots,T-1\\},$

where the arithmetic operations are treated mod $T$. There is an algorithm
that computes

$\sum_{(i,j,k)\in C_{T}}x_{i}y_{j}z_{k}$

in $O(n\log^{2}T)$ time. For convenience, we zero-index into each vector.

###### Proof.

We first rewrite the sum as

$\displaystyle\sum_{(i,j,k)\in C_{T}}x_{i}y_{j}z_{k}$
$\displaystyle=\sum_{k}z_{-k}\sum_{(i,j):i+j=k,i-j\in\\{0,\ldots,T-1\\}}x_{i}y_{j}$
$\displaystyle=\sum_{k}z_{-k}\sum_{i:2i\in\\{k,\ldots,k+T-1\\}}x_{i}y_{k-i}.$

Note that $\\{i\in\mathbb{Z}/T:2i\in[k,k+T-1]\\}$ is a union of two intervals
$I_{k}$ and $J_{k}$ in $\mathbb{Z}/T$ which are given by

$\displaystyle I_{k}$
$\displaystyle=\left[\left\lceil\frac{k}{2}\right\rceil,\left\lfloor\frac{k+T-1}{2}\right\rfloor\right],\,\,$
$\displaystyle J_{k}$
$\displaystyle=\left[\left\lceil\frac{k+T}{2}\right\rceil,\left\lfloor\frac{k+2T-1}{2}\right\rfloor\right].$

We now split the inner sum over $I_{k}$ and $J_{k}.$ We also split the outer
sum over even and odd values of $k$ so that the intervals shift by one with
each term. This gives four sums to compute, each of which is similar. For the
first sum, we wish to compute

$\displaystyle\sum_{k^{\prime}:0\leq 2k^{\prime}<T}z_{-2k^{\prime}}$
$\displaystyle\sum_{i\in I_{2k^{\prime}}}x_{i}y_{k-i}$
$\displaystyle=\sum_{k^{\prime}:0\leq 2k^{\prime}<T}z_{-2k^{\prime}}\sum_{i\in
k^{\prime}+I_{0}}x_{i}y_{2k^{\prime}-i}.$

We will show how to compute all of the inner sums quickly. Define the
“stride-$2$” circulant matrix $Y\in\mathbb{R}^{n\times n}$ by
$Y_{kj}=y_{2k-j}.$ Define $Y^{\prime}_{kj}$ by $Y^{\prime}_{kj}=Y_{kj}$ for
$j\in k+I_{0}$ and $0$ otherwise. The sums we wish to evaluate are precisely
the entries of $Y^{\prime}x$ by construction.

Now we show how to evaluate $Y^{\prime}x.$ By increasing the dimension of
$Y^{\prime}$ and $x$ by at most $T$ (and continuing the circulant pattern), we
may assume without loss of generality that $n$ is a multiple of $|I_{0}|.$

The support of $Y^{\prime}$ contains a “stripe” of width $\left|I_{0}\right|$
running parallel to the diagonal. This stripe decomposes into a union of
$2n/T$ right triangles with legs of length $|I_{0}|,$ and with $n/T$ in each
of two orientations. More formally, by permuting rows, we may assume that
$I_{0}=[0,a)$ where $a=|I_{0}|.$ Then the triangles described are the
restrictions of $Y^{\prime}$ to sets of the form
$[ra,(r+1)a]\times[sa,(s+1)a].$ Each such triangular submatrix is of the form
considered in Lemma 2.3, so admits $O(|I_{0}|\log^{2}|I_{0}|)=O(T\log^{2}T)$
matrix-vector multiplication. Since there are $2n/T$ such submatrices of
$Y^{\prime}$, matrix vector multiplication with $Y^{\prime}$ runs in
$O(n\log^{2}T)$ time. This allows the sum in question to be computed in the
same time.

Each of the other three sums can be evaluated similarly and gives the same
runtime. So the runtime stated in the lemma follows. ∎

We also collect a couple of basic facts that follow from the definition of a
$p$-sample.

######  Proposition 2.6.

Let $S$ be a $p$-sample for $x$, and let $L$ be a subset of indices.

  1. 1.

For all $i$ in $L$, we have $S\cap L=\\{i\\}$ with probability at least
$\frac{1}{2}p(1-2p|L|).$

  2. 2.

With probability at least $\frac{1}{2}|L|p(1-2p|L|)$, $|S\cap L|=1.$

###### Proof.

Consider a fixed index $i$ in $L.$ With probability at least $p/2$ we have
$i\in S.$ Conditioned on $i\in S$, we have $\operatorname{Pr\/}(j\in S)\leq
2p$ for all $j\neq i.$ Taking a union bound over all $j\in L$ gives
$\operatorname{Pr\/}(S\cap L=\\{i\\})\geq p/2(1-2p|L|).$ The events $S\cap
L=\\{i\\}$ are mutually exclusive as $i$ ranges over $L$, so the proposition
follows. ∎

##  3 Constructing an $\ell_{0}$ sampler.

Before giving the $\ell_{0}$-sampler we give two basic constructions for
$1$-sparse recovery and singleton testing using Kahtri-Rao structured
measurements. We want Kahtri-Rao measurements here, so that we can apply our
random sign flips modewise.

######  Lemma 3.1.

Let $X\in\mathbb{R}^{n^{q}}.$ Suppose that $\sigma^{(i)}_{j}$ for $i\in[N]$
and $j\in[q]$ are random sign vectors. For $N\geq c\log
n\log\frac{1}{\delta}$, the measurements
$\left\langle\sigma^{(i)}_{1}\otimes\ldots\otimes\sigma^{(i)}_{q},X\right\rangle$
are sufficient to

  1. 1.

recover $X$ with $1-\delta$ probability if $X$ is $1$-sparse111$X$ being
$1$-sparse means that $X$ has at most one nonzero entry.

  2. 2.

to distinguish between $X$ being $1$-sparse and $X$ having at least two
nonzero entries with at least $1-\delta$ probability.

###### Proof.

First note that if $X$ is $1$-sparse then all measurements will be equal in
absolute value.

Now, suppose that $X$ is $1$-sparse with $k$ in $\textnormal{supp}(X)$. Say
that the “sign pattern” for $k$ is the sequence of signs of
$\left\langle\sigma^{(i)}_{1}\otimes\ldots\otimes\sigma^{(i)}_{q},X\right\rangle$
for $i\in[N].$ We can recover $k$ as long as its sign pattern is unique among
indices, up to inverting all the signs (since an entry of $X$ could be either
positive or negative). The entries of each measurement are pairwise
independent, so the probability that $k_{1}$ and $k_{2}$ have the same sign on
a given measurement is $1/2.$ The probability that they have the same sign on
each of $N$ measurements is $1/2^{N}$, so the probability that some other
index has the same sign pattern as $k$, up to inverting signs, is at most
$2*n^{q}/2^{N}$ which is at most $\delta$ for
$N\geq\log\frac{2n^{q}}{\delta}.$ So if $X$ is $1$-sparse, then we recover $X$
with $1-\delta$ probability.

To test if $X$ is $1$-sparse we first run the $1$-sparse recovery algorithm to
obtain a $1$-sparse candidate $Y$ for $X$ (or if the recovery fails, then we
immediately return False). Now we would like to test if $X-Y$ is $0.$ To do
this, we simply take additional sign measurements in the form of those in the
lemma, and test if they all give $0.$ Clearly if $X-Y=0$ then this will occur.
Now consider the case where $X-Y\neq 0.$ Consider the natural compression
$\mathbb{R}^{\otimes q}\rightarrow\mathbb{R}^{\otimes(q-1)}$ obtained by
applying $\sigma_{q}$ along the $q$-th mode. If $X-Y$ is non-zero then it has
a 1-dimensional fiber $v$ along mode $q$ which is non-zero. It is clear that
$\left\langle v,\sigma_{q}\right\rangle$ is non-zero with probability at least
$1/2.$ Therefor applying this inductively shows that each measurement
$\left\langle\sigma^{(i)}_{1}\otimes\ldots\otimes\sigma^{(i)}_{q},X-Y\right\rangle$
is non-zero with probability at least $1/2^{q}.$ So
$O(2^{q}\log\frac{1}{\delta})$ additional measurements suffice. ∎

######  Remark 3.2.

The above algorithm shows that there exists a collection of sign measurements
that allows for $1$-sparse recovery with probability $1.$ Thus such a
collection of signs only needs to be sampled once, and it will be valid with
probability $1-\delta.$ For the equality testing step, the randomness is
essential. However, conditioned on having a no-failure $1$-sparse recovery
algorithm, we note that it yields no false negatives. That is, if $X$ is
$1$-sparse, then the tester will always accept.

Our next lemma for $\ell_{0}$ sampling shows that if we succeed in isolating a
single element of the support of $X$ at an $O(1/|\textnormal{supp}(X)|)$
sampling level, then that element is close to uniform from the support.

######  Lemma 3.3.

Let $S$ be a $p$-sample, where $p\leq\frac{1}{4|\textnormal{supp}(X)|}$. Then
for all $i\in\textnormal{supp}(X),$

$\operatorname{Pr\/}(i\in
S\big{|}|S\cap\textnormal{supp}(X)|=1)\in\left[\frac{0.25}{|\textnormal{supp}(x)|},\frac{4}{|\textnormal{supp}(X)|}\right].$

###### Proof.

See the supplementary. ∎

We are now ready to give our construction of an $\ell_{0}$ sampler.

######  Theorem 3.4.

For $q=2,3$ there there is a linear sketch of $X\in\mathbb{R}^{n^{q}}$ with
sketching dimension $m=O(\log\frac{1}{\delta}\log^{2}n(\log\log
n+\log\frac{1}{\delta}))$ space, and a sampling algorithm that succeeds with
probability $1-\delta$, which, conditioned on succeeding, outputs an index $i$
in $\textnormal{supp}(X).$ Conditioned on succeeding, the probability that the
algorithm outputs $i\in\textnormal{supp}(X)$ is in
$[\frac{0.5}{|\textnormal{supp}(X)|},\frac{2}{|\textnormal{supp}(X)|}].$ The
algorithm also returns the value of $X_{i}.$ Moveover, the entries of the
sketching matrix can be taken to be in $\\{0,+1,-1\\}.$

For $q=2$, our sketch can be applied to rank-one tensors in $O(mn)$ time, and
when $q=3$ can be applied in $O(mn\log^{2}n)$ time.

###### Proof.

Our sketch consists of $O(\log n^{q})$ $p$-samples with sampling probabilities
$p=1,\frac{1}{2},\frac{1}{4}\ldots,\frac{1}{n^{q}}$. Additionally, we take $k$
independent such samples for each $p$ which we denote by
$S^{(p)}_{1},\ldots,S^{(p)}_{k}.$

There is some $p$ in
$[\frac{1}{8|\textnormal{supp}(X)|},\frac{1}{4|\textnormal{supp}(X)|}.$ For
this value of $p$ and fixed $i$, we have that
$\left|S^{(p)}_{i}\cap\textnormal{supp}(X)\right|=1$ with probability at least
$1/8$ by Proposition 2.6. Thus, the probability that this holds for some $i$
is at least $1-(7/8)^{k},$ which is at least $1-\delta$ for
$k\geq\log\frac{1}{\delta}.$

Our $\ell_{0}$ sampling algorithm is to to iterate through the samples in
increasing order of $p$, and within each level to iterate through the
$S^{(p)}_{i}$’s in increasing order of $i.$ For each $(p,i)$ we run the
singleton tester from Lemma 3.1. If the singleton tester accepts, then we use
the $1$-sparse recovery scheme to return the appropriate index and value. As
we have seen, with probability at least $1-\delta$, the algorithm will return
for some $p\leq\frac{1}{4|\textnormal{supp}(X)|}.$ We also use the
optimization of Remark 3.2, so that we can assume the $1$-sparse tester never
yields false negatives.

As long as this occurs, the algorithm will be successful by the previous Lemma
3.3 unless the first $S_{i}^{(p)}$ for which the singleton accepts is not
$1$-sparse. Call this event $E.$ If $E$ occurs, then the singleton tester must
fail on one of $O(k\log n)$ trials, so we would like to make its failure
probability at most $O(\delta/(k\log n).$ By Lemma 3.1 this can be
accomplished using a singleton-testing sketch of size $O(\log n\log(k\log
n/\delta))$ which is $O(\log n\log\log n+\log n\log\frac{1}{\delta})$ when
$k=O(\log\frac{1}{\delta}).$

Now, we have $k$ singleton sketches for each of $\log n$ sampling
probabilities $p$, resulting in a total sketching dimension of
$O(k\log^{2}n(\log\log
n+\log\frac{1}{\delta}))=O(\log\frac{1}{\delta}\log^{2}n(\log\log
n+\log\frac{1}{\delta}))$ ∎

###  3.1 Conclusion and Future Work

We gave a new approach to constructing sampling-based sketches which can be
quickly applied to rank-one tensors with at most $3$ modes. To demonstrate the
utility of this approach we showed how to speed up sketching for $\ell_{0}$
estimation, and for constructing $\ell_{1}$ embeddings.

A number of intriguing questions remain. For example, are there other
constructions of $p$-samples that apply to higher mode tensors? Additionally,
our sketches require $O(n)$ time to construct each entry of the sketch. While
this can be faster than the $O(n^{2})$ or $O(n^{3})$ required to expand a
rank-one tensor, it seems reasonable to hope for sketches that takes closer to
$O(1)$ time per entry. This is the case for CountSketch based constructions
that work by sampling and then hashing entries of the sample into buckets.
Unfortunately in our setting it is not clear how to efficiently compose our
$p$-samples with a CountSketch.

Finally, it would interesting to extend our ideas to other problems where
sampling-based sketches have been applied.

##  4 Experiments

We evaluate the correctness of the $\ell_{0}$ samplers. Our $\ell_{0}$-sampler
is theoretically guaranteed to output a uniformly random entry of the support,
up to some constant factor. That is, the probability we output a fixed entry
of the support of a tensor $X$ is
$[\frac{c_{1}}{|\textnormal{supp}(X)|},\frac{c_{2}}{|\textnormal{supp}(X)|}]$
for some absolute constants $c_{1}$ and $c_{2}.$ In order to keep our analysis
simple, the constants $c_{1}$ and $c_{2}$ gotten from unwinding our proof are
more extreme than necessary. We remedy this by empirically showing that our
sampler is in fact much closer to uniform.

We choose our $\ell_{0}$ sampler to have $10$ buckets at each sampling level.
For an $N\times N\times N$ tensor $X$ our sampling rates begin at $1/N^{3}$
and increase in powers of $5.$ While these parameters are somewhat differnt
from what we chose theoretically, we find that they give good practical
results.

In all experiments we work with $3$-mode tensors. We consider two different
support structures for our experiments. The first support shape we consider,
we call the “disjoint rectangle” support. In this model, the support consists
of two boxes $B_{1}$ and $B_{2}$ of dimensions $(x_{1},x_{2},x_{3})$ and
$(y_{1},y_{2},y_{3}),$ such that

Recall that the first step of $p$-sample construction is to permute each mode
separately, so up to permutation there is only one configuration of rectangles
to consider. Moreover, up to this permutation symmetry, all entries of $B_{1}$
are equivalent to one another, and all entries of $B_{2}$ are equivalent to
one another. Thus to test for uniformity, it suffices to check that the number
of samples from $B_{1}$ is approximately $|B_{1}|/(|B_{1}|+|B_{2}|).$

The second support shape that we consider is a rectangle $B$ of dimensions
$(x_{1},x_{2},x_{3})$, along with an additional $x_{1}x_{2}x_{3}$ random
entries (sampled from outside the rectangle) which we call $R$, the random
component of the support. Note that we choose the size of $R$ to be the same
as the size of $B.$ This is mainly for convenience, so that our $\ell_{0}$
tester is successful if approximately half the samples lie in $R.$

In each of experiments we use the $p$-sample construction that we previously
discussed to build an $\ell_{0}$ tester. To make the experiments run faster,
we assume access to a perfect singleton tester. While such a singleton tester
does not exist in practice, it can be approximated arbitrarily well using
ideas that we previously discussed. Here our aim is to understand how well our
$p$-sample construction works as a proxy for a truly uniform sample.

All of our experiments suggest that our $\ell_{0}$ sampling procedure behaves
very nearly perfectly on the tensors described above. All experiments are
accurate to within a few percent of what one would expect for uniform
sampling. See section C in the appendix for tables showing our experimental
data.

The code for the experiments, along with additional implementations of our
sketch are available at https://github.com/wswartworth/tensorSampling.

## References

  * [AMS96] Noga Alon, Yossi Matias and Mario Szegedy  “The space complexity of approximating the frequency moments”  In _Proceedings of the twenty-eighth annual ACM symposium on Theory of computing_ , 1996, pp. 20–29
  * [KS99] Thomas Kailath and Ali H Sayed  “Fast reliable algorithms for matrices with structure”  SIAM, 1999
  * [Ind06] Piotr Indyk  “Stable distributions, pseudorandom generators, embeddings, and data stream computation”  In _Journal of the ACM (JACM)_ 53.3 ACM New York, NY, USA, 2006, pp. 307–323
  * [IM08] Piotr Indyk and Andrew McGregor  “Declaring independence via the sketching of sketches.”  In _SODA_ 8, 2008, pp. 737–745
  * [VZ12] Elad Verbin and Qin Zhang  “Rademacher-sketch: A dimensionality-reducing embedding for sum-product norms, with an application to earth-mover distance”  In _International Colloquium on Automata, Languages, and Programming_ , 2012, pp. 834–845  Springer
  * [PP13] Ninh Pham and Rasmus Pagh  “Fast and scalable polynomial kernels via explicit feature maps”  In _Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining_ , 2013, pp. 239–247
  * [CW14] Kenneth L Clarkson and David P Woodruff  “Sketching for M-estimators: A unified approach to robust regression”  In _Proceedings of the twenty-sixth annual ACM-SIAM symposium on Discrete algorithms_ , 2014, pp. 921–939  SIAM
  * [CF14] Graham Cormode and Donatella Firmani  “A unifying framework for 0-sampling algorithms”  In _Distributed and Parallel Databases_ 32 Springer, 2014, pp. 315–335
  * [Woo+14] David P Woodruff  “Sketching as a tool for numerical linear algebra”  In _Foundations and Trends® in Theoretical Computer Science_ 10.1–2 Now Publishers, Inc., 2014, pp. 1–157
  * [WZ18] David P Woodruff and Qin Zhang  “Distributed statistical estimation of matrix products with applications”  In _Proceedings of the 37th ACM SIGMOD-SIGACT-SIGAI Symposium on Principles of Database Systems_ , 2018, pp. 383–394
  * [AK19] Thomas D Ahle and Jakob BT Knudsen  “Almost optimal tensor sketch”  In _arXiv preprint arXiv:1909.01821_ , 2019
  * [MSW19] Michela Meister, Tamas Sarlos and David Woodruff  “Tight dimensionality reduction for sketching low degree polynomial kernels”  In _Advances in Neural Information Processing Systems_ 32, 2019
  * [Ahl+20] Thomas D Ahle et al.  “Oblivious sketching of high-degree polynomial kernels”  In _Proceedings of the Fourteenth Annual ACM-SIAM Symposium on Discrete Algorithms_ , 2020, pp. 141–160  SIAM
  * [LWY21] Yi Li, David Woodruff and Taisuke Yasuda  “Exponentially improved dimensionality reduction for l1: Subspace embeddings and independence testing”  In _Conference on Learning Theory_ , 2021, pp. 3111–3195  PMLR
  * [MOW21] Alexander Munteanu, Simon Omlor and David Woodruff  “Oblivious sketching for logistic regression”  In _International Conference on Machine Learning_ , 2021, pp. 7861–7871  PMLR

##  Appendix A Construction of an $\ell_{1}$ Embedding

We design a sketch for $\ell_{1}$ embeddings similar to [VZ12]. However we
show how to design the sketch in such a way so as to allow fast application to
rank-one tensors.

We describe the general structure of our sketch here. We will later choose
parameters in order to obtain the desired guarantees.

#### Description of sketch.

Our sketch $\mathbf{S}$ consists of a series of $\ell$ sampling levels
numbered $0,\ldots,\ell-1$. At each level $h$, there are two parameters:
$T_{h}$ and $p_{h}.$ $T_{h}$ is the number of buckets at level $h$, and
$p_{h}$ is the sampling probability at level $h.$ Each bucket in level $h$ is
formed by taking a random sign combination of coordinates corresponding to a
$(p_{h}/T_{h})$-sample. So if $S$ is a $(p_{h}/T_{h})$-sample, a bucket in
level $h$ takes the form

$\frac{1}{p_{h}}\sum_{i\in S}\sigma_{i}x_{i},$

where $x\in\mathbb{R}^{n}$ is the vector that we sketch, and
$\sigma\in\\{+1,-1\\}^{n}$. We will choose $\sigma$ so that $\sigma_{i}$ is a
Kronecker product of $q$ random sign vectors, each of which is $4$-wise
independent. We use the notation $B_{1}^{(h)},\ldots,B_{T_{h}}^{(h)}$ to
represent the buckets at sampling level $h.$ We will allow $p_{h}$ to be
larger than $1$, as long as $p_{h}/T_{h}\leq 1.$

####  $\ell_{1}$ sketch analysis

We first give the following technical lemma, which is given here to shorten
the proof of the following lemma.

######  Lemma A.1.

Let $x\in\mathbb{R}^{d}$ have entries in $[0,\alpha].$ Let
$X_{1},\ldots,X_{k}\sim X$ be i.i.d. with
$\operatorname{Pr\/}(X=\left|x_{j}\right|)\geq p$ for all $j.$ Then for $k\geq
c\log\frac{1}{\delta}\max(\frac{\alpha}{p\left\|x\right\|_{1}},1)$, we have
$\frac{1}{k}(X_{1}+\ldots+X_{k})\geq\frac{1}{4}p\left\|x\right\|_{1}$ with
probability at least $1-\delta.$

###### Proof.

It suffices to consider the case where $\operatorname{Pr\/}(X=x_{j})=p$ for
all $j$ and is $0$ with the remaining probability, since this majorizes the
random variable in the lemma. Note that
$\operatorname{E\/}(X_{1}+\ldots+X_{k})=pk\left\|x\right\|_{1}$. We also have

$\operatorname{E\/}((X_{1}+\ldots+X_{k})^{2})=pk\left\|x\right\|_{2}^{2}+k(k-1)p^{2}\left\|x\right\|_{1}^{2},$

So

$\text{Var}(X_{1}+\ldots+X_{k})\leq pk\left\|x\right\|_{2}^{2}\leq
pk\alpha\left\|x\right\|_{1}.$

For $k\geq\frac{16\alpha}{p\left\|x\right\|_{1}}$ we have
$\sqrt{pk\alpha\left\|x\right\|_{1}}\leq\frac{1}{4}pk\left\|x\right\|_{1}$,
and so by Chebyshev’s inequality,

$\operatorname{Pr\/}(X_{1}+\ldots+X_{k}\leq\frac{1}{2}pk\left\|x\right\|_{1})\leq\frac{1}{4}.$

If instead we have
$k\geq(24\log(1/\delta))\cdot\max(\frac{16\alpha}{p\left\|x\right\|_{1}},1)$,
then group the sum into $24\log(1/\delta)$ terms (each group of size at least
one), each of which is at least half their expectation with $3/4$ probability.
By a Chernoff bound, at least half the groups will be at least this large with
$1-\delta$ probability, which means that the entire sum will be at least $1/4$
its expectation, as desired. ∎

Our next task is to lower bound the contribution of a given sampling level.

######  Lemma A.2.

Suppose that the number of modes $q$ is at most $3$. Consider a sampling level
of our sketch as described above, with sampling probability $p$ and $T$
buckets. Suppose that $\left\|x\right\|_{1}=1$. Fix $m$ and $\gamma$ and let
$L$ be the set of coordinates of $x$ with magnitude in $(m,\gamma m]$. Then
for $(p/T)|L|\leq 1/4,$ $T\geq c\log(1/\delta)$ and $p\geq
c\log(1/\delta)\frac{\gamma m}{\left\|x_{L}\right\|_{1}}$, we have

$\left|B_{1}\right|+\ldots+\left|B_{T}\right|\geq c\left\|x_{L}\right\|_{1}$

with probability at least $1-\delta.$

###### Proof.

For convenience, set $p^{\prime}=p/T.$ Consider a fixed bucket $B_{k}$ formed
from $p^{\prime}$-sample $S_{k}$ and random signs $\sigma_{k}$. Also consider
fixed $j\in L.$ The probability that $j\in S_{k}$ is at least $p^{\prime}/2$
by definition of the sampler.

Recall that $\sigma_{k}$ is a product of $q\in\\{2,3\\}$ $4$-wise independent
sign vectors, say $\tau_{1},\tau_{2},\tau_{3}$ We claim that conditioned on
$j\in S_{k}$, we have $\left|B_{k}\right|\geq\frac{c}{p}\left|x_{j}\right|$
with constant probability. To see this, consider first sketching by $\tau_{1}$
to compress along the first mode. For the fiber $y$ of $x$ that contains $j$,
$\left|\left\langle\tau_{1},y\right\rangle\right|$ is at least
$c\left\|y\right\|_{2}\geq c\left\|x_{j}\right\|_{2}$ with constant
probability, since $\tau_{1}$ has $4$-wise independent signs 222This follows
for example from the classic analysis of the AMS sketch[AMS96]. Taking the
average of $O(1)$ such sign measurements gives a $1\pm(1/3)$ multiplicative
approximation to the $\ell_{2}$ norm of $y$ with $2/3$ probability say (by a
Chebyshev bound). Thus, each measurement individually must be
$c\left\|y\right\|_{2}$ with constant probability.. Then iteratively apply the
same argument to sketching by $\tau_{2}$ and $\tau_{3}$ along the remaining
modes.

Using the assumption that $p^{\prime}|L|\leq 1/4,$ we have that for $i\in L$,
$\operatorname{Pr\/}(|S_{k}\cap L|=\\{i\\})\geq p^{\prime}/4$ by 2.6. Applying
the previous Lemma A.1 with $p^{\prime}/8$ gives that

$\frac{1}{T}(p\left|B_{1}\right|+\ldots+p\left|B_{T}\right|)\geq\frac{1}{4}(p^{\prime}/8)\left\|x\right\|_{1}$

with probability at least $1-\delta$ when $T\geq
c\log\frac{1}{\delta}\max(\frac{\gamma
m}{(p^{\prime}/8)\left\|x_{L}\right\|_{1}},1)$. Recalling that
$p^{\prime}=p/T,$ this rearranges to the desired bound.

∎

###  A.1 No Dilation

The following proposition will help to bound the contribution of the small
coordinates to a given sampling level. This will be the key place where we
take advantage of cancellation. Unlike earlier analyses such as [VZ12], we do
not have independence, so instead we will show how to get away with a Markov
bound (earlier analyses could have applied a similar argument).

######  Proposition A.3.

Suppose that $x\in\mathbb{R}^{n}$ with $\left\|x\right\|_{1}\leq 1$ and
$\left\|x\right\|_{\infty}\leq\alpha.$ Let
$\sigma_{1},\ldots,\sigma_{n}\in\\{0,-1,1\\}$ be such that
$\operatorname{E\/}(\sigma_{i}\sigma_{j})=0$ for all $i\neq j$ and
$\operatorname{E\/}(\left|\sigma_{i}\right|)\leq q.$ Then

$\operatorname{E\/}(\left|\sigma_{1}y_{1}+\ldots+\sigma_{n}y_{n}\right|)\leq\min(\sqrt{\alpha
q},q).$

###### Proof.

By Jensen’s inequality,

$\displaystyle\operatorname{E\/}(\left|\sigma_{1}y_{1}+\ldots+\sigma_{n}y_{n}\right|)^{2}$
$\displaystyle\leq\operatorname{E\/}\left|\sigma_{1}y_{1}+\ldots+\sigma_{n}y_{n}\right|^{2}$
$\displaystyle=\operatorname{E\/}(\sigma_{1}^{2}y_{1}^{2}+\ldots+\sigma_{n}^{2}y_{n}^{2})$
$\displaystyle\leq qy_{1}^{2}+\ldots+qy_{n}^{2}$
$\displaystyle=q\left\|y\right\|_{2}^{2}.$

Note that
$\left\|y\right\|_{2}^{2}\leq\left\|y\right\|_{1}\left\|y\right\|_{\infty}\leq\alpha$,
so the first part of the bound follows. For the second half, we have

$\displaystyle\operatorname{E\/}(\left|\sigma_{1}y_{1}+\ldots+\sigma_{n}y_{n}\right|)$
$\displaystyle\leq\operatorname{E\/}(\left|\sigma_{1}\right|\left|y_{1}\right|+\ldots+\left|\sigma_{n}\right|\left|y_{n}\right|)$
$\displaystyle\leq q\left\|y\right\|_{1}\leq q.$

∎

######  Lemma A.4.

Consider applying our $\ell_{1}$-embedding sketch $\mathbf{S}$ to a vector
$x.$ Then with with failure probability at most,
$p_{1}/\beta_{1}+\ldots+p_{\ell}/\beta_{\ell}+1/10$, we have the bound

$\left\|\mathbf{S}x\right\|_{1}\leq
10\left(\sum_{i}\left\|x_{[\alpha_{i},\beta_{i}]}\right\|_{1}+\sum_{i}\min(\sqrt{T_{i}\alpha_{i}/p_{i}},1)\right)$

###### Proof.

Let $E_{i}$ be the event of sampling a coordinate larger than $\beta_{i}$ in
layer $i$. There are at most $1/\beta_{i}$ such coordinates, so the
probability that $E_{i}$ occurs is at most $p/\beta_{i}$ by a union bound.

Now let $x_{\leq\alpha_{i}}$ be $x$ where all coordinates larger in magnitude
than $\alpha_{i}$ are zeroed out. The previous lemma bounds the expected
contribution of $x_{\leq\alpha_{i}}$ to layer $i$ by
$\min(\sqrt{T\alpha_{i}/q},1)$ (after rescaling).

Now let $x_{[\alpha_{i},\beta_{i}]}$ contain all coordinates of $x$ between
$\alpha_{i}$ and $\beta_{i}$ with the remaining coordinates zeroed out. The
expected contribution of $x_{[\alpha_{i},\beta_{i}]}$ to layer $i$ is
$\left\|x_{[\alpha_{i},\beta_{i}]}\right\|_{1}$ since the layers are unbiased.

Thus with failure probability at most,
$p_{1}/\beta_{1}+\ldots+p_{\ell}/\beta_{\ell}+1/10$, we have the bound

$10\left(\sum_{i}\left\|x_{[\alpha_{i},\beta_{i}]}\right\|_{1}+\sum_{i}\min(\sqrt{T_{i}\alpha_{i}/p_{i}},1)\right)$

on the total $\ell_{1}$ norm of the sketch. ∎

As a corollary, a standard net argument, such as that given in [CW14] for
example, extends our embedding result for a single vector to a subspace.

######  Corollary A.5.

Let $A\in\mathbb{R}^{n\times d}$ be a matrix. There is an oblivious sketch
$\mathbf{S}$ with sketching dimension $m=O(d^{2}\log n+\log^{4}n)$ such that

$c_{1}\left\|Av\right\|_{1}\leq\left\|\mathbf{S}Av\right\|_{1}\leq
c_{2}\left\|Av\right\|_{1},$

where $0<c_{1}<c_{2}$ are absolute constants. Moreover all entries of
$\mathbf{S}$ can be taken to be in $\\{0,+1,-1\\}$, and $\mathbf{S}$ can be
applied to a tensor in $(\mathbb{R}^{k})^{\otimes 2}$ in $O(mk)$ time and can
be applied to a tensor in $(\mathbb{R}^{k})^{\otimes 3}$ in $O(mk\log^{2}k)$
time.

#### Combining the bounds

######  Theorem A.6.

There is an $O(1)$-distortion $\ell_{1}$ embedding sketch $\mathbf{S}$
$\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ with sketching dimension
$m=O(\log^{4}n+\log^{2}(1/\delta)\log n)$ that satisfies
$\left\|\mathbf{S}x\right\|_{1}\leq c\left\|x\right\|_{1}$ with constant
probability, and satisfies $\left\|\mathbf{S}x\right\|_{1}\geq
c\left\|x\right\|_{1}$ with probability at least $1-\delta.$

Moreover our sketch can be applied to rank one tensors in
$\mathbb{R}^{k}\otimes\mathbb{R}^{k}$ in $O(mk)$ time, and to rank one tensors
in $\mathbb{R}^{k}\otimes\mathbb{R}^{k}\otimes\mathbb{R}^{k}$ in
$O(mk\log^{2}k)$ time.

###### Proof.

For $h\geq 0$, we set $p_{h}=q^{h}$ for a fixed $q$ and set $T_{h}=T$ for a
fixed value of $T$, both to be determined. For notational convenience we will
also define a layer $-1$ of the sketch with $p_{-1}=q^{-1}$ and $T_{-1}=T.$
While $q^{-1}$ is larger than $1,$ $T$ will be chosen so that $q^{-1}/T<1$.

#### Contraction Bound.

Divide the coordinates of $x$ into levels $L_{i}$ for $i\geq 0$, with $L_{i}$
consisting of the coordinates of $x$ in $(q^{i+1},q^{i}]$ and all other
coordinates replaced with $0.$ We say that a level $L_{i}$ is heavy if
$\left\|x_{L_{i}}\right\|_{1}\geq 1/(10\ell).$

Suppose that $L_{i}$ is heavy. We will choose parameters so that level $i-1$
of the sketch, preserves a constant fraction of $L_{i}$’s mass. Setting
$p=q^{i-1}$ in Lemma A.2, the conditions we need are $q^{i-1}|L_{i}|/T\leq
1/4$, $T\geq c\log(1/\delta)$ and $q^{i-1}\geq
c\log(1/\delta)\frac{q^{i}}{\left\|x_{L_{i}}\right\|_{1}}.$ Under the
assumption that $L_{i}$ is heavy, we have
$\left\|x_{L_{i}}\right\|_{1}\geq\frac{1}{10\ell}$. Also
$\left|L_{i}\right|\leq 1/q^{i+1}$ since $\left\|x\right\|_{1}=1.$ Therefore
it suffices to choose parameters so that the following relations hold:

  1. 1.

$\frac{1}{q}\geq c(\log\frac{1}{\delta})\ell$

  2. 2.

$T\geq\frac{c}{q^{2}}$

(Note that the condition $T\geq c\log\frac{1}{\delta}$ is redundant in light
of the two listed.)

#### Dilation Bound.

Next we bound the dilation of the sketch on $x$ with $\left\|x\right\|_{1}=1.$
In the dilation bound above, we choose $\alpha_{i}=p^{i+3}$ and
$\beta_{i}=p^{i-1}.$

Then the first sum in the dilation bound Lemma A.4, is bounded by $4$ since
each coordinate of $x$ appears in the sum at most $4$ times. For the second
sum, we have

$\sum_{i=1}^{\ell}\min(\sqrt{T_{i}\alpha_{i}/p_{i}},1)\leq\sum_{i=1}^{\ell}\sqrt{Tq^{i+3}/q^{i}}=\sum_{i=1}^{\ell}\sqrt{Tq^{3}}.$

We will choose parameters so that $Tq^{3}<1/\ell^{2}$, so that this sum will
be bounded by $1.$ Then

$\sum_{i=0}^{\ell}\min(\sqrt{T_{i}\alpha_{i}/p_{i}},1)\leq 1+1=2.$

We therefore get an $O(1)$ dilation bound with failure probability at most
$\frac{1}{10}+\sum_{i=1}^{\ell}\frac{q^{i}}{q^{i-1}}\leq q\ell.$ Thus we must
choose $q\leq c/\ell.$

#### Choosing Parameters.

Taking stock of our constraints on parameters, we have

  1. 1.

$\frac{1}{q}\geq c(\log\frac{1}{\delta})\ell$

  2. 2.

$T\geq\frac{c}{q^{2}}$

  3. 3.

$Tq^{3}<1/\ell^{2}$

  4. 4.

$q\leq\frac{c}{\ell}$.

These can be satisfied by setting
$\frac{1}{q}=c\max(\log\frac{1}{\delta}\ell,\ell^{2})$ and $T=c/q^{2}.$ Since
$\ell\leq\log n$ this yields an overall bound of
$T=O(\log^{4}n+\log^{2}(1/\delta)\log^{2}n).$ Since there are $T$ buckets in
each of $\ell$ levels, this gives the desired bound on the sketching
dimension.

Finally the time bounds follow by using our construction of $p$-samples that
admit fast summation: we first apply the $\sigma_{i}$’s given in our
construction, which takes $O(n)$ time. We then compute the sum of the
resulting tensor entries over a $p$-sample that admits $\tilde{O}(n)$
summation. This takes $O(n)$ time in the two mode case, and $O(n\log^{2}n)$
time in the three mode case. ∎

##  Appendix B Proof of Lemma 2.9

Let $i\in\textnormal{supp}(X)$ By 2.6, we have

$\operatorname{Pr\/}(S\cap\textnormal{supp}(X)=\\{i\\})\geq\frac{1}{2}p(1-2p|\textnormal{supp}(X)|)\geq\frac{1}{4}p.$

From the definition of a $p$-sample we have

$\operatorname{Pr\/}(S\cap\textnormal{supp}(X)=\\{i\\})\leq\operatorname{Pr\/}(i\in
S)\leq p.$

Also from 2.6,

$\displaystyle\operatorname{Pr\/}(|S\cap\textnormal{supp}(X)|=1)$
$\displaystyle\geq\frac{1}{2}|\textnormal{supp}(X)|p(1-2p|\textnormal{supp}(X)|)$
$\displaystyle\geq\frac{1}{4}\left|\textnormal{supp}(X)\right|p.$

By taking a union bound over $i\in\textnormal{supp}(X),$ we have

$\operatorname{Pr\/}(|S\cap\textnormal{supp}(X)|=1)\leq|\textnormal{supp}(X)|p.$

It follows from this that

$\displaystyle\operatorname{Pr\/}(i\in S\big{|}|S\cap\textnormal{supp}(X)|=1)$
$\displaystyle=\frac{\operatorname{Pr\/}(S\cap\textnormal{supp}(X)=\\{i\\})}{\operatorname{Pr\/}(|S\cap\textnormal{supp}(X)|=1)}$
$\displaystyle\geq\frac{(1/4)p}{|\textnormal{supp}(X)p|}=\frac{1}{4|\textnormal{supp}(X)|}$

and

$\displaystyle\operatorname{Pr\/}(i\in S\big{|}|S\cap\textnormal{supp}(X)|=1)$
$\displaystyle\leq\frac{p}{(1/4)|\textnormal{supp}(X)|p}$
$\displaystyle=\frac{4}{|\textnormal{supp}(X)|}.$

##  Appendix C Experiment Data

All experiments are carried out on a tensor of shape $40\times 40\times 40$
for $1000$ trials. For our $\ell_{0}$ sampler, we take our levels to have
sampling probabilities increasing in powers of $5.$ We additionally use $10$
buckets per level. The following tables give the raw data for our experiments.

Rectangle dimensions | Fraction in first rectangle | Expected Fraction | Failures
---|---|---|---
((1, 1, 20), (1, 1, 1)) | 0.9650 | 0.9524 | 0
((1, 10, 20), (1, 1, 1)) | 0.9970 | 0.9950 | 0
((1, 20, 20), (1, 1, 1)) | 0.9980 | 0.9975 | 1
((20, 20, 20), (1, 1, 1)) | 1.0000 | 0.9999 | 1
((1, 1, 20), (10, 10, 10)) | 0.0200 | 0.0196 | 0
((1, 10, 20), (10, 10, 10)) | 0.1610 | 0.1667 | 0
((1, 20, 20), (10, 10, 10)) | 0.2970 | 0.2857 | 0
((20, 20, 20), (10, 10, 10)) | 0.8959 | 0.8889 | 1
((1, 1, 20), (20, 20, 20)) | 0.0010 | 0.0025 | 3
((1, 10, 20), (20, 20, 20)) | 0.0160 | 0.0244 | 2
((1, 20, 20), (20, 20, 20)) | 0.0500 | 0.0476 | 0
((20, 20, 20), (20, 20, 20)) | 0.5130 | 0.5000 | 0
Table 1: In each experiment, the support consists of two disjoint boxes. The dimensions of the two boxes are given in the leftmost column. The table shows the fraction of points that our $\ell_{0}$ sampler chooses in the first rectangle, as well as the fraction that we expect for a perfect $\ell_{0}$ sampler. The rightmost column records the number of times our $\ell_{0}$ sampler failed. Table 2: Each row of the table corresponds to an experiment run on a support shape consisting of a box of dimension specified by the left column of the table, as well as an equal number of additional uniformly sampled entries. We run $1000$ trials with our $\ell_{0}$ sampler, and record the fraction of samples from successful runs of the sampler that are in the rectangle. The rightmost column records the number of times the sampler failed. Rectangle shape | Fraction in rectangle | Failures
---|---|---
(1, 1, 1) | 0.47 | 5
(1, 1, 3) | 0.51 | 0
(1, 1, 9) | 0.52 | 2
(1, 1, 27) | 0.52 | 0
(1, 3, 1) | 0.52 | 1
(1, 3, 3) | 0.51 | 4
(1, 3, 9) | 0.50 | 0
(1, 9, 1) | 0.51 | 3
(1, 9, 3) | 0.51 | 0
(1, 9, 9) | 0.50 | 0
(3, 1, 1) | 0.49 | 1
(3, 1, 3) | 0.50 | 1
(3, 1, 9) | 0.52 | 0
(3, 3, 1) | 0.49 | 0
(3, 3, 3) | 0.53 | 0
(3, 3, 9) | 0.51 | 0
(3, 9, 1) | 0.51 | 0
(3, 9, 3) | 0.49 | 0
(3, 9, 9) | 0.51 | 0
(9, 1, 1) | 0.53 | 1
(9, 1, 3) | 0.50 | 0
(9, 1, 9) | 0.51 | 0
(9, 3, 1) | 0.51 | 0
(9, 3, 3) | 0.49 | 0
(9, 3, 9) | 0.53 | 0
(9, 9, 1) | 0.54 | 0
(9, 9, 3) | 0.50 | 1
(9, 9, 9) | 0.52 | 1
(1, 3, 27) | 0.49 | 1
(1, 9, 27) | 0.50 | 0
(1, 27, 1) | 0.51 | 0
(1, 27, 3) | 0.50 | 0
(1, 27, 9) | 0.50 | 0
(1, 27, 27) | 0.48 | 0
(3, 1, 27) | 0.52 | 1
(3, 3, 27) | 0.51 | 1
(3, 9, 27) | 0.53 | 0
(3, 27, 1) | 0.51 | 0
(3, 27, 3) | 0.50 | 0
(3, 27, 9) | 0.53 | 2
(3, 27, 27) | 0.47 | 3
(9, 1, 27) | 0.48 | 1
(9, 3, 27) | 0.53 | 0
(9, 9, 27) | 0.51 | 1
(9, 27, 1) | 0.50 | 1
(9, 27, 3) | 0.47 | 1
(9, 27, 9) | 0.50 | 1
(9, 27, 27) | 0.49 | 0
(27, 1, 1) | 0.53 | 0
(27, 1, 3) | 0.54 | 0
(27, 1, 9) | 0.52 | 2
(27, 1, 27) | 0.50 | 2
(27, 3, 1) | 0.52 | 0
(27, 3, 3) | 0.49 | 2
(27, 3, 9) | 0.50 | 0
(27, 3, 27) | 0.50 | 0
(27, 9, 1) | 0.50 | 1
(27, 9, 3) | 0.50 | 0
(27, 9, 9) | 0.51 | 0
(27, 9, 27) | 0.53 | 0
(27, 27, 1) | 0.47 | 0
(27, 27, 3) | 0.51 | 1
(27, 27, 9) | 0.48 | 0
(27, 27, 27) | 0.48 | 0


# Detector Control System and Efficiency Performance for CMS RPC at GIF++

M. Gul,11footnotetext: Corresponding author. G. Gonzalez
Blanco,22footnotetext: Corresponding author. A. Cimmino  S. Crucy  A. Fagot
A. A. O. Rios  M. Tytgat  N. Zaganidis  S. Aly  Y. Assran  A. Radi  A. Sayed
G. Singh  M. Abbrescia  G. Iaselli  M. Maggi  G. Pugliese  P. Verwilligen  W.
V. Doninck  S. Colafranceschi  A. Sharma  L. Benussi  S. Bianco  D. Piccolo
F. Primavera  V. Bhatnagar  R. Kumari  A. Mehta  J. Singh  A. Ahmad  M. I.
Asghar  S. Muhammad  I. M. Awan  H. R. Hoorani  W. Ahmed  H. Shahzad  M. A.
Shah  S. W. Cho  S. Y. Choi  B. Hong  M. H. Kang  K. S. Lee  J. H. Lim  S. K.
Park  M. S. Kim  M. Goutzvitz  G. Grenier  F. Lagarde  F. Lagarde  C. U.
Estrada  I. Pedraza  C. B. Severiano  S. Carrillo Moreno  F. Vazquez Valencia
L. M. Pant  S. Buontempo  N. Cavallo  M. Esposito  F. Fabozzi  G. Lanza  L.
Lista  S. Meola  M. Merola  I. Orso  P. Paolucci  F. Thyssen  A. Braghieri  A.
Magnani  P. Montagna  C. Riccardi  P. Salvini  I. Vai  P. Vitulo  Y. Ban  S.
J. Qian  M. Choi  Y. Choi  J. Goh  D. Kim  A. Aleksandrov  R. Hadjiiska  P.
Iaydjiev  M. Rodozov  S. Stoykova  G. Sultanov  M. Vutova  A. Dimitrov  L.
Litov  B. Pavlov  P. Petkov  D. Lomidze  I. Bagaturia  C. Avila  A. Cabrera
J. C. Sanabria  I. Crotty  and J. Vaitkus

###### Abstract

In the framework of the High Luminosity LHC upgrade program, the CMS muon
group built several different RPC prototypes that are now under test at the
new CERN Gamma Irradiation Facility (GIF++). A dedicated Detector Control
System (DCS) has been developed using the WinCC-OA tool to control and monitor
these prototype detectors and to store the measured parameters data.
Preliminary efficiency studies that set the base performance measurements of
CMS RPC for starting aging studies are also presented.

##  1 Introduction

The High Luminosity Large Hadron Collider (HL-LHC) machine will induce higher
background radiations compared to the current operating conditions. It is
important to study the performance and stability of the currently installed
and future detectors in a high radiation environment. Focused on these
requirements, the CERN Engineering- (EN) and Physics- Department (PH) made a
joint project, the Gamma Irradiation Facility (GIF++) [1]. GIF++ is the new
CERN irradiation facility located in the North Area of the CERN Super Proton
Synchrotron (SPS). It is a unique place where high energy ($\sim$100 GeV)
charged particles (mainly muons) are combined with a high flux of gamma
radiation (662 keV) produced by 13.9 TBq 137Cs source [2]. An attenuator
system is installed to vary the gamma flux on the two sides of the source
independently. A schematic overview of the GIF++ is shown in figure 1.
The Compact Muon solenoid (CMS) is one of the two general-purpose detectors at
the LHC [3], which uses Resistive Plate Chambers (RPCs) along with other
detectors [4] for the muon detection. In April 2015, the RPC detectors were
installed at the GIF++ to study the performance at a radiation dose equivalent
at 3000 fb-1 of the CMS. A dedicated control system has been built to control
these detectors and archive the relevant parameters using the WinCC-OA (PVSS)
Supervisory Control And Data Acquisition (SCADA) system [5]. The system
controls high voltage and low voltage supplies and monitors temperature,
pressure and humidity of both the RPC gas and the environment. The source
status and attenuator values are accessed through the Data Interchange
Protocol (DIP), published centrally by the Engineering Department. The RPC gas
supply is controlled and monitored by an external WinCC-OA project, that
shares relevant parameters with this project. All relevant parameters are
archived in a Structured Query Language (SQL) database (DB) for offline
analysis.
One of the features of the GIF++ RPC DCS system is accessing the source status
and attenuator values. Based on this information, RPC performance parameters
(efficiency, working point, cluster size and resistivity) are measured. To
retrieve the data from the database, a specific algorithm has been developed
to synchronize the detector parameters (current and voltage) with the external
parameters (temperature, pressure and humidity). It enables to monitor
precisely the effect of external parameters on the detector.

Figure 1: An overview of the GIF++.

##  2 Motivation

At HL-LHC, RPCs and associated electronics will operate at much higher
luminosity ($\sim$ 5 $\times$ 1034 cm-2 s-1) while accumulating large
radiation dose expected from more than 3000 fb-1. A precise understanding of
the ageing of detector’s material and currently installed electronics is
needed in intensively high radiation. The GIF++ gives the same conditions
which will be present at CMS during HL-LHC. The RPCs are irradiated with
photons, produced by a 137Cs radiation source of 13.9 TBq. At the same time, a
high-momentum-particles beam extracted from the SPS, mainly consist of muons,
is used to study the performance of the detectors by measuring and comparing
their efficiency. The two main gas components of the RPC detector are C2H2F4
and SF6 which are going to be gradually banned due to greenhouse effects. An
R&D of different gases mixture is under study at GIF++ to replace the two
gases by eco-friendly gases. When a good candidate for a new RPC gas mixture
is identified, the study will be further carried out at GIF++ to include long-
term assessment of the radiation tolerance of the chambers operated with the
new gas mixture. Further details can be found in [6].

##  3 WinCC-OA

Large experiments at CERN use WinCC-OA as a tool to develop control systems.
It allows for the description of a device in terms of a data point, with data
point elements representing its parameters. These can then be addressed
directly to write to and read from the corresponding device. Parameters of
interest can then be stored by WinCC-OA in an internal database for offline
analyses. WinCC-OA provides the facility to build a graphical user interface
(GUI).

##  4 The CMS RPC DCS project at GIF++

The CMS RPC DCS at GIF++ has been developed by using WinCC-OA 3.11 and
extended with the standard Joint Control Project (JCOP) framework [7]. The
JCOP framework provides extra functions such as standardized Finite State
Machine (FSM), the additional Graphical User Interface (GUI), the alarm
handlers and the ORACLE database interface [8]. The high voltage and low
voltage power supplies used in the RPC GIF++ setup consist of CAEN SY1527
mainframe modules as well as CAEN EASY modules, with additional ADC modules
used to read out gas and environmental sensors (pressure, temperature and
humidity). The project has access to the hardware registered through an Object
Linking and Embedding (OLE) for Process Control (OPC) server provided by CAEN
using the OPC protocol [9]. The project controls the high voltage and low
voltage system through the OPC protocol. The environmental and gas sensors (
for pressure, temperature and humidity) are also readout via the OPC protocol.
The source status and attenuator values are available centrally via the Data
Interchange Protocol (DIP). The project has been designed as a distributed one
in order to be able to communicate with other projects and to read valuable
information. Communication has been established with the central GIF++ DCS,
such that the information from the gas system, like flow rates are readable.
The Finite State machine (FSM) hierarchy of the project is based on the naming
convention of the trolley, where the detectors are installed. Each trolley has
six sections and each section accommodates one detector. Currently three CMS
RPCs trolleys are installed in the GIF++. Trolley 1 (RPC Consolidation) is
equipped with spare RPCs, trolley 2 with small glass RPCs and trolley 3 with
prototypes of improved RPCs. Detailed information of the trolleys are given
here [10]. A schematic overview of the DCS project is shown in figure 2.

Figure 2: DCS project overview.

###  4.1 High and Low Voltage System

The high and low voltage system is controlled and monitored through the CAEN
OPC server. Each gap of a chamber is independently connected to a single high
voltage channel which improves the granularity of control. The RPC front-end
electronics requires digital and analog power supplies [11]. Each low voltage
line has been shared between two front-end-boards (FEBs) for digital as well
as for analog. The low voltage boards are installed in the CAEN main frame and
controlled by DCS.

###  4.2 Environmental and Gas Parameters

The performance of RPCs strongly depends on the temperature and pressure of
the environment. Hence, it is important to measure the environmental
parameters (temperature, pressure and humidity) at different locations. The
applied voltage is corrected for the environmental temperature and pressure in
order to include its effects. This procedure is described in detail in [12].
Figure 3a gives a plot for environmental temperature, pressure and humidity.
The environmental and gas sensors (temperature, pressure and humidity) are
readable through ADC (analog-to-digital converter) board which is installed in
the EASY crate. The JCOP framework gives the opportunity to convert online the
ADC counts into physical values. The trending feature provides a comparison
among different sensors located at different positions.

###  4.3 High Voltage Scan and Stability Test

The project has been designed for R&D of detectors, hence, should be able to
perform high voltage scanning or stability tests. For high voltage scanning, a
separate branch has been incorporated in the Finite State Machine (FSM) tree,
where the user operates each detector independently. The stability test runs
for long time. Based on the requirements, a dedicated manager is used to apply
the stability script and restart it automatically. A high voltage scanning
plot for one of the CMS RPC chambers (RE3) at GIF++ is shown in figure 3b.

|
---|---
($\mathbf{a}$) | ($\mathbf{b}$)
Figure 3: (a): Environmental temperature (oC), pressure (mb) and humidity (%).
Time is on x-axis while red, blue and green lines show the values of pressure,
temperature and humidity respectively on y-axis. Pressure and temperature
values used for operating voltage correction of RPCs. (b): A high voltage scan
for CMS RE3 chamber. The x-axis shows time while the y-axis shows the voltage
(V) and current ($\mu$A) values. The red, blue and green lines are voltages
while cyan, brown and orange lines are the corresponding current values for
bottom, top narrow and top wide gaps respectively.

###  4.4 FSM and GUI

The JCOP framework provides FSM toolkits in WinCC-OA based on State Machine
Interface (SMI++). It offers an easy, robust and safe way to control the full
detector through the definition of a finite number of states, transitions and
actions (ON, OFF, STANDBY, Ramping Up, Ramping Down). All the DCS hierarchy
nodes are implemented through the FSM mechanism, shown in figure 4.
WinCC-OA provides a user friendly Graphic User Interface (GUI) panel- an
intuitive tool to control, monitor and operate the detectors in safe mode. It
provides flexibility to combine text, graphical objects and synoptic diagrams.
GUI panels can be used to see the online behavior of the detector in the form
of plots, tables and histograms.

Figure 4: FSM main tree and high voltage scan panel using GUI.

###  4.5 DataBase

To study the behavior of the detector over time and to make offline data
analysis, it is necessary to store all the important parameters in a database.
WinCC-OA’s internal SQL database is used in this project. The stored data can
be extracted for the offline analysis using a GUI.

##  5 Detector performance study

###  5.1 RPCs setup at GIF++ and Method for efficiency study

Four RPC chambers (two RE2 and two RE4 [11]) are placed parallel to each other
in vertical position. Several high voltage scans were performed using
different radiation levels (starting from the absence of radiation source), to
define the optimal operating voltage of each chamber, called working point
(WP). The efficiency (E) for different background radiation levels is
calculated using the formula:

E = NRPC/NTRACK.

The muon track is reconstructed using three RPC planes and extrapolated to the
RPC under test, looking to the closest cluster (a strip or set of continuous
strips). NTRACK corresponds to the number of muons passing through the
reference three RPC detectors at the same time. NRPC corresponds to the number
of fired clusters in the chamber under test. The dependency of the efficiency
E with respect to the effective high voltage HVeff [13] can be fitted using
the sigmoidal curve described by the subsequent formula:

$\langle$E$\rangle$ = Emax / (1 + exp (-$\lambda$ (HVeff \- HV50%))

where Emax is the maximum efficiency reached by chambers at
HV$\rightarrow\infty$, $\lambda$ is proportional to the slope of sigmoid at
flex point and HV50% is the high voltage at which a chamber reaches 50% of its
maximum efficiency. Working point of a chamber is defined by the formula HVwp
= HVknee \+ 150V, where HVknee is the voltage at which efficiency is 95% of
the maximal one.

###  5.2 Efficiency Results

In figure 5a, the efficiency as a function of HVeff for different radiation
levels is shown. The maximum efficiency decreases with the amount of radiation
received by the detector. In figure 5b, the maximum efficiency as a function
of the gamma hits rate is presented for four RPC chambers. The RPCs were
placed parallel to each other, RPC-1 being the closest to the source and RPC-4
being the most distant. The radiation doze depends on the distance between the
RPCs and the gamma source, RPCs 3 and 4 received a smaller dose as compared to
RPCs 1 and 2.

|
---|---
($\mathbf{a}$) | ($\mathbf{b}$)
Figure 5: (a): Efficiency vs HVeff for different gamma attenuator factors.
(b): Maximum efficiency vs gamma rate at HVwp for four RPCs.

##  6 Conclusion

A DCS project for CMS RPCs has successfully been implemented and tested in the
CERN GIF++. Since June, 2015 the project is running in a stable state,
operating the detectors and archiving the data. The hardware integrated in the
project, fully controls high voltage scanning and stability tests. The
environmental and gas sensors are included and used for Temperature/Pressure
corrections. Gas flow-meters are read through central DCS at GIF++ and the
data are used to study the behavior of different gases. All useful parameters
are archived in the internal database for offline analysis. As the project is
designed for detector R&D studies, any new hardware can be added easily and
safely.
The performance of the CMS RPC chambers at GIF++ has been studied and compared
at different radiation levels. At a rate of 600 Hz-cm-2, the Effmax of the
chamber was 95%.

## Acknowledgments

We wish to congratulate our colleagues in the CERN Engineering- (EN) and
Physics- Department (PH) for successful operation of the GIF++. We thank the
technical and administrative staff at CERN, other CMS institutes and RPC
group. Many thanks to ATLAS colleague Marino Romano for his technical support
to develop the project.

## References

  * [1] M. R. Jakel et al., _CERN GIF++: A new irradiation facility to test large-area particle detectors for the high-luminosity LHC program,_ PoS TIPP 2014 102 (2014).
  * [2] B. A. Gonzalez et al., _Ageing Studies on the First Resistive-MicroMeGaS Quadruplet at GIF++,_ EDP Sciences, 2015.
  * [3] The CMS Collaboration, _The CMS experiment at the CERN LHC,_ JINST 3, S08004 (2008), doi:10.1088/1748-0221/3/08/S08004.
  * [4] The CMS Collaboration, _The performance of the CMS muon detector in proton-proton collisions at $\sqrt{s}$ = 7 TeV at the LHC, _JINST 8 (2013) P11002. doi:10.1088/1748-0221/8/11/P11002.
  * [5] CERN EN-Department, _WinCC Open Architecture (OA) Service,_ wikis.web.cern.ch/wikis/display/EN/WinCC-OA+Service.
  * [6] The CMS Collaboration, _TECHNICAL PROPOSAL FOR THE PHASE-II UPGRADE OF THE COMPACT MUON SOLENOID,_ CERN-LHCC-2015-010.
  * [7] O. Holme et al., _THE JCOP FRAMEWORK,_ 10th ICALEPCS Int. Conf. on Accelerator & Large Expt. Physics Control Systems. Geneva, 10 - 14 Oct 2005 , WE2.1-6O (2005).
  * [8] G. Polese et al., _The Detector Control Systems for the CMS Resistive Plate Chamber,_ J. Phys. Conf. Ser. 219 (2010) 022019. doi:10.1088/1742-6596/219/2/022019.
  * [9] O. Passalacqua et al., _Integrating OPC Data into GSN Infrastructures,_ IADIS International Conference APPLIED COMPUTING 2008, arXiv:0808.0055v1.
  * [10] S. C. Moreno et al., _CMS RPC preliminary results for aging studies at new CERN GIF++ facility,_ proceedings from the XIII RPC worksop and related detectors, UGent Belgium, February, 2016.
  * [11] The CMS Collaboration, _The Muon Project Technical Design Report,_ CERN/LHC 97-32, CMS TDR 3, 15 December, 1997.
  * [12] M. Abbrescia et al., _Operation, performance and upgrade of the CMS Resistive Plate Chamber system at LHC,_ Nucl. Instr. and Meth. A 732, 195 (2013).
  * [13] F. Thyssen et at., _Performance of the Resistive Plate Chambers in the CMS experiment_ JINST 7 (2012) C01104. doi:10.1088/1748-0221/7/01/C01104.