# Polylogarithmic Sketches for Clustering

Moses Charikar  Erik Waingarten

###### Abstract

Given $n$ points in $\ell_{p}^{d}$, we consider the problem of partitioning
points into $k$ clusters with associated centers. The cost of a clustering is
the sum of $p^{\text{th}}$ powers of distances of points to their cluster
centers. For $p\in[1,2]$, we design sketches of size
$\mathrm{poly}(\log(nd),k,1/\epsilon)$ such that the cost of the optimal
clustering can be estimated to within factor $1+\epsilon$, despite the fact
that the compressed representation does not contain enough information to
recover the cluster centers or the partition into clusters. This leads to a
streaming algorithm for estimating the clustering cost with space
$\mathrm{poly}(\log(nd),k,1/\epsilon)$. We also obtain a distributed memory
algorithm, where the $n$ points are arbitrarily partitioned amongst $m$
machines, each of which sends information to a central party who then computes
an approximation of the clustering cost. Prior to this work, no such streaming
or distributed-memory algorithm was known with sublinear dependence on $d$ for
$p\in[1,2)$.

###### Contents

  1. 1 Introduction
    1. 1.1 Our results
    2. 1.2 Technical Overview
  2. 2 Sketching Median Costs
    1. 2.1 Statement of Main Lemma
    2. 2.2 Sketch for Optimizing a Single Coordinate
    3. 2.3 Grouping Single Coordinate Sketches
    4. 2.4 Proof of Lemma 2.1
    5. 2.5 Storing the randomness
  3. 3 Streaming $(k,p)$-Clustering Costs in $\ell_{p}$
    1. 3.1 Preliminaries: Coreset of [HV20] and $\ell_{p}^{p}$-median sketch
    2. 3.2 Proof of Theorem 6
  4. A Proof of (9)
  5. B Sketching Medoid Costs

##  1 Introduction

Given a large number of high-dimensional points, is it possible to compress
the raw representation into a very compact sketch so that we can understand
how clusterable the data is from just this highly compressed representation?
Given $n$ points in $d$ dimensions, we consider the problem of approximating
the cost of clustering them into $k$ clusters from a compressed representation
whose size is polylogarithmic in both $n$ and $d$.

For $n,d\in\mathbbm{N}$, let $P=\\{x_{1},\dots,x_{n}\\}\in\mathbbm{R}^{d}$ be
any set of points with polynomially bounded entries (i.e., all coordinates may
be represented with $O(\log(nd))$ bits). The $(k,z)$-clustering problem in
$\ell_{p}^{d}$, asks to partition $P$ into at most $k$ clusters
$C_{1},\dots,C_{k}$ so as to minimize

$\displaystyle\sum_{\ell=1}^{k}\min_{c_{\ell}\in\mathbbm{R}^{d}}\sum_{x\in
C_{\ell}}\|x-c_{\ell}\|_{p}^{z}.$ (1)

The problem is a generalization of the $k$-means and $k$-median problem; in
particular, in Euclidean space $(p=2)$, $z=2$ corresponds to $k$-means, and
$z=1$ to $k$-median.

We note that the raw representation of the dataset uses $O(nd\log(nd))$ bits,
and that any algorithm which outputs optimal cluster centers
$c_{1},\dots,c_{k}\in\mathbbm{R}^{d}$, or the optimal clustering
$C_{1},\dots,C_{k}$ must utilize $\Omega(kd)$, or $\Omega(n\log k)$ bits of
space, respectively. Hence, such algorithms cannot decrease the dependency on
both $n$ and $d$ simultaneously. However, this does not rule out an
exponential compression, from $O(nd\log(nd))$ bits to $\mathrm{polylog}(nd)$
bits (for constant $k$ and $z$), for algorithms which approximate the optimal
clustering _cost_ , which only needs $O(\log(nd))$ bits. In this work, we show
that it is possible to design sketches of size
$\mathrm{poly}(\log(nd),1/\epsilon)$ bits which $\epsilon$-approximate the
optimal clustering _cost_ , despite the fact that we do not have enough
information to compute the clusters nor the cluster centers which achieve such
cost.

Our results fit into a line of prior work on approximating the cost of
optimization problems without necessarily computing an optimal solution. These
have been investigated before for various problems and in various contexts,
including estimating minimum spanning tree [Ind04, CEF+05, CS09], minimum cost
matchings and Earth Mover’s Distance [Ind04, AIK08, ABIW09, BI14], minimum
vertex cover and maximum matching [PR07, NO08, YYI09, ORRR12, CKK20, KMNFT20,
Beh21] and model-fit [KV18, KBV20, BGLT20]. Specifically for clustering
problems, the value of the clustering cost is an important statistic; used,
for example, in the “elbow method” for determining the number of clusters
needed. We will show these sketches may be efficiently maintained on a stream
as well as for distributed-memory models, implying $\mathrm{polylog}(nd)$-bit
algorithms for these models of computation.

We start by reviewing a set of techniques in the literature to either reduce
the dependence on the data set size or the dependence on the dimension.

#### Coresets.

The coreset technique is a “dataset compression” mechanism, aiming to reduce
the dependency on $n$. From the $n$ points $P\subset\mathbbm{R}^{d}$, an
algorithm computes a much smaller (weighted) set of points
$S\subset\mathbbm{R}^{d}$, $w\colon S\to\mathbbm{R}_{\geq 0}$, such that the
cost of clustering the weighted points $S,w$ approximates that of $P$.
Following a long line of work [BHPI02, HPM04, AHPV05, Che09, LS10, FL11,
FSS13, BFL16, SW18, HV20, CASS21], the best coreset constructions for
$(k,z)$-clustering in $\ell_{p}$ achieve sizes
$\tilde{O}(k/\epsilon^{4})\cdot\min\\{1/\epsilon^{z-2},k\\}$ for a
$(1\pm\epsilon)$-approximation. The ensuing streaming and distributed-memory
algorithms maintain a coreset of the input; these algorithms find
(approximately) optimal centers $c_{1},\dots,c_{k}\in\mathbbm{R}^{d}$ and use
space complexity
$d\cdot\tilde{O}(k/\epsilon^{4})\cdot\min\\{1/\epsilon^{z-2},k\\}\cdot\mathrm{polylog}(n)$.111The
$\log n$-factors arise from utilizing the “merge-and-reduce” framework for
maintaining coresets on a stream [BS80, AHPV05], and the fact the coreset
constructions are randomized.

#### Dimension Reduction and Sketching.

In addition to constructing coresets, an algorithm may seek to optimize the
dependence on $d$. There is a large body of work studying (oblivious)
dimensionality reduction and sketching, where strong compression results are
known for computing distances [AMS99, KOR00, SS02, BYJKS04, CCF04, CM05, IW05,
Ind06, JW09, KNW10, AKO11, AKR15, BBCY17]. For example, for $p\in[1,2]$ there
exists a (randomized) sketch
$\mathbf{sk}\colon\mathbbm{R}^{d}\to\mathbbm{R}^{t}$ with $t$ much smaller
than $d$ such that, for any two vectors $x,y\in\mathbbm{R}^{d}$, an algorithm
can approximate $\|x-y\|_{p}$ from $\mathbf{sk}(x)$ and $\mathbf{sk}(y)$ with
high probability. While these results are encouraging, leveraging such
sketches for distance computation in order to compress entire optimization
problems like (1) is highly nontrivial. The challenge is that (1) implicitly
considers distances among infinitely many vectors, and we need to rule out the
possibility of spurious low cost solutions in the “sketched” space which do
not have an analog in the original space. In particular, prior to this work,
no streaming or distributed-memory algorithm was known which could reduce the
dependence on $d$ for $p\in[1,2)$.

There is one setting, of Euclidean space ($p=2$), where one can sketch vectors
while preserving (1). A sequence of works [BZD10, CEM+15, BBCA+19, MMR19] show
that applying a Johnson-Lindenstrauss [JL84] map
$\boldsymbol{\Pi}\colon\mathbbm{R}^{d}\to\mathbbm{R}^{t}$ with
$t=O(z^{4}\log(k/\epsilon)/\epsilon^{2})$, sketches Euclidean vectors to
$O(t\log(nd))$ bits and preserves (1) up to $1\pm\epsilon$. We emphasize that
Euclidean space $p=2$ is special in this regard, because the Johnson-
Lindenstrauss map achieves _dimension reduction_ , a property known not to
hold in $\ell_{1}$ [BC05, LN04, ACNN11]. In particular, $d$-dimensional
vectors $x\in\mathbbm{R}^{d}$ in Euclidean space are sketched to vectors
$\boldsymbol{\Pi}(x)\in\mathbbm{R}^{t}$ in Euclidean space, i.e., one
estimates $\|x-y\|_{2}$ by $\|\boldsymbol{\Pi}(x)-\boldsymbol{\Pi}(y)\|_{2}$.
Thus the optimization problem (1) for $d$-dimensional Euclidean space reduces
to the same optimization problem for a much smaller dimensional Euclidean
space. This can therefore be composed with known coreset constructions.
Importantly, the “sketched” space inherits all geometric properties of
Euclidean spaces, a key aspect of prior works, and the reason they do not
extend beyond Euclidean space. The technical challenge in applying sketches
for $\ell_{p}$ when $p\neq 2$ is that the “sketched” space is non-
geometric.222For example, the sketched space for $\ell_{p}$ with $p\neq 2$
does not satisfy the triangle inequality: it is not the case that for any
$x,y,z\in\mathbbm{R}^{d}$, the estimate of $(\mathbf{sk}(x),\mathbf{sk}(y))$
plus the estimate of $(\mathbf{sk}(y),\mathbf{sk}(z))$ is less than the
estimate of $(\mathbf{sk}(x),\mathbf{sk}(z))$. On the other hand, for $p=2$,
the estimates of $(\mathbf{sk}(x),\mathbf{sk}(y))$,
$(\mathbf{sk}(y),\mathbf{sk}(z))$, and $(\mathbf{sk}(x),\mathbf{sk}(z))$ are
$\|\mathbf{sk}(x)-\mathbf{sk}(y)\|_{2}$,
$\|\mathbf{sk}(y)-\mathbf{sk}(z)\|_{2}$, and
$\|\mathbf{sk}(x)-\mathbf{sk}(z)\|_{2}$, so the triangle inequality does hold
in the sketched space.

###  1.1 Our results

We give a streaming and distributed-memory algorithm for $(k,p)$-clustering in
$\ell_{p}$ with space complexity $\mathrm{poly}(\log(nd),k,1/\epsilon)$ bits.

######  Theorem 1 (Streaming $(k,p)$-Clustering in $\ell_{p}$).

For $p\in[1,2]$, there exists an insertion-only streaming algorithm which
processes a set of $n$ points $x_{1},\dots,x_{n}\in\mathbbm{R}^{d}$ utilizing
$\mathrm{poly}(\log(nd),k,1/\epsilon)$ bits which outputs a parameter
$\boldsymbol{\eta}\in\mathbbm{R}$ satisfying

$(1-\epsilon)\min_{\begin{subarray}{c}C_{1},\dots,C_{k}\\\ \text{partition
$[n]$}\end{subarray}}\sum_{\ell=1}^{k}\min_{c_{\ell}\in\mathbbm{R}^{d}}\sum_{i\in
C_{\ell}}\|x_{i}-c_{\ell}\|_{p}^{p}\leq\boldsymbol{\eta}\leq(1+\epsilon)\min_{\begin{subarray}{c}C_{1},\dots,C_{k}\\\
\text{partition
$[n]$}\end{subarray}}\sum_{\ell=1}^{k}\min_{c_{\ell}\in\mathbbm{R}^{d}}\sum_{i\in
C_{\ell}}\|x_{i}-c_{\ell}\|_{p}^{p}$

with probability at least $0.9$.

######  Theorem 2 (Distributed-Memory $(k,p)$-Clustering in $\ell_{p}$).

For $p\in[1,2]$, there exists a public-coin protocol where $m$ machines
receive an arbitrary partition of $n$ points
$x_{1},\dots,x_{n}\in\mathbbm{R}^{d}$, each communicates
$\mathrm{poly}(\log(md),k,1/\epsilon)$ bits to a central authority who outputs
a parameter $\boldsymbol{\eta}\in\mathbbm{R}$ satisfying

$(1-\epsilon)\min_{\begin{subarray}{c}C_{1},\dots,C_{k}\\\ \text{partition
$[n]$}\end{subarray}}\sum_{\ell=1}^{k}\min_{c_{\ell}\in\mathbbm{R}^{d}}\sum_{i\in
C_{\ell}}\|x_{i}-c_{\ell}\|_{p}^{p}\leq\boldsymbol{\eta}\leq(1+\epsilon)\min_{\begin{subarray}{c}C_{1},\dots,C_{k}\\\
\text{partition
$[n]$}\end{subarray}}\sum_{\ell=1}^{k}\min_{c_{\ell}\in\mathbbm{R}^{d}}\sum_{i\in
C_{\ell}}\|x_{i}-c_{\ell}\|_{p}^{p}$

with probability at least $0.9$.

Both algorithms will follow from applying a coreset and compressing the
representation of the coreset points into sketches to recover single-cluster
cost. Specifically, the bottleneck for our algorithm will be estimating the
cost of $(k,p)$-clustering in $\ell_{p}$ for $k=1$. We give a linear sketch
such that given a set of points $x_{1},\dots,x_{n}\in\mathbbm{R}^{d}$, one may
approximate the $\ell_{p}^{p}$-median cost:

$\min_{y\in\mathbbm{R}^{d}}\sum_{i=1}^{n}\|x_{i}-y\|_{p}^{p}.$

Most of the technical work will be devoted to sketching this
“$\ell_{p}^{p}$-median cost” objective. Then, the streaming and distributed-
memory algorithm will evaluate the sum of $\ell_{p}^{p}$-median costs for all
possible partitions of the coreset points into $k$ parts. The following
theorem gives a linear sketch for approximating the $\ell_{p}^{p}$-median
cost.333A related although different work is that of approximating the
$\ell_{p}^{p}$-median (for instance, see Appendix F of [BIRW16]). An
$\ell_{p}^{p}$-median is a vector in $\mathbbm{R}^{d}$ which means the sketch
outputs $d$ numbers; however, we will sketch the $\ell_{p}^{p}$-median cost,
which is a real number. Hence, our sketch will use
$\mathrm{poly}(\log(nd),1/\epsilon)$ space, as opposed to $\Omega(d)$ space
needed to describe an $\ell_{p}^{p}$-median.

######  Theorem 3 ($\ell_{p}^{p}$-Median Sketch).

For $p\in[1,2]$, there exists a linear sketch which processes a set of $n$
points $x_{1},\dots,x_{n}\in\mathbbm{R}^{d}$ into a vector $\mathbbm{R}^{t}$
with $t=\mathrm{poly}(\log(nd),1/\epsilon)$ and outputs a parameter
$\boldsymbol{\eta}\in\mathbbm{R}$ satsifying

$(1-\epsilon)\min_{y\in\mathbbm{R}^{d}}\sum_{i=1}^{n}\|x_{i}-y\|_{p}^{p}\leq\boldsymbol{\eta}\leq(1+\epsilon)\min_{y\in\mathbbm{R}^{d}}\sum_{i=1}^{n}\|x_{i}-y\|_{p}^{p}$

with probability at least $0.9$.

There are a few important remarks to make:

  * •

The requirement that $p\leq 2$ is necessary for the exponential compression we
desire. For $p>2$, there are strong lower bounds for sketching distances which
show that such sketches require $\Omega(d^{1-2/p})$ space [BYJKS04]. For
$p<1$, we are not aware of small coresets.

  * •

The focus of this work is on optimizing the space complexity of the sketch,
and while we do not explicitly specify the running time of the sketching and
streaming algorithms, a naive implementation runs in time
$(k\log(nd)/\epsilon)^{(k\log n/\epsilon)^{O(1)}}$. The exponential factor is
due to the fact that we evaluate the cost of all possible partitions of the
$(k\log(n)/\epsilon)^{O(1)}$-coreset points into $k$ clusters. One could
alleviate the exponential dependence to $(k/\epsilon)^{O(1)}$ (as opposed to
$(k\log n/\epsilon)^{O(1)}$) by running more sophisticated approximation
algorithms [BHPI02, KSS04] on the sketched representation of the
coreset.444The one subtlety is that the algorithm should be implemented
without explicitly considering the $d$-dimensional representation of the
points. Instead, it should only use the sketches of Theorem 3. We note that a
super-polynomial dependence on $k$ should is unavoidable, because
$(1\pm\epsilon)$-approximations for $(k,z)$-clustering problems, for _non-
constant_ $k$, are NP-hard [ACKS15, LSW17, CAK19].

  * •

It would be interesting to generalize Theorem 1 to dynamic streams. The reason
our algorithm works in the insertion-only model is that we utilize the coreset
of [HV20] with the merge-and-reduce framework [BS80, AHPV05] which do not
support deletions. While there exist dynamic coreset constructions for the
streaming model [BFL+17, HSYZ19], our use of coresets is not entirely black-
box. Other dynamic coresets, like [HK20], focus on update time and do not
optimize the space complexity. We must ensure that the algorithm for
constructing the coreset does not utilize the $d$-dimensional representation
of the dataset points. The coreset construction of [HV20] only consider
distances between the dataset points, so it suffices for us to only maintain a
sketch of the dataset points.

  * •

The fact that $z=p$ in our theorems above is a consequence of our techniques.
It is unclear to us whether this assumption is necessary, although our
approach hinges on the fact $\ell_{p}^{p}$ is additive over the $d$
coordinates. We leave this as a problem for future work.

A similar, yet importantly different notion of $(k,z)$-clustering considers
_medoid_ cost, where the centers of the $k$ clusters $c_{1},\dots,c_{k}$ are
restricted to be dataset points. While seemingly similar to the
$(k,z)$-clustering objective where centers are unrestricted, these two are
qualitatively very different from a sketching perspective. In Appendix B, we
show that while a two-pass sketching algorithm may $\epsilon$-approximate the
medoid cost, $\epsilon$-approximations for one-pass sketching algorithm
require polynomial space.

###  1.2 Technical Overview

We give an overview of Theorem 3. Once that is established, combining the
$\ell_{p}^{p}$-median sketch with coresets, thereby establishing Theorems 1
and 2 is (relatively) straight-forward (for more details, see Section 3).
Recall that for $p\in[1,2]$, we will process $n$ points
$x_{1},\dots,x_{n}\in\mathbbm{R}^{d}$, and aim to return an approximation to
the $\ell_{p}^{p}$-median cost:

$\displaystyle\min_{c\in\mathbbm{R}^{d}}\sum_{i=1}^{n}\|x_{i}-c\|_{p}^{p}.$
(2)

It will be useful to assume that the points are centered, i.e.,
$\sum_{i=1}^{n}x_{i}=0\in\mathbbm{R}^{d}$ (we can enforce this because our
sketches will be linear). The approach will come from the fact that the above
optimization problem decomposes into a sum of $d$ independent optimizations,
one for each coordinate, and (2) seeks to evaluate the sum. Specifically, we
may write

$\displaystyle\min_{c\in\mathbbm{R}^{d}}\sum_{i=1}^{n}\|x_{i}-c\|_{p}^{p}=\sum_{j=1}^{d}\min_{c_{j}\in\mathbbm{R}}\sum_{i=1}^{n}|x_{ij}-c_{j}|^{p}.$

Furthermore, for any fixed $j\in[d]$, estimating

$\displaystyle\min_{c_{j}\in\mathbbm{R}}\sum_{i=1}^{n}|x_{ij}-c_{j}|^{p}$ (3)

is much more amenable to $\ell_{p}$-sketching. Specifically, we let
$x_{\cdot,j}\in\mathbbm{R}^{n}$ be the vector containing the $j$-th
coordinates of all $n$ points, and ${\bf 1}\in\mathbbm{R}^{n}$ be the all-1’s
vector. Then, the quantity
$\sum_{i=1}^{n}|x_{ij}-c_{j}|^{p}=\|x_{\cdot,j}-c_{j}{\bf 1}\|_{p}^{p}$, and
since the $\ell_{p}$-sketches are linear, an algorithm may maintain
$\mathbf{sk}(x_{\cdot,j})\in\mathbbm{R}^{t}$ (for $t=\mathrm{poly}(\log(nd))$)
and after processing, could iterate through various values of
$c_{j}\in\mathbbm{R}$ to evaluate

$\left(\sum_{i=1}^{n}|x_{ij}-c_{j}|^{p}\right)^{1/p}=\|x_{\cdot,j}-c_{j}{\bf
1}\|_{p}\approx_{1\pm\epsilon}\text{ estimate of
}(\mathbf{sk}(x_{\cdot,j}),\mathbf{sk}(c_{j}{\bf 1})),$

and output the smallest value of $c_{j}\in\mathbbm{R}$ found. In order to
guarantee an $(1\pm\epsilon)$-approximation of (3), only
$\mathrm{poly}(1/\epsilon)$ values of $c_{j}$ need to be tried (since first
evaluating $c_{j}=0$ will specify the range where the optimal $c_{j}$ may
lie). A simple union bound implies that for any fixed $j\in[d]$, we can
prepare a small sketch $\mathbf{sk}(x_{\cdot,j})$ from which we can
approximate (3).

In summary, we want to estimate the sum of $d$ minimization problems. Even
though each of the $d$ problems could be solved independently with a linear
sketch, we do not want to process $d$ linear sketches (as this increases
space). In addition, we do not know which of the $d$ minimizations will
significantly affect the sum; hence, if we only (uniformly) sampled few
$\boldsymbol{j}_{1},\dots,\boldsymbol{j}_{t}\sim[d]$ and only processed $t\ll
d$ sketches along the sampled dimensions, the variance of the estimator may be
too large, making it completely useless. The technique we will use was
recently developed in [CJLW22], building on [AKO10, JW18], under the name
“$\ell_{p}$-sampling with meta-data.” In this paper, we further develop the
ideas, and apply them to sketches for clustering in a simple and modular way.
We refer the reader to Remark 4 (following this technical overview), where we
expand on the comparison to [CJLW22].

The goal is to approximate the sum of the $d$ minimization problems by
_importance sampling_ (see Chapter 9 of [Owe13]). While importance sampling is
a well-known technique, it’s use in (one-pass) linear sketching algorithms is
counter-intuitive, and we are not aware of any linear sketches which use
importance sampling in the literature, expect for this and the recent work of
[CJLW22, CJK+22]. Importance sampling will aim to estimate (2) by sampling
with respect to an alternate distribution $\mathcal{D}$. In particular, (2)
may be re-written as

$\displaystyle d\cdot\mathop{{\bf
E}\/}_{\boldsymbol{j}\sim[d]}\left[\min_{c_{\boldsymbol{j}}\in\mathbbm{R}}\sum_{i=1}^{n}|x_{i\boldsymbol{j}}-c_{\boldsymbol{j}}|^{p}\right]=d\mathop{{\bf
E}\/}_{\boldsymbol{j}\sim\mathcal{D}}\left[\mathbf{Y}_{\boldsymbol{j}}\right]\qquad\text{where}\qquad\mathbf{Y}_{\boldsymbol{j}}\stackrel{{\scriptstyle\rm
def}}{{=}}\min_{c_{\boldsymbol{j}}\in\mathbbm{R}}\sum_{i=1}^{n}|x_{i\boldsymbol{j}}-c_{\boldsymbol{j}}|^{p}\cdot\dfrac{1}{\operatorname{{\bf
Pr}}_{\mathcal{D}}[\boldsymbol{j}]},$ (4)

where $\mathcal{D}$ is a distribution chosen so the variance of the random
variable $\mathbf{Y}_{\boldsymbol{j}}$ for $\boldsymbol{j}\sim\mathcal{D}$ is
bounded. Once the variance of the random variable is bounded, only a few
samples are needed to estimate its expectation in (4). In general, the
alternate distribution $\mathcal{D}$ depends on the data in order to decrease
the variance; for instance, coordinates $j\in[d]$ whose value of (3) is higher
should be sampled more often. Hence, importance sampling inherently interacts
with the data in a two-stage process: 1) first, it samples
$\boldsymbol{j}\sim\mathcal{D}$ (where the distribution is data-dependent),
and 2) second, it evaluates $\mathbf{Y}_{\boldsymbol{j}}$ by using (3) and
$\mathop{{\bf Pr}\/}_{\mathcal{D}}[\boldsymbol{j}]$ for the value
$\boldsymbol{j}\in[d]$ specified in the first step.

In a two-pass algorithm, the two steps may be implemented sequentially. A
_sampling_ sketch, like that of [JW18], is used to sample
$\boldsymbol{j}\sim\mathcal{D}$ in the first pass. In the second pass, the
algorithm knows the value of the sampled $\boldsymbol{j}$, so it maintains a
sketch $\mathbf{sk}(x_{\cdot,\boldsymbol{j}})$ of size $t$ and a sketch
$\mathbf{sk}^{\prime}(\mathop{{\bf Pr}\/}_{\mathcal{D}}[\boldsymbol{j}])$ of
size $t^{\prime}$ (to estimate $\operatorname{{\bf
Pr}}_{\mathcal{D}}[\boldsymbol{j}]$) from which it can evaluate the random
variable $\mathbf{Y}_{\boldsymbol{j}}$. The counter-intuitive aspect is that,
in this case, we will perform both steps in one-pass:

  * •

We will use an $\ell_{p}$-sampling sketch of [JW18] to sample from an
importance sampling distribution $\mathcal{D}$, and

  * •

Concurrently, we prepare $2d$ linear sketches: $d$ sketches
$\mathbf{sk}(x_{\cdot,j})$ to evaluate (3), one for each $j\in[d]$, and $d$
sketches $\mathbf{sk}^{\prime}(\operatorname{{\bf Pr}}_{\mathcal{D}}[j])$ to
evaluate $\operatorname{{\bf Pr}}_{\mathcal{D}}[j]$, one for each $j\in[d]$.
The non-trivial part is to _sketch the sketches_. by compressing the $2d$
linear sketches into a $O(\mathrm{polylog}(nd))$-bit Count-Min data structure
[CM05].

The guarantee will be that the $\ell_{p}$-sampling sketch of [JW18] generates
a sample $\boldsymbol{j}\sim\mathcal{D}$, and the Count-Min data structure can
recover an approximation

$\widehat{\mathbf{sk}_{1}}\approx\mathbf{sk}(x_{\cdot,\boldsymbol{j}})\qquad\text{and}\qquad\widehat{\mathbf{sk}_{2}}\approx\mathbf{sk}^{\prime}(\operatorname{{\bf
Pr}}_{\mathcal{D}}[\boldsymbol{j}]).$

Furthermore, the sketch evaluation algorithm, which executes on the
approximation $\widehat{\mathbf{sk}_{1}}$ and $\widehat{\mathbf{sk}_{2}}$,
should be able to recover $(1\pm\epsilon)$-approximations to (3) and
$\operatorname{{\bf Pr}}_{\mathcal{D}}[\boldsymbol{j}]$, so that the ratio of
the two is a $(1\pm 2\epsilon)$-approximation to
$\mathbf{Y}_{\boldsymbol{j}}$.

While the above plan provides a general recipe for importance sampling, the
idea of “sketching the sketches” may not be applied in a black-box manner.
First, the alternate distribution $\mathcal{D}$ should admit a sampling
sketch. Second, the sketch evaluation algorithm for $\mathbf{sk}(x_{\cdot,j})$
and $\mathbf{sk}^{\prime}(\operatorname{{\bf Pr}}_{\mathcal{D}}[j])$ should be
robust to the errors introduced by the Count-Min compression. Bounding the
errors introduced by the Count-Min data structure, and ensuring that the
approximate sketches $\widehat{\mathbf{sk}_{1}}$ and
$\widehat{\mathbf{sk}_{2}}$ constitutes the bulk of the technical work.
Specifically for us, the plan is executed as follows: when
$\sum_{i=1}^{n}x_{i}=0\in\mathbbm{R}^{d}$, every $j$ satisfies (see Appendix
A)

$\displaystyle\dfrac{\min_{c_{j}\in\mathbbm{R}}\sum_{i=1}^{n}|x_{ij}-c_{j}|^{p}}{\|x_{\cdot,j}\|_{p}^{p}}\in[2^{-p},1].$
(5)

Hence, we will let $\mathcal{D}$ be the distribution supported on $[d]$ given
by setting, for each $j\in[d]$,

$\displaystyle\mathop{{\bf
Pr}\/}_{\boldsymbol{j}\sim\mathcal{D}}\left[\boldsymbol{j}=j\right]$
$\displaystyle=\frac{\|x_{\cdot,j}\|_{p}^{p}}{Z}\qquad\text{where}\qquad
Z=\sum_{j=1}^{d}\|x_{\cdot,j}\|_{p}^{p}=\sum_{i=1}^{n}\sum_{j=1}^{d}|x_{ij}|^{p}.$

Note that (5) implies the variance of $\mathbf{Y}_{\boldsymbol{j}}$ for
$\boldsymbol{j}\sim\mathcal{D}$ is appropriately bounded. Furthermore, since
$\mathcal{D}$ is an $\ell_{p}$-sampling distribution, the $\ell_{p}$-sampling
sketches of [JW18] are useful for sampling $\boldsymbol{j}\sim\mathcal{D}$.
Finally, the approach of [JW18] is particularly suited for bounding the errors
incurred by Count-Min on $\widehat{\mathbf{sk}_{1}}$ and
$\widehat{\mathbf{sk}_{2}}$, which we overview below.

At a high level, the $\ell_{p}$-sampling sketch of [JW18] generates a sample
$\boldsymbol{j}$ from $[d]$ by identifying a _heavy hitter_ in a random
scaling of the vector specifying the sampling probabilities. In particular,
the algorithm generates
$\boldsymbol{u}_{1},\dots,\boldsymbol{u}_{d}\sim\mathrm{Exp}(1)$ and
identifies an entry $j\in[d]$ in the vector

$\left(\frac{\|x_{\cdot,1}\|_{p}}{\boldsymbol{u}_{1}^{1/p}},\frac{\|x_{\cdot,2}\|_{p}}{\boldsymbol{u}_{2}^{1/p}},\dots,\frac{\|x_{\cdot,d-1}\|_{p}}{\boldsymbol{u}_{d-1}^{1/p}},\frac{\|x_{\cdot,d}\|_{p}}{\boldsymbol{u}_{d}^{1/p}}\right)\in\mathbbm{R}^{d},$

whose value satisfies

$\displaystyle\dfrac{\|x_{\cdot,j}\|_{p}}{\boldsymbol{u}_{j}^{1/p}}\gtrsim\left(\sum_{j^{\prime}=1}^{d}\frac{\|x_{\cdot,j^{\prime}}\|_{p}^{p}}{\boldsymbol{u}_{j^{\prime}}}\right)^{1/p},$
(6)

and is the largest among those heavy hitters. For the coordinate
$\boldsymbol{j}\in[d]$ recovered by the $\ell_{p}$-sampling sketch [JW18], the
inequality (6) gives a lower bound on how large
$1/\boldsymbol{u}_{\boldsymbol{j}}^{1/p}$ will be. In particular, by applying
the same transformation to the vector of _sketches_ ,

$\displaystyle\left(\dfrac{\mathbf{sk}(x_{\cdot,1})}{\boldsymbol{u}_{1}^{1/p}},\dots,\dfrac{\mathbf{sk}(x_{\cdot,d})}{\boldsymbol{u}_{d}^{1/p}}\right)\in(\mathbbm{R}^{t})^{d}\qquad\text{and}\qquad\left(\dfrac{\mathbf{sk}^{\prime}(\operatorname{{\bf
Pr}}_{\mathcal{D}}[1])}{\boldsymbol{u}_{1}^{1/p}},\dots,\dfrac{\mathbf{sk}^{\prime}(\operatorname{{\bf
Pr}}_{\mathcal{D}}[d])}{\boldsymbol{u}_{d}^{1/p}}\right)\in(\mathbbm{R}^{t^{\prime}})^{d},$
(7)

the $t$ and $t^{\prime}$ coordinates corresponding to the sketches
$\mathbf{sk}(x_{\cdot,\boldsymbol{j}})\in\mathbbm{R}^{t}$ and
$\mathbf{sk}^{\prime}(\operatorname{{\bf
Pr}}_{\mathcal{D}}[\boldsymbol{j}])\in\mathbbm{R}^{t^{\prime}}$ will be heavy
hitters of those vectors as well. Namely, with only
$\mathrm{poly}(\log(nd),1/\epsilon)$-bits, the Count-Min data structure will
recover the entries of $\mathbf{sk}(x_{\cdot,\boldsymbol{j}})$ and
$\mathbf{sk}^{\prime}(\operatorname{{\bf Pr}}_{\mathcal{D}}[\boldsymbol{j}])$
up to a small additive error, proportional to the $\ell_{1}$-norm of (7). We
know the distribution of sketched vectors (7) (since these are simply
$\ell_{p}$-sketches [Ind06]), so we will be able to bound the additive error
and show that the sketch evaluation algorithms of $\widehat{\mathbf{sk}_{1}}$
and $\widehat{\mathbf{sk}_{2}}$ return the desired
$(1\pm\epsilon)$-approximations.

######  Remark 4 (Comparison to [CJLW22]).

The technique, “$\ell_{p}$-sampling with meta-data”, arises in [CJLW22] in the
following context. They seek a linear sketch
$\mathbf{sk}\colon\mathbbm{R}^{d}\to\mathbbm{R}^{t}$ which can process a
vector $y\in\mathbbm{R}^{d}$ and evaluate a weighted $\ell_{1}$-norm,
$\sum_{i=1}^{d}w_{i}(y)\cdot|y_{i}|$, where the weights
$w_{1}(y),\dots,w_{d}(y)\in\mathbbm{R}_{\geq 0}$ are themselves dependent on
the vector $y$. This arises as an algorithmic step in streaming algorithms for
geometric minimum spanning tree and the earth-mover’s distance. Mapping the
above formulation to our setting, we want to evaluate a weighted
$\ell_{1}$-norm as well, where the $i$-th weight corresponds to $\mathop{{\bf
Pr}\/}_{\boldsymbol{j}\sim\mathcal{D}}[\boldsymbol{j}=i]$, and the $i$-th
value seek to sum is $\mathbf{Y}_{i}$ (as in (4)). The perspective of this
technique as importance sampling (as presented in this work) is new. Indeed,
the appropriate setting of weights is only apparent once one multiplies and
divides the contribution of the $j$-th coordinate by $\|x_{\cdot,j}\|_{p}^{p}$
to define $\mathcal{D}$.

##  2 Sketching Median Costs

###  2.1 Statement of Main Lemma

######  Theorem 5.

Fix $n,d\in\mathbbm{N}$, as well as $p\in[1,2]$ and $\epsilon,\delta\in(0,1)$.
There exists a linear sketch using $\mathrm{poly}(\log
d,1/\epsilon,\log(1/\delta))$ space which processes a set of $n$ points
$x_{1},\dots,x_{n}\in\mathbbm{R}^{d}$, and outputs a parameter
$\boldsymbol{\eta}\in\mathbbm{R}$ which satisfies

$\min_{y\in\mathbbm{R}^{d}}\sum_{i=1}^{n}\|y-x_{i}\|_{p}^{p}\leq\boldsymbol{\eta}\leq(1+\epsilon)\min_{y\in\mathbbm{R}^{d}}\sum_{i=1}^{n}\|y-x_{i}\|_{p}^{p}$

with probability at least $1-\delta$.

We work with the following representation of a linear sketch. The processed
set of $n$ points in $\mathbbm{R}^{d}$ are stacked to form a vector
$x\in\mathbbm{R}^{nd}$. A linear sketch using space $s$ is a distribution
$\mathcal{M}$ supported on $s\times(nd)$ matrices. The theorem states that for
any fixed $x_{1},\dots,x_{n}\in\mathbbm{R}^{d}$, with probability $1-\delta$
over the draw of $\mathbf{S}\sim\mathcal{M}$, an algorithm with access to the
vector $\mathbf{S}x\in\mathbbm{R}^{s}$ and $\mathbf{S}$ can output
$\boldsymbol{\eta}$ satisfying the above guarantees.

Linear sketches of the above form imply efficient streaming algorithms, albeit
with some subtleties. It is useful to first assume that the streaming
algorithm can store its randomness for free (we will address this in
Subsection 2.5) so that it knows the matrix $\mathbf{S}$. In particular, since
$\mathbf{S}\in\mathbbm{R}^{s\times nd}$ acts on the vector
$x\in\mathbbm{R}^{nd}$ which vertically stacks
$x_{1},\dots,x_{n}\in\mathbbm{R}^{d}$, the columns of $\mathbf{S}$ may be
broken up into $n$ groups of size $d$, so

$\mathbf{S}=\left[\begin{array}[]{cccc}\mathbf{S}_{1}&\mathbf{S}_{2}&\dots&\mathbf{S}_{n}\end{array}\right],\qquad\text{and}\qquad\mathbf{S}x=\sum_{i=1}^{n}\mathbf{S}_{i}x_{i}.$

In the _insertion-only_ model, an algorithm would process the points one-at-a-
time, and at time-step $j$, maintain
$\sum_{i=1}^{j}\mathbf{S}_{i}x_{i}\in\mathbbm{R}^{s}$. In the _turnstile_
model of streaming, there is a subtlety in the implementation; namely, as the
algorithm receives insertions and deletions of points in $\mathbbm{R}^{d}$, it
must know which index $i\in[n]$ it is considering. The reason is that the
algorithm should know which of the sub-matrix $\mathbf{S}_{i}$ to update the
point with.

For our application of the $\ell_{p}^{p}$-median sketch to $(k,p)$-clustering
in $\ell_{p}$, we consider a weighted $\ell_{p}^{p}$-median. Namely, for
points $x_{1},\dots,x_{n}\in\mathbbm{R}^{d}$ and weights
$\lambda_{1},\dots,\lambda_{n}\in[0,1]$ with $\sum_{i=1}^{n}\lambda_{i}=1$,
the $\ell_{p}^{p}$-median cost with respect to weights
$\lambda_{1},\dots,\lambda_{n}$ is

$\min_{y\in\mathbbm{R}^{d}}\sum_{i=1}^{n}\lambda_{i}\|y-x_{i}\|_{p}^{p}.$

It is useful to first consider of $\lambda_{1}=\dots=\lambda_{n}=1/n$. For
general weights, the sketch will receive as input
$\mathbf{S}=[\mathbf{S}_{1},\dots,\mathbf{S}_{n}]\in\mathbbm{R}^{s\times(nd)}$,
the vector
$\sum_{i=1}^{n}\lambda_{i}^{1/p}\mathbf{S}_{i}x_{i}\in\mathbbm{R}^{s}$, and
the weights $\lambda_{1},\dots,\lambda_{n}$.

#### Centering Points

There is a straight-forward way to process the points so as to assume they are
centered. Specifically, the average point may be subtracted from every point
by applying a linear map, and since our sketch is linear, subtracting the
average point may be incorporated into the sketch. For weights
$\lambda_{1},\dots,\lambda_{n}\in[0,1]$ satisfying
$\sum_{i=1}^{n}\lambda_{i}=1$, we consider the linear map

$(x_{1},\dots,x_{n})\mathop{\mapsto}\left(x_{1}-\sum_{i=1}^{n}\lambda_{i}x_{i},\dots,x_{n}-\sum_{i=1}^{n}\lambda_{i}x_{i}\right)\in\mathbbm{R}^{nd}.$

Hence, we assume, without loss of generality, that the points
$x_{1},\dots,x_{n}\in\mathbbm{R}^{d}$ satisfy

$\displaystyle\sum_{i=1}^{n}\lambda_{i}x_{i}=0\in\mathbbm{R}^{d}.$ (8)

The centering is useful for deriving the following set of inequalities, which
will be useful for our sketching procedures. Suppose we denote
$y\in\mathbbm{R}^{d}$ as the point which minimizes
$\sum_{i=1}^{n}\lambda_{i}\|y-x_{i}\|_{p}^{p}$. Then, for every $j\in[d]$,

$\displaystyle\sum_{i=1}^{n}\lambda_{i}|y_{j}-x_{ij}|^{p}\leq\sum_{i=1}^{n}\lambda_{i}|x_{ij}|^{p}\leq
2^{p}\sum_{i=1}^{n}\lambda_{i}|y_{j}-x_{ij}|^{p}.$

Importantly for us, every $j\in[d]$ satisfies

$\displaystyle
2^{-p}\leq\dfrac{\min_{y_{j}\in\mathbbm{R}}\sum_{i=1}^{n}\lambda_{i}|y_{j}-x_{ij}|^{p}}{\sum_{i=1}^{n}\lambda_{i}|x_{ij}|^{p}}\leq
1.$ (9)

We let $\mathcal{D}$ be the distribution supported on $[d]$ given by setting,
for each $j\in[d]$,

$\displaystyle\mathop{{\bf
Pr}\/}_{\boldsymbol{j}\sim\mathcal{D}}\left[\boldsymbol{j}=j\right]$
$\displaystyle=\frac{1}{Z}\sum_{i=1}^{n}\lambda_{i}|x_{ij}|^{p}\qquad\text{where}\qquad
Z=\sum_{j^{\prime}=1}^{d}\sum_{i=1}^{n}\lambda_{i}|x_{ij^{\prime}}|^{p}=\sum_{i=1}^{n}\lambda_{i}\|x_{i}\|_{p}^{p}.$

Then, the quantity we want to estimate may be equivalently re-written as:

$\displaystyle\sum_{j=1}^{d}\min_{y_{j}\in\mathbbm{R}}\sum_{i=1}^{n}\lambda_{i}|y_{j}-x_{ij}|^{p}$
$\displaystyle=Z\cdot\mathop{{\bf
E}\/}_{\boldsymbol{j}\sim\mathcal{D}}\left[\frac{\min_{y_{\boldsymbol{j}}\in\mathbbm{R}}\sum_{i=1}^{n}\lambda_{i}|y_{\boldsymbol{j}}-x_{i\boldsymbol{j}}|^{p}}{\sum_{i=1}^{n}\lambda_{i}|x_{i\boldsymbol{j}}|^{p}}\right],$
(10)

where the value within the expectation is bounded between $2^{-p}$ and $1$.
Furthermore, the quantity $Z$ will be sketched with an $\ell_{p}$-sketch, and
a sample $\boldsymbol{j}\sim\mathcal{D}$ will be drawn with an
$\ell_{p}$-sampling sketch. Hence, the plan is to produce
$t=O(1/\epsilon^{2})$ samples of
$\boldsymbol{j}_{1},\dots,\boldsymbol{j}_{t}\sim\mathcal{D}$, and produce a
sketch to evaluate the numerator inside the expectation, and the denominator
inside the expectation. Taking an empirical average of the samples to estimate
the expectation, and multiplying it by the estimate of $Z$ will give the
desired estimator.

######  Lemma 2.1 (Main Lemma).

For any $n,d\in\mathbbm{N}$, $p\in[1,2]$ and $\epsilon,\delta\in(0,1)$, let
$s=\mathrm{poly}(\log d,1/\epsilon,1/\delta)$.555See (16) for the specific
polynomial bounds. There exists a distribution $\mathcal{S}$ over
$s\times(nd)$ matrices, and an algorithm such that for any $n$ vectors
$x_{1},\dots,x_{n}\in\mathbbm{R}^{d}$ and any
$\lambda_{1},\dots,\lambda_{n}\in[0,1]$ with $\sum_{i=1}^{n}\lambda_{i}=1$ and
$\sum_{i=1}^{n}\lambda_{i}x_{i}=0$, the following occurs:

  * •

We sample $\mathbf{S}=[\mathbf{S}_{1},\dots,\mathbf{S}_{n}]\sim\mathcal{S}$,
and we give the algorithm as input $\mathbf{S}$,
$\sum_{i=1}^{n}\mathbf{S}_{i}(\lambda_{i}^{1/p}x_{i})$, and
$\lambda_{1},\dots,\lambda_{n}$.

  * •

The algorithm outputs a tuple of three numbers
$(\boldsymbol{j},\boldsymbol{\alpha},\boldsymbol{\beta})\in[d]\times\mathbbm{R}_{\geq
0}\times\mathbbm{R}_{\geq 0}$. With probability at least $1-\delta$ over the
draw of $\mathbf{S}\sim\mathcal{S}$, we have the following two inequalities:

$\displaystyle(1-\epsilon)\left(\sum_{i=1}^{n}\lambda_{i}|x_{i\boldsymbol{j}}|^{p}\right)^{1/p}$
$\displaystyle\leq\boldsymbol{\alpha}\leq(1+\epsilon)\left(\sum_{i=1}^{n}\lambda_{i}|x_{i\boldsymbol{j}}|^{p}\right)^{1/p},$
$\displaystyle(1-\epsilon)\min_{z\in\mathbbm{R}}\left(\sum_{i=1}^{n}\lambda_{i}|x_{i\boldsymbol{j}}-z|^{p}\right)^{1/p}$
$\displaystyle\leq\boldsymbol{\beta}\leq(1+\epsilon)\min_{z\in\mathbbm{R}}\left(\sum_{i=1}^{n}\lambda_{i}|x_{i\boldsymbol{j}}-z|^{p}\right)^{1/p}.$

  * •

Furthermore, the distribution of the random variable $\boldsymbol{j}$ is
$\epsilon 2^{-p}$-close in total variation distance to $\mathcal{D}$.

Proof: [Proof of Theorem 5 assuming Lemma 2.1] Given Lemma 2.1, the proof of
Theorem 5 is straight-forward. We fix $\lambda_{1}=\dots=\lambda_{n}=1/n$, and
we first handle the centering. We will utilize Lemma 2.1 which requires
vectors $x_{1},\dots,x_{n}$ to satisfy $\sum_{i=1}^{n}\lambda_{i}x_{i}=0$;
hence, we sketch the vectors $x_{1}^{\prime},\dots,x_{n}^{\prime}$ given by
$x_{i}^{\prime}=x_{i}-\sum_{h=1}^{n}\lambda_{h}x_{h}$, which are now centered.
By linearity, this is equivalent to maintaining the vector

$\sum_{i=1}^{n}\lambda_{i}^{1/p}\mathbf{S}_{i}(x_{i}-\sum_{h=1}^{n}\lambda_{h}x_{h})=\sum_{i=1}^{n}\left(\lambda_{i}^{1/p}\mathbf{S}_{i}-\lambda_{i}\sum_{h=1}^{n}\lambda_{h}^{1/p}\mathbf{S}_{h}\right)x_{i}\in\mathbbm{R}^{s}.$

We take $t=\omega(1/\epsilon^{2})$ independent sketches from Lemma 2.1 with
accuracy parameter $\epsilon/2$ and error probability $\delta=o(1/t)$. This,
in turn, gives us $t$ independent samples
$(\boldsymbol{j}_{1},\boldsymbol{\alpha}_{1},\boldsymbol{\beta}_{1}),\dots,(\boldsymbol{j}_{t},\boldsymbol{\alpha}_{1},\boldsymbol{\beta}_{t})$.
By taking a union bound over the $t$ executions of Lemma 2.1, with high
probability, every $\boldsymbol{\alpha}_{1},\dots,\boldsymbol{\alpha}_{t}$ and
$\boldsymbol{\beta}_{1},\dots,\boldsymbol{\beta}_{t}$ satisfy

$\boldsymbol{\alpha}_{\ell}^{p}\approx_{(1+\epsilon
p/2)}\sum_{i=1}^{n}\lambda_{i}|x_{i\boldsymbol{k}}|^{p}\qquad\text{and}\qquad\boldsymbol{\beta}_{\ell}^{p}\approx_{(1+\epsilon
p/2)}\min_{z\in\mathbbm{R}}\sum_{i=1}^{n}\lambda_{i}|x_{i\boldsymbol{k}}-z|^{p},$

and $\boldsymbol{j}_{1},\dots,\boldsymbol{j}_{t}$ are independent draws from a
distribution $\mathcal{D}^{\prime}$ which is $\epsilon 2^{-p}$-close to
$\mathcal{D}$. For estimating $Z$, we use an $\ell_{p}$-sketch to accuracy
$\epsilon/2$ and failure probability $\delta=o(1)$. For example, the sketch
for $Z$ may proceed by applying an $\ell_{p}$-sketch [Ind06] to the stacked
vector $x^{\prime}\in\mathbbm{R}^{nd}$ where

$x_{ij}^{\prime}=\lambda_{i}^{1/p}\cdot x_{ij},$

so that the $\ell_{p}$ norm of $x^{\prime}$ is exactly $Z^{1/p}$. Let
$\widehat{\mathbf{Z}}$ be the estimate for the $Z$. For our estimate
$\boldsymbol{\eta}$ that we will output, we set

$\boldsymbol{\eta}=\widehat{\mathbf{Z}}\cdot\frac{1}{t}\sum_{\ell=1}^{t}\mathrm{minmax}\left\\{2^{-p},\left(\frac{\boldsymbol{\beta}_{\ell}}{\boldsymbol{\alpha}_{\ell}}\right)^{p},1\right\\},$

where $\mathrm{minmax}(l,x,u)$ is $l$ if $x\leq l$, $u$ if $u\geq x$, and $x$
otherwise. To see why our estimator approximates (10), we have
$\widehat{\mathbf{Z}}$ is a $(1\pm\epsilon/2)$-approximation of $Z$. The
latter quantity is the empirical average of $t$ i.i.d random variables, each
of which is bounded by $2^{-p}$ and $1$. In particular, we have that with
probability at least $1-o(1)$, Chebyshev’s inequality, and the conditions of
$\boldsymbol{\beta}_{\ell}$ and $\boldsymbol{\alpha}_{\ell}$,

$\mathop{{\bf
E}\/}_{\boldsymbol{j}\sim\mathcal{D}^{\prime}}\left[\dfrac{\min_{z\in\mathbbm{R}}\sum_{i=1}^{n}\lambda_{i}|x_{i\boldsymbol{j}}-z|^{p}}{\sum_{i=1}^{n}\lambda_{i}|x_{i\boldsymbol{j}}|^{p}}\right]\approx_{(1+2\epsilon
p)}\frac{1}{t}\sum_{\ell=1}^{t}\mathrm{minmax}\left\\{2^{-p},\left(\frac{\boldsymbol{\beta}_{\ell}}{\boldsymbol{\alpha}_{\ell}}\right)^{p},1\right\\}.$

It remains to show that

$\mathop{{\bf
E}\/}_{\boldsymbol{j}\sim\mathcal{D}^{\prime}}\left[\dfrac{\min_{z\in\mathbbm{R}}\sum_{i=1}^{n}\lambda_{i}|x_{i\boldsymbol{j}}-z|^{p}}{\sum_{i=1}^{n}\lambda_{i}|x_{i\boldsymbol{j}}|^{p}}\right]\approx_{(1\pm\epsilon)}\mathop{{\bf
E}\/}_{\boldsymbol{j}\sim\mathcal{D}}\left[\dfrac{\min_{z\in\mathbbm{R}}\sum_{i=1}^{n}\lambda_{i}|x_{i\boldsymbol{j}}-z|^{p}}{\sum_{i=1}^{n}\lambda_{i}|x_{i\boldsymbol{j}}|^{p}}\right].$

This follows from two facts: (1) $\mathcal{D}^{\prime}$ and $\mathcal{D}$ are
$\epsilon 2^{-p}$ close, since the random variable is at most $1$, the
expectations are off by at most an additive $\epsilon 2^{-p}$-factor, and (2)
both quantities above are the average of random variables which are at least
$2^{-p}$, so an additive $\epsilon 2^{-p}$ error is less than a multiplicative
$(1\pm\epsilon)$-error.

The above gives an estimate which is a $1\pm\epsilon$-approximation with
probability $1-o(1)$, in order to boost the probability of success to
$1-\delta$, we simply repeat $O(\log(1/\delta))$ times and output the median
estimate.

The remainder of the section is organized as follows. We give in the (next)
Subsection 2.2, the necessary sketches for obtaining $\boldsymbol{\alpha}$ and
$\boldsymbol{\beta}$ for a fixed coordinate $j$. Then, in the following
Subsection 2.3, we show how we combine various sketches from Subsection 2.2
for different $j\in[d]$ to obtain $\boldsymbol{\alpha}$ and
$\boldsymbol{\beta}$ up to some additive error. Finally, the proof of Lemma
2.1 appears in Subsection 2.4, where we apply a randomized transformation to
the input so that the additive error from Subsection 2.3 is a multiplicative
error for the specific sampled $\boldsymbol{j}$.

###  2.2 Sketch for Optimizing a Single Coordinate

In this subsection, we give linear sketches which are useful for optimizing
over a single coordinate. Specifically, given the $n$ vectors
$x_{1},\dots,x_{n}\in\mathbbm{R}^{d}$ and $j\in[d]$, we consider the $k$-th
coordinate of the $n$ vectors $x_{1j},x_{2j},\dots,x_{nj}\in\mathbbm{R}$.
Hence, the linear sketches in this section will act on vectors in
$\mathbbm{R}^{n}$, corresponding to the $j$-th coordinates of the points, and
will give approximations to

$\sum_{i=1}^{n}\lambda_{j}|x_{ij}|^{p}\quad\text{(Corollary~{}\ref{cor:cross-
ps})}\quad\text{and}\qquad\min_{y_{j}\in\mathbbm{R}}\sum_{i=1}^{n}\lambda_{j}|y_{j}-x_{ij}|^{p}.\quad\text{(Lemma~{}\ref{lem:single-
coord-opt})}$

The lemma statements also consider an additive error term,
$\mathrm{err}\in\mathbbm{R}_{\geq 0}$, which will be necessary when combining
these sketches in Subsection 2.3; however, it may be helpful to consider
$\mathrm{err}=0$ on first reading.

######  Lemma 2.2.

For any $n\in\mathbbm{N}$, $p\in[1,2]$ and $\epsilon,\delta\in(0,1)$, let
$s=O(\log(1/\delta)/\epsilon^{2})$. There exists a distribution $\mathcal{M}$
over $s\times n$ matrices, and an algorithm such that for any
$x\in\mathbbm{R}^{n}$, $\lambda_{1},\dots,\lambda_{n}\in[0,1]$, and
$y\in\mathbbm{R}$, the following occurs:

  * •

We sample $\mathbf{S}\sim\mathcal{M}$ as well as a random vector
$\boldsymbol{\chi}=(\boldsymbol{\chi}_{1},\dots,\boldsymbol{\chi}_{s})\in\mathbbm{R}^{s}$
where each is an i.i.d $p$-stable random variable. For any
$\mathrm{err}\in\mathbbm{R}_{\geq 0}$, we give the algorithm as input
$\mathbf{S}$, $\mathbf{S}(\lambda^{1/p}\circ
x)+\mathrm{err}\cdot\boldsymbol{\chi}$, the parameters
$\lambda_{1},\dots,\lambda_{n}$, and $y$.666The notation $\lambda^{1/p}\circ
x\in\mathbbm{R}^{n}$ denotes the Hadamard product, where $(\lambda^{1/p}\circ
x)_{i}=\lambda_{i}^{1/p}\cdot x_{i}$.

  * •

The algorithm outputs a parameter
$\widehat{\boldsymbol{\eta}}\in\mathbbm{R}_{\geq 0}$, which depends on
$\mathbf{S},\mathbf{S}(\lambda^{1/p}\circ
x)+\mathrm{err}\cdot\boldsymbol{\chi}$, the parameters
$\lambda_{1},\dots,\lambda_{n}$, and $y$ which satisfies with probability at
least $1-\delta$ over $\mathbf{S}$ and $\boldsymbol{\chi}$,

$\displaystyle(1-\epsilon)\left(\sum_{i=1}^{n}\lambda_{i}|x_{i}-y|^{p}+\mathrm{err}^{p}\right)^{1/p}\leq\widehat{\boldsymbol{\eta}}\leq(1+\epsilon)\left(\sum_{i=1}^{n}\lambda_{i}|x_{i}-y|^{p}+\mathrm{err}^{p}\right)^{1/p}.$
(11)

Furthermore, for every $j\in[s]$, the random variables
$(\mathbf{S}(\lambda^{1/p}\circ x))_{j}\in\mathbbm{R}$ are distributed as
$\|\lambda^{1/p}\circ x\|_{p}\cdot\boldsymbol{\chi}_{j}$, where
$\boldsymbol{\chi}_{j}$ are independent, $p$-stable random variables.

Proof: We notice that this simply corresponds to an $\ell_{p}$-sketch of the
vector $z\in\mathbbm{R}^{n}$, which is given by letting each
$z_{i}=\lambda^{1/p}\circ(x_{i}-y)$, so that the $\ell_{p}$-sketch of [Ind06]
would accomplish this task. Since the algorithm receives $\mathbf{S}$,
$\mathbf{S}(\lambda^{1/p}\circ x)+\mathrm{err}\cdot\chi$ and
$y\in\mathbbm{R}$, the algorithm may compute $\mathbf{S}(\lambda^{1/p}\circ
y\cdot{\bf 1})$, where ${\bf 1}\in\mathbbm{R}^{n}$ is an all-1’s vector, and
evaluate the sketch

$\mathbf{S}(\lambda^{1/p}\circ
x)+\mathrm{err}\cdot\boldsymbol{\chi}-\mathbf{S}(\lambda^{1/p}\circ y\cdot{\bf
1})=\mathbf{S}(\lambda^{1/p}\circ(x-y\cdot{\bf
1}))+\mathrm{err}\cdot\boldsymbol{\chi}$

by linearity. Furthermore, note that the error simply corresponds to an
$\ell_{p}$-sketch of the vector $z^{\prime}\in\mathbbm{R}^{n+1}$ which sets
$z^{\prime}_{i}=z_{i}$ for $i\neq n+1$ and $z_{n+1}^{\prime}=\mathrm{err}$.

######  Corollary 2.3.

For any $n\in\mathbbm{N}$, $p\in[1,2]$ and $\epsilon,\delta\in(0,1)$, let
$s=O(\log(1/\delta)/\epsilon^{2})$. There exists a distribution $\mathcal{M}$
over $s\times n$ matrices, and an algorithm such that for any
$x\in\mathbbm{R}^{n}$, and any $\lambda_{1},\dots,\lambda_{n}\in[0,1]$, the
following occurs:

  * •

We sample $\mathbf{S}\sim\mathcal{M}$ and a random vector
$\boldsymbol{\chi}=(\boldsymbol{\chi}_{1},\dots,\boldsymbol{\chi}_{s})\in\mathbbm{R}^{s}$
of i.i.d $p$-stable random variables. For any
$\mathrm{err}\in\mathbbm{R}_{\geq 0}$. We give the algorithm as input
$\mathbf{S}$ and $\mathbf{S}(\lambda^{1/p}\circ
x)+\mathrm{err}\cdot\boldsymbol{\chi}$.

  * •

With probability at least $1-\delta$ over $\mathbf{S}$ and
$\boldsymbol{\chi}$, the algorithm outputs a parameter
$\widehat{{\boldsymbol{\gamma}}}\in\mathbbm{R}_{\geq 0}$, which depends on
$\mathbf{S}$, and $\mathbf{S}(\lambda^{1/p}\circ x)$ which satisfies

$\displaystyle(1-\epsilon)\left(\sum_{i=1}^{n}\lambda_{i}|x_{i}|^{p}+\mathrm{err}^{p}\right)^{1/p}\leq\widehat{{\boldsymbol{\gamma}}}\leq(1+\epsilon)\left(\sum_{i=1}^{n}\lambda_{i}|x_{i}|^{p}+\mathrm{err}^{p}\right)^{1/p}.$

Furthermore, for every $j\in[s]$, the random variables
$(\mathbf{S}(\lambda^{1/p}\circ x))_{j}\in\mathbbm{R}$ are independent and
distributed as $\|\lambda^{1/p}\circ x\|_{p}\cdot\boldsymbol{\chi}_{j}$, where
$\boldsymbol{\chi}_{j}$ are independent, $p$-stable random variables.

Proof: We apply Lemma 2.2 to the vector $x\in\mathbbm{R}^{n}$ with $y=0$.

######  Lemma 2.4.

For any $n\in\mathbbm{N}$, $p\in[1,2]$ and $\epsilon,\delta\in(0,1)$, let
$s=O(\log(1/(\epsilon\delta))/\epsilon^{2})$. There exists a distribution
$\mathcal{M}$ over $s\times n$ matrices, and an algorithm such that for any
$x\in\mathbbm{R}^{n}$ and any $\lambda_{1},\dots,\lambda_{n}\in[0,1]$ with
$\sum_{i=1}^{n}\lambda_{i}=1$, whenever $\sum_{i=1}^{n}\lambda_{i}x_{i}=0$,
the following occurs:

  * •

We sample $\mathbf{S}\sim\mathcal{M}$ and a random vector
$\boldsymbol{\chi}=(\boldsymbol{\chi}_{1},\dots,\boldsymbol{\chi}_{s})$ of
i.i.d $p$-stable random variables. For any $\mathrm{err}\in\mathbbm{R}_{\geq
0}$. We give the algorithm as input $\mathbf{S}$,
$\mathbf{S}(\lambda^{1/p}\circ x)+\mathrm{err}\cdot\boldsymbol{\chi}$, the
parameters $\lambda_{1},\dots,\lambda_{n}$, and a parameter
$\gamma\in\mathbbm{R}_{\geq 0}$ satisfying

$\displaystyle(1-\epsilon)\left(\sum_{i=1}^{n}\lambda_{i}|x_{i}|^{p}\right)^{1/p}\leq\gamma\leq(1+\epsilon)\left(\sum_{i=1}^{n}\lambda_{i}|x_{i}|^{p}\right)^{1/p}.$

  * •

The algorithm outputs a parameter
$\widehat{\boldsymbol{\beta}}\in\mathbbm{R}_{\geq 0}$ which satisfies

$\displaystyle(1-\epsilon)\min_{z\in\mathbbm{R}}\left(\sum_{i=1}^{n}\lambda_{i}|x_{i}-z|^{p}+\mathrm{err}^{p}\right)^{1/p}\leq\widehat{\boldsymbol{\beta}}\leq(1+\epsilon)\min_{z\in\mathbbm{R}}\left(\sum_{i=1}^{n}\lambda_{i}|x_{i}-z|^{p}+\mathrm{err}^{p}\right)^{1/p}.$
(12)

Furthermore, for every $j\in[s]$, the random variables
$(\mathbf{S}(\lambda^{1/p}\circ x))_{j}$ are independent and distributed as
$\|\lambda^{1/p}\circ x\|_{p}\cdot\boldsymbol{\chi}_{j}$, where
$\boldsymbol{\chi}_{j}$ are independent, $p$-stable random variables.

Proof: We will utilize the sketch from Lemma 2.2, while varying the $y$’s to
find the minimum. Specifically, let $t=16\cdot 2^{p}/\epsilon$, and let the
distribution $\mathcal{M}$ be the same as that of Lemma 2.2 instantiated with
error probability $1-t\delta$ and accuracy parameter $\epsilon/2$. We
discretize the interval $[-4\gamma,4\gamma]$ into $t$, evenly-spaced out
points $y_{1},\dots,y_{t}\subset[-4\gamma,4\gamma]$ such that
$y_{\ell+1}-y_{\ell}=8\gamma/t$. We utilize the algorithm in Lemma 2.2 to
obtain estimates
$\widehat{\boldsymbol{\eta}}_{1},\dots,\widehat{\boldsymbol{\eta}}_{t}$
satisfying (11) with $y_{1},\dots,y_{t}$, respectively. Then, we output

$\widehat{\boldsymbol{\beta}}=\min_{\ell\in[t]}\widehat{\boldsymbol{\eta}}_{\ell}.$

Since we amplified the error probability to less than $t\delta$, we may
assume, by a union bound, that all estimates
$\\{\boldsymbol{\eta}_{\ell}\\}_{\ell\in[t]}$ satisfy (11) with $y_{\ell}$
with probability at least $1-\delta$. First, for any $\ell\in[t]$,

$\displaystyle\min_{z\in\mathbbm{R}}\left(\sum_{i=1}^{n}\lambda_{i}|x_{i}-z|^{p}\right)^{1/p}\leq\left(\sum_{i=1}^{n}\lambda_{i}|x_{i}-y_{\ell}|^{p}\right)^{1/p},$

and therefore, the lower bound in (12) is implied by (11). To prove the upper
bound in (12), denote $z\in\mathbbm{R}$ as the true minimizer of
$(\sum_{i=1}^{n}\lambda_{i}|x_{i}-z|^{p})^{1/p}$. By the fact
$\sum_{i=1}^{n}\lambda_{i}=1$ and the triangle inequality, we have

$\displaystyle|z|=\left(\sum_{i=1}^{n}\lambda_{i}|z|^{p}\right)^{1/p}\leq\left(\sum_{i=1}^{n}\lambda_{i}|x_{i}-z|^{p}\right)^{1/p}+\left(\sum_{i=1}^{n}\lambda_{i}|x_{i}|^{p}\right)^{1/p}\leq
2\left(\sum_{i=1}^{n}\lambda_{i}|x_{i}|^{p}\right)^{1/p},$

so $|z|\leq 2(1+\epsilon)\gamma\leq 4\gamma$, and thus
$z\in[-4\gamma,4\gamma]$. Let $\ell\in[t]$ be such that $|y_{t}-z|\leq
4\gamma/t$. Then, again by the triangle inequality and the fact
$\sum_{i=1}^{n}\lambda_{i}x_{i}=0$,

$\displaystyle\left(\sum_{i=1}^{n}\lambda_{i}|x_{i}-y_{t}|^{p}\right)^{1/p}\leq\left(\sum_{i=1}^{n}\lambda_{i}|x_{i}-z|^{p}\right)^{1/p}+4\gamma/t\leq\left(1+4(1+\epsilon)2^{p}/t\right)\left(\sum_{i=1}^{n}\lambda_{i}|x_{i}-z|^{p}\right)^{1/p}.$

By the setting of $t$,
$(\sum_{i=1}^{n}\lambda_{i}|x_{i}-y_{t}|^{p})^{1/p}\leq(1+\epsilon/2)(\sum_{i=1}^{n}\lambda_{i}|x_{i}-z|^{p})^{1/p}$,
and by (11), we obtain the desired upper bound.

###  2.3 Grouping Single Coordinate Sketches

In this subsection, we show how to compress $d$ linear sketches (one for each
coordinate) from Subsection 2.2. In the lemma that follows, the parameter
$m\in\mathbbm{N}$ should be considered the sketch size of the sketches in
Subsection 2.2, and the linear sketch will take the $d$ sketches from
Subsection 2.2 (represented as a vector $\mathbbm{R}^{dm}$). Each of the $d$
linear sketches have each coordinate of $\mathbbm{R}^{m}$ distributed as an
i.i.d scaled $p$-stable random variable (specified by the last sentence in
Corollary 2.3 and Lemma 2.4). Thus, we write the $d$ sketches as
$\boldsymbol{\Psi}_{1}v_{1},\dots,\boldsymbol{\Psi}_{d}v_{d}\in\mathbbm{R}^{m}$,
where $v_{j}\in\mathbbm{R}$ is a scaling, and
$\boldsymbol{\Psi}_{1},\dots,\boldsymbol{\Psi}_{d}\in\mathbbm{R}^{m\times n}$
are i.i.d $p$-stable matrices.

######  Lemma 2.5 ($p$-stable Sketch Compression via Count-Min).

Let $d,m\in\mathbbm{N}$, $\epsilon,\delta\in(0,1)$, and let
$t=O(\log(d/\delta))$. There exists a distribution $\mathcal{C}$ over
$(10tm/\epsilon^{p})\times(dm)$ matrices, and an algorithm such that for any
$v\in\mathbbm{R}^{d}$, the following occurs:

  * •

We sample $\mathbf{C}\sim\mathcal{C}$ and a $(dm)\times d$ matrix
$\boldsymbol{\Psi}$, where
$\boldsymbol{\Psi}=\mathrm{diag}(\boldsymbol{\Psi}_{1},\dots,\boldsymbol{\Psi}_{d})$,
and each
$\boldsymbol{\Psi}_{j}=(\boldsymbol{\chi}_{j1},\dots,\boldsymbol{\chi}_{jm})\in\mathbbm{R}^{m}$
are independent $p$-stable random vectors.777Hence, the vector
$\boldsymbol{\Psi}v\in\mathbbm{R}^{dm}$ is given by vertically stacking $d$
vectors of the form $v_{j}\boldsymbol{\Psi}_{j}\in\mathbbm{R}^{m}$. The
algorithm receives as input $\mathbf{C}$ and $\mathbf{C}\boldsymbol{\Psi}v$.

  * •

The algorithm outputs, for each $j\in[d]$, a sequence of $t$ vectors
$\widehat{\boldsymbol{z}}_{j}^{(1)},\dots,\widehat{\boldsymbol{z}}_{j}^{(t)}\in\mathbbm{R}^{m}$
which satisfy, for each $t^{\prime}\in[t]$,

$\displaystyle\widehat{\boldsymbol{z}}_{j}^{(t^{\prime})}=v_{j}\boldsymbol{\Psi}_{j}+\boldsymbol{\mathrm{err}}_{j}^{(t^{\prime})}\cdot\boldsymbol{\chi}_{j}^{(t^{\prime})}.$
(13)

where $\boldsymbol{\chi}_{j}^{(t^{\prime})}\in\mathbbm{R}^{m}$ is a vector of
independent $p$-stable random variables, and
$\boldsymbol{\mathrm{err}}_{j}^{(t^{\prime})}\in\mathbbm{R}_{\geq 0}$ only
depends on $\mathbf{C}$. With probability at least $1-\delta$ over
$\mathbf{C}$, for every $j\in[d]$

$\displaystyle\left|\left\\{t^{\prime}\in[t]:\boldsymbol{\mathrm{err}}_{j}^{(t^{\prime})}\leq\epsilon\|v\|_{p}\right\\}\right|\geq
t/2.$ (14)

Proof: The matrix $\mathbf{C}$ is a Count-Min matrix which given a vector
$u\in\mathbbm{R}^{dm}$ given by vertically stacking $d$ vectors
$u_{1},\dots,u_{d}\in\mathbbm{R}^{m}$ repeats the following process: for each
$t^{\prime}\in[t]$, we sample a hash function
$\boldsymbol{h}_{t^{\prime}}\colon[d]\to[10/\epsilon^{p}]$, and for each
$\ell\in[10/\epsilon^{p}]$ store the vector

$\boldsymbol{b}_{t^{\prime},\ell}\stackrel{{\scriptstyle\rm
def}}{{=}}\sum_{j\in[d]}{\bf 1}\\{\boldsymbol{h}_{t^{\prime}}(j)=\ell\\}\cdot
u_{j}\in\mathbbm{R}^{m}.$

In particular, the output $\mathbf{C}u$ consists of stacking $t\cdot
10/\epsilon^{p}$ vectors
$(\boldsymbol{b}_{t^{\prime},\ell}\in\mathbbm{R}^{m}:t^{\prime}\in[t],\ell\in[10/\epsilon^{p}])$,
which gives the desired bound of $10mt/\epsilon^{p}$ on the output dimension
of $\mathbf{C}$. For each $j\in[d]$ and $t^{\prime}\in[t]$, the algorithm lets
$\boldsymbol{\ell}=\boldsymbol{h}_{t^{\prime}}(j)$ and sets

$\widehat{\boldsymbol{z}}_{j}^{(t^{\prime})}=\boldsymbol{b}_{t^{\prime},\ell}=v_{j}\boldsymbol{\Psi}_{j}+\sum_{j^{\prime}\in[d]\setminus\\{j\\}}{\bf
1}\\{\boldsymbol{h}_{t^{\prime}}(j^{\prime})=\boldsymbol{\ell}\\}\cdot
v_{j^{\prime}}\cdot\boldsymbol{\Psi}_{j^{\prime}}.$

We now apply the $p$-stability property to the right-most summand, to notice
that

$\boldsymbol{\mathrm{err}}_{j}^{(t^{\prime})}=\left(\sum_{j^{\prime}\in[d]\setminus\\{j\\}}{\bf
1}\\{\boldsymbol{h}_{t^{\prime}}(j^{\prime})=\boldsymbol{\ell}\\}\cdot
v_{j^{\prime}}^{p}\right)^{1/p},$

which only depends on $\mathbf{C}$. Furthermore, the inner most summand is at
most
$\epsilon^{p}/10\sum_{j^{\prime}\in[d]\setminus\\{j\\}}v_{j^{\prime}}^{p}$ in
expectation. By Markov’s inequality, each
$\boldsymbol{\mathrm{err}}_{j}^{(t^{\prime})}\leq\epsilon\|v\|_{p}$ with
probability at least $9/10$. Since $t=O(\log(d/\delta))$, the probability that
(14) is not satisfied for each $j\in[d]$ is at most $\delta/d$ by a Chernoff
bound, so that a union bound gives the desired guarantees.

The above lemma allows us to compress $d$ many $p$-stable sketches into
$O(\log(d/\delta)/\epsilon^{p})$ many $p$-stable sketches, albeit with some
error. Since the $p$-stable sketches that we will use (from Corollary 2.3 and
Lemma 2.4) are exactly of the form $\boldsymbol{\Psi}v$ for some vector $v$,
Lemma 2.5 will allow us to compress them. Namely, we will consider $d$
sketches from Corollary 2.3 and Lemma 2.4 and utilize Lemma 2.5; for each
$j\in[d]$, we will be able to recover $t$ noisy versions of the sketch of
Corollary 2.3 and Lemma 2.4 for coordinate $j$. Importantly, the noise is of
the form an error times a $p$-stable random variable, and these are the kinds
of errors that Corollary 2.3 and Lemma 2.4 can easily handle.

######  Lemma 2.6 ($p$-stable Sketch Recovery for Sample).

For $n,d\in\mathbbm{N}$, $p\in[1,2]$ and $\epsilon,\delta\in(0,1)$, let
$s=O(\log^{2}(d/\delta)/\epsilon^{2+p})$. There exists a distribution
$\mathcal{R}$ over $s\times(nd)$ matrices, and an algorithm such that for any
vectors $y_{1},\dots,y_{n}\in\mathbbm{R}^{d}$ and weights
$\lambda_{1},\dots,\lambda_{n}\in[0,1]$ with $\sum_{i=1}^{n}\lambda_{i}=1$,
the following occurs with probability at least $1-\delta$:

  * •

We sample $\mathbf{S}=[\mathbf{S}_{1},\dots,\mathbf{S}_{n}]\sim\mathcal{R}$
and we give the algorithm as input $\mathbf{S}$, and the vector
$\sum_{i=1}^{n}\mathbf{S}_{i}(\lambda_{i}^{1/p}y_{i})\in\mathbbm{R}^{s}$.

  * •

The algorithm outputs $d$ numbers
$\boldsymbol{\alpha}_{1},\dots,\boldsymbol{\alpha}_{d}\in\mathbbm{R}_{\geq 0}$
such that each $j\in[d]$ satisfies

$(1-\epsilon)\left(\sum_{i=1}^{n}\lambda_{i}|y_{ij}|^{p}\right)^{1/p}-\mathrm{err}\leq\boldsymbol{\alpha}_{j}\leq(1+\epsilon)\left(\sum_{i=1}^{n}\lambda_{i}|y_{ij}|^{p}\right)^{1/p}+\mathrm{err},$

where $\mathrm{err}\in\mathbbm{R}_{\geq 0}$ is an additive error satisfying

$\mathrm{err}\leq\epsilon\left(\sum_{j=1}^{d}\sum_{i=1}^{n}\lambda_{i}|y_{ij}|^{p}\right)^{1/p}.$

Proof: We combine Corollary 2.3 and Lemma 2.5. We describe the distribution
$\mathcal{R}$ over $s\times(nd)$ matrices by giving a procedure for sampling
$\mathbf{S}\sim\mathcal{R}$. $\mathbf{S}$ can be naturally expressed as a
concatenation of matrices $\mathbf{S}=[\mathbf{S}_{1},\dots,\mathbf{S}_{n}]$.

  * •

We let $\mathcal{M}$ be the distribution over $s_{0}\times n$ matrices of
Corollary 2.3 with error probability at most $\delta/(2dt)$ and accuracy
$\epsilon$ (so that we may union bound over $d$ sketches later) so that
$s_{0}=O(\log(dt/\delta)/\epsilon^{2})$. We take $d$ independent samples
$\mathbf{S}_{1}^{\prime},\dots,\mathbf{S}_{d}^{\prime}\sim\mathcal{M}$.

  * •

For each $j\in[d]$, we let $P_{j}$ be the $n\times(nd)$ matrix where given the
vector $y^{\prime}\in\mathbbm{R}^{nd}$ given by vertically stacking
$\lambda_{1}^{1/p}y_{1},\dots,\lambda_{n}^{1/p}y_{n}\in\mathbbm{R}^{d}$, sets
$y_{\cdot,j}^{\prime}=P_{k}j^{\prime}$, where
$y_{\cdot,j}^{\prime}=\lambda^{1/p}\circ(y_{i,j})_{i\in[n]}\in\mathbbm{R}^{n}$.
Let $P$ be the $(nd)\times(nd)$ matrix which stacks these matrices vertically.

  * •

We sample $\mathbf{C}\sim\mathcal{C}$ as in Lemma 2.5 with $m=s_{0}$, where we
set the accuracy parameter $\epsilon/2$ and the failure probability
$\delta/2$. We let

$\mathbf{S}=\mathbf{C}\cdot\mathrm{diag}(\mathbf{S}_{1}^{\prime},\dots,\mathbf{S}_{d}^{\prime})\cdot
P.$

Intuitively, we will apply our sketch $\mathbf{S}$ on the vector the matrix
$\mathbf{S}$ may be interpreted as first applying $d$ sketches of Corollary
2.3 to the vectors $(\lambda^{1/p}\circ y_{\cdot,1}),\dots,(\lambda^{1/p}\circ
y_{\cdot,d})\in\mathbbm{R}^{n}$, and then applying $\mathbf{C}$ from Lemma
2.5. The algorithm for producing the estimates
$\boldsymbol{\alpha}_{1},\dots,\boldsymbol{\alpha}_{d}$ proceeds by applying
the algorithm of Lemma 2.5 to obtain, for each $j\in[d]$ a sequence of $t$
vectors
$\widehat{\boldsymbol{z}}_{j}^{(1)},\dots,\widehat{\boldsymbol{z}}_{j}^{(t)}\in\mathbbm{R}^{s_{0}}$.
We apply the algorithm of Corollary 2.3 to each of the $t$ vectors to obtain
estimates
$\boldsymbol{\alpha}_{j}^{(1)},\dots,\boldsymbol{\alpha}_{j}^{(t)}\in\mathbbm{R}_{\geq
0}$, and we let
$\boldsymbol{\alpha}_{j}=\mathrm{median}\\{\boldsymbol{\alpha}_{j}^{(t^{\prime})}:t^{\prime}\in[t]\\}$.

To see why this works, consider the collection of $d$ vectors

$\boldsymbol{z}_{j}=\mathbf{S}_{j}^{\prime}(\lambda^{1/p}\circ
y_{\cdot,j})\in\mathbbm{R}^{s_{0}},$

and notice that by Corollary 2.3, every $j\in[d]$ and $\ell\in[s]$,
$\boldsymbol{z}_{j,\ell}\sim(\sum_{i=1}^{n}\lambda_{i}|y_{i,j}|^{p})^{1/p}\cdot\boldsymbol{\chi}_{j,\ell}$,
where $\boldsymbol{\chi}_{k,\ell}$ are independent, $p$-stable random
variables. Indeed, if we write $v\in\mathbbm{R}^{d}$ as the vector which sets

$v_{j}=\left(\sum_{i=1}^{n}\lambda_{i}|y_{i,j}|^{p}\right)^{1/p},$

then vertically stacking the vectors
$\boldsymbol{z}_{1},\dots,\boldsymbol{z}_{d}\in\mathbbm{R}^{s_{0}}$ gives a
vector which is equivalently distributed as $\boldsymbol{\Psi}v$, where
$\boldsymbol{\Psi}$ is the matrix from Lemma 2.5. In particular, with
probability at least $1-\delta/2$, the algorithm of Lemma 2.5 outputs $dt$
vectors
$(\widehat{\boldsymbol{z}}_{j}^{(t^{\prime})}:j\in[d],t^{\prime}\in[t])$ which
satisfy

$\displaystyle\widehat{\boldsymbol{z}}_{j}^{(t^{\prime})}=\mathbf{S}_{k}^{\prime}(\lambda^{1/p}\circ
y_{\cdot,j})+\boldsymbol{\mathrm{err}}_{j}^{(t^{\prime})}\cdot\boldsymbol{\chi}_{j,t^{\prime}}.$
(15)

Hence, with probability at least $1-\delta/(2dt)$, the algorithm of Corollary
2.3 applied to $\widehat{\boldsymbol{z}}_{j}^{(t^{\prime})}$ outputs an
estimate $\boldsymbol{\alpha}_{j}^{(t^{\prime})}$ satisfying
$(1-\epsilon)(v_{j}^{p}+(\boldsymbol{\mathrm{err}}_{j}^{(t^{\prime})})^{p})^{1/p}\leq\boldsymbol{\alpha}_{j}^{(t^{\prime})}\leq(1+\epsilon)(v_{j}^{p}+(\boldsymbol{\mathrm{err}}_{j}^{(t^{\prime})})^{p})^{1/p}$,
and therefore, we have that each $\boldsymbol{\alpha}_{j}^{(t^{\prime})}$
satisfies

$(1-\epsilon)v_{j}-\boldsymbol{\mathrm{err}}_{j}^{(t^{\prime})}\leq\boldsymbol{\alpha}_{j}^{(t^{\prime})}\leq(1+\epsilon)v_{j}+2\cdot\boldsymbol{\mathrm{err}}_{j}^{(t^{\prime})}.$

Since at least $t/2$ of $t^{\prime}\in[t]$ satisfies
$\boldsymbol{\mathrm{err}}_{j}^{(t^{\prime})}\leq\epsilon/2\cdot\|v\|_{p}$,
the median $\boldsymbol{\alpha}_{j}^{(t^{\prime})}$ satisfies the desired
error guarantee. Applying a union bound over all $dt$ applications of
Corollary 2.3 and Lemma 2.5 gives the desired guarantees.

######  Lemma 2.7 ($p$-stable Sketch Recovery for Optimizer).

For $n,d\in\mathbbm{N}$, $p\in[1,2]$ and $\epsilon,\delta\in(0,1)$, let
$s=O(\log(d/\delta)\cdot\log(\log d/(\epsilon\delta))/\epsilon^{2+p})$. There
exists a distribution $\mathcal{O}$ over $s\times(nd)$ matrices, and an
algorithm such that for any vectors $y_{1},\dots,y_{n}\in\mathbbm{R}^{d}$ and
any set of weights $\lambda_{1},\dots,\lambda_{n}\in[0,1]$ where
$\sum_{i=1}^{n}\lambda_{i}=0$, whenever $\sum_{i=1}^{n}\lambda_{i}y_{i}=0$,
the following occurs with probability at least $1-\delta$:

  * •

We sample $\mathbf{S}=[\mathbf{S}_{1},\dots,\mathbf{S}_{n}]\sim\mathcal{O}$
and we give the algorithm as input $\mathbf{S}$,
$\sum_{i=1}^{n}\mathbf{S}_{i}(\lambda_{i}^{1/p}y_{i})$, the parameters
$\lambda_{1},\dots,\lambda_{n}$, an index $j_{0}\in[d]$, and a parameter
$\gamma\in\mathbbm{R}_{\geq 0}$ satisfying

$(1-\epsilon)\left(\sum_{i=1}^{n}\lambda_{i}|y_{ij_{0}}|^{p}\right)^{1/p}\leq\gamma\leq(1+\epsilon)\left(\sum_{i=1}^{n}\lambda_{i}|y_{ij_{0}}|^{p}\right)^{1/p}.$

  * •

The algorithm outputs a parameter
$\widehat{\boldsymbol{\beta}}\in\mathbbm{R}_{\geq 0}$ which satisfies

$(1-\epsilon)\min_{z\in\mathbbm{R}}\left(\sum_{i=1}^{n}\lambda_{i}|y_{ij_{0}}-z|^{p}\right)^{1/p}-\mathrm{err}\leq\widehat{\boldsymbol{\beta}}\leq(1+\epsilon)\min_{z\in\mathbbm{R}}\left(\sum_{i=1}^{n}\lambda_{i}|y_{ij_{0}}-z|^{p}\right)^{1/p}+\mathrm{err},$

where $\mathrm{err}\in\mathbbm{R}_{\geq 0}$ is an additive error satisfying

$\mathrm{err}\leq\epsilon\left(\sum_{j=1}^{d}\sum_{i=1}^{n}\lambda_{i}|y_{i,j}|^{p}\right)^{1/p}.$

Proof: The proof follows similarly to that of Lemma 2.6; the only difference
is that instead of using the sketch of Corollary 2.3, we use the sketch of
Lemma 2.4. For completeness, we describe the distribution $\mathcal{O}$ over
$s\times(nd)$ matrices by giving a procedure for sampling
$\mathbf{S}\sim\mathcal{O}$:

  * •

We let $\mathcal{M}$ be the distribution over $s_{0}\times n$ matrices from
Lemma 2.4 with accuracy $\epsilon$ and failure probability $\delta/(2t)$,
where $s_{0}=O(\log(t/(\epsilon\delta))/\epsilon^{2})$. We take $d$
independent samples
$\mathbf{S}_{1}^{\prime},\dots,\mathbf{S}_{d}^{\prime}\sim\mathcal{M}$. Note
that even though we take $d$ independent samples, we will only require that
the sketch $t$ evaluations of the $\mathbf{S}_{j_{0}}\sim\mathcal{M}$ succeed
(hence, we amplify the error probability to $\delta/(2t)$, as opposed to
$\delta/(2td)$ as in Lemma 2.6).

  * •

We sample $\mathbf{C}\sim\mathcal{C}$ as in Lemma 2.5 with $m=s_{0}$, where we
set the accuracy parameter $\epsilon/2$ and failure probability $\delta/2$.
Recalling the definition of $P$ (see Item 2 in the proof of Lemma 2.6, we let

$\mathbf{S}=\mathbf{C}\cdot\mathrm{diag}(\mathbf{S}_{1}^{\prime},\dots,\mathbf{S}_{d}^{\prime})\cdot
P.$

Similarly to the proof of Lemma 2.6, $\mathbf{S}$ may be interpreted as
applying the sketch of Lemma 2.5 to $d$ vectors in $\mathbbm{R}^{s_{0}}$, each
$j\in[d]$ of which is an independent sketch
$\mathbf{S}_{j}^{\prime}(\lambda^{1/p}\circ y_{\cdot
j})\in\mathbbm{R}^{s_{0}}$, where $\mathbf{S}_{j}^{\prime}\sim\mathcal{M}$ is
the sketch of Lemma 2.4. Again, we consider the collection of $d$ vectors
$\boldsymbol{z}_{j}=\mathbf{S}_{j}^{\prime}(\lambda^{1/p}\circ
y_{\cdot,j})\in\mathbbm{R}^{s_{0}}$, for all $j\in[d]$, and by Lemma 2.4,
every $j\in[d]$ has $\boldsymbol{z}_{j}\sim\|\lambda^{1/p}\circ
y_{\cdot,j}\|_{p}\cdot\boldsymbol{\Psi}_{j}\in\mathbbm{R}^{s_{0}}$, where
$\boldsymbol{\Psi}_{j}$ is an independent, $p$-stable random vector. Writing
$v\in\mathbbm{R}^{d}$ by $v_{j}=\|\lambda^{1/p}\circ y_{\cdot,j}\|_{p}$, and
we apply the algorithm of Lemma 2.5 and focus on the $t$ vectors
$\widehat{\boldsymbol{z}}_{j_{0}}^{(1)},\dots,\widehat{\boldsymbol{z}}_{j_{0}}^{(t)}\in\mathbbm{R}^{s_{0}}$
which satisfy

$\widehat{\boldsymbol{z}}_{j_{0}}^{(t^{\prime})}=\mathbf{S}_{j_{0}}^{\prime}(\lambda^{1/p}\circ
y_{\cdot,j_{0}})+\boldsymbol{\mathrm{err}}_{j_{0}}^{(t^{\prime})}\cdot\boldsymbol{\chi}_{j,t^{\prime}}.$

We apply the algorithm of Lemma 2.4 to each of the vectors
$\widehat{\boldsymbol{z}}_{j_{0}}\in\mathbbm{R}^{s}$, while giving as input
the parameter $\gamma$ to obtain the estimate
$\widehat{\boldsymbol{\beta}}^{(1)},\dots,\widehat{\boldsymbol{\beta}}^{(t)}$.
Then, we set
$\widehat{\boldsymbol{\beta}}=\mathrm{median}\\{\widehat{\boldsymbol{\beta}}^{(t^{\prime})}:t^{\prime}\in[t]\\}$.
Similarly to the proof of Lemma 2.6, $\widehat{\boldsymbol{\beta}}$ provides
the desired approximation guarantees.

###  2.4 Proof of Lemma 2.1

The distribution over matrices $\mathcal{S}$ will be given by utilizing the
sketches in Lemmas 2.6 and 2.7. In particular, let $\mathcal{R}$ denote the
distribution over $s_{1}\times(nd)$ matrices from Lemma 2.6 where
$s_{1}=O(\log^{2}(d/(\delta_{1}\epsilon_{1}))/\epsilon_{1}^{2+p})$ using
failure probability at most $\delta_{1}=\delta\epsilon 2^{p}/3$ and accuracy
$\epsilon_{1}$, and let $\mathcal{O}$ denote the distribution over
$s_{2}\times(nd)$ matrices from Lemma 2.7 where
$s_{2}=O(\log(d/\delta_{2})\log(\log
d/(\epsilon_{2}\delta_{2}))/\epsilon_{2}^{2+p})$ using failure probability at
most $\delta_{2}=\delta/3$ and accuracy $\epsilon_{2}$. Specifically, we let

$\epsilon_{1}=c_{1}\cdot\epsilon^{1+2/p}\cdot\delta^{1+1/p}\qquad\text{and}\qquad\epsilon_{2}=c_{2}\cdot\epsilon^{1+1/p}\cdot\delta^{1/p},$

for small enough constant $c_{1},c_{2}>0$. Thus, the total space

$\displaystyle
s=s_{1}+s_{2}=O\left(\frac{\log^{2}(d/(\epsilon\delta))}{\epsilon^{4+p+4/p}\delta^{3+p+2/p}}\right)+O\left(\frac{\log(d/\delta)\log(\log
d/(\epsilon\delta))}{\epsilon^{3+p+2/p}}\right).$ (16)

A draw $\mathbf{S}=[\mathbf{S}_{1},\dots,\mathbf{S}_{n}]\sim\mathcal{S}$ is an
$(s_{1}+s_{2})\times(nd)$ matrix which is generated as follows:

  * •

We first sample a collection
$\boldsymbol{u}_{1},\dots,\boldsymbol{u}_{d}\sim\mathrm{Exp}(1)$ independently
at random. We let
$\mathbf{U}=\mathrm{diag}(1/\boldsymbol{u}_{1}^{1/p},\dots,1/\boldsymbol{u}_{d}^{1/p})$
be the $d\times d$ matrix with entries $1/\boldsymbol{u}_{j}^{1/p}$ across the
diagonal, and let $\mathbf{V}=\mathrm{diag}(\mathbf{U},\dots,\mathbf{U})$
denote the $(nd)\times(nd)$ matrix placing $n$ copies of the matrix
$\mathbf{U}$ in the diagonals.

  * •

Then, we sample
$\mathbf{S}^{(1)}=[\mathbf{S}^{(1)}_{1},\dots,\mathbf{S}^{(1)}_{n}]\sim\mathcal{R}$
and
$\mathbf{S}^{(2)}=[\mathbf{S}^{(2)}_{1},\dots,\mathbf{S}^{(2)}_{n}]\sim\mathcal{O}$.
The final matrix is given by

$\mathbf{S}=\left[\begin{array}[]{c}\mathbf{S}^{(1)}\\\
\mathbf{S}^{(2)}\end{array}\right]\cdot\mathbf{V}.$

Intuitively, the above operation consists of keeping two linear sketches: one
of $\sum_{i=1}^{n}\mathbf{S}^{(1)}_{i}(\lambda_{i}^{1/p}\boldsymbol{y}_{i})$,
and one of
$\sum_{i=1}^{n}\mathbf{S}^{(2)}_{i}(\lambda_{i}^{1/p}\boldsymbol{y}_{i})$,
where $\boldsymbol{y}_{1},\dots,\boldsymbol{y}_{n}\in\mathbbm{R}^{d}$ are
given by letting, for each $j\in[d]$,

$\displaystyle\boldsymbol{y}_{ij}$
$\displaystyle=\frac{1}{\boldsymbol{u}_{j}^{1/p}}\cdot x_{ij}.$ (17)

We note that $\sum_{i=1}^{n}\lambda_{i}x_{i}=0$ implies
$\sum_{i=1}^{n}\lambda_{i}\boldsymbol{y}_{i}=0$, as we simply re-scaled all
coordinates of $x_{\cdot,j}$ by the same value. We now specify the description
of the algorithm for outputting the triple
$(\widehat{\boldsymbol{j}},\widehat{\boldsymbol{\alpha}},\widehat{\boldsymbol{\beta}})$,
when given as input $\mathbf{S}$ and
$\sum_{i=1}^{n}\mathbf{S}_{i}(\lambda_{i}^{1/p}x_{i})$. We will proceed in the
following way:

  1. 1.

We will first apply the algorithm of Lemma 2.6 with the matrix
$\mathbf{S}^{(1)}$ and the vector
$\sum_{i=1}^{n}\mathbf{S}^{(1)}_{i}(\lambda_{i}^{1/p}\boldsymbol{y}_{i})$.
Notice that this algorithm outputs a sequence of $d$ numbers,
$\boldsymbol{\alpha}_{1},\dots,\boldsymbol{\alpha}_{d}\in\mathbbm{R}_{\geq
0}$. We will set

$\displaystyle\widehat{\boldsymbol{j}}\stackrel{{\scriptstyle\rm
def}}{{=}}\mathop{\mathrm{argmax}}_{j\in[d]}\boldsymbol{\alpha}_{j}\qquad\text{and}\qquad\widehat{\boldsymbol{\alpha}}\stackrel{{\scriptstyle\rm
def}}{{=}}\boldsymbol{\alpha}_{\widehat{\boldsymbol{j}}}\cdot\boldsymbol{u}_{\widehat{\boldsymbol{j}}}^{1/p}.$
(18)

  2. 2.

Then, we will apply the algorithm of Lemma 2.7 with the matrix
$\mathbf{S}^{(2)}$ and the vector
$\sum_{i=1}^{n}\mathbf{S}^{(2)}(\lambda_{i}^{1/p}\boldsymbol{y}_{i})$,
parameters $\lambda_{1},\dots,\lambda_{n}$, and the input
$\widehat{\boldsymbol{j}}$ and $\gamma=\widehat{\boldsymbol{\alpha}}$. The
algorithm produces an output $\boldsymbol{\beta}$, and we let

$\widehat{\boldsymbol{\beta}}\stackrel{{\scriptstyle\rm
def}}{{=}}\boldsymbol{\beta}\cdot\boldsymbol{u}_{\boldsymbol{j}}^{1/p}.$

We now consider the following event $\mathcal{E}$, which is a function of the
random variables $\boldsymbol{u}_{1},\dots,\boldsymbol{u}_{d}$ and
$x_{1},\dots,x_{n}\in\mathbbm{R}^{d}$. The goal is to show that (1) the event
$\mathcal{E}$ is satisfied for a random
$\boldsymbol{u}_{1},\dots,\boldsymbol{u}_{d}\sim\mathrm{Exp}(1)$ with high
probability, and (2) if the event $\mathcal{E}$ is satisfied, then
$(\widehat{\boldsymbol{j}},\widehat{\boldsymbol{\alpha}},\widehat{\boldsymbol{\beta}})$
satisfies the conditions of Lemma 2.1.

######  Definition 2.8.

Consider a fixed setting of $u_{1},\dots,u_{d}\in\mathbbm{R}_{>0}$, and let

$w_{j}\stackrel{{\scriptstyle\rm
def}}{{=}}\frac{1}{u_{j}}\cdot\sum_{i=1}^{n}\lambda_{i}|x_{ij}|^{p},\qquad
j^{*}\stackrel{{\scriptstyle\rm
def}}{{=}}\mathop{\mathrm{argmax}}_{j\in[d]}w_{j},\qquad\text{and}\qquad
j^{**}\stackrel{{\scriptstyle\rm
def}}{{=}}\mathop{\mathrm{argmax}}_{j\in[d]\setminus\\{j^{*}\\}}w_{j}.$

Finally, we also let

$\eta\stackrel{{\scriptstyle\rm def}}{{=}}\epsilon\delta 2^{-p}/1000.$

We say the event $\mathcal{E}$ is satisfied for the setting of
$u_{1},\dots,u_{d}$ when the following conditions hold:

$\displaystyle w_{j^{*}}$ $\displaystyle\geq\eta\sum_{j=1}^{d}w_{j}$ (19)
$\displaystyle w_{j^{*}}$ $\displaystyle\geq(1+\eta)w_{j^{**}}.$ (20)

######  Claim 2.9.

Consider a random setting of
$\boldsymbol{u}_{1},\dots,\boldsymbol{u}_{d}\sim\mathrm{Exp}(1)$, and let
$\boldsymbol{j}^{*}\in[d]$ be the random variable (which depends on
$\boldsymbol{u}_{1},\dots,\boldsymbol{u}_{d}$) which sets $\boldsymbol{j}^{*}$
to the maximum $\boldsymbol{w}_{j}$, as in Definition 2.8. Then,
$\boldsymbol{j}^{*}$ is drawn from $\mathcal{D}$.

Proof: The proof is a simple identity of the exponential distribution. See
Fact 5.3 of [CJLW22] for the complete computation.

######  Claim 2.10.

If we sample $\boldsymbol{u}_{1},\dots,\boldsymbol{u}_{d}\sim\mathrm{Exp}(1)$,
the event $\mathcal{E}$ is satisfied with probability at least
$1-\epsilon\delta 2^{-p}/3$.

Proof: This follows from Lemma 5.5 of [CJLW22]. Namely, consider the vector
$z\in\mathbbm{R}^{d}$ where $z_{j}=\sum_{i=1}^{n}\lambda_{i}|x_{ij}|^{p}$.
Lemma 5.5 implies $|z_{j^{*}}|/\boldsymbol{u}_{j^{*}}\geq\eta\|z\|_{1}$ and
$|z_{j^{*}}|/\boldsymbol{u}_{j^{*}}\geq(1+\eta)|z_{j^{**}}|/\boldsymbol{u}_{j^{**}}$
with probability at least $1-4\eta$, giving us (19) and (20).

######  Claim 2.11.

Let $u_{1},\dots,u_{d}\in\mathbbm{R}_{>0}$ be a fixed setting where event
$\mathcal{E}$ is satisfied. Then, let $y_{1},\dots,y_{n}$ be given from
$x_{1},\dots,x_{n}$ by (17), and run the algorithm of Lemma 2.6 while given
access to $\mathbf{S}^{(1)}$ and
$\sum_{i=1}^{n}\mathbf{S}^{(1)}_{i}(\lambda_{i}^{1/p}y_{i})$, to output
$\boldsymbol{\alpha}_{1},\dots,\boldsymbol{\alpha}_{d}$. Then, with
probability at least $1-\delta_{1}$, we have

$\displaystyle\mathop{\mathrm{argmax}}_{j\in[d]}\boldsymbol{\alpha}_{j}$
$\displaystyle=j^{*}\qquad\text{and}$
$\displaystyle(1-\epsilon_{2})\left(\sum_{i=1}^{n}\lambda_{i}|y_{ij^{*}}|^{p}\right)^{1/p}$
$\displaystyle\leq\boldsymbol{\alpha}_{j^{*}}\leq(1+\epsilon_{2})\left(\sum_{i=1}^{n}\lambda_{i}|y_{ij^{*}}|^{p}\right)^{1/p}$

Proof: We first note that applying Lemma 2.6, as well as condition (19), with
probability at least $1-\delta_{1}$, the algorithm outputs
$\boldsymbol{\alpha}_{1},\dots,\boldsymbol{\alpha}_{d}$ satisfying

$\displaystyle(1-\epsilon_{1})w_{j}^{1/p}-\frac{\epsilon_{1}}{\eta^{1/p}}\cdot
w_{j^{*}}^{1/p}\leq\boldsymbol{\alpha}_{j}\leq(1+\epsilon_{1})w_{j}^{1/p}+\frac{\epsilon_{1}}{\eta^{1/p}}\cdot
w_{j^{*}}^{1/p}.$ (21)

In other words, we use (20) to say

$\displaystyle\left(1-\epsilon_{1}-\frac{\epsilon_{1}}{\eta^{1/p}}\right)\cdot
w_{j^{*}}^{1/p}\leq\boldsymbol{\alpha}_{j^{*}}\qquad\text{and, for all $j\neq
j^{*}$,}\quad\boldsymbol{\alpha}_{j}\leq\boldsymbol{\alpha}_{j^{**}}$
$\displaystyle\leq(1+\epsilon_{1})w_{j^{**}}^{1/p}+\frac{\epsilon_{1}}{\eta^{1/p}}\cdot
w_{j^{*}}^{1/p}$
$\displaystyle\leq\left(\frac{1}{1+\eta/(2p)}+\epsilon_{1}+\frac{\epsilon_{1}}{\eta^{1/p}}\right)w_{j^{*}}^{1/p}.$

Hence, the fact
$\mathop{\mathrm{argmax}}_{j\in[d]}\boldsymbol{\alpha}_{j}=j^{*}$ follows from
the fact

$\frac{1}{1+\eta/(2p)}+\epsilon_{1}+\frac{\epsilon_{1}}{\eta^{1/p}}\leq
1-\epsilon_{1}-\frac{\epsilon_{1}}{\eta^{1/p}}$

when $\epsilon_{1}\leq c\eta^{1+1/p}/p$, for a small constant $c>0$. The
condition on $\boldsymbol{\alpha}_{j^{*}}$ then follows from (21) since
$\epsilon_{1}\left(1+1/\eta^{1/p}\right)\leq\epsilon_{2}$.

######  Claim 2.12.

Let $u_{1},\dots,u_{d}\in\mathbbm{R}_{>0}$ be a fixed setting where
$\mathcal{E}$ is satisfied. Then, let $y_{1},\dots,y_{n}$ be given from
$x_{1},\dots,x_{n}$ according to (17), and notice
$\sum_{i=1}^{n}\lambda_{i}y_{i}=0$. Run the algorithm of Lemma 2.7 while given
access to $\mathbf{S}^{(2)}$,
$\sum_{i=1}^{n}\mathbf{S}^{(2)}_{i}(\lambda_{i}^{1/p}y_{i})$, the parameters
$\lambda_{1},\dots,\lambda_{n}$,$j^{*}$, and a parameter
$\gamma\in\mathbbm{R}_{>0}$ (which is set to $\widehat{\boldsymbol{\alpha}}$)
satisfying

$(1-\epsilon_{2})\left(\sum_{i=1}^{n}\lambda_{i}|y_{ij^{*}}|^{p}\right)^{1/p}\leq\gamma\leq(1+\epsilon_{2})\left(\sum_{i=1}^{n}\lambda_{i}|y_{ij^{*}}|^{p}\right)^{1/p}.$

With probability at least $1-\delta_{2}$, the algorithm outputs a parameter
$\boldsymbol{\beta}\in\mathbbm{R}_{\geq 0}$ which satisfies

$(1-\epsilon)\min_{z\in\mathbbm{R}}\left(\sum_{i=1}^{n}\lambda_{i}|y_{ij^{*}}-z|^{p}\right)^{1/p}\leq\boldsymbol{\beta}\leq(1+\epsilon)\min_{z\in\mathbbm{R}}\left(\sum_{i=1}^{n}\lambda_{i}|y_{ij^{*}}-z|^{p}\right)^{1/p}$

Proof: The proof proceeds similarly to above. Namely, the algorithm of Lemma
2.7 outputs a parameter $\widehat{\boldsymbol{\beta}}$ which is a
multiplicative $(1\pm\epsilon_{2})$-approximation to
$\min_{z\in\mathbbm{R}}(\sum_{i=1}^{n}|y_{ij^{*}}-z|^{p})^{1/p}$. with an
additive error of
$\epsilon_{2}(\sum_{j=1}^{d}\sum_{i=1}^{n}\lambda_{i}|y_{i,j}|^{p})^{1/p}$.
Hence, we will upper bound the additive error.

$\displaystyle\sum_{j=1}^{d}\sum_{i=1}^{n}\lambda_{i}|y_{i,j}|^{p}$
$\displaystyle=\sum_{j=1}^{d}w_{j}\leq\frac{1}{\eta}\cdot w_{j^{*}},$

and the right-most expression is exactly
$\frac{1}{\eta}\sum_{i=1}^{n}\lambda_{i}|y_{ij^{*}}|^{p}$. Specifically, the
additive error incurred is at most

$\frac{\epsilon_{2}}{\eta^{1/p}}\left(\sum_{i=1}^{n}\lambda_{i}|y_{ij^{*}}|^{p}\right)^{1/p}\leq\frac{2\epsilon_{2}}{\eta^{1/p}}\cdot\min_{z\in\mathbbm{R}}\left(\sum_{i=1}^{n}\lambda_{i}|y_{ij^{*}}-z|^{p}\right)^{1/p}.$

In particular, since $\epsilon_{2}+2\epsilon_{2}/\eta^{1/p}\leq\epsilon$, we
obtain the desired claim.

We note that combining Claims 2.9, 2.10, 2.11, and 2.12 will give the desired
lemma. With probability $1-\epsilon\delta 2^{-p}/3$, event $\mathcal{E}$ is
satisfied by Claim 2.10, and hence Claim 2.11 imply that in the triple
$(\widehat{\boldsymbol{k}},\widehat{\boldsymbol{\alpha}},\widehat{\boldsymbol{\beta}})$,
we have $\widehat{\boldsymbol{j}}=j^{*}$ with probability at least
$1-\delta_{1}=1-\epsilon\delta 2^{-p}/3$. Furthermore, we always have, by
definition of $\boldsymbol{y}$,

$\displaystyle\boldsymbol{u}_{j^{*}}^{1/p}\left(\sum_{i=1}^{n}\lambda_{i}|\boldsymbol{y}_{ij^{*}}|^{p}\right)^{1/p}$
$\displaystyle=\left(\sum_{i=1}^{n}\lambda_{i}|x_{ij^{*}}|^{p}\right)^{1/p}$
$\displaystyle\boldsymbol{u}_{j^{*}}^{1/p}\min_{z\in\mathbbm{R}}\left(\sum_{i=1}^{n}\lambda_{i}|\boldsymbol{y}_{ij^{*}}-z|^{p}\right)^{1/p}$
$\displaystyle=\min_{z\in\mathbbm{R}}\left(\sum_{i=1}^{n}\lambda_{i}|x_{ij^{*}}-z|^{p}\right)^{1/p}.$

In particular, Claim 2.11 implies $\widehat{\boldsymbol{\alpha}}$ is a
$(1\pm\epsilon_{2})$-multiplicative approximation. This establishes the
approximation requirement of $\widehat{\boldsymbol{\alpha}}$, but it also
allows us to utilize it in Claim 2.12, which implies the desired approximation
on $\widehat{\boldsymbol{\beta}}$. Finally, the probability that
$\widehat{\boldsymbol{j}}=j^{*}$ is at least the probability that event
$\mathcal{E}$ is satisfied, and the sketch of Lemma 2.6 succeeds, which both
occur with probability at least $1-\epsilon 2^{-p}$ by the setting of the
failure probabilities on these events. By Claim 2.9, the variation distance of
the random variable $\widehat{\boldsymbol{j}}$ from $\mathcal{D}$ is at most
$\epsilon 2^{-p}$.

###  2.5 Storing the randomness

The previous section described the $\ell_{p}^{p}$-median sketch assuming
infinite precision on the random variables, as well as complete access to the
randomness used. Here, we remark on how to: 1) utilize bounded precision on
the entries of $\mathbf{S}$ (to avoid storing infinite-precision real
numbers), and 2) utilize Nisan’s pseudorandom generator in order to maintain
$\mathbf{S}$ implicitly. This technique, first used in [Ind06] for the
$\ell_{p}$-sketch, is by now standard in the streaming literature.

#### Bounded Precision

We assume that the input vectors $x_{1},\dots,x_{n}\in\mathbbm{R}^{d}$ and the
weights $\lambda_{1},\dots,\lambda_{n}$ have polynomially bounded entries
(i.e., every number lies $[-\mathrm{poly}(nd),\mathrm{poly}(nd)]$ and is
specified by $O(\log(nd))$ bits). The only infinite-precision numbers are the
$p$-stable random variables and the exponential random variables, so we will
discretize them so they require only $O(\log(nd))$ random bits to generate
them. Similarly to Claim 3 of [Ind06] our algorithm’s output, as a function of
the $p$-stable and exponential random variables, always considers polynomially
bounded entries and has bounded derivatives in all but a
$1/\mathrm{poly}(nd)$-fraction of the probability space. With probability
$1/\mathrm{poly}(nd)$, the $p$-stable random variables may become too large or
exponentials too small. Furthermore, for any fixed setting of the $p$-stable
random variables, a discontinuity in the algorithm’s output occurs when the
exponential random variables cause $\boldsymbol{k}$ may suddenly switch, when
the maximum of $\boldsymbol{\alpha}_{k}$ in (18) is not unique. When we
discretize the exponential distribution, the probability this occurs is at
most $1/\mathrm{poly}(nd)$.

#### Using Nisan’s Pseudorandom Generator

Instead of utilizing $snd\cdot O(\log(nd))$ random bits to generate the matrix
$\mathbf{S}$, we show that we may utilize Nisan’s pseudorandom generator. This
technique is standard when the sketch is linear: we may re-order the elements
of the stream so as to update the algorithm with a finite-state machine which
reads the random bits in $O(\log(nd))$-sized chunks. The minor subtlety is
that algorithm of Lemma 2.1 requires access to $\mathbf{S}$ in order to
produce its output $(\boldsymbol{j},\boldsymbol{\alpha},\boldsymbol{\beta})$,
and hence refers back to previous random bits.

There is a simple fix: for each fixed setting of the input
$x_{1},\dots,x_{n}\in(\\{0,1\\}^{O(\log(nd))})^{nd}$ and weights
$\lambda_{1},\dots,\lambda_{n}\in(\\{0,1\\}^{O(\log(nd))})^{n}$, and
$j^{\prime}\in[d]$ and $y\in\mathbbm{R}$ specified by $O(\log(nd))$ bits,
there exists a low-space finite state machine which reads the randomness in
$O(\log(nd))$ chunks and checks whether $\boldsymbol{j}=j^{\prime}$, and if
so, checks that $\boldsymbol{\alpha}$ satisfies the guarantees of Lemma 2.1,
and accurately estimates $\boldsymbol{\beta}$ if $z=y$. The finite state
machine receives its input in the order of the coordinates; this way, it can
temporarily store $\boldsymbol{u}_{j}$ to compute $\boldsymbol{\alpha}_{j}$
and find (18). Furthermore, while processing coordinate $j^{\prime}$, it
checks the cost of using $y$ as a center. Since Nisan’s pseudorandom generator
fools every such finite state machine, and we try at most
$O(d)\cdot\mathrm{poly}(\log d,1/\epsilon)$ possible coordinates $j^{\prime}$
and centers $y$, a union bound will show that the output of the algorithm
under true random bits and under Nisan’s pseudorandom generator differs by at
most $1/\mathrm{poly}(nd/\epsilon)$ in total variation distance.

##  3 Streaming $(k,p)$-Clustering Costs in $\ell_{p}$

We give the following application of our $\ell_{p}^{p}$-median sketch to
estimating the $(k,p)$-clustering in $\ell_{p}$ spaces on a stream for
$p\in[1,2]$. The space complexity of our streaming algorithm will be
$\mathrm{poly}(\log(nd),k,1/\epsilon)$, and will work in insertion-only
streams.

######  Theorem 6.

Fix $n,d,k\in\mathbbm{N}$, as well as $p\in[1,2]$ and $\epsilon\in(0,1)$.
There exists an insertion-only streaming algorithm which processes a set of
$n$ points $x_{1},\dots,x_{n}\in\mathbbm{R}^{d}$ utilizing
$\mathrm{poly}(\log(nd),k,1/\epsilon)$ space, and outputs a parameter
$\boldsymbol{\eta}\in\mathbbm{R}$ which satisfies

$(1-\epsilon)\min_{\begin{subarray}{c}C_{1},\dots,C_{k}\\\ \text{partition
$[n]$}\end{subarray}}\sum_{\ell=1}^{k}\min_{c_{\ell}\in\mathbbm{R}^{d}}\sum_{i\in
C_{\ell}}\|x_{i}-c_{\ell}\|_{p}^{p}\leq\boldsymbol{\eta}\leq(1+\epsilon)\min_{\begin{subarray}{c}C_{1},\dots,C_{k}\\\
\text{partition
$[n]$}\end{subarray}}\sum_{\ell=1}^{k}\min_{c_{\ell}\in\mathbbm{R}^{d}}\sum_{i\in
C_{\ell}}\|x_{i}-c_{\ell}\|_{p}^{p}$

with probability at least $0.9$.

###  3.1 Preliminaries: Coreset of [HV20] and $\ell_{p}^{p}$-median sketch

We will utilize the following result of [HV20].

######  Definition 3.1 (Strong Coresets).

Let $x_{1},\dots,x_{n}\in\mathbbm{R}^{d}$ be a set of points and $p\in[1,2]$.
The set $S\subset[n]$ and $w\colon S\to\mathbbm{R}_{\geq 0}$ is an
$\epsilon$-_strong coreset_ for $(k,z)$-clustering in $\ell_{p}$ if for any
set of $k$ points $c_{1},\dots,c_{\ell}\in\mathbbm{R}^{d}$,

$\displaystyle\sum_{i=1}^{n}\min_{\ell\in[k]}\|x_{i}-c_{\ell}\|_{p}^{z}\approx_{1\pm\epsilon}\sum_{i\in
S}w(i)\min_{\ell\in[k]}\|x_{i}-c_{\ell}\|_{p}^{z}.$

######  Corollary 3.2 (Corollary 5.18 of [HV20]).

For any $p\in[1,2]$, there exists a randomized algorithm that, given
$O(1)$-approximate $\ell_{p}$-distances between $n$ points
$x_{1},\dots,x_{n}\in\mathbbm{R}^{d}$, an integer $k\in\mathbbm{N}$,
$z\in[1,2]$, and $\epsilon,\delta\in(0,1/2)$, outputs a (random) subset
$\mathbf{S}\subset[n]$ of $\mathrm{poly}(k\log(1/\delta)/\epsilon)$ points and
a set of (random) weights $\boldsymbol{w}\colon\mathbf{S}\to\mathbbm{R}_{\geq
0}$ which is a strong coreset for $(k,p)$-clustering in $\ell_{p}$ with
probability at least $1-\delta$.888Corollary 5.18 is stated with $z\geq 2$;
however, the proof works for any $z\in[1,2]$ as well.

It is important that the algorithm of Corollary 3.2 only needs (approximate)
access to distances between the $n$ points, as even storing a single point
requires $\Omega(d)$ space. In particular, the $d$-dimensional representation
of the coreset points is never stored. Even though [HV20] does not explicitly
mention this in Corollary 5.18, it can be verified by inspecting Algorithm 1.
One additional fact is needed in Step 2 of Algorithm 1 of [HV20]: a
$O(1)$-approximate center set $C^{*}$ may be found using a subset of the
dataset points when $z\geq 1$ is not too large.

######  Fact 3.3.

Consider any $k\in\mathbbm{N}$, $z\geq 1$ and $p\in[1,\infty)$, as well as
points $x_{1},\dots,x_{n}\in\mathbbm{R}^{d}$. There exists a set of centers
$c_{1}^{*},\dots,c_{k}^{*}\in\\{x_{1},\dots,x_{n}\\}$ which give a
$2^{z}$-approximation to the $(k,z)$-clustering in $\ell_{p}$.

Proof: Let $c_{1},\dots,c_{k}\in\mathbbm{R}^{d}$ denote the optimal centers
for $(k,z)$-clustering $x_{1},\dots,x_{n}$, and let $C_{1},\dots,C_{k}$ denote
the partition of $[n]$ from the center set $c_{1},\dots,c_{k}$, i.e., for
$\ell\in[k]$,

$C_{\ell}=\left\\{i\in[n]:\forall\ell^{\prime}\in[k],\|x_{i}-c_{\ell}\|_{p}\leq\|x_{i}-c_{\ell^{\prime}}\|_{p}\right\\},$

with ties broken arbitrarily. For each $\ell\in[k]$, let
$\boldsymbol{i}_{\ell}\sim C_{\ell}$ be a uniformly random index from
$C_{\ell}$, and let $\boldsymbol{y}_{\ell}=x_{\boldsymbol{i}_{\ell}}$. Then,
we show that the expected cost of using $\boldsymbol{y}_{\ell}$ as a center
for cluster $C_{\ell}$ cannot increase the cost significantly,

$\displaystyle\mathop{{\bf E}\/}_{\boldsymbol{i}_{\ell}\sim
C_{\ell}}\left[\sum_{i\in
C_{\ell}}\|x_{i}-\boldsymbol{y}_{\ell}\|_{p}^{z}\right]$
$\displaystyle\leq\mathop{{\bf E}\/}_{\boldsymbol{i}_{\ell}\sim
C_{\ell}}\left[\sum_{i\in
C_{\ell}}\left(\|x_{i}-c_{\ell}\|_{p}+\|\boldsymbol{y}_{\ell}-c_{\ell}\|_{p}\right)^{z}\right]\qquad\text{(H\"{o}lder
Inequality)}$ $\displaystyle\leq 2^{z-1}\sum_{i\in
C_{\ell}}\|x_{i}-c_{\ell}\|_{p}^{z}+\mathop{{\bf
E}\/}_{\boldsymbol{i}_{\ell}\sim C_{\ell}}\left[|C_{\ell}|\cdot
2^{z-1}\|\boldsymbol{y}_{\ell}-c_{\ell}\|_{p}^{z}\right]=2^{z}\sum_{i\in
C_{\ell}}\|x_{i}-c_{\ell}\|_{p}^{z}.$

Hence, for each cluster $\ell\in[k]$, there exists a dataset point
$c_{\ell}^{*}$ such that, utilizing $c_{\ell}^{*}$ instead of the optimal
center $c_{\ell}$ incurs at most a factor of $2^{z}$ in the cost of that
cluster, and therefore, the $(k,z)$-clustering cost can be at most $2^{z}$
times more than with the optimal centers $c_{1},\dots,c_{k}$.

We now re-state Theorem 5 for the case that the weights
$\lambda_{1},\dots,\lambda_{n}$ are unknown at time of processing.

######  Lemma 3.4.

Fix $m,d\in\mathbbm{N}$, as well as $p\in[1,2]$ and $\epsilon,\delta\in(0,1)$,
let $s=\mathrm{poly}(\log d,1/\epsilon,\log(1/\delta))$. There exists a
distribution $\mathcal{S}$ over $s\times(md)$ matrices and an algorithm which
satisfies the following.

  * •

For any set of points $x_{1},\dots,x_{m}\in\mathbbm{R}^{d}$, and any set of
weights $\lambda_{1},\dots,\lambda_{n}\in[0,1]$ with
$\sum_{i=1}^{n}\lambda_{i}=1$, we sample
$\mathbf{S}=[\mathbf{S}_{1},\dots,\mathbf{S}_{n}]\sim\mathcal{S}$ and give the
algorithm as input the matrix $\mathbf{S}$, the weights
$\lambda_{1},\dots,\lambda_{m}$, and the vectors

$\mathbf{S}_{j}x_{i}\in\mathbbm{R}^{s}\qquad\text{ for all $i,j\in[m]$}.$

  * •

With probability $1-\delta$ over the draw of $\mathbf{S}$, the algorithm
outputs a parameter $\boldsymbol{\eta}\in\mathbbm{R}_{\geq 0}$ which satisfies

$\min_{y\in\mathbbm{R}^{d}}\sum_{i=1}^{n}\lambda_{i}\|y-x_{i}\|_{p}^{p}\leq\boldsymbol{\eta}\leq(1+\epsilon)\min_{y\in\mathbbm{R}^{d}}\sum_{i=1}^{n}\lambda_{i}\|y-x_{i}\|_{p}^{p}.$

Proof: The above is a simple consequence of Lemma 2.1. As in the proof of
Theorem 5, the algorithm utilizes the algorithm of Lemma 2.1 with the vector

$\displaystyle\sum_{i=1}^{m}\lambda_{i}^{1/p}\mathbf{S}_{i}x_{i}-\sum_{i=1}^{m}\sum_{j=1}^{m}\lambda_{i}\lambda_{j}^{1/p}\mathbf{S}_{j}x_{i}$
$\displaystyle=\sum_{i=1}^{m}\lambda_{i}^{1/p}\mathbf{S}_{i}\left(x_{i}-\sum_{j=1}^{m}\lambda_{j}x_{j}\right),$

which it can compute given access to $\mathbf{S}_{j}x_{i}\in\mathbbm{R}^{d}$
for all $i,j\in[m]$. The collection of vectors
$x_{i}-\sum_{j=1}^{m}\lambda_{j}x_{j}$ is centered with respect to
$\lambda_{1},\dots,\lambda_{m}$, and is exactly the vector that Lemma 2.1
expects to see.

###  3.2 Proof of Theorem 6

Let $M$ be an upper bound on the size of the coreset given from running the
merge-and-reduce framework for coresets in the streaming model [BS80] with $n$
points, where we use the coreset construction of Corollary 3.2. Notice that
$M=\mathrm{poly}(k\log(n)/\epsilon)$, where the extra $\mathrm{poly}(\log
n)$-factors come from the merge-and-reduce framework. For each $m\in[M]$, we
instantiate the distribution $\mathcal{S}_{m}$ given by Lemma 3.4 with
parameters $m,d$, $p$, $\epsilon$, and $\delta=k^{-M}/(100k)$. We sample
$\mathbf{S}^{(1)},\dots,\mathbf{S}^{(M)}$, where each
$\mathbf{S}^{(m)}=[\mathbf{S}_{1}^{(m)},\dots,\mathbf{S}_{m}^{(m)}]\sim\mathcal{S}_{m}$.
Note that for each $m\in[M]$, $\mathbf{S}^{(m)}$ is an $s\times(md)$ matrix
with $s=\mathrm{poly}(k\log(n)/\epsilon)$.

When a new vector $x_{i}\in\mathbbm{R}^{d}$ comes in the stream, we generate
the following:

  * •

An $\ell_{p}$-sketch of $x_{i}\in\mathbbm{R}^{d}$ for computing
$O(1)$-approximate $\ell_{p}$-distances which succeeds with probability
$1-o(1/n^{2})$ [Ind06]. This sketch is used to compute distances between
dataset points, as these are needed in order to compute the coreset of
Corollary 3.2.

  * •

For each $m\in[M]$ and $j\in[m]$, we store the vector
$\mathbf{S}_{j}^{(m)}x_{i}\in\mathbbm{R}^{s}$. This prepares the
$\ell_{p}^{p}$-median sketch of Lemma 3.4 to be evaluated on a cluster of size
$m$.

While processing the stream, we maintain a coreset $S\subset[n]$ with weights
$w\colon S\to\mathbbm{R}_{\geq 0}$ by using the merge-and-reduce framework of
Bentley and Saxe [BS80] for the streaming model. We maintain the
$\ell_{p}$-sketch of $x_{i}$ and the $O(M^{2})$ vectors
$\mathbf{S}_{j}^{(m)}x_{i}$ for all $m\in[M]$ and $j\in[m]$ whenever $i\in S$.
If $i\notin S$, we no longer maintain the $\ell_{p}$-sketch or the vectors
$\mathbf{S}_{j}^{(m)}x_{i}$. This ensures that the total space complexity is
$\mathrm{poly}(\log(nd),k,1/\epsilon)$ because we store: 1) $O(M^{3})$ vectors
in $\mathbbm{R}^{s}$, since there are $M$ points in the coreset, and each
$i\in S$ stores the $O(M^{2})$ vectors
$\mathbf{S}^{(m)}_{j}x_{i}\in\mathbbm{R}^{s}$, 2) the weights $w(i)$ for each
$i\in S$, and 3) the $\ell_{p}$-sketches for the points $x_{i}$.

When we finish processing the stream, we will evaluate the weighted cost of
every clustering of $S$ coreset points. For each clustering, and for each
cluster in the clustering, we will utilize the (weighted)
$\ell_{p}^{p}$-median sketch in order to estimate its cost. First, notice that
it suffices to compute the cost of all clustering of the coreset:

$\displaystyle\min_{\begin{subarray}{c}C_{1},\dots,C_{k}\\\ \text{partition
$[n]$}\end{subarray}}\sum_{\ell=1}^{k}\min_{c_{\ell}\in\mathbbm{R}^{d}}\sum_{i\in
C_{\ell}}\|x_{i}-c_{\ell}\|_{p}^{p}$
$\displaystyle=\min_{c_{1},\dots,c_{k}\in\mathbbm{R}^{d}}\sum_{i=1}^{n}\min_{\ell\in[k]}\|x_{i}-c_{\ell}\|_{p}^{p}$
$\displaystyle\approx_{1\pm\epsilon}\min_{c_{1},\dots,c_{k}\in\mathbbm{R}^{d}}\sum_{i\in
S}w(i)\min_{\ell\in[k]}\|x_{i}-c_{\ell}\|_{p}^{p}\quad\text{(def.~{}\ref{def:strong-
coreset})}$ $\displaystyle=\min_{\begin{subarray}{c}C_{1},\dots,C_{k}\\\
\text{partition
$S$}\end{subarray}}\sum_{\ell=1}^{k}\min_{c_{\ell}\in\mathbbm{R}^{d}}\sum_{i\in
C_{\ell}}w(i)\|x_{i}-c_{\ell}\|_{p}^{p}.$

Consider a fixed partition $C_{1},\dots,C_{k}$ of $S$. We evaluate the
weighted cost of clustering $S$ with $C_{1},\dots,C_{k}$ by evaluating the
cost of each cluster $C_{\ell}$ for $\ell\in[k]$ and summing these costs up.
Consider $C_{\ell}=\\{i_{1},\dots,i_{m}\\}\subset S$, and for $t\in[m]$, let

$\lambda_{t}=\frac{w(i_{t})}{\sum_{j=1}^{t}w(i_{j})}.$

Recall that we stored the vectors
$\mathbf{S}^{(m)}_{j}x_{i_{t}}\in\mathbbm{R}^{s}$, for all $j\in[m]$, so that
we may utilize the algorithm of Lemma 3.4 to obtain an
$(1\pm\epsilon)$-approximation to

$\boldsymbol{\eta}\approx_{1\pm\epsilon}\min_{c_{\ell}\in\mathbbm{R}^{d}}\sum_{t=1}^{m}\lambda_{t}\|x_{i_{t}}-c_{\ell}\|_{p}^{p}.$

Then, $\boldsymbol{\eta}\sum_{t=1}^{m}w(i_{t})$ gives us a $(1\pm
O(\epsilon))$-approximation to the cost of the cluster $C_{\ell}$. Since the
$\ell_{p}^{p}$-median sketch succeeds with very high probability, we may union
bound over all evaluations of at most $k$ cluster costs, for all $O(k^{M})$
possible clusterings of $S$. Hence, outputting the smallest cost gives us the
desired approximation.

##  Appendix A Proof of (9)

Here, we show that, for $y\in\mathbbm{R}^{d}$ being the minimizer of
$\sum_{i=1}^{n}\lambda_{I}\|y-x_{i}\|_{p}^{p}$ and
$\sum_{i=1}^{n}\lambda_{i}x_{i}=0$, the quantity
$\sum_{i=1}^{n}\lambda_{i}|y_{j}-x_{ij}|^{p}$ is up to a constant factor, the
same as $\sum_{i=1}^{n}\lambda_{i}|x_{ij}|^{p}$. Formally, for every
$j\in[d]$,

$\displaystyle\sum_{i=1}^{n}\lambda_{i}|y_{j}-x_{ij}|^{p}$
$\displaystyle\leq\sum_{i=1}^{n}\lambda_{i}|x_{ij}|^{p}\quad\text{(cost of
$y_{j}$ smaller than $0$)}$
$\displaystyle=\sum_{i=1}^{n}\lambda_{i}\left|x_{ij}-\sum_{i^{\prime}=1}^{n}\lambda_{i^{\prime}}x_{i^{\prime}j}\right|^{p}\quad\text{($x_{1j},\dots,x_{nj}$
are centered)}$
$\displaystyle\leq\sum_{i=1}^{n}\lambda_{i}\left(\sum_{i^{\prime}=1}^{n}\lambda_{i^{\prime}}|x_{ij}-x_{i^{\prime}j}|\right)^{p}\quad\text{(Triangle
Inequality)}$
$\displaystyle\leq\sum_{i=1}^{n}\sum_{i^{\prime}=1}^{n}\lambda_{i}\lambda_{i^{\prime}}|x_{ij}-x_{i^{\prime}j}|^{p},\quad\text{(Jensen's
Inequality)}$

Furthermore,

$\displaystyle\sum_{i=1}^{n}\sum_{i^{\prime}=1}^{n}\lambda_{i}\lambda_{i^{\prime}}|x_{ij}-x_{i^{\prime}j}|^{p}$
$\displaystyle=\sum_{i=1}^{n}\sum_{i^{\prime}=1}^{n}\lambda_{i}\lambda_{i^{\prime}}|x_{ij}-y_{j}+y_{j}-x_{i^{\prime}j}|^{p}$
$\displaystyle\leq\sum_{i=1}^{n}\sum_{i^{\prime}=1}^{n}\lambda_{i}\lambda_{i^{\prime}}\left(|y_{j}-x_{ij}|+|y_{j}-x_{i^{\prime}j}|\right)^{p}\quad\text{(Triangle
Inequality)}$
$\displaystyle\leq\sum_{i=1}^{n}\sum_{i^{\prime}=1}^{n}2^{p-1}\lambda_{i}\lambda_{i^{\prime}}\left(|y_{j}-x_{ij}|^{p}+|y_{j}-x_{i^{\prime}j}|^{p}\right)\quad\text{(H\"{o}lder
Inequality)}$
$\displaystyle=2^{p}\sum_{i=1}^{n}\lambda_{i}|y_{j}-x_{ij}|^{p}.$

##  Appendix B Sketching Medoid Costs

Another notion of centrality of a set of points is the medoid: for a metric
space $(X,d_{X})$ and a set of $n$ points $x_{1},\dots,x_{n}\in X$, the medoid
of the set of points is

$\mathop{\mathrm{argmin}}_{y\in\\{x_{1},\dots,x_{n}\\}}\sum_{i=1}^{n}d_{X}(x_{i},y).$

Analogously, we may define the $\ell_{p}^{p}$-medoid cost of a set of points
$x_{1},\dots,x_{n}\in\mathbbm{R}^{d}$ as

$\min_{y\in\\{x_{1},\dots,x_{n}\\}}\sum_{i=1}^{n}\|x_{i}-y\|_{p}^{p}.$

Similarly to the case of $\ell_{p}^{p}$-medians, the hope is to leverage
sketching algorithms for $\ell_{p}$-norms in order to sketch the
$\ell_{p}^{p}$-medoid cost. However, as we now show, the $\ell_{p}^{p}$-median
and $\ell_{p}^{p}$-medoid are very different from the sketching perspective.

######  Theorem 7 (Two-Pass Streaming Algorithm for $\ell_{p}^{p}$-Medoid).

Fix $n,d\in\mathbbm{N}$, as well as $p\in[1,2]$ and $\epsilon\in(0,1)$. There
exists a two-pass, insertion-only streaming algorithm using space
$\mathrm{poly}(\log(nd),1/\epsilon)$ which processes a set of $n$ points
$P=\\{x_{1},\dots,x_{n}\\}\in\mathbbm{R}^{d}$, and outputs a parameter
$\boldsymbol{\eta}\in\mathbbm{R}$ which satisfies

$\displaystyle\min_{z\in
P}\sum_{i=1}^{n}\|x_{i}-z\|_{p}^{p}\leq\boldsymbol{\eta}\leq(1+\epsilon)\min_{z\in
P}\sum_{i=1}^{n}\|x_{i}-z\|_{p}^{p}$

with probability at least $1-o(1)$.

The streaming algorithm will use the fact that there exists a linear sketch
for $\ell_{p}$ norms.

######  Theorem 8 (Sketching $\ell_{p}$ [Ind06]).

For any $m\in\mathbbm{N}$, $p\in[1,2]$, and $\epsilon,\delta\in(0,1)$, let
$s=O(\log(1/\delta)/\epsilon^{2})$. There exists a distribution $\mathcal{M}$
over $s\times m$ matrices, and an algorithm which takes as input a vector in
$\mathbbm{R}^{s}$. For any $x\in\mathbbm{R}^{m}$, with probability at least
$1-\delta$ over $\mathbf{S}\sim\mathcal{M}$, the algorithm on input
$\mathbf{S}x$, outputs $\boldsymbol{\eta}\in\mathbbm{R}_{\geq 0}$ satisfying

$\displaystyle(1-\epsilon)\|x\|_{p}\leq\boldsymbol{\eta}\leq(1+\epsilon)\|x\|_{p}.$

Proof: We let $\mathcal{M}$ be the distribution of Theorem 8 with accuracy
parameter $\epsilon$ and failure probability $o(1/n)$. The sketch proceeds in
the following way:

  * •

We sample $\mathbf{S}\sim\mathcal{M}$. In the first pass, we consider the
vector $x\in\mathbbm{R}^{nd}$ which stacks the $n$ vectors
$x_{1},\dots,x_{n}\in\mathbbm{R}^{d}$ and we maintain $\mathbf{S}x$.

  * •

In the second pass, we process the vectors
$x_{1},\dots,x_{n}\in\mathbbm{R}^{d}$. For $i=1,\dots,n$, we let
$y_{i}\in\mathbbm{R}^{nd}$ be the vector which vertically stacks $x_{i}$ for
$n$ repetitions. We query the algorithm on $\mathbf{S}(x-y_{i})$ to obtain an
estimate $\boldsymbol{\eta}_{i}$, and we maintain the minimum such
$\boldsymbol{\eta}_{i}$.

We note that by the amplification of the success probability to $o(1/n)$, we
may union bound over the $n$ queries to the streaming algorithm in the second
pass.The bounds on the value of $\boldsymbol{\eta}$ then follows from the fact
that $\boldsymbol{\eta}_{i}$ is an $(1\pm\epsilon)$-approximation to the
$\ell_{p}^{p}$-medoid cost which utilizes $x_{i}$ as the medoid.

######  Theorem 9.

[Lower Bound for One-Pass $\ell_{p}^{p}$-Medoids] Fix $p\in[1,\infty)$ and
$\epsilon\geq 10/n$. For large enough $n\in\mathbbm{N}$, a one-pass,
insertion-only streaming algorithm which processes a set $P$ of $n$ points in
$\\{0,1\\}^{2n}$ and outputs a number $\boldsymbol{\eta}\in\mathbbm{R}_{\geq
0}$ satisfying

$\min_{z\in P}\sum_{x\in
P}\|x-z\|_{p}^{p}\leq\boldsymbol{\eta}\leq(2-\epsilon)\cdot\min_{z\in
P}\sum_{x\in P}\|x-z\|_{p}^{p}$

with probability at least $9/10$ must use $\Omega(\epsilon n)$ bits of space.

In particular, this implies that any streaming algorithm which can approximate
the medoid cost in one pass and $\mathrm{polylog}(nd)$ space achieves
approximation no better than $2$.

The proof follows a reduction to one-way communication complexity of indexing.

######  Definition B.1 (Indexing Communication Problem).

The Indexing communication problem is a two-party one-way communication game
parametrized by $m\in\mathbbm{N}$.

  * •

Alice receives a bit-string $y\in\\{0,1\\}^{m}$, and Bob receives an index
$i\in[m]$.

  * •

Alice and Bob may use public-randomness, and Alice must produce a single
message to Bob so that Bob outputs $y_{i}$ with probability at least $9/10$.

######  Theorem 10 (One-way Communication of Indexing).

Any one-way communication protocol for the Indexing problem uses $\Omega(m)$
bits of communication.

Proof: [Proof of Theorem 9] We devise a protocol which utilizes a streaming
algorithm for approximating the medoid in order to solve the indexing problem.
Specifically, let $\alpha=\epsilon/10$ and suppose Alice receives a bit-string
$y\in\\{0,1\\}^{\alpha n}$, where $\alpha n\in\mathbbm{N}$. She considers the
subset

$P_{A}(y)=\\{e_{i}\in\\{0,1\\}^{2n}:y_{i}=1\\}.$

She process the stream of points $P_{A}(y)$ and communicates the contents of
the memory and the number of points in $k=P_{A}(y)$. Bob, upon seeing
$i\in[\alpha n]$ and $k\leq\alpha n$, needs to determine whether $e_{i}\in
P_{A}(y)$ or $e_{i}\notin P_{A}(y)$. Bob considers the subset

$P_{B}(i,k)=\left\\{e_{i}+e_{j}:j\in\\{\alpha n+1,\dots,\alpha
n+n-k\\}\right\\},$

so that $|P_{B}(i,k)|+|P_{A}(y)|=n$. Bob inserts the points $P_{B}(i,k)$ into
the stream, which results in the subset of points $P=P_{B}(i,k)\cup P_{A}(y)$
of $n$ points in $\\{0,1\\}^{2n}$. The crucial point is that Bob may
distinguish the two cases by considering the estimate $\boldsymbol{\eta}$ of
the medoid cost that the streaming algorithm produces. In particular,

  * •

If $e_{i}\in P_{A}(y)$, then, by choosing a candidate center
$z^{\prime}=e_{i}$, we have

$\min_{z\in P}\sum_{x\in P}\|x-z\|_{p}^{p}\leq\sum_{x\in
P}\|e_{i}-x\|_{p}^{p}\leq(k-1)\cdot 2+(n-k)\leq(1+2\alpha)n.$

In this case, the algorithm produces an estimate $\boldsymbol{\eta}$ which is
at most $(1+2\alpha)(2-\epsilon)n<2n-2$.

  * •

On the other hand, if $e_{i}\notin P_{A}(y)$, then every two points non-equal
$x,x^{\prime}\in P$ satisfy $\|x-x^{\prime}\|_{p}^{p}\geq 2$, so that

$\min_{z\in P}\sum_{x\in P}\|x-z\|_{p}^{p}\geq 2\cdot(n-1).$

Specifically, if the streaming algorithm was able to return a
$(2-\epsilon)$-approximation, Bob could distinguish between the two cases.
Since Bob communicates $O(\log_{2}(\epsilon n))\ll\epsilon n$ bits to encode
the number $k$, the space complexity of the streaming algorithm must be
$\Omega(\epsilon n)$.

## References

  * [ABIW09] Alexandr Andoni, Khanh Do Ba, Piotr Indyk, and David Woodruff.  Efficient sketches for earth-mover distance, with applications.  In Proceedings of the 50th Annual IEEE Symposium on Foundations of Computer Science (FOCS ’2009), 2009.
  * [ACKS15] Pranjal Awasthi, Moses Charikar, Ravishankar Krishnaswamy, and Ali Kemal Sinop.  The hardness of approximation of euclidean k-means.  In 31st International Symposium on Computational Geometry (SoCG 2015). Schloss Dagstuhl-Leibniz-Zentrum fuer Informatik, 2015.
  * [ACNN11] Alexandr Andoni, Moses Charikar, Ofer Neiman, and Huy L. Nguyen.  Near linear lower bound for dimension reduction in l1.  In Proceedings of the 52nd Annual IEEE Symposium on Foundations of Computer Science (FOCS ’2011), 2011.
  * [AHPV05] Pankaj K. Agarwal, Sariel Har-Peled, and Kasturi R. Varadarajan.  Geometric approximation via coresets.  Combinatorial and computational geometry, 2005.
  * [AIK08] Alexandr Andoni, Piotr Indyk, and Robert Krauthgamer.  Earth mover distance over high-dimensional spaces.  In Proceedings of the 19th ACM-SIAM Symposium on Discrete Algorithms (SODA ’2008), pages 343–352, 2008.
  * [AKO10] Alexandr Andoni, Robert Krauthgamer, and Krzysztof Onak.  Polylogarithmic approximation for edit distance and the asymmetric query complexity.  In Proceedings of the 51st Annual IEEE Symposium on Foundations of Computer Science (FOCS ’2010), 2010.
  * [AKO11] Alexandr Andoni, Robert Krauthgamer, and Krzysztof Onak.  Streaming algorithms from precision sampling.  In Proceedings of the 52nd Annual IEEE Symposium on Foundations of Computer Science (FOCS ’2011), 2011.
  * [AKR15] Alexandr Andoni, Robert Krauthgamer, and Ilya Razenshteyn.  Sketching and embedding are equivalent for norms.  In Proceedings of the 47th ACM Symposium on the Theory of Computing (STOC ’2015), pages 479–488, 2015.  Available as arXiv:1411.2577.
  * [AMS99] Noga Alon, Yossi Matias, and Mario Szegedy.  The space complexity of approximating the frequency moments.  Journal of Computer and System Sciences, 58(1):137–147, 1999.
  * [BBCA+19] Luca Becchetti, Marc Bury, Vincent Cohen-Addad, Fabrizio Grandoni, and Chris Schwiegelshohn.  Oblivious dimension reduction for $k$-means: Beyond subspaces and the johnson-lindenstrauss lemma.  In Proceedings of the 51th ACM Symposium on the Theory of Computing (STOC ’2019), 2019.
  * [BBCY17] Jaroslaw Blasiok, Vladimir Braverman, Stephen R. Chestnut, and Robert Krauthgamerand Lin F. Yang.  Streaming symmetric norms via measure concentration.  In Proceedings of the 50th ACM Symposium on the Theory of Computing (STOC ’2017), 2017.
  * [BC05] Bo Brinkman and Moses Charikar.  On the impossibility of dimension reduction in l1.  Journal of the ACM, 52(5):766–788, 2005.
  * [Beh21] Soheil Behnezhad.  Time-optimal sublinear algorithms for matching and vertex cover.  In Proceedings of the 62nd Annual IEEE Symposium on Foundations of Computer Science (FOCS ’2021), 2021.
  * [BFL16] Vladimir Braverman, Dan Feldman, and Harry Lang.  New frameworks for offline and streaming coreset constructions.  arXiv preprint arXiv:1612.00889, 2016.
  * [BFL+17] Vladimir Braverman, Gereon Frahling, Harry Lang, Christian Sohler, and Lin F Yang.  Clustering high dimensional dynamic data streams.  In Proceedings of the 34th International Conference on Machine Learning (ICML ’2017), 2017.
  * [BGLT20] Guy Blanc, Neha Gupta, Jane Lange, and Li-Yang Tan.  Estimating decision tree learnability with polylogarithmic sample complexity.  In Proceedings of Advances in Neural Information Processing Systems 33 (NeurIPS ’2020), 2020.
  * [BHPI02] Mihai Badoiu, Sariel Har-Peled, and Piotr Indyk.  Approximate clustering via core-sets.  In Proceedings of the 34th ACM Symposium on the Theory of Computing (STOC ’2002), 2002.
  * [BI14] Arturs Bačkurs and Piotr Indyk.  Better embeddings for planar earth-mover distance over sparse sets.  In Proceedings of the 41st International Colloquium on Automata, Languages and Programming (ICALP ’2014), 2014.
  * [BIRW16] Arturs Backurs, Piotr Indyk, Ilya Razenshteyn, and David P. Woodruff.  Nearly-optimal bounds for sparse recovery in generic norms, with applications to $k$-median sketching.  In Proceedings of the 27th ACM-SIAM Symposium on Discrete Algorithms (SODA ’2016), pages 318–337, 2016.  Available as arXiv:1504.01076.
  * [BS80] Jon Louis Bentley and James B Saxe.  Decomposable searching problems i. static-to-dynamic transformation.  Journal of Algorithms, 1(4):301–358, 1980.
  * [BYJKS04] Ziv Bar-Yossef, T.S. Jayram, Ravi Kumar, and D. Sivakumar.  An information statistics approach to data stream and communication complexity.  Journal of Computer and System Sciences, 68(4):702–732, 2004.
  * [BZD10] Christos Boutsidis, Anastasios Zouzias, and Petros Drineas.  Random projections for $k$-means clustering.  In Proceedings of Advances in Neural Information Processing Systems 23 (NeurIPS ’2010), 2010.
  * [CAK19] Vincent Cohen-Addad and CS Karthik.  Inapproximability of clustering in lp metrics.  In 2019 IEEE 60th Annual Symposium on Foundations of Computer Science (FOCS), pages 519–539. IEEE, 2019.
  * [CASS21] Vincent Cohen-Addad, David Saulpic, and Chris Schwiegelshohn.  A new coreset framework for clustering.  In Proceedings of the 53rd ACM Symposium on the Theory of Computing (STOC ’2021), 2021.
  * [CCF04] Moses Charikar, Kevin Chen, and Martin Farach-Colton.  Finding frequent items in data streams.  Theoretical Computer Science, 312(1):3–15, 2004.
  * [CEF+05] Artur Czumaj, Funda Engün, Lance Fortnow, Avner Magen, Ilan Newman, Ronitt Rubinfeld, and Christian Sohler.  Approximating the weight of the euclidean minimum spanning tree in sublinear time.  SIAM Journal on Computing, 2005.
  * [CEM+15] Michael B. Cohen, Sam Elder, Cameron Musco, Christopher Musco, and Mădălina Persu.  Dimensionality reduction for $k$-means clustering and low rank approximation.  In Proceedings of the 47th ACM Symposium on the Theory of Computing (STOC ’2015), 2015.
  * [Che09] Ke Chen.  On coresets for k-median and k-means clustering in metric and euclidean spaces and their applications.  SIAM Journal on Computing, 39(3):923–947, 2009.
  * [CJK+22] Artur Czumaj, Shaofeng H.-C. Jiang, Robert Krauthgamer, Pavel Veselý, and Mingwei Yang.  Streaming facility location in high dimension via new geometric hashing.  arXiv preprint arXiv:2204.02095, 2022.
  * [CJLW22] Xi Chen, Rajesh Jayaram, Amit Levi, and Erik Waingarten.  New streaming algorithms for high dimensional emd and mst.  In Proceedings of the 54th ACM Symposium on the Theory of Computing (STOC ’2022), 2022.
  * [CKK20] Yu Chen, Sampath Kannan, and Sanjeev Khanna.  Sublinear algorithms and lower bounds for metric tsp cost estimation.  In 47th International Colloquium on Automata, Languages, and Programming (ICALP 2020). Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2020.
  * [CM05] Graham Cormode and Shan Muthukrishnan.  An improved data stream summary: the count-min sketch and its applications.  Journal of Algorithms, 2005.
  * [CS09] Artur Czumaj and Chirstian Sohler.  Estimating the weight of metric minimum spanning trees in sublinear time.  SIAM Journal on Computing, 2009.
  * [FL11] Dan Feldman and Michael Langberg.  A unified framework for approximating and clustering data.  In Proceedings of the 43rd ACM Symposium on the Theory of Computing (STOC ’2011), 2011.
  * [FSS13] Dan Feldman, Melanie Schmidt, and Christian Sohler.  Turning big data into tiny data: constant-size coresets for k-means, pca and projective clustering.  In Proceedings of the 24th ACM-SIAM Symposium on Discrete Algorithms (SODA ’2013), 2013.
  * [HK20] Monika Henzinger and Sagar Kale.  Fully-dynamic coresets.  In 28th Annual European Symposium on Algorithms (ESA 2020). Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2020.
  * [HPM04] Sariel Har-Peled and Soham Mazumdar.  On coresets for $k$-means and $k$-median clustering.  In Proceedings of the 36th ACM Symposium on the Theory of Computing (STOC ’2004), 2004.
  * [HSYZ19] Wei Hu, Zhao Song, Lin F. Yang, and Peilin Zhong.  Nearly optimal dynamic $k$-means clustering for high-dimensional data.  arXiv preprint arXiv:1802.00459, 2019.
  * [HV20] Lingxiao Huang and Nisheeth K. Vishnoi.  Coresets for clustering in euclidean spaces: importance sampling is nearly optimal.  In Proceedings of the 52nd ACM Symposium on the Theory of Computing (STOC ’2020), 2020.
  * [Ind04] Piotr Indyk.  Algorithms for dynamic geometric problems over data streams.  In Proceedings of the 36th ACM Symposium on the Theory of Computing (STOC ’2004), 2004.
  * [Ind06] Piotr Indyk.  Stable distributions, pseudorandom generators, embeddings, and data stream computation.  Journal of the ACM, 53(3):307–323, 2006.
  * [IW05] Piotr Indyk and David Woodruff.  Optimal approximations of the frequency moments of data streams.  In Proceedings of the 37th ACM Symposium on the Theory of Computing (STOC ’2005), pages 202–208, 2005.
  * [JL84] William Johnson and Joram Lindenstrauss.  Extensions of Lipschitz mappings into a Hilbert space.  In Conference in modern analysis and probability (New Haven, Conn., 1982), volume 26 of Contemporary Mathematics, pages 189–206. 1984\.
  * [JW09] Thathachar S. Jayram and David Woodruff.  The data stream complexity of cascaded norms.  In Proceedings of the 50th Annual IEEE Symposium on Foundations of Computer Science (FOCS ’2009), 2009.
  * [JW18] Rajesh Jayaram and David Woodruff.  Perfect lp sampling in a data stream.  In Proceedings of the 59th Annual IEEE Symposium on Foundations of Computer Science (FOCS ’2018), 2018.
  * [KBV20] Weihao Kong, Emma Brunskill, and Gregory Valiant.  Sublinear optimal policy value estimation in contextual bandits.  In Proceedings of the 23rd International Conference on Artificial Intelligence and Statistics (AISTATS ’2020), 2020.
  * [KMNFT20] Michael Kapralov, Slobodan Mitrović, Ashkan Norouzi-Fard, and Jakab Tardos.  Space efficient approximation to maximum matching size from uniform edge samples.  In Proceedings of the Fourteenth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 1753–1772. SIAM, 2020.
  * [KNW10] Daniel M. Kane, Jelani Nelson, and David P. Woodruff.  On the exact space complexity of sketching and streaming small norms.  In Proceedings of the 21st ACM-SIAM Symposium on Discrete Algorithms (SODA ’2010), 2010.
  * [KOR00] Eyal Kushilevitz, Rafail Ostrovsky, and Yuval Rabani.  Efficient search for approximate nearest neighbor in high dimensional spaces.  SIAM Journal on Computing, 30(2):457–474, 2000.
  * [KSS04] Amit Kumar, Yogish Sabharwal, and Sandeep Sen.  A simple linear time $(1+\epsilon)$-approximation algorithm for $k$-means clustering in any dimension.  In Proceedings of the 45th Annual IEEE Symposium on Foundations of Computer Science (FOCS ’2004), 2004.
  * [KV18] Weihao Kong and Gregory Valiant.  Estimating learnability in the sublinear data regime.  In Proceedings of Advances in Neural Information Processing Systems 31 (NeurIPS ’2018), pages 5455–5464, 2018.
  * [LN04] James R. Lee and Assaf Naor.  Embedding the diamond graph in lp and dimension reduction in l1.  Geometric and Functional Analysis, 14(4):745–747, 2004.
  * [LS10] Michael Langberg and Leonard J. Schulman.  Universal epsilon-approximators for integrals.  In Proceedings of the 21st ACM-SIAM Symposium on Discrete Algorithms (SODA ’2010), 2010.
  * [LSW17] Euiwoong Lee, Melanie Schmidt, and John Wright.  Improved and simplified inapproximability for k-means.  Information Processing Letters, 120:40–43, 2017.
  * [MMR19] Konstantin Makarychev, Yuri Makarychev, and Ilya Razenshteyn.  Performance of johnson-lindenstrauss transform for $k$-means and $k$-medians clustering.  In Proceedings of the 51th ACM Symposium on the Theory of Computing (STOC ’2019), 2019.
  * [NO08] Huy N Nguyen and Krzysztof Onak.  Constant-time approximation algorithms via local improvements.  In 2008 49th Annual IEEE Symposium on Foundations of Computer Science, pages 327–336. IEEE, 2008.
  * [ORRR12] Krzysztof Onak, Dana Ron, Michal Rosen, and Ronitt Rubinfeld.  A near-optimal sublinear-time algorithm for approximating the minimum vertex cover size.  In Proceedings of the twenty-third annual ACM-SIAM symposium on Discrete Algorithms, pages 1123–1131. SIAM, 2012.
  * [Owe13] Art B. Owen.  Monte Carlo theory, methods, and examples.  2013\.
  * [PR07] Michal Parnas and Dana Ron.  Approximating the minimum vertex cover in sublinear time and a connection to distributed algorithms.  Theoretical Computer Science, 381(1-3):183–196, 2007.
  * [SS02] Michael Saks and Xiaodong Sun.  Space lower bounds for distance approximation in the data stream model.  In Proceedings of the 34th ACM Symposium on the Theory of Computing (STOC ’2002), 2002.
  * [SW18] Christian Sohler and David Woodruff.  Strong coresets for $k$-median and subspace approximation.  In Proceedings of the 59th Annual IEEE Symposium on Foundations of Computer Science (FOCS ’2018), 2018.
  * [YYI09] Yuichi Yoshida, Masaki Yamamoto, and Hiro Ito.  An improved constant-time approximation algorithm for maximum~ matchings.  In Proceedings of the forty-first annual ACM symposium on Theory of computing, pages 225–234, 2009.