New Frameworks for Offline and
Streaming Coreset Constructions

arXiv:1612.00889v1 [] 2 Dec 2016

Vladimir Braverman

∗

Dan Feldman†

Harry Lang‡

Abstract
Let P be a set (called points), Q be a set (called queries) and a function f : P ×Q →
[0, ∞) (called cost). For an error P
parameter  > 0, a set S ⊆ P with
P a weight function
w : P → [0, ∞) is an ε-coreset if s∈S w(s)f (s, q) approximates p∈P f (p, q) up to a
multiplicative factor of 1 ± ε for every given query q ∈ Q. Coresets are used to solve
fundamental problems in machine learning of streaming and distributed data.
We construct coresets for the k-means clustering of n input points, both in an
arbitrary metric space and d-dimensional Euclidean space. For Euclidean space, we
present the first coreset whose size is simultaneously independent of both d and n. In
particular, this is the first coreset of size o(n) for a stream of n sparse points in a
d ≥ n dimensional space (e.g. adjacency matrices of graphs). We also provide the first
generalizations of such coresets for handling outliers. For arbitrary metric spaces, we
improve the dependence on k to k log k and present a matching lower bound.
For M -estimator clustering (special cases include the well-known k-median and
k-means clustering), we introduce a new technique for converting an offline coreset
construction to the streaming setting. Our method yields streaming coreset algorithms
requiring the storage of O(S + k log n) points, where S is the size of the offline coreset.
In comparison, the previous state-of-the-art was the merge-and-reduce technique that
required O(S log2a+1 n) points, where a is the exponent in the offline construction’s
dependence on −1 . For example, combining our offline and streaming results, we
produce a streaming metric k-means coreset algorithm using O(−2 k log k log n) points
of storage. The previous state-of-the-art required O(−4 k log k log6 n) points.

∗

Department of Computer Science, Johns Hopkins University. This material is based upon work supported
in part by the National Science Foundation under Grant No. 1447639, by the Google Faculty Award and
by DARPA grant N660001-1-2-4014. Its contents are solely the responsibility of the authors and do not
represent the official view of DARPA or the Department of Defense.
†
Department of Computer Science, University of Haifa
‡
Department of Mathematics, Johns Hopkins University. This material is based upon work supported by
the Franco-American Fulbright Commission.

1

Introduction

In the algorithmic field of computer science, we usually have an optimization problem at
hand and a state-of-the-art or a straight-forward exhaustive search algorithm that solves it.
The challenge is then to suggest a new algorithm with a better running time, storage or
other feature. A different and less traditional approach is to use data reduction, which is a
compression of the input data in some sense, and to run the (possibly inefficient) existing
algorithm on the compressed data. In this case, the problem of solving the problem at hand
reduced to the computing a problem-dependent compression such that:
1. an existing algorithm that solves the optimization problem on the original (complete)
data, will yield a good approximate solution for the original data when applied on the
compressed data.
2. The time and space needed for constructing the compression and running the optimization algorithm on the coreset will be better than simply solving the problem on
the complete data.
There are many approaches for obtaining such a provable data reduction for different
problems and from different fields, such as using uniform sampling, random projections (i.e.,
the Johnson-Lindenstrauss Lemma), compressed sensing, sketches or PCA.
In this paper we focus on a specific type of a reduced data set, called coreset (or coreset) that was originated in computational geometry, but now applied in other fields such as
computer vision and machine learning. Our paper is organized into basic sections: results for
maintaining coresets over data streams (Section 3) and results for offline coresets (Sections 46). We briefly introduce both of these topics the remainder of this section. Many of our
results, along with comparison to prior works, are summarized in Table 1 in Section 2.
In the Appendix A (Section 7) we summarize the merge-and-reduce technique that is used
in previous approaches [Che09a, HPM04, HPK07, AMR+ 12, FL11]. In the Appendix B (Section 8) we provide an alternative framework that generalizes our main result (Theorem 3.1),
applying to a wide-array of constructions although giving a weaker bound.

1.1

Streaming Results

In the streaming model of computation, the input arrives sequentially. This differs from the
standard model where the algorithm is given free access to the entire input. Given a memory
that is linear in the size of the input, these models are evidently equivalent; therefore the
goal of a streaming algorithm is to perform the computation using a sublinear amount of
memory.
Our stream consists of n elements p1 , . . . , pn . In the streaming model (or more specifically
the insertion-only streaming model, since points that arrive will never be deleted), we attempt
to compute our solution using o(n) memory. Sometimes the algorithm will be allowed to
pass over the stream multiple times, resulting in another parameter called the number of
passes. All of our algorithms use polylog(n) memory and require only a single pass.
1

Prior to the current work, the merge-and-reduce technique due to Har-Peled and Mazumdar [HPM04] and Bentley and Sax [BS80] was used to maintain a coreset on an insertion-only
stream. For a summary of this technique, see Section 7 in the Appendix. In this paper we
introduce an alternative technique that reduces the multiplicative overhead from log2a+1 n
to log n (here, a is the offline construction’s dependence on 1/). While our method is not
as general as merge-and-reduce (it requires that the function in question satisfies more than
just the “merge” and “reduce” properties, defined in Section 7), it is general enough to apply
to all M -estimators. For the special case of our coreset offline construction for M -estimators
(introduced in Section 6), we use a more tailored method that causes this to be log n additive
overhead. Therefore our streaming space complexity matches our offline space complexity,
both of which improve upon the state-of-the-art.
The offline coreset construction of [FL11] has the following structure: first, a bicriterion approximation is computed. Second, points are sampled according a distribution that
depends only on the distances between points of the input and their assigned bicriterion
centers. This suggests a two-pass streaming algorithm (which we later combine into a single
pass): in the first pass, construct a bicriterion using an algorithm such as [BMO+ 11]. In
the second pass, sample according to the bicriterion found in the first pass. This provides a
two-pass algorithm for a coreset using O(−2 k log k log n)-space. Our contribution is showing
how these two passes can be combined into a single-pass. Using the algorithm of [BMO+ 11]
to output O(k log n) centers at any time, we show that this is sufficient to carry out the
sampling (originally in the second pass) in parallel without re-reading the stream. Our main
lemma (Lemma 3.7) shows that the bicriterion, rather than just providing “central” points to
concentrate the sampling, actually can be thought of as a proof of the importance of points
for the coreset (technically, a bound on the “sensitivity” that we define at the beginning of
Section 3). Moreover, the importance of points is non-increasing as the stream progresses,
so we can maintain a sample in the streaming setting without using any additional space.

1.2

Offline Results

The name coreset was suggested by Agarwal, Har-Peled, and Varadarajan in [AHPV04] as a
small subset S of points for a given input set P , such that any shape from a given family that
covers S will also cover P , after expanding the shape by a factor of (1 + ε). In particular,
the smallest shape that covers S will be a good approximation for the smallest shape that
covers P . For approximating different cost functions, e.g. the sum of distances to a given
shape, we expect that the total weight of the sample will be similar to the number n of input
points. Hence, in their seminal work [HPM04], Har-Peled and Mazumdar used multiplicative
weights for each point in S, such that the weighted sum of distances from S to a given shape
from the family will approximate its sum of distances from the original data. In [HPM04]
each shape in the family was actually a set of k points, and the application was the classic
k-means problem.
In this paper, we are given an input set P (called points), a family (set) Q of items,
called queries and a function f : P → [0, ∞) that is called a cost function. A coreset is then
a subset S of P , that is associated with a non-negative weight function u : S → [0, ∞) such
2

P
that, for every given
query
q
∈
Q,
the
sum
of
original
costs
p∈P f (p, q) is approximated by
P
the weighted sum p∈S u(p)f (p, q) of costs in S up to a multiplicative factor, i.e.,
X
X
X
(1 − ε)
f (p, q) ≤
u(p)f (p, q) ≤ (1 + ε)
f (p, q)
p∈P

p∈S

p∈P

While our framework is general we demonstrate it on the k-means problem and its variant.
There are at least three reasons for this: (i) This is a fundamental problem in both computer
science and machine learning, (ii) This is probably the most common clustering technique
that used in practice, (iii) Many other clustering and non-clustering problems can be reduced
to k-means; e.g. Mixture of Gaussians, Bregman Clustering, or DP-means [FFK11, LBK15,
BLK15]. In this context we suggest offline coreset constructions for k-clustering queries, that
can be constructed in a streaming fashion using our streaming approach, and are:
• of size linear in k (for d > log k) and arbitrary metric space of dimension d. Current
coresets that are subset of the input have size at least cubic in k [LS10]. This is by
reducing the total sensitivity to O(1) without introducing negative weights that might
be conditioned on the queries as in [FL11].
• of size independent of d for the Euclidean case of k-means (squared distances). This
is particular useful for sparse input set of points where d ≥ n, such as in adjacency
matrices of graphs, document-term, or image-object matrices. Recent coreset for sparse
k-means of [BF15] is of size exponential in 1/ε and thus turn to O(n) when used when
the merge-and-reduce tree (where ε is replaced by O(ε/ log(n)). The result of [FSS13]
for k-means is exponential in k/ε and fails with constant probability, so also cannot
be used with streaming. Another result of [FSS13] suggests a coreset type set for kmeans of size O(k/ε) but which is based on projections that loss the sparsity of the
data. Similar sparsity loss occurs with other projection-type compression methods e.g.
in [CEM+ 15]. Nevertheless, we use the technique in [FSS13] to bound the dimension
of the k-means problem by O(k/ε).
• of size independent of d for the Euclidean case, and non-squared distances, using weak
coresets. These coresets can be used to approximates the optimal solution, but not
every set of k centers. Unlike the weak coresets in [FL11, FMS07], we can use any
existing heuristic on these coresets, as explained in Section 6.1.
• Robust to outliers. This is since the general pseudo-metric definition we used (inspired
by [FS12], support m-estimators which is a tool for handling outliers [Tyl87]. Unlike
in [FS12] our coresets are linear (and not exponential) in k, and also independent of n
(and not logarithmic in n).

2

Related Work

The following table summarizes previous work along with our current results. By far,
the most widely-studied problems in this class have been the k-median and k-means func3

tions. In general, the extension to arbitrary M -estimators is non-trivial; the first such result
was [FS12]. Our approach naturally lends itself to this extension. M -estimators are highly
important for noisy data or data with outliers. As one example, Huber’s estimator is widely
used in the statistics community [HHR11, Hub81]. It was written that “this estimator is so
satisfactory that it has been recommended for almost all situations” [Zha11]. Our results
work not only for Huber’s estimator but for all M -estimators, such as the Cauchy and Tukey
biweight functions which are also well-used functions.
Note that in the below table, Õ notation is used to write in terms of d, , k, and log n
(therefore hiding factors of log log n but not log n).
Problem
Offline Size
Streaming Size
−d
Euclidean k-means
O(k log n)
O(k−d logd+2 n)
Euclidean k-means
O(k 3 −(d+1) )
O(k 3 −(d+1) logd+2 n)
Euclidean k-means
O(dk 2 −2 log n)
O(dk 2 −2 log8 n)
Euclidean k-means
O(dk log k−4 )
O(dk log k−4 log5 n)
Õ((d/)O(d) k log n)
Õ((d/)O(d) k logO(d) n)
Euclidean k-means
−2
−2
Euclidean k-means
O( k log k min(k/, d)) O( k log k min( k , d) + k log n)
Metric k-means
O(−2 k 2 log2 n)
O(−2 k 2 log8 n)
Metric k-means
O(−4 k log k log n)
O(−4 k log k log6 n)
−2
Metric k-means
O( k log k log n)
O(−2 k log k log n)
Õ(dk 2 −2 log n)
Euclidean k-median
O(dk 2 −2 log8 n)
Euclidean k-median
O(k−d log n)
O(k−d logd+2 n)
Euclidean k-median
O(k 2 −d )
O(k 2 −(d) logd+1 n)
Euclidean k-median
O(d−2 k log k)
O(d−2 k log k log3 n)
Euclidean k-median
O(d−2 k log k)
O(d−2 k log k + k log n)
2
Metric k-median
O(k 2 −2 log n)
O(k 2 −2 log8 n)
Metric k-median
O(−2 k log k log n)
O(−2 k log k log4 n)
Metric k-median
O(−2 k log k log n)
O(−2 k log k log n)
2
Euclidean M -estimator
O(−2 k O(k) d2 log n)
O(−2 k O(k) d2 log5 n)
Euclidean M -estimator
O(d−2 k log k)
O(d−2 k log k + k log n)
Metric M -estimator
O(−2 k O(k) log4 n)
O(−2 k O(k) log7 n)
−2
Metric M -estimator
O( k log k log n)
O(−2 k log k log n)

Paper
[HPM04]
[HPK07]
[Che09a]
[FL11]
[AMR+ 12]
**
[Che09b]
[FL11]
**
[Che09a]
[HPM04]
[HPK07]
[FL11]
**
[Che09a]
[FL11]
**
[FS12]
**
[FS12]
**

Table 1: Summary of Related Work

Framework. A generic framework for coreset construction was suggested in [FL11]. The
main technique is a reduction from coreset to ε-approximations, that can be computed using
non-uniform sampling. The distribution of the sampling is based on the importance of each
point (in some well defined sense), and the size of the coreset depends on the sum of these
importance levels. This term of importance appeared in the literature as leverage score (in
the context of low-rank approximation, see [PKB14] and references therein), or, for the case
of k-clustering, sensitivity [LS10]. The proof of many previous coreset constructions were
4

significantly simplified by using this framework, and maybe more importantly, the size of
these coresets was sometimes significantly reduced; see [FL11] for references. Many of these
coresets size can be further improve by our improved framework, as explained below.
The size of the coreset in [FL11] depends quadratically on the sum of sensitivities, called
total sensitivity [LS10]. In this paper, we reduce the size of the coreset that are constructed
by this framework to be only near-linear (t log t) in the total sensitivity t. In addition, we
generalize and significantly simplify the notation and results from this framework.
k-means.
In the k-means problem we wish to compute a set k of centers (points) in
some metric space, such that the sum of squared distances to the input points is minimized,
where each input point is assigned to its nearest center. The corresponding coreset is a
positively weighted subset of points that approximates this cost to every given set of k
centers. First deterministic coresets of size exponential in d were first suggested by HarPeled and Mazumdar in [HPM04]. The first coreset construction of size polynomial in d was
suggested by Ke-Chen in [Che09a] using several sets of uniform sampling.
The state-of-the-art is the result of Schulman and Langberg [LS10] who suggested a
coreset of size O(d2 k 3 /ε2 ) for k-means in the Euclidean case based on non-uniform sampling.
The distribution is similar to the distribution in [FMS07] over the input points, however
in [FMS07] the goal was to have weaker version coresets that can be used to solve the
optimal solution, but are of size independent of d.
Some kind of coreset for k-means of size near-linear in k was suggested in [FL11]. However, unlike the definition of this paper, the multiplicative weights of some of the points in
this coreset were (i) negative, and (ii) depends on the query, i.e., instead of a weight w(p) > 0
for an input point p, as in this paper, the weight is w(p, C) ∈ R where C is the set of queries.
While exhaustive search was suggested to compute a PTAS for the coreset, it is not clear how
to compute existing algorithms or heuristics (what is actually done in practice) on such a
coreset. On the contrary, generalizing an existing approximation algorithm for the k-means
problem to handle positively weights is easy, and public implementations are not hard to
find (e.g. in Matlab and Python).
k-mean for handling outliers.
Coresets for k-means and its variants that handle
outliers via m-estimators were suggested in [FS12], and also inspired our paper. The size
of these coresets is exponential in k and also depend on logn. For comparison, we suggest
similar coreset of size near-linear in k, and independent of n. PTAS for handling exactly
m-outliers was suggested in [Che08] but with no coreset or streaming version.
Streaming.
The metric results of [Che09a, FL11] and Euclidean results of [Che09a,
HPM04, HPK07, FL11] that rely on merge-and-reduce have already been mentioned. A
summary of these results appears in the tables below. For the specific case of Euclidean
space, a more diverse set of stronger results is known. In particular, coreset constructions
are known that do not begin with a bicriterion solution, and whose streaming variant does
not rely on merge-and-reduce [AMR+ 12]. With the additional assumption in Euclidean
5

space that the points lie on a discrete grid {1, . . . , ∆}d , alternative techniques are known for
k-means and other problems, even when the stream allows the deletion of points [FS05].

3

Streaming Algorithm

We present a streaming algorithm for constructing a coreset for metric k-means clustering
that requires the storage of O(−2 k log n) points. The previous state-of-the-art [FL11] required the storage of O(−4 k log k log6 n) points. In this section we assume the correctness
of our offline algorithm, which is proven in Section 6.
More generally, our technique works for M -estimators, a general class of clustering objectives that includes the well-known k-median and k-means functions as special cases. Other
special cases include the Cauchy functions, the Tukey functions, and the Lp norms. Our
method combines a streaming bicriterion algorithm [BMO+ 11] and a batch coreset construction [FL11] to create a streaming coreset algorithm. The space requirements are combined
addivitely, therefore ensuring no overhead.
The streaming algorithm of [BMO+ 11] provides a bicriterion solution using O(k log n)
space. Our new offline construction of Section 6 requires O(−2 k log k log n) space. Therefore
our main theorem yields a streaming algorithm that combines these spaces additively, therefore requiring O(−2 k log k log n) space while maintaining a coreset for k-means clustering.
The previous state-of-the-art framework that works for the metric variant (other methods
are known to improve upon this for the special case of Euclidean space) was the mergeand-reduce technique [BS80] that yields a streaming algorithm requiring O(−4 k log k log6 n)
space, incurring an overhead of Θ(log5 n) over the offline coreset size. In comparison, our
framework incurs no overhead. The additional improvement in our space is due the improved
offline construction given in Sections 5-7.
We now state our main theorem. The result is stated in full generality: the k-median
clustering in a ρ-metric space (see Definition 6.1). Note that metric k-means clustering
corresponds to setting ρ = 2. Also, the probability of success 1 − δ typically has one of two
meanings: that the construction succeeds at the end of the stream (a weaker result), or that
the construction succeeds at every intermediate point of the stream (a stronger result). Our
theorem gives the stronger result, maintaining a valid coreset at every point of the stream.
Theorem 3.1 (Main Theorem). There exists an insertion-only streaming algorithm that
maintains a (k, )-coreset for k-median clustering in a ρ-metric space, requires the storage
of O(ρ2 −2 k log(ρk) log(n) log(1/δ)) points, has poly(k, log n, ρ, , log(1/δ)) worst-case update
time, and succeeds at every point of the stream with probability 1 − δ.
Our method can be applied to the coreset constructions of [HPM04, HPK07, Che09a,
FL11] with a multiplicative overhead of O(log n). Our second theorem is a more generally
applicable technique; it applies to all constructions that first compute a bicriterion solution
and then sample points according to the bicriterion solution. The constructions of [HPM04,
HPK07, Che09a, FL11] follow this outline, and we are unaware of any constructions which do

6

not. The theorem yields immediate corollaries as well as reducing certain streaming coreset
problems to that of constructing an offline coreset.
Theorem 3.2. Given an offline algorithm that constructs a (k, )-coreset consisting of S =
S(n, k, , δ) points with probability 1 − δ by sampling points based on a bicriterion solution,
there exists a streaming algorithm requiring the storage of O(S log n) points that maintains
a (k, )-coreset on an insertion-only stream.
Proof Sketch. The known offline coreset constructions start with a bicriterion solution of
O(k) points. We modify the algorithm of [BMO+ 11] to output O(k log n) centers; this is
trivial since the final step of the algorithm of [BMO+ 11] is to take the O(k log n) centers
stored in memory and reduce them to exactly k centers to provide a solution. Our first
modification to the original algorithm is thus to simply remove this final step, but we must
also keep a datastructure storing log(1/) intermediate states of these O(k log n) centers.
See Section 8 for a precise description of our modification and the sampling method, applied
to the construction of [FL11] as an example (but equally applicable to [HPM04, HPK07,
Che09a]). As the high-level idea, since the bicriterion given to the offline construction
consists of O(k log n) centers instead of exactly k, the number of additional points taken
in the coreset increases by a factor of O(log n).
Two important corollaries include:
1. Using the result of [FL11], we obtain a streaming algorithm that maintains a (k, )coreset with negative weights for metric k-median requiring the storage of O(−2 k log n)
points.
2. Given a O(k · poly(, log n, log(1/δ))) point (k, )-coreset, we would obtain a streaming
algorithm that maintains a (k, )-coreset (with only positive weights) for metric kmedian requiring the storage of O(k · poly(, log n, log(1/δ))) points. This differs from
Theorem 3.1 in that the dependence on k is linear instead of O(k log k).

3.1

Definitions

We begin by defining a ρ-metric space, which is defined in full as Definition 6.1. Briefly, let
X be a set. If D : X × X → [0, ∞) is a symmetric function such that for every x, z ∈ X we
have that D(x, z) ≤ ρ(D(x, y) + D(y, z)) for every y ∈ X, when we call (X, D) a ρ-metric
space. Note that this is a weakening of the triangle inequality, and at ρ = 1 we recover the
definition of a metric space. All M -estimators can be re-cast for a certain constant value of
ρ, and k-means is obtained with ρ = 2. This generality is therefore useful and working in
this language allows us to naturally generalize our results to any M -estimator.
The k-median problem is, given an input set P and an integer k ≥ 1, to find a set C of
k points that minimizes:
X
min D(p, c)
p∈P

c∈C

7

We use OPTk (P ) to denote this minimal value. As this is NP-Hard to compute, we settle
for an approximation. The notion of a bicriterion approximation is well-known; we state a
definition that suits our needs while also fitting into the definition of previous works.
Definition 3.3 ((α, β)-approximation). An (α, β)-approximation forP
the k-median clustering
of a multiset P is a map π : P → B for some set B such that
p∈P w(p)D(p, π(p)) ≤
αOPTk (P ) and |B| ≤ βk.
We now define a coreset:
Definition 3.4 ((k, )-coreset). A (k, )-coreset of a multiset P is a weighted
Pset (S, v) with
k
non-negative weight functionPv such that for every Z ∈ X we have (1 − ) p∈P D(p, Z) ≤
P
s∈S v(s)D(s, Z) ≤ (1 + )
p∈P D(p, Z).
Coresets with arbitrary weight functions (i.e. with negative weights allowed) have been
considered [FL11, etc]. However, computing approximate solutions on these coresets in
polynomial-time remains a challenge, so we restrict our definition to non-negative weight
functions. This ensures that an approximate solution can be quickly produced. This implies
a PTAS for Euclidean space and a polynomial-time 2γ(1 + )-approximation for general
metric spaces (where γ is the best polynomial-time approximation factor for the problem
in the batch setting). This factor of 2γ(1 + ) is well-known in the literature, see [COP03,
BMO+ 11, GMM+ 03] for details.

3.2

Constant-Approximation Algorithm

Let Pi denote the prefix of the stream {p1 , . . . , pi }. The entire stream is then Pn . Consider
the moment when the first i points have arrived, meaning that the prefix Pi is the current set
of arrived points. The algorithm A of [BMO+ 11] provides an (O(1), O(log n))-approximation
of Pi in the following sense. Define f0 : ∅ → ∅ as the null map, and define Bi = image(fi ).
Upon receiving point pi , algorithm A defines a map fi : Bi−1 ∪ {pi } → Bi . We define
πi : Pi → Bi by πi (pj ) = fi (fi−1 (. . . (fj (pj )) . . .)) for each 1 ≤ j ≤ i. These mappings have
an essential gaurantee stated in the following lemma.
Theorem 3.5 ([BMO+ 11]). For every 1 ≤ i ≤ n, after receiving Pi , Algorithm A(k, n, δ)
definesP
a function fi such that with probability 1 − δ, using the above definition of πi , the
bound
p∈Pi D(p, πi (p)) ≤ αOPTk (Pi ) holds. The algorithm deterministically requires the
storage of O(k(log n + log(1/δ))) points.

3.3

Offline Coreset Construction

We briefly describe the offline coreset construction. The proof of correctness can be found
in Sections 5 and 6. It is this construction that we will maintain in the streaming setting.
The sensitivity of a point p ∈ P is defined as:
D(p, Z)
q∈P D(q, Z)

s(p) = max P
Z∈X k

8

0
Notice that 0 ≤
P s(p) ≤ 1. We give an upper bound s (p) ∈0 [s(p), 1]. Define
Pthe total
0
sensitivity t = p∈P s(p). Likewise, we give an upper bound t ≥ t where t = p∈P s0 (p)
and will show that t = O(k). The sampling probability distribution at point p is set to
s0 (p)/t0 . We take an i.i.d. sample from P of size m for any m ≥ ct0 −2 (log n log t0 + log(1/δ))
where c is a constant.
Let R be the union of these m i.i.d. samples, and then define a weight function v : R →
[0, ∞) where v(r) = (|R|s0 (r))−1 . It is proven as one of our main theorems (Theorem 6.6)
that the weighted set (R, v) is a (k, )-coreset for P .

3.4

Bounding the Sensitivity

Consider the prefix Pi which is the input after the first i points have arrived. Using Algorithm
A we obtain an (α, β)-approximation πi where α = O(1) and β = O(log n). Recall that Bi
is the image of this approximation, i.e. Bi = image(πi (Pi )).
Running an offline (γ, λ)-approximation algorithm on Bi , we obtain a multiset Ci of at
most λk distinct points. Let p0 denote the element of Ci nearest to πi (p) (this is the element
of Ci P
that p gets mapped to when we pass from Pi → Bi → Ci ). The following lemma implies
that p∈P w(p)D(p, p0 ) ≤ ᾱOPT(P ) where ᾱ = ρα + 2ρ2 γ(α + 1). This is an observation
used widely in the literature [COP03, BMO+ 11], but we include a proof for completeness.
Lemma 3.6. Let B be a (α, β)-approximation of P , and let C be a (γ, λ)-approximation of
B. Then C is a (ρα + 2ρ2 γ(α + 1), λ)-approximation of A.
Proof. Let π : P → B be the (α, β)-approximation of P and let t : B → C be the (γ, λ)approximation
P of B. In the following, all sums
P will be taken over all p ∈ P . The hypotheses
∗
state that
D(p, π(p)) ≤ αOPT
P Let P be∗ an
P(P ) and∗ D(π(p), t(π(p))) 1≤ γ OPT(B).
D(π(p), P ) ≤
optimal
clustering of P , that is
D(p, P ) = OPT(P ). Then 2 OPT(B) ≤
P
ρ (D(π(p), p) + D(p, P ∗ )) ≤ ρ(α + 1)OPT(P ). The factor of 12 comes from the fact that
OPT(B) is defined using centers restricted to B (see [GMM+ 03] for details). We now write
P
P
D(p, t(π(p))) ≤ ρ (D(p, π(p))+D(π(p), t(π(p)))) ≤ (ρα+2ρ2 γ(α+1))OPT(P ) as desired.
We now prove the following lemma which gives us our sampling probability s0 (p). Recall
that for the construction to succeed, the sampling probability s0 (p) must be at least the
sensitivity s(p) (defined in the previous subsection). Since we focus on a single iteration, we
drop subscripts and write C = Ci and P = Pi . Let p 7→ p0 be an (ᾱ, λ)-approximation of P .
Define P (p) = {q ∈ P : q 0 = p0 } to be the cluster containing p.
Lemma 3.7. Let the map p 7→ p0 define an (ᾱ, λ)-approximation for the k-median clustering
of P . For every point p ∈ P :
ρᾱD(p, p0 )
ρ2 (ᾱ + 1)
P
s(p) ≤
+
0
|P (p)|
q∈P D(q, q )

9

Proof. For an arbitrary Z ∈ X k we need to provide a uniform bound for
ρD(p, p0 )
ρD(p0 , Z)
D(p, Z)
P
P
≤
+
q∈P D(q, Z)
q∈P D(q, Z)
q∈P D(q, Z)

P

ρD(p0 , Z)
ᾱρD(p, p0 )
P
+
≤P
0
q∈P D(q, q )
q∈P D(q, Z)

(1)

P
P
where the second inequality holds because q∈P D(q, q 0 ) ≤ ᾱOPT(P ) ≤ q∈P D(q, Z). To
bound the last term, recall that q 0 = p0 for all q ∈ P (p) so:
X
X
D(p0 , Z)|P (p)| =
D(p0 , Z) =
D(q 0 , Z)
q∈P (p)

≤ρ

q∈P (p)

X

(D(q 0 , q) + D(q, Z))

q∈P (p)

≤ρ

X

D(q 0 , q) + ρ

q∈P

≤ ρᾱ

X

D(q, Z)

q∈P (p)

X

D(q, Z) + ρ

q∈P

X

D(q, Z)

q∈P (p)

≤ ρ(ᾱ + 1)

X

D(q, Z)

q∈P

Dividing by |P (p)|

P

q∈P

D(q, Z) gives
ρ(ᾱ + 1)
D(p0 , Z)
≤
|P (p)|
q∈P D(q, Z)

P

Substituting this in (16) yields the desired result.
We therefore define our upper bound s0 (p) as in the
An immediate but extremely
P lemma.
0
0
important consequence of Lemma 3.7 is that t = p∈P s (p) = ρᾱ + ρ2 (ᾱ + 1)k ≤ 3ρ2 ᾱk.
This can be seen by directly summing the formula given in the lemma.

3.5

Streaming Algorithm

We now state Algorithm 1, which we then prove maintains a coreset. To use Lemma 3.7
0
to
P determine0 s (p), we will compute the cluster sizes |P (p)| and estimate the clustering cost
q∈P D(q, q ). We must bound the clustering cost from below because we need an upperbound of s(p). On Line 9, L is an estimate of the cost of clustering P to the centers C. On
Line 4, c is the absolute constant used in Theorem 6.6.
Algorithm 1 outputs (R, v) such that point p is sampled with probability xs0 (p) where x
is defined on Line 4. For each p that has arrived, the value of s0 (p) is non-increasing (notice
that it is defined as the minimum of itself and a new value on Line 13), so it is possible to
10

Algorithm 1: Input: stream of n points in a ρ-metric space,  > 0, k ∈ N, maximum
failure probability δ > 0
1 Initilization:
2 R ← ∅
0
2
3 t ← ρᾱ + ρ (ᾱ + 1)k
−2
0
4 x ← 2c (log n log t + log(1/δ))
5 Initialize A(k, n, δ)
6 Update Process: after receiving pi
7 (Bi , fi ) ← A.update(pi )
8 C ← an (γ, λ)-approximation of Bi
P
−1
−1
9 L ← D(pi , C) + ᾱ (1 + )
r∈R v(r)D(r, C)
10 for r ∈ R do
11
πi (s) ← fi (πi−1 (r))
0)
2
12
z(r) ← ρᾱD(p,p
+ ρ|P(ᾱ+1)
L
(p)|
13
s0 (r) ← min(s0 (r), z(r))
14 for r ∈ R do
15
if u(r) > xs0 (r) then
16
Delete r from R
17 u(pi ) ← uniform random from [0, 1)
0
18 if u(pi ) ≤ xs (pi ) then
19
Add pi to R
20
πi (pi ) ← fi (pi )
21 Query Process:
22 for each r ∈ R do
23
v(r) ← (|R|s0 (r))−1

11

maintain this in the streaming setting since once the deletion condition on Line 15 becomes
satisfied, it remains satisfied forever.
We now proceed with the proof of Theorem 1. Since the probability of storing point p is
xs0 (p), the expected space of Algorithm 1 is xt0 . By Lemma 3.7 that implies t0 ≤ 3ρ2 ᾱk, we
then bound the expected space as 2c−2 (log n log t + log(1/δ))(3ρ2 ᾱk). Simplifying notation
by defining an absolute constant c̃ (a function of c and ᾱ), we write this expected space as
c̃ρ2 −2 k log(ρk) log n log(1/δ). By a Chernoff bound, the high-probability gaurantee follows
by replacing c̃ with 2c̃.
Proof of Theorem 1. The proof of Theorem 1 can be divided into the following pieces:
1. P
For correctness (to satisfy the bound given in Lemma 3.7, we must show that L ≤
0
p∈Pi D(p, p ). For space, it is important that L is not too small. In particular, the
space grows as 1/L. We show that L is a -approximation of the true cost.
2. The value of |P (p)| can be computed exactly for every p. This is needed on Line 12.
3. The construction of Algorithm 1 that samples p with probability xs0 (p) can be processed
to be identical to the offline construction of Subsection 3.3 that takes an i.i.d. sample
of size xt0 from the distribution where p is given sampling probability s0 (p)/t0 .
1: To lower bound the clustering cost, inductively assume that we have a (k, )-coreset
Si−1 of Pi−1 . Note that pi is a (k, )-coreset of itself (in fact a (k, 0)-coreset of itself),
so Si−1 ∪ {pi } is a (k, )-coreset of Pi . Let L be the cost of clustering Si−1 ∪ {pi } to C.
Therefore the cost of clustering Pi to C is in the interval [(1 − )L, (1 + )L]. Recall that the
upper bound on s(p) from Lemma 3.7 is:
ρᾱD(p, p0 )
ρ2 (ᾱ + 1)
P
+
0
|P (p)|
q∈P D(q, q )
P
0
0
0
By using L in place of the true cost
p∈Pi D(p, p ) for defining s (p), the value of t =

P
1+
0
2
4
p∈Pi s (p) increases to at most ρᾱ 1− +ρ (ᾱ+1)λk = O(ρ k). Here there is no dependence
1+
is bounded by an absolute constant.
on  since we assume  ≤ 1/2, so 1−
2: Computing |P (p)| is straightforward. Define w(b) = |{p ∈ P
P : π(p) = b}| and then let
h : B → C be the (γ, λ)-approximate clustering. Then |P (p)| = b∈h−1 (p0 ) w(b).
3: In Algorithm 1, we sample point p with probability s0 (p) to maintain a set M of
non-deterministic size. We now argue that this can be converted to the desired coreset,
where an i.i.d. sample of size m is taken from the distribution s0 /t0 . First, by a Chernoff bound, |R| ≥ E[|R|]/2 with probability 1 − exp(−E[|R|]/8). Since E[|R|] = xt0 =
(2c−2 log(ρk) log n log(1/δ)) · (ρᾱ + ρ2 (ᾱ + 1)k) = Ω(log(n)), we have that |R| ≥ E[|R|]/2
with probability 1 − O(1/n). Then by the union bound, this inequality holds true at each of
the n iterations of receiving a point throughout the entire stream. Recall that for the offline
coreset construction outlined in Subsection 3.3 to hold, we need an i.i.d. sample of at least
m = ct0 −2 (log n log t0 + log(1/δ)). By Lemma 3.7 t0 = ρᾱ + ρ2 (ᾱ + 1)k, and so by plugging
in values we see that E[|R|] = xt0 ≥ 2m. Having that |R| ≥ m with probability 1 − O(1/n),
12

it is well-known (see [DDH+ 07] for example) that this can be converted to the required i.i.d.
sample.

4
4.1

Preliminaries for Offline Coreset Construction
Query space

In our framework, as in [FL11], we are given a finite input set P of items that are called
points, a (usually infinite) set Q of items that are called queries, and a cost function f that
maps each pair of a point in P and a query in Q to a non-negative number f (p, q). The cost
of the set P to this query is the sum over all the costs,
X
f¯(P, q) :=
f (p, q).
p∈P

More generally, each input point might be given a positive multiplicative weight w(p) > 0,
and the overall cost of each point is then reweighed, so that
X
w(p)f (p, q).
f¯(P, w, q) =
p∈P

The tuple (P, w, f, Q) is thus define our input problem and we call it a query space. In the
case that the points are unweighted we can simply define w(p) = 1 for every p ∈ P .
However, for the following sections, it might help to scale the weights so that their sum
is 1. In this case, we can think of the weights as a given distribution over the input points,
and the cost f¯(P, w, q) is the expected value of f (p, q) for a point that is sampled at random
from P . For the unweighted case we thus have w(p) = 1/n for each point, and
1X
f (p, q)
f¯(P, w, q) =
n p∈P
is the average cost per point. The cost function f is usually also scaled as will be explained
later, to have values f (p, q) between 0 and 1, or −1 to 1.
Example 1: Consider the following problem where the input is a set P of n points in
d
R . Given a ball B = B(c, r) of radius r that is centered in c ∈ Rd , we wish to compute the
fraction of input points that are covered by this ball, |P|P∩B|
. More generally, the query is a
|
set of k balls, and we wish to compute the fraction of points in P that are covered by the
union of these balls. In this case, each input point p ∈ P has a weight w(p) = 1/n, the set
Q of queries is the union over every k balls in Rd ,

Q = B(c1 , r1 ) ∪ · · · ∪ B(ck , rk ) | ∀i ∈ [k] : ci ∈ Rd , ri ≥ 0 ,
(2)
and the cost f (p, q) for a query q = {B1 ∪ · · · ∪ Bk } ∈ Q is
Peither 1/n if p is inside one of the
balls of q, and 0 otherwise. The overall cost f¯(P, w, q) = p∈P f (p, q) is thus the fraction of
points of P that are covered by the union of these k balls.
13

The motivation of defining query spaces is usually to solve some related optimization
problem, such as the query q that minimizes the cost f¯(P, w, q). In this case, the requirement
to approximate every query in Q is too strong, and we may want to approximate only the
optimal query in some sense. To this end, we may wish to replace the set Q by a function
that assigns a different set of queries Q(S) for each subset S of P . For the correctness of
our results, we require that this function Q will be monotonic in the following sense: if T is
a subset of S then Q(T ) must be a subset of Q(S). If we wish to have a single set Q(P ) of
queries as above, we can simply define Q(S) := Q(P ) for every subset S of P , so that the
desired monotonic property will hold.
Example 2: For the set Q in (2), define Q(S) to be the set of balls in Rd such that the
center of each ball is a point in S. More generally, we can require that the center of each
ball will be spanned by (i.e., linear combination of) at most 10 points from S.
We now conclude with the formal definitions for the above discussion.
Definition 4.1 (weighted set). Let S be a subset of some set P and w : P → [0, ∞) be a
function. The pair (S, w) is called a weighted set.
Definition 4.2 (query space). [FL11]. Let Q be a function that maps every set S ⊆ P to
a corresponding set Q(S), such that Q(T ) ⊆ Q(S) for every T ⊆ S. Let f : P × Q(P ) → R
be a cost function. The tuple (P, w, Q, f ) is called a query space.
We denote the cost of a query q ∈ Q(P ) by
X
w(p)f (p, q).
f (P, w, q) :=
p∈P

4.2

(ε, ν)-Approximation

Consider the query space (P, w, Q, f ) in Example 1 for k = 1, and suppose that we wish to
compute a set S ⊆ P , such that, for every given ball B, the fraction of points in P that are
covered by B, are approximately the same as the fraction of the points that are covered in
S, up to a given small additive percentage error ε > 0. That is, for every ball B,
|P ∩ B| |S ∩ B|
−
≤ ε.
|P |
|S|
By defining the weight u(p) = 1/|S| for every p ∈ S, this implies that for every query
(ball) q = B,
f¯(P, w, q) − f¯(S, u, q) =

X 1
X 1
|P ∩ B| |S ∩ B|
· f (p, q) −
· f (p, q) =
−
≤ ε.
|P |
|S|
|P |
|S|
p∈P
p∈S

A weighted set (S, u) that satisfies the last inequality for every query q ∈ Q(S) is called an
ε-approximation for the query space (P, w, f, Q). Note that the above example assumes that
the maximum answer to a query q is f (p, q) ≤ 1. Otherwise, the error guaranteed by an
ε-approximation for a query q is ε maxp∈P |f (p, q)|.
14

The above inequalities implies that if a ball covers a fraction of at least ε points from P
(i.e., at least εn points), then it must cover at least one point of P . If we only ask for this
(weaker) property from S then S is called an ε-net. To obtain the new results of this paper,
we use a tool that generalizes the notion of ε-approximation and ε-net, but less common in
the literature, and is known as (ε, ν)-approximation [LLS01].
By letting a = f¯(P, w, q) and b = f¯(S, u, q), an ε-approximation implies that |a − b| ≤ ε
for every query q, and ε-net implies that b > 0 if a ≥ ε. Following [LLS00], we define below
a distance function that maps a positive real |a − b|ν for each pair of positive real numbers
a and b. A specific value ν will imply |a − b|ν ≤ ε, i.e., that S is an ε-approximation, and a
different value of ν will imply that S is an ε-net for P . This is formalized as follows.
Definition 4.3 ((ε, ν)-approximation [LLS00]). Let ν > 0. For every a, b ≥ 0, we define the
distance function
|a − b|
.
|a − b|ν =
a+b+ν
P
Let (P, w, Q, f ) be a query space such that p∈P w(p) = 1 and f : P → [0, 1]. For ε > 0, a
weighted set (S, u) is an (ε, ν)-approximation for this query space, if for every q ∈ Q(S) we
have
f (P, w, q) − f (S, u, q) ν ≤ ε
(3)
Corollary 4.4 ( [LLS00, HP11]). Let a, b ≥ 0 and τ, ν > 0 such that |a − b|ν ≤ τ . Let ε > 0.
Then
(i) (ε-sample/approximation). If τ = ε/4 and ν = 1/4 then
|a − b| ≤ ε.
(ii) (ε-net). If τ = 1/4 and ν = ε then
(a ≥ ε

⇒

b > 0).

(iii) (Relative ε-approximation). Put µ > 0. If ν = µ/2 and τ = ε/9 then
(1) a ≥ µ
(2) a < µ

4.3

⇒
⇒

(1 − ε)a ≤ b ≤ (1 + ε)a.
b ≤ (1 + ε)µ.

Constructing ε-Approximations

Unlike the notion of coresets in the next section, the idea of ε-approximations is known
for decades [VC71, Mat89]. In particular, unlike coresets, ε-approximations and (ε, ν)approximations in general, can be constructed using simple uniform random sampling of
the points. The size of the sample depends linearly on the complexity of the queries in the
sense soon to be defined in this section.
Intuitively, and in most practical cases, including the examples in this paper, this complexity is roughly the number of parameters that are needed to define a single query. For
15

example, a ball in Rd can be defined by d + 1 parameters: its center c ∈ Rd , which is defined
using d numbers and the radius r > 0 which is an additional number. A query of k balls is
similarly defined by k(d + 1) parameters. However, by the set theory we have |R| = |Rm | for
every integer c ≥ 1, which means that we can encode every m integers to a single integer.
We can thus always reduce the number of parameters that are needed to define a query from
m to 1 by redefining our cost function f and re-encoding the set of queries.
There are also natural examples of query spaces whose cost function f is defined by
one parameter, but the size of the sampling needed for obtaining an ε-approximation is
unbounded, e.g., the query space where f (p, q) = sign(sin(pq)) and P = Q = Rd , where
sign(x) = 1 if x > 0 and 0 otherwise; see details e.g. in [?]. Hence, a more involved
definition of complexity is needed as follows.
While the number of subsets from a given set of n points is 2n , its can be easily verified
that the number of subsets that can be covered by a ball in Rd is roughly nO(d) . The exponent
O(d) is called the VC-dimension of the family of balls in Rd . The following definition is a
simple generalization by [FL11] for query spaces where the query set is a function and not a
single set. The original definition of pseudo-dimension can be found e.g. in [LLS00] and is
very similar to the definition of VC-dimension given by [VC71] as well as many other similar
measures for the complexity of a family of shapes.
Definition 4.5 (dimension [FL11]). For a query space (P, w, Q, f ) and r ∈ [0, ∞) we define
range(q, r) = {p ∈ P | w(p) · f (p, q) ≤ r} .
The dimension of (P, w, Q, f ) is the smallest integer d such that for every S ⊆ P we have
{range(q, r) | q ∈ Q(S), r ∈ [0, ∞)} ≤ |S|d .
The main motivation of the above definition, is that it tells us how many samples we
need to take uniformly at random from the input set P , to get an (ε, ν)-approximation (S, u)
as follows.
Theorem
4.6 ([VC71, LLS00, FL11]). Let (P, w, Q, f ) be a query space of dimension d such
P
that p∈P w(p) = 1 and f : P → [0, 1], and let ε, δ, ν > 0. Let S be a random sample from
P , where every point p ∈ P is sampled independently with probability w(p). Assign a weight
1
for every p ∈ S. If
u(p) = |S|
1
|S| ∈ Ω(1) · 2
εν


 
 
1
1
d log
+ log
ν
δ

then, with probability at least 1 − δ, (S, u) is an (ε, ν)-approximation for the query space
(P, w, Q, f ).

16

5

Improved Coreset Framework

5.1

Improved (ε, ν)-approximations

In this section we show the first technical result of this paper:
P that the additive error
ε maxp∈P |f (p, q)| in (3) can be replaced by ε, and the assumption i=1 wi = 1 in Theorem 4.6
can be removed. The condition for this to work is that the total importance value t, as
defined below, is small. More precisely, the required sample size will be near-linear with
t. Also, the uniform sampling in Theorem 4.6 will be replaced by non-uniform sampling
with a distribution that is proportional to the importance of each point. This result is
essentially a generalization and significant simplification of the framework in [FL11] that
used ε-approximations for constructing coresets. Maybe more importantly: using the idea
of (ε, ν)-approximations we are able to show that the sample size is near linear in t while
in [FL11] it is quadratic in t. For some applications, such an improvement
means turning a
√
theoretical result into a practical result, especially when t is close to n.
We define the importance of a point as the maximum absolute weighted cost w(p)f (p, q)
of a point p, over all the possible queries q, i.e,
s(p) := w(p) max |f (p, q)|,
q∈Q(P )

and hope that this sum is small (say, constant or log n), in P
other words, that not all the
points are very important. More precisely, if the sum t = p∈P s(p) of these costs is 1,
then we prove below that a new query space (P, w0 , Q, f 0 ) can be constructed, such that an
(ε, ν)-approximation (S, u) for the new query space would imply
|f (P, w, q)| − |f (S, u, q)|

ν

≤ ε.

(4)

That is, the additive error in (3) is reduced as desired.
The new query space (P, w0 , Q, f 0 ) is essentially a re-scaling of the original weights by
their importance, w0 (p) = s(p). To make sure that the cost of each query will still be the
same, we need to define f 0 such that w(p)f (p, q) = w0 (p)f 0 (p, q). This implies f 0 (p, q) :=
w(p)f (p, q)/s(p). While the new cost f¯0 (P, w0 , q) is the same as the old one f¯(P, w, q) for
every query q, the maximum value of |f 0 (p, q)| is 1, by definition of s(p), even if |f (p, q)| is
arbitrarily large. Hence, the additive error ε|f 0 (p, q)| in (3) reduced to ε in (4).
More generally, an (ε, ν/t)-approximation for (P, w0 , Q, f 0 ) would yield (4). Using the
uniform sample construction of Theorem 4.6, this implies that to get (4) we need to increase
the sample size by a factor that is nearly linear in t.
Theorem 5.1.
• Let (P, w, Q, f ) be a query space where f (p, q) ≤ 1 for every p ∈ P and q ∈ Q(P ).
• Let s : P → (0, ∞) such that s(p) ≥ w(p) maxq∈Q(P ) f (p, q).
P
• Let t = p∈P s(p).
17

• Let w0 : P → [0, 1] such that w0 (p) := s(p)/t.
• Let f 0 : P × Q(P ) → [0, 1] be defined as f 0 (p, q) := w(p) ·

f (p, q)
.
s(p)

• Let (S, u) be an (ε, ν/t)-approximation for (P, w0 , Q, f 0 ).
• Let u0 (p) = u(p) ·

w(p)
w0 (p)

for every p ∈ S.

Then for every q ∈ Q(S),
|f (P, w, q) − f (S, u0 , q)|ν ≤ ε.

(5)

Proof. Put q ∈ Q(S).
|f (P, w, q) − f (S, u0 , q)|ν =

X

w(p)f (p, q) −

p∈P

= t

X s(p)

X

t

· w(p) ·

X

· f (p, q)

(6)
ν

f (p, q)
f (p, q)
u(p) · w(p) ·
−t
s(p)
s(p)
p∈S
X

w0 (p) · f 0 (p, q) − t

p∈P

=

s(p)
t

p∈S

p∈P

= t

X u(p)w(p)

X

u(p) · f 0 (p, q)

p∈S

w0 (p) · f 0 (p, q) −

p∈P

X

(7)
ν

u(p) · f 0 (p, q)

p∈S

ν

(8)
ν/t

≤ ε,

(9)

where (6) is by the definition of u0 , (7) is by the definition of f 0 , (8) follows since
|ta − tb|ν =

|a − b|
|ta − tb|
=
= |a − b|ν/t ,
ta + tb + ν
a + b + ν/t

(10)

and (9) is by (3) and the definition of S in the sixth bullet of the statement.
Plugging Theorem 4.6 in Theorem 5.1 yields our main technical result that would imply
smaller coresets in the next sections.
Theorem 5.2. Let (P, w, Q, f ) be a query space of dimension d, where f is non-negative,
and let ε, δ, ν > 0. Let S be a random sample from P , where every point p ∈ P is sampled
tw(p)
independently with probability s(p)/t. Assign a weight u0 (p) = s(p)·|S|
for every p ∈ S. If

 
 
t
1
t
d log
+ log
|S| ∈ Ω(1) · 2
εν
ν
δ
then, with probability at least 1 − δ,
|f (P, w, q) − f (S, u0 , q)|ν ≤ ε.
18

0
0
Proof. By Theorem
P 5.10 it suffices to prove that S is an (ε, ν/t)-approximation for (P, w , Q, f ).
Indeed, since
p∈P w (p) = 1 and S is a random sample where p ∈ P is sampled with
0
probability w (p) = s(p)/t, the weighted set (S, u) is, with probability at least 1 − δ, an
(ε, ν/t)-approximation for (P, w0 , Q, f 0 ) by Theorem 4.6.

While Theorem 4.6 suggests to construct an ε-approximation simply by taking a uniform
random sample from P , Theorem 5.2 requires us to take non-uniform sample where the
distribution
P is defined by the importance s(·). Bounding these importances such that the
sum t = p∈P s(p) of importances will be small raise a new optimization problem that, as
we will see in the next sections, might be not easy at all to solve. So what did we gain by
proving Theorem 4.6 (or the framework of [FL11] in general)?
Most of the existing papers for constructing coresets essentially had to bound the total
importance of the related problem anyway. However, without Theorem 4.6, their proofs also
had to deal with complicated terms that involved ε and include sophisticated probability
arguments. Essentially each paper had to re-prove in some sense the very involved proofs
in [VC71, LLS00] that were researched for decades, as well as the mathematics behind the
usage of Definition 4.3. Beside the complicated proofs, the final bounds, the dependency on
ε and δ, as well as the coreset construction algorithm, were usually sub-optimal compared
to Theorem 4.6. On the contrary, bounding the total importance allows us to focus on a
deterministic results (no δ involved) and the terms s(p) can be approximated up to constant
factors (and not (1 + ε)-factors). This is demonstrated in Section 6.

5.2

Improved Coreset Framework

While the result in the previous section can be used to improve the quality of (ε, ν)approximation in general, its main motivation is to construct the more recent type of data
structures that is sometimes called coresets.
As explained in Theorem 4.6, ε-approximation and (ε, ν)-approximation in general, can
be easily computed using uniform random sampling. While ε-approximations are useful for
hitting sets, where we wish to know how much points are covered by a set of shapes, they
are less relevant for shape fitting, where we wish to approximate the sum of distances of the
input points to a given query shape or model. The main reason is that in shape fitting the
maximum contribution of a point to the overall cost (sum of covered points) is bounded by 1,
while in general its distance to the shape is unbounded. Using the notation of the previous
section: the importance of each point is the same, so Theorem 5.2 yields uniform sample S
and the same error as Theorem 4.6.
Example: In the Euclidean k-means problem the input is a set P of n points in Rd .
For simplicity, assume that the set is unweighted, that is w(p) = 1/n for every p ∈ P . A
query in this context is a set of points (centers) in Rd so Q(P ) is the family (set) of k balls
in Rd . We denote the squared Euclidean distance from a point p ∈ P to a center c ∈ Rd by
D(p, c) = kp − ck22 . The distance from p to its closest center in C = {c1 , · · · , ck } is then
D(p, C) := min D(p, c) = min kp − ck22 .
c∈C

c∈C

19

By defining g(p, C) = D(p, C) we obtain the query space (P, w, g, Q), where ḡ(P, w, C) is the
average squared distances to a given (query) set C of k centers. Our goal is to compute a
weighted subset (S, u) such that ḡ(P, w, C) will approximate the average squared distances
ḡ(S, u, C) for every set of centers.
Suppose that (S, u) is an ε-approximation of (P, w, Q, g). Then
|ḡ(P, w, C) − ḡ(S, u, C)| ≤ ε max g(p, C).
p∈P

(11)

That is,
X
1X
D(p, C) −
u(p)D(p, C) ≤ ε max D(p, C).
p∈P
n p∈P
p∈P

(12)

In other words, the additive error depends on the maximum distance between a point to a
given center, which can be arbitrary large, unless we assume, say, that both the input points
and centers are inside the unit cube. Theorem 5.2 would not improve this bound by itself,
since the importance of each point, maxc∈Q(P ) D(p, C) is unbounded.
In this section we wish to compute a weighted subset (S, u) that will approximate the
average distance ḡ(P, w, C) for every query, up to a multiplicative factor of 1 ± ε without
further assumptions. Such a set is sometimes called a coreset as follows.
Definition 5.3 (ε-coreset). For an error parameter ε ∈ (0, 1), the weighted set (S, u) is an
ε-coreset for the query space (P, w, Q, g) if S ⊆ P and for every q ∈ Q(S),
(1 − ε)ḡ(P, w, q) ≤ ḡ(S, u, q) ≤ (1 + ε)ḡ(P, w, q).
An equivalent definition of an ε-coreset, is that the additive error ε maxp∈P g(p, C) in (11)
is replaced by εḡ(P, w, q). That is, for every q ∈ Q(S)
|ḡ(P, w, q) − ḡ(S, u, q)| ≤ εḡ(P, w, q).
ε-coreset implies that not only the average, but also the sum of squared distances (or sum
of costs, in general) are preserved up to 1 ± ε. Note also that simply multiplying the weight
of each point in S by (1 − ε) would yield a one sided error,
ḡ(P, w, q) ≤ ḡ(S, u, q)| ≤

2ε
1+ε
· ḡ(P, w, q) = 1 +
· ḡ(P, w).
1−ε
1−ε

If we assume in addition that ε ∈ (0, 1 − 2/c) for some c > 2, then 1 − ε ≥ 2/c and thus
ḡ(P, w, q) ≤ ḡ(S, u, q)| ≤ (1 + cε)ḡ(P, w).
Hence, an ε/c-coreset (S, u) implies,
ḡ(P, w, q) ≤ ḡ(S, u, q) ≤ (1 + ε)ḡ(P, w, q).

(13)

For example, if ε ∈ (0, 1/2), then an (ε/4)-coreset would yield (13), by substituting c = 4.
20

The main observation for getting the desired multiplicative approximation is that a multiplicative approximation of εḡ(P, w, q), can be turned into an additive approximation of ε
by replacing g(p, q) with its scaled version
f (p, q) =

g(p, q)
.
ḡ(P, w, q)

To get a coreset for g as in (13) it suffices to have an ε-approximation for f . In addition, for
many problems, while the importance s(p) of a point in p is unbounded with respect to g, it
is bounded with respect to f . We formalize this observation as follows.
Corollary 5.4. Let (P, w, Q, g) be a query space for some non-negative function g, and
define f : P × Q(P ) → R such that for every p ∈ P and q ∈ Q(P ) we have
f (p, q) =

g(p, q)
.
ḡ(P, w, q)



1
Let (S, u ) be a ε/4,
-approximation for (P, w0 , Q, f 0 ) as defined in Theorem 5.1. Then
4t
(S, u0 ) is an ε-coreset for (P, w, Q, g), i.e., for every q ∈ Q(S)
0

|ḡ(P, w, q) − ḡ(S, u, q)| ≤ εḡ(P, w, q).
Proof. Put q ∈ Q(S), τ = ε/4 and ν = 1/4. Applying Theorem 5.1 with ν yields
|f¯(P, w, q) − f¯(S, u0 , q)|ν ≤ τ.
Since f¯(P, w, q) = 1, this implies
1−

ḡ(S, u0 , q)
ḡ(P, w, q)

= |f¯(P, w, q) − f¯(S, u0 , q)|ν ≤ τ.
ν

Substituting a = 1, and b = f¯(S, u0 , q) in Corollary 4.4(i) yields
1−

ḡ(S, u0 , q)
≤ ε.
ḡ(P, w, q)

Multiplying by ḡ(P, w, q) yields an ε-coreset as
|ḡ(P, w, q) − ḡ(S, u0 , q)| ≤ εḡ(P, w, q).

Combining Corollary 5.4 and Theorem 4.6 yields the following theorem.

21

Theorem 5.5. Let (P, w, Q, g) be a query space, where g is a non-negative function. Let
s : P → (0, ∞) such that
w(p)g(p, q)
s(p) ≥ max P
,
q∈Q(P )
p∈P w(p)g(p, q)
P
and t = p∈P s(p). Let d be the dimension of (P, w0 , Q, f 0 ) as defined in Theorem 5.4. Let
c ≥ 1 be a sufficiently large constant, and let S be a random sample of

 
ct
1
|S| ≥ 2 d log t + log
,
ε
δ
points from P , such that for every p ∈ P and q ∈ S we have p = q with probability s(p)/t.
t·w(p)
Let u0 (p) = s(p)|S|
for every p ∈ S. Then, with probability at least 1 − δ, (S, u0 ) is an ε-coreset
for (P, w, Q, g), i.e.,
∀Q ∈ Q(S) : (1 − ε)g(P, w, q) ≤ g(S, u0 , q) ≤ (1 + ε)g(P, w, q).

6

Application: Smaller and Generalized Coreset for kMeans

Definition 6.1 (ρ-metric space). Let X be a set, and D : X 2 → [0, ∞) be a function. Let
ρ > 0. The pair (X, D) is a ρ-metric space if for every (x, y, z) ∈ X 3 we have
D(x, z) ≤ ρ(D(x, y) + D(y, z)).
For C ⊆ X we denote D(x, C) := minc∈C D(x, c) assuming that such minimum exists.
Moreover, for ε > 0, the pair (X, D) is a (ψ, ε)-metric space if for every (x, y, z) ∈ X 3 we
have
D(x, y)
+ εD(y, z)).
|D(x, z) − D(y, z)| ≤
ψ
Note that for every x, y ∈ X, and a center cy ∈ C that is closest to y, we have
D(x, C) = min D(x, c) ≤ D(x, cy ) ≤ ρ(D(x, y) + D(y, cy )) = ρ(D(x, y) + D(y, C)).
c∈C

(14)

A simple way to decide whether (D, X) is indeed a ρ-metric space is to use the following
bound, known as Log-Log Lipschitz, that is usually easier to compute. The following lemma
is very similar to [FS12, Lemma 2.1], where in the proof of (i) the constant 4 that appeared
there is replaced here by 1/ε.
Lemma 6.2 (Lemma 2.1(ii) in [FS12]). Let D̃ : [0, ∞) → [0, ∞) be a monotonic nondecreasing function that satisfies the following (Log-Log Lipschitz) condition: there is r > 0
such that for every x > 0 and ∆ > 1 we have
D̃(∆x) ≤ ∆r D̃(x).
22

Let (X, dist) be a metric space, and let D(x, y) = D̃(dist(x, y)) for every x, y ∈ X. Then
(X, D) is a ρ-metric space for ρ = max {2r−1 , 1}, and a (ψ, ε)-metric space for every ε ∈ (0, 1)
and ψ = (r/ε)r .
For example, consider a metric space (X, dist) and the function D̃(x) = x2 that corresponds to the squared distance D(p, q) = D̃(dist(p, q)) = (dist(p, q))2 . Note that (X, D) is
not a metric space since the triangle inequality does not hold. However, for every x > 0 and
∆>1
D̃(x∆) = (x∆)2 = ∆2 x2 = ∆2 D̃(x).
Hence, by Lemma 6.2 r = 2, ρ = 2r−1 = 2, (D, X) is a 2-metric space and a (ψ, ε)-metric
space for ψ = (2/ε)2 .
To compute a coreset for a problem we need to decide what are the important points,
or more formally, to use Theorem 5.5 we need to bound the importance s(p) of each point
p ∈ P . To do this, we usually need to solve another optimization problem that is usually
related to computing the query with the minimal cost. For example, the bound of importance
of a point in the k-mean problem, as will be defined later, is based on the optimal solution
for the k-means problem. Unfortunately, this optimization problem is usually hard and the
main motivation for constructing the coreset in the first place.
There are two leeways from this chicken-and-egg problem:
(i) Use the merge-and-reduce approach that reduces the problem of computing a coreset
n
small weighted
for a large set of n items to the problem of computing coresets for 2|S|
(core)sets. Each input coreset 2|S| is reduced to a coreset of |S| and merged with
another such coreset, where 2|S| is the minimum size of input set that can be reduced
to half using the given coreset construction. In this case, if this coreset construction
takes time f (|S|) than, since there are such O(n) constructions, the overall running
time will then be O(n) · f (|S|).
(ii) For problems such as k-means, it is NP-hard to compute the optimal solution, even
for a small set of n = O(k) points. Instead of computing an optimal solution, usually a constant factor approximation suffices for computing the importance of each
point. Since for many problems such an approximation is also unknown or too slow to
compute, (α, β)-approximation, (or bi-criteria, or bicriterion), can be used instead as
explained below.
Suppose that the optimal solution for the k-means problem on a set is OP Tk . That is,
there is a set C of k centers whose cost (sum of squared distances to the input points of P )
is OP Tk , and there is no such set of smaller cost. Then a set C̃ is called an α-approximation
if its cost is at most α · OP T . However, for many problems, even a rougher approximation
would do: instead of using k centers to approximate OP Tk by a factor of α, we use βk
centers, where each input point may be assigned for the nearest center. Note that the cost is
still compared to OP Tk and not to OP Tβk . We define (α, β)-approximation formally below.
For our purposes later, we generalize this common definition of (α, β)-approximation, and
23

allow a point to be assigned to a different center than its nearest one, as long as the overall
cost is small.
Definition 6.3 ((α, β)-approximation.). Let (X, D) be a ρ-metric space, and k ≥ 1 be an
integer. Let (P, w, Q, g) be a query space such that P ⊆ X, Q(P ) = {C ⊆ X | |C| ≤ k}, and
g : P × Q(P ) be a function such that
g(p, C) = D(p, C) := min D(p, c).
c∈C

(15)

Let α, β ≥ 0, B ⊆ X such that |B| ≤ βk and B : P → B such that
X
w(p)g(p, {B(p)}) ≤ α min ḡ(P, w, q).
q∈Q(P )

p∈P

Then B is called an (α, β)-approximation for (P, w, Q, g).
For every b ∈ B, we denote by Pb = {p ∈ P | B(p) = b} the points that are mapped to
the center b. We also denote p0 = B(p) for every p ∈ P .
One of the main tools in our novel streaming algorithm is also a technique to update an
(α, β)-approximation. However, due to memory limitations, our streaming algorithm cannot
attach each point to its nearest center, but still the distances to the approximated centers
is bounded. We thus generalize the definition of an (α, β)-approximation, which is usually a
set B of size βk, to a function B : P → B that assigns each input point to a (not necessarily
closest) center in B, while the overall cost is still bounded.
Since an (α, β)-approximation yields a weaker result compared to a PTAS or α-approximation,
it can usually be computed very quickly. Indeed, a very general framework for constructing
(α, β)-approximation to any query space with a small VC-dimension is suggested in [FL11]
where α = 1 + ε and β = O(log n).
Reducing (α, β)-approximation to an α = O(1) approximation. The size of the
coreset usually depends on α and β. However, if they are reasonably small (e.g. polynomial
in log n), we can reduce the approximation factor and number of centers in few phases as
follows: (i) Compute an (α, β)-approximation for small (but maybe not constant) α and β.
(ii) Compute an ε-coreset for ε = 1/2 using this approximation. (iii) Compute an O(1) factor
approximation on the coreset. Since the coreset is small, such an approximation algorithm
can run inefficiently, say, in polynomial time if the coreset is of size (log n)O(1) . The resulting
O(1) approximation for the coreset is also an O(1) approximation for the original set, by
the definition of the ε = 1/2 coreset. (iv) Recompute the coreset for the complete (original)
data using the O(1) approximation instead of the (α, β)-approximation to obtain a coreset
of size independent of both n and d.
Assumption 6.4. In what follows we assume that:
• P is a set that is contained in X where (X, D) is a ρ-metric space and (φ, ε)-metric
space as defined in (6.1).
24

Algorithm 2: Coreset(P, w, B, m)
Input:
A weighted set (P, w) where P ⊆ X and (D, X) is a ρ-metric space,
(α, β)-approximation B : P → B,
and sample size m ≥ 1.
Output:
A pair (S, u) that satisfies Theorem 6.6.
1
2
3

4

5
6
7
9

for each b ∈ B do
Set Pb ← {p ∈ P | B(p) = b}
for each b ∈ B and p ∈ Pb do
w(p)
w(p)D(p, B(p))
P
.
+
Set Prob(p) ← P
2 q∈P w(q)D(q, B(q)) 2|B| q∈Pb w(q)
Pick a sample S of at least m points from P such that for each q ∈ S and p ∈ P we
have q = p with probability Prob(p).
for each p ∈ P do
w(p)
Set u(p) ← |S|·Prob(p)
.
Set u(p) ← 0 for each p ∈ P \ S. /* Used only in the analysis.
*/
return (S, u)

• (P, w, Q, g) is a query space as defined in (15).
• dk denotes the dimension of f 0 as defined in Corollary 5.5.
• c is a sufficiently large constant that can be determined from the proofs of the theorems.
• B is an (α, β)-approximation for (P, w, Q, g) as in Definition 6.3.
• We are given an error parameter ε ∈ (0, 1).
• We are given a maximum probability of failure δ ∈ (0, 1).
We begin with the following claim, that is a simplified and generalized version of a similar
claim in [LS10].
Lemma 6.5. For every b ∈ B, and p ∈ Pb have
w(p)D(p, C)
ραw(p)D(p, p0 )
ρ2 (α + 1)
P
≤P
+
.
0
q∈P w(q)D(q, C)
q∈P w(q)D(q, q )
q∈Pb w(q)

P

Proof. Put p ∈ P and b ∈ B such that p ∈ Pb . We need to bound
w(p)D(p, C)
ρw(p)D(p, p0 )
ρw(p)D(p0 , C)
≤P
+P
q∈P w(q)D(q, C)
q∈P w(q)D(q, C)
q∈P w(q)D(q, C)

P

αρw(p)D(p, p0 )
ρw(p)D(p0 , C)
P
P
≤
+
,
0
q∈P w(q)D(q, q )
q∈P w(q)D(q, C)
25

(16)

where the first inequality holds by (14), and the second inequality holds since B is an (α, β)approximation.
To bound the last term, we sum the inequality D(p0 , C) ≤ ρ(D(p0 , q) + D(q, C)) over
every q ∈ Pb to obtain
X
X
X
D(p0 , C)
w(q) =
w(q)D(p0 , C) ≤
w(q) · ρ(D(p0 , q) + D(q, C))
q∈Pb

q∈Pb

=ρ

q∈Pb

X

w(q)D(q 0 , q) + ρ

q∈Pb

≤ ρα

X

X

w(q)D(q, C)

q∈Pb

w(q)D(q, C) + ρ

q∈Pb

≤ ρ(α + 1)

X

w(q)D(q, C)

q∈Pb

X

w(q)D(q, C).

q∈P

Dividing by

P

q∈Pb

w(q) ·

P

q∈P

w(q)D(q, C) yields

ρ(α + 1)
D(p0 , C)
≤P
.
q∈P w(q)D(q, C)
q∈Pb w(q)

P

Substituting this in (16) yields the desired result
ραw(p)D(p, p0 )
ρ2 (α + 1)
w(p)D(p, C)
P
≤P
+
.
0
q∈P w(q)D(q, C)
q∈P w(q)D(q, q )
q∈Pb w(q)

P

Our first theorem suggests a coreset for k-means of size near-quadratic in k and quadratic
in ε, based on our improved framework and last lemma. Existing work [LS10] for obtaining
such coresets with only positive weights requires size cubic in k.
Theorem 6.6. Under Assumption 6.4, let t = k · ρ2 (α + 1)β. Let (S, u) be the output of a
call to algorithm Coreset(P, w, B, m), where

 
1
ct
.
m ≥ 2 dk log t + log
ε
δ
Then, with probability at least 1 − δ, (S, u) is an ε-coreset of size m for (P, w, Q, g).
Proof. Let f : P → [0, ∞) such that for every C ∈ Q(P ),
D(p, C)
,
q∈P w(q)D(q, C)

f (p, C) = P
and define

ραw(p)D(p, p0 )
ρ2 (α + 1)
P
+
.
s(p) = P
0
q∈P w(q)D(q, q )
q∈Pb w(q)
26

By Lemma 6.5,
X
XX
XX
max w(p)|f (p, C)| =
max w(p)f (p, C) ≤
s(p)
p∈P

C∈Q(P )

b∈B p∈Pb

C∈Q(P )

b∈B p∈Pb

2

= ρα + |B| · ρ (α + 1) = ρ(α + |B|α + |B|) ∈ O(ρ2 (α + 1)βk).
Applying Theorem 5.5 with the query space (P, w, Q, g) then yields the desired bound.
Our second theorem in this section suggests a coreset for k-means of size near-linear in
k by combining new observations with our improved framework.
Theorem 6.7. Under Assumption 6.4, let t = α/φ. Let (S, u) be the output of a call to
algorithm Coreset(P, w, B, m), where

 
ck(t + β)
1
m≥
d log t + log(βk) + log
.
2
ε
δ
Then, with probability at least 1 − δ, (S, u) is an ε-coreset of size m for (P, w, Q, g).
Proof. Let



2D(p, p0 )
0
H = p ∈ P | |D(p, C) − D(p , C)| ≤
.
φ
We need to bound by O(ε) the expression
P
p∈P (w(p) − u(p)) · D(p, C)
P
q∈P D(q, C)
P
0
p∈H (w(p) − u(p)) · (D(p, C) − D(p , C))
P
≤
q∈P D(q, C)
P
X q∈P w(q)D(b, C)
P b
+
·
r∈P w(r)D(r, C)
b∈B
X (w(p) − u(p)) · D(p, C)
X (w(p) − u(p)) · (D(p, C) − D(p0 , C))
P
P
−
w(q)D(b,
C)
q∈P
q∈Pb w(q)D(b, C)
b
p∈P
p∈P ∩H
b

(17)

(18)
(19)
!
.

(20)

b

Put b ∈ B and p ∈ Pb . Then,
X (w(p) − u(p)) · D(p, C)
X (w(p) − u(p)) · (D(p, C) − D(p0 , C))
P
P
−
q∈Pb w(q)D(b, C)
q∈Pb w(q)D(b, C)
p∈P
p∈P ∩H
b

b

X w(p) − u(p)
P
=
q∈Pb w(q)
p∈P

(21)

b

+

X (w(p) − u(p)) · (D(p, C) − D(b, C))
P
.
q∈Pb w(q)D(b, C)

(22)

p∈Pb \H

We now prove that, with probability at least 1 − cδ, each of the expressions (18), (21)
and (22) is bounded by 2ε.
27

Bound on (18):
every p ∈ P ,

Let h(p, C) =

D(p,C)−D(p0 ,C)
P
q∈P w(q)D(q,C)

if p ∈ H and h(p, C) = 0 otherwise. For

2αw(p)D(p, p0 )
2αProb(p)
2w(p)D(p, p0 )
P
P
≤
≤
.
w(p) · |h(p, C)| ≤
0
φ q∈P w(q)D(q, C)
φ q∈P w(q)D(q, q )
φ
Hence, using t = 2α/φ in Theorem 5.5, with probability at least 1 − δ,
P

p∈H (w(p)

− u(p)) · (D(p, C) − D(p0 , C))
X
P
(w(p) − u(p))h(p, C) = |h̄(P, w, q)−h̄(S, u, q)| ≤ ε,
=
q∈P D(q, C)
p∈H
(23)

Bound on (21):

Let I(p, b) = 1/

P

q∈Pb

w(q) if p ∈ Pb and I(p, b) = 0 otherwise. We have

max w(p)I(p, b) ≤ 2|B|Prob(p),
p∈P

P
where Prob(p) is defined in Line 3 of Algorithm 2. Hence, p∈P max w(p)I(p, b) ≤ 2|B|.
Also, each point in S is sampled with probability proportional to 2|B|Prob(p) and for c ≥ 2,

 



2|B|
1
2|B|
|B|
|S| ≥ 2
2 log(2|B|) + log
= 2
log(|B|) + log
.
ε
δ
ε
δ
Using d = 1 and replacing δ by δ/|B| in Theorem 5.5 yields that S is an ε-coreset for
(P, w, {b} , I), with probability at least 1 − δ/|B|. That is
|

X w(p) − u(p)
¯ w, b) − I(S,
¯ u, b)| ≤ ε.
P
| = |I(P,
w(q)
q∈P
b
p∈P

(24)

b

By the union bound, with probability at least 1 − δ, the last inequality holds for every b ∈ B
simultaneously.
Bound on (22):

Since (D, X) is a ρ-metric,

D(p, b)
|D(p, C) − D(b, C)| ≤
+ εD(b, C) ≤ max
φ
Hence




2D(p, b)
, 2εD(b, C) .
φ

X (w(p) − u(p)) · (D(p, C) − D(b, C))
X w(p) + u(p)
P
P
≤ 2ε
.
w(q)D(b,
C)
q∈Pb
q∈Pb w(q)

p∈Pb \H

(25)

p∈Pb \H

The last expression is bounded using (24), as
X w(p) + u(p)
X (1 + ε)w(p)
P
P
≤
= (1 + ε).
q∈Pb w(q)
q∈Pb w(q)

p∈Pb \H

p∈Pb \H

28

(26)

Bounding (17):
least 1 − 10δ,
P

p∈P (w(p)

P

By combining the above bounds we obtain that with probability at

− u(p)) · D(p, C)

q∈P

D(q, C)

≤ ε + (ε + 2ε(1 + ε))

X
b∈B

P
w(q)D(b, C)
Pq∈Pb
≤ cε.
r∈P w(r)D(r, C)

Replacing ε with ε/c, and δ with δ/c then proves the theorem.

The values of ρ and ψ. Algorithm 2 can be applied to compute a coreset for any given
variant of the k-means/median problem given a set P and a ρ or (ρ, ε)-metric (X, D). The
only difference is the size of the required coreset. The parameters ρ and φ can usually be
computed easily using Lemma 6.2. For example, in the case of distances to the power of
r ≥ 1, the value of ρ is roughly 2r and the value of φ is roughly εr . For most common
m-estimators the values are similar, when r is some constant.
The values of α and β can be α = β = O(1) for k-means and all its variants, by using
the generic algorithm for computing (α, β)-approximation in [FL11] with α = β = O(log n),
and then use the technique for reducing α and β. For bounding the approximation of the
bi-criteria approximation in [FL11] only the pseudo-dimension of the problem is required,
which is usually easy to compute as explained below.
Dimension d. Unlike the total sensitivity t, the dimension d for numerous problems was
already computed in many papers in computational geometry and machine learning (in the
context of PAC learning). This includes reductions and connections to similar notions such
as the shattering dimension or the VC-dimension of a set. General techniques for computing
the dimension of a set based on number of parameters to define a query, or number of
operations that are needed to answer a query can be found in the book [AB99].
Dimension of the k-means problem and its variants. Note that, unlike sensitivity,
the dimension is less dependent on the exact type of distance function. For example, using
Euclidean distance or Euclidean distance to the power of 3 as a variant for the k-means
clustering problem does not change the dimension. This is because set of ranges for both of
these problems is the same: subsets of the input points that can be covered by k balls. It is
easy to compute the dimension for the k-means problem (the query space (P, w, Q, g) in our
paper, as well as the modified query space for the function f . These bounds can be found
in [FL11]. In addition, [FL11] provide a simple reduction that shows that the dimension for
k centers is the same as the dimension of 1 center multiplied by k.
In short, for the Euclidean space, the dimension of the k-means/median problem is O(dk),
and for metric spaces (graph) the dimension is O(k log n).

29

Smaller coreset for k-means queries.
Consider the k-means queries in Rd , i.e., the
P
cost is the sum of squared distances
p∈P D(p, C) over every point in P to its nearest
d
center in a given set C of k points in R . It was proven that projecting P onto an O(k/ε)dimensional subspace that minimizes its sum of squared distances, known as the low-rank
approximation of P , would preserve this sum, to any set of k centers, up to a factor of 1 ± ε.
This is in some sense an ε-coreset for P of size n and that is not subset of the input, but of
low-dimensionality. In particular, this result implies that there is a set of centers (known as
centroid set [HPK07]) that is contained in a O(k/ε)-dimensional space, such that every set
of k centers in Rd can be replaced by a set of k centers in the centroid set, that would yield
the same cost up to a factor of 1 ± ε. In particular, this implies that the dimension of the
k-means problem can be reduced to O(k/ε) instead of dk, i.e., independent of d. Combining
this result with our paper yields the first coreset for k-means of size independent of d that
is subset of the input (in particular, preserve the sparsity of the input points) that also
supports streaming.
Weak coresets of size independent of d. For the non-Euclidean case or non-squared
distances it seems that it is impossible to obtain coreset of size independent of d. However,
coreset as defined above (sometimes called strong coreset) approximates every query in
the set of queries, while the main application and motivation for constructing coreset is to
compute the optimal query or its approximation. A weak coreset is a small set that can
be used to give such an approximation. The exact definition of weak coreset also changes
from paper to paper. In particular, a weak coreset for k-means was suggested in [FMS07].
However, to extract the approximated solution from the coreset we must run an exhaustive
search and cannot use existing algorithms or heuristics as in the case of strong coreset.
In this paper, following [FL11], we use a simple and general definition of weak coreset,
that is also more practical. Instead of defining a unique (static, global) set Q of queries,
we define Q to be a function that maps every subset (potential coreset) S of P to a set
of queries. A weak coreset S needs only to approximate the queries in Q(S). It turns out
that for many case the (generalized definition) of dimension for such a query space is much
smaller compared to the traditional case where Q(S) = Q(P ) is the same for every subset
S ⊆ P.
To be able to use this property in our existing proofs, we require a monotonicity property:
that the set Q(T ) of queries that are assigned to a subset T ⊂ S must be contained in Q(S).
If we can prove that for every S ⊆ P , the set Q(S) contains a (1 + ε)-approximation to both
the optimal solution of S and P , then we can extract such an approximation from S. For
k-means, k-median and their variants, it was proven in [SV07] that the optimal k-centers of a
set S can be approximated by such a set that is spanned by O(k/ε) points in S. By defining
Q(S) to be the union of the optimal center of P , with all the centers that are spanned by
O(k/ε) points in S, we get that S is a weak coreset. Note that the definition of Q can be
explicit (without knowing the optimal center of P ) and is needed only to bound its dimension
as in Definition 4.5.

30

Extracting a (1 + ε)-approximation from the coreset. It was proven in [FL11] that
for problems such as k-median such weak coresets have dimension O(k/ε), i.e., independent
of d. Unlike [FL11], we suggest here a very simple way to extract the approximated solution
from the coreset: compute the weighted coreset (S, u) (of size independent of d) as defined in
Algorithm 2, and then use any given algorithm to compute a (1 + ε) approximation set C of
k centers on the coreset (or any other trusted heuristic that we hope computes such a set C).
Since C is not necessarily spanned by few points in S, it may not be a good approximation for
the original set P and we should not return it. However,the proof in [SV07] is constructive
and shows a near-linear time algorithm that using such an approximated solution set C, we
can compute another (1 + O(ε))-approximation C 0 that has the additional property that C 0
is spanned by O(k/ε) points in S. Hence, C 0 is both a near-optimal solution to S and in
Q(S), so it must be a (1 + ε)-approximation for the optimal solution of P .

7

Appendix A: Merge and Reduce Tree

We now briefly introduce the previous technique for maintaining coresets in the streaming
setting due to Har-Peled and Mazumdar [HPM04] and Bentley and Sax [BS80]. In this
method, a merge-and-reduce tree is built by using an offline coreset construction as a blackbox. Previously merge-and-reduce was the only known technique for building a streaming
coreset for metric k-median, and it relies solely on the following two properties:
1. Merge: The union of (k, )-coresets is a (k, )-coreset.
2. Reduce: A (k, )-coreset of a (k, δ)-coreset is a (k,  + δ)-coreset.
The merge-and-reduce tree works as follows. There are buckets Bi for i ≥ 0. In each step,
the bucket B0 takes in a segment of O(1) points from the stream. Then the tree works like
counting in binary: whenever buckets B0 to Bi−1 are full, these i buckets are merged and
then reduced by taking a (k, log n )-coreset and storing the result in Bi .
Let s be the space of offline construction, which depends on  as −a . At the end of
the stream, O(log n) buckets have been used and each bucket uses O(s loga n) space; this
incurs a multiplicative overhead of Θ(loga+1 n) in the storage requirement. The second factor
comes from using the accuracy parameter log n , which is necessary by Property 2 since the
construction will be compounded O(log n) times. Due to this compounding, the runtime is
multiplied by a factor of O(log n).

8

Appendix B: General Streaming Reduction

We present a general technique for converting an offline coreset construction to a streaming coreset construction with O(log n) overhead. Given a ρ-metric space (X, D) (recall
Definition 6.1), we build a query space (P, w, g, Q) in the same way as in Definition 6.3:
Q(P ) = {C ⊆ X | |C| ≤ k} and g(p, C) = D(p, C) := minc∈C D(p, c). Here k is a positive
integer that denotes how many centers may be used for the clustering.
31

Our bicriterion algorithm is an adjustment to the algorithms of [BMO+ 11] and [CCF]
with the following important difference: our bicriterion is online, so we do not delete and
reassign centers as is [BMO+ 11]. This “online” property is critical for the algorithm to work
and is one of the main technical ideas. Although a fully online bicriterion can require linear
space, we maintain a division of the stream P into a prefix R and a suffix P \ R such that
our bicriterion is online on the suffix P \ R and the prefix R can be largely ignored. To
maintain this property of being online for the suffix, we incur only a modest space increase
from O(k log n) to O(log( 1 )k log n).
After having an online bicriterion (it is further explained below why this property is essential), the offline coreset algorithms perform non-uniform sampling procedure with carefully
chosen probabilities that are defined by the bicriterion. Equipped with our new bicriterion algorithm, implementing the sampling procedure is rather straightforward computation
which is explained in Section 8.2. As a result, we can implement any sampling-based coreset
algorithm for k-median without merge-and-reduce and in one pass. As such it is applicable
to several coreset constructions (such as k-means and other M -estimators). In addition,
we believe that our methods will work with other objective functions as well such as (k, j)subspace, and we hope that future work will investigate these directions.
Many clustering algorithms [COP03, GMM+ 03, BMO+ 11] maintain a weighted set (B, u)
of points (which are selected using a facility-location algorithm). Upon arrival of an update
(p, w(p)) from the stream, this update is added to the set by u(p) ← u(p) + w(p).
In these algorithms, only a single operation is performed on (B, u) which we call MOVE.
For two points p, p0 ∈ B with weights u(p) and u(p0 ), the function MOVE(p, p0 ) does the
following: u(p0 ) ← u(p0 ) + u(p) and u(p) ← 0. This essentially moves the weight at location
p to the location of p0 . The motivation for building (B, u) will be compression; (B, u) will
be maintained over the stream P in such a way that |B| = O(log |P |).
Throughout this section, we will assume for ease of exposition that each point in the
stream is from a distinct location. This simplifies the analysis, allowing B to take only a
single value for each input point. The algorithm works for general inputs without modification, requiring only a more careful notation for the analysis. Additionally, we state the
algorithm for unweighted input (where w(p) = 1 for all p ∈ P ) and the parameter n is
defined as |P |. We still include w(p) throughout the algorithm and analysis, as the analysis
generalizes to weighted inputs where the parameter n is replaced by the sum of all weights
(after normalizing the minimally weighted point to 1).

8.1

An Algorithm for building a coreset

First, let us describe how the algorithm of [BMO+ 11] works. We will modify this algorithm
as part of our coreset construction. In this summary, we alter the presentation from that
of [BMO+ 11] to more fluidly transition to our modified version but the algorithm remains the
same. [BMO+ 11] operates in phases i ≥ 1. This means that the algorithm maintains a phase
number i (used internally by the algorithm), beginning in phase i = 1. As the stream arrives,
the algorithm may decide to increment the phase number. Let (Ri , wi ) denote the prefix of
the input received before the end of phase i, and let OPTk (Ri ) denote the minimal value of
32

Figure 1: A bicriterion approximation. Although connections are not optimal, the sum of
all connections costs is O(OPT).
ḡ(Ri , wi , C) over every C ∈ Q. When phase i + 1 begins, a value Li+1 is declared on Line 27
as a lower-bound for the current value of OPTk (Ri+1 ). The algorithm has computed (Mi , ui ),
which we inductively assume is a bicriterion approximation for Ri (more precisely, a map
B : Ri → Mi such that ∪p∈Ri (B(p), w(p)) = (Mi , ui ). However, to maintain polylogarithmicspace the algorithm pushes (Mi , ui ) to the beginning of the stream and restarts the bicriterion
construction. This means that the algorithm, at this point, restarts by viewing the stream as
(P, w − wi + ui ) (i.e. replacing (Ri , wi ) with (Mi , ui )). Continuing in this way, the algorithm
maintains a bicriterion (Mi+1 , ui+1 ) for (Ri+1 , wi+1 − wi + ui ) (which is also a bicriterion for
(Ri+1 , wi+1 ) by Theorem 8.2) until the next phase change is triggered.
Now we explain our modifications to [BMO+ 11] (see Algorithm 3). The first step is
that the bicriterion our algorithm builds must be determined “online” in the following sense:
upon receiving a point (x, w(x)), the value of B(x) must be determined (and never be altered)
before receiving the next point from the stream.
This generalization is necessary for the following reason. Suppose we connect an element
p to a center b1 . Later in the stream, we open a new center b2 that becomes the closest
center to p. However, using polylogarithmic space, we have already deleted b1 and/or p from
memory and the state of our algorithm is identical to the case where b1 remains the closest
center to p. Therefore the connections must be immutable, and this results in non-optimal
connections.
Definition 8.1. An online [α, β]-bicriterion is an algorithm that maintains a bicriterion
(B, t) over a stream X, and operates online in the following sense. Upon arrival of each
point p, it is immediately decided whether p is added to B, and then t(p) is determined. Both
of these decisions are permanent.
Upon receiving an update (p, w(p)) from the stream, the algorithm may call MOVE(p, p0 )
for the nearest p0 ∈ B to p. In the analysis we use the function B : P → B that maps
each point p to its immediate location after this initial move (either p itself, or p0 ). If future
moves are performed on p, this does not change the value B(p). B is not stored by the
algorithm due to space constraints; only the value of B(p) for the most recent point p is used
by the algorithm, and older values are used in the analysis only. We will show that B is a
(O(1), O(log n))-approximation that allows us to maintain a coreset over the stream.
33

We now state Algorithm 3. φ and γ are constants (dependent on ρ) used in the analysis
that are defined as in [BMO+ 11]. Each point has a flag that is either raised or lowered (given
by the F lag function). All points have their flag initially lowered, given on Line 1. A lowered
flag shows that the point is being read for the first time (being received from the stream),
and a raised flag shows that the point is being re-read by the algorithm (having been stored
in memory).
On Line 22, (Mi , ui ) is the weighted set (B, u) as it exists at the end of phase i. We
define the cost of MOVE(p, p0 ) to be w(p)D(p, p0 ). The value Ki is therefore the total cost of
all moves performed in phase i.
At a phase change, the set (B, u) is pushed onto the beginning of the stream with all
points having a raised flag. This means that during the next |B| iterations of the outer-loop
where a point (x, w(x)) is received we actually receive a point from memory. We continue
to process the stream after these |B| points are read.
The following theorem summarizes the guarantees of this algorithm. We note that, as
stated, the algorithm’s runtime of O(nk log n) is not correct - in fact the number of phases
may be arbitrarily large. However, using the same technique as detailed in Section 3.3
n
) while requiring O(k 2 log2 n)
of [BMO+ 11] the number of phases can be bounded to O( k log
n
time per phase.
Theorem 8.2 ([BMO+ 11]). Let Algorithm 3 process a stream (P, w) of at most n points.
Let Ri denote the prefix of the stream received before the end of phase i, and let Ki and
Li be their final values from the algorithm (which are never modified after phase i). With
probability at least 1 − n1 , the following statements all hold after processing each point:
P
ρφγ
1.
i Ki ≤ φ−ρ OPTk (Ri )
2. The total runtime is O(nk log n)
3. For every phase i, Li ≤ OPTk (Ri ) ≤ φLi ≤ Li+1
4. At the execution of Line 22, Mi consists of O(k log n) points
Part 1 of the preceding theorem, which bounds the total cost of all moves performed
by the algorithm, also serves as an upper bound on the Earth-Mover distance between the
stream (P, w) and the maintained set (B, u).
Definition 8.3 (Earth-Mover Distance). Let (A, wA ) and (B, wB ) be weighted sets in (X, D).
Morever, let them be of equal weight in the sense that Σa∈A w(a) = Σb∈B w(b). Define a
“movement” from (A, wA ) to (B, wB ) to be a weighted set (Y, v) in (X × X, D) such that
∪(a,b)∈Y (a, v(a, b)) = (A, wA ) and ∪(a,b)∈Y (b, v(a, b)) = (B, wB ). Then dEM (A, wA , B, wB ) is
the minimum of Σ(a,b)∈Y v((a, b))D(a, b) over all movements Y from (A, wA ) to (B, wB ).
Another way to view the preceding definition is in terms of probability distributions.
Over all joint distributions over A × B with marginal distributions (A, wA ) and (B, wB ), we
seek the minimum possible cost (as defined) for any such joint distribution.
34

Algorithm 3: Input: integer k, ρ-metric space (X, D), stream (P, w) of n weighted
points from (X, D). Output: B(x) after receiving each point x, and a weighted set
(Mi , ui ) after each phase i ≥ 1
1 F lag(x) = 0 for all x ∈ P
2 L1 ← minimum D(x, y) for any x, y in the first k distinct points
3 i ← 1
4 K1 ← 0
5 B ← ∅
6 u(x) ← 0 for all x
7 for each point (x, w(x)) received do
8
u(x) ← u(x) + w(x)
9
y ← argminy∈Q D(x, y) (break ties arbitrarily)
, 1}, otherwise I ← 0
10
I ← 1 with probability min{ Liw(x)D(x,y)
/k(1+log2 n)
11
if I then
12
if F lag(x) = 0 then
13
B(x) ← x
14
B ← B ∪ {x}
15
else
16
Ki ← Ki + w(x)D(x, y)
17
u(y) ← u(y) + u(x) /* Step 1 of MOVE(x, y)
*/
18
u(x) ← 0 /* Step 2 of MOVE(x, y)
*/
19
if F lag(x) = 0 then
20
B(x) ← y
21
if Ki > γLi or |B| > (γ − 1)(1 + log2 n)k then
22
(Mi , ui ) ← (B, u)
23
F lag(b) ← 1 for all b ∈ B
24
Push (B, u) onto the stream (P, w) before the next point to read
25
B←∅
26
u(x) ← 0 for all x
27
Li+1 ← φLi
28
q←0
29
Ki+1 ← 0
30
i←i+1

35

From now on we will write a weighted set A instead of (A, wA ) when the meaning is clear.
If we start with a set A0 and apply n operations of MOVE until it becomes the set An , we can
provide an upper-bound for dEM (A0 , An ) by summing dEM (Ai , Ai+1 ) for 0 ≤ i < n. This is
a direct application of the triangle inequality. And if Ai+1 is obtained from Ai by applying
MOVE(p, p0 ), then dEM (Ai , Ai+1 ) = D(p, p0 )w(p), the cost of this move.
The Earth-Mover distance is important for clustering problems for the following reason. For any weighted sets (P, w) and (B, u) and query C, |ḡ(P, w, C) − ḡ(B, u, C)| ≤
dEM (P, w, B, u). This is immediate from a repeated application of the triangle-inequality
(proofs are found in Theorem 2.3 of [GMM+ 03] as well as in [COP03, BMO+ 11, Guh09]).

Theorem 8.4. There exists an algorithm stores O(log 1 k log n) points and maintains for
every prefix (R, w0 ) of the stream (P, w): (1) a weighted set (M, u) such that dEM (M, u, R, w0 ) ≤
OPTk (R), and (2) a (O(1), O(log( 1 ) log n))-bicriterion B for (R, w0 ). Morever, this bicriterion B is computed online in the sense that B(x) is determined upon receiving x from the
stream.
Proof. The algorithm will consist of running Algorithm 3 and storing certain information
for the λ most recent phases (where λ depends on ).
Let i be the current phase number, and let (R, w0 ) be the points received so far. We remind
the reader that when the meaning is clear we suppress notation and write R for the weighted
ργ
)e. The prefix R and the set
set (R, w0 ), and likewise for (M, u). Define λ = 2 + dlogφ ( (φ−ρ)
M in the statement of the theorem will be Ri−λ and Mi−λ . To upper bound dEM (R, M ),
which by definition is the minimum cost of any movement from R to M, we note that one
movement
is the set of moves carried out by the algorithm through phase i − λ, whose cost is
Pi−λ
ρφγ
i−λ
K
.
Part
1 of Theorem 8.2 then shows that dEM (R, M ) ≤ Σj=1
OPTk (Ri−λ ).
Kj ≤ φ−ρ
j
j=1
1−λ
Combining Statements 3 and 4 of the theorem shows that OPTk (Ri−λ ) ≤ φ OPTk (R).
ρφγ 1−λ
Therefore dEM (R, M ) ≤ φ−ρ
φ OPTk (R) ≤ OPTk (R) as desired.
As for the second statement of the theorem, the algorithm defines the map B and we are
currently interested in the restriction of B to R \ Ri−λ . B maps to at most O(k log n) points
per phase - this is guaranteed by Statement 5 (a direct result the termination condition on
Line 17). Over the last λ phases, this then maps to O(λk log n) = O(( 1 )k log n) points.
Therefore we have that β = O(( 1 ) log n). And the value of α is immediate from Statement 1
of Theorem 8.2, since B incurs only a subset of the costs as Kj (the subset that comes from
the new portion of the stream Rj \ Rj−1 ).
As a final note, we do not store B(P ) in memory (it may be linear in space). For
algorithms in the following sections it is only required to know it’s value for the most recent
point received; previous values are used only in the analysis.

8.2

Maintaining a coreset over the stream

In the previous section, we presented an algorithm that maintains an (α, β)-approximation of
the stream (this is given by the function B : P → B). In this section we show how we can use
this approximation to carry out the coreset construction of Algorithm 2 on an insertion-only
36

stream. In the offline construction, a sample S of m points is taken from X according to a
distribution where point p sampled with probability depending on D(p, B(p)), |B|, and nB(p)
(the total weight of points connected to the center B(p), which is written as Σq∈Pb w(q) where
b = B(p)) - the specific formula for the probability written below (and comes from Line 3 of
Algorithm 2). All three of these quantities can easily be maintained over the stream using
Algorithm 4 since B is online (i.e. B(p) never changes).
P rob(p) =

w(p)
w(p)D(p, B(p))
P
P
+
2 q∈P w(q)D(q, B(q)) 2|B| q∈Pb w(q)

Upon receiving a point p, we assign r(p) a uniform random number in the interval (0, 1).
This is the threshold for keeping p in our sample S - we keep the point p if and only if
r(p) < P rob(p). For each point p, P rob(p) is non-increasing as the stream progresses (this is
immediate from the formula and from the fact that clusters never decrease in size since B is
online). Therefore after receiving point p, we update P rob(s) for each s ∈ S and delete any
such s that drop below their threshold: r(s) ≥ P rob(s). Once a point crosses the threshold,
it may be deleted since P rob(s) is non-increasing and so it will remain below the threshold
at all future times. In this way, the construction exactly matches the output as if the offline
Algorithm 2 had been used.
Algorithm 4: Input: integer k, ρ-metric space (X, D), stream P of points n from
(X, D), online bicriterion B, sample size m ≥ 1. Output: a coreset S
1 S ← ∅
2 for each point p that arrives do
3
Assign r(p) a uniform random in (0, 1)
4
S ← S ∪ {p}
5
Compute Prob(p) according to Line 3 of Algorithm 2
6
for each point s ∈ S do
7
Update Prob(s)
8
if Prob(s) ≤ r(s) then
9
Delete s from S
10
Update u(s) according to Line 6 of Algorithm 2
Algorithm 3 provides the function B and a weighted set Mi−λ . Beginning in phase i − λ,
we will begin running Algorithm 4. This outputs the sample S for phases i − λ + 1 until the
current phase i. The following theorem shows that Mi−λ ∪ Si is a coreset for the stream. Of
course, we need to do this construction for each of the λ = O(log( 1 )) most recent phases, so
the space gets multiplied by this factor.
We return to definition 6.2 of an r Log-Log Lipschitz function D̃. Given a metric space
(X, dist), the space (X, D) where D = D̃(dist) is a ρ-metric space for ρ = max{2r−1 , 1}. It
is well-known that most M -estimators can be recast as a Log-Log Lipschitz function for a
low constant value of r. For example, k-means has r = 2.

37

Theorem 8.5. There exists a single-pass streaming algorithm requiring O(−O(1) k log n(log n+
log k + log 1δ )) space that maintains a (k, )-coreset for the k-median clustering of an r-LogLog-Lipschitz function D̃ on a stream of at most n elements with probability at least 1 − δ.
Proof. Running Algorithm 2 (producing the Mi ) and Algorithm 3 (producing the Si ) in
parallel, we will show that Mi−λ ∪ Si is the desired coreset of the stream.
We already have by Theorem 8.4 that dEM (Mi−λ , Ri−λ ) ≤ OPTk (P ) ≤ ḡ(P, w, C) for
any set C of size k. Also, by Theorem 6.6 we have the Si is an (k, )-coreset for P \
Ri−λ . This is because (in the statement of Theorem 6.6) we are required to carry out the

1
d
log
t
+
log(βk)
+
log
where t = α/ψ.
construction of Algorithm 2 using m ≥ ck(t+β)
2
ε
δ
It was shown in Lemma 6.2 that ψ = (/r)r . We are using the (O(1), O(log( 1 ) log n))approximation B, and in [FL11] it is shown that d = O(log n) for k-median in a ρ-metric
space. So therefore the minimal possible value of m thatsatisfies the hypotheses of the
−r
n)
theorem is O( k( ε+log
log n log −r + log(k log n) + log 1δ ). Simplifying notation, this is
2
O(−O(1) k log n(log n + log k + log 1δ )).
We write COST(A, C) to denote ḡ(A, wA , C) to simplify notation; by A ∪ B we mean
(A ∪ B, wA + wB ). Taking the union, we get that |COST(P, C) − COST(Mi−λ ∪ Si , C)| ≤
|COST(Ri−λ , C) − COST(Mi−λ , C)| + |COST(P \ Ri−λ , C) − COST(Si , C)|. The first term is
upper-bounded by the Earth-Mover distance, and the second term is upper-bounded by 
since Si is a (k, )-coreset for P \ Ri−λ . So therefore Mi−λ ∪ Si is a 2-coreset for the stream
P , and the proof is complete after rescaling .
The previous theorem has, as special cases, streaming coresets for k-median, k-means,
Lp , Cauchy estimators, Tukey estimators, etc. This the first algorithm that does not use
merge-and-reduce for building coresets over streams for any of these problems. Moreover,
the constant in the exponent for  is small, for the example of k-median the  dependence is
−3 log(1/).

References
[AB99]

M. Anthony and P. L. Bartlett. Neural Network Learning: Theoretical Foundations. Cambridge University Press, Cambridge, England, 1999.

[AHPV04] P. K. Agarwal, S. Har-Peled, and K. R. Varadarajan. Approximating extent
measures of points. Journal of the ACM, 51(4):606–635, 2004.
[AMR+ 12] Marcel R. Ackermann, Marcus Märtens, Christoph Raupach, Kamil Swierkot,
Christiane Lammersen, and Christian Sohler. Streamkm++: A clustering algorithm for data streams. J. Exp. Algorithmics, 17:2.4:2.1–2.4:2.30, May 2012.
[BF15]

Artem Barger and Dan Feldman. k-means for streaming and distributed big
sparse data. SDM’16 and arXiv preprint arXiv:1511.08990, 2015.

38

[BLK15]

Olivier Bachem, Mario Lucic, and Andreas Krause. Coresets for nonparametric estimation-the case of dp-means. In Proceedings of the 32nd International
Conference on Machine Learning (ICML-15), pages 209–217, 2015.

[BMO+ 11] Vladimir Braverman, Adam Meyerson, Rafail Ostrovsky, Alan Roytman,
Michael Shindler, and Brian Tagiku. Streaming k-means on well-clusterable
data. In Proceedings of the Twenty-second Annual ACM-SIAM Symposium on
Discrete Algorithms, SODA ’11, pages 26–40. SIAM, 2011.
[BS80]

Jon Louis Bentley and James B Saxe. Decomposable searching problems i.
static-to-dynamic transformation. Journal of Algorithms, 1(4):301–358, 1980.

[CEM+ 15] Michael B Cohen, Sam Elder, Cameron Musco, Christopher Musco, and
Madalina Persu. Dimensionality reduction for k-means clustering and low rank
approximation. In Proceedings of the Forty-Seventh Annual ACM on Symposium
on Theory of Computing, pages 163–172. ACM, 2015.
[Che08]

Ke Chen. A constant factor approximation algorithm for k -median clustering
with outliers. In Shang-Hua Teng, editor, SODA, pages 826–835. SIAM, 2008.

[Che09a]

Ke Chen. On coresets for k-median and k-means clustering in metric and euclidean spaces and their applications. SIAM Journal on Computing, 39(3):923–
947, 2009.

[Che09b]

Ke Chen. On coresets for k-median and k-means clustering in metric and euclidean spaces and their applications. SIAM J. Comput., 39(3):923–947, August
2009.

[COP03]

Moses Charikar, Liadan O’Callaghan, and Rina Panigrahy. Better streaming
algorithms for clustering problems. In Proceedings of the Thirty-fifth Annual
ACM Symposium on Theory of Computing, STOC ’03, pages 30–39, New York,
NY, USA, 2003. ACM.

[DDH+ 07] Anirban Dasgupta, Petros Drineas, Boulos Harb, Ravi Kumar, and Michael W.
Mahoney.
Sampling algorithms and coresets for lp regression.
CoRR,
abs/0707.1714, 2007.
[FFK11]

Dan Feldman, Matthew Faulkner, and Andreas Krause. Scalable training of
mixture models via coresets. In Advances in Neural Information Processing
Systems, pages 2142–2150, 2011.

[FL11]

Dan Feldman and Michael Langberg. A unified framework for approximating
and clustering data. In Proceedings of the Forty-third Annual ACM Symposium
on Theory of Computing, STOC ’11, pages 569–578, New York, NY, USA, 2011.
ACM.

39

[FMS07]

D. Feldman, M. Monemizadeh, and C. Sohler. A PTAS for k-means clustering
based on weak coresets. In SoCG, 2007.

[FS05]

G. Frahling and C. Sohler. Coresets in dynamic geometric data streams. In Proc.
37th Annu. ACM Symp. on Theory of Computing (STOC), pages 209–217, 2005.

[FS12]

Dan Feldman and Leonard J Schulman. Data reduction for weighted and outlierresistant clustering. In Proceedings of the twenty-third annual ACM-SIAM symposium on Discrete Algorithms, pages 1343–1354. SIAM, 2012.

[FSS13]

Dan Feldman, Melanie Schmidt, and Christian Sohler. Turning big data into
tiny data: Constant-size coresets for k-means, pca and projective clustering. In
Proceedings of the Twenty-Fourth Annual ACM-SIAM Symposium on Discrete
Algorithms, pages 1434–1453. SIAM, 2013.

[GMM+ 03] Sudipto Guha, Adam Meyerson, Nina Mishra, Rajeev Motwani, and Liadan
O’Callaghan. Clustering data streams: Theory and practice. IEEE Trans. on
Knowl. and Data Eng., 15(3):515–528, March 2003.
[Guh09]

Sudipto Guha. Tight results for clustering and summarizing data streams. In
Proceedings of the 12th International Conference on Database Theory, ICDT ’09,
pages 268–275, New York, NY, USA, 2009. ACM.

[HHR11]

Frank Hampel, Christian Hennig, and Elvezio Ronchetti. A smoothing principle
for the huber and other location m-estimators. Computational Statistics & Data
Analysis, 55(1):324 – 337, 2011.

[HP11]

Sariel Har-Peled. Geometric approximation algorithms, volume 173. American
mathematical society Providence, 2011.

[HPK07]

S. Har-Peled and A. Kushal. Smaller coresets for k-median and k-means clustering. Discrete Comput. Geom., 37(1):3–19, 2007.

[HPM04]

S. Har-Peled and S. Mazumdar. On coresets for k-means and k-median clustering. In STOC, 2004.

[Hub81]

P. J. Huber. Robust statistics. 1981.

[LBK15]

Mario Lucic, Olivier Bachem, and Andreas Krause. Strong coresets for hard
and soft bregman clustering with applications to exponential family mixtures.
CoRR, abs/1508.05243, 2015.

[LLS00]

Yi Li, Philip M. Long, and Aravind Srinivasan. Improved bounds on the sample
complexity of learning. In SODA, pages 309–318, 2000.

[LLS01]

Y. Li, P. M. Long, and A. Srinivasan. Improved bounds on the sample complexity
of learning. Journal of Computer and System Sciences (JCSS), 62, 2001.
40

[LS10]

M. Langberg and L. J. Schulman. Universal ε approximators for integrals. To appear in proceedings of ACM-SIAM Symposium on Discrete Algorithms (SODA),
2010.

[Mat89]

Jiřı́ Matoušek. Construction of epsilon nets. In Proceedings of the fifth annual
symposium on Computational geometry, pages 1–10. ACM, 1989.

[PKB14]

Dimitris Papailiopoulos, Anastasios Kyrillidis, and Christos Boutsidis. Provable
deterministic leverage score sampling. In Proceedings of the 20th ACM SIGKDD
international conference on Knowledge discovery and data mining, pages 997–
1006. ACM, 2014.

[SV07]

Nariankadu D Shyamalkumar and Kasturi Varadarajan. Efficient subspace approximation algorithms. In SODA, volume 7, pages 532–540. Citeseer, 2007.

[Tyl87]

David E Tyler. A distribution-free m-estimator of multivariate scatter. The
Annals of Statistics, pages 234–251, 1987.

[VC71]

V. N. Vapnik and A. Y. Chervonenkis. On the uniform convergence of relative
frequencies of events to their probabilities. Theory Prob. Appl., 16:264–280,
1971.

[Zha11]

Z. Zhang.
M-estimators.
http://research. microsoft.com/enus/um/people/zhang/INRIA/ Publis/Tutorial-Estim/node20.html, [accessed
July 2011].

41

