# Experimental Evaluation of Fully Dynamic $k$-Means via Coresets

Monika Henzinger Institute of Science and Technology Austria (ISTA),
Klosterneuburg, Austria David Saulpic11footnotemark: 1 Leonhard Sidl
University of Vienna, Austria. Work done while at ISTA.

###### Abstract

For a set of points in $\mathbb{R}^{d}$, the Euclidean $k$-means problems
consists of finding $k$ centers such that the sum of distances squared from
each data point to its closest center is minimized. Coresets are one the main
tools developed recently to solve this problem in a big data context. They
allow to compress the initial dataset while preserving its structure: running
any algorithm on the coreset provides a guarantee almost equivalent to running
it on the full data. In this work, we study coresets in a fully-dynamic
setting: points are added and deleted with the goal to efficiently maintain a
coreset with which a $k$-means solution can be computed. Based on an algorithm
from Henzinger and Kale [ESA’20], we present an efficient and practical
implementation of a fully dynamic coreset algorithm, that improves the running
time by up to a factor of 20 compared to our non-optimized implementation of
the algorithm by Henzinger and Kale, without sacrificing more than 7% on the
quality of the $k$-means solution.

##  1 Introduction

One of the most fundamental tools of data analysis is clustering, where groups
of points that are similar, or ”close” in some metric space are identified.
This task can be formalised using the $k$-means problem: find $k$ centers,
such that the sum of distances squared from each point to its closest center
is minimized. For a given set of centers, this sum is called the _cost_ of the
solution. A solution with low cost can be used to represent the original data,
with only a small loss in precision, by replacing each point with its closest
center. A solution also results in a valid clustering of the instance, where a
cluster is defined as the set of all points assigned to the same center.

The $k$-means problem has numerous applications, especially for instances with
a very large number of input points [32, 2]. Various algorithms to solve the
problem with different emphasis – running in near-linear time, as [18], or
using minimal few memory, as [6] – have been described.

The recently most studied strategy to find a $k$-means solution on big data
problems is to use _coresets_ : for a parameter $\varepsilon>0$, a set C is an
$\varepsilon$-coreset of the input $P$ if evaluating the cost of any candidate
solution on $C$ gives the same cost as evaluating it on the full input, up to
a factor of $(1\pm\varepsilon)$.

This line of work developed coresets with a size _independent of the original
dataset_ , which can be computed equally fast as a solution to $k$-means[17,
30]. Additionally, coresets are particularly useful when the available memory
is restricted, such as in streaming or distributed settings. We add to the
versatility of coresets by focusing on maintaining a coreset in the _fully
dynamic_ model, where points can be inserted or deleted from the dataset at
every time step.

Algorithms specifically designed for the dynamic setting are beneficial in any
applications, e.g. social networks or web-search queries, where the dataset is
evolving over time, and one wishes to maintain a good clustering. The naive
alternative to using a dynamic algorithm is to run a static algorithm after
each update; however this can become inefficient, and the solution may even be
outdated as soon as it is computed. To cope with this, Henzinger and Kale [29]
presented an algorithm that can efficiently maintain an $\varepsilon$-coreset.
More precisely, the complexity to deal with any update is
$O(k^{2}d\operatorname*{polylog}(n))$, and the $\varepsilon$-coreset obtained
has optimal size (essentially $O(k\varepsilon^{-4})$).

Once a coreset can be maintained, it also becomes very efficient to update a
solution to $k$-means: For instance, running the popular $k$-means++ algorithm
of [1] would result in a $O(\log k)$ approximation, with a complexity at every
time step of $O(k^{2}d\operatorname*{polylog}(n))$. Alternatively, algorithms
that provide a constant factor approximation with the same running time can be
used (see e.g. [40, 33]).

##### Our Contribution.

In this work, we explore the practical characteristics of the algorithm
described in [29] for maintaining both a coreset and a $k$-means
solution.111https://git.ista.ac.at/lsidl/dynamiccoreset We evaluate the
practical relevance of several algorithmic details required for strict
theoretical guarantees. Furthermore, we propose and test several, partly
heuristic, optimizations that improve the theoretical running time of the
dynamic algorithm further. After calibrating all parameters on small-scale
datasets from [44], we test our optimized algorithm on large-scale real-world
datasets. Even though our algorithm is in theory versatile enough to be used
in any metric space, we focus on euclidean datasets, as they are most common.

We found that an implementation that closely follows the algorithm of [29]
already vastly outperforms our baseline, despite its heavy details required
for theoretical guarantees. Combined with our optimizations, this results in
an algorithm running more than $550$ times faster on real-world dataset than
recomputing from scratch – which is the only baseline algorithm we are aware
of.

For the incremental setting, where all operations are insertions, our
optimizations lead to an algorithm provably maintaining a coreset with running
time $O(kd\log(n))$ per insertion, where $n$ is the size of the dataset. In
the opposite decremental setting, with only deletions, our algorithm also runs
in time $O(kd\log n)$, but without a theoretical guarantee on its quality.
Combining those two guarantees, our algorithm maintains a coreset under both
insertion and deletion in amortized time $O(kd\operatorname*{polylog}n)$,
improving [29] by a factor $k$.

###  1.1 Related Work

##### On Fully-Dynamic Clustering

The study of the $k$-means problem dates back to Lloyd’s work [35] on
quantization of a continuous signal. Since then, a lot of theoretical research
was dedicated to understanding the cost function. In summary, the problem is
NP-hard, even when $d=2$ [36] or $k=2$ [20]. However, it is possible to
compute a $(1+\varepsilon)$-approximation [14] in low dimensional space, and a
$5.91$-approximation in general spaces [13]. An orthogonal line of work, more
focused on the practical aspect, attempts to find solutions that are as good
as possible while minimizing running time: the most celebrated result – which
had a large influence on the practical success of $k$-means – is the
$k$-means++ algorithm of [1], a very simple algorithm that computes a $O(\log
k)$-approximation in time $O(nkd)$.

The fully dynamic setting received lots of attention in recent years, with
works of [34, 6, 25] for $k$-means, but also [15, 4] for the related facility
location problem or [12, 3] for $k$-center. However, only few of those works
have been carefully implemented, and mostly for the $k$-center problem [11,
10]. For dynamic _graph_ algorithms there is, however, a large body of
empirical work, see [27] for a survey.

##### On Coresets

As mentioned above, the coreset paradigm is an efficient way to solve
clustering in big data, as they allow to “turn big data into tiny data” [24],
on which the clustering problem can be efficiently solved. Introduced by [28],
the research on $k$-means coresets is booming (see e.g. [22, 31, 19, 16] and
references therein), which culminated in optimal coresets for $k$-means
clustering of size $O(\varepsilon^{-2}k\min(\varepsilon^{-2},\sqrt{k}))$ [17,
30], that can be computed as fast as any constant-factor approximation
algorithm. Using the $k$-means++ algorithm, the running time is therefore
$O(nkd)$. This progress is spreading to related problems as well, such as
capacitated clustering [5], clustering of lines [38] or clustering with
missing values [7]. A more comprehensive introduction to the coreset
literature can be found in [21].

###  1.2 Organisation of the Paper

We first introduce the relevant definitions and standard algorithm. In Section
3, we present the algorithm from [29] and the optimization we introduced. We
present the result of our experiments in Section 4.2.

##  2 Preliminaries

###  2.1 Definitions

We call a weighted set a set $X\subset\mathbb{R}^{d}$ with an associated
weight function $w:X\rightarrow\mathbb{R}^{+}$.

A _solution_ to $k$-means is any $k$-tuple of points in $\mathbb{R}^{d}$. The
(Euclidean) $k$-means problem is defined as follows. Given a set of points $X$
with weights $w$ and an integer $k$, the goal is to find the solution
$S^{*}\subset\mathbb{R}^{d}$ that minimizes
$\textsc{cost}_{w}(S,X):=\sum_{x\in X}w(x)\min_{s\in
S}\operatorname*{dist}^{2}(x,s)$. A solution $S$ is an $\alpha$-approximation
if it satisfies
$\textsc{cost}_{w}(S,X)\leq\alpha\cdot\textsc{cost}_{w}(S^{*},X)$. An
$(\alpha,\beta)$-bicriteria approximation is a set of $\beta k$ many centers,
with cost at most $\alpha\textsc{cost}_{w}(S^{*},X)$.

######  Definition 2.1.

For $\varepsilon\geq 0$, an $\varepsilon$-coreset for $k$-means is a set
$C\subseteq U$ with weights $w^{\prime}$ such that, for any candidate solution
$S$,

$\textsc{cost}_{w^{\prime}}(S,C)=(1\pm\varepsilon)\textsc{cost}_{w}(S,X).$

In the dynamic model, the input changes and the algorithm has to efficiently
maintain a solution to $k$-means. More precisely, for any set
$X\subset\mathbb{R}^{d}$ with an associated weight function
$w:X\rightarrow\mathbb{R}^{+}$, the algorithm should support operations
$Insert(x)$ that inserts point $x$ to the current dataset.222Our algorithm is
actually more general and supports insertion of weighted points. Additionally,
$Delete(x)$ removes $x$ from $X$. After each operation, the algorithm should
return an $\varepsilon$-coreset. We note that one can then run any (static)
approximation algorithm to extract a good solution from the coreset.

###  2.2 Static Coresets

Schwiegelshohn and Sheikh-Omar [42] present a comprehensive description of
coreset algorithms, with an experimental evaluation. As mentioned previously,
$\varepsilon$-coresets with optimal size
$O(\varepsilon^{-2}k\min(\varepsilon^{2},k^{1/2}))$ can be constructed in time
$\tilde{O}(|X|kd)$ [17, 30].

The result of [42] is that the _Sensitivity Sampling_ algorithm from Feldman
and Langberg [23] performs best. This algorithm has a theoretical running time
of $\tilde{O}(|X|kd)$, and produces an $\varepsilon$-coreset of size
essentially $\frac{kd}{\varepsilon^{4}}$ many points [23]. One of the strength
of this algorithm is its simplicity: to compute an $\varepsilon$-coreset, it
first computes an $(O(1),O(1))$-bicriteria approximation, and then samples
points (essentially) proportionate to their cost in the bicriteria solution.
The initial bicriteria approximation can be computed using the $k$-means++
algorithm. This algorithm works in rounds, each round sampling one point
proportionate to the squared distance to the points already sampled. This is
guaranteed to be a $O(1)$-approximation after $2k$ iterations [43]. Following
this, one iteration of _Lloyd’s Algorithm_ is performed: each center is
replaced by the optimal center for its cluster, which is the center of mass of
the cluster. While this does not formally improve the quality of
approximation, is results in a significantly better solution in practice. Note
that, [23] showed that, when (roughly) $\frac{kd}{\varepsilon^{4}}$ points are
sampled using sensitivity sampling, they form an $\varepsilon$-coreset (with
each point weighted by the inverse of its sampling probability). We give the
pseudo-code of the algorithms in Appendix B

##  3 The Dynamic Algorithms

###  3.1 Description of the base algorithm.

A thorough description of the dynamic framework can be found in the original
paper [29]. We provide here a detailed overview (without proofs) in order to
explain our heuristic improvements. The algorithm relies on two fundamental
properties of coresets, that allow to compose coresets together:

######  Property 3.1.

Let $C_{1}$ be an $\varepsilon$-coreset for $X_{1}$, and $C_{2}$ be an
$\varepsilon$-coresets for $X_{2}$. Then, $C_{1}\cup C_{2}$ is an
$\varepsilon$-coreset of $X_{1}\cup X_{2}$.

######  Property 3.2.

Let $C_{1}$ be an $\varepsilon$-coreset for $X$, and $C_{2}$ be a
$\delta$-coresets for $C_{1}$. Then, $C_{2}$ is an
$(\varepsilon+\delta+\varepsilon\delta)$-coreset of $X$.

The high-level idea of the algorithm is to use a combination of those two
properties. To compute an $\varepsilon$-coreset of a set $X$, one can split
$X$ into two halves $X_{1},X_{2}$, and compute $\varepsilon/2$-coresets
$C_{1},C_{2}$ for those. By Property 3.1, the union of $C_{1}$ and $C_{2}$ is
an $\varepsilon/2$-coreset for $X$. To reduce the size further, let $C$ be an
$\varepsilon/2$-coreset of $C_{1}\cup C_{2}$: Property 3.2 ensures that $C$ is
an $\varepsilon$-coreset of $X$. The reason for this “two-level” scheme is
that $C_{1}\cup C_{2}$ has size much smaller than $X$, and, thus, maintaining
a coreset of it will be more efficient than maintaining it for $X$. More
specifically, when the dataset is updated by an insert or remove operation,
half of the work is already done: if the update happens to be in $X_{1}$, then
$C_{2}$ is still an $\varepsilon/2$-coreset for $X_{2}$ and does not need to
be recomputed. Leveraging this idea recursively allows to dynamically maintain
a coreset with essentially $O(\log n)$ applications of Property 3.2.

More formally, the dataset is hierarchically decomposed in a bottom-up fashion
as follows: the algorithm maintains a partition of the dataset $X$ into parts
of size $1$, that are recursively assembled in a binary tree structure. A tree
node $v$ _represents_ the set of points contained in the leaves of the subtree
rooted at $v$, and we denote by $p(v)$ the parent of $v$, and by
$\gamma_{1}(v)$, $\gamma_{2}(v)$ the children of $v$ in the tree.

The algorithm ensures that, after processing an update, the following
invariant holds: for each tree node $v$ at height $\ell$, there is a set
$C(v)$ that is a $\frac{\varepsilon\ell}{\lceil\log n\rceil}$-coreset of the
points represented by $v$. This ensures the correctness of the algorithm: for
the root $r$, the set $C(r)$ is an $\varepsilon$-coreset of the current
dataset. This invariant is maintained using the following operations for point
insertions and deletions. At first we assume that the number of points in the
dataset is about $n$.

##### Insert($x$)

: to insert the point $x$, the algorithm adds $x$ to the first empty leaf
$v_{0}$ of the tree. To restore the invariant the algorithm needs to update
every set $C(v)$, for all ancestor $v$ of $v_{0}$ in the tree. For this, let
$v_{1},...,v_{\log n}$ such that $v_{i}=p(v_{i-1})$ be the sequence of
ancestors of $v_{0}$ ordered by height. The sets $C(v)$ are updated in this
order: first, let $C(v_{0})$ be an $\varepsilon/\lceil\log n\rceil$-coreset
for $v_{0}$. Then, for $i=1$ to $\log n$, let $C_{1}$ and $C_{2}$ be the
$\frac{\varepsilon(i-1)}{\lceil\log n\rceil}$-coresets of $\gamma_{1}(v_{i})$
and $\gamma_{2}(v_{i})$: let $C(v_{i})$ be an $\varepsilon/\lceil\log
n\rceil$-coreset of $C_{1}\cup C_{2}$. Using Properties 3.1 and 3.2, this is
indeed an $\frac{\varepsilon i}{\lceil\log n\rceil}$-coreset of the points
represented by $v_{i}$, and therefore the invariant is satisfied at the end of
this procedure.

##### Delete($x$)

: the algorithm identifies (with e.g. a hashtable) in which leaf $x$ is
stored, and removes it. Afterwards, all coresets of parents of the leaf
involved are updated similarly as for an insertion, in order to maintain the
invariant.

Details of how to delete and insert leafs are given in Appendix A.

In case the number of points in the dataset is not roughly $n$ but varies a
lot, the work of the dynamic algorithm is partitioned into phases. Each phase
ends after the number of points in the dataset changed by $50\%$ compared to
the beginning of the phase. At the beginning of a phase, each set $C(v)$ is
recomputed in order to maintain the invariant regarding the core-set size – as
the value of $\lceil\log n\rceil$ changes between phases, the precision of the
coreset in each tree node also has to change. The recomputation time can be
amortized over the updates that caused the change in $\lceil\log n\rceil$. By
performing the recomputation “spread out” over subsequence updates, the time
complexity can even be turned into a worst-case bound (see [29]), but as this
is not the focus of our work, we did not implement this.

This algorithm uses a static coreset construction as a blackbox. As explained
in Section 2.2, the state-of-the-art is the algorithm from [23], with a
bicriteria solution computed with $k$-means++. We use this algorithm both for
the dynamic algorithm and as a baseline. We note that better results may be
achievable with a different coreset construction: however, we expect the
relative comparison between the different dynamic algorithms to stay alike.

###  3.2 Running time optimizations.

We propose several improvements to the implementation of the previous
algorithm, some are quite standard and some are new heuristics. Our first
optimization is to compress the lower levels of the tree: instead of each leaf
representing a single point, we introduce a parameter $s$, and make each leaf
represent between $s/2$ and $s$ points – except for one special leaf that may
contain fewer than $s/2$ points. All insertions are performed in the special
leaf, and whenever it reaches size $s$, it is turned into a “normal” leaf, and
a new special one is created. This requires to slightly change the insert and
delete procedures. Furthermore, instead of computing a precise
$\varepsilon/\log n$-coreset, we parameterize our experiments directly by the
size of the coreset $s$. Taking $s\approx\frac{kd\log n}{\varepsilon^{2}}$
ensures that we compute an $\varepsilon$-coreset with some good probability.
In practice much smaller coreset size work just as well, as our experiments
show.

Our main approach for optimization is to perform lazy updates in order to
group the recomputations together: we show how to maintain the coreset while
avoiding some recomputations. Our improvement for insertion provably maintains
a coreset, while the deletion case is only heuristic.

###  3.3 Optimizing Insertions.

For insertions, each tree node has an _epoch_ : an epoch starts (i) either
after $s$ points have been inserted into the set represented by the node, (ii)
when a point is removed from the node or (iii) when the coreset of the node is
recomputed for another reason (e.g. beginning of a new phase). If a new epoch
is started for one node, the same is done for all ancestors of this node.
During an epoch, the tree node maintains two coresets: one for the points
present at the beginning of the epoch, and one for newly inserted points.
Since there are at most $s$ inserted points, it is very easy to maintain the
latter: it merely consists of all inserted points, with their initial weights.
The former is constant, as a new epoch is started if a point is removed.
Property 3.1 ensures that the union of those two sets is a coreset for the
whole set represented by the tree node. Therefore, the algorithm stays
correct. Furthermore, at the start of a new epoch, a coreset of size $s$ is
recomputed for all points that are currently represented by the node; and the
coreset for newly inserted points is set to $\emptyset$.

Although this optimization has no effect in the worst-case scenario, it
significantly improves the running time in insertion-only streams, as shown in
the next lemma. This hints that this optimization is beneficial when several
points are added in a row.

######  Lemma 3.1.

In a stream of $n$ insertion, the algorithm described above has total running
time $\tilde{O}(nkd)$, and maintains a valid $\varepsilon$-coreset at each
time step.

###### Proof.

Running the static coreset algorithm on a set of size $O(s)$, to produce a
coreset with size $s$, has running time $O(skd)$. To recompute the coreset of
a node at the beginning of an epoch, the algorithm collects all coresets
maintained by the children of this node – there are $4$ of them, $2$ for each
children – and then computes a coreset of size $s$ from the union of those $4$
coresets. This takes time $O(skd)$, as explained above. Furthermore, any tree
node has at most $O(\log n)$ many ancestors in a tree. Therefore, if a new
epoch is started at a node, $O(\log n)$ coreset are recomputed, and the total
complexity is $O(skd\log n)$. Now, if every inserted point gives a budget
$O(kd\log n)$ to each of its $O(\log n)$ ancestors, then, when a new epoch
starts from one node, its budget is enough to pay for the total recomputation.
Therefore, the overall complexity after $n$ insertions is at most the total
budget, which is $\tilde{O}(nkd)$.

The set maintained by the algorithm is an $\varepsilon$-coreset: the proof
follows from Property 3.1 and the fact that, during a given epoch, each node
stores a $0$-coreset for the points inserted during that epoch. ∎

###  3.4 Optimizing Deletions

We propose the following heuristic. When $p$ is deleted, the algorithm first
checks if $p$ is part of the _final_ coreset. If this is not the case, the
point is marked for removal, but no further changes are made. Otherwise, the
algorithm removes $p$, and all other points currently marked for removal,
using the non-optimized algorithm. When recomputing the coreset of a node due
to an insertion, all points marked and represented by this node are deleted
(and unmarked). We note that, when the final coreset has size much smaller
than $n$, many points may be deleted before a recomputation is triggered – if
points are deleted randomly, one needs to delete $O(n/s)$ points before
hitting one that is in the final coreset. However, those deleted points still
influence the coreset computation (by being influential at lower levels of the
tree). In order to balance the running time improvement and the quality of the
coreset, we introduce a cut-off $\Delta$, and proceed to recomputation also
when $\Delta\cdot n$ of the total points are marked to be removed. We present
in Section 4.2.1 different experiments to optimized $\Delta$ in two particular
cases:

$\bullet$ if points are deleted in the same order as they are inserted – e.g.,
a sliding windows – then it is best to set $\Delta\approx s/n$. In theory,
this would lead to a speed-up of the order $s$: indeed, $s$ points inserted
consecutively are stored in the same tree leaf by our algorithm. Therefore,
instead of recomputing coresets for this leaf and its parents $s$ times, our
optimization recompute it only once.

$\bullet$ if points are deleted randomly, then they come from potentially very
different leaves and the previous argument does not apply. However, random
deletions have a weaker affect on the $k$-means cost, and it appears possible
to take a much larger $\Delta$, e.g., $\Delta=3\%$. This ensure that, in
average, the $k$-means cost changes by only $3\%$, which is the order of error
we expect from our algorithm. We show experimentally that such a large
$\Delta$ indeed yields to a large speed-up.

Those results hint how to pick $\Delta$ according to the expected evolution of
the dataset.

###  3.5 Flattening the Tree

Our algorithm is heavily reliant on a balanced binary-tree datastructure,
which has to be recomputed regularly. Similar tree-like datastructures are
commonly used in practice and can often be made more efficient by fixing the
number of levels of the tree. Such strategies may increase the work done at
each level, but can reduce overhead significantly. We tested different ways of
flattening the tree into a _shallow tree_.

We suggest to use a complete $g$-ary tree with a fixed height $h$. To find the
optimal degree $g$ of the internal nodes, we balance the work done in the
leafs with the work in the remaining nodes as follows: there are $g^{h}$ many
leaves, each of those containing $n/g^{h}$ many points. Computing a coreset on
a leaf costs therefore essentially $k\cdot n/g^{h}$. On a tree node of arity
$g$, the algorithm collects $g$ coresets of size $s$ and computes a coreset of
those: therefore, the work done is $k\cdot gs$. To balance this work, we want
$g$ such that $k\cdot n/g^{h}=k\cdot gs$, i.e., $g=\sqrt[h+1]{n/s}$.

This optimization requires a priori knowledge of the size of the dataset $n$:
to cope with it, the algorithm can work in phases, and recompute from scratch
every time the value of $n$ changes by a large factor. We simplified the
setting for our experiments and ensured the number of points stayed close to a
fixed value. Those experiments show that, the shallow tree may provide a
speed-up, however highly dependent on the size of the dataset and the
parameters $k$ and $s$: typically, for $h=2,k=50,s=90$, the algorithm is
faster for $n=20.000$ but slower for $n=60.000$; this trend is reversed for
$h=5$. Precise results are shown in Section D.3. It appears that the
improvement is uncertain, as it depends on the size of the dataset. Therefore,
we find it preferable to have an adaptive height, as in our original
algorithm. However, this experience shows the potential benefits of an a
priory knowledge on the data set size.

###  3.6 Space Requirement and Data Structure.

In the algorithms as described above, the total memory size is $2nd$: at each
node of the tree, $s$ points are stored (which takes memory $sd$), and a
binary tree with $n/s$ leaves contains $2n/s$ nodes in total. However, we use
the following approach to reduce the memory requirement to $nd$, i.e., improve
it by a factor 2. First, instead of storing the $d$-dimensional points in each
node, one can store pointers to those input points.333More precisely, we store
each input point in a hashtable with keys the identifier of each point, in
order to ensure fast deletions of the points. This works as long as coreset
points are input points. This may not always be the case, as centers from the
bicriteria approximation algorithm also are coreset points. One can modify the
used bicriteria algorithm to ensure that all those centers are input points by
removing the Lloyd step of the $k$-means++ algorithm. The drawback is that the
quality of the bicriteria approximation worsen in practice, which in turn may
worsen the quality of the coreset. The weight (which fits into a single
variable) still need to be stored explicitly in each tree node. With those
improvements,the memory requirement is only $nd+n$ ($nd$ for the original
points and $n$ for all the pointers). We tested this in Section 4.2.7.

##  4 Experiments

###  4.1 Experimental Setup

####  4.1.1 Performance metrics.

In all experiments, we measured (1) the distortion of the coreset, (2) the
quality of the $k$-means solution computed on the coreset, and (3) the running
time of maintaining both the coreset and a $k$-means solution. All experiments
were repeated five times and the average result is reported. For the running
time measurements, we timed the computation of the coreset itself and the
computation of the $k$-means solution with the C++ library Chrono. (We did not
include the time needed to evaluate the distortion and quality.) Note that we
use a logarithmic scale whenever we plot running time results.

##### Coreset distortion:

To evaluate how well a coreset $C$ represents the original datapoints of
dataset $X$, one can consider a set $\mathcal{S}$ of candidate solutions, and
evaluate the _distortion_ of the coreset on each of those solutions as
follows:

$\displaystyle
D_{C}=\max_{S\in\mathcal{S}}\left(\max\left(\frac{\textsc{cost}(S,X)}{\textsc{cost}(S,C)},\frac{\textsc{cost}(S,C)}{\textsc{cost}(S,X)}\right)\right)-1$

If $\mathcal{S}$ is the set of all possible solutions, this definition ensures
that $C$ is a $D_{C}$-coreset for $X$. However, it is impossible to enumerate
efficiently over all those solutions – and [42] showed the co-NP-hardness of
checking whether a given set is a coreset. To cope with this, [42] proposed to
consider $\mathcal{S}$ to be a single solution, computed via $k$-means++
algorithm on the coreset $C$. We slightly extend $\mathcal{S}$ by also adding
a solution computed by $k$-means++ on the full dataset $X$. Note that [42]
showed that adding solutions generated uniformly at random under some natural
distributions to $\mathcal{S}$ was pointless, as for those random solutions
the quality of the coreset is always very good.

#####  $k$-means quality:

As our coreset can be used to dynamically compute a solution $S_{C}$ to
$k$-means, we compare the cost of $S_{C}$ to a solution $S_{X}$ computed on
the full dataset as follows:

$Q_{\textsc{$k$-means}}=\textsc{cost}(S_{X},X)/\textsc{cost}(S_{C},X).$

To get the solution $S_{X}$, we use again the $k$-means++ algorithm.

####  4.1.2 Data sets and Update Sequences.

Since most data sets available for clustering are not intrinsically dynamic,
we followed different (standard) strategies to simulate a dynamical behaviour,
starting from an ordered data set,

  1. 1.

An _insertions-only_ data set can be generated by inserting each point after
another in the same order as they occur in the static dataset.

  2. 2.

A _sliding window_ of size $t$ can be used to include removal operations.
First $t$ points are inserted without any point removals. In the next phase
(called the _window_), operations alternate between removing the oldest point
in the dynamic dataset and inserting the next points of the static dataset.

  3. 3.

A _random window_ with insertion probability $\pi$ is created by inserting the
next point of the static dataset with a probability of $\pi$ or removing a
random existing point with a probability of $1-\pi$.

  4. 4.

A _snake window_ of size $t$ is created by chaining multiple random windows
together. First, a random window with insertion probability $0.9$ is
constructed, until the dynamic dataset contains $t$ points. Then, the
insertion probability is reduced to $\pi=0.1$ and operations are added, until
$0.2\cdot t$ points remain in the dynamic data set. This pattern can be
repeated, until all points of the static data set are inserted.

These different datasets allow us to highlight various aspects of the
algorithm, namely the impact of insertion (for insertion-only), the impact of
making deletions in order (for sliding windows) versus random deletions
(random window), and the impact of varying the size of the dataset (snake
window).

_Synthetic data set._ To optimize the parameters of the algorithms and show
details of the running time behaviour, we used a subset of a synthetic data
set with two dimensions, called _Birch_ [44], which consists of $100.000$
points having $100$ random Gaussian clusters with random size in
$\mathbb{R}^{2}$. For this dataset, we order the points as follows: All points
in the same cluster are placed consecutively in the order. Then we created and
tested both an insertions-only as well as a sliding window dataset of size
$20000$ with a $k$ of $10$. Due to the ordering and the fact that the number
of ground truth clusters in the data set is larger than $k$, the optimal
$k$-means solution changes over time. Finally, to test the algorithm in a more
challenging setting, we also tested a snake window of size $20.000$ with
$80.000$ operations applied in total. We call the latter _Birch-snake_.

_Real-world data sets._ With the parameters chosen on the synthetic data sets
we then tested our algorithms on larger real-world data sets from the UCI
Database [37]. We used the datasets (1) _Taxi_ [41], which stores the start
coordinates of taxi rides in two dimensions, (2) _Twitter_ [8, 9], which
contains location information of tweets in two dimensions, and (3) _Census_
[39] dataset, which has 68 dimensions. For all these data sets we removed
duplicate points and reordered the points randomly.

We constructed dynamic data sets for Twitter and Census by introducing a snake
window of size $t$, where $t\in\\{0.5\text{e}6,1\text{e}6,1.5\text{e}6\\}$. To
ensure that the size of the dataset stays roughly constant (and be able to
have an indication of the running time for a given data set size), we modified
the description above and kept the insertion probability at $0.1$ until
$0.95\cdot t$ points remain in the dataset. To avoid side effects due to the
construction of the dataset, we start the measurements only after the first
$t$ operations. We then measured the performance of the next $500.000$
operations.

For the data set Taxi, we created a dynamic dataset using a standard sliding
windows of size $t\in\\{0.5\text{e}6,1\text{e}6\\}$. Measurements were again
taken during $500.000$ operations after the first $t$ points had been
inserted.

See Table 2 in Appendix for a listing of all the update sequences we created.
The number of clusters $k$ was set to 10. Recall that theory recomments to
choose $s$ linear in $k\log n$. Thus, we use $s=5k$ for the small synthetic
dataset and $s=50k$ for the real-world data sets.

##### Baseline algorithms.

We are not aware of any other implementations of dynamic coreset algorithms
and, thus, we compare our algorithm with four simple baselines:

$\bullet$ static, that computes a coreset from scratch after each update using
Sensitivity Sampling in $\mathcal{O}(|X|kd)$

$\bullet$ random, that chooses after each update a uniform sample of the
dataset in $\mathcal{O}(s)$

$\bullet$ only k-means, which computes a $k$-means solution directly on the
data in $\mathcal{O}(|X|kd)$. Obviously, this does not allow to compare corset
distortion. The running time cannot be directly compared to other algorithms,
since no coreset is created by this algorithm. However, if the ultimate goal
is to compute a $k$-means solution, this baseline still gives a worthwile
comparison.

$\bullet$ ILP, which calculates a $k$-means solution (centers are required to
be part of the dataset) directly on the whole data set. We used the commercial
solver Gurobi [26]. This algorithm was at least 10.000 times slower, with a
cost close to the one computed by our algorithm. We defer the related
discussion to Appendix D.4

Every experiment was performed on an Intel(R) Core(TM) i5-1235U CPU (4.4 GHz)
with 16 GiB of RAM running on Ubuntu 22.04.2 LTS with Kernel
5.19.0-45-generic. All implementations are in C++, compiled using g++ version
11.3.0 with the optimization flag -O3.

###  4.2 Analysis of the algorithms and calibration of parameters on small
dataset

In this section, we analyze and compare the implemented algorithms using the
small dataset _Birch_.

Note that the running times of the coreset algorithms include the time to
compute a $k$-means solution on top of the coreset, while the running time of
only k-means only measures the time for computing a $k$-means solution.

(a) Insertion only dataset with a total of $20000$ operations

(b) Snake Window with $80000$ operations and size $20000$

Figure 1: Running time comparison of all coreset algorithms. The static data
set Birch with $k=10$ and $s=50$ was used. The dynamic data sets were created
as insert only and snake window as described in Section 4.1.2. The plot shows
the running time of the algorithm after each update to the dataset. The
running times are averaged over 20 consecutive operations.

####  4.2.1 Optimizing Deletions

We first determine the parameter $\Delta$ needed for the optimized deletion
algorithm from Section 3.4. Figure 7 presents the running time of the
optimized dynamic algorithm for different values of $\Delta$. For sliding
windows, i.e. when points are removed in the same order as inserted, the best
running time is achieved for $\Delta=2s/n$. The reason is explained in Section
3.4: points inserted consecutively are stored in the same leaf, leading to a
huge reduction in recomputation time, as recomputation happens only once
$\Delta\cdot n$ points have been deleted, instead of after each deletion.
444Taking $2s/n$ instead of $s/n$ allows to save a bit more, as the two leaves
concerned by the deletions are consecutive in the tree, and therefore share
lots of ancestor, for which only one recomputation is needed with
$\Delta=2s/n$, while $2$ are required for $\Delta=s/n$. As expected, the cost
of the $k$-means solution rises linearly with $\Delta$, increasing by a factor
$1.02$ for $\Delta=2s/n$ compared to $\Delta=0$. The plots are presented in
Figure 7 of the appendix.

Using a dataset where points are removed randomly, a much larger $\Delta$ is
required to obtain the best possible speedup. However, the cost is also less
affected, varying by less than $2\%$ for all tested $\Delta$. The reason is
that removing points uniformly at random from all inserted clusters allows a
$k$-means solution to stay accurate for more operations. Exact results are
show in the appendix in Figure 7(b). Thus, as the impact on the $k$-means cost
is small, we choose $\Delta=2s/n$ in the subsequent experiments where
deletions appear in order, and $\Delta=0.03$ otherwise.

####  4.2.2 Running time comparison for only insertions.

The behavior of the algorithms when only inserting points can be seen in
Figure 1(a). The experiments confirm the theoretical running times of the
different algorithms as described in Section 4.1.2. Both in theory and in the
experiments, the static coreset and calculating the $k$-means solution
directly has a running time that is linear with the input size, while the
dynamic algorithms have a logarithmic dependency.

####  4.2.3 Running time comparison for sliding window.

We evaluated the behavior of the algorithms on a sliding window. The running
time of all tested algorithms stayed constant in the window, except for the
optimized dynamic which pays every $2s=100$ deletions for recomputing, and
achieve an average speedup of $20$ compared to the base dynamic coreset. The
static coreset and only $k$-means were slower that the optimized dynamic
coreset by factors $164$ and $790$ respectively. The exact results are shown
in the appendix (Figure 6). We note, that the data points are removed in the
same order as they were inserted, resulting in optimal conditions for the
optimized dynamic algorithm.

(a) Running time comparison

(b) Distortion of the coresets

Figure 2: Impact of _k_. The Birch-snake dataset was used with a varying _k_.
The coreset size was set to $5k$ and time measurements are averages over all
update operations after the first $20.000$ points have been inserted. The
experiments were repeated five times and the median was used.

####  4.2.4 Running time comparison for Snake Window.

The snake-windows dataset represents a setting where both the size and points
in dynamic data change drastically. The running time of the different
algorithms can be seen in Figure 1(b). As the number of points in the dynamic
data set periodically decrease and increase in $4000$ operation intervals,
both the running time of the static coreset and $k$-means on the total data
set exhibit both a periodic pattern. This is not the case for the dynamic
coreset algorithms. The main reason for the lack of symmetry is that randomly
removing and inserting points leads to an imbalance in the tree, where each
leaf no longer stores the exact same number of points. This slightly reduces
the efficiency of both dynamic coreset algorithms. The spikes in the optimized
dynamic coreset algorithm indicate when recomputations happened.

####  4.2.5 Impact of Coreset Size $s$.

The effect of the coreset size on running time and $k$-means cost was tested
with $k=10$ and $s\in\\{40,...,280\\}$. We present the exact results in the
appendix (Figure 5).

The first observation is that the increase in running time with $s$ is faster
for the dynamic coresets than for the static coreset. This is due to the fact
that the dependency in $s$ of the static algorithm is $sk$ (to compute one
$k$-means solution on the coreset of size $s$), while the dynamic algorithm
needs to compute $\log n$ coreset of size $s$ – hence running time $sk\log n$.
With the used settings, the base dynamic coreset is already slower than the
static coreset at $s=350$. However, the number of points in the dynamic data
set used is relatively small ($20.000$ points). For experiments with larger
datasets, see Section 4.3.

Looking at the distortion of the different coresets, selecting points randomly
is beneficial when the size of the coreset $s$ grows and starts competing with
the other algorithms when $s=150$, but can be very poor for small values of
$s$.

####  4.2.6 Impact of number of clusters $k$.

The effect of $k$ on running time and coreset distortion is plotted in Figure
2. For each value of $k\in\\{10,...,280\\}$, a point corresponds to the
average running time (respectively cost) of each algorithm between operation
$20.000$ and $80.000$ of the dataset – to avoid side effects due to the
initial insertions.

We observe that the value of $k$ has slightly positive impact on the quality
of all the algorithms, but their relative order remains the same. However, for
running time, the non-optimized dynamic algorithm performs poorly for large
values of $k$: For $k$ roughly 100 or above it is slower than the static
algorithm. This is due to the fact that the dynamic implementation is quite
intricate, and its running time scales with $sk\approx k^{2}$: Thus it cannot
compete with the simple coreset construction for large values of $k$. However,
this effect is mitigated with our optimized implementation, which is faster
than the static algorithm for all values of $k$ we experimented with and also
has a better coreset quality than the non-optimized dynamic algorithm.

####  4.2.7 Reduction of memory overhead.

We also tested a variant of the algorithm without the Lloyd step in the
bicriteria algorithm, as described in Section 3.6. We use $k$ of 10 and $s=50$
with the Birch-snake dataset. The results are the following: as expected, the
memory footprint of the tree was drastically reduced. However, to process
deletions, the coreset is stored in a hashtable (unordered map), and, thus,
the memory reduction was only $8\%$. Furthermore, the running time was reduced
by $~{}11\%$, and the cost of the $k$-means solution increases by $5\%$. We
decided to privilege quality over memory in the remaining, and did not use
this variant for our experiments on large data sets.

Table 1: Results of experiments with large datasets. Speedup of the optimized dynamic coreset over each of the other algorithms, and also $k$-means quality ($Q_{k}$) and coreset distortion ($D_{C}$) for each algorithm. ODyn is the optimized dynamic coreset, Dyn the dynamic coreset, Stat the static coreset, Rand the random coreset and KM is only $k$-means. The datasets have size $t$ and are created with $k=10$, $s=50k$. Input | Twitter | Taxi | Census
---|---|---|---
t $[e6]$ | 0.5 | 1.0 | 1.5 | 0.5 | 1.0 | 0.5 | 1.5
Speed-up | Dyn | 11 | 12 | 12 | 20 | 23 | 2.5 | 2.59
Stat | 240 | 447 | 580 | 270 | 550 | 158 | 431
Rand | 0.039 | 0.051 | 0.049 | 0.037 | 0.054 | 0.02 | 0.02
KM | 850 | 1663 | 2100 | 1100 | 2000 | 2028 | 5717
$Q_{k}$ | ODyn | 0.96 | 0.96 | 0.95 | 0.83 | 0.85 | 0.96 | 0.95
Dyn | 0.95 | 0.96 | 0.95 | 0.83 | 0.84 | 0.95 | 0.95
Stat | 0.95 | 0.95 | 0.94 | 0.89 | 0.89 | 0.97 | 0.97
Rand | 0.91 | 0.92 | 0.92 | 0.46 | 0.16 | 0.96 | 0.96
$D_{C}$ | ODyn | 0.62 | 0.62 | 0.65 | 0.91 | 1.3 | 0.66 | 0.71
Dyn | 0.62 | 0.64 | 0.67 | 0.91 | 1.3 | 0.66 | 0.75
Stat | 0.11 | 0.10 | 0.10 | 0.37 | 0.57 | 0.09 | 0.09
Rand | 0.12 | 0.20 | 0.21 | 3.8 | 22 | 0.07 | 0.08

####  4.2.8 Shallow Tree.

We implemented shallow trees of different height $h$ as described in Section
3.5. To give maximum advantage to the shallow tree algorithm, we use a random
sliding window with insertion probability $\pi=0.5$ and size 20000, which
keeps the data size roughly constant: no restructuring is necessary and the
number of points in each leaf is predetermined. The Brich dataset with
$k=10,50$ and $s=50,90$ was used. To verify that the degree
$g=\sqrt[h+1]{n/s}$ predicted by theory is indeed correct, we compared the
running time when $h=1$ with different values for $g$ using the parameters
shown in Table 2. In Figure 8, we show that the observed minimum lies very
close to the theoretical optimum of $g=\sqrt[h+1]{n/s}=20$. Therefore, we used
the theoretical values for $g$ in all further experiments.

When using the optimal $g$, the shallow tree of height one could compete with
the base dynamic coreset. The speedup is however dependent on many factors
like the height of the tree, the size of the coreset, $k$ and $s$. We show the
speedup of shallow trees with different heights $h$ and optimal $g$ in Figure
9(b) in the appendix. As expected, there exists one optimal set of $g$ and $h$
for each combination of $n$, $k$, $s$ and changing either parameter would
require a complete restructuring of the tree and a recalculation of all
coresets. Especially when the size of the dataset varies over time, frequent
restructuring of the tree would be necessary: this would yield essentially our
main algorithm. Therefore, we did not use a shallow tree for our experiments
on large data sets

###  4.3 Experiments on large data sets

We present an evaluation of our optimized dynamic algorithm on real-world
dataset. We used the datasets described in Section 4.1.2: Twitter and Census
with a snake window and $\Delta=0.03$ (as determined above), and Taxi with a
sliding window and $\Delta=2s/n$. We chose $k=10$ and $k=25$ (the latter is
presented in Appendix D.1 and shows results for the Taxi and the Twitter data
set) and $s=50k$.

We show in Table 1 the speedup of the optimized dynamic algorithm compared to
the baselines, together with the respective coreset distortion and $k$-means
qualities.

For $k=10$ (resp. 25) the speedup of the optimized dynamic coreset improves
with larger input sizes and lies between 2.5 and 23 (resp. 10 and 19) compared
to the dynamic coreset, between 158 and 580 (resp. 94 and 367) compared to the
static coreset and between 850 and 5717 (resp. 541 and 1725) for only
$k$-means.

Both for $k=10$ and $k=25$ the performance results are similar to the previous
data set: For both dynamic coresets, the standard and the optimized ones, the
distortion is worse than for the static coreset algorithm. However, the
quality of the solution remains high, where the dynamic coresets only incur a
small quality penalty. _Thus, due to the large speedups that it achieves, the
optimized dynamic coreset algorithm is a clear improvement over the static
one_. The running time improvement over the only-$k$-means algorithm is even
larger, but the quality of the solution is somewhat worse.

The performance of the random coreset is very dependent on the structure of
the used dataset. The fact that the solution quality of the random coreset on
the Twitter and Census datasets are comparable to the results of only
$k$-means indicates that the clusters in those datasets have roughly the same
size and exhibit no special structure. This is, however, not the case for the
Taxi dataset, which indicates more structured clusters. _It is, thus, not
advisable to use the a random coreset without extensive a priory knowledge of
the dataset._

##  5 Conclusion

We empirically evaluated the theoretically best algorithm for maintaining a
coreset for $k$-means in a fully dynamic setting given in [29]. We show that
our optimized variant of the algorithm offers a speed-up of 2 to 23
(increasing with the data set size), without loss in the quality, compared to
the straightforward implementation of the dynamic algorithm in [29]. With our
algorithm we can maintain a coreset as well as a $k$-means solution with a
quality comparable to executing a _static_ algorithm at every time step, while
achieving a speedup of about two orders of magnitude. Note that this was the
state-of-the-art prior to our implementation.

We provide various strategies to optimize even further if some basic
statistics of the dataset are known. First one is to tune the co called
_deletion cutoff_ $\Delta$ from our optimized dynamic coreset. While a
$\Delta$ of $2s/n$ is sufficient for cases where points are deleted roughly in
the same order as they are inserted, larger values of up to $0.03$ are
beneficial when deletions occur randomly.

Second, a shallow tree can be used when the size of the data set stays roughly
constant over time. Together with a priory knowledge of the size of the data
set, the optimal arity of the shallow tree can be calculated and fixed over
the computations to achieve a significant speedup. Without precise knowledge
of these parameters, we recommend using the optimized dynamic algorithm
instead.

Third, we present a strategy for reducing memory overhead with only a slight
penalty in solution quality. This reduction in quality if however additionally
compensated by a lower running time.

A natural question following our work is to extend those experiments to
$k$-median clustering. Although theory predicts exactly the same behaviour, to
our knowledge it has not been confirmed in practice – the practical evaluation
of coresets from [42] has also been performed only with respect to the
$k$-means objective.

Another exciting question rises from our experiments: while both dynamic
algorithms generate worse coresets than the static algorithm, surprisingly,
the quality of the resulting $k$-means solution is great. The difference with
the static algorithm is at most $1\%$ for the dataset census (Table 1), even
when the distortion of the dynamic coreset is $0.72$ compared to $0.09$ for
the static coreset (as shown in Table 1). Understanding and explaining this
gap seems a promising avenue: coreset construction may be much better “when it
matters”.

## Acknowledgements.

This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(Grant agreement No. 101019564 “The Design of Modern Fully Dynamic Data
Structures (MoDynStruct)” and the Austrian Science Fund (FWF) project Z 422-N,
project “Static and Dynamic Hierarchical Graph Decompositions”, I 5982-N, and
project “Fast Algorithms for a Reactive Network Layer (ReactNet)”, P 33775-N,
with additional funding from the netidee SCIENCE Stiftung, 2020–2024.

D. Sauplic has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement No 101034413.

## References

  * [1] David Arthur and Sergei Vassilvitskii.  k-means++: the advantages of careful seeding.  In Nikhil Bansal, Kirk Pruhs, and Clifford Stein, editors, Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA 2007, New Orleans, Louisiana, USA, January 7-9, 2007, pages 1027–1035. SIAM, 2007.
  * [2] Bahman Bahmani, Benjamin Moseley, Andrea Vattani, Ravi Kumar, and Sergei Vassilvitskii.  Scalable k-means++.  Proc. VLDB Endow., 5(7):622–633, 2012.
  * [3] MohammadHossein Bateni, Hossein Esfandiari, Hendrik Fichtenberger, Monika Henzinger, Rajesh Jayaram, Vahab Mirrokni, and Andreas Wiese.  Optimal fully dynamic _k_ -center clustering for adaptive and oblivious adversaries.  In Nikhil Bansal and Viswanath Nagarajan, editors, Proceedings of the 2023 ACM-SIAM Symposium on Discrete Algorithms, SODA 2023, Florence, Italy, January 22-25, 2023, pages 2677–2727. SIAM, 2023.
  * [4] Sayan Bhattacharya, Silvio Lattanzi, and Nikos Parotsidis.  Efficient and stable fully dynamic facility location.  In NeurIPS, 2022.
  * [5] Vladimir Braverman, Vincent Cohen-Addad, Shaofeng H.-C. Jiang, Robert Krauthgamer, Chris Schwiegelshohn, Mads Bech Toftrup, and Xuan Wu.  The power of uniform sampling for coresets.  In 63rd IEEE Annual Symposium on Foundations of Computer Science, FOCS 2022, Denver, CO, USA, October 31 - November 3, 2022, pages 462–473. IEEE, 2022.
  * [6] Vladimir Braverman, Gereon Frahling, Harry Lang, Christian Sohler, and Lin F. Yang.  Clustering high dimensional dynamic data streams.  In Doina Precup and Yee Whye Teh, editors, Proceedings of the 34th International Conference on Machine Learning, ICML 2017, Sydney, NSW, Australia, 6-11 August 2017, volume 70 of Proceedings of Machine Learning Research, pages 576–585. PMLR, 2017.
  * [7] Vladimir Braverman, Shaofeng H.-C. Jiang, Robert Krauthgamer, and Xuan Wu.  Coresets for clustering with missing values.  In Marc’Aurelio Ranzato, Alina Beygelzimer, Yann N. Dauphin, Percy Liang, and Jennifer Wortman Vaughan, editors, Advances in Neural Information Processing Systems 34: Annual Conference on Neural Information Processing Systems 2021, NeurIPS 2021, December 6-14, 2021, virtual, pages 17360–17372, 2021.
  * [8] T-H. Hubert Chan, Arnaud Guerqin, and Mauro Sozio.  Fully dynamic k-center clustering.  In Proceedings of the 2018 World Wide Web Conference, WWW ’18, page 579–587, Republic and Canton of Geneva, CHE, 2018. International World Wide Web Conferences Steering Committee.
  * [9] T-H. Hubert Chan, Arnaud Guerqin, and Mauro Sozio.  Fully dynamic k-center clustering github, 2018.
  * [10] T.-H. Hubert Chan, Arnaud Guerquin, Shuguang Hu, and Mauro Sozio.  Fully dynamic $k$-center clustering with improved memory efficiency.  IEEE Trans. Knowl. Data Eng., 34(7):3255–3266, 2022.
  * [11] T.-H. Hubert Chan, Arnaud Guerquin, and Mauro Sozio.  Fully dynamic _k_ -center clustering.  In Pierre-Antoine Champin, Fabien Gandon, Mounia Lalmas, and Panagiotis G. Ipeirotis, editors, Proceedings of the 2018 World Wide Web Conference on World Wide Web, WWW 2018, Lyon, France, April 23-27, 2018, pages 579–587. ACM, 2018.
  * [12] T.-H. Hubert Chan, Silvio Lattanzi, Mauro Sozio, and Bo Wang.  Fully dynamic k-center clustering with outliers.  In Yong Zhang, Dongjing Miao, and Rolf H. Möhring, editors, Computing and Combinatorics - 28th International Conference, COCOON 2022, Shenzhen, China, October 22-24, 2022, Proceedings, volume 13595 of Lecture Notes in Computer Science, pages 150–161. Springer, 2022.
  * [13] Vincent Cohen-Addad, Hossein Esfandiari, Vahab S. Mirrokni, and Shyam Narayanan.  Improved approximations for euclidean _k_ -means and _k_ -median, via nested quasi-independent sets.  In Stefano Leonardi and Anupam Gupta, editors, STOC ’22: 54th Annual ACM SIGACT Symposium on Theory of Computing, Rome, Italy, June 20 \- 24, 2022, pages 1621–1628. ACM, 2022.
  * [14] Vincent Cohen-Addad, Andreas Emil Feldmann, and David Saulpic.  Near-linear time approximation schemes for clustering in doubling metrics.  In J. ACM, volume 68, 2021.
  * [15] Vincent Cohen-Addad, Niklas Hjuler, Nikos Parotsidis, David Saulpic, and Chris Schwiegelshohn.  Fully dynamic consistent facility location.  In Hanna M. Wallach, Hugo Larochelle, Alina Beygelzimer, Florence d’Alché-Buc, Emily B. Fox, and Roman Garnett, editors, Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada, pages 3250–3260, 2019.
  * [16] Vincent Cohen-Addad, Kasper Green Larsen, David Saulpic, and Chris Schwiegelshohn.  Towards optimal lower bounds for k-median and k-means coresets.  In Stefano Leonardi and Anupam Gupta, editors, STOC ’22: 54th Annual ACM SIGACT Symposium on Theory of Computing, Rome, Italy, June 20 \- 24, 2022, pages 1038–1051. ACM, 2022.
  * [17] Vincent Cohen-Addad, Kasper Green Larsen, David Saulpic, Chris Schwiegelshohn, and Omar Ali Sheikh-Omar.  Improved coresets for euclidean k-means.  In NeurIPS, 2022.
  * [18] Vincent Cohen-Addad, Silvio Lattanzi, Ashkan Norouzi-Fard, Christian Sohler, and Ola Svensson.  Fast and accurate $k$-means++ via rejection sampling.  In Hugo Larochelle, Marc’Aurelio Ranzato, Raia Hadsell, Maria-Florina Balcan, and Hsuan-Tien Lin, editors, Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual, 2020.
  * [19] Vincent Cohen-Addad, David Saulpic, and Chris Schwiegelshohn.  A new coreset framework for clustering.  In Samir Khuller and Virginia Vassilevska Williams, editors, STOC ’21: 53rd Annual ACM SIGACT Symposium on Theory of Computing, Virtual Event, Italy, June 21-25, 2021, pages 169–182. ACM, 2021.
  * [20] Sanjoy Dasgupta and Yoav Freund.  Random projection trees for vector quantization.  IEEE Transactions on Information Theory, 55(7):3229–3242, 2009\.
  * [21] Dan Feldman.  Introduction to core-sets: an updated survey, 2020.
  * [22] Dan Feldman and Michael Langberg.  A unified framework for approximating and clustering data.  In Lance Fortnow and Salil P. Vadhan, editors, Proceedings of the 43rd ACM Symposium on Theory of Computing, STOC 2011, San Jose, CA, USA, 6-8 June 2011, pages 569–578. ACM, 2011.
  * [23] 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, page 569–578, New York, NY, USA, 2011. Association for Computing Machinery.
  * [24] Dan Feldman, Melanie Schmidt, and Christian Sohler.  Turning big data into tiny data: Constant-size coresets for k-means, pca, and projective clustering.  SIAM J. Comput., 49(3):601–657, 2020.
  * [25] Hendrik Fichtenberger, Silvio Lattanzi, Ashkan Norouzi-Fard, and Ola Svensson.  Consistent k-clustering for general metrics.  In Dániel Marx, editor, Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms, SODA 2021, Virtual Conference, January 10 \- 13, 2021, pages 2660–2678. SIAM, 2021.
  * [26] Gurobi Optimization, LLC.  Gurobi Optimizer Reference Manual, 2023.
  * [27] Kathrin Hanauer, Monika Henzinger, and Christian Schulz.  Recent advances in fully dynamic graph algorithms - A quick reference guide.  ACM J. Exp. Algorithmics, 27:1.11:1–1.11:45, 2022.
  * [28] Sariel Har-Peled and Soham Mazumdar.  On coresets for k-means and k-median clustering.  In Proceedings of the 36th Annual ACM Symposium on Theory of Computing, Chicago, IL, USA, June 13-16, 2004, pages 291–300, 2004.
  * [29] Monika Henzinger and Sagar Kale.  Fully-Dynamic Coresets.  In Fabrizio Grandoni, Grzegorz Herman, and Peter Sanders, editors, 28th Annual European Symposium on Algorithms (ESA 2020), volume 173 of Leibniz International Proceedings in Informatics (LIPIcs), pages 57:1–57:21, Dagstuhl, Germany, 2020. Schloss Dagstuhl–Leibniz-Zentrum für Informatik.
  * [30] Lingxiao Huang, Jian Li, and Xuan Wu.  Towards optimal coreset construction for (k, z)-clustering: Breaking the quadratic dependency on k.  CoRR, abs/2211.11923, 2022.
  * [31] Lingxiao Huang and Nisheeth K. Vishnoi.  Coresets for clustering in euclidean spaces: importance sampling is nearly optimal.  In Konstantin Makarychev, Yury Makarychev, Madhur Tulsiani, Gautam Kamath, and Julia Chuzhoy, editors, Proccedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing, STOC 2020, Chicago, IL, USA, June 22-26, 2020, pages 1416–1429. ACM, 2020.
  * [32] Abiodun M Ikotun, Absalom E Ezugwu, Laith Abualigah, Belal Abuhaija, and Jia Heming.  K-means clustering algorithms: A comprehensive review, variants analysis, and advances in the era of big data.  Information Sciences, 2022.
  * [33] Silvio Lattanzi and Christian Sohler.  A better k-means++ algorithm via local search.  In International Conference on Machine Learning, pages 3662–3671. PMLR, 2019.
  * [34] Silvio Lattanzi and Sergei Vassilvitskii.  Consistent k-clustering.  In Doina Precup and Yee Whye Teh, editors, Proceedings of the 34th International Conference on Machine Learning, ICML 2017, Sydney, NSW, Australia, 6-11 August 2017, volume 70 of Proceedings of Machine Learning Research, pages 1975–1984. PMLR, 2017.
  * [35] Stuart Lloyd.  Least squares quantization in pcm.  IEEE transactions on information theory, 28(2):129–137, 1982.
  * [36] Meena Mahajan, Prajakta Nimbhorkar, and Kasturi R. Varadarajan.  The planar k-means problem is np-hard.  Theor. Comput. Sci., 442:13–21, 2012.
  * [37] Kolby Nottingham Markelle Kelly, Rachel Longjohn.  The uci machine learning repository.  https://archive.ics.uci.edu.
  * [38] Yair Marom and Dan Feldman.  k-means clustering of lines for big data.  In Hanna M. Wallach, Hugo Larochelle, Alina Beygelzimer, Florence d’Alché-Buc, Emily B. Fox, and Roman Garnett, editors, Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada, pages 12797–12806, 2019.
  * [39] Meek, Thiesson, and Heckerman.  US Census Data (1990).  UCI Machine Learning Repository.  DOI: https://doi.org/10.24432/C5VP42.
  * [40] Ramgopal R Mettu and C Greg Plaxton.  Optimal time bounds for approximate clustering.  Machine Learning, 56:35–60, 2004.
  * [41] Luis Moreira-Matias, Michel Ferreira, Joao Mendes-Moreira, J. J., and L. L.  Taxi Service Trajectory - Prediction Challenge, ECML PKDD 2015.  UCI Machine Learning Repository, 2015.  DOI: https://doi.org/10.24432/C55W25.
  * [42] Chris Schwiegelshohn and Omar Ali Sheikh-Omar.  An empirical evaluation of $k$-means coresets, 2022.
  * [43] Dennis Wei.  A constant-factor bi-criteria approximation guarantee for k-means++.  Advances in neural information processing systems, 29, 2016.
  * [44] T. Zhang, R. Ramakrishnan, and M. Livny.  Birch: A new data clustering algorithm and its applications.  Data Mining and Knowledge Discovery, 1(2):141–182, 1997.

##  Appendix A Illustrations of the Algorithms

This section gives an illustrative overview of the procedures to insert a leaf
(Figure 3) and remove a leaf (Figure 4) from the used tree structure.

Figure 3: Insertion of a new Node. Nodes are represented as circles and leafs
are labeled $a$ to $f$. New nodes are marked with an asterisk. To insert a new
leaf, a new child node (∗) is added to the leftmost leaf with the lowest depth
($c$) and then swapped with its parent node ($c$). The new leaf ($f^{*}$) is
then added as a sibling of $c$. Figure 4: Deletion of a Node $e$. Nodes are
represented as circles and leafs are labeled $a$ to $e$. To remove leaf $e$,
it is first swapped with the rightmost leaf of the tree with the largest depth
($b$). In the second step, the sibling of $e$ is swapped with its parent node
and made a new leaf by deleting all its children.

##  Appendix B Pseudo-codes

Here, we introduce the algorithms described in the main text as pseudo code.
The include the $k$-means++ algorithm with input weights (Algorithm 1) and the
algorithm used to select point for the coreset of Sensitivity Sampling
(Algorithm 2). We furthermore describe how to insert and remove points from
the dynamic coreset and its underlying tree structure (Algorithm 4, Algorithm
5 and helper function Algorithm 3).

Algorithm 1 $k$-means++$(X,w,n)$

1:A set $X$ with a weight function $w$, and an integer $k^{\prime}$

2:A set $S$ with $k^{\prime}$ centers

3:Sample a point $x$ randomly from $X$

4:Initialize $S\leftarrow x$

5:while $i<k^{\prime}$ do

6: Sample a point $x$ according to the distribution
$w(x)\textsc{cost}(x,S)/\textsc{cost}_{w}(X,S)$

7: Update $S\leftarrow S\cup\\{x\\}$

8:end while

9:Let $S^{\prime}\leftarrow\emptyset$

10:for Each cluster $C_{s}$ induced by $s\in S$ do

11: $S^{\prime}\leftarrow S^{\prime}\cup\\{\frac{1}{|C_{s}}\sum_{p\in
C_{s}}p\\}$

12:end for

13:return $S^{\prime}$

Algorithm 2 Selecting Coreset Points

1:A set $X$ with a weights function $w$, and integer $k$, Number of Coreset
Points to select $s$

2:A list of points $C$ and weights $w^{\prime}$, such that the multiset
defined by $C$ with weights $w^{\prime}$ is an $\varepsilon$-coreset.

3:$S\leftarrow\textsc{$k$-means{++}}(X,w,2k)$

4:For any point $x$, let $w(S_{x})$ be the total weight of its cluster in
solution $S$.

5:Define the sensitivity of $x$
$s(x)=w(x)\textsc{cost}(x,S)/\textsc{cost}_{w}(X,S)+\frac{1}{w(S_{x})}$

6:for $i<s$ do

7: Sample a point $x$ proportionate to $s(x)$

8: $C[i]\leftarrow x$,
$w^{\prime}[i]\leftarrow\frac{\sum_{x^{\prime}}s(x^{\prime})}{(s-k)s(x)}$

9:end for

10:for Each center $s\in S$ with cluster $S_{s}$: do

11: Add $s$ to $C$ with weight $(1+\varepsilon)w(S_{s})-w^{\prime}(S_{s})$

12:end for

13:return $C,w^{\prime}$

Algorithm 3 UpdateInnerCoreset($v,k,s$)

1:A tree node $v$, the cardinality of the internal coresets $s$ and the number
of centers $k$. We are also given an algorithm Coreset($X,w,k,s$), that, for
any set of points $X$ with weights $w$, returns a weighted coreset of size $s$
for $k$-means on $X$.

2:if $v$ is a leaf then

3: $C(v),w^{\prime}(v)\leftarrow$Coreset($C(v),w^{\prime}(v),k,s$)

4:else

5: $C_{init}\leftarrow C(c_{1}(v))\cup C(c_{2}(v))$

6: $w\leftarrow w(c_{1}(v))\cup w(C_{2}(v))$

7: if $|C_{init}|\leq s$ then

8: $C(v)\leftarrow C_{init}$, $w^{\prime}(v)\leftarrow w$

9: else

10: $C(v),w(v)\leftarrow$Coreset$(C_{init},w,k,s)$

11: end if

12:end if

13:UpdateInnerCoreset($p(v),k,s$)

Algorithm 4 Insert(x, w, s, k)

1:A point $x$ with weight $w$ to insert into the tree $t$. The cardinality of
the coresets $s$ and the number of centers $k$.

2:We store $i$, which points to a leaf, which contains less than $s$ points or
$NULL$ if no such leaf exists.

3:if $i==NULL$ then

4: Let $l$ be the leftmost leaf with the lowest depth in tree $t$

5: Let $a,b$ be new nodes

6: $c_{1}(l)\leftarrow a$, $c_{2}(l)\leftarrow b$

7: $p(a)\leftarrow l$, $p(b)\leftarrow l$

8: $C(a)\leftarrow C(l),w^{\prime}(a)\leftarrow w^{\prime}(l)$

9: $i\leftarrow b$

10: $UpdateInnerCoreset(l,k,s)$

11:end if

12:$C(i)\leftarrow C(i)\cup x$

13:$w^{\prime}(i)\leftarrow w^{\prime}(i)\cup w$

14:$UpdateInnerCoreset(p(i),k,s)$

15:if $|C(i)|\geq s$ then

16: $i\leftarrow NULL$

17:end if

Algorithm 5 Remove(x, s, k, v)

1:A point $x$ stored in leaf $v$ to remove from the tree $t$. The cardinality
of the coresets $s$ and the number of centers $k$.

2:We store $i$, which points to a leaf, which contains less than $s$ points or
$NULL$ if no such leaf exists.

3:$C(v)\leftarrow C(v)\setminus x$

4:$UpdateInnerCoreset(p(v),k,s)$

5:if $|C(v)|\leq s/2$ then

6: Let $r$ be the rightmost leaf with the highest depth in tree $t$ and $x$
its sibling

7: $C(i)\leftarrow C(i)\cup C(v)$, $w^{\prime}(i)\leftarrow w^{\prime}(i)\cup
w^{\prime}(v)$

8: $C(v)\leftarrow C(r)$, $w^{\prime}(v)\leftarrow w^{\prime}(r)$

9: $swap(s,p(x))$

10: Make $x$ a leaf

11: $UpdateInnerCoreset(x,k,s)$

12: $UpdateInnerCoreset(i,k,s)$

13: if $|C(i)|\geq s$ then

14: $i\leftarrow NULL$

15: end if

16:end if

##  Appendix C Experimental Parameters

Table 2: Overview of experimental parameters Experiment | Dataset | Updates | size | k | s | $\Delta$
---|---|---|---|---|---|---
Optimizing deletions | Birch | sliding window | 10000 | 10 | 50 | 0-0.14
Running time only inserts | Birch | only inserts | 20000 | 10 | 50 | –
Running time sliding window | Birch | sliding Window | 20000 | 10 | 50 | $2s/n$
Running time snake window | Birch | snake window | 20000 | 10 | 50 | 0.03
ILP | Birch | Only inserts | 300 | 10 | 20 | –
Impact of s | Birch | snake window | 20000 | 10 | 40-280 | 0.03
Impact of k | Birch | snake window | 20000 | 10-280 | $5k$ | 0.03
Shallow Tree impact of $g$ | Birch | sliding window | 20000 | 10 | 50 | –
Shallow Tree impact of $h$ | Birch | sliding window | 20000-80000 | 10,50 | 50,90 | –
Memory reduction | Birch | snake window | 20000 | 10 | 50 | –
Large Twitter | Twitter | snake window | 0.5e6-1.5e6 | 10 | 500 | 0.03
Large Taxi | Taxi | sliding window | 0.5e6-1.0e6 | 10 | 500 | 2s/n
Large Census | Census | snake window | 0.5e6-1.5e6 | 10 | 500 | 0.03

##  Appendix D Additional experimental Data

This section includes additional experimental data that is not part of the
main body.

###  D.1 Large experiments.

Here, we give the experimental results using large datasets and $k=25$.

Table 3: Results of experiments with large datasets. Speedup of the optimized dynamic coreset over each of the other algorithms, and also $k$-means quality ($Q_{k}$) and coreset distortion ($D_{C}$) for each algorithm. ODyn stands for optimized dynamic coreset, Dyn for dynamic coreset, Stat for static coreset, Rand for random coreset and KM for only $k$-means. The datasets have size $t$ and are extracted from the static datasets Twitter and Taxi with $k=25$, $s=50k$ as described in Section 4.1.2. Input | Twitter | Taxi
---|---|---
t $[e6]$ | 0.5 | 1.0 | 1.5 | 0.5 | 1.0
Speed-up | Dyn | 10 | 11 | 12 | 18 | 19
Stat | 106 | 229 | 367 | 94.0 | 176
Rand | 0.02 | 0.04 | 0.04 | 0.02 | 0.03
KM | 541 | 1090 | 1725 | 709 | 1275
$Q_{k}$ | ODyn | 0.94 | 0.94 | 0.95 | 0.88 | 0.81
Dyn | 0.95 | 0.95 | 0.94 | 0.88 | 0.81
Stat | 0.95 | 0.94 | 0.93 | 0.91 | 0.87
Rand | 0.87 | 0.89 | 0.86 | 0.16 | 0.05
$D_{C}$ | ODyn | 0.53 | 0.53 | 0.58 | 0.53 | 0.85
Dyn | 0.54 | 0.54 | 0.60 | 0.57 | 0.86
Stat | 0.10 | 0.11 | 0.11 | 0.23 | 0.39
Rand | 0.20 | 0.18 | 0.22 | 12.6 | 67.1

###  D.2 Small experiments.

Here, we show the experiments on small datasets used to optimizes the
parameters of the tested algorithms. Figure 5 contains the impact of the
parameter $s$ on the dynamic algorithms. We show the detailed running time of
the various algorithms in a sliding widow in Figure 6.

(a) Running time comparison

(b) Distortion of the coresets

Figure 5: Impact of _s_. The Birch-snake dataset with the parameters given in
table 2 was used with a varying _s_. Measurements are averaged over all update
operations inside of the window. The experiments were repeated five times and
the median was used. Figure 6: Running time comparison of all coreset
algorithms using a sliding window. The static data set Birch with $k=10$ and
$s=50$ was used. The plot shows the running time of the algorithm after each
update to the dataset.

(a) Running time of the sliding Window with ordered deletions

(b) Running time of the sliding Window with random deletions

(c) Cost of the $k$-means solution of the sliding window with ordered
deletions

(d) Cost of the $k$-means solution of the sliding window with random deletions

Figure 7: Impact of $\Delta$ on the remove optimization. The Birch dataset
with a sliding window of size $10000$ with the parameters given in table 2 was
used. Points were removed in the same order as inserted or in random order.
The cost of the optimal $k$-means++ solution on the coreset and running time
is averaged over all steps in the sliding window.

###  D.3 Shallow Tree.

Figure 7 contains details of experiments to determine the best values for the
remove cut-off $\Delta$ while Figure 8 and Figure 9(b) show the impact of the
two parameters $g$ and $h$ on the shallow tree.

The difference between $k=10$ and $k=50$ can be explained by the fact that our
algorithm needs to sometimes merge nodes in the tree (to keep each leaf large
enough), while the shallow tree does not. This feature could be incorporated
in our algorithm as well. However, this is very specific to the sliding window
case, where the size of the dataset stays fixed.

Figure 8: Impact of _g_ on the running time of the shallow tree of height one.
All measurements were taken in a sliding window of size $20000$ using
parameters given in table 2. The experiments were repeated five times and the
median was used.

(a) Impact of _h_ on the running time of the shallow tree of height one with
$k=10$ and $s=50$.

(b) Impact of _h_ on the running time of the shallow tree of height one with
$k=50$ and $s=90$.

Figure 9: Impact of _h_ on the running time of the shallow tree with optimal
$g$. All measurements were taken in a sliding window of different sizes using
parameters given in table 2.

###  D.4 ILP-based algorithm.

We tested the quality of this algorithm on a small dataset containing only
insertions (see Table 2), running the ILP from scratch at every iteration. As
its running time is by several orders of magnitude larger than any other
algorithm, we introduced a time limit of $100.000$ times the running time of
the dynamic coreset algorithm. This limit was already reached for very small
input sizes of around $300$ points. See Figure 10 for a comparison on the
results up to this data size. The solutions found by the ILP have slightly
better cost than the dynamic corset algorithm, but only when the ILP fully
converges before the time limit. For larger inputs the solver could not
produce solutions of a similar quality. ILP was therefore not further
investigated. Instead, we used the cost achieved as $k$-means++ algorithm as
the benchmark for quality.

(a) Running time comparison

(b) Cost of the $k$-means solutions on different coresets

Figure 10: Comparison with ILP. An insertion only data set from the first 300
points of the Birch-snake dataset was used with parameters given in table 2.
The cutoff for the running time of the ILP was set to $1e5$ times the running
time of the dynamic coreset. Measurements are averages over all update
operations. The experiments were repeated three times and the median was used.