# Streaming Euclidean $k$-median and $k$-means with $o(\log n)$ Space

Vincent Cohen-Addad Google Research. E-mail<EMAIL_ADDRESS>David P.
Woodruff Carnegie Mellon University. E-mail<EMAIL_ADDRESS>Supported in part by a Simons Investigator Award.  Samson Zhou Texas A&M
University. E-mail<EMAIL_ADDRESS>Work done in part while at Rice
University and UC Berkeley.

###### Abstract

We consider the classic Euclidean $k$-median and $k$-means objective on data
streams, where the goal is to provide a $(1+\varepsilon)$-approximation to the
optimal $k$-median or $k$-means solution, while using as little memory as
possible. Over the last 20 years, clustering in data streams has received a
tremendous amount of attention and has been the test-bed for a large variety
of new techniques, including coresets, the merge-and-reduce framework,
bicriteria approximation, sensitivity sampling, and so on. Despite this
intense effort to obtain smaller sketches for these problems, all known
techniques require storing at least $\Omega(\log(n\Delta))$ words of memory,
where $n$ is the size of the input and $\Delta$ is the aspect ratio. A natural
question is if one can beat this logarithmic dependence on $n$ and $\Delta$.
In this paper, we break this barrier by first giving an insertion-only
streaming algorithm that achieves a $(1+\varepsilon)$-approximation to the
more general $(k,z)$-clustering problem, using
$\tilde{\mathcal{O}}\left(\frac{dk}{\varepsilon^{2}}\right)\cdot(2^{z\log
z})\cdot\min\left(\frac{1}{\varepsilon^{z}},k\right)\cdot\text{poly}(\log\log(n\Delta))$
words of memory. Our techniques can also be used to achieve two-pass
algorithms for $k$-median and $k$-means clustering on dynamic streams using
$\tilde{\mathcal{O}}\left(\frac{1}{\varepsilon^{2}}\right)\cdot\text{poly}(d,k,\log\log(n\Delta))$
words of memory.

##  1 Introduction

Clustering is the problem of partitioning an input dataset to identify and
extract important underlying structural information, so that ideally, points
in the same cluster have similar properties and points in different clusters
have dissimilar properties. Various formulations of clustering have become
important cornerstones in combinatorial optimization, computational geometry,
computer vision, data science, and machine learning, while classic
formulations of clustering such as the $k$-median and $k$-means problems have
been studied since the 1950s [S+56, Mac67].

More generally, $k$-median and $k$-means are specific parameterizations of the
$(k,z)$-clustering problem, which is one of the most commonly studied center-
based clustering objectives, where the input data elements belong to a metric
space, every cluster is induced by a corresponding center point, and the cost
of each data element is a fixed function of the distance between the data
point and the corresponding cluster center. Formally, the input is a set $X$
of $n$ points $x_{1},\ldots,x_{n}$ along with a distance function dist, a
cluster parameter $k>0$, and a positive integer exponent $z>0$. The
$(k,z)$-clustering objective is defined to be

$\min_{C\subset\mathbb{R}^{d},|C|=k}\sum_{i=1}^{n}\min_{c\in
C}{\text{dist}}(x_{i},c)^{z}.$

When $z=1$ and $z=2$, the problem is known as $k$-median clustering and
$k$-means clustering, respectively.

#### Clustering in the streaming model.

As modern datasets have significantly increased in size, often consisting of
hundreds of millions of points, attention has shifted to large-scale
computational models that require neither the storage of the dataset nor
multiple passes over the data, both of which could require expensive memory or
runtime usage. Perhaps the simplest of these models is the insertion-only
streaming model, where in our setting, the points $x_{1},\ldots,x_{n}$ of $X$
arrive sequentially, and the goal is to output a clustering whose cost is a
$(1+\varepsilon)$-multiplicative approximation of the optimal clustering of
$X$ while using space sublinear in $n$, for some input accuracy parameter
$\varepsilon>0$. Although it is not possible to explicitly output the label
for each point in the dataset using sublinear space, streaming algorithms
often instead implicitly define the clusters by outputting the center of each
cluster, so that each point is assigned to the cluster with the closest
center.

Most streaming algorithms for $(k,z)$-clustering follow the same template.
Over the course of the stream, they first construct a coreset, which is a
small weighted subset of points from the original dataset (see Definition 1.10
for a formal definition). After the stream terminates, they then compute an
optimal or near-optimal $(k,z)$-clustering of the coreset to produce the
cluster centers.

In an offline setting, where the input dataset $X$ is fully accessible at all
times, the best known coreset constructions sample and store
$\tilde{\mathcal{O}}\left(\frac{k}{\varepsilon^{2}}\right)\cdot\min\left(k,\frac{1}{\varepsilon^{z}}\right)$111Here
we use
$\tilde{\mathcal{O}}\left(f\left(n,d,k,\Delta,\frac{1}{\varepsilon}\right)\right)$
to denote
$\mathcal{O}\left(f\left(n,d,k,\Delta,\frac{1}{\varepsilon}\right)\right)\cdot\operatorname*{polylog}\left(f\left(n,d,k,\Delta,\frac{1}{\varepsilon}\right)\right)$.
weighted points of the input dataset [CLSS22]. Although the sampled points
could have weight $\operatorname*{poly}(n,k,d,\Delta)$, where $\Delta$ is the
aspect ratio, these weights can be represented by $\mathcal{O}\left(1\right)$
words of space, where a single word is $\Theta(\log(nd\Delta))$ bits of space,
i.e., the space necessary to represent a single coordinate of an input point
along with the corresponding weight. Thus it follows that there exists a
succinct representation of the input dataset with size _independent_ of the
number $n$ of input points.

Therefore, a natural question to ask is:

> _Is it possible to perform approximate $(k,z)$-clustering in the streaming
> model using $o_{k,d,\varepsilon}(\log n)$ words of memory?_

More specifically, the best offline coreset constructions use
$\tilde{\mathcal{O}}\left(\frac{kd}{\varepsilon^{2}}\min\left(k,\frac{1}{\varepsilon^{z}}\right)\log(n\Delta)\right)$
bits of space to store the weighted subset of the input points.

> _Can $(k,z)$-clustering be performed in the streaming model with the same
> space complexity?_

At first glance, the answer might seem unlikely. The simplest and perhaps most
standard approach for insertion-only streams is the merge-and-reduce
framework, which first partitions the stream into a number of blocks, creating
a coreset on-the-fly for each block. The framework can be viewed as a binary
tree, where each node contains a coreset of the merger of the coresets
contained at the children nodes, and the root node is a coreset for the entire
stream. However, the main pitfall is that because each step incurs distortion
through the construction of a new coreset, the coresets for each block require
accuracy $\left(1+\mathcal{O}\left(\frac{\varepsilon}{\log n}\right)\right)$,
which results in space
$\tilde{\mathcal{O}}\left(\frac{kd}{\varepsilon^{2}}\min\left(k,\frac{1}{\varepsilon^{z}}\right)\log(n\Delta)\log^{3}n\right)$.
Moreover, all known streaming algorithms implicitly leverage analysis that
requires correctness at all times of the stream, thereby incurring
polylogarithmic dependencies in the number $n$ of input points.

On the other hand, we only ask for correctness at the end of the stream, so we
do not inherently require correctness at all times. Moreover, there exist
special structured problems in the streaming model where correctness at all
times does not require additional polylogarithmic overhead [BCIW16, BDN17,
BCI+17, Bla18].

###  1.1 Our Contributions

In this paper, we answer the question in the affirmative for insertion-only
streams and in the negative for insertion-deletion streams. As a warm-up, we
first give a simple algorithm for $(k,z)$-clustering on insertion-only streams
that uses $o(\log n)$ words of space and nearly matches the same space
complexity as the best offline coreset constructions.

######  Theorem 1.1.

Given an accuracy parameter $\varepsilon\in(0,1)$, an integer $k>0$ for the
number of clusters, a constant $z\geq 1$, and a data stream consisting of
points $X\subseteq[\Delta]^{d}$ with $X=\\{x_{1},\ldots,x_{n}\\}$, there
exists a one-pass streaming algorithm that uses
$\tilde{\mathcal{O}}\left(\frac{dk}{\varepsilon^{2}}\right)\cdot(2^{z\log
z})\cdot\min\left(\frac{1}{\varepsilon^{z}},k\right)\cdot\operatorname*{poly}(\log\log(n\Delta))$
words of memory and outputs a $(1+\varepsilon)$-approximation to
$(k,z)$-clustering at all times in the stream.

Theorem 1.1 is the first result to achieve $(k,z)$-clustering on insertion-
only streams using $o(\log n)$ words of space. That is, even for constant $d$,
$k$, $\varepsilon$, and $z$, existing works use $\operatorname*{poly}(\log n)$
words of space, while Theorem 1.1 achieves $\operatorname*{poly}(\log\log n)$
words of space and in fact, its guarantees apply for all times in the stream.
Moreover, we match or improve previous algorithms, e.g., [HM04, HK07, Che09,
FL11, BFLR19], for insertion-only streams in terms of dependencies on $d$,
$k$, and $\frac{1}{\varepsilon}$ (see Figure 1 for a comprehensive summary).
Furthermore, we match the best known _offline_ coreset constructions [CLSS22]
up to lower-order terms in the bit complexity.

Despite its simplicity, we shall use Theorem 1.1 as an opportunity to present
several important structural properties that we will build upon.

#### Online sensitivity sampling.

First, we define the online sensitivities for $(k,z)$-clustering as follows:

######  Definition 1.2 (Online sensitivity for $(k,z)$-clustering).

Let $x_{1},\ldots,x_{n}$ be a sequence of points with metric dist and let
$X_{t}=x_{1},\ldots,x_{t}$ for all $t\in[n]$. The _online sensitivity_ of
$x_{t}$, denoted $\sigma_{t}$, is

$\sigma_{t}:=\max_{C:|C|\leq
k}\frac{\operatorname*{cost}(x_{t},C)}{\operatorname*{cost}(X_{t},C)}=\max_{C:|C|\leq
k}\frac{{\text{dist}}(x_{t},C)^{z}}{\sum_{i=1}^{t}{\text{dist}}(x_{i},C)^{z}}.$

Sampling elements of the input with probability proportional to approximations
of their sensitivities and online sensitivities, i.e., sensitivity sampling,
and its variants have commonly been used in many other applications [DDH+09,
VX12, CP15, CMM17, BDM+20, HV20, CSS21, MMM+22, CLSS22, MMM+23, WY23b],
because it is one of the most intuitive algorithms for acquiring a
representative subset of the input – each element of the input is sampled with
probability proportional to a quantity that informally measures how
“important” the element is.

The sample complexity of sensitivity sampling is proportional to the total
sensitivity, that is, the sum of the sensitivities of each element. Although
in some cases analyzing the total sensitivity can be quite involved [DDH+09,
VX12, CP15, CMM17, BDM+20, WY23b], we give a simple proof upper bounding the
sum of the online sensitivities:

######  Theorem 1.3.

(Upper bound on sum of online sensitivities) Let
$X=\\{x_{1},\ldots,x_{n}\\}\subset[\Delta]^{d}$ and for each point $x_{t}$
with $t\in[n]$, let $\sigma_{t}$ denote its online sensitivity for
$(k,z)$-clustering for any $z\geq 1$. Then
$\sum_{t=1}^{n}\sigma_{t}=\mathcal{O}\left(2^{2z}k\log^{2}(nd\Delta)\right).$

We leverage Theorem 1.3 to analyze the space complexity of online sensitivity
sampling in Theorem 3.5, which provides a separate but sub-optimal approach
for $(k,z)$-clustering on insertion-only streams. However, we utilize online
sensitivity sampling as a crucial subroutine en route to achieving our main
result in Theorem 1.1. We summarize these results in Figure 1.

Streaming algorithm | Words of Memory
---|---
[HK07], $z\in\\{1,2\\}$ | $\tilde{\mathcal{O}}\left(\frac{dk^{1+z}}{\varepsilon^{\mathcal{O}\left(d\right)}}\log^{d+z}n\right)$
[HM04], $z\in\\{1,2\\}$ | $\tilde{\mathcal{O}}\left(\frac{dk}{\varepsilon^{d}}\log^{2d+2}n\right)$
[Che09], $z\in\\{1,2\\}$ | $\tilde{\mathcal{O}}\left(\frac{d^{2}k^{2}}{\varepsilon^{2}}\log^{8}n\right)$
[FL11], $z\in\\{1,2\\}$ | $\tilde{\mathcal{O}}\left(\frac{d^{2}k}{\varepsilon^{2z}}\log^{1+2z}n\right)$
Sensitivity and rejection sampling [BFLR19] | $\tilde{\mathcal{O}}\left(\frac{d^{2}k^{2}}{\varepsilon^{2}}\log n\right)$
Online sensitivity sampling, i.e., Theorem 3.5 | $\tilde{\mathcal{O}}\left(\frac{d^{2}k^{2}}{\varepsilon^{2}}\log n\right)$
Merge-and-reduce with coreset of [CLSS22] | $\tilde{\mathcal{O}}\left(\frac{dk}{\varepsilon^{2}}\log^{4}n\right)\cdot\min\left(\frac{1}{\varepsilon^{z}},k\right)$
This work, i.e., Theorem 1.1 | $\tilde{\mathcal{O}}\left(\frac{dk}{\varepsilon^{2}}\right)\cdot\min\left(\frac{1}{\varepsilon^{z}},k\right)\cdot\text{poly}(\log\log n)$
Fig. 1: Table of $(k,z)$-clustering algorithms on insertion-only streams. We
summarize existing results with $z=\mathcal{O}\left(1\right)$ and
$\Delta=\operatorname*{poly}(n)$ for the purpose of presentation.

#### Upper bounds, lower bounds and separations for dynamic streams.

A natural follow-up question to ask is whether our results extend to dynamic
streams, where points may be inserted and deleted in updates to the stream. We
first show that a streaming algorithm that uses a single pass on a insertion-
deletion stream cannot even provide a $2$-approximation to the _cost_ of an
optimal $(k,z)$-clustering using $o(\log n)$ words of space. Note that by
comparison, our algorithm in Theorem 1.1 actually provides a strong coreset at
the end of an insertion-only stream and thus not only outputs a set of $k$
near-optimal centers, but also a
$(1+\mathcal{O}\left(\varepsilon\right))$-approximation to the cost induced by
those $k$ centers and by extension, a $(1+\varepsilon)$-approximation to the
optimal $(k,z)$-clustering cost. Since the following statement, Theorem 1.4,
shows that any algorithm that provides a $2$-approximation on dynamic streams
requires $\Omega(\log^{2}n)$ bits, i.e., $\Omega(\log n)$ words of space, then
we cannot expect to achieve such a result for dynamic streams.

######  Theorem 1.4.

Let $z$ be a constant. Then even for $k=1$, any algorithm that with
probability at least $\frac{2}{3}$, simultaneously outputs a $2$-approximation
to the optimal $(k,z)$-clustering cost at all times of a dynamic stream of
length $n$ for points in $\\{0,1\\}^{d}$ with $d=\Omega(\log n)$ must use
$\Omega(\log^{2}n)$ bits of space.

At first glance, Theorem 1.4 may not seem like a separation from Theorem 1.1
because the former requires the assumption that $d=\Omega(\log n)$ and the
latter has a linear dependence on $d$. However, we show that if the goal is
just to output an estimation to the _cost_ of the optimal $(k,z)$-clustering
at the end of an insertion-only stream, then our upper bounds do not require a
linear dependence on $d$:

######  Theorem 1.5.

Given an accuracy parameter $\varepsilon\in(0,1)$, an integer $k>0$ for the
number of clusters, a constant $z>0$, and a data stream consisting of points
$X\subseteq[\Delta]^{d}$ with $X=\\{x_{1},\ldots,x_{n}\\}$, there exists a
one-pass streaming algorithm that uses
$\tilde{\mathcal{O}}\left(\frac{k}{\varepsilon^{2}}\right)\min\left(\frac{1}{\varepsilon^{z}},k\right)\cdot\operatorname*{poly}(\log\log(n\Delta))$
words of memory and outputs a $(1+\varepsilon)$-approximation to the cost of
the optimal $(k,z)$-clustering at all times in the stream.

Combined with Theorem 1.4, Theorem 1.5 provides a separation for the
$(k,z)$-clustering problem between insertion-only and dynamic streams.

We also show that any one-pass dynamic streaming algorithm that computes a
constant factor-approximation to the clustering cost from a weighted sample of
the input points must use $\Omega(\log^{2}n)$ bits of space.

######  Theorem 1.6.

Any one-pass dynamic streaming algorithm that computes a $2$-approximation to
the clustering cost from a weighted sample of the input points with
probability at least $\frac{2}{3}$ must use $\Omega(\log^{2}n)$ bits of space.

On the other hand, we show that the $\Omega(\log^{2}n)$ space barrier can be
broken on dynamic streams if we permit algorithms an additional pass over the
stream.

######  Theorem 1.7.

There exists a two-pass dynamic streaming algorithm that outputs a
$(1+\varepsilon)$-coreset for $k$-median and $k$-means clustering, with
probability at least $\frac{2}{3}$, and uses
$\tilde{\mathcal{O}}\left(\frac{1}{\varepsilon^{2}}\right)\cdot\operatorname*{poly}(d,k,\log\log(n\Delta))$
words of space. We remark that Theorem 1.7 can be generalized to
$(k,z)$-clustering for all $z\in[1,2]$.

#### Wasserstein-$z$ bicriteria embedding.

Along the way, we introduce a quadtree embedding technique that achieves a
structured bicriteria approximation to $(k,z)$-clustering and is key to our
two-pass dynamic streaming algorithm. Although quadtree embeddings are popular
for $k$-median clustering in big data models, generalizations to
$(k,z)$-clustering are not known for a single quadtree and in fact, there are
simple examples on a line that show that the expected distortion between the
squared distances of a set of $n$ points, i.e., $z=2$, and the estimated
distance by a quadtree is $\Omega(n)$, e.g., see Appendix A. Recently,
[CLN+20] overcame this barrier in the offline setting by considering multiple
quadtrees and taking the minimum estimated distance for pairs of points across
the quadtrees, but this approach does not seem to work for the streaming
setting, because we will no longer be able to estimate these distances after
embedding the quadtree into $L_{1}$.

We instead develop a quadtree embedding technique for $(k,z)$-clustering that
uses a bicriteria approximation, which often suffices for downstream
applications. Indeed, we crucially use our embedding as a subroutine toward
our $(1+\varepsilon)$-approximate streaming algorithm for $(k,z)$-clustering
in Theorem 1.1. Our embedding has implications to the Wasserstein-$z$
distance, which is a distance between probability measures and corresponds to
the important earth mover distance for $z=1$. We use $\operatorname*{WASSZZ}$
to denote the $z$-th power of the Wasserstein-$z$ distance and formally recall
its definition in Section 1.6. Our embedding then has the following
guarantees:

######  Theorem 1.8.

Let $\mu,\nu\in\mathbb{R}^{[\Delta]^{d}}$ be probability measures such that
$\nu$ has support at most $k$ on a set $C\subset[\Delta]^{d}$. Then there
exists a quadtree embedding $W_{s}$ parameterized by a random shift parameter
$s$ and explicit mappings $\psi,\phi$ such that with probability at least
$0.99$, $\psi(\nu)$ is a probability mass with support at most
$\mathcal{O}\left(k\right)$ on the set $\phi(C)$ and
$\|W_{s}(\mu-\psi(\nu))\|_{1}\leq\mathcal{O}\left(d^{1+0.5z}\log^{z-1}\Delta\right)\cdot\operatorname*{WASSZZ}(\mu,\nu).$

In fact, the set $\phi(C)$ includes the set $C$, so $\psi(\nu)$ is a
probability measure that includes the support of $\nu$. Thus, our embedding
can be viewed as a sketch that produces a bicriteria approximation for the
Wasserstein-$z$ distance.

###  1.2 Technical Overview: Initial Challenges

In this section, we describe the intuition behind our main algorithm and how
it overcomes significant barriers for previous techniques.

#### Merge-and-reduce does not work.

A standard approach for $(k,z)$-clustering on datasets of $\mathbb{R}^{d}$ on
insertion-only streams is the merge-and-reduce framework, due to [BS80, HM04].
Given a coreset construction algorithm for $(k,z)$-clustering with size
$S(n,d,k,\varepsilon,\delta)$ where $\varepsilon$ is the desired accuracy and
$\delta$ is the failure probability, the merge-and-reduce framework first
partitions the stream into consecutive blocks of size
$S(n,d,k,\varepsilon^{\prime},\delta^{\prime})$, where
$\varepsilon^{\prime}=\frac{\varepsilon}{\mathcal{O}\left(\log n\right)}$ and
$\delta^{\prime}=\frac{\delta}{\operatorname*{poly}(n)}$. A coreset with
accuracy $(1+\varepsilon^{\prime})$ and failure probability $\delta^{\prime}$
is then computed for each block, so that each coreset also uses space
$S(n,d,k,\varepsilon^{\prime},\delta^{\prime})$. We can then view these
coresets as the leaves of a binary tree of height $\mathcal{O}\left(\log
n\right)$, where each node in the tree at depth $t$ denotes a coreset computed
from two coresets at depth $t+1$. See Figure 2 for an example of the merge-
and-reduce framework.

The coreset at the root node represents a coreset for the entire stream, and
since each node only requires the coresets of its children nodes, this process
can be done on-the-fly. Because there is a $(1+\varepsilon^{\prime})$
multiplicative loss of accuracy at each level, then the coreset at the root
node has accuracy $(1+\varepsilon^{\prime})^{\mathcal{O}\left(\log
n\right)}=(1+\varepsilon)$, as desired. However, all known coreset
constructions use $\Omega\left(\frac{1}{\varepsilon^{\prime}}\right)$ space,
and thus for $\varepsilon^{\prime}=\frac{\varepsilon}{\mathcal{O}\left(\log
n\right)}$, the merge-and-reduce approach would not work for our goal because
it incurs extraneous $\log n$ factors.

Stream:Depth 3:$C_{3,1}$$C_{3,2}$$C_{3,3}$$C_{3,4}$Depth
2:$C_{2,1}$$C_{2,2}$Depth 1:$C_{1,1}$ Fig. 2: Merge and reduce framework. Each
coreset $C_{t,i}$ at depth $t$ is a coreset of the merger of the two children
coresets $C_{t+1,2i-1}$ and $C_{t+1,2i}$ at depth $t+1$.

#### Quadtrees are not enough.

Another approach for $(k,z)$-clustering on dynamic streams is placing the
input space into a quadtree [Che09], which is a tree of randomly shifted grids
on $[\Delta]^{d}$. In the quadtree, the coarsest grid contains a single cell
representing the entire input space of $[\Delta]^{d}$ and each subsequent
level of the quadtree corresponds to a refinement of the previous grid with
smaller grid cells, e.g., by splitting each grid cell in the previous level
into $2^{d}$ smaller grid cells. It then suffices to identify the finest level
that contains a number of non-empty cells below a certain threshold, as well
as the non-empty cells and the number of points in those non-empty cells at
that level, e.g., using sparse recovery [BFL+17, SYZ18]. However, because the
number of cells in the quadtree is $2^{\mathcal{O}\left(d\right)}$, then the
frequency vector of the number of points in each cell has dimension
$2^{\mathcal{O}\left(d\right)}$ and thus recovery of the non-empty cells and
the number of points in those cells often requires $\mathcal{O}\left(d\right)$
space, which, combined with subsequent union bounds over nets with size
$\left(\frac{n\Delta}{\varepsilon}\right)^{\mathcal{O}\left(kd\right)}$,
results in $\operatorname*{poly}(d,\log(n\Delta))$ dependencies prohibitive
for our goal and in fact even the $\operatorname*{poly}(d)$ dependencies
cannot be removed through known dimensionality reduction techniques [MMR19,
ISZ21, CW22b, CW22a] without incurring extraneous
$\operatorname*{poly}\left(\frac{1}{\varepsilon}\right)$ factors.

#### Offline sensitivity sampling is not enough.

Another approach would be to adapt recent methods for offline coreset
constructions to the streaming model. Unfortunately, many of the optimal or
near-optimal constructions [HV20, CSS21, CLSS22] use multiple stages of
sampling, which seems inherently infeasible for one-pass algorithms in the
streaming model.

Luckily, the sensitivity framework [FL11, FS12, BFL+21] for coreset
construction seems amenable to adaptation to data streams. The (offline)
_sensitivity_ of each point $x_{1},\ldots,x_{n}\in\mathbb{R}^{d}$ of a set $X$
is a quantity that informally measures the “importance” of that point.
Formally, the sensitivity of $x_{t}$ for the $(k,z)$-clustering problem is
defined as

$\max_{C\subset\mathbb{R}^{d}:|C|\leq
k}\frac{\operatorname*{cost}(x_{t},C)}{\operatorname*{cost}(X,C)}=\max_{C\subset\mathbb{R}^{d}:|C|\leq
k}\frac{{\text{dist}}(x_{t},C)^{z}}{\sum_{t=1}^{n}{\text{dist}}(x_{t},C)^{z}}.$

Traditionally, the sensitivity sampling framework independently samples a
fixed number of points with replacement, so that each point is sampled with
probability proportional to (some approximation to) its sensitivity. However,
it can be shown that sampling each point independently without replacement, so
that the total number of sampled points is a random variable, is also a valid
coreset construction algorithm. Recent works have considered online variants
of sensitivity sampling for other problems such as subspace embedding and low-
rank approximation [CMP20, BDM+20, WY23b].

#### Online sensitivity sampling.

We thus consider online sensitivity sampling for $(k,z)$-clustering. The
argument of correctness is relatively straightforward. We first fix a set $C$
of $k$ centers and show that the expectation of the cost of clustering with
$C$ for the sampled points is an unbiased estimator of the clustering cost of
$C$ with respect to the input set $X$. We then upper bound the variance of the
cost of clustering the coreset with $C$, which allows us to apply a standard
martingale argument that shows concentration, i.e., the coreset approximately
preserves the clustering cost with respect to $C$. Although there is an
arbitrary number of subsets of $\mathbb{R}^{d}$ of size $k$, it is well-known
that to achieve an approximately optimal clustering, it suffices to only show
correctness on a net of size
$\left(\frac{n}{\varepsilon}\right)^{\mathcal{O}\left(kd\right)}$, and to
adjust the probability of failure and apply a union bound.

Since each point is sampled with probability proportional to its online
sensitivity, the total number of points sampled is proportional to the sum of
the online sensitivities of the points. We thus upper bound the total online
sensitivity as follows. We first break the stream into $\mathcal{O}\left(\log
n\right)$ blocks where the cost of the optimal $(k,z)$-clustering doubles. For
each block, we show that the online sensitivity of the points that arrive in a
block that ends at time $t$ can be charged to either the cost of an optimal
clustering $K_{t}$ for the stream up to $t$ or to the number of points in some
cluster of $K_{t}$.

More specifically, because the online sensitivity of a point $x_{t}$ is
defined as

$\max_{C:\,C\subset\mathbb{R}^{d},\,|C|=k}\frac{{\text{dist}}^{z}(x_{t},C)}{\sum_{i=1}^{t}{\text{dist}}^{z}(x_{i},C)},$

then we can use the generalized triangle inequality to upper bound the online
sensitivity of $x_{t}$ by

$\frac{{\text{dist}}^{z}(x_{t},C)}{\sum_{i=1}^{t}{\text{dist}}^{z}(x_{i},C)}\leq\frac{2^{z-1}{\text{dist}}^{z}(x_{t},\pi(x_{t}))}{\sum_{i=1}^{t}{\text{dist}}^{z}(x_{i},C)}+\frac{2^{z-1}{\text{dist}}^{z}(\pi(x_{t}),C)}{\sum_{i=1}^{t}{\text{dist}}^{z}(x_{i},C)},$

where $\pi$ is the mapping to the closest point in $K_{t}$. We can further
lower bound the denominator $\sum_{i=1}^{t}{\text{dist}}^{z}(x_{i},C)$ by
$2\sum_{i=1}^{t}{\text{dist}}^{z}(x_{i},\pi(x_{i}))$ using the optimality of
$K_{t}$, therefore upper bounding the sum of the first term across all $t$. We
can also use a charging argument to upper bound the second term by
$\mathcal{O}\left(\frac{1}{|S_{t}|}\right)$, where $S_{t}$ is the subset of
$\\{x_{1},\ldots,x_{t}\\}$ mapped to $\pi(x_{t})$ at time $t$. Since there are
$k$ possible clusters for $\pi(x_{t})$,
$\sum_{i=1}^{n}\frac{1}{i}=\mathcal{O}\left(\log n\right)$, and the stream is
partitioned into $\mathcal{O}\left(\log n\right)$ blocks, this informally
gives our $\mathcal{O}\left(k\log^{2}n\right)$ upper bound for the total
online sensitivity, i.e., Theorem 1.3.

###  1.3 Technical Overview: Algorithmic Intuition for Insertion-Only Streams

Despite these efforts, online sensitivity sampling is insufficient for our
ultimate goal because it has an
$\mathcal{O}\left(dk^{2}\log^{2}(n\Delta)\right)$ dependency, which is
prohibitive in both the $k$ and $\log(nd\Delta)$ factors. However, a crucial
observation is that the points sampled by online sensitivity sampling form a
stream ${\mathcal{S}}^{\prime}$ of length
$\operatorname*{poly}\left(k,d,\log(nd\Delta),\frac{1}{\varepsilon}\right)$ of
weighted points that forms a coreset of the input points. Thus if we could
somehow access ${\mathcal{S}}^{\prime}$ in a single pass without storing all
the points of ${\mathcal{S}}^{\prime}$, then the length of
${\mathcal{S}}^{\prime}$ is now small enough for us to simply run a merge-and-
reduce algorithm on ${\mathcal{S}}^{\prime}$.

Unfortunately, online sensitivity sampling requires the storage of the entire
set of sampled points to compute the online sensitivities of future points. A
natural idea would be to use a data structure that could give “good”
approximations to the online sensitivities without using prohibitively large
space. However, if we require correctness at all times, it seems likely that
any analysis that essentially requires a union bound over all times $t\in[n]$
would result in prohibitive $\log(n\Delta)$ factors. Thus, we instead settle
for a data structure that could give “good” approximations to the online
sensitivities “most” of the time, without using prohibitively large space.

Our main insight is that any coreset to the underlying dataset at time $t-1$
precisely gives this guarantee at time $t$. That is, suppose we have a
$(1+\varepsilon)$-coreset at time $t-1$. We can use the coreset to compute
approximations to the online sensitivity of $x_{t}$, which we can then use to
sample $x_{t}$ into a conceptual stream ${\mathcal{S}}^{\prime}$. Because
${\mathcal{S}}^{\prime}$ is formed by online sensitivity sampling, then at all
times $t\in[n]$, the optimal clustering to the set of weighted points in
${\mathcal{S}}^{\prime}$ at time $t$ is a $(1+\varepsilon)$-approximation to
the optimal clustering for the first $t$ points of the original stream. We can
then apply a standard merge-and-reduce algorithm on ${\mathcal{S}}^{\prime}$
to obtain a $(1+\varepsilon)$-coreset to ${\mathcal{S}}^{\prime}$, which
translates to a $(1+\mathcal{O}\left(\varepsilon\right))$-coreset for the
original stream and guarantees that we now have a coreset for time $t$. We can
then iterate on this process to perform online sensitivity sampling at time
$t+1$.

Our analysis for online sensitivity sampling shows that with high probability,
the number of sampled points into ${\mathcal{S}}^{\prime}$ will be
$\operatorname*{poly}\left(k,d,\log(nd\Delta),\frac{1}{\varepsilon}\right)$.
Crucially, since merge-and-reduce requires space that is polylogarithmic in
the length of the input stream, then running merge-and-reduce on
${\mathcal{S}}^{\prime}$ will use $o(\log n)$ words of space. We summarize our
approach in Figure 3.

Bounds on sum ofonline sensitivitiesOnline sensitivitysamplingSmall
implicitstream ${\mathcal{S}}$’Merge-and-reduceframeworkGood estimates
ofonline sensitivities Fig. 3: Flowchart of our simple algorithm and analysis
for $(k,z)$-clustering on insertion-only streams. .

###  1.4 Technical Overview: Algorithmic Intuition for Dynamic Streams

Unfortunately, merge-and-reduce approaches generally do not seem to work for
insertion-deletion, i.e., dynamic streams. Similarly, online sensitivity
sampling and its variants do not seem to immediately work for dynamic streams.
Consider the following definition for the sensitivity of a point
$x\in[\Delta]^{d}$ with respect to a dataset $X=\\{x_{1},\ldots,x_{n}\\}$,
defined as

$\max_{C:\,C\subset\mathbb{R}^{d},\,|C|=k}\frac{{\text{dist}}^{z}(x,C)}{\sum_{i=1}^{n}{\text{dist}}^{z}(x_{i},C)}.$

The issue for dynamic streams is twofold, both related to the possible
subsequent deletion of points: 1) sampled points can later be removed from the
stream, 2) the online sensitivity of a point $x$ at some time in the stream
may be less than the sensitivity of the point $x$ with respect to $X$, because
the data stream may later remove additional points. That is, the online
sensitivity of a point at the time it arrives in the data stream can be either
significantly higher or significantly lower than its sensitivity at the end of
the stream. For example, let $X$ be the dataset at the end of the stream and
let $Y$ be an additional set of points so that $X\cup Y$ is the set of points
when $x$ arrives in the stream. Then the sensitivity of $x$ is defined with
respect to $X$, while the online sensitivity of $x$ is by definition, the
sensitivity of $x$ with respect to $X\cup Y$, and thus less than the true
sensitivity of $x$ with respect to $X$. Although we do not see how to overcome
these issues in a single pass over the data stream and indeed we give
impossibility results for $o(\log n)$ space (in words of memory) algorithms
for one-pass dynamic streams, fortunately it seems these issues can be
resolved using a second pass.

We first show that an approximately optimal clustering can be used to compute
approximations of the sensitivities of the points. Therefore, our hope is to
use a data structure that could give “good” approximations to the optimal
clustering, without using prohibitively large space. To that end, our
algorithm uses two main steps: 1) we first embed the Wasserstein-$z$ distance
into $L_{1}$ and then 2) we use a streaming algorithm to estimate the $L_{1}$
distance. For the ease of discussion, we first consider $z=1$, i.e.,
$k$-median clustering, so that the Wasserstein-$z$ distance corresponds to the
Earth Mover Distance (EMD).

#### EMD sketch for approximate $k$-median clustering cost.

To embed $\operatorname*{EMD}$ on $[\Delta]^{d}$ into $L_{1}$, we first
generalize a standard quadtree consisting of refinements of a randomly shifted
grid, i.e., [IT03, BIRW16]. Suppose without loss of generality that
$\Delta=2^{\ell}$ for some non-negative integer $\ell$. Given a shift
parameter $s=(s_{1},\ldots,s_{d})\in\mathbb{Z}^{d}$ and
$t\in\\{0,1,\ldots,\ell+1\\}$, we define the linear map $G_{s,t}$ from
$\mathbb{R}^{[\Delta]^{d}}$ into $L_{1}$ by first defining the grid
${\mathcal{G}}_{s,t}$ over $\mathbb{Z}^{d}$ with side length $2^{t}$, so that
$s=(s_{1},\ldots,s_{d})$ lies on one of the corners of the grid. For
$\mu\in\mathbb{R}^{[\Delta]^{d}}$, we define $G_{s,t}\mu$ as the frequency
vector whose coordinates correspond to the total mass in each cell/hypercube
of the grid ${\mathcal{G}}_{s,t}$. Then we define the mapping $G_{s}\mu$ to be
the concatenation of the vectors

$G_{s}\mu=(G_{s,0}\mu)\circ(2\cdot G_{s,1}\mu)\circ\ldots\circ(2^{t}\cdot
G_{s,t}\mu)\circ\ldots\circ(2^{\ell}\cdot G_{s,\ell}\mu).$

The main intuition is that for a mass vector $\mu$ and an assignment vector
$\nu$ on $[\Delta]^{d}$, the cost to the optimal transport induced by some
mass in $\mu$ and $\nu$ corresponds to the finest grid the mass does not
appear in the same cell. Due to the scaling of each grid, the $L_{1}$ mass
attributed to the frequency vector $G_{s}(\mu-\nu)$ is proportional to the
cost needed to move the mass from $\mu$ to $\nu$. Thus the grid embedding not
only gives a good approximation to EMD, but also enjoys the property that it
can only overestimate the EMD, which we shall utilize in the analysis.

We also remark that compared to other grid embeddings, our data structure
achieves worse approximation guarantees but crucially uses less space. In
fact, the worse approximation guarantees ultimately only means that we sample
more points into the implicit stream ${\mathcal{S}}^{\prime}$ in such a way
that the downstream guarantees of running a merge-and-reduce algorithm on
${\mathcal{S}}^{\prime}$ will not be affected. Given the embedding of
$\operatorname*{EMD}$ into $L_{1}$, we then use a streaming algorithm to
estimate the $L_{1}$ distance. The guarantee of the $L_{1}$ streaming
algorithm ensures that after the first pass, the accuracy of the estimated
costs of the optimal clusterings ensures that we have a good approximation to
the sensitivity of each point.

#### Putting things together for $k$-median.

Given two passes over the stream, we can use the first pass over the data
stream to maintain the sketch of the EMD embedding. Importantly, the EMD
embedding and the corresponding $L_{1}$ sketch are both linear sketches, and
so they can handle both insertions and deletions.

We then use the second pass over the data stream to perform sensitivity
sampling, since we can now compute the sensitivity of each point with respect
to the final data set $X$ at the end of the stream. Unfortunately, sensitivity
sampling is still incompatible with the dynamic setting, since sampled points
could be subsequently deleted. Thus, we instead subsample from the universe
with probability proportional to each point’s sensitivity, so that if a
universe element is subsampled, all updates to the element throughout the
stream are reported. We then use a sparse recovery scheme to track the updates
to the sampled universe elements. Since there are known linear sketches for
sparse recovery, then our algorithm can handle both insertions and deletions
while simulating sensitivity sampling. Thus at the end of the stream, our
algorithm outputs a weighted subset of the input points that forms a
$(1+\varepsilon)$-coreset of the underlying point set for $k$-median
clustering.

###  1.5 Technical Overview: Embedding for $(k,z)$-Clustering

To extend our approach to $(k,z)$-clustering, it seems crucial to develop an
analog for EMD sketching for the Wasserstein-$z$ distance. Unfortunately
despite substantial efforts, such an embedding is not known for a single
quadtree, because the distortion incurred by a splitting hyperplane in the
quadtree is exponential in $z$, while the probability that two points are
split by a hyperplane is inversely linear in their distance. In fact, there
are simple examples on a line that show that the expected distortion between
the squared distances of a set of $n$ points, i.e., $z=2$, and the estimated
distance by a quadtree is $\Omega(n)$. See Appendix A for one such example.
Hence while $k$-median, i.e., $z=1$ is ideal for the quadtree approach, larger
values of $z$ can have unbounded distortion.

Recently, [CLN+20] overcame this barrier in the offline setting by considering
multiple quadtrees and taking the minimum estimated distance for pairs of
points across the quadtrees. This approach does not seem to work for the
streaming setting, because we can only separately embed each of the quadtrees
into $L_{1}$ and then we will only have access to the estimated clustering
costs by each of the quadtrees. However, the minimum estimate of the sum of
the squared distances is not equal to the sum of the minimum estimated squared
distances; the former is what the embedding into $L_{1}$ would give, but the
latter is what the approach of [CLN+20] requires.

#### Wasserstein-$z$ embedding preliminaries.

We instead develop a quadtree embedding technique that achieves a bicriteria
approximation for $(k,z)$-clustering. As in the EMD sketch, we select a shift
$s=(s_{1},\ldots,s_{d})\in[\Delta]^{d}$ uniformly at random and then add $s$
to each of the input points $x_{1},\ldots,x_{n}\in[\Delta]^{d}$. We create a
quadtree so that the root of the tree represents a grid with side length
$2\Delta$. We then partition the grid into $2^{d}$ smaller hypercubes with
side length $\Delta$, so that for each of these hypercubes that contains a
shifted input point, we create a node representing the hypercube and add the
node as a child of the root note in the tree, using an edge with weight
$\sqrt{d}\Delta$. We repeat this procedure until every node representing a
hypercube contains at most a single point, resulting in a tree where all
leaves contain a single point and are all at the same height, which is at most
$\mathcal{O}\left(\log\Delta\right)$.

We embed the input points $X$ into the quadtree. As it will be useful to index
from the leaves of the tree, we define the $t$-th level of the quadtree to
correspond with the hypergrid with cells of length $2^{t}$. For
$\mu\in\mathbb{R}^{[\Delta]^{d}}$, we define $W_{s,t}\mu$ to be the frequency
vector over the hypercubes of the hypergrid ${\mathcal{G}}_{s,t}$ at level $t$
that counts the total mass in each hypercube. We define the mapping $W_{s}\mu$
to be the concatenation of the vectors

$W_{s}\mu=(W_{s,0}\mu)\circ((2\sqrt{d})^{z}\cdot
W_{s,1}\mu)\circ\ldots\circ((2^{z}\sqrt{d})^{z}\cdot
W_{s,t}\mu)\circ\ldots\circ((2^{\ell}\sqrt{d})^{z}\cdot W_{s,\ell}\mu).$

We remark that up to this point, the approach is the same as previous quadtree
embeddings.

#### Wasserstein-$z$ embedding through bicriteria approximation.

Now, for a query point $q$, we say that $q$ is _bad at level $i$_ if there
exists a hyperplane of the quadtree decomposition of length $2^{i}$ that has
distance less than $\frac{2^{i}}{d\log\Delta}$ from $q$. Otherwise, we say
that $q$ is good at level $i$. Observe that if $q$ is good at level $i$, but
$q$ and $x\in X$ are separated at level $i$ but not $i+1$, then
$\|q-x\|_{2}^{z}>\frac{2^{iz}}{d^{z}\log^{z}\Delta}$, since all hyperplanes at
level $i$ are at least distance $\frac{2^{i}}{d^{2}\log^{2}\Delta}$ from $q$.

Moreover since $q$ and $x$ are not separated at level $i+1$, then the incurred
estimated cost for $q$ and $x$ in the quadtree is most $(2^{i}\sqrt{d})^{z}$,
so the distortion will be at most $d^{0.5z}\log^{z}\Delta$.

On the other hand, if $q$ is bad at level $i$, then by definition, it has
distance less than $\frac{2^{i}}{d\log\Delta}$ from a hyperplane of the
quadtree decomposition of length $2^{i}$ and so $q$ and $x$ incur estimated
cost $2^{iz}d$ by the quadtree, then the distortion could be significantly
larger. To address this issue, we define a mapping $\phi$ to subsets $S_{q}$
of $\mathbb{R}^{d}$ for a query point $q$ as follows. We first add $q$ to
$S_{q}$. If $q$ is bad at level $i$, then consider each hyperplane $H$ of the
quadtree of length $2^{i}$ that is too close to $q$, i.e.,
${\text{dist}}(q,H)\leq\frac{2^{i}}{d\log\Delta}$. We create a copy
$q^{(H)}_{i}$ corresponding to the projection of $q$ onto $H$, so that
$q^{(H)}_{i}$ is responsible for serving the points assigned to $q$ that are
on the other side of $H$. We add $q^{(H)}_{i}$ to $S_{q}$ and proceed top-
down, repeatedly adding points to $S_{q}$ as necessary. See Figure 4 for an
example of this process.

Note that if $q$ is bad at level $i$, it could be too close to multiple
hyperplanes of the quadtree with length $2^{i}$; in this case, we add a point
$q^{(H)}_{i}$ to $S_{q}$ for each hyperplane $H$ for which
${\text{dist}}(q,H)\leq\frac{2^{i}}{d^{2}\log^{2}\Delta}$. That is, $q$ could
induce multiple points $q^{(H)}_{i}$ to be added to $S_{q}$.

$\triangle$$\times$$\triangle$$\triangle$$\triangle$$\triangle$ Fig. 4:
Example of quadtree for Wasserstein-$z$ embedding. Black circles are the input
points and black triangles denote the query centers. One of the four queries
(top left cell) is too close to splitting hyperplanes, which could cause too
much distortion, since the query is mapped to the cell centers denoted by the
black X. Hence, we create an additional center on the cell boundary, marked by
the red triangle.

By construction, the resulting mapping $\phi(C):=\cup_{j=1}^{k}\phi(C_{j})$
for a set of $k$ centers $C=\\{C_{1},\ldots,C_{k}\\}$ will induce an estimated
cost by the quadtree that has small distortion, i.e.,
$d^{0.5z}\log^{z}\Delta$, from the actual cost. With more careful fine-grained
analysis, we show that the expected distortion is at most
$\mathcal{O}\left(d^{1+0.5z}\log^{z-1}\Delta\right)$. To ensure the bicriteria
guarantee, it remains to show that $\phi(C)$ contains at most
$\mathcal{O}\left(k\right)$ centers in expectation.

Consider a level $i$ and a center $C_{j}$ with $j\in[k]$. Let ${\mathcal{H}}$
be the set of hyperplanes of the hypergrid at level $i$. By a union bound over
the $d$ dimensions, the probability that
${\text{dist}}(C_{j},{\mathcal{H}})<\frac{2^{i}}{d\log\Delta}$ is at most
$\mathcal{O}\left(\frac{d}{d\log\Delta}\right)=\mathcal{O}\left(\frac{1}{\log\Delta}\right)$.
Hence, summing up over all the $\mathcal{O}\left(\log\Delta\right)$ levels,
the expected number of center copies that a center $C_{j}$ can generate into
$\phi(C)$ is at most $\mathcal{O}\left(1\right)$. Finally, summing up over
$j\in[k]$, the expected number of generated points is at most
$\mathcal{O}\left(k\right)$ and thus the expected size of $\phi(C)$ is at most
$k+\mathcal{O}\left(k\right)=\mathcal{O}\left(k\right)$.

#### Putting things together for $(k,z)$-clustering.

To adapt our two-pass dynamic streaming $k$-median clustering algorithm for
$k$-means clustering, it is instructive to consider which steps break down.
Recall that the first pass of our $k$-median algorithm was used to set up the
EMD embedding, which subsequently gave a $\mathcal{O}\left(\log
k+\log\log\Delta\right)$ approximation to the sensitivity of each point in the
second pass, so that we could perform sensitivity sampling through sparse
recovery in the second pass.

For $k$-means clustering, the natural approach would be to replace the EMD
embedding with our Wasserstein-$z$ embedding. However, our Wasserstein-$z$
embedding may distort the estimated sensitivity by a factor of
$\mathcal{O}\left(d^{1+0.5z}\log^{z-1}\Delta\right)$, in which case
sensitivity sampling in the second pass could require sampling
$\Omega(\log^{z-1}\Delta)$ points. Unfortunately, even for $z=2$, this no
longer gives $\mathcal{O}\left(\log(n\Delta)\right)$ total bits of space if we
store each sampled point using $\mathcal{O}\left(\log(n\Delta)\right)$ bits of
space.

To overcome this issue, we store an approximation of its _offset_ from each of
the centers instead of storing an explicit representation of each point. That
is, for a point $x$ and a set $C^{\prime}$ of $\mathcal{O}\left(k\right)$
centers, suppose $c^{\prime}(x)$ is the closest center of $C^{\prime}$ to $x$.
Let the offset from $x$ to $c^{\prime}(x)$ be defined by
$y^{\prime}=x-c^{\prime}(x)$. Suppose we round each coordinate of $y^{\prime}$
to a power of
$\left(1+\operatorname*{poly}\left(\varepsilon,\frac{1}{d},\frac{1}{\log\Delta}\right)\right)$
to form a vector $y$.

Let $Y$ be the set of all points of $X$ rounded in this manner. We show that
the cost of any clustering $C$ on $Y$ is a $(1+\varepsilon)$-approximation to
the cost of the clustering $C$ and $X$, so it suffices to instead consider
$(k,z)$-clustering on $Y$. The main insight from considering $Y$ is that due
to the points all being rounded to a power of
$\left(1+\operatorname*{poly}\left(\varepsilon,\frac{1}{d},\frac{1}{\log\Delta}\right)\right)$
away from a point in $C^{\prime}$, then there is a compact representation of
each point in $Y$ that does not require $\Omega(d\log\Delta)$ bits to
represent. In fact, to represent each point $y$, it suffices to store the
identity of $c^{\prime}(x)$ as well as the _exponents_ of the offsets, which
only requires $\mathcal{O}\left(d\log\frac{\log\Delta}{\varepsilon}\right)$
bits per sample. Hence for $k$-means clustering, the algorithm still uses
$o(\log^{2}(n\Delta))$ total bits of space.

###  1.6 Preliminaries

For an integer $n>0$, we use $[n]$ to denote the set $\\{1,\ldots,n\\}$. We
use $\operatorname*{poly}(n)$ to denote a fixed polynomial in $n$ and
$\operatorname*{polylog}(n)$ to denote $\operatorname*{poly}(\log n)$. If an
event occurs with probability at least $1-\frac{1}{\operatorname*{poly}(n)}$,
we say the event occurs with high probability.

For vectors $x,y\in\mathbb{R}^{d}$, we use ${\text{dist}}(x,y)$ to denote the
Euclidean distance $\|x-y\|_{2}$, so that
$\|x-y\|_{2}^{2}=\sum_{i=1}^{d}(x_{i}-y_{i})^{2}$. More generally, we define
the $L_{z}$ norm of $x$ by $\|x\|_{z}$ so that
$\|x\|_{z}^{z}=\sum_{i=1}^{d}x_{i}^{z}$. For a set $S$, we use
${\text{dist}}(x,S)$ to denote $\min_{y\in S}{\text{dist}}(x,y)$ and
similarly, for sets $X$ and $S$, we define ${\text{dist}}(X,S)=\min_{x\in
X,y\in S}{\text{dist}}(x,y)$. For a fixed $z\geq 1$ and sets
$X,C\subset\mathbb{R}^{d}$ with $X=\\{x_{1},\ldots,x_{n}\\}$ we use
$\operatorname*{cost}(X,C)$ to denote
$\sum_{i=1}^{n}{\text{dist}}(x_{i},C)^{z}$.

######  Fact 1.9 (Generalized triangle inequality).

For any $z\geq 1$ and $x,y,z\in\mathbb{R}^{d}$, we have

${\text{dist}}(x,y)^{z}\leq
2^{z-1}({\text{dist}}(x,w)^{z}+{\text{dist}}(w,y)^{z}).$

Let $x,y\in\mathbb{R}^{n}$ with $n=[\Delta]^{d}$ be two vectors such that
$\|x\|_{1}=\|y\|_{1}$. Let $F(x,y)$ be the family of functions
$f:[\Delta]^{d}\times[\Delta]^{d}\to\mathbb{R}^{\geq 0}$ so that for any
$i\in[\Delta]^{d}$, we have $\sum_{j\in[\Delta]^{d}}f(i,j)=x_{i}$ and for any
$j\in[\Delta]^{d}$, we have $\sum_{i\in[\Delta]_{d}}f(i,j)=y_{j}$. Then we
define the earth mover distance (EMD) between $x$ and $y$ by

$\operatorname*{EMD}(x,y)=\min_{f\in
F}\sum_{i,j\in[\Delta]^{d}}f(i,j)\|i-j\|_{2}.$

For general vectors $x,y\in\mathbb{R}^{n}$, we define

$\operatorname*{EMD}(x,y)=\inf_{x^{\prime}\preceq x,y^{\prime}\preceq
y,\|x^{\prime}\|_{1}=\|y^{\prime}\|_{=}1}\operatorname*{EMD}(x^{\prime},y^{\prime})+d\Delta(\|x-x^{\prime}\|_{1}+\|y-y^{\prime}\|_{1}).$

More generally, we define the quantity $\operatorname*{WASSZZ}(x,y)$ for
$x,y\in\mathbb{R}^{n}$ with $n=[\Delta]^{d}$ and $\|x\|_{1}=\|y\|_{1}$ by

$\operatorname*{WASSZZ}(x,y)=\min_{f\in
F}\sum_{i,j\in[\Delta]^{d}}f(i,j)\|i-j\|_{2}^{z}$

and for general vectors $x,y\in\mathbb{R}^{n}$ by

$\operatorname*{WASSZZ}(x,y)=\inf_{x^{\prime}\preceq x,y^{\prime}\preceq
y,\|x^{\prime}\|_{1}=\|y^{\prime}\|_{=}1}\operatorname*{WASSZZ}(x^{\prime},y^{\prime})+d^{z}\Delta^{z}(\|x-x^{\prime}\|_{1}+\|y-y^{\prime}\|_{1}).$

We remark that the quantity $\operatorname*{WASSZZ}(x,y)$ corresponds to the
$z$-th power of the Wasserstein-$z$ distance between $x$ and $y$; we use
$\operatorname*{WASSD}(x,y)$ to denote the Wasserstein-$z$ distance between
$x$ and $y$, in general omitting the dependency on $z$ when it is clear from
context.

We also define $x^{+}=\frac{|x|+x}{2}$ to be the vector containing the
positive entries of $x$ and $x^{-}=x-x^{+}$ to be the vector containing the
negative entries of $x$. Then we define

$\|x\|_{\operatorname*{EMD}}=\operatorname*{EMD}(x^{+},x^{-}).$

We now formally define the notion of a (strong) coreset.

######  Definition 1.10 (Coreset).

Given an approximation parameter $\varepsilon>0$, and a set $X$ of points
$x_{1},\ldots,x_{n}\in\mathbb{R}^{d}$ with distance function dist, a _coreset_
for $(k,z)$ clustering is a subset $S$ of weighted points of $X$ with weight
function $w$ such that for any set $C$ of $k$ points, we have

$(1-\varepsilon)\sum_{t=1}^{n}{\text{dist}}(x_{t},C)^{z}\leq\sum_{q\in
S}w(q){\text{dist}}(q,S)^{z}\leq(1+\varepsilon)\sum_{t=1}^{n}{\text{dist}}(x_{t},C)^{z}.$

Currently, the state-of-the-art coreset construction for $(k,z)$-clustering is
the following:

######  Theorem 1.11.

[CLSS22] Given an accuracy parameter $\varepsilon\in(0,1)$, there exists a
coreset construction for $(k,z)$-clustering that samples
$\tilde{\mathcal{O}}\left(\frac{k}{\varepsilon^{2}}2^{z\log
z}\cdot\min\left(\frac{1}{\varepsilon^{z}},k\right)\right)$ weighted points.

Thus, merge-and-reduce on a stream of length $n$ using the coreset
construction of [CLSS22] in Theorem 1.11 offers the following guarantee:

######  Lemma 1.12.

There exists a one-pass streaming algorithm on insertion-only data streams of
length $n$ that, with probability at least $0.99$, simultaneously outputs a
$(1+\varepsilon)$-approximation to $(k,z)$-clustering at all times of the
stream, using $\tilde{\mathcal{O}}\left(\frac{k}{\varepsilon^{2}}2^{z\log
z}\cdot\left(\frac{1}{\varepsilon^{z}},k\right)\right)\cdot\operatorname*{polylog}(n)$
bits of space.

We recall the following property of the Johnson-Lindenstrauss (JL)
transformation.

######  Theorem 1.13 (Johnson-Lindenstrauss lemma).

[JL84] Let $X\subset\mathbb{R}^{d}$ be a set of $n$ points and
$m=\mathcal{O}\left(\frac{1}{\varepsilon^{2}}\log n\right)$. There exists a
family of random linear maps ${\mathfrak{M}}:\mathbb{R}^{d}\to\mathbb{R}^{m}$
such that with probability at least $0.99$ over the choice of
${\mathfrak{M}}$, we simultaneously have for all $x,y\in X$,

$(1-\varepsilon)\|x-y\|_{2}\leq\|{\mathfrak{M}}x-{\mathfrak{M}}y\|_{2}\leq(1+\varepsilon)\|x-y\|_{2}.$

We also recall the following concentration inequalities:

######  Theorem 1.14 (Bernstein’s concentration inequality).

[Ber27] Let $X_{1},\ldots,X_{m}$ be independent random variables with
$\mathbb{E}\left[X_{i}\right]<\infty$ and $X_{i}\geq 0$ for all $i\in[m]$. Let
$X=\sum_{i=1}^{m}X_{i}$ and let $\gamma>0$. Then

$\mathbf{Pr}\left[X\leq\mathbb{E}\left[X\right]-\gamma\right]\leq\exp\left(\frac{-\gamma^{2}}{2\sum_{i}\mathbb{E}\left[X_{i}^{2}\right]}\right).$

Moreover, if $X_{i}-\mathbb{E}\left[X_{i}\right]\leq\Delta$ for all $i\in[m]$,
then for
$\sigma_{i}^{2}:=\mathbb{E}\left[X_{i}^{2}\right]-\mathbb{E}\left[X_{i}\right]^{2}$,

$\mathbf{Pr}\left[X\geq\mathbb{E}\left[X\right]+\gamma\right]\leq\exp\left(\frac{-\gamma^{2}}{2\sum_{i}\sigma_{i}^{2}+2\gamma\Delta/3}\right).$

######  Theorem 1.15 (Freedman’s inequality).

[Fre75] Suppose $Y_{0},Y_{1},\ldots,Y_{n}$ is a scalar martingale with
difference sequence $X_{1},\ldots,X_{n}$. Specifically, we initiate $Y_{0}=0$
and set $Y_{i}=Y_{i-1}+X_{i}$ for all $i\in[n[$ Let $R\geq|X_{t}|$ for all
$t\in[n]$ with high probability. We define the predictable quadratic variation
process of the martingale by
$w_{k}:=\sum_{t=1}^{k}\underset{t-1}{\mathbb{E}}\left[X_{t}^{2}\right]$, for
$k\in[n]$. Then for all $\varepsilon\geq 0$ and $\sigma^{2}>0$, and every
$k\in[n]$,

$\mathbf{Pr}\left[\max_{t\in[k]}|Y_{t}|>\varepsilon\text{ and
}w_{k}\leq\sigma^{2}\right]\leq
2\exp\left(-\frac{\varepsilon^{2}/2}{\sigma^{2}+R\varepsilon/3}\right).$

##  2 EMD Sketch

We first describe the embedding of $\operatorname*{EMD}$ on $[\Delta]^{d}$
into $L_{1}$. Suppose without loss of generality that $\Delta=2^{\ell}$ for
some non-negative integer $\ell$.

Given $s=(s_{1},\ldots,s_{d})\in\mathbb{Z}^{d}$ and
$t\in\\{0,1,\ldots,\ell\\}$, we define the linear map $G_{s,t}$ from
$\mathbb{R}^{[\Delta]^{d}}$ into $L_{1}$ by first defining the grid
${\mathcal{G}}_{s,t}$ over $\mathbb{Z}^{d}$ with side length $2^{t}$, so that
$s=(s_{1},\ldots,s_{d})$ lies on one of the corners of the grid. Then for
$\mu\in\mathbb{R}^{[\Delta]^{d}}$, we define $G_{s,t}\mu$ as the frequency
vector over the hypercubes of the grid ${\mathcal{G}}_{s,t}$ that counts the
total mass in each hypercube. To avoid ambiguity, we say that each cell of a
grid of length $2^{t}$ has closed boundaries on one side and open boundaries
on the other side, e.g., a cell that contains a point $(x_{1},\ldots,x_{d})$
may contain points $(y_{1},\ldots,y_{d})$ where $y_{i}<x_{i}+2^{k}$, but does
not contain any points $(y_{1},\ldots,y_{d})$ where $y_{i}\geq x_{i}+2^{k}$.
Then we define the mapping $G_{s}\mu$ to be the concatenation of the vectors

$G_{s}\mu=(G_{s,0}\mu)\circ(2\cdot G_{s,1}\mu)\circ\ldots\circ(2^{t}\cdot
G_{s,t}\mu)\circ\ldots\circ(2^{\ell}\cdot G_{s,\ell}\mu).$

We first require the following structural property: The following is analogous
to Claim 6.1 in [BIRW16], generalized to dimension $d$.

######  Lemma 2.1.

Let $X_{1},\ldots,X_{n}$ be non-negative random variables not necessarily
independent, such that for every $i\in[n]$ and $d,t>0$,

$\mathbf{Pr}\left[X_{i}\geq t\right]\leq\frac{\gamma d}{t},$

for an absolute constant $\gamma>0$. Let $S=\sum_{i\in[n]}\alpha_{i}X_{i}$ for
a set of coefficients $\alpha_{1},\ldots,\alpha_{n}\geq 0$ with
$\sum_{i\in[n]}\alpha_{i}=1$. Then there exists an absolute constant $B>0$
such that for every $\delta\in(0,1)$,

$\mathbf{Pr}\left[S\geq\frac{B\gamma
d}{\delta}\left(H(\alpha)+\log\frac{2\gamma
d}{\delta}\right)\right]\leq\delta,$

where $H(\alpha)$ is the entropy of the distribution over $[n]$ induced by
$\alpha_{1},\ldots,\alpha_{n}$.

###### Proof.

Let $T_{1},\ldots,T_{n}\geq 0$ be non-negative parameters to be fixed. We
define ${\mathcal{E}}$ to be the event that $X_{i}\leq T_{i}$ simultaneously
for all $i\in[n]$. By the assumption that $\mathbf{Pr}\left[X_{i}\geq
T_{i}\right]\leq\frac{\gamma d}{T_{i}}$ for all $i\in[n]$ and a union bound,
we have that

$\mathbf{Pr}\left[\neg{\mathcal{E}}\right]\leq\sum_{i\in[n]}\frac{\gamma
d}{T_{i}}.$

Moreover, there exists an absolute constant $B$ such that
$\mathbb{E}\left[X_{i}\mid{\mathcal{E}}\right]\leq B\gamma d\log T_{i}$. Hence
by Markov’s inequality,

$\mathbf{Pr}\left[S\leq\frac{B\gamma
d}{\delta}\left(\sum_{i\in[n]}\alpha_{i}\log
T_{i}\right)\mid{\mathcal{E}}\right]\geq 1-\frac{\delta}{2}.$

By setting $T_{i}=\frac{2\gamma d}{\alpha_{i}\delta}$, we have that

$\mathbf{Pr}\left[\neg{\mathcal{E}}\right]\leq\sum_{i\in[n]}\frac{\gamma
d}{T_{i}}\leq\sum_{i\in[n]}\frac{\alpha_{i}\delta}{2}=\frac{\delta}{2}$

and

$\sum_{i\in[n]}\alpha_{i}\log T_{i}=H(\alpha)+\log\frac{2\gamma d}{\delta}.$

∎

We next generalize a statement lower bounding the contraction of the EMD
embedding into $L_{1}$ by [BIRW16, IT03] from $[\Delta]^{2}$ to
$[\Delta]^{d}$. Namely, [BIRW16, IT03] showed that for $[\Delta]^{2}$, the
contraction is upper bounded by a constant amount. Here we show that for
$[\Delta]^{d}$, the contraction is at most $\mathcal{O}\left(\sqrt{d}\right)$.

######  Lemma 2.2.

For every $\mu,\nu\in\operatorname*{EMD}_{[\Delta]^{d}}$ and every
$s=(s_{1},\ldots,s_{d})\in\mathbb{Z}^{d}$, it holds that
$\|\mu-\nu\|_{\operatorname*{EMD}}\leq\mathcal{O}\left(\sqrt{d}\right)\cdot\|G_{s}(\mu-\nu)\|_{1}$.

###### Proof.

Suppose without loss of generality that $\mu$ and $\nu$ have non-intersecting
support, since otherwise the intersection can be subtracted from both $\mu$
and $\nu$. Moreover, suppose without loss of generality that $\mu$ and $\nu$
are probability measures, i.e., $\|\mu\|_{1}=\|\nu\|_{1}=1$.

We first consider the matching induced by pairing mass contained in the same
cells of grids, starting with ${\mathcal{G}}_{s,0}$ and then iterating through
${\mathcal{G}}_{s,1},\ldots$. Since $\mu$ and $\nu$ have non-intersecting
support and all points in $[\Delta]^{d}$ fall in separate cells in
${\mathcal{G}}_{s,0}$, then we have $\|G_{s,0}(\mu-\nu)\|_{1}=2$. Then there
is at most $1-\frac{1}{2}\|G_{s,1}(\mu-\nu)\|_{1}$ mass that can be matched
together in ${\mathcal{G}}_{s,1}$, which induces cost at most
$2\sqrt{d}(1-\frac{1}{2}\|G_{s,1}(\mu-\nu)\|_{1})$.

Next, at most
$\frac{1}{2}\|G_{s,1}(\mu-\nu)\|_{1}-\frac{1}{2}\|G_{s,2}(\mu-\nu)\|_{1}$ mass
can be matched together in ${\mathcal{G}}_{s,2}$, inducing cost at most

$4\sqrt{d}\left(\frac{1}{2}\|G_{s,1}(\mu-\nu)\|_{1}-\frac{1}{2}\|G_{s,2}(\mu-\nu)\|_{1}\right).$

More generally, the cost induced by grid ${\mathcal{G}}_{s,i}$ can be at most

$2^{i}\sqrt{d}\left(\frac{1}{2}\|G_{s,1}(\mu-\nu)\|_{1}-\frac{1}{2}\|G_{s,2}(\mu-\nu)\|_{1}\right).$

Hence, we have

$\displaystyle\|\mu-\nu\|_{\operatorname*{EMD}}$
$\displaystyle\leq\frac{\sqrt{d}}{2}\sum_{i=0}^{\ell-1}2^{i}(\|G_{s,i}(\mu-\nu)\|_{1}-\|G_{s,i+1}(\mu-\nu)\|_{1})$
$\displaystyle\leq\frac{\sqrt{d}}{2}\|G_{s}(\mu-\nu)\|_{1},$

since

$G_{s}\mu=(G_{s,0})\mu\circ(2\cdot G_{s,1}\mu)\circ\ldots\circ(2^{t}\cdot
G_{s,t}\mu)\circ\ldots\circ(2^{\ell}\cdot G_{s,\ell}\mu).$

∎

Before lower bounding the dilation of the EMD embedding, we first define the
following concept to quantify the granularity at which the randomly shifted
grid intersects a specific $L_{1}$ ball.

######  Definition 2.3.

Let $x\in\mathbb{R}^{d}$, $R\in(0,2\Delta]$, and let
$s=(s_{1},\ldots,s_{d})\in[\Delta]^{d}$ be sampled uniformly at random. Let

${\mathcal{A}}_{x,R}(s):=\frac{\min\left(\\{2^{t}\,\mid\,{\mathcal{G}}_{s,t}\text{
does not cut }B_{L_{1}}(x,R)\\}\right)}{R}.$

Intuitively, ${\mathcal{A}}_{x,R}(s)$ is the side length of the finest grid
among the $\ell+1$ grids ${\mathcal{G}}_{s,0},\ldots,{\mathcal{G}}_{s,\ell}$
that does not intersect the $L_{1}$ ball $B_{L_{1}}(x,R)$.

We next generalize Lemma 6.2 in [BIRW16] to $d$-dimensional space.
Intuitively, the following statement shows that the probability that coarser
grids, i.e., grids containing cells with larger side lengths, are linearly
less likely to intersect with a fixed $L_{1}$ ball.

######  Lemma 2.4.

There exists a universal constant $\gamma>0$ such that for every
$x\in\mathbb{R}^{d}$, $R\in(0,2\Delta]$ and $T>0$,

$\underset{s\in[\Delta]^{d}}{\mathbf{Pr}}\left[{\mathcal{A}}_{x,R}(s)\geq
T\right]\leq\frac{\gamma d}{T}.$

###### Proof.

Observe that by definition, ${\mathcal{A}}_{x,R}(s)\geq T$ if and only if the
grid with the largest side length less than $RT$ intersects with the ball
$B_{L_{1}^{d}}(x,R)$. Note that for each of the $d$ dimensions, the
probability that $B_{L_{1}^{d}}(x,R)$ intersects with the grid of length $RT$
is at most $\frac{\gamma}{T}$ for some absolute constant $\gamma>0$. Thus by a
union bound, we have that
$\underset{s\in[\Delta]^{d}}{\mathbf{Pr}}\left[{\mathcal{A}}_{x,R}(s)\geq
T\right]\leq\frac{\gamma d}{T}$. ∎

Similarly, we generalize Lemma 6.3 in [BIRW16] to $[\Delta]^{d}$. Intuitively,
the following statement relates the embedding of the difference of two
elementary vectors to the quantity ${\mathcal{A}}_{u,R}(s)$ representing the
side length of the finest grid that does not cut the $L_{1}$ ball
$B_{L_{1}}(u,R)$.

######  Lemma 2.5.

Let $x,y\in[\Delta]^{d}$ and let $s=(s_{1},\ldots,s_{d})\in[\Delta]^{d}$ be
sampled uniformly at random. Let $e_{x}$ and $e_{y}$ be the basis vectors
corresponding to $x$ and $y$ respectively. Then there exists an absolute
constant $\gamma>0$ such that $\|G_{s}(e_{x}-e_{y})\|_{1}\geq\gamma\cdot
R\cdot{\mathcal{A}}_{u,R}(s)$ for every $u\in\mathbb{R}^{2}$ and
$R\in(0,2\Delta]$ such that the ball $B_{L_{1}}(u,R)$ contains both $x$ and
$y$.

###### Proof.

Observe that by the definition of ${\mathcal{A}}_{u,R}$, no grid with side
length at least $R\cdot{\mathcal{A}}_{u,R}(s)$ contributes to
$\|G_{s}(e_{x}-e_{y})\|_{1}$. On the other hand, all grids with side length
less than $R\cdot{\mathcal{A}}_{u,R}(s)$ contribute a geometric series toward
$\|G_{s}(e_{x}-e_{y})\|_{1}$ with total sum at most $\gamma\cdot
R\cdot{\mathcal{A}}_{u,R}(s)$ for some constant $\gamma>0$. ∎

The following statement is the generalization of Claim 6.2 in [BIRW16] to
$[\Delta]^{d}$. The main point is to we show that we can upper bound
$\|G_{s}\mu\|_{1}$ in a way that the assumptions of Lemma 2.1 can be applied.
In the following statement, we say a random variable $X$ statistically
dominates a random variable $Y$ if for all real numbers $t$, we have
$\mathbf{Pr}\left[X\geq t\right]\geq\mathbf{Pr}\left[Y\geq t\right]$.

######  Lemma 2.6.

Let $\mu$ and $\nu$ be two probability measures over $[\Delta]^{d}$ and
suppose the optimal transportation of $\mu$ to $\nu$ consists of moving mass
$w_{i}$ from the point $x_{i}\in[\Delta]^{d}$ to the point
$y_{i}\in[\Delta]^{d}$, for all $i\in[k]$. Let
$\\{B_{j}=B_{L_{1}}(u_{j},R_{j})\\}_{j=1}^{q}$ be a collection of $L_{1}$
balls in $d$-dimensional space such that for every $i\in[k]$, there exists
$f(i)\in[q]$ such that both $x_{i}$ and $y_{i}$ are contained within
$B_{f(i)}$. For each $j\in[q]$, let

$\widetilde{w_{j}}=\sum_{i:f(i)=j}w_{i}.$

Let $s=(s_{1},\ldots,s_{d})\in[\Delta]^{d}$ be sampled uniformly at random.
Then there exists an absolute constant $\gamma>9$ such that the random
variable

$\|G_{s}(\mu-\nu)\|_{1}\leq\sum_{i\in[k]}w_{i}\|G_{s}(e_{x_{i}}-e_{y_{i}})\|_{1}$

is statistically dominated by $S=\sum_{j\in[k]}\widetilde{w_{j}}\cdot
R_{j}\cdot X_{i}$ for some non-negative and not necessarily independent random
variables $X_{1},\ldots,X_{k}$ with the property for every $i$ and $t>0$,

$\mathbf{Pr}\left[X_{i}\geq T\right]\leq\frac{\gamma d}{T}.$

###### Proof.

The claim immediately follows by applying Lemma 2.4, Lemma 2.5, and the
triangle inequality. ∎

The following statement is the generalization of Claim 6.3 in [BIRW16] to
$[\Delta]^{d}$. Intuitively, we upper bound the $L_{1}$ norm of the embedding
$G_{s}(\mu-\nu)$ in terms of the entropy $H(\alpha)$ of the coefficient vector
$\alpha$.

######  Lemma 2.7.

Assuming the notation and conditions from Lemma 2.6,

$\mathbf{Pr}\left[\|G_{s}(\mu-\nu)\|_{1}\leq\mathcal{O}\left(d\right)\cdot
H(\alpha)\cdot T\right]\geq 0.999,$

where $T=\sum_{j=1}^{q}\widetilde{w_{j}}R_{j}=\sum_{i=1}^{k}w_{i}R^{f(i)}$ and
$\alpha_{j}=\frac{\widetilde{w_{j}}R_{j}}{T}$ for $j\in[q]$.

###### Proof.

The claim immediately holds by Lemma 2.6 and Lemma 2.1. ∎

We now upper bound the dilation of the EMD embedding into $L_{1}$. Namely, we
show that the embedding into $L_{1}$ can be at most a multiplicative
$\mathcal{O}\left(d(\log k+\log\log\Delta)\right)$ factor larger than the EMD.

######  Lemma 2.8.

Suppose that $\mu$ and $\nu$ are two probability measures over $[\Delta]^{d}$
and in addition, $\mu$ has support of size at most $k$. Then

$\mathbf{Pr}\left[\|G_{s}(\mu-\nu)\|_{1}\leq\mathcal{O}\left(d(\log
k+\log\log\Delta)\right)\cdot\|\mu-\nu\|_{\operatorname*{EMD}}\right]\geq
0.999.$

###### Proof.

Let the support of $\mu$ be $x_{1},\ldots,x_{k}\in[\Delta]^{d}$. Let
${\mathcal{B}}$ be the family of $\mathcal{O}\left(k\log\Delta\right)$ balls
$\\{B(x_{i},2^{j-1})\\}_{i\in[k],j\in[2+\log\Delta]}$.

Consider the optimal transportation from $\mu$ to $\nu$. Each edge of length
$L$ in this transportation is contained in a ball of radius
$\mathcal{O}\left(L\right)$. Hence, by taking
$T\leq\mathcal{O}\left(1\right)\cdot\|\mu-\nu\|_{\operatorname*{EMD}}$ in
Lemma 2.7, we have

$\mathbf{Pr}\left[\|G_{s}(\mu-\nu)\|_{1}\leq\mathcal{O}\left(t\log
t\right)\cdot
H(\alpha)\cdot\|\mu-\nu\|_{\operatorname*{EMD}}\right]\leq\frac{1}{t},$

and it remains to upper bound $H(\alpha)$. Since the support of $\alpha$ has
size at most $\mathcal{O}\left(k\log\Delta\right)$, then
$H(\alpha)\leq\mathcal{O}\left(\log(k\log\Delta)\right)\leq\mathcal{O}\left(\log
k+\log\log\Delta\right)$. Therefore,

$\mathbf{Pr}\left[\|G_{s}(\mu-\nu)\|_{1}\leq\mathcal{O}\left(d\right)\cdot\mathcal{O}\left(\log
k+\log\log\Delta\right)\cdot\|\mu-\nu\|_{\operatorname*{EMD}}\right]\geq
0.999.$

∎

In summary, we have the following bounds on the contraction and the dilation
of the EMD embedding into $L_{1}$.

######  Theorem 2.9.

For every $\mu,\nu\in\operatorname*{EMD}_{[\Delta]^{d}}$ and every
$s=(s_{1},\ldots,s_{d})\in\mathbb{Z}^{d}$, it holds that

$\|\mu-\nu\|_{\operatorname*{EMD}}\leq\mathcal{O}\left(\sqrt{d}\right)\cdot\|G_{s}(\mu-\nu)\|_{1}.$

In addition, if $\mu$ and $\nu$ are two probability measures over
$[\Delta]^{d}$ and $\mu$ has support of size at most $k$. Then

$\mathbf{Pr}\left[\|G_{s}(\mu-\nu)\|_{1}\leq\mathcal{O}\left(d\right)\cdot\mathcal{O}\left(\log
k+\log\log\Delta\right)\cdot\|\mu-\nu\|_{\operatorname*{EMD}}\right]\geq
0.999,$

###### Proof.

The first part follows from Lemma 2.2 and the second part follows from Lemma
2.8. ∎

We now require simple generalizations of statements from [BIRW16] to the
Wasserstein-$z$ distance, as follows:

######  Lemma 2.10 (Claim 5.1 in [BIRW16]).

Let $\mu$ be a fixed $k$-sparse probability measure over $[\Delta]^{d}$ and
let $R>0$ be a fixed radius. For every $k$-sparse $\mu^{\prime}$ such that
$\operatorname*{WASSD}(\mu,\mu^{\prime})\leq R$, there exists an optimal flow
between $\mu$ and $\mu^{\prime}$ that is supported on at most $2k$ pairs of
points.

Algorithm 1 $\operatorname*{WASSD}$ Net Construction

1:Radius $R>0$, parameters $d$, $\Delta$, $k$

2:Net ${\mathcal{M}}$ on $k$-sparse probability measures with respect to
$\operatorname*{WASSD}$

3:$m_{0}\leftarrow\frac{R}{100d\Delta k}$, ${\mathcal{M}}\leftarrow\emptyset$

4:for $c:\operatorname*{supp}\mu\to\mathbb{Z}_{>0}$ such that
$\sum_{(x_{1},\ldots,x_{d})\in\operatorname*{supp}\mu}c(x_{1},\ldots,x_{d})\leq
2k$ do

5:
${\mathcal{I}}\leftarrow\\{(i,x_{1},\ldots,x_{d})\,\mid\,(x_{1},\ldots,x_{d})\in\operatorname*{supp}\mu,i\in[c(x_{1},\ldots,x_{d})]$

6: for $\ell:{\mathcal{I}}\to\\{1,1.01,1.01^{2},\ldots,2\Delta\\}$ do

7: for $(i,x_{1},\ldots,x_{d})\in{\mathcal{I}}$ and for all
$p(i,x_{1},\ldots,x_{d})\in{\textsc{BuildNet}}((x_{1},\ldots,x_{d}),\ell(i,x_{1},\ldots,x_{d}))$
do

8: for
$m:{\mathcal{I}}\to\\{0,m_{0},1.01m_{0},1.01^{2}m_{0},\ldots,\min(1,R)\\}$ do

9: if for every $(x_{1},\ldots,x_{d})\in\operatorname*{supp}\mu$, it holds
that
$\sum_{i:(i,x_{1},\ldots,x_{d})\in{\mathcal{I}}}m(i,x_{1},\ldots,x_{d})\leq\mu(x_{1},\ldots,x_{d})$
then

10: Let $\mu^{\prime}$ be a measure over $[\Delta]^{d}$ that is identically
zero

11: for $(x_{1},\ldots,x_{d})\in\operatorname*{supp}\mu$ do

12: $s\leftarrow 0$

13: for $i:(i,x_{1},\ldots,x_{d})\in{\mathcal{I}}$ do

14: $s\leftarrow s+m(i,x_{1},\ldots,x_{d})$

15:
$\mu^{\prime}(p(i,x_{1},\ldots,x_{d}))\leftarrow\mu^{\prime}(p(i,x_{1},\ldots,x_{d}))+m(i,x_{1},\ldots,x_{d})$

16:
$\mu^{\prime}(x_{1},\ldots,x_{d})\leftarrow\mu^{\prime}(x_{1},\ldots,x_{d})+\mu(x_{1},\ldots,x_{d})-s$

17: ${\mathcal{M}}\leftarrow{\mathcal{M}}\cup\\{\mu^{\prime}\\}$

######  Lemma 2.11.

The doubling dimension of the set of $k$-sparse probability measures over
$[\Delta]^{d}$ under $\operatorname*{WASSD}$ is
$\mathcal{O}\left(dk\log\log\Delta\right)$.

###### Proof.

Let $\mu$ be a fixed $k$-sparse probability measure over $[\Delta]^{d}$ and
let $R>0$ be a fixed radius. Let $B_{\operatorname*{WASSD}}(\mu,R)$ denote the
set of points $\nu$ such that $\operatorname*{WASSD}(\mu,\nu)\leq R$. We show
that $B_{\operatorname*{WASSD}}(\mu,R)$ can be covered with
$(\log\Delta)^{\mathcal{O}\left(kd\right)}$ $\operatorname*{WASSD}$-balls with
radius $\frac{R}{2}$ that are centered at $k$-sparse measures. Note that it
suffices to cover $B_{\operatorname*{WASSD}}(\mu,R)$ with
$(\log\Delta)^{\mathcal{O}\left(kd\right)}$ $\operatorname*{WASSD}$-balls with
radius $\frac{R}{4}$ that are centered at arbitrary probability measures.

In Algorithm 1, we construct a set of measures ${\mathcal{M}}$ that form the
centers of balls with radius $\frac{R}{4}$ whose union together covers
$B_{\operatorname*{WASSD}}(\mu,R)$. The high-level approach is to first
enumerate over all possible topologies of the optimal flow, then enumerate
over all possible lengths of the corresponding edges, then enumerate over all
possible supports, and finally, enumerate over all possible masses that are
transported over the edges.

We first assume access to a subroutine ${\textsc{BuildNet}}(p,r)$ that returns
a $\frac{r}{100}$-net of $B_{L_{z}}(p,r)\cap[\Delta]^{d}$, i.e., a net over
the $L_{z}$ ball of radius $r$ centered that $p$ restricted to the points in
$[\Delta]^{d}$. It follows that
$|{\mathcal{M}}|\leq(\log\Delta)^{\mathcal{O}\left(dk\right)}$ and thus, the
runtime of Algorithm 1 is also at most
$(\log\Delta)^{\mathcal{O}\left(dk\right)}$. It remains to show that for every
$k$-sparse probability measure $\mu^{\prime}$ such that
$\operatorname*{WASSD}(\mu,\mu^{\prime})\leq R$, there exists
$\mu^{\prime\prime}\in{\mathcal{M}}$ with
$\operatorname*{WASSD}(\mu^{\prime\prime},\mu^{\prime})\leq\frac{R}{4}$.

By Lemma 2.10, there is an optimal flow between $\mu$ and $\mu^{\prime}$
containing at most $2k$ edges. Thus, in the outer for loop of Algorithm 1,
there exists at least one guess of $c(x_{1},\ldots,x_{d})$ that indeed
corresponds to the number of outgoing edges from $(x_{1},\ldots,x_{d})$ in the
optimal flow. In the second for loop of Algorithm 1 that guesses the flow for
each edge, there is at least one choice that guesses all flows within a
multiplicative factor of $1.01$. Hence, there exists a measure $\tilde{\mu}$
such that

  1. (1)

$\operatorname*{supp}(\tilde{\mu})\in\operatorname*{supp}(\mu)\cup\\{p(i,x_{1},\ldots,x_{d})\\}_{(i,x_{1},\ldots,x_{d})\in{\mathcal{I}}}$

  2. (2)

$\operatorname*{WASSD}(\mu^{\prime},\tilde{\mu})\leq\frac{R}{50}$

  3. (3)

There exists a flow between $\mu$ and $\tilde{\mu}$ with cost at most $1.02R$
and transports mass from each
$(x_{1},\ldots,x_{d})\in\operatorname*{supp}(\mu)$ to some point
$\\{(x_{1},\ldots,x_{d})\\}\cup\\{p(i,x_{1},\ldots,x_{d})\\}_{(i,x_{1},\ldots,x_{d})\in{\mathcal{I}}}$.

Consequently, Algorithm 1 will guess $\operatorname*{supp}(\tilde{\mu})$ and
then the measure at the support. We show that some guess $\mu^{\prime\prime}$
by Algorithm 1 will satisfy
$\operatorname*{supp}(\mu^{\prime\prime})\subseteq\operatorname*{supp}(\tilde{\mu})$
and $\operatorname*{WASSD}(\tilde{\mu},\mu^{\prime\prime})\leq\frac{R}{25}$.

To that end, we first define $\mu^{\prime\prime}$. Consider a fixed
$(x_{1},\ldots,x_{d})\in\operatorname*{supp}(\mu)$ and the corresponding
multi-set
$\\{(x_{1},\ldots,x_{d})\\}\cup\\{p(i,x_{1},\ldots,x_{d})\,\mid\,(i,x_{1},\ldots,x_{d})\in{\mathcal{I}}\\}$.
Round down the mass of $\tilde{\mu}$ at the coordinates
$\\{p(i,x_{1},\ldots,x_{d})\\}$ to the closest element of
$\\{0,m_{0},1.01m_{0},1.01^{2}m_{0},\ldots,\min(1,R)\\}$ and let
$\mu^{\prime\prime}$ be the resulting measure. We set
$\mu^{\prime\prime}((x_{1},\ldots,x_{d}))=\sum_{i:(i,x_{1},\ldots,x_{d})\in{\mathcal{I}})}(\tilde{\mu}(p(i,x_{1},\ldots,x_{d}))-\mu^{\prime\prime}(p(i,x_{1},\ldots,x_{d})))$
and emphasize that $\mu^{\prime\prime}$ is one of the measures guessed by
Algorithm 1.

We now show
$\operatorname*{WASSD}(\tilde{\mu},\mu^{\prime\prime})\leq\frac{R}{25}$, by
upper bounding the following two terms:

  1. (1)

We first upper bound the contribution from
$(i,x_{1},\ldots,x_{d})\in{\mathcal{I}}$ for which
$\tilde{\mu}(p(i,x_{1},\ldots,x_{d}))<m_{0}$, so that
$\mu^{\prime\prime}(p(i,x_{1},\ldots,x_{d}))=0$. There exist at most $2k$
elements $(i,x_{1},\ldots,x_{d})\in{\mathcal{I}}$ and these elements can be
rerouted with cost at most $2\Delta\sqrt{d}km_{0}\leq\frac{R}{50}$.

  2. (2)

We next upper bound the contribution from
$(i,x_{1},\ldots,x_{d})\in{\mathcal{I}}$ for which
$\tilde{\mu}(p(i,x_{1},\ldots,x_{d}))\geq m_{0}$, so that
$\mu^{\prime\prime}(p(i,x_{1},\ldots,x_{d}))=0$ is within a multiplicative
$1.01$ factor of $\tilde{\mu}(p(i,x_{1},\ldots,x_{d}))=0$. Hence, the total
contribution of such elements is at most
$0.01\operatorname*{WASSD}(\mu,\tilde{\mu})\leq\frac{R}{50}$.

Therefore,
$\operatorname*{WASSD}(\mu^{\prime}-\mu^{\prime\prime})\leq\operatorname*{WASSD}(\mu^{\prime}-\tilde{\mu})+\operatorname*{WASSD}(\tilde{\mu}-\mu^{\prime\prime})\leq\frac{R}{50}+\frac{R}{50}+\frac{R}{50}<\frac{R}{4}$.
∎

We can use the results for the doubling dimension of the set of $k$-sparse
Wasserstein-$z$ (and in particular, earth mover distance for $z=1$)
probability measures as follows.

######  Lemma 2.12 (Lemma 3.1 in [BIRW16]).

Let ${\mathcal{M}}=(X,\rho)$ be a $K$-quasi-metric space and $Y\subseteq X$
have doubling dimension $d$. Suppose there exists a sketch of size $s$ that
approximates distances between points of $X$ and $Y$ with distortion $D$ with
probability at least $\frac{2}{3}$. Then for every
$\varepsilon\in\left(0,\frac{1}{2}\right)$, $\lambda\in(0,\Lambda)$, and
$y_{0}\in Y$, one can sketch points of $X$ with sketch size
$\mathcal{O}\left(s(d\log(DK/\varepsilon)+\log\log(\Lambda/\lambda))\right)$,
so that from this sketch for $x\in X$ with $\rho(x,y_{0})\leq\Lambda$, with
probability at least $\frac{2}{3}$, we can recover a point $y^{\prime}\in Y$
such that

$\rho(x,y^{\prime})\leq\max((1+\varepsilon)DK\cdot\rho(x,Y),\lambda).$

Lemma 2.12 avoids a brute force search over the entire metric space by
cleverly utilizing the doubling dimension to perform an efficient top-down
search. In particular, it first partitions the search space into balls of
radius $R$ and identifies the ball $B$ that contains the smallest distance
from $X$. It then partitions $B$ into smaller balls of radius $\frac{R}{2}$
and iterates to identify $y^{\prime}$.

##  3 Online Sensitivity Sampling and Insertion-Only Streams

In this section, we present our simple algorithm for $(k,z)$-clustering in the
insertion-only model. Along the way, we first show in Section 3.1 important
properties of online sensitivity sampling that we believe could be of
independent interest. Namely, we show that online sensitivity sampling gives a
$(1+\varepsilon)$-coreset to the underlying dataset for $(k,z)$-clustering. We
also upper bound the sum of the online sensitivities, which allows us to upper
bound the number of samples procured by online sensitivity sampling. Finally,
we give our algorithm for $(k,z)$-clustering on insertion-only streams in
Section 3.2.

###  3.1 Online Sensitivity Sampling

In this section, we show that online sensitivity sampling can be used to
achieve a $(1+\varepsilon)$-coreset for $(k,z)$-clustering. We show that an
approximately optimal clustering can be used to compute approximations of the
online sensitivities. We remark that the same statement can also be used to
compute approximations of the sensitivities of query points.

######  Lemma 3.1.

Let $S\subset\mathbb{R}^{d}$ with $|S|=k$ satisfy
$\operatorname*{cost}(X,S)\leq\gamma{\text{OPT}}$ for some $\gamma\geq 1$,
where OPT denotes the cost of an optimal $(k,z)$-clustering. Let
$Z\leq\operatorname*{cost}(X,S)\leq\beta Z$. Let $\pi:X\to S$ be the mapping
from $X$ to $S$ that induces the clustering cost $\operatorname*{cost}(X,S)$
and let
$Z^{\prime}\leq\operatorname*{cost}(C,\pi(X))\leq\beta^{\prime}Z^{\prime}$.
Then for any $C\subset\mathbb{R}^{d}$ with $|C|=k$,

$\frac{\operatorname*{cost}(x,C)}{\operatorname*{cost}(X,C)}\leq\frac{2^{z}\cdot
4\gamma\operatorname*{cost}(x,C)}{Z^{\prime}+Z}\leq\frac{2^{z}\cdot
4\beta\beta^{\prime}\gamma\operatorname*{cost}(x,C)}{\operatorname*{cost}(X,C)}.$

###### Proof.

Let $S\subset\mathbb{R}^{d}$ with $|S|=k$ satisfy
$\operatorname*{cost}(X,S)\leq\gamma{\text{OPT}}$ for some $\gamma\geq 1$. Let
$\pi:X\to S$ be the mapping from $X$ to $S$ that induces the clustering cost
$\operatorname*{cost}(X,S)$. Since
$\gamma\operatorname*{cost}(X,C)\geq\gamma{\text{OPT}}\geq\operatorname*{cost}(X,S)$,

$\displaystyle\frac{\operatorname*{cost}(x,C)}{\operatorname*{cost}(X,C)}$
$\displaystyle=\frac{4\gamma\operatorname*{cost}(x,C)}{4\gamma\operatorname*{cost}(X,C)}=\frac{4\gamma\operatorname*{cost}(x,C)}{2\gamma\operatorname*{cost}(X,C)+2\gamma\operatorname*{cost}(X,C)}\leq\frac{4\gamma\operatorname*{cost}(x,C)}{2\gamma\operatorname*{cost}(X,C)+2\operatorname*{cost}(X,S)}$

Note that $\operatorname*{cost}(X,S)=\operatorname*{cost}(X,\pi(X))$, so that

$\displaystyle\frac{\operatorname*{cost}(x,C)}{\operatorname*{cost}(X,C)}$
$\displaystyle\leq\frac{4\gamma\operatorname*{cost}(x,C)}{2\gamma\operatorname*{cost}(X,C)+2\operatorname*{cost}(X,\pi(X))}$
$\displaystyle\leq\frac{4\gamma\operatorname*{cost}(x,C)}{\operatorname*{cost}(X,C)+\operatorname*{cost}(X,\pi(X))+\operatorname*{cost}(X,\pi(X))},$

since $\gamma\geq 1$. By the generalized triangle inequality, i.e., Fact 1.9,
we then have
$(\operatorname*{cost}(X,C)+\operatorname*{cost}(X,\pi(X))\geq\frac{1}{2^{z}}\operatorname*{cost}(C,\pi(X))$,
so that

$\displaystyle\frac{\operatorname*{cost}(x,C)}{\operatorname*{cost}(X,C)}$
$\displaystyle\leq\frac{4\gamma\operatorname*{cost}(x,C)}{\frac{1}{2^{z}}\operatorname*{cost}(C,\pi(X))+\operatorname*{cost}(X,\pi(X))}$
$\displaystyle\leq\frac{2^{z}\cdot
4\gamma\operatorname*{cost}(x,C)}{\operatorname*{cost}(C,\pi(X))+\operatorname*{cost}(X,\pi(X))}\leq\frac{2^{z}\cdot
4\gamma\operatorname*{cost}(x,C)}{\operatorname*{cost}(X,C)}.$

Thus, for $\operatorname*{cost}(X,S)=\operatorname*{cost}(X,\pi(X)$ and
$Z\leq\operatorname*{cost}(X,S)\leq\beta Z$, we have

$\displaystyle\frac{\operatorname*{cost}(x,C)}{\operatorname*{cost}(X,C)}$
$\displaystyle\leq\frac{2^{z}\cdot
4\gamma\operatorname*{cost}(x,C)}{\operatorname*{cost}(C,\pi(X))+Z}\leq\frac{2^{z}\cdot
4\beta\gamma\operatorname*{cost}(x,C)}{\operatorname*{cost}(X,C)}.$

Finally, since
$Z^{\prime}\leq\operatorname*{cost}(C,\pi(X))\leq\beta^{\prime}Z^{\prime}$,
then we have

$\frac{\operatorname*{cost}(x,C)}{\operatorname*{cost}(X,C)}\leq\frac{2^{z}\cdot
4\gamma\operatorname*{cost}(x,C)}{Z^{\prime}+Z}\leq\frac{2^{z}\cdot
4\beta\beta^{\prime}\gamma\operatorname*{cost}(x,C)}{\operatorname*{cost}(X,C)}.$

∎

We next bound the sum of the online sensitivities for $(k,z)$-clustering. See
1.3

###### Proof.

For each $t\in[n]$, we define $X_{t}=x_{1},\ldots,x_{t}$ and
${\text{OPT}}(X_{t})$ to be the objective of an optimal $(k,z)$ clustering of
$X_{t}$. Let $t_{0}=1,t_{1},\ldots,t_{y}$ be a subsequence of $1,\ldots,n$
such that ${\text{OPT}}(X_{t_{i}})>2{\text{OPT}}(X_{t_{i-1}})$ and either
$t_{i}-1=t_{i-1}$ or ${\text{OPT}}(X_{t_{i}-1})\leq
2{\text{OPT}}(X_{t_{i-1}})$ for all $i\in[y]$. Note that
$y=\mathcal{O}\left(\log(nd\Delta)\right)$ for constant $z$ since each of the
$n$ points can contribute at most $d\cdot\Delta^{z}$ cost. Hence the sum of
the online sensitivities at times $t_{i}$ where $t_{i}-1=t_{i-1}$ is at most
$\mathcal{O}\left(\log(nd\Delta)\right)$. We consider the sum of online
sensitivities at the remaining times.

Let $i\in[y]$ be fixed and consider a set $K_{t_{i}}$ of $k$ centers such that
${\text{cost}}(X_{t_{i}},K_{t_{i}})={\text{OPT}}(X_{t_{i}})$. Let
$\pi:X_{i}\to K_{i}$ be a mapping from each point of $X_{t_{i}}$ to its
closest center in $K_{t_{i}}$ so that $\pi(x)=\operatorname*{argmin}_{p\in
K_{t_{i}}}{\text{dist}}(x,p)$ for all $x\in X_{t_{i}}$. By the generalized
triangle inequality, i.e., Fact 1.9, we have that for $t\in(t_{i-1},t_{i}]$
and for any set of $C$ centers, i.e., $C\subseteq[\Delta]^{d}$ with $|C|=k$,

$\frac{{\text{dist}}^{z}(x_{t},C)}{\sum_{j=1}^{t}{\text{dist}}^{z}(x_{j},C)}\leq\frac{2^{z-1}{\text{dist}}^{z}(x_{t},\pi(x_{t}))}{\sum_{j=1}^{t}{\text{dist}}^{z}(x_{j},C)}+\frac{2^{z-1}{\text{dist}}^{z}(\pi(x_{t}),C)}{\sum_{j=1}^{t}{\text{dist}}^{z}(x_{j},C)}.$

We upper bound the first term of the right-hand side by lower bounding the
denominator. Due to the optimality of $K_{t_{i}}$, we have

$\sum_{j=1}^{t}{\text{dist}}^{z}(x_{j},C)\geq\sum_{j=1}^{t_{i-1}}{\text{dist}}^{z}(x_{j},C)\geq\sum_{j=1}^{t_{i-1}}{\text{dist}}^{z}(x_{j},K_{t_{i-1}})={\text{OPT}}(X_{t_{i-1}})>\frac{1}{2}{\text{OPT}}(X_{t_{i}}).$

Since $t\in(t_{i-1},t_{i}]$ and
${\text{OPT}}(X_{t_{i}})=\sum_{j=1}^{t_{i}}{\text{dist}}^{z}(x_{j},\pi(x_{j}))\geq\sum_{j=1}^{t}{\text{dist}}^{z}(x_{j},\pi(x_{j}))$,
then

$\frac{2^{z-1}{\text{dist}}^{z}(x_{t},\pi(x_{t}))}{\sum_{j=1}^{t}{\text{dist}}^{z}(x_{j},C)}\leq\frac{2^{z}{\text{dist}}^{z}(x_{t},\pi(x_{t}))}{\sum_{j=1}^{t}{\text{dist}}^{z}(x_{j},\pi(x_{j}))}.$

To bound the second term of the right-hand side, let $S_{t}$ denote the subset
of $X_{t}$ mapped to $\pi(x_{t})$ at time $t_{i}$, i.e.,
$S_{t}=\pi^{-1}(\pi(x_{t}))\cap X_{t}$. Then

${\text{dist}}^{z}(\pi(x_{t}),C)\cdot|S_{t}|=\sum_{p\in
S_{t}}{\text{dist}}^{z}(\pi(x_{t}),C).$

By the generalized triangle inequality, i.e., Fact 1.9,

$\displaystyle{\text{dist}}^{z}(\pi(x_{t}),C)\cdot|S_{t}|$
$\displaystyle=\sum_{p\in S_{t}}{\text{dist}}^{z}(\pi(x_{t}),C)$
$\displaystyle\leq 2^{z-1}\sum_{p\in
S_{t}}\left[{\text{dist}}^{z}(\pi(p),p)+{\text{dist}}^{z}(p,C)\right]$
$\displaystyle\leq 2^{z-1}\sum_{p\in
X_{t}}\left[{\text{dist}}^{z}(\pi(p),p)+{\text{dist}}^{z}(p,C)\right],$

since $S_{t}=\pi^{-1}(\pi(x_{t}))\cap X_{t}$ is a subset of $X_{t}$. Because
$X_{t}\subseteq X_{t_{i}}$, then

$\sum_{p\in X_{t}}{\text{dist}}^{z}(\pi(p),p)\leq\sum_{p\in
X_{t_{i}}}{\text{dist}}^{z}(\pi(p),p)={\text{OPT}}(X_{t_{i}})\leq
2{\text{OPT}}(X_{t}).$

By the optimality of ${\text{OPT}}(X_{t})$,

$\sum_{p\in X_{t}}{\text{dist}}^{z}(\pi(p),p)\leq 2\sum_{p\in
X_{t}}{\text{dist}}^{z}(p,C)=2\sum_{j=1}^{t}{\text{dist}}^{z}(x_{j},C).$

Hence,

${\text{dist}}^{z}(\pi(x_{t}),C)\cdot|S_{t}|\leq 3\cdot
2^{z-1}\sum_{j=1}^{t}{\text{dist}}^{z}(x_{j},C),$

so that

$\frac{{\text{dist}}^{z}(\pi(x_{t}),C)}{\sum_{j=1}^{t}{\text{dist}}^{z}(x_{j},C)}\leq\frac{3\cdot
2^{z-1}}{|S_{t}|}.$

Putting things together,

$\displaystyle\frac{{\text{dist}}^{z}(x_{t},C)}{\sum_{j=1}^{t}{\text{dist}}^{z}(x_{j},C)}$
$\displaystyle\leq\frac{2^{z-1}{\text{dist}}^{z}(x_{t},\pi(x_{t}))}{\sum_{j=1}^{t}{\text{dist}}^{z}(x_{j},C)}+\frac{2^{z-1}{\text{dist}}^{z}(\pi(x_{t}),C)}{\sum_{j=1}^{t}{\text{dist}}^{z}(x_{j},C)}$
$\displaystyle\leq\frac{2^{z}{\text{dist}}^{z}(x_{t},\pi(x_{t}))}{\sum_{j=1}^{t}{\text{dist}}^{z}(x_{j},\pi(x_{j}))}+\frac{3\cdot
2^{2z-1}}{|S_{t}|}.$

Therefore,

$\displaystyle\sum_{t=t_{i-1}+1}^{t_{i}}\sigma_{t}$
$\displaystyle\leq\sum_{t=t_{i-1}+1}^{t_{i}}\left(\frac{2^{z}{\text{dist}}^{z}(x_{t},\pi(x_{t}))}{\sum_{j=1}^{t}{\text{dist}}^{z}(x_{j},\pi(x_{j}))}+\frac{3\cdot
2^{2z-1}}{|S_{t}|}\right)$
$\displaystyle=\sum_{t=t_{i-1}+1}^{t_{i}}\frac{2^{z}{\text{dist}}^{z}(x_{t},\pi(x_{t}))}{\sum_{j=1}^{t}{\text{dist}}^{z}(x_{j},\pi(x_{j}))}+\sum_{t=t_{i-1}+1}^{t_{i}}\frac{3\cdot
2^{2z-1}}{|S_{t}|}$ $\displaystyle\leq
2^{z}+\sum_{t=t_{i-1}+1}^{t_{i}}\frac{3\cdot 2^{2z-1}}{|S_{t}|}.$

Since $S_{t}=\pi^{-1}(\pi(x_{t}))\cap X_{t}$ and $\pi^{-1}$ can map onto one
of $k$ different sets of points, then

$\sum_{t=t_{i-1}+1}^{t_{i}}\frac{1}{|S_{t}|}\leq k\sum_{t=1}^{n}\frac{1}{t}.$

Hence,

$\displaystyle\sum_{t=t_{i-1}+1}^{t_{i}}\sigma_{t}$ $\displaystyle\leq
2^{z}+(3k\cdot 2^{2z-1})\sum_{t=1}^{n}\frac{1}{t}=\mathcal{O}\left(2^{2z}k\log
n\right).$

Since the sequence $t_{0},\ldots,t_{y}$ satisfies
$y=\mathcal{O}\left(\log(nd\Delta)\right)$, then

$\sum_{t=1}^{n}\sigma_{t}=\sum_{i=1}^{y}\sum_{t=t_{i-1}+1}^{t_{i}}\sigma_{t}=\mathcal{O}\left(2^{2z}k\log^{2}(nd\Delta)\right).$

∎

While the input points are from $[\Delta]^{d}$, the optimal cluster centers
may not lie in $[\Delta]^{d}$. Although there is an arbitrary number of
subsets of $\mathbb{R}^{d}$ of size $k$, we now recall the well-known fact
that to achieve an approximately optimal clustering, it suffices to only show
correctness on a sufficiently-sized net. For completeness, we include the
proof here.

######  Lemma 3.2.

Let $X\subset[\Delta]^{d}$ and let $z\geq 1$ be a constant. Then there exists
a set $S$ of size
$|S|=\left(\frac{n\Delta}{\varepsilon}\right)^{\mathcal{O}\left(kd\right)}$,
such that
$(1-\varepsilon)\operatorname*{cost}(A,C)\leq\operatorname*{cost}(X,C)\leq(1+\varepsilon)\operatorname*{cost}(A,C)$
for any $C\in S$, implies
$(1-\varepsilon)\operatorname*{cost}(A,C)\leq\operatorname*{cost}(X,C)\leq(1+\varepsilon)\operatorname*{cost}(A,C)$
for any set $C\subset\mathbb{R}^{d}$ with $|C|=k$.

###### Proof.

Note that since $X\subset[\Delta]^{d}$, then either the optimal
$(k,z)$-clustering is zero or at least $\Omega\left(\frac{1}{2^{z}}\right)$.
Moreover, note that moving a point or a center by a distance $D$ can only
change the clustering cost by $n(\Delta d)^{z}$. Hence, it suffices to create
a $D$-net of $[\Delta]^{d}$ for
$D=\mathcal{O}\left(\frac{\varepsilon}{n(\Delta d)^{z}}\right)$, which has
size $\left(\frac{n\Delta}{\varepsilon}\right)^{\mathcal{O}\left(d\right)}$
for constant $z\geq 1$. The claim then follows from the observation that it
suffices to choose $k$ points from this net for any set of $k$ queries. ∎

We now show the correctness of online sensitivity sampling for
$(k,z)$-clustering.

######  Lemma 3.3.

Let ${\mathcal{P}}$ be a process that for each $t\in[n]$, samples each point
$x_{t}$ with probability $p_{t}\geq\min(1,\gamma\cdot\sigma_{t})$, where
$\sigma_{t}$ is the online sensitivity of $x_{t}$ and
$\gamma=\mathcal{O}\left(\frac{dk}{\varepsilon^{2}}\log\frac{n\Delta}{\varepsilon}\right)$
and gives the point weight $\frac{1}{p_{t}}$ if $x_{t}$ is sampled. Then with
high probability, i.e., $1-\frac{1}{\operatorname*{poly}(n)}$, ${\mathcal{P}}$
is a $(1+\varepsilon)$-coreset for $(k,z)$-clustering.

###### Proof.

Let $C$ be a set of $k$ centers. Let $x_{1},\ldots,x_{n}\in\mathbb{R}^{d}$ be
a sequence of points and for $t\in[n]$, let $X_{t}=\\{x_{1},\ldots,x_{t}\\}$.
Let $Y_{0}=0,Y_{1},\ldots,Y_{n}$ be a martingale with difference sequence
$Z_{0},\ldots,Z_{n}$, where for $t\geq 1$, we define $Z_{t}=0$ if
$|Y_{t-1}|>\varepsilon\operatorname*{cost}(C,X_{t-1})$. Otherwise, for
$|Y_{t-1}|\leq\varepsilon\operatorname*{cost}(C,X_{t-1})$, we define

$X_{t}=\begin{cases}\left(\frac{1}{p_{t}}-1\right)\operatorname*{cost}(C,x_{t})\qquad&\text{if
}x_{t}\text{ is sampled}\\\
-\operatorname*{cost}(C,x_{t})\qquad&\text{otherwise}\end{cases}.$

Note that $\mathbb{E}\left[X_{t}\right]=0$ and thus $\mathbb{E}\left[Y_{t}\mid
Y_{1},\ldots,Y_{t-1}\right]=Y_{t-1}$, so that $Y_{0},\ldots,Y_{t}$ is a valid
martingale sequence. Furthermore, by the construction of $Y_{t}$, we have that
$Y_{t}=\operatorname*{cost}(C,A_{t})-\operatorname*{cost}(C,X_{t})$, where
$A_{t}$ is the coreset at time $t$.

Observe that if $p_{t}=1$, then $X_{t}=0$. Otherwise, we have that
$\mathbb{E}\left[X_{t}^{2}\right]\leq\frac{1}{p_{t}}\operatorname*{cost}(C,x_{t})^{2}$.
Let $\gamma=\frac{100}{\varepsilon^{2}}\sigma_{t}\log\frac{1}{\delta}$ for
some parameter $\delta$ that we will set, which will ultimately give the
claimed setting for $\gamma$. Since $p_{t}\geq\gamma\cdot\sigma_{t}$, then
$p_{t}\geq\frac{100\log\frac{1}{\delta}}{\varepsilon^{2}}\sigma_{t}$, then
$p_{t}\geq\frac{100\log\frac{1}{\delta}}{\varepsilon^{2}}\frac{\operatorname*{cost}(C,x_{t})}{\operatorname*{cost}(C,X_{t})}$,
so that
$\mathbb{E}\left[X_{t}^{2}\right]\leq\frac{\varepsilon^{2}}{100\log\frac{1}{\delta}}\operatorname*{cost}(C,x_{t})\cdot\operatorname*{cost}(C,X_{t})$.
Hence,

$\sum_{i=1}^{t}\mathbb{E}\left[X_{i}^{2}\right]\leq\sum_{i=1}^{t}\frac{\varepsilon^{2}}{100}\operatorname*{cost}(C,x_{i})\cdot\operatorname*{cost}(C,X_{i})\leq\frac{\varepsilon^{2}}{100\log\frac{1}{\delta}}\operatorname*{cost}(C,X_{t})^{2}.$

Furthermore, we have
$X_{i}\leq\frac{\varepsilon^{2}}{100\log\frac{1}{\delta}}\operatorname*{cost}(C,X_{i})$
for all $i\in[t]$. Then by Freedman’s inequality, i.e., Theorem 1.15,

$\mathbf{Pr}\left[\max_{i\in[t]}|Y_{i}|\geq\varepsilon\operatorname*{cost}(C,X_{t})\right]\leq
2\exp\left(-\frac{\varepsilon\operatorname*{cost}(C,X_{t})^{2}/2}{\frac{\varepsilon^{2}}{100\log\frac{1}{\delta}}\operatorname*{cost}(C,X_{t})^{2}+\frac{\varepsilon^{3}}{300\log\frac{1}{\delta}}\operatorname*{cost}(C,X_{t})^{2}}\right),$

and so in particular,

$\mathbf{Pr}\left[|\operatorname*{cost}(C,A_{t})-\operatorname*{cost}(C,X_{t})|\leq\varepsilon\operatorname*{cost}(C,X_{t})\right]\geq
1-\delta.$

By Lemma 3.2 for each time $t\in[n]$, it suffices to union bound over a net of
size $\left(\frac{n\Delta}{\varepsilon}\right)^{\mathcal{O}\left(kd\right)}$
to ensure that

$|\operatorname*{cost}(C,A_{t})-\operatorname*{cost}(C,X_{t})|\leq
2\varepsilon\operatorname*{cost}(C,X_{t}),$

for any set $C\subset\mathbb{R}^{d}$ with $|C|=k$. We further union bound over
all $t\in[n]$ to ensure that

$|\operatorname*{cost}(C,A_{t})-\operatorname*{cost}(C,X_{t})|\leq
2\varepsilon\operatorname*{cost}(C,X_{t}),$

for any set $C\subset\mathbb{R}^{d}$ with $|C|=k$ and for any $t\in[n]$. Thus
by setting
$\log\frac{1}{\delta}=\mathcal{O}\left(kd\log\frac{n\Delta}{\varepsilon}\right)$,
we have that with high probability, ${\mathcal{P}}$ is a coreset for
$(k,z)$-clustering at all times in the stream. ∎

######  Lemma 3.4.

Let
$\gamma=\mathcal{O}\left(\frac{dk}{\varepsilon^{2}}\log\frac{n\Delta}{\varepsilon}\right)$
and $\tau>1$ be some parameter. Let ${\mathcal{P}}$ be a process that for each
$t\in[n]$, samples each point $x_{t}$ with probability $p_{t}$, where

$p_{t}\leq\min(1,\tau\gamma\sigma_{t}),$

for the online sensitivity $\sigma_{t}$ of $x_{t}$. Then with high
probability, ${\mathcal{P}}$ contains at most $\mathcal{O}\left(\frac{\tau
dk^{2}}{\varepsilon^{2}}\log\frac{n\Delta}{\varepsilon}\right)$ points.

###### Proof.

For each $t\in[n]$, let $Y_{t}$ denote the indicator random variable for
whether $x_{t}$ is sampled by ${\mathcal{P}}$, i.e., $Y_{t}=1$ if $x_{t}$ is
sampled by ${\mathcal{P}}$ and $Y_{t}=0$ otherwise, so that
$Y:=\sum_{t=1}^{n}Y_{t}$ is the total number of sampled points. We have that
$\mathbb{E}\left[Y_{t}\right]=p_{t}$ for each $t\in[n]$ and thus

$\mathbb{E}\left[Y\right]=\sum_{t=1}^{n}\mathbb{E}\left[Y_{t}\right]=\sum_{t=1}^{n}p_{t}\leq\tau\gamma\cdot\sum_{t=1}^{n}\sigma_{t}.$

By Theorem 1.3, we have

$\sum_{t=1}^{n}\sigma_{t}=\mathcal{O}\left(2^{2z}k\log^{2}(nd\Delta)\right),$

so that

$\mathbb{E}\left[Y\right]\leq\tau\gamma\cdot\mathcal{O}\left(2^{2z}k\log^{2}(nd\Delta)\right).$

Since $p_{t}\in[0,1]$, we have $\mathbb{E}\left[Y_{t}^{2}\right]\leq p_{t}$
and so
$\sum_{t=1}^{n}\mathbb{E}\left[Y_{t}^{2}\right]\leq\tau\gamma\cdot\sum_{t=1}^{n}\sigma_{t}$.
Hence by Freedman’s inequality, i.e., Theorem 1.15, we have that for
$\gamma=\mathcal{O}\left(\frac{dk}{\varepsilon^{2}}\log\frac{n\Delta}{\varepsilon}\right)$,
$Y=\mathcal{O}\left(\frac{\tau
dk^{2}}{\varepsilon^{2}}\log\frac{n\Delta}{\varepsilon}\right)$ with high
probability. ∎

######  Theorem 3.5 (Online sensitivity sampling).

Let ${\mathcal{P}}$ be a process that for each $t\in[n]$, samples each point
$x_{t}$ with probability $p_{t}\geq\min(1,\gamma\cdot\sigma_{t})$, where
$\sigma_{t}$ is the online sensitivity of $x_{t}$ and
$\gamma=\mathcal{O}\left(\frac{dk}{\varepsilon^{2}}\log\frac{n\Delta}{\varepsilon}\right)$
and gives the point weight $\frac{1}{p_{t}}$ if $x_{t}$ is sampled. Then with
high probability, ${\mathcal{P}}$ is a coreset for $(k,z)$-clustering that
contains at most
$\mathcal{O}\left(\frac{dk^{2}}{\varepsilon^{2}}\log\frac{n\Delta}{\varepsilon}\right)$
points.

###### Proof.

We argue by induction at each time $t\in[n]$. At time $t=1$, the online
sensitivity of the first point is $1$ and thus the first point is sampled, and
so ${\mathcal{P}}$ is a coreset of $x_{1}$.

Now suppose that at some time $t\in[n]$, ${\mathcal{P}}$ is a coreset of
$X_{t}=\\{x_{1},\ldots,x_{t}\\}$. Then by the definition of coreset, the cost
of any set $C$ of $k$ centers will be preserved by ${\mathcal{P}}$. Hence, the
coreset will output a $(1+\varepsilon)$-approximation $q^{\prime}_{t+1}$ to
the online sensitivity of the next point $x_{x+1}$. Then by setting
$q_{t+1}=2q^{\prime}_{t+1}$, we have that $2\sigma_{t+1}>q_{t+1}>\sigma_{t+1}$
and so by setting $p_{t+1}=\gamma\cdot q_{t+1}$ for
$\gamma=\mathcal{O}\left(\frac{dk}{\varepsilon^{2}}\log\frac{n}{\varepsilon}\right)$,
the conditions of Lemma 3.3 are satisfied, so that with high probability,
${\mathcal{P}}$ remains a coreset at time $t+1$. It follows by induction that
${\mathcal{P}}$ is a coreset at all times.

Moreover, it follows that each point $x_{t}$ at time $t\in[n]$ is sampled with
probability at most $2\gamma\sigma_{t}$. Hence by Lemma 3.4, we have that with
high probability, ${\mathcal{P}}$ contains
$\mathcal{O}\left(\frac{dk^{2}}{\varepsilon^{2}}\log\frac{n\Delta}{\varepsilon}\right)$
points. ∎

###  3.2 $(k,z)$-Clustering on Insertion-Only Streams

In this section, we present a simple algorithm for $(k,z)$-clustering on
insertion-only streams using $o(\log n)$ words of space. The algorithm
proceeds as follows. At all times, we maintain a data structure that with high
probability, produces a $(1+\varepsilon)$-coreset $C_{t}$ to the point sets
$X_{t}=\\{x_{1},\ldots,x_{t}\\}$ at all times $t\in[n]$. Conditioned on the
correctness of the coreset $C_{t-1}$ at time $t-1$, we use $C_{t-1}$ to
approximate the online sensitivity of point $x_{t}$, which in turn, produces a
probability $p_{t}$ of sampling $x_{t}$ into a secondary stream
${\mathcal{S}}^{\prime}$. We do not maintain ${\mathcal{S}}^{\prime}$, but
rather feed it into a coreset construction algorithm, which then produces
$C_{t}$, completing the process for time $t$.

See 1.1

###### Proof.

For $t\in[n]$, let ${\mathcal{E}}_{t}$ be the event that the data structure
$C_{t}$ at time $t$ is a $(1+\varepsilon)$-coreset of the data set
$X_{t}=\\{x_{1},\ldots,x_{t}\\}$. For all $t\in[n]$, let
${\mathcal{S}}^{\prime}_{t}$ be the stream consisting of the weighted points
sampled by sensitivity sampling. Observe that conditioned on
${\mathcal{E}}_{t-1}$, we can compute a $(1+\varepsilon)$-approximation
$q_{t}$ to the online sensitivity $\sigma_{t}$ of the point $x_{t}$. Thus by
setting $p_{t}=\min(1,2\gamma q_{t})$ for
$\gamma=\mathcal{O}\left(\frac{dk}{\varepsilon^{2}}\log\frac{n\Delta}{\varepsilon}\right)$,
then by Lemma 3.3, we have that ${\mathcal{S}}^{\prime}_{t}$ is a
$(1+\varepsilon)$-coreset of $X_{t}$. Let ${\mathcal{F}}_{t}$ be the event
that ${\mathcal{S}}^{\prime}_{t}$ is a $(1+\varepsilon)$-coreset of $X_{t}$.
By applying a coreset construction with failure probability
$\frac{1}{\operatorname*{poly}(n)}$ and a rescaling of $\varepsilon$, it
follows that $C_{t}$ is a $(1+\varepsilon)$-coreset of $X_{t}$, conditioned on
${\mathcal{E}}_{t-1}$ and ${\mathcal{F}}_{t}$. By Lemma 3.4, we have that
$|{\mathcal{S}}^{\prime}_{t}|\leq\operatorname*{poly}\left(d,k,\frac{1}{\varepsilon},\log
n,\log\Delta\right)$. Hence, we have that
$\sum_{t=1}^{n}\mathbf{Pr}\left[\not{\mathcal{F}}_{t}\,\mid\,{\mathcal{E}}_{t-1}\right]\leq\frac{1}{\operatorname*{poly}(n)}$
and
$\sum_{t=1}^{n}\mathbf{Pr}\left[\not{\mathcal{E}}_{t}\,\mid\,{\mathcal{F}}_{t},{\mathcal{E}}_{t-1}\right]\leq\frac{1}{\operatorname*{poly}\left(d,k,\frac{1}{\varepsilon},\log
n,\log\Delta\right)}$. Thus by a union bound, we have that with high
probability, $C_{t}$ is a $(1+\varepsilon)$-coreset for $X_{t}$ at all times
$t\in[n]$. Finally by Theorem 1.11, we have that our algorithm uses
$\tilde{\mathcal{O}}\left(\frac{dk}{\varepsilon^{2}}\right)\cdot(2^{z\log
z})\cdot\min\left(\frac{1}{\varepsilon^{z}},k\right)\cdot\operatorname*{poly}(\log\log(n\Delta))$
words of memory. ∎

##  4 Two-Pass Dynamic Streaming Algorithm for $k$-Median Clustering

In this section, we present a two-pass algorithm that uses $o(\log n)$ words
of space for $k$-median clustering, using the techniques that we have built up
over the previous sections. Our algorithm outputs a weighted subset of the
input points, which forms a $(1+\varepsilon)$-coreset of the underlying point
set. In Section 6, we show it is impossible to achieve this in $o(\log n)$
words of space using a single pass over the data.

We first require preliminaries on $L_{1}$ norm estimation in an insertion-only
stream and on the notion of consistent clustering. We then incorporate these
ideas into our analysis along with structural results regarding sensitivity
sampling similar to those that were given in the previous section.

###  4.1 $L_{1}$ Sketch

We start by describing a streaming algorithm for $L_{1}$-approximation. In
this setting, the input is a stream of elements $u_{1},\ldots,u_{m}\in[n]$
that implicitly defines an underlying frequency vector $f\in\mathbb{R}^{n}$,
so that each update $u_{t}=(c_{t},\Delta_{t})$ changes coordinate
$c_{t}\in[n]$ by $\Delta_{t}\in[-M,M]$ for $M=\operatorname*{poly}(n)$. More
formally, for all $i\in[n]$,

$f_{i}=|\sum_{t\in[m]:c_{t}=i}\Delta_{t}.$

The problem is then to output a $(1+\varepsilon)$-approximation to

$\|f\|_{1}=\sum_{i=1}^{n}|f_{i}|=|f_{1}|+\ldots+|f_{n}|.$

We define the general family of $p$-stable distributions, though we shall
ultimately only require Cauchy random variables, i.e., the $p$-stable
distribution with $p=1$.

######  Definition 4.1 ($p$-stable distribution).

[Zol89] For $0<p\leq 2$, there exists a distribution ${\mathcal{D}}_{p}$
called the $p$-stable distribution ${\mathcal{D}}_{p}$, such that for any
positive integer $n$ with $Z_{1},\ldots,Z_{n}\sim{\mathcal{D}}_{p}$ and vector
$x\in\mathbb{R}^{n}$, then $\sum_{i=1}^{n}Z_{i}x_{i}\sim Z\cdot\|x\|_{p}$ for
$Z\sim{\mathcal{D}}_{p}$. The $p$-stable distribution has probability density
function $f(x)=\Theta\left(\frac{1}{1+|x|^{1+p}}\right)$ for $p\leq 2$ and is
equal to the Cauchy distribution for $p=1$ (and also the normal distribution
for $p=2$).

For completeness, [Nol03] presented a now standardized method for sampling a
random variable $X$ from the $p$-stable distribution, by generating $\theta$
uniformly at random from the interval
$\left[-\frac{\pi}{2},\frac{\pi}{2}\right]$, $r$ uniformly at random from the
interval $[0,1]$, and outputting

$X=f(r,\theta)=\frac{\sin(p\theta)}{\cos^{1/p}(\theta)}\cdot\left(\frac{\cos(\theta(1-p))}{\log\frac{1}{r}}\right)^{\frac{1}{p}-1}.$

The $p$-stable distribution is utilized in an $L_{1}$ estimation algorithm of
[Ind06]. The algorithm creates
$\ell=\mathcal{O}\left(\frac{1}{\varepsilon^{2}}\left(\log\frac{1}{\varepsilon\delta}+\log\log
m\right)\right)$ vectors $v^{(1)},\ldots,v^{(\ell)}\in\mathbb{R}^{n}$. For
each $i\in[\ell]$, the vector $v^{(i)}$ consists of independently generated
Cauchy random variables and the algorithm maintains $Z_{i}=\langle
v^{(i)},f\rangle$ throughout the stream. At the end of the stream, the
algorithm outputs $\operatorname*{median}_{i\in[\ell]}Z_{i}$ as the estimation
for $\|f\|_{1}$.

We have the following guarantees for the algorithm of $p$-stable sketch:

######  Theorem 4.2.

[Ind06] For $p\in(0,2]$ and any dynamic stream of length $m$ on a universe of
size $n$, there exists an algorithm that with probability at least $1-\delta$,
provides a $(1+\varepsilon)$-approximation to the $L_{p}$ norm of the stream,
using $\mathcal{O}\left(\frac{1}{\varepsilon^{2}}\log
nm\log\frac{1}{\delta}\right)$ bits of space.

Putting the $L_{1}$ sketch and the EMD embedding together, we have:

######  Lemma 4.3.

There exists an algorithm that uses $\mathcal{O}\left(kd(\log
k+\log\log\Delta)\right)$ words of space and outputs a clustering with
estimated cost $Z$ such that with probability at least $0.98$,

${\text{OPT}}\leq Z\leq\mathcal{O}\left(d^{1.5}(\log
k+\log\log\Delta)\right){\text{OPT}},$

where OPT is the optimal $k$-median clustering cost at the end of the stream.

###  4.2 Dynamic Streaming Algorithm

We now describe our two pass streaming algorithm for insertion-deletion
streams. In the first pass over the data stream, our algorithm maintains a
sketch of the EMD embedding, as in Algorithm 2.

Algorithm 2 EMD Sketch

1:Input points $x_{1},\ldots,x_{n}\in[\Delta]^{d}$

2:Data structure to support EMD queries

3:Initiate random shift $s$

4:Create grid embedding $G_{s}$

5:Let $\mu_{t}$ be a probability distribution and $\nu$ be a $k$-sparse
probability distribution

6:return a constant-factor approximation to $\|G_{s}(\mu_{t}-\nu)\|_{1}$ using
an $L_{1}$ sketch, i.e., Theorem 4.2

Because the data structure is a linear sketch, it is amenable to insertions
and deletions. We would like to use the EMD sketch to perform sensitivity
sampling over the points, where we define the sensitivity of a point
$x\in\mathbb{R}^{d}$ with respect to a set $X\subset\mathbb{R}^{d}$ by

$\varphi_{X}(x):=\max_{C\subset\mathbb{R}^{d}:|C|\leq
k}\frac{\operatorname*{cost}(x,C)}{\operatorname*{cost}(X,C)}=\max_{C\subset\mathbb{R}^{d}:|C|\leq
k}\frac{{\text{dist}}(x,C)^{z}}{\sum_{t=1}^{n}{\text{dist}}(x_{t},C)^{z}}.$

Throughout the section, the set $X$ will be the set of points at the end of
the stream and thus we omit the subscript $X$ for clarity. We show that that
an approximation to the sensitivity of each point can be computed in Algorithm
3.

Algorithm 3 Approximation of Sensitivity for $k$-Median Clustering

1:Input dataset $X\subset[\Delta]^{d}$

2:Approximate sensitivity for $x$ for any query $x\in[\Delta]^{d}$

3:$m\leftarrow\tilde{\mathcal{O}}\left(kd\log\log\Delta\right)$

4:Initialize an EMD sketch, i.e., $m$ instances of $L_{1}$ sketches on random
EMD embeddings

5:Update the EMD sketch with $X$ $\triangleright$ First pass

6:Use the EMD sketch to find a near-optimal set $S$ of $k$ centers

7:Let $\widetilde{\operatorname*{cost}(X,S)}$ be the estimated cost of $S$ on
$X$ $\triangleright$ Lemma 4.3

8:Let ${\mathcal{M}}$ be the net from Lemma 3.2

9:$q(x)\leftarrow 0$

10:for $C\in{\mathcal{M}}$ do

11:
$q(x)\leftarrow\max\left(q(x),\frac{8\gamma\operatorname*{cost}(x,C)}{\operatorname*{cost}(S,C)+\widetilde{\operatorname*{cost}(X,S)}}\right)$
$\triangleright$ Lemma 3.1

12:return $q(x)$

Unfortunately, sensitivity sampling is incompatible with the dynamic setting,
where sampled points could be subsequently deleted. Instead, in the second
pass over the data stream, our algorithm uses the EMD sketch to subsample from
the universe with probability proportional to each point’s sensitivity. Thus
if a universe element is subsampled, all updates to the element throughout the
stream are reported. The universe element is then fed into a sparse recovery
scheme, which is a linear sketch that can handle both insertions and
deletions. We present the algorithm in full in Algorithm 4.

Algorithm 4 Two-pass Dynamic Streaming Algorithm for $k$-Median

1:Stream of updates of length $m=\operatorname*{poly}(n)$ to coordinates of
$[\Delta]^{d}$, approximation parameter $\varepsilon\in(0,1)$, number of
clusters $k$, parameter $z\geq 1$

2:$(1+\varepsilon)$-coreset for $(k,z)$-median clustering

3:First pass: implicitly compute an approximate sensitivity $q(x)$ for all
$x\in[\Delta]^{d}$ $\triangleright$ Algorithm 3

4:$s\leftarrow\operatorname*{poly}\left(k,d,\frac{1}{\varepsilon^{2}},\log\log
n,\log\log\Delta\right)$

5:$p(x)=\min\left(1,\frac{k^{2}d}{\varepsilon^{2}}\log k\cdot q(x)\right)$

6:Initialize an $100s$-sparse recovery algorithm $A$ (see Lemma 4.5)

7:Draw hash function $h$ from a family of hash functions where for all
$x\in\mathbb{R}^{d}$, $\mathbf{Pr}\left[h(x)=1\right]=p(x)$ and $h(x)=0$
otherwise

8:for each update in the second pass do

9: if the update is to $x\in[\Delta]^{d}$ with $h(x)=1$ then

10: Update $A$ with the corresponding update to $x$

11:return the points $x$ output by $A$ weighted by $\frac{1}{q(x)}$

We first claim that with constant probability, we can obtain a bounded
approximation to the sensitivity of each point.

######  Lemma 4.4.

There exists an algorithm that uses
$\mathcal{O}\left(kd^{2}\log\log\Delta(\log k+\log\log n)\right)$ words of
space and outputs $q(x)$ for all points $x\in\mathbb{R}^{d}$ in the stream,
such that with probability at least $0.98$,

$\varphi(x)\leq q(x)\leq\mathcal{O}\left(d^{1.5}(\log
k+\log\log\Delta)\right)\cdot\varphi(x),$

where $\varphi(x)$ is the sensitivity of $x$ with respect to a point set
$X\subset\mathbb{R}^{d}$ for the $k$-median problem.

###### Proof.

Let $\widetilde{Z}$ be the smallest estimated cost by the EMD sketch scaled by
$\mathcal{O}\left(\sqrt{d}\right)$, across all sets of $k$ centers in the
doubling dimension net. Let OPT be the optimal $k$-median clustering cost of
the dataset at the end of the stream. Let ${\mathcal{E}}$ be the event that

${\text{OPT}}\leq\widetilde{Z}\leq\mathcal{O}\left(d^{1.5}(\log
k+\log\log\Delta)\right){\text{OPT}},$

so that $\mathbf{Pr}\left[{\mathcal{E}}\right]$ with probability at least
$0.98$, by Lemma 4.3.

Recall that for a point $x$, its sensitivity is defined by

$\varphi(x)=\max_{C:\,C\subset\mathbb{R}^{d},\,|C|=k}\frac{\operatorname*{cost}(x,C)}{\operatorname*{cost}(X,C)},$

where $X$ is the set of the points remaining at the end of the stream.

By Lemma 3.1, a constant factor approximation to a clustering that achieves a
clustering cost that is an $\mathcal{O}\left(d^{1.5}(\log
k+\log\log\Delta)\right)$-approximation to OPT can be used to compute a
$\mathcal{O}\left(d^{1.5}(\log k+\log\log\Delta)\right)$-approximation to
$\frac{\operatorname*{cost}(x,C)}{\operatorname*{cost}(X,C)}$ for any set $C$
of $k$ centers. Thus conditioned on ${\mathcal{E}}$ the estimate
$\widetilde{Z}$ of the EMD sketch can be used to compute a
$\mathcal{O}\left(d^{1.5}(\log k+\log\log\Delta)\right)$-approximation to
$\varphi(x)$ for all $x\in[\Delta]^{d}$. Hence, the EMD sketch can only be
used to compute $q(x)$ such that with probability at least $0.98$,

$\varphi(x)\leq q(x)\leq\mathcal{O}\left(d^{1.5}(\log
k+\log\log\Delta)\right)\cdot\varphi(x),$

simultaneously for all $t\in[n]$.

The space complexity results from maintaining
$\tilde{\mathcal{O}}\left(kd\log\log\Delta\right)$ instances of the EMD
sketch, i.e., the $L_{1}$ sketch on the EMD embedding. Note that since the
grid of the EMD embedding has size $\Delta^{\mathcal{O}\left(d\right)}$, then
each $L_{1}$ sketch uses $\mathcal{O}\left(d\log\Delta\right)$ bits of space.
Hence, the total space is
$\tilde{\mathcal{O}}\left(kd^{2}\log\log\Delta\right)$ words of space. ∎

We then recall a well-known deterministic algorithm for sparse recovery.

######  Lemma 4.5.

There exists a deterministic algorithm that recovers the non-zero coordinates
of a $k$-sparse vector of length $n$ in a dynamic stream, using
$\mathcal{O}\left(k\log n\right)$ bits of space.

We now show that if we perform sensitivity sampling on the universe, then the
number of sampled elements in the stream is upper bounded by $100s$, so that
the sparse recovery scheme used in Algorithm 4 will properly recover all
sampled points.

######  Lemma 4.6.

Let ${\mathcal{E}}$ be the event that for all $x\in X$, we have

$\varphi(x)\leq q(x)\leq\mathcal{O}\left(d^{1.5}(\log
k+\log\log\Delta)\right)\cdot\varphi(x).$

Then conditioned on ${\mathcal{E}}$, the set of nonzero points induced by the
input to $A$ has sparsity at most $100s$ with probability at least $0.99$.

###### Proof.

Although as many as $n$ points may be input to $A$, the nonzero points at the
end of the stream form exactly $X$. Because $A$ is a sparse recovery
algorithm, it suffices to consider the points $x$ of $X$ with $h(x)=1$.
Observe that conditioned on ${\mathcal{E}}$, the expected number of sampled
points is $\sum_{x\in X}p(x)\leq s$. Thus by Markov’s inequality, the set of
nonzero points induced by the input to $A$ has sparsity at most $100s$ with
probability at least $0.99$. ∎

We next recall the following statement showing that sensitivity sampling
results in a $(1+\varepsilon)$-coreset of the underlying dataset.

######  Theorem 4.7 (Theorem B.9 in [WY23a]).

For $X\subset[\Delta]^{d}$ of $n$ points, let $\varphi(x)$ be the sensitivity
of $x\in X$, let $q(x)\geq\varphi(x)$, and let
$p(x)\geq\min\left(1,\frac{k^{2}d}{\varepsilon^{2}}\log k\cdot q(x)\right)$.
Then with probability at least $0.99$, independently sampling each point $x\in
X$ with probability $p(x)$ and reweighting each sampled point by
$\frac{1}{p(x)}$ gives a $(1+\varepsilon)$-coreset.

Importantly, the proof of Theorem 4.7 in [WY23a] relies on a simple
Bernstein’s inequality, i.e. Theorem 1.14 and a union bound over a net whose
size is determined by the VC dimension. Thus by applying the following version
of Bernstein’s inequality for $t$-wise independent random variables with
$t=\tilde{\mathcal{O}}\left(\frac{d^{2}k^{2}}{\varepsilon^{2}}\right)\cdot\operatorname*{poly}(\log\log(n\Delta))$,
the same result holds.

######  Theorem 4.8.

[Sko22] Let $X_{1},\ldots,X_{n}$ be $t$-wise independent random variables with
$|X_{i}-\mathbb{E}\left[X_{i}\right]|\leq 1$. Let $X=\sum_{i=1}^{n}X_{i}$ and
$V=\sum_{i=1}^{n}\operatorname*{Var}{X_{i}}$. Then there exists an absolute
constant $C>0$ such that for $d\leq k$ with
$\log\frac{d}{V}<\max\left(\frac{d}{n},2\right)$,

$\mathbf{Pr}\left[|X-\mathbb{E}\left[X\right]|>\gamma\right]\leq
C\cdot\left(\frac{\sqrt{dV}}{\gamma}\right)^{d}.$

In particular, note that we can normalize each independent random variable by
the total cost of $X$ with respect to a fixed set $C$ of $k$ centers, so that
$|X_{i}-\mathbb{E}\left[X_{i}\right]|\leq 1$. The main point of Theorem 4.8 is
that we can now derandomize the hash function $h$ by instead using a $t$-wise
independent hash function on a universe of size $\operatorname*{poly}(n)$,
which can be generated and stored using $\mathcal{O}\left(t\log n\right)$ bits
of space [WC81].

We now justify the correctness of Algorithm 4.

######  Lemma 4.9.

Let $Z$ be the set of weight samples output by $A$ at the end of the stream.
With probability at least $0.96$, $Z$ is a $(1+\varepsilon)$-coreset for $X$
for $k$-median clustering.

###### Proof.

By Lemma 4.4, we have that with probability at least $0.98$,

$\varphi(x)\leq q(x)\leq\mathcal{O}\left(d^{1.5}(\log
k+\log\log\Delta)\right)\cdot\varphi(x).$

Let $Z_{0}$ be the set of weighted points obtained by sampling each $x\in X$
with probability $p(x)$ and reweighting by $\frac{1}{p(x)}$, where
$p(x)=\min\left(1,\frac{k^{2}d}{\varepsilon^{2}}\log k\cdot q(x)\right)$ By
Theorem 4.7, $Z_{0}$ will be a $(1+\varepsilon)$-coreset of $X$ with
probability at least $0.99$. By Lemma 4.6, $Z_{0}$ will have at most $100m$
points with probability at least $0.99$, in which case these points will be
recovered by $A$, by Lemma 4.5. Thus by a union bound, $Z$ is a
$(1+\varepsilon)$-coreset for $X$ for $k$-median clustering with probability
at least $0.96$. ∎

Putting things together, we obtain the following guarantees of Algorithm 4.

######  Theorem 4.10.

Given an accuracy parameter $\varepsilon\in(0,1)$, an integer $k>0$ for the
number of clusters, and a dynamic stream of length $m=\operatorname*{poly}(n)$
defining a set of points $X\subseteq[\Delta]^{d}$ with
$X=\\{x_{1},\ldots,x_{n}\\}$, there exists a one-pass streaming algorithm that
uses
$\tilde{\mathcal{O}}\left(\frac{d^{2.5}k^{2}}{\varepsilon^{2}}\right)\cdot\operatorname*{poly}(\log\log(n\Delta))$
words of memory and outputs a $(1+\varepsilon)$-approximation to $k$-median
clustering.

###### Proof.

Consider Algorithm 4. Correctness follows from Lemma 4.9. To analyze the space
complexity, observe that the first pass sets up the EMD sketch, which by Lemma
4.4, uses a total of $\tilde{\mathcal{O}}\left(kd^{2}\log\log\Delta\right)$
words of space.

Observe that by Lemma 4.6, the sparse recovery data structure in the second
pass uses space
$\tilde{\mathcal{O}}\left(\frac{d^{2.5}k^{2}}{\varepsilon^{2}}\right)\cdot\operatorname*{poly}(\log\log(n\Delta))$
after derandomization. Therefore, the total space used by Algorithm 4 is
$\tilde{\mathcal{O}}\left(\frac{d^{2.5}k^{2}}{\varepsilon^{2}}\right)\cdot\operatorname*{poly}(\log\log(n\Delta))$.
∎

##  5 Wasserstein-$z$ Embedding and $(k,z)$-Clustering

In this section, we present our Wasserstein-$z$ embedding and how to
subsequently utilize it for $(k,z)$-clustering on dynamic streams. We first
present the embedding in Section 5.1. The algorithm and analysis for
$(k,z)$-clustering are presented in Section 5.2.

###  5.1 Wasserstein-$z$ Embedding

We first describe our adaptation of the standard quadtree embedding technique
to achieving a bicriteria approximation for $(k,z)$-clustering. As in the EMD
sketch, we select a shift $s=(s_{1},\ldots,s_{d})\in[\Delta]^{d}$ uniformly at
random. For the sake of presentation, it will be easier to add $s$ to each of
the input points $x_{1},\ldots,x_{n}\in[\Delta]^{d}$ in this case, rather than
shifting the underlying geometric data structure by $s$ as in the case of the
EMD sketch. We then create a tree in the following manner.

Assume without loss of generality that $\Delta$ is a power of two. The root of
the tree represents a grid with side length $2\Delta$ on the non-negative
orthant, so that the origin lies on a corner of the grid. Because
$s,x_{1},\ldots,x_{n}\in[\Delta]^{d}$, then all shifted points are in this
grid with side length $2\Delta$. We thus say that all points belong to the
root node of the tree.

We then partition the grid into $2^{d}$ smaller hypercubes with side length
$\Delta$. For each of these hypercubes that contains a shifted input point, we
create a node representing the hypercube and add the node as a child of the
root node in the tree, using an edge with weight $\sqrt{d}\Delta$. We define
these edges as having height zero, connecting the root node and the children
nodes at height one. Observe that $\sqrt{d}\Delta$ is the maximum distance
between any two points in the parent cube, since the original input points
satisfy $x_{1},\ldots,x_{n}\in[\Delta]^{d}$. Moreover, since there are only
$n$ input points, then the root node has at most $n$ children nodes at height
one. We repeat this procedure until every node representing a hypercube
contains at most a single point, resulting in a tree where all leaves contain
a single point and are all at the same height, which is at most
$\mathcal{O}\left(\log\Delta\right)$.

As it will be useful to index from the leaves of the tree, we will also define
the $t$-th level of the quadtree to be the hypergrid with length $2^{t}$, as
opposed to height $t$ corresponding to the hypergrid with length
$\frac{\Delta}{2^{t-1}}$, recalling that we assume without loss of generality
$\Delta$ is a power of two.

We first embed the input points $X$ into the quadtree. That is, for
$\mu\in\mathbb{R}^{[\Delta]^{d}}$, we define $W_{s,t}\mu$ as the frequency
vector over the hypercubes of the hypergrid ${\mathcal{G}}_{s,t}$ at level $t$
that counts the total mass in each hypercube. To avoid ambiguity, we again say
that each cell of a grid of length $2^{t}$ has closed boundaries on one side
and open boundaries on the other side, e.g., a cell that contains a point
$(x_{1},\ldots,x_{d})$ may contain points $(y_{1},\ldots,y_{d})$ where
$y_{i}<x_{i}+2^{k}$, but does not contain any points $(y_{1},\ldots,y_{d})$
where $y_{i}\geq x_{i}+2^{k}$. As before, we define the mapping $W_{s}\mu$ to
be the concatenation of the vectors

$W_{s}\mu=(W_{s,0}\mu)\circ((2\sqrt{d})^{z}\cdot
W_{s,1}\mu)\circ\ldots\circ((2^{z}\sqrt{d})^{z}\cdot
W_{s,t}\mu)\circ\ldots\circ((2^{\ell}\sqrt{d})^{z}\cdot W_{s,\ell}\mu).$

Now, for a query point $q$, we say that $q$ is _bad at level $i$_ if there
exists a hyperplane of the quadtree decomposition of length $2^{i}$ that has
distance less than $\frac{2^{i}}{d\log\Delta}$ from $q$. Otherwise, we say
that $q$ is good at level $i$. Observe that if $q$ is good at level $i$, but
$q$ and $x\in X$ are first separated at level $i$, then
$\|q-x\|_{2}^{z}>\frac{2^{iz}}{d^{z}\log^{z}\Delta}$ and the incurred
estimated cost for $q$ and $x$ in the quadtree is most $(2^{i}\sqrt{d})^{z}$,
so the distortion will be at most $d^{0.5z}\log^{z}\Delta$.

On the other hand, if $q$ is bad at level $i$, then by definition, it has
distance less than $\frac{2^{i}}{d\log\Delta}$ from a hyperplane of the
quadtree decomposition of length $2^{i}$ and so $q$ and $x$ incur estimated
cost $2^{iz}d$ by the quadtree, then the distortion could be significantly
larger.

Thus for a query point $q$, we define a mapping $\phi$ to subsets $S_{q}$ of
$\mathbb{R}^{d}$ as follows. That is, we set $\phi(q)=S_{q}$ and describe a
well-defined iterative process to construct $S_{q}$.

We initialize $S_{q}=\\{q\\}$, i.e., we first add $q$ to $S_{q}$. If $q$ is
bad at level $i$, then for each hyperplane $H$ of the quadtree of length
$2^{i}$ that is too close to $q$, i.e.,
${\text{dist}}(q,H)\leq\frac{2^{i}}{d\log\Delta}$, we create a copy
$q^{(H)}_{i}$ corresponding to the projection of $q$ onto $H$ (or the closest
point on the projection line on the other side of $H$), so that $q^{(H)}_{i}$
is responsible for serving the points assigned to $q$ that are on the other
side of $H$. We add the point $q^{(H)}_{i}$ to $S_{q}$ for each hyperplane $H$
for which ${\text{dist}}(q,H)\leq\frac{2^{i}}{d^{2}\log\Delta}$, noting that
$q$ could induce multiple points $q^{(H)}_{i}$ to be added to $S_{q}$. We then
proceed top-down across the levels of $i$, repeatedly adding points to $S_{q}$
as necessary. See Figure 4 for an example of this process.

We first show that $\phi(C)$ contains at most $\mathcal{O}\left(k\right)$
centers in expectation.

######  Lemma 5.1.

Given a set $C=\\{C_{1},\ldots,C_{k}\\}$ of $k$ centers, $\phi(C)$ contains at
most $\mathcal{O}\left(k\right)$ centers in expectation over the choice of the
randomized shift parameter $s\in[\Delta]^{d}$.

###### Proof.

Consider a level $i$ and a center $C_{j}$ with $j\in[k]$. Let ${\mathcal{H}}$
be the set of hyperplanes of the hypergrid at level $i$. The probability that
${\text{dist}}(C_{j},{\mathcal{H}})<\frac{2^{i}}{d\log\Delta}$ is at most

$\mathcal{O}\left(\frac{d}{d\log\Delta}\right)=\mathcal{O}\left(\frac{1}{\log\Delta}\right).$

Hence, summing up over all the $\mathcal{O}\left(\log\Delta\right)$ levels,
the expected number of center copies that a center $C_{j}$ can generate is at
most $\mathcal{O}\left(1\right)$. Finally, summing up over $j\in[k]$, then the
expected number of center copies that $C$ can generate, i.e., the number of
additional points in $\phi(C)$ besides $C$, is at most
$\mathcal{O}\left(k\right)$. ∎

By construction, the resulting mapping $\phi(C)$ for a set of $k$ centers
$C=\\{C_{1},\ldots,C_{k}\\}$ will induce an estimated cost by the quadtree
that has small distortion from the actual cost. We now formalize this
statement.

######  Lemma 5.2.

Let $\mu,\nu\in\mathbb{R}^{[\Delta]^{d}}$ be probability measures so that
$\nu$ has support at most $k$ on a set $C$. Then with probability at least
$0.99$, there exists a probability measure $\nu^{\prime}$ with support at most
$\mathcal{O}\left(k\right)$ on the set $\phi(C)$ such that

$\|W_{s}(\mu-\nu^{\prime})\|_{1}\leq\mathcal{O}\left(d^{1+0.5z}\log^{z-1}\Delta\right)\cdot\operatorname*{WASSZZ}(\mu,\nu).$

###### Proof.

Let $C=\\{C_{1},\ldots,C_{k}\\}$ be the support of $\nu$ and let
$x\in[\Delta]^{d}$ be any point in the support of $\mu$. Consider each level
$t$ for which $x$ and $C_{j}$ are separated by the quadtree. For a fixed axis
$j\in[d]$, let ${\text{Proj}}(x)$ denote the projection of a point $x$ onto
the fixed axis and similarly, let
${\text{Proj}}(C)=\cup_{i\in[k]}{\text{Proj}}(C_{i})$. Let
${\text{dist}}({\text{Proj}}(C),{\text{Proj}}(x))=2^{t}\cdot\alpha$ for some
parameter $\alpha\geq 0$.

Observe that if $\alpha\geq 1$, then $C_{j}$ and $x$ will be split at level
$t$ along the direction of axis $j$. This is because $C_{j}$ and $x$ will be
in different cells regardless of the random choice of the hyperplane chosen
along the direction of axis $j$ of the quadtree of length $2^{t}$.

On the other hand, if $\alpha=\mathcal{O}\left(\frac{1}{d\log\Delta}\right)$,
then $C_{j}$ and $x$ will not be split at level $t$ along the direction of
axis $j$. This is because we create a point $C^{\prime}_{j}$ in $\phi(C)$
corresponding to the projection of $C_{j}$ onto the hyperplane along the
direction of axis $j$ of the quadtree of length $2^{t}$, so that
$C^{\prime}_{j}$ and $x$ are not separated by level $j$.

Otherwise, the probability is $\alpha$ that $C_{j}$ and $x$ are split at level
$t$ along the direction of axis $j$, which causes the contribution to
$\|W_{s}(\mu-\nu^{\prime})\|_{1}$ by $\|(2^{t}\sqrt{d})^{z}\cdot
W_{s,t}(\mu-\nu^{\prime})\|_{1}$ to be $(2^{t}\sqrt{d})^{z}$. Hence, with
probability $\alpha$, the distortion between ${\text{dist}}(C,x)$ is at most
$\left(\frac{\sqrt{d}}{\alpha}\right)^{z}$. Therefore, the expected distortion
to the contribution of the cost caused by the direction of axis $j$ is at most
$\mathcal{O}\left(d^{0.5z}\log^{z-1}\Delta\right)$. Summing across all $d$
directions, we have that the expected distortion to the contribution of the
cost is at most $\mathcal{O}\left(d^{1+0.5z}\log^{z-1}\Delta\right)$.

In summary, the distortion to the cost is at most
$\mathcal{O}\left(d^{1+0.5z}\log^{z-1}\Delta\right)$ with probability at least
$0.9999$ by Markov’s inequality. Similarly, by Lemma 5.1 and Markov’s
inequality, it follows that the number of centers in $\nu^{\prime}$ resulting
from $\phi(C)$ is at most $\mathcal{O}\left(k\right)$. Therefore, the desired
claim follows. ∎

We define the mapping $\psi(W_{s},\mu,\nu)$ for probability measures
$\mu,\nu\in\mathbb{R}^{[\Delta]^{d}}$ so that $\nu$ has support at most $k$ on
a set $C$. We define the output of the mapping
$\nu^{\prime}=\psi(W_{s},\mu,\nu)$ as follows. Suppose that mass $m_{i}$ is
moved between $\mu_{i}$ and $\nu_{i}$ in the optimal Wasserstein-$z$
transport. If there is a separating plane $H$ between $\mu_{i}$ and $\nu_{i}$
in $\phi(C)$, then we add mass $m_{i}$ to $\nu^{\prime}_{i}$ at the coordinate
that corresponds to the point that is the projection of $\nu_{i}$ onto $H$.
Otherwise, we add mass $m_{i}$ to $\nu^{\prime}_{i}$ at the coordinate that
corresponds to the point represented by $\nu_{i}$. In other words,
$\nu^{\prime}=\psi(W_{s},\mu,\nu)$ is the set $\phi(C)$ with the appropriate
mass in the optimal Wasserstein-$z$ transport. When the parameters of $W_{s}$
and $\mu$ are clear from context, we abuse notation and use $\psi(\nu)$ to
denote $\psi(W_{s},\mu,\nu)$.

Then from the definition of $\psi$ replacing $\nu^{\prime}$ in Lemma 5.2, we
immediately have: See 1.8

###  5.2 $(k,z)$-Clustering

In this section, we give a two-pass dynamic streaming algorithm for
$(k,z)$-clustering. Recall that for $k$-median, the first pass of our
algorithm in Section 4 was used to set up the EMD embedding, which
subsequently gave a $\mathcal{O}\left(\log k+\log\log\Delta\right)$
approximation to the sensitivity of each point in the second pass. Therefore,
we could perform sensitivity sampling through sparse recovery. We first
describe in Algorithm 5 the subroutine that can be used to estimate the
sensitivity of each point in the stream.

Algorithm 5 Approximation of Sensitivity for $(k,z)$-Clustering

1:Input dataset $X\subset[\Delta]^{d}$

2:Approximate sensitivity for $x$ for any query $x\in[\Delta]^{d}$

3:$m\leftarrow\tilde{\mathcal{O}}\left(kd\log\log\Delta\right)$

4:Initialize a Wasserstein-$z$ sketch, i.e., $m$ instances of $L_{1}$ sketches
on random Wasserstein-$z$ embeddings

5:Update the Wasserstein-$z$ sketch with $X$ $\triangleright$ First pass

6:Use the Wasserstein-$z$ sketch to find a near-optimal set $S$ of $k$ centers

7:Let $\widetilde{\operatorname*{cost}(X,S)}$ be the estimated cost of $S$ on
$X$ $\triangleright$ Lemma 4.3

8:Let ${\mathcal{M}}$ be the net from Lemma 3.2

9:$q(x)\leftarrow 0$

10:for $C\in{\mathcal{M}}$ do

11: $q(x)\leftarrow\max\left(q(x),\frac{2^{z}\cdot
4\gamma\operatorname*{cost}(x,C)}{\operatorname*{cost}(S,C)+\widetilde{\operatorname*{cost}(X,S)}}\right)$
$\triangleright$ Lemma 3.1

12:return $q(x)$

Unfortunately, the Wasserstein-$z$ embedding gives a distortion of
$\mathcal{O}\left(d^{1+0.5z}\log^{z-1}\Delta\right)$, so that sensitivity
sampling in the second pass could incur at least an extra $\log^{z-1}\Delta$
factor. Even for $z=2$, this no longer gives
$\mathcal{O}\left(\log(n\Delta)\right)$ total bits of space if we further use
$\mathcal{O}\left(\log(n\Delta)\right)$ bits of space to represent each
sampled point. Thus, instead of storing an explicit representation of each
point, we store an approximation of its _offset_ from each of the centers.

In particular, for a point $x$ and a set $C^{\prime}$ of
$\mathcal{O}\left(k\right)$ centers, let $c^{\prime}(x)$ be the closest center
of $C^{\prime}$ to $x$. We define the offset $y^{\prime}=x-c^{\prime}(x)$ and
round each coordinate of $y^{\prime}$ to a power of
$\left(1+\operatorname*{poly}\left(\varepsilon,\frac{1}{d},\frac{1}{\log\Delta}\right)\right)$
to form a vector $y$. Then to represent $y$, it suffices to store the identity
of $c^{\prime}(x)$ as well as the _exponents_ of the offsets, which only
requires $\mathcal{O}\left(d\log\frac{\log\Delta}{\varepsilon}\right)$ bits
per sample. Hence for $k$-means clustering, the algorithm still uses
$o(\log^{2}(n\Delta))$ total bits of space. We give the algorithm in full in
Algorithm 6.

Algorithm 6 Two-pass Dynamic Streaming Algorithm for $(k,z)$-Clustering

1:Stream of updates of length $m=\operatorname*{poly}(n)$ to coordinates of
$[\Delta]^{d}$, approximation parameter $\varepsilon\in(0,1)$, number of
clusters $k$, parameter $z\geq 1$

2:$(1+\varepsilon)$-coreset for $(k,z)$-median clustering

3:Initialize a Wasserstein-$z$ sketch $W$ for $(k,z)$-clustering

4:On the first pass, pass each stream update to $W$

5:After the first pass, use $W$ to implicitly define $q(x)$ for all
$x\in[\Delta]^{d}$

6:$s\leftarrow\log^{z-1}\Delta\cdot\operatorname*{poly}\left(k,d,\frac{1}{\varepsilon^{2}},\log\log
n,\log\log\Delta\right)$

7:$p(x)=\min\left(1,\frac{k^{2}d}{\varepsilon^{2}}\log k\cdot q(x)\right)$

8:Use $G$ to find a set $C^{\prime}$ of $\mathcal{O}\left(k\right)$ centers
with cost
$\mathcal{O}\left(2^{z+2}d^{1.5z}\log^{2z}\Delta\right)\cdot{\text{OPT}}$

9:Initialize an $100s$-sparse recovery algorithm $A$ (see Lemma 4.5)

10:Draw hash function $h$ from a family of hash functions where for all
$x\in\mathbb{R}^{d}$, $\mathbf{Pr}\left[h(x)=1\right]=p(x)$ and $h(x)=0$
otherwise

11:for each update in the second pass do