# Min-Sum Clustering (with Outliers)

Sandip Banerjee The Hebrew University of Jerusalem<EMAIL_ADDRESS>Research supported in part by Yuval Rabani’s NSFC-ISF grant 2553-17. Rafail
Ostrovsky University of California, Los Angeles<EMAIL_ADDRESS>Supported
in part by DARPA under Cooperative Agreement No: HR0011-20-2-0025, NSF Grant
CNS-2001096, US-Israel BSF grant 2015782, Google Faculty Award, JP Morgan
Faculty Award, IBM Faculty Research Award, Xerox Faculty Research Award, OKAWA
Foundation Research Award, B. John Garrick Foundation Award, Teradata Research
Award, and Lockheed-Martin Corporation Research Award. The views and
conclusions contained herein are those of the authors and should not be
interpreted as necessarily representing the official policies, either
expressed or implied, of DARPA, the Department of Defense, or the U.S.
Government. The U.S. Government is authorized to reproduce and distribute
reprints for governmental purposes not withstanding any copyright annotation
therein. Yuval Rabani The Hebrew University of Jerusalem,
<EMAIL_ADDRESS>Research supported in part by NSFC-ISF grant 2553-17
and by NSF-BSF grant 2018687.

###### Abstract

We give a constant factor polynomial time pseudo-approximation algorithm for
min-sum clustering with or without outliers. The algorithm is allowed to
exclude an arbitrarily small constant fraction of the points. For instance, we
show how to compute a solution that clusters 98% of the input data points and
pays no more than a constant factor times the optimal solution that clusters
99% of the input data points. More generally, we give the following bicriteria
approximation: For any $\epsilon>0$, for any instance with $n$ input points
and for any positive integer $n^{\prime}\leq n$, we compute in polynomial time
a clustering of at least $(1-\epsilon)n^{\prime}$ points of cost at most a
constant factor greater than the optimal cost of clustering $n^{\prime}$
points. The approximation guarantee grows with $\frac{1}{\epsilon}$. Our
results apply to instances of points in real space endowed with squared
Euclidean distance, as well as to points in a metric space, where the number
of clusters, and also the dimension if relevant, is arbitrary (part of the
input, not an absolute constant).

##  1 Introduction

We consider min-sum $k$-clustering. This is the problem of partitioning an
input dataset of $n$ points into $k$ clusters with the objective of minimizing
the sum of intra-cluster pairwise distances. We consider primarily the
prevalent setting of instances of points in ${\mathbb{R}}^{d}$ endowed with a
distance function equal to the squared Euclidean distance (henceforth refered
to as the $\ell_{2}^{2}$ case). Our results apply also to the case of
instances of points endowed with an explicit metric (henceforth refered to as
the metric case). Note that we consider $k$ (and $d$, if relevant) to be part
of the input, rather than an absolute constant. In these and similar cases we
give polynomial time approximation algorithms that cluster all but a
negligible constant fraction of outliers at a cost that is at most a constant
factor larger than the optimum clustering. More specifically, for any
$\epsilon>0$, if the optimum we compete against is required to cluster any
number $n^{\prime}\leq n$ of points, our algorithm clusters at least
$(1-\epsilon)n^{\prime}$ points and at most $n^{\prime}$ points, and pays a
constant factor more than the optimum for $n^{\prime}$ points. The constant
depends on $\epsilon$.

Clustering in general is a fundamental question in unsupervised learning. The
question originated in the social sciences and now is widely applicable in
data analysis and machine learning, in areas including bioinformatics,
computer vision, pattern recognition, signal processing, fraud/spam/fake news
filtering, and market/population segmentation. Clustering is also a list of
fundamental discrete optimization problems in computational geometry that have
been studied for decades by theoreticians, in particular (but not exclusively)
as simple non-convex targets of machine learning. Some clustering problems,
notably centroid-based criteria such as $k$-means, have been studied
extensively. We currently have a fairly tight analysis of their complexity in
the worst case (e.g. [ADHP09, MNV09, ANFSW17, CC19]) and under a wide range of
restrictive conditions: low dimension (e.g., [HPK07, CAKM16, FRS16]) fixed $k$
(e.g. [KSS05, FMS07, Che09, FL11]), various notions of stability (e.g.
[ORSS06, ABS10, KK10, AMR11, CAS17]), restrictive models of computation (e.g.,
[BO07, AJM09, BMO+11, BFL16]), etc., as well as practically appealing
heuristics (e.g., Lloyd’s iteration, local search) and supportive theoretical
justification (e.g., some of the afore-mentioned papers and also [AV07,
JG12]).

Theoretical understanding of density-based clustering criteria, and in
particular min-sum clustering, is far less developed. There are clearly cases
in practice where, for instance, min-sum clustering coincides far better with
the intuitive clustering objective than standard centroid-based criteria. A
simple illustrative example is the case of separating two concentric dense
rings of points in the plane. Moreover, min-sum clustering satisfies
Kleinberg’s consistency axiom, whereas a fairly large class of centroid-based
criteria including $k$-means and $k$-median do not satisfy this axiom [Kle02,
ZBD09].

However, the state-of-the-art for computing min-sum clustering remains
inferior to alternatives. Min-sum $k$-clustering is NP-hard in the
$\ell_{2}^{2}$ case (e.g., using arguments from [ADHP09]), and also for the
metric case (see [FK01]), even for $k=2$. In the $\ell_{2}^{2}$ case, it can
be solved in polynomial time if both $d$ and $k$ are absolute constants
[IKI94]. In the metric case with arbitrary $k$, approximating min-sum
clustering to within a factor better than $1.415$ is NP-hard [GI03, CCL19]. If
$k$ is a fixed constant, the problem admits a PTAS, both in the $\ell_{2}^{2}$
case and in the metric case [FdlVKKR03]; see also [GH98, Ind99, Mat00, Sch00]
for previous work in this vein. If $k=o(\log n/\log\log n)$, then there is a
constant factor approximation algorithm for the $\ell_{2}^{2}$ case [CS07]. In
the metric case, assuming that $k=o(\log n/\log\log n)$ and the instance
satisfies a certain clusterability/stability condition, a partition close to
optimal can be computed in polynomial time [BBG09, BB09] (see also [VBR+11]
for some applications and experimental results in this vein). We note that
practical applications often require $k$ in the thousands, so the above
restrictions on $k$ are unrealistic in those cases.

In the worst case, and under no restrictions on the instance, the best known
approximation guarantee known is an $O(\log n)$ approximation algorithm
[BFSS19] for the metric case. This improves upon a slightly worse and much
earlier guarantee [BCR01]. In both papers, the factor is derived from
representing the input metric space approximately as a convex combination of
hierarchically separeted tree (HST) metrics [Bar96, Bar98, FRT03]. This incurs
logarithmic distortion, which is asymptotically tight in the worst case. In
HST metrics, min-sum clustering can be approximated to within a constant
factor. Thus, a fundamental challenge of the study of min-sum clustering is to
eliminate the gap between the hardness of approximation lower bound of $1.415$
and the approxiomation guarantee upper bound of $O(\log n)$. We show that a
constant factor approximation is possible, if one is willing to regard as
outliers a small fraction of the input dataset. For the $\ell_{2}^{2}$ case,
we are not aware of any previous non-trivial guarantee for $k\gg\log
n/\log\log n$.

Our results are derived using a reduction from min-sum clustering to a
centroid-based criterion with (soft) capacity constraints. This can be done
exactly in the $\ell_{2}^{2}$ case, and approximately in the metric case,
though to get polynomial time algorithms we use an approximation in both
cases. This reduction underlies also some of the above-mentioned previous
results on min-sum clustering. The outcome of this reduction is a $k$-median
or $k$-means problem with non-uniform capacities. If we are aiming for a
constant factor approximation then we can afford to violate the capacities by
a constant factor. There are nice results on approximating $k$-median with
non-uniform capacities, for instance [DL16]. Unfortunately, these results do
not seem applicable here, because their input is a metric space. The
reduction, even for the metric case, does not generate a metric instance of
capacitated $k$-median (the triangle inequality is violated unboundedly).
Nevertheless, we do draw some ideas from this literature.

Our min-sum clustering algorithm is based on the the well-trodden path of
using the primal-dual schema repeatedly to search for a good Lagrange
multiplier in lagrangian relaxation of the problem (see [JV01] for the origin
of this method). The dual program has a variable for every data point, and a
constraint for every possible cluster. The dual ascent process requires
detection of constraints that become tight. In our case, this is a non-trivial
problem, which we solve only approximately. As usual, the dual values are used
to “buy” the opening of the clusters that become tight, and we have to contend
with points contributing simultaneously to multiple clusters. This is done, as
usual, by creating a conflict graph among the tight clusters and choosing an
independent set in this graph. However, in our case there are unsual
complications. The connection cost is a distance (not a metric in the
$\ell_{2}^{2}$ case, but this is a minor concern) multiplied by the
cardinality of the cluster. If there is a conflict between a large cluster and
some small clusters, we have the following dilemma. If we open the large
cluster, the unclustered points in the small clusters may lack dual “money” to
connect to the large cluster; they can only afford the distance multiplied by
the cardinality of their small cluster. If, on the other hand, we open (some
or all of) the small clusters, assigning the unclustered points in the large
cluster to those small clusters might inflate their cardinality by a super-
constant factor, leaving all points with insufficient funds to connect to the
inflated clusters.

We resolve this dilemma as follows (using in part some ideas from [CR05]). We
open larger clusters first, so if a cluster is not opened, it is smaller than
the conflicting cluster that blocked it. Unclustered points are not assigned
to the blocking cluster, but rather aggregated around each blocking cluster to
form their own clusters of appropriate cardinality. We use approximate
cardinality, in scales which are powers of a constant $b$. As we require the
Lagrange multiplier preserving (LMP) property, we must have sufficient “funds”
to pay the opening costs in full (but can setttle for paying just a fraction
of the connection cost). This is possible if in a scale of, say, $b^{p}$ we
have, say, at least $b^{2+p}$ unclustered points in clusters of scaled
cardinality $b^{p}$ (each set of roughly this size can afford to open its own
cluster). If a blocking cluster is blocking fewer points in this scale, we
can’t afford to cluster them and must discard them as outliers. This is the
primary source of the excess outliers.

As usual, the search for a good Lagrange multiplier may end with two integer
primal solutions, one with fewer than $k$ clusters and one with more than $k$
clusters, whose convex combination is a feasible fractional bipoint solution
to the $k$-clustering problem. In our case, as we already may have to give up
on some outliers, we can simply output either the $<k$ solution or the $k$
largest clusters in the $>k$ solution. We point out that these extra outliers
can be avoided by using a more sophisticated “rounding” of the bipoint
solution, but given our loss in the primal-dual phase, it would not improve
meaningfully our guarantees.

The above description sums up the algorithm in the case that $n^{\prime}=n$.
Our result also extends to the case that the optimal solution is also allowed
to discard some outliers (we may have to discard some more). The main
additional issue in the case $n^{\prime}<n$ is that in the primal-dual phase
we may open a cluster that brings the number of clustered points from below
$n^{\prime}$ to above $n^{\prime}$. In this case, some points in this last
cluster need to be discarded, but then the remaining clustered points might
have insufficient “funds” to open the last cluster. If we have many clusters,
we can afford to eliminate the smallest cluster, declaring its points as
outliers, and use the dual values of the points in that cluster to pay for
opening the last one. If there is a small number of clusters, we may assume
that the primal-dual phase opened less than $k$ clusters (to ensure this
property, if $k$ is a small constant, we employ the known PTAS; thus we can
assume that $k$ is large). Our approach in this case draws from [AS16], where
a similar issue is addressed in the case of the sum-of-radii $k$-clustering
problem. Though the questions are quite different, we use a similar idea of
computing a (slightly) non-Lagrange multiplier preserving approximation to the
lagrangian relaxation. The LMP property is regularly used in the argument that
the bipoint solution is both feasible and cheap; the approach we adopt
requires an extra argument to bound the cost of a bipoint solution that
incorporates a non-LMP solution.

The rest of the paper is organized as follows. Section 2 introduces some basic
definitions and claims. Section 3 describes the algorithm. Section 4 analyzes
the algorithm. For conciseness, the paper presents the $\ell_{2}^{2}$ case.
Our main result is Theorem 8. The metric case is essentially identical, and is
briefly explained in Theorem 9. We note that we made no effort to optimize the
constant factor guarantees, throughout the paper.

##  2 Definitions and Preliminary Claims

Consider an instance of min-sum clustering that is defined by a set of points
$X\subset{\mathbb{R}}^{d}$ and a target number of clusters $k\in{\mathbb{N}}$.
Let $n=|X|$. The cost of a cluster $Y\subset X$ is

$\operatorname{cost}(Y)=\frac{1}{2}\cdot\sum_{x,y\in Y}\|x-y\|_{2}^{2}.$

The center of mass (or mean) of $Y$ is
$\operatorname{cm}(Y)=\frac{1}{|Y|}\cdot\sum_{x\in Y}x$. The following
proposition is a well-known fact (for instance, see [IKI94]).

######  Proposition 1.

The following assertions hold for every finite set $Y\subset{\mathbb{R}}^{d}$.

  1. 1.

The center of mass $\operatorname{cm}(Y)$ is the unique minimizer of
$\sum_{x\in Y}\|x-y\|_{2}^{2}$ over $y\in{\mathbb{R}}^{d}$.

  2. 2.

$\operatorname{cost}(Y)=|Y|\cdot\sum_{x\in
Y}\|x-\operatorname{cm}(Y)\|_{2}^{2}$.

A min-sum $k$-clustering of $X$ is a partition of $X$ into $k$ disjoint
subsets $X_{1},X_{2},\dots,X_{k}$ that minimizes over all possible
partititions

$\sum_{i=1}^{k}\operatorname{cost}(X_{i})=\sum_{i=1}^{k}|X_{i}|\cdot\sum_{x\in
X_{i}}\|x-\operatorname{cm}(X_{i})\|_{2}^{2}.$

In the version allowing outliers, we are given a target $n^{\prime}\leq n$ of
the number of points to cluster, and we are required that
$\left|\bigcup_{i=1}^{k}X_{i}\right|\geq n^{\prime}$. Clearly, the version
without outliers is a special case of the version with outliers with
$n^{\prime}=n$. Let $\operatorname{opt}(X,n^{\prime},k)$ denote the optimal
min-sum cost of clustering $n^{\prime}$ points in $X$ into $k$ clusters.
Formally, we can express the goal as a problem of optimizing an exponential
size integer program:

$\begin{array}[]{lll}\hbox{minimize}&\sum_{Y\subset
X}\operatorname{cost}(Y)\cdot z_{Y}&\\\ \hbox{s.t.}&\sum_{Y\ni
x}z_{Y}+w_{x}\geq 1&\forall x\in X\\\ &\sum_{Y\subset X}z_{Y}\leq k&\\\
&\sum_{x\in X}w_{x}\leq n-n^{\prime}&\\\ &z_{Y}\in\\{0,1\\}&\forall Y\subset
X\\\ &w_{x}\in\\{0,1\\}&\forall x\in X.\end{array}$ (1)

Fix $b\in{\mathbb{N}}$, $b>1$. For $i\in{\mathbb{N}}$, let
$\operatorname{floor}_{b}(i)=b^{\lfloor\log_{b}i\rfloor}$. For $Y\subset X$,
let $\operatorname{ctr}(Y)$ be a reference point that we set for now as
$\operatorname{ctr}(Y)=\operatorname{cm}(Y)$. Define

$\operatorname{cost}_{b}(Y)=\operatorname{floor}_{b}(|Y|)\cdot\sum_{y\in
Y}\|y-\operatorname{ctr}(Y)\|_{2}^{2}.$

In other words, (assuming $\operatorname{ctr}(Y)=\operatorname{cm}(Y)$) we
revise $\operatorname{cost}(Y)$ by rounding $|Y|$ down to the nearest power of
$b$. Thus,
$\frac{1}{b}\cdot\operatorname{cost}(Y)<\operatorname{cost}_{b}(Y)\leq\operatorname{cost}(Y)$.
We relax the integer program (1) as follows ($b$ to be determined later):

$\begin{array}[]{lll}\hbox{minimize}&\sum_{Y\subset
X}\operatorname{cost}_{b}(Y)\cdot z_{Y}&\\\ \hbox{s.t.}&\sum_{Y\ni
x}z_{Y}+w_{x}\geq 1&\forall x\in X\\\ &\sum_{Y\subset X}z_{Y}\leq k&\\\
&\sum_{x\in X}w_{x}\leq n-n^{\prime}&\\\ &z_{Y}\geq 0&\forall Y\subset X\\\
&w_{x}\geq 0&\forall x\in X.\end{array}$ (2)

Then, following a well-traveled path, we lagrangify the constraint on the
number of clusters to get the following lagrangian relaxation ($\lambda$
denotes the unknown Lagrange multiplier).

$\begin{array}[]{lll}\hbox{minimize}&\sum_{Y\subset
X}\operatorname{cost}_{b}(Y)\cdot z_{Y}+\lambda\cdot\left(\sum_{Y\subset
X}z_{Y}-k\right)&\\\ \hbox{s.t.}&\sum_{Y\ni x}z_{Y}+w_{x}\geq 1&\forall x\in
X\\\ &\sum_{x\in X}w_{x}\leq n-n^{\prime}&\\\ &z_{Y}\geq 0&\forall Y\subset
X\\\ &w_{x}\geq 0&\forall x\in X.\end{array}$ (3)

For fixed $\lambda$, this is a linear program, and its dual is:

$\begin{array}[]{lll}\hbox{maximize}&\sum_{x\in
X}\alpha_{x}-\gamma\cdot(n-n^{\prime})-\lambda\cdot k&\\\
\hbox{s.t.}&\sum_{x\in
Y}\alpha_{x}\leq\lambda+\operatorname{cost}_{b}(Y)&\forall Y\subset X\\\
&0\leq\alpha_{x}\leq\gamma&\forall x\in X.\end{array}$ (4)

Notice that the linear program (3) can be interpreted as a relaxation of the
“facility location” version of the problem, with $\lambda$-uniform cluster
opening costs.

######  Lemma 2.

For any $\lambda$, the optimal value of the linear program (4) is a lower
bound on the optimal value of the integer program (1).

###### Proof.

Consider any optimal solution $(z,w)$ to the integer program (1). Notice that
we may assume that $\sum_{Y\subset X}z_{Y}=k$, otherwise we can split some
clusters to get exactly $k$ of them. Splitting clusters cannot increase the
cost of the solution. This is also a feasible solution to the linear program
(3). Moreover, the Lagrange term $\lambda\cdot\left(\sum_{Y\subset
X}z_{Y}-k\right)$ zeroes out, and $\sum_{Y\subset
X}\operatorname{cost}_{b}(Y)\cdot z_{Y}\leq\sum_{Y\subset
X}\operatorname{cost}(Y)\cdot z_{Y}$. By weak duality, the value of any
feasible solution to the dual program (4) is a lower bound on the value of any
feasible solution to the linear program (3). ∎

An obvious issue with the dual program (4) is that the number of constraints
is exponential in $n$. We want to construct a dual solution by growing the
dual variables, however, it is not clear how to detect new tight dual
constraints without enumerating over the $\exp(n)$ number of constraints. We
now address this issue. First consider the following fact.

######  Proposition 3.

Let $Y$ be a finite set of points in ${\mathbb{R}}^{d}$. There exists $y\in Y$
such that $\sum_{x\in Y}\|x-y\|_{2}^{2}\leq 2\cdot\sum_{x\in
Y}\|x-\operatorname{cm}(Y)\|_{2}^{2}$. (We note that the factor of $2$ can be
improved to $1+\epsilon$, for any $\epsilon>0$, using the center of mass of
$O(1/\epsilon^{2})$ points in $Y$, e.g. [IKI94, FdlVKKR03].)

Proof sketch. Notice that for every $y\in{\mathbb{R}}^{d}$,

$\sum_{x\in Y}\|x-y\|_{2}^{2}\leq\sum_{x\in
Y}\|x-\operatorname{cm}(Y)\|_{2}^{2}+|Y|\cdot\|y-\operatorname{cm}(Y)\|_{2}^{2}$

(see, e.g. [ORSS06]). Thus, by picking $y\in Y$ that minimizes
$\|y-\operatorname{cm}(Y)\|_{2}^{2}$, the proposition follows. ∎

An immediate consequence of Proposition 3 is that ${\cal F}=X$ is a set of $n$
points in ${\mathbb{R}}^{d}$, such that for every $Y\subset X$ there exists a
point $c_{Y}\in{\cal F}$ such that

$\sum_{x\in Y}\|x-\operatorname{cm}(Y)\|_{2}^{2}\leq\sum_{x\in
Y}\|x-c_{Y}\|_{2}^{2}\leq 2\cdot\sum_{x\in
Y}\|x-\operatorname{cm}(Y)\|_{2}^{2}.$ (5)

(We can improve the factor of $2$ to any constant $1+\epsilon$ by increasing
the size of ${\cal F}$ to $n^{O(1/\epsilon^{2})}$.) Now, given ${\cal F}$, set
initially $\operatorname{ctr}(Y)=c_{Y}$ for all $Y\subset X$. Notice that this
puts $\operatorname{cost}_{b}(Y)=\operatorname{floor}_{b}(|Y|)\cdot\sum_{y\in
Y}\|y-c_{Y}\|_{2}^{2}$. We consider the following revised dual program.

$\begin{array}[]{lll}\hbox{maximize}&\sum_{x\in
X}\alpha_{x}-\gamma\cdot(n-n^{\prime})-\lambda\cdot k&\\\
\hbox{s.t.}&\sum_{x\in
Y}\alpha_{x}\leq\lambda+\operatorname{floor}_{b}(|Y|)\cdot\sum_{x\in
Y}\|x-y\|_{2}^{2}&\forall Y\subset X,\ \forall y\in Y\\\
&0\leq\alpha_{x}\leq\gamma&\forall x\in X.\end{array}$ (6)

######  Lemma 4.

For any $\lambda$, the optimal value of the linear program (6) is at most
twice the optimal value of the integer program (1).

###### Proof.

The dual of the linear program (6) is

$\begin{array}[]{ll}\min&\left\\{\sum_{Y\subset X}\sum_{y\in
Y}\operatorname{floor}_{b}(|Y|)\cdot\sum_{x\in Y}\|x-y\|_{2}^{2}\cdot
z_{Y,y}+\lambda\cdot\left(\sum_{Y\subset X}\sum_{y\in
Y}z_{Y,y}-k\right):\right.\\\ &\left.\forall x\in X,\ \sum_{Y\ni x}\sum_{y\in
Y}z_{Y,y}+w_{x}\geq 1\wedge\sum_{x\in X}w_{x}\leq n-n^{\prime}\wedge z,w\geq
0\right\\}\end{array}$ (7)

Consider an optimal clustering of any $n^{\prime}$ points in $X$ into $k$
disjoint clusters $Y_{1},Y_{2},\dots,Y_{k}$. For all $Y\subset X$, set
$z_{Y,y}$ to be the indicator that $Y$ is a cluster in this list and
$y=c_{Y}$. Also, for all $x\in X$ set $w_{x}$ to be the indicator that $x$ is
not clustered. Clearly, this is a feasible solution to the linear program (7),
so its value is an upper bound on the optimal value of the linear program (6).
The Lagrange term vanishes as there are exactly $k$ non-zero values $Z_{Y,y}$.
Thus, the upper bound is

$\sum_{j=1}^{k}\operatorname{floor}_{b}(|Y_{j}|)\cdot\sum_{x\in
Y_{j}}\|x-c_{Y_{j}}\|_{2}^{2}\leq
2\cdot\sum_{j=1}^{k}\operatorname{floor}_{b}(|Y_{j}|)\cdot\sum_{x\in
Y_{j}}\|x-\operatorname{cm}(Y_{j})\|_{2}^{2}\leq
2\cdot\sum_{j=1}^{k}\operatorname{cost}(Y_{j}),$

where the first inequality uses Equation (5). ∎

In the primal-dual procedure, there is an active set
$\operatorname{active}\subset X$ of points for which it is safe to raise the
dual variable $\alpha_{x}$ for all $x\in\operatorname{active}$. We need to
detect when a new dual constraint becomes tight and requires the removal of
the points that are involved from $\operatorname{active}$. This can be done in
polynomial time for the revised dual program (6) as follows. For every $y\in
X$ and for every $j\in\\{0,1,2,\dots,\log_{b}\operatorname{floor}_{b}(n)\\}$,
we check if there exists $Y\subset X$ that satisfies ($i$) $y\in Y$; ($ii$)
$Y\cap\operatorname{active}\neq\emptyset$; ($iii$)
$\log_{b}\operatorname{floor}_{b}(|Y|)=j$; ($iv$) $\sum_{x\in
Y}\alpha_{x}\geq\lambda+b^{j}\cdot\sum_{x\in Y}\|x-y\|_{2}^{2}$. In order to
do this, consider the set of points $C_{y,j}=\\{x\in X:\ \alpha_{x}\geq
b^{j}\cdot\|x-y\|_{2}^{2}\\}$, and sort $C_{y,j}$ by nonincreasing order of
$\alpha_{x}-b^{j}\cdot\|x-y\|_{2}^{2}$.

######  Lemma 5.

There exists a choice of $Y,y,j$ that satisfies ($i$)–($iv$) iff there exists
a choice of $y,j$ such that $|C_{y,j}|\geq b^{j}$ and
$C_{y,j}\cap\operatorname{active}\neq\emptyset$ and the first
$\min\\{|C_{y,j}|,b^{j+1}-1\\}$ points in the above order that include $y$ and
at least one point from $\operatorname{active}$ are a set that satisfies
($i$)–($iv$).

###### Proof.

Clearly the existence of $y,j$ such that $C_{y,j}$ has the listed properties
implies the existence of $Y,y,j$ that satisfy ($i$)–($iv$). As for the other
direction, consider $Y,y,j$ that satisfy ($i$)–($iv$). Clearly $y\in C_{y,j}$.
Suppose that there exists a point $x\in Y\setminus C_{y,j}$. Then, putting
$Y^{\prime}=Y\setminus\\{x\\}$ and
$j^{\prime}=\log_{b}\operatorname{floor}_{b}(|Y^{\prime}|)\leq j$, we have
that $Y^{\prime},y,j^{\prime}$ also satisfy ($i$)–($iv$). Thus, we may assume
that $Y\subset C_{y,j}$. Now, the choice in the lemma of a subset of $C_{y,j}$
maximizes $\sum_{x\in Y}\left(\alpha_{x}-b^{j}\cdot\|x-y\|_{2}^{2}\right)$
subject to the conditions ($i$)–($iii$). Thus, this subset also satisfies
($iv$). ∎

There are $O(n\log n)$ pairs $y,j$. Listing and sorting each $C_{y,j}$ takes
at most $O(n\log n)$ operations. Listing the candidate $Y\subset C_{y,j}$ and
checking it takes $O(|C_{y,j}|)$ operations. Thus, finding a new tight
constraint can be done in polynomial time. (Trivially, we can discretize the
increase of the dual variables and/or use binary search to find the increase
that causes a new constraint to become tight. As we’re dealing with squared
Euclidean distance, if the input consists of finite precision rational
numbers, then all computed values are finite precision rational numbers.)

##  3 The Algorithm

We now describe the following three-phase primal-dual algorithm (see Algorithm
1 on page 1) that can be used to solve the facility location version of min-
sum clustering. In addition to the pointset $X$, the cluster opening cost
$\lambda$, and the target number of points $n^{\prime}$, the algorithm also
gets a (sufficiently large, TBD) parameter $b$ that governs the excess number
of discarded outliers in its output. Throughout the algorithm, sets of points
$Y\subset X$ will maintain values $\operatorname{card}_{b}(Y)$ and
$\operatorname{ctr}(Y)$. Clearly, we cannot do this explicitly and efficiently
for every set $Y\subset X$. We use Lemma 5 and its consequences to implement
the operations that we need, without storing explicitly these values for more
than $n$ sets. This affects only the first phase of the algorithm. For $x\in
X$ and $Y\subset X$, we denote throughout the paper
$d(x,Y)=b^{\operatorname{card}_{b}(Y)}\cdot\|x-\operatorname{ctr}(Y)\|_{2}^{2}$.
This is interpreted according to the relevant values of
$\operatorname{card}_{b}(Y)$ and $\operatorname{ctr}(Y)$.

Phase 1 constructs a dual solution and collects candidate clusters. During
phase 1, a point $x$ is either active or inactive. Initially, for all $x\in
X$, we set $\alpha_{x}$ to $0$, and we set $x$ to be active. The set of
candidate clusters $\operatorname{preclusters}$ is empty. We raise all active
$x$ at a uniform rate, and pause to change the status of points and clusters
if one of the following events happens.

  * •

There exists an active $x\in X$ and a cluster $Y\in\operatorname{preclusters}$
such that $\alpha_{x}\geq d(x,Y)$. In this case, replace $Y$ by $Y\cup\\{x\\}$
in $\operatorname{preclusters}$. The new cluster in
$\operatorname{preclusters}$ inherits the $\operatorname{card}_{b}$ and
$\operatorname{ctr}$ values from $Y$. Also set $x$ to be inactive.

  * •

There exists $Y\subset X$ that contains an active point and $y\in Y$ such that
the dual constraint associated with the pair $Y,y$ is tight. Explicitly,

$\sum_{x\in Y}\alpha_{x}\geq\lambda+\operatorname{cost}_{b}(Y),$

where we set
$\operatorname{card}_{b}(Y)=\log_{b}\operatorname{floor}_{b}(|Y|)$ and
$\operatorname{ctr}(Y)=y$. In this case, add an inclusion-wise minimal such
$Y$ to $\operatorname{preclusters}$ and set all $x\in Y$ to be inactive (and
set $\operatorname{card}_{b}(Y)$ and $\operatorname{ctr}(Y)$ as stated above).

The first phase ends as soon as the number of active $x\in X$ drops to
$n-n^{\prime}$ or lower. If this number drops below $n-n^{\prime}$, we do not
add the last cluster $Y_{\operatorname{last}}$ to
$\operatorname{preclusters}$, but keep it separately. Note that each new tight
constraint causes at least one point $x\in X$ to become inactive, hence the
number of sets $Y$ that require storing explicitly their parameters
$\operatorname{card}_{b}(Y)$ and $\operatorname{ctr}(Y)$ is at most
$n^{\prime}\leq n$.

In phase 2, we trim the set of candidate clusters and assign points uniquely
to the clusters in the trimmed list, as follows. Note that we need the
parameters $\operatorname{card}_{b}$ and $\operatorname{ctr}$ only for
clusters for which these values were stored explicitly in phase 1. Define a
conflict graph on the clusters in $\operatorname{preclusters}$. Two clusters
$Y_{1},Y_{2}\in\operatorname{preclusters}$ are connected by an edge in the
conflict graph iff there exists $x\in Y_{1}\cap Y_{2}$ such that
$\alpha_{x}>\max\\{d(x,Y_{1}),d(x,Y_{2})\\}$. In other words, the edge
$\\{Y_{1},Y_{2}\\}$ indicates that there is $x\in Y_{1}\cap Y_{2}$ that
contributes to the opening cost $\lambda$ of both $Y_{1}$ and $Y_{2}$. Next,
take a lexicographically maximal independent set ${\cal I}$ in the conflict
graph, ordering $\operatorname{preclusters}$ by non-increasing order of
$\operatorname{card}_{b}(Y)$, breaking ties arbitrarily. We group the points
clustered in $\operatorname{preclusters}$ into meta-clusters of the form
$(Y,Y^{\prime})$, where $Y\in{\cal I}$ indicates the meta-cluster, and
$Y^{\prime}$ is a set of points. (Thus, the entire meta-cluster associated
with $Y$ is $\cup_{(Y,Y^{\prime})\in\operatorname{metaclusters}}Y^{\prime}$.)
In particular, for $Y\in{\cal I}$, we put
$(Y,Y)\in\operatorname{metaclusters}$. Any remaining points in
$\operatorname{preclusters}$ are added as follows. If $Y^{\prime}\not\in{\cal
I}$, then let $Y^{\prime\prime}$ be the set of remaining points in $\\{x\in
Y^{\prime}:\ \alpha_{x}=\max_{y\in Y^{\prime}}\alpha_{y}\\}$, and let
$Y\in{\cal I}$ be such that $Y$ precedes $Y^{\prime}$ in the order on
$\operatorname{preclusters}$ and $\\{Y,Y^{\prime}\\}$ is an edge. Add
$(Y,Y^{\prime\prime})$ to $\operatorname{metaclusters}$, with
$\operatorname{card}_{b}(Y^{\prime\prime})=\operatorname{card}_{b}(Y^{\prime})$
and $\operatorname{ctr}(Y^{\prime\prime})=\operatorname{ctr}(Y)$. Finally, if
fewer than $n^{\prime}$ points were thus assigned to meta-clusters, add
$(Y_{\operatorname{last}},Y)$ to $\operatorname{metaclusters}$, where $Y$ is a
set of previously unclustered points from $Y_{\operatorname{last}}$ of the
cardinality needed to complete the number of clustered point to $n^{\prime}$.
(Notice that at least $n^{\prime}$ points are clustered in
$\operatorname{preclusters}\cup\\{Y_{\operatorname{last}}\\}$, so this is
possible.)

Phase 3 determines the final output clustering of the points. For every meta-
cluster $(Y,\cdot)$ and for every integer $p\leq\operatorname{card}_{b}(Y)$,
let $n_{Y,p}$ denote the number of points $x\in X$ such that there exists
$(Y,Y^{\prime})\in\operatorname{metaclusters}$ with $Y^{\prime}\ni x$ and
$\operatorname{card}_{b}(Y^{\prime})=p$. We open clusters as follows. For
$p=\operatorname{card}_{b}(Y)$, we open
$\max\\{1,\left\lfloor\frac{n_{Y,p-2}+n_{Y,p-1}+n_{Y,p}}{b^{2+p}}\right\rfloor\\}$
clusters and assign all the points counted in $n_{Y,p-2},n_{Y,p-1},n_{Y,p}$ to
these clusters, as evenly as possible.

######  Lemma 6.

The number of points in each such cluster is at most $2b^{2+p}$, and if $Y\neq
Y_{\operatorname{last}}$ then this number is at least $b^{p}$.

###### Proof.

If we open one cluster, then clearly $n_{Y,p-2}+n_{Y,p-1}+n_{Y,p}<2b^{2+p}$.
If we open $s>1$ clusters, then we must have $sb^{2+p}\leq
n_{Y,p-2}+n_{Y,p-1}+n_{Y,p}<(s+1)b^{2+p}$. Thus, the number of points in each
cluster is between $b^{2+p}$ and $(1+1/s)b^{2+p}$. Clearly, if $Y\neq
Y_{\operatorname{last}}$, then $(Y,Y)\in\operatorname{metaclusters}$, and by
the definition of $p=\operatorname{card}_{b}(Y)$, it holds that $|Y|\geq
b^{p}$. ∎

For $p<\operatorname{card}_{b}(Y)-2$, we open
$\left\lfloor\frac{n_{Y,p}}{b^{2+p}}\right\rfloor$ clusters. If this number is
at least $1$, we assign all the points counted in $n_{Y,p}$ to these clusters,
as evenly as possible. If this number is $0$, we discard all the points
counted in $n_{Y,p}$ as outliers.

######  Lemma 7.

In this step, if no cluster is opened then the number of points that are
discarded is less than $b^{2+p}$, and otherwise the number of points in each
cluster is at least $b^{2+p}$ and less than $2b^{2+p}$.

###### Proof.

The assertion is trivial. ∎

We are now ready to define our min-sum $k$-clustering algorithm (see Algorithm
2 on page 2). If $k\leq\frac{4}{\epsilon}$, we can run a PTAS or a constant
factor approximation for fixed $k$ (for instance [FdlVKKR03, CS07]).111These
papers consider only the case without outliers. The PTAS in [FdlVKKR03]
enumerates over cluster sizes and approximate cluster centers, then computes
an optimal assignment of the data points to the approximate centers, given the
corresponding cluster sizes. Clearly, the algorithm can be adapted trivially
to handle the case with outliers by modifying the target sum of cluster sizes.
Otherwise, our algorithm follows the general schema of the lagrangian
relaxation method. Let $\delta>0$ be determined later. We run the procedure
PrimalDual on various values of $\lambda$, and if the smallest returned
cluster has at most $\frac{\epsilon}{3}\cdot n^{\prime}$ points, we remove
this cluster. Using binary search on the Lagrange multiplier $\lambda$, we
find two values $\lambda_{1}<\lambda_{2}$, with
$\lambda_{2}-\lambda_{1}<\delta$, that satisfy the following property. The
above process (running PrimalDual, then removing the smallest cluster if it’s
sufficiently small) returns $k_{1}>k$ clusters for $\lambda=\lambda_{1}$, and
$k_{2}\leq k$ clusters for $\lambda=\lambda_{2}$. If
$\frac{k-k_{2}}{k_{1}-k_{2}}\geq 1-\frac{\epsilon}{4}$, we output the $k$
largest clusters in the solution for $\lambda_{1}$, and otherwise we output
the solution for $\lambda_{2}$.

######  Theorem 8.

The excution of procedure MinSumClustering$(X,k,n^{\prime},\epsilon)$ computes
a clustering of $X^{\prime}\subset X$ into $k$ clusters such that
$|X^{\prime}|\in[(1-\epsilon)n^{\prime},n^{\prime}]$, and the total cost of
the clustering of $X^{\prime}$ is at most
$O\left(\frac{1}{\epsilon^{3}}\right)\cdot\operatorname{opt}(X,n^{\prime},k)$.
The time complexity of this computation is
$\operatorname{poly}(n,\log(1/\epsilon),\log\Delta)$, where $\Delta$ is the
ratio of largest to non-zero smallest $\|\cdot\|_{2}^{2}$ distance in $X$.

###### Proof.

The performance guarantee is an immediate consequence of Corollary 11 below.
The running time is a straightforward analysis of the code. ∎

######  Theorem 9.

The same claim applies to instances of points in a metric space
$(X,\operatorname{dist})$, with $\operatorname{dist}$ replacing
$\|\cdot\|_{2}^{2}$ in the code and in the claim.

Proof sketch. The $\|\cdot\|_{2}^{2}$ distance can be replaced by any metric
distance $\operatorname{dist}$ in all claims starting from Proposition 3. The
proofs sometime require minor changes. In particular, in Lemma 13, the factor
$\frac{1}{9}$ can be improved to $\frac{1}{3}$ on account of the triangle
inequality, and this improves all the other constants that depend on it. ∎

##  4 Proofs

In this section we analyze the min-sum $k$-clustering algorithm. The analysis
builds on the following guarantees of the primal-dual schema.

######  Theorem 10.

For every $\epsilon\in(0,1]$ there exists a constant $c=c_{\epsilon}$ such
that the following holds. Let $\operatorname{clusters}$ be the output of
procedure PrimalDual($X,\lambda,n^{\prime},b$), and let $\alpha$ be the dual
solution computed during the execution of this procedure. Set
$\gamma=\max_{x\in X}\alpha_{x}$. Then,

  1. 1.

$(\alpha,\gamma)$ is a feasible solution to the dual program (6).

  2. 2.

$\sum_{Y\in\operatorname{clusters}}|Y|\in[(1-\frac{\epsilon}{3})n^{\prime},n^{\prime}]$.

  3. 3.

$c\cdot\left(|\operatorname{clusters}|-1\right)\cdot\lambda+\sum_{Y\in\operatorname{clusters}}\operatorname{cost}(Y)\leq
c\cdot\sum_{Y\in\operatorname{clusters}}\sum_{x\in Y}\alpha_{x}$.

######  Corollary 11.

Let $\operatorname{clusters}$ be the output of procedure MinSumClustering.
Then, the following assertions hold:

  1. 1.

$|\operatorname{clusters}|\leq k$.

  2. 2.

$\sum_{Y\in\operatorname{clusters}}|Y|\in[(1-\epsilon)n^{\prime},n^{\prime}]$.

  3. 3.

$\sum_{Y\in\operatorname{clusters}}\operatorname{cost}(Y)\leq\frac{8(c+1)}{\epsilon}\cdot\operatorname{opt}(X,n^{\prime},k)$.

###### Proof.

The first assertion follows directly from the definition of the procedure.

For the second assertion, let $\lambda_{i},\operatorname{clusters}_{i}$ be the
values that determine the output of the procedure. By Theorem 10,
$\sum_{Y\in\operatorname{clusters}_{i}}|Y|\geq(1-\frac{\epsilon}{3})n^{\prime}$.
If $Y_{\min,\lambda_{i}}$ is removed from $\operatorname{clusters}_{i}$, then
$|Y_{\min,\lambda_{i}}|\leq\frac{\epsilon}{3}\cdot n^{\prime}$. Thus, if $i=2$
then clearly the assertion holds. If $i=1$, then $\rho_{1}\geq
1-\frac{\epsilon}{4}$. Therefore,

$k_{1}\leq\frac{1}{\rho_{1}}\cdot
k\leq\left(1+\frac{\epsilon}{4-\epsilon}\right)\cdot
k\leq\left(1+\frac{\epsilon}{3}\right)\cdot k.$

Thus, the procedure removes from the output at most a fraction of
$\frac{\epsilon}{3}$ of the clusters in $\operatorname{clusters}_{1}$. As the
removed clusters are the smallest, they contain at most
$\frac{\epsilon}{3}\cdot n^{\prime}$ points.

As for the third assertion, consider the two solutions
$\operatorname{clusters}_{1},\operatorname{clusters}_{2}$ that are used to
determine the procedure’s output. For $i=1,2,$, let
$\alpha_{i},\operatorname{preclusters}_{i}$ be the output of PrimalDualPhase1
during the computation of $\operatorname{clusters}_{i}$. Put
$\gamma_{i}=\max_{x\in X}\alpha_{i,x}$. Let

$(\alpha,\gamma)=\rho_{1}(\alpha_{1},\gamma_{1})+(1-\rho_{1})(\alpha_{2},\gamma_{2}).$

Clearly, $(\alpha,\gamma)$ is a feasible solution to the dual LP (6) with the
constant $\lambda=\rho_{1}\lambda_{1}+(1-\rho_{1})\lambda_{2}$. Notice that
there are exactly $n-n^{\prime}$ points that are not included in
$\operatorname{preclusters}_{i}$. Each point $x\in X$ which is not included in
$\operatorname{preclusters}_{i}$ has $\alpha_{i,x}=\gamma_{i}$. (Notice that
all the points that are excluded are active. This is true even for points that
are discarded from the last tight cluster that gets included in
$\operatorname{preclusters}_{i}$.) So, the value of the solution
$(\alpha,\gamma)$ is

$\displaystyle 2\cdot\operatorname{opt}(X,n^{\prime},k)$ $\displaystyle\geq$
$\displaystyle\sum_{x\in X}\alpha_{x}-(n-n^{\prime})\gamma-\lambda k$
$\displaystyle=$ $\displaystyle\rho_{1}\cdot\left(\sum_{x\in
X}\alpha_{1,x}-(n-n^{\prime})\gamma_{1}-\lambda
k_{1}\right)+(1-\rho_{1})\cdot\left(\sum_{x\in
X}\alpha_{2,x}-(n-n^{\prime})\gamma_{2}-\lambda k_{2}\right)$ $\displaystyle=$
$\displaystyle\rho_{1}\cdot\left(\sum_{Y\in\operatorname{preclusters}_{1}}\sum_{x\in
Y}\alpha_{1,x}-\lambda
k_{1}\right)+(1-\rho_{1})\cdot\left(\sum_{Y\in\operatorname{preclusters}_{2}}\sum_{x\in
Y}\alpha_{2,x}-\lambda k_{2}\right)$ $\displaystyle\geq$
$\displaystyle\rho_{1}\cdot\left(\sum_{Y\in\operatorname{preclusters}_{1}}\sum_{x\in
Y}\alpha_{1,x}-\lambda_{1}k_{1}\right)+(1-\rho_{1})\cdot\left(\sum_{Y\in\operatorname{preclusters}_{2}}\sum_{x\in
Y}\alpha_{2,x}-\lambda_{2}k_{2}\right)$
$\displaystyle-\delta\cdot(k_{1}+k_{2}),$

where the first inequality follows from Lemma 4, and the first equality uses
the fact that $k=\rho_{1}k_{1}+(1-\rho_{1})k_{2}$.

For $i=1,2$ consider the final value of $\operatorname{clusters}_{i}$. By
definition, $k_{i}$ is one less than the number of clusters returned from
procedure PrimalDual, so by Theorem 10,

$\sum_{Y\in\operatorname{clusters}_{i}}\operatorname{cost}(Y)\leq
c\cdot\left(\sum_{Y\in\operatorname{clusters}_{i}}\sum_{x\in
Y}\alpha_{x}-k_{i}\cdot\lambda_{i}\right).$

In particular, the right-hand side is non-negative. Notice that if
$\rho_{1}\geq 1-\frac{\epsilon}{4}$, then clearly
$\rho_{1}>\frac{\epsilon}{4}$ and the cost of the output clustering is at most
$\sum_{Y\in\operatorname{clusters}_{1}}\operatorname{cost}(Y)$. Similarly, if
$\rho_{1}<1-\frac{\epsilon}{4}$, then $1-\rho_{1}>\frac{\epsilon}{4}$ and the
cost of the output clustering is at most
$\sum_{Y\in\operatorname{clusters}_{2}}\operatorname{cost}(Y)$. Either way, we
get that the cost of the clustering is at most
$\frac{8c}{\epsilon}\cdot\operatorname{opt}(X,n^{\prime},k)+\frac{4\delta}{\epsilon}\cdot(k_{1}+k_{2})\leq\frac{8(c+1)}{\epsilon}\cdot\operatorname{opt}(X,n^{\prime},k)$.
∎

We now proceed to analyze the primal-dual algorithm and to prove Theorem 10.
The notation follows Algorithm 1.

######  Lemma 12.

At the end of executing procedure PrimalDualPhase1, for every
$Y\in\operatorname{preclusters}$ and for every $x\in Y$, we have that
$\alpha_{x}\geq d(x,Y)$.

###### Proof.

When $Y$ is added to $\operatorname{preclusters}$ then there exists $y\in Y$
and $j=\log_{b}\operatorname{floor}_{b}(|Y|)$ such that $Y\subset C_{y,j}$. We
set $\operatorname{card}_{b}(Y)=j$ and $\operatorname{ctr}(Y)=y$, so by the
definition of $C_{y,j}$ the lemma holds. If a point $x$ is later added to $Y$,
the condition for doing it is that $\alpha_{x}\geq d(x,Y)$. ∎

######  Lemma 13.

At the end of executing procedure PrimalDualPhase2, the following assertions
hold:

  1. 1.

$\sum_{(Y,Y^{\prime})\in\operatorname{metaclusters}}|Y^{\prime}|=n^{\prime}$.

  2. 2.

For every $(Y,Y^{\prime})\in\operatorname{metaclusters}$ and for every $x\in
Y^{\prime}$, we have that $\alpha_{x}\geq\frac{1}{9}\cdot d(x,Y^{\prime})$.

###### Proof.

The first assertion holds as the points that are clustered in
$\operatorname{metaclusters}$ are all the points that are clustered in
$\operatorname{preclusters}$ plus some of the points clustered in
$Y_{\operatorname{last}}$. The number of such points is at most $n^{\prime}$
without $Y_{\operatorname{last}}$, and at least $n^{\prime}$ with
$Y_{\operatorname{last}}$. The algorithm takes from $Y_{\operatorname{last}}$
exactly the number of points needed to complete the number in
$\operatorname{preclusters}$ to $n^{\prime}$.

For the second assertion, consider
$(Y,Y^{\prime})\in\operatorname{metaclusters}$ and $x\in Y^{\prime}$. If
$Y=Y^{\prime}$, then Lemma 12 guarantees the assertion. Otherwise, consider
$Y^{\prime\prime}\in\operatorname{preclusters}$ that caused $(Y,Y^{\prime})$
to be added to $\operatorname{metaclusters}$. In particuar, $Y^{\prime}\subset
Y^{\prime\prime}$,
$\operatorname{card}_{b}(Y^{\prime})=\operatorname{card}_{b}(Y^{\prime\prime})\leq\operatorname{card}_{b}(Y)$,
and there exists $z\in Y\cap Y^{\prime\prime}$ such that
$\alpha_{z}>\max\\{d(z,Y),d(z,Y^{\prime\prime})\\}$. By the choice of
$Y^{\prime}$ in the algorithm, $\alpha_{x}=\max_{x^{\prime}\in
Y^{\prime\prime}}\alpha_{x^{\prime}}$, so it must be that
$\alpha_{z}\leq\alpha_{x}$. Notice that by Lemma 12,

$\alpha_{x}\geq
d(x,Y^{\prime\prime})=b^{\operatorname{card}_{b}(Y^{\prime})}\cdot\|x-\operatorname{ctr}(Y^{\prime\prime})\|_{2}^{2}.$

Also

$\displaystyle\alpha_{z}$ $\displaystyle>$
$\displaystyle\max\\{b^{\operatorname{card}_{b}(Y)}\cdot\|z-\operatorname{ctr}(Y)\|_{2}^{2},b^{\operatorname{card}_{b}(Y^{\prime\prime})}\cdot\|z-\operatorname{ctr}(Y^{\prime\prime})\|_{2}^{2}\\}$
$\displaystyle\geq$ $\displaystyle
b^{\operatorname{card}_{b}(Y^{\prime})}\cdot\max\\{\|z-\operatorname{ctr}(Y)\|_{2}^{2},\|z-\operatorname{ctr}(Y^{\prime\prime})\|_{2}^{2}\\}.$

Thus,

$\displaystyle d(x,Y^{\prime})$ $\displaystyle=$ $\displaystyle
b^{\operatorname{card}_{b}(Y^{\prime})}\cdot\|x-\operatorname{ctr}(Y)\|_{2}^{2}$
$\displaystyle\leq$ $\displaystyle
b^{\operatorname{card}_{b}(Y^{\prime})}\cdot\left(\|x-\operatorname{ctr}(Y^{\prime\prime})\|_{2}+\|z-\operatorname{ctr}(Y^{\prime\prime})\|_{2}+\|z-\operatorname{ctr}(Y)\|_{2}\right)^{2}$
$\displaystyle\leq$ $\displaystyle 9\cdot\alpha_{x}.$

This completes the proof. ∎

Let $\operatorname{clusters}$ be the output of procedure PrimalDualPhase3.
Recall that every cluster $Z\in\operatorname{clusters}$ is derived in some
iteration indexed by $(Y,\cdot)\in\operatorname{metaclusters}$. It holds that
either $Z\subset Y_{\max}$ or $Z\subset Y_{p}$ for some
$p<\operatorname{card}_{b}(Y)-2$. Notice that in the latter case, $Y\neq
Y_{\operatorname{last}}$. We will set implicitly $\operatorname{card}_{b}(Z)$
as follows.

$\operatorname{card}_{b}(Z)=\left\\{\begin{array}[]{ll}\operatorname{card}_{b}(Y)&Z\subset
Y_{\max},\\\ p&Z\subset Y_{p},\
p<\operatorname{card}_{b}(Y)-2.\end{array}\right.$

We will also set implicitly $\operatorname{ctr}(Z)=\operatorname{ctr}(Y)$.

######  Lemma 14.

If $Z\subset Y_{p}$ for some $p<\operatorname{card}_{b}(Y)-2$, then
$\sum_{x\in Z}\alpha_{x}\geq b\cdot\lambda$. The same is true if $Z\subset
Y_{\max}$, $Y\neq Y_{\operatorname{last}}$, and $|Y_{\max}|\geq
b^{2+\operatorname{card}_{b}(Y)}$.

###### Proof.

Consider $x\in Z\subset Y_{p}$, $p<\operatorname{card}_{b}(Y)-2$. There is a
pair $(Y,Y^{\prime})\in\operatorname{metaclusters}$ such that
$p=\operatorname{card}_{b}(Y^{\prime})<\operatorname{card}_{b}(Y)-2$, and
$x\in Y^{\prime}$. Moreover, there is
$Y^{\prime\prime}\in\operatorname{preclusters}$ such that $Y^{\prime}\subset
Y^{\prime\prime}$ and $\operatorname{card}_{b}(Y^{\prime\prime})=p$ and
$\alpha_{x}=\max_{y\in Y^{\prime\prime}}\alpha_{y}$. By the definition of
$\operatorname{card}_{b}$, we have that $|Y^{\prime\prime}|<b^{1+p}$.
Therefore, $\alpha_{x}>\frac{\lambda}{b^{1+p}}$. By Lemma 7, $|Z|\geq
b^{2+p}$, hence the conclusion.

A similar argument applies to $Z\subset Y_{\max}$, assuming that
$|Y_{\max}|\geq b^{2+\operatorname{card}_{b}(Y)}$. In this case, if $Z\supset
Y$ then we have $\sum_{x\in Y}\alpha_{x}\geq\lambda$. As
$|Y|<b^{1+\operatorname{card}_{b}(Y)}$, $Z$ also contains more than
$b^{2+\operatorname{card}_{b}(Y)}-b^{1+\operatorname{card}_{b}(Y)}=b^{2+\operatorname{card}_{b}(Y)}\cdot\left(1-\frac{1}{b}\right)$
points from pairs $(Y,Y^{\prime})\in\operatorname{metaclusters}$,
$Y^{\prime}\neq Y$. By the argument for $Y_{p}$, for each such point $x$ we
have
$\alpha_{x}>\frac{\lambda}{b^{1+\operatorname{card}_{b}(Y^{\prime})}}>\frac{\lambda}{b^{1+\operatorname{card}_{b}(Y)}}$.
Overall, we get that $\sum_{x\in
Z}\alpha_{x}>\lambda+b^{2+\operatorname{card}_{b}(Y)}\cdot\left(1-\frac{1}{b}\right)\cdot\frac{\lambda}{b^{1+\operatorname{card}_{b}(Y)}}=b\cdot\lambda$.
If $Z$ does not contain $Y$, then the argument for $Y_{p}$ holds. ∎

######  Lemma 15.

For every $Z\in\operatorname{clusters}$, $\operatorname{cost}(Z)\leq
2b^{2}\cdot\sum_{x\in Z}d(x,Z)$.

###### Proof.

We have $\operatorname{cost}(Z)=|Z|\cdot\sum_{x\in
Z}\|x-\operatorname{cm}(Z)\|_{2}^{2}\leq|Z|\cdot\sum_{x\in
Z}\|x-\operatorname{ctr}(Y)\|_{2}^{2}$. By Lemmas 6 and 7, $|Z|\leq
2b^{2}\cdot\operatorname{card}_{b}(Z)$. ∎

Proof of Theorem 10. First, consider the feasibility of $(\alpha,\gamma)$.
Clearly, $\gamma$ is set in the theorem to satisfy the constraints that
include it. Regarding the constraints that involve only $\alpha$, we prove
that they are satisfied throughout the execution of procedure
PrimalDualPhase1. The proof is by induction on the number of inactive points.
Clearly, the initial $\alpha$ is feasible. Now, suppose that $\alpha$ is
feasible for some number of inactive points, and consider the next step when
this number increases and a set $A\subset\operatorname{active}$ is removed
from $\operatorname{active}$ (we will use $\operatorname{active}$ here to
denote the set before the removal of $A$). Let $\alpha^{\prime}$ denote the
values of the dual variables just before $A$ is removed from
$\operatorname{active}$. If there exist $Y\subset X$ and $y\in Y$ such that
the constraint for the pair $Y,y$ is violated, then clearly
$Y\cap\operatorname{active}\neq\emptyset$, otherwise the same constraint would
have been violated by the solution $\alpha$, as $\alpha$ and $\alpha^{\prime}$
differ only on $\operatorname{active}$. But then there is some intermediate
value $\alpha^{\prime\prime}$ such that $\alpha^{\prime\prime}_{x}=\alpha_{x}$
for all $x\not\in\operatorname{active}$ and
$\alpha_{x}\leq\alpha^{\prime\prime}_{x}<\alpha^{\prime}_{x}$ for all
$x\in\operatorname{active}$, which causes this constraint (or another
constraint involving active points) to become tight. Therefore, at least one
point would have been removed from $\operatorname{active}$ before we reach the
values $\alpha^{\prime}$, in contradiction with our assumptions.

Next, consider the number of points clustered in the output
$\operatorname{clusters}$ of procedure PrimalDual. Clearly, procedure
PrimalDualPhase2 clusters in $\operatorname{metaclusters}$ exactly
$n^{\prime}$ points. Some of these points are discarded by procedure
PrimalDualPhase3. Consider some $(Y,\cdot)\in\operatorname{metaclusters}$. By
Lemma 7, the number of points discarded from these clusters is less than

$\sum_{p<\operatorname{card}_{b}(Y)-2}b^{2+p}=\frac{b^{\operatorname{card}_{b}(Y)}-b^{2}}{b-1}.$

On the other hand, all the points in $Y_{\max}$ are clustered in
$\operatorname{clusters}$, as $Y\neq Y_{\operatorname{last}}$ in this case.
Clearly, the number of points in $Y_{\max}$ is at least $|Y|\geq
b^{\operatorname{card}_{b}(Y)}$. Thus, less than $\frac{1}{b-1}\cdot
n^{\prime}\leq\epsilon\cdot n^{\prime}$ points get discarded.

Finally, let’s consider the cost of the clustering. Let
$Z\in\operatorname{clusters}$ be a cluster that satisfies the conditions of
Lemma 14. Then,

$\displaystyle\sum_{x\in Z}\alpha_{x}-\lambda$ $\displaystyle>$
$\displaystyle\frac{b-1}{b}\cdot\sum_{x\in Z}\alpha_{x}$ $\displaystyle\geq$
$\displaystyle\frac{b-1}{9b}\cdot\sum_{x\in Z}d(x,Z)$ $\displaystyle\geq$
$\displaystyle\frac{b-1}{18b^{3}}\cdot\operatorname{cost}(Z),$

where the first inequality follows from Lemma 14, the second inequality
follows from Lemma 13, and the third inequality follows from Lemma 15. The
remaining clusters are sets $Y_{\max}$ with
$|Y_{\max}|<b^{2+\operatorname{card}_{b}(Y)}$ and a subset of
$Y_{\operatorname{last}}$. Consider a cluster
$Y_{\max}\in\operatorname{clusters}$. We have that

$\sum_{x\in Y_{\max}}\alpha_{x}-\lambda=\sum_{x\in
Y}\alpha_{x}-\lambda+\sum_{x\in Y_{\max}\setminus Y}\alpha_{x}\geq\sum_{x\in
Y}d(x,Y)+\frac{1}{9}\cdot\sum_{x\in Y_{\max}\setminus
Y}d(x,Y)\geq\frac{1}{18b^{2}}\cdot\operatorname{cost}(Y_{\max}).$

Finally, if there’s a cluster $Z\subset Y_{\operatorname{last}}$, then

$\sum_{x\in Z}\alpha_{x}\geq\frac{1}{9}\cdot\sum_{x\in
Z}d(x,Y_{\max})\geq\frac{1}{18b^{2}}\cdot\operatorname{cost}(Z).$

Thus, we can set
$c=c_{\epsilon}=\frac{18b^{3}}{b-1}=\frac{18(1+\epsilon)^{3}}{\epsilon^{2}}$.
∎

1:procedure PrimalDual($X,\lambda,n^{\prime},b$)

2:
$\alpha,\operatorname{preclusters},Y_{\operatorname{last}}\leftarrow$PrimalDualPhase1($X,\lambda,n^{\prime},b$)

3:
$\operatorname{metaclusters}\leftarrow$PrimalDualPhase2($X,n^{\prime},b,\alpha,\operatorname{preclusters},Y_{\operatorname{last}}$)

4:
$\operatorname{clusters}\leftarrow$PrimalDualPhase3($X,b,\operatorname{metaclusters}$)

5: return $\operatorname{clusters}$

6:end procedure

7:

8:procedure PrimalDualPhase1($X,\lambda,n^{\prime},b$)

9: $\operatorname{active},\operatorname{preclusters}\leftarrow X,\emptyset$

10: $\alpha_{x}\leftarrow 0$ for all $x\in X$

11: while $|\operatorname{active}|>n-n^{\prime}$ do

12: raise $\alpha_{x}$ at a uniform rate for all
$x\in\operatorname{active}$$\triangleright$ stop raising when one of the
following two cases happens

13: if $\exists x\in\operatorname{active}$ and
$Y\in\operatorname{preclusters}$ such that $\alpha_{x}\geq d(x,Y)$ then

14:
$\operatorname{card}_{b}(Y\cup\\{x\\}),\operatorname{ctr}(Y\cup\\{x\\})\leftarrow\operatorname{card}_{b}(Y),\operatorname{ctr}(Y)$

15:
$\operatorname{preclusters},\operatorname{active}\leftarrow\operatorname{preclusters}\setminus\\{Y\\}\cup\\{Y\cup\\{x\\}\\},\operatorname{active}\setminus\\{x\\}$

16: else if $\exists Y\subset X$ and $y\in Y$ such that
$Y\cap\operatorname{active}\neq\emptyset$ the dual constraint for $Y,y$ is
tight then

17:
$\operatorname{card}_{b}(Y),\operatorname{ctr}(Y)\leftarrow\log_{b}\operatorname{floor}_{b}(|Y|),y$

18: if $|\operatorname{active}\setminus Y|<n-n^{\prime}$ then $\triangleright$
use an inclusion-wise minimal such $Y$

19: return $\alpha,\operatorname{preclusters},Y$

20: else

21:
$\operatorname{preclusters},\operatorname{active}\leftarrow\operatorname{preclusters}\cup\\{Y\\},\operatorname{active}\setminus
Y$

22: end if

23: end if

24: end while

25: return $\alpha,\operatorname{preclusters},\emptyset$

26:end procedure

27:

28:procedure
PrimalDualPhase2($X,n^{\prime},b,\alpha,\operatorname{preclusters},Y_{\operatorname{last}}$)

29: $\operatorname{active},\operatorname{metaclusters}\leftarrow\\{x\in X:\
x\in Y\in\operatorname{preclusters}\vee x\in
Y_{\operatorname{last}}\\},\emptyset$

30: for $Y\in\operatorname{preclusters}$, by order of nonicreasing
$\operatorname{card}_{b}(Y)$ do

31: if $\exists(Y^{\prime},Y^{\prime})\in\operatorname{metaclusters}$ with
$x\in Y\cap Y^{\prime}$ and
$\alpha_{x}>\max\left\\{d(x,Y),d(x,Y^{\prime})\right\\}$ then

32:
$Y^{\prime\prime},\operatorname{card}_{b}(Y^{\prime\prime}),\operatorname{ctr}(Y^{\prime\prime})\leftarrow\\{x\in
Y\cap\operatorname{active}:\ \alpha_{x}=\max_{y\in
Y}\alpha_{y}\\},\operatorname{card}_{b}(Y),\operatorname{ctr}(Y^{\prime})$

33:
$\operatorname{metaclusters}\leftarrow\operatorname{metaclusters}\cup\\{(Y^{\prime},Y^{\prime\prime})\\}$

34: $\operatorname{active}\leftarrow\operatorname{active}\setminus
Y^{\prime\prime}$

35: else

36: remove each $x\in Y$ from any $Y^{\prime\prime}\ni x$ with
$(Y^{\prime},Y^{\prime\prime})\in\operatorname{metaclusters}$ $\triangleright$
$\alpha_{x}\leq d(x,Y^{\prime\prime})$;
$\operatorname{card}_{b}(Y^{\prime\prime}),\operatorname{ctr}(Y^{\prime\prime})$
don’t change

37:
$\operatorname{metaclusters}\leftarrow\operatorname{metaclusters}\cup\\{(Y,Y)\\}$

38: $\operatorname{active}\leftarrow\operatorname{active}\setminus Y$

39: end if

40: end for

41:
$Y,\operatorname{card}_{b}(Y),\operatorname{ctr}(Y)\leftarrow\\{|\operatorname{active}|-n+n^{\prime}$
points in
$Y_{\operatorname{last}}\cap\operatorname{active}\\},\operatorname{card}_{b}(Y_{\operatorname{last}}),\operatorname{ctr}(Y_{\operatorname{last}})$

42: if $Y\neq\emptyset$ then

43:
$\operatorname{metaclusters}\leftarrow\operatorname{metaclusters}\cup\\{(Y_{\operatorname{last}},Y)\\}$

44: end if

45: return $\operatorname{metaclusters}$

46:end procedure

47:

48:procedure PrimalDualPhase3($X,b,\operatorname{metaclusters}$)

49: $\operatorname{clusters}\leftarrow\emptyset$

50: for $(Y,\cdot)\in\operatorname{metaclusters}$ do

51: $Y_{\max}\leftarrow\\{x\in X:\ \exists Y^{\prime}\ni x\hbox{ such that
}(Y,Y^{\prime})\in\operatorname{metaclusters}\wedge\operatorname{card}_{b}(Y^{\prime})\geq\operatorname{card}_{b}(Y)-2\\}$

52: $\operatorname{clusters}\leftarrow\operatorname{clusters}\cup$
Partition($Y_{\max},\max\\{1,\left\lfloor|Y_{\max}|/b^{2+\operatorname{card}_{b}(Y)}\right\rfloor\\}$)

53: for $p<\operatorname{card}_{b}(Y)-2$ do

54: $Y_{p}\leftarrow\\{x\in X:\ \exists Y^{\prime}\ni x\hbox{ such that
}(Y,Y^{\prime})\in\operatorname{metaclusters}\wedge\operatorname{card}_{b}(Y^{\prime})=p\\}$

55: if $|Y_{p}|\geq b^{2+p}$ then

56: $\operatorname{clusters}\leftarrow\operatorname{clusters}\cup$
Partition($Y_{p},\left\lfloor|Y_{p}|/b^{2+p}\right\rfloor$)

57: end if

58: end for

59: end for

60: return $\operatorname{clusters}$

61:end procedure

62:

63:procedure Partition($S,m$)$\triangleright$ $m\geq 1$

64: partition $S$ as evenly as possible into $m$ disjoint subsets
$S_{1},S_{2},\dots,S_{m}$

65: return $\\{S_{1},S_{2},\dots,S_{m}\\}$

66:end procedure

Algorithm 1 Algorithm Primal-Dual

1:procedure MinSumClustering($X,k,n^{\prime},\epsilon$)

2: $\lambda_{1}\leftarrow 0$, $\lambda_{2}\leftarrow\sum_{x,y\in
X}\|x-y\|_{2}^{2}$

3: $\operatorname{clusters}_{1}\leftarrow\\{\\{x\\}:\ x\in X\\}$,
$\operatorname{clusters}_{2}\leftarrow X$

4: $b\leftarrow\frac{1+\epsilon}{\epsilon}$

5: $\delta\leftarrow\frac{2}{(n+k)\lambda_{2}}$$\triangleright$ we need
$\delta\leq\frac{2}{(n+k)\operatorname{opt}(X,n^{\prime},k)}$

6: while $\lambda_{2}-\lambda_{1}>\delta$ do

7: $\lambda=\frac{1}{2}\cdot(\lambda_{1}+\lambda_{2})$

8:
$\operatorname{clusters}\leftarrow\hbox{PrimalDual}(X,\lambda,n^{\prime},b)$

9: $Y_{\min,\lambda}\leftarrow$ smallest cluster in $\operatorname{clusters}$

10: $k^{\prime}\leftarrow|\operatorname{clusters}|-1$

11: if $|Y_{\min,\lambda}|\leq\frac{\epsilon}{3}\cdot n^{\prime}$ then

12:
$\operatorname{clusters}\leftarrow\operatorname{clusters}\setminus\\{Y_{\min,\lambda}\\}$

13: end if

14: if $k^{\prime}>k$ then

15:
$\lambda_{1},\operatorname{clusters}_{1},k_{1}\leftarrow\lambda,\operatorname{clusters},k^{\prime}$

16: else$\triangleright$ $k^{\prime}\leq k$

17:
$\lambda_{2},\operatorname{clusters}_{2},k_{2}\leftarrow\lambda,\operatorname{clusters},k^{\prime}$

18: end if

19: end while

20: $\rho_{1}\leftarrow\frac{k-k_{2}}{k_{1}-k_{2}}$ $\triangleright$
$k_{1}>k\geq k_{2}\geq 0$

21: if $\rho_{1}\geq 1-\frac{\epsilon}{4}$ then

22: return $\\{k$ largest sets in $\operatorname{clusters}_{1}\\}$

23: else

24: return $\operatorname{clusters}_{2}$ $\triangleright$ If
$|\operatorname{clusters}_{2}|<k$, split clusters arbitrarily to get exactly
$k$

25: end if

26:end procedure

Algorithm 2 Algorithm Min-Sum-Clustering

## References

  * [ABS10] P. Awasthi, A. Blum, and O. Sheffet.  Stability yields a PTAS for $k$-median and $k$-means clustering.  In Proc. of the 51st Ann. IEEE Symp. on Foundations of Computer Science, pages 309–318, 2010.
  * [ADHP09] D. Aloise, A. Deshpande, P. Hansen, and P. Popat.  NP-hardness of Euclidean sum-of-squares clustering.  Machine Learning, 75(2):245–248, May 2009.
  * [AJM09] N. Ailon, R. Jaiswal, and C. Monteleoni.  Streaming $k$-means approximation.  In Proc. of the 23rd Ann. Conf. on Neural Information Processing Systems, pages 10–18, 2009.
  * [AMR11] D. Arthur, B. Manthey, and H. Röglin.  Smoothed analysis of the $k$-means method.  J. ACM, 58(5):19:1–19:31, October 2011.
  * [ANFSW17] S. Ahmadian, A. Norouzi-Fard, O. Svensson, and J. Ward.  Better guarantees for $k$-means and Euclidean $k$-median by primal-dual algorithms.  Proc. of the 58th Ann. IEEE Symp. on Foundations of Computer Science, pages 61–72, 2017.
  * [AS16] S. Ahmadian and C. Swamy.  Approximation algorithms for clustering problems with lower bounds and outliers.  In Proc. of the 43rd Int’l Colloq. on Automata, Languages, and Programming, pages 69:1–69:15, 2016.
  * [AV07] D. Arthur and S. Vassilvitskii.  $k$-means++: The advantages of careful seeding.  In Proc. of the 18th Ann. ACM-SIAM Symp. on Discrete Algorithms, pages 1027–1035, 2007.
  * [Bar96] Y. Bartal.  Probabilistic approximation of metric spaces and its algorithmic applications.  In Proc. of the 37th Ann. IEEE Symp. on Foundations of Computer Science, page 184, 1996.
  * [Bar98] Y. Bartal.  On approximating arbitrary metrices by tree metrics.  In Proc. of the 30th Ann. ACM Symp. on Theory of Computing, pages 161–168, 1998.
  * [BB09] M.-F. Balcan and M. Braverman.  Finding low error clusterings.  In COLT 2009 - The 22nd Conference on Learning Theory, 2009.
  * [BBG09] M.-F. Balcan, A. Blum, and A. Gupta.  Approximate clustering without the approximation.  In Proc. of the 20th Ann. ACM-SIAM Symp. on Discrete Algorithms, pages 1068–1077, 2009.
  * [BCR01] Y. Bartal, M. Charikar, and D. Raz.  Approximating min-sum $k$-clustering in metric spaces.  In Proc. of the 33rd Ann. ACM Symp. on Theory of Computing, pages 11–20, 2001.
  * [BFL16] V. Braverman, D. Feldman, and H. Lang.  New frameworks for offline and streaming coreset constructions.  CoRR, abs/1612.00889, 2016.
  * [BFSS19] B. Behsaz, Z. Friggstad, M. R. Salavatipour, and R. Sivakumar.  Approximation algorithms for min-sum $k$-clustering and balanced $k$-median.  Algorithmica, 81(3):1006–1030, 2019.
  * [BMO+11] V. Braverman, A. Meyerson, R. Ostrovsky, A. Roytman, M. Shindler, and B. Tagiku.  Streaming $k$-means on well-clusterable data.  In Proc. of the 22nd Ann. ACM-SIAM Symp. on Discrete Algorithms, pages 26–40, 2011.
  * [BO07] P. Bunn and R. Ostrovsky.  Secure two-party $k$-means clustering.  In Proc. of the 14th Ann. ACM Conf. on Computer and Communications Security, pages 486–497, 2007.
  * [CAKM16] V. Cohen-Addad, P. N. Klein, and C. Mathieu.  Local search yields approximation schemes for $k$-means and $k$-median in Euclidean and minor-free metrics.  Proc. of the 57th Ann. IEEE Symp. on Foundations of Computer Science, pages 353–364, 2016.
  * [CAS17] V. Cohen-Addad and C. Schwiegelshohn.  On the local structure of stable clustering instances.  pages 49–60, 2017.
  * [CC19] V. Cohen-Addad and Karthik C.S.  Inapproximability of clustering in $l_{p}$ metrics.  In Proc. of the 60th Ann. IEEE Symp. on Foundations of Computer Science, pages 519–539, 2019.
  * [CCL19] V. Cohen-Addad, Karthik C. S., and E. Lee.  On approximability of $k$-means, $k$-median, and $k$-minsum clustering, 2019.
  * [Che09] K. Chen.  On coresets for $k$-median and $k$-means clustering in metric and Euclidean spaces and their applications.  SIAM J. Comput., 39:923–947, 2009.
  * [CR05] J. Chuzhoy and Y. Rabani.  Approximating $k$-median with non-uniform capacities.  In Proc. of the 16th Ann. ACM-SIAM Symp. on Discrete Algorithms, page 952–958, 2005.
  * [CS07] A. Czumaj and C. Sohler.  Small space representations for metric min-sum $k$-clustering and their applications.  In Proc. of the 24th Ann. Conf. on Theoretical Aspects of Computer Science, pages 536–548, 2007.
  * [DL16] H. G. Demirci and S. Li.  Constant approximation for capacitated $k$-median with $(1+\epsilon)$-capacity violation.  ArXiv, abs/1603.02324, 2016.
  * [FdlVKKR03] W. Fernandez de la Vega, M. Karpinski, C. Kenyon, and Y. Rabani.  Approximation schemes for clustering problems.  In Proc. of the 35th Ann. ACM Symp. on Theory of Computing, pages 50–58, 2003.
  * [FK01] W. Fernandez de la Vega and C. Kenyon.  A randomized approximation scheme for metric MAX-CUT.  Journal of Computer and System Sciences, 63(4):531–541, 2001.
  * [FL11] D. Feldman and M. Langberg.  A unified framework for approximating and clustering data.  In Proc. of the 43rd Ann. ACM Symp. on Theory of Computing, pages 569–578, 2011.
  * [FMS07] D. Feldman, M. Monemizadeh, and C. Sohler.  A PTAS for $k$-means clustering based on weak coresets.  In Proc. of the 23rd Ann. Symp. on Computational Geometry, pages 11–18, 2007.
  * [FRS16] Z. Friggstad, M. Rezapour, and M. R. Salavatipour.  Local search yields a PTAS for $k$-means in doubling metrics.  Proc. of the 57th Ann. IEEE Symp. on Foundations of Computer Science, pages 365–374, 2016.
  * [FRT03] J. Fakcharoenphol, S. Rao, and K. Talwar.  A tight bound on approximating arbitrary metrics by tree metrics.  In Proc. of the 35th Ann. ACM Symp. on Theory of Computing, pages 448–455, 2003.
  * [GH98] N. Guttmann-Beck and R. Hassin.  Approximation algorithms for min-sum $p$-clustering.  Discret. Appl. Math., 89(1-3):125–142, 1998.
  * [GI03] V. Guruswami and P. Indyk.  Embeddings and non-approximability of geometric problems.  In Proc. of the 14th Ann. ACM-SIAM Symp. on Discrete Algorithms, pages 537–538, 2003.
  * [HPK07] S. Har-Peled and A. Kushal.  Smaller coresets for $k$-median and $k$-means clustering.  Discrete Comput. Geom., 37(1):3–19, January 2007.
  * [IKI94] M. Inaba, N. Katoh, and H. Imai.  Applications of weighted Voronoi diagrams and randomization to variance-based $k$-clustering.  In Proc. of the 10th Ann. Symp. on Computational Geometry, pages 332–339, 1994.
  * [Ind99] P. Indyk.  A sublinear time approximation scheme for clustering in metric spaces.  In Proc. of th 40th Ann. IEEE Symp. on Foundations of Computer Science, pages 154–159, 1999.
  * [JG12] R. Jaiswal and N. Garg.  Analysis of k-means++ for separable data.  In Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques, pages 591–602, 2012.
  * [JV01] K. Jain and V. V. Vazirani.  Approximation algorithms for metric facility location and $k$-median problems using the primal-dual schema and Lagrangian relaxation.  J. ACM, 48(2):274–296, 2001.
  * [KK10] A. Kumar and R. Kannan.  Clustering with spectral norm and the $k$-means algorithm.  In Proc. of the 51st Ann. IEEE Symp. on Foundations of Computer Science, pages 299–308, 2010.
  * [Kle02] J. Kleinberg.  An impossibility theorem for clustering.  In Proc. of the 15th Int’l Conf. on Neural Information Processing Systems, pages 463–470, 2002.
  * [KSS05] A. Kumar, Y. Sabharwal, and S. Sen.  Linear time algorithms for clustering problems in any dimensions.  In Proc. of the 32nd Int’l Conf. on Automata, Languages and Programming, pages 1374–1385, 2005.
  * [Mat00] J. Matoušek.  On approximate geometric $k$-clustering.  Discrete & Computational Geometry, 24(1):61–84, Jan 2000.
  * [MNV09] M. Mahajan, P. Nimbhorkar, and K. Varadarajan.  The planar $k$-means problem is NP-hard.  In Proc. of the 3rd Int’l Workshop on Algorithms and Computation, pages 274–285, 2009.
  * [ORSS06] R. Ostrovsky, Y. Rabani, L. J. Schulman, and C. Swamy.  The effectiveness of Lloyd-type methods for the $k$-means problem.  In Proc. of the 47th Ann. IEEE Symp. on Foundations of Computer Science, pages 165–176, 2006.
  * [Sch00] L. J. Schulman.  Clustering for edge-cost minimization.  In Proc. of the 32nd Ann. ACM Symp. on Theory of Computing, pages 547–555, 2000.
  * [VBR+11] K. Voevodski, M.-F. Balcan, H. Röglin, S.-H. Teng, and Y. Xia.  Min-sum clustering of protein sequences with limited distance information.  In Proc. of the 1st Int’l Conf. on Similarity-Based Pattern Recognition, pages 192–206, 2011.
  * [ZBD09] R. B. Zadeh and S. Ben-David.  A uniqueness theorem for clustering.  In Proc. of the 25th Ann. Conf. on Uncertainty in Artificial Intelligence, pages 639–646, 2009.


11institutetext:  Bio-Inspired System Design Lab, Department of Biomedical
Engineering, Amirkabir University of Technology (Tehran Polytechnic), Tehran
15914, Iran
22institutetext: Department of Electrical Engineering, Amirkabir University of
Technology (Tehran Polytechnic), Tehran 15914, Iran
correspondingauthor: 22email<EMAIL_ADDRESS>
# AVIS Engine: High Performance Networking Layer for Simulation Applications

Amir Mohammad Zarif Shahsavan Nejad 11 0000-0003-2969-1573 Amir Mahdi Zarif
Shahsavan Nejad 11 Amirali Setayeshi 22 0009-0000-2894-7640 Soroush
Sadeghnejad 11 0000-0001-8286-7756

###### Abstract

Autonomous vehicles are one of the most popular and also fast-growing
technologies in the world. As we go further, there are still a lot of
challenges that are unsolved and may cause problems in the future when it
comes to testing in real world. Simulations on the other hand have always had
a huge impact in the fields of science, technology, physics, etc. The
simulation also powers real-world Autonomous Vehicles nowadays. Therefore, We
have built an Autonomous Vehicle Simulation Software - called AVIS Engine -
that provides tools and features that help develop autonomous vehicles in
various environments. AVIS Engine features an advanced input and output system
for the vehicle and includes a traffic system and vehicle sensor system which
can be communicated using the fast networking system and ROS Bridge.

###### Keywords:

Simulation Computer Networking Optimization Artificial Intelligent

##  1 Introduction

Artificial Intelligence (AI) is one of the biggest topics in the world today
and plays an important role in people’s lives. Hence, it should be more and
more intelligent and safe to be able play its role perfectly. Simulations have
always been a part of robotics and AI such as industrial robots, surgical and
medical robots [25], humanoid robots, etc. which can become extremely complex
and difficult as we progress [8] so is the RoboCup goal for 2050. RoboCup also
involved simulators for wheeled robot since the beginnings [12, 16] and also
made the simulation to more robotics research groups till this day [12].
Autonomous vehicles (AVs) are one of the applications of AI in the real world
and are designed to help people travel safer, faster, easier, and cheaper than
non-electric and non-autonomous vehicles. AVs also promise to be more
environmentally friendly, as they are electrically powered and produce fewer
pollutants than internal combustion engine vehicles [21].

In most cases, many experiments are required to approximately guarantee the
safety of these modern AVs. The experiments can also be expensive. Testing and
experimenting in simulation environments is a solution to build a much more
secure AV that costs less compared to building it in the real world.

Simulating different objects and machines can be very beneficial because it is
easy to study them. On the other hand, real-world research can be very
expensive. For example, in the field of AVs, experimenting with a real car is
much more expensive and time consuming than a simulation given the cost of the
tools and equipment needed to build one AV. Simulations can be beneficial in
studying algorithms, training machine learning models, finding different
solutions and approaches to solving a problem, and identifying "driving"
variables. [20]

For instance, obstacle avoidance or path planning in the real physical form of
the vehicle can be extremely difficult and also expensive, and certainly
errors can have some unavoidable consequences. Usually, testing algorithms in
simulator can be very time-saving and also cheaper than other solutions.
Resetting and quickly testing other conditions without cost is one of the
other advantages of using simulators, which can be useful in machine learning
and deep learning solutions.

All different scenarios and conditions should be considered when a technology
is based on AI, and one solution to this is the use of simulations. Some
researchers use videos to develop lane recognition keeping algorithms, and
there is no other input or output from the vehicle. But a simulator can
actually process those inputs and outputs. Nowadays, the trend of AI is
growing so fast and there are so many things that are related to AI in some
way and have been extended to an application of AI.

There are other simulators for autonomous vehicles that also aim to support
and advance research in this area, including Apollo [5], AutonoViSim [6],
CARLA [11], Force-based, GTAV, SimMobilityST [17], SUMMIT [9], SUMO [19] and
TORCS [2]. Table 1 shows a comparison between the existing simulators and
platforms.

Table 1: A Comparison between existing simulators (Derived from [9]) Simulator | Real-world Maps | Unregulated Behaviors | Dense Traffic | Realistic Visuals
---|---|---|---|---
Apollo | $\times$ | $\times$ | $\times$ |
AutonoViSim | $\times$ |  |  |
CARLA | $\times$ | $\times$ | $\times$ |
Force-based | $\times$ | $\times$ |  |
GTAV | $\times$ | $\times$ | $\times$ |
SimMobilityST |  | $\times$ |  | $\times$
SUMMIT |  |  |  |
SUMO | $\times$ | $\times$ |  |
AVIS Engine (Ours) |  |  |  |

We have developed simulation software for AVs called AVIS Engine (Autonomous
Vehicles Intelligent Simulation Software)111AVIS Engine Website:
https://avisengine.com that allows researchers to explore get their hands on
the research and development of AVs in a way that makes the process truly
feasible and fast. CARLA and Microsoft AirSim [26] that are notable in the AV
simulation. Both Microsoft AirSim and CARLA are powered by the Unreal Engine.
Both simulators provide powerful realistic 3D simulations. AVIS Engine uses
Unity 3D as its main engine. Both Unreal Engine and Unity 3D are suitable for
high-fidelity rendering in applications such as robotics simulations. [13]).
The AVIS Engine can also be used to experiment with the sound of the
environment, i.e. it provides multiple listeners placed around the vehicle as
microphones that pick up the sound of any object in the scene such as other
vehicles, pedestrians, ambulances and other criteria.

The AVIS Engine includes most of the features required for the AI development
of the AVs are helpful for researchers and developers developing AVs and AI.

##  2 AVIS Engine

Figure 1: A screenshot of the main urban environment in AVIS Engine. The panel
at the top of the screen is used for configuration and visualization of the
necessary data of the vehicle

AVIS Engine is a robust and fast simulation software for autonomous vehicles
that extremely simplifies the development of AVs and is quite easy to use.
Figure 1 shows the main interface of the software. The simulator uses its own
messaging system that relies on the Transmission Control Protocol (TCP). The
program transmits data and commands over a TCP connection where the simulator
is the server peer of the connection and the program is the client peer of the
connection [bib8].

In order to ensure reliable, practical, and secure data transfer and to allow
smooth communication between the server (simulator) and the client (script),
we have chosen TCP instead of UDP for the implementation of this system. In
general, clients can be written in any technology or programming language that
supports TCP. So far, AVIS Engine communicates with Python, C++, ROS (Robotic
Operating System) and MATLAB [1] and is flexible to support other technologies
in the future.

The simulator itself is based on Unity [3] and the C# programming language,
which enabled the development of powerful cross-platform simulation software.
The simulator includes and offers a wide range of features and key values, to
name a few:

###  2.1 High Performance

AVIS Engine is compatible with low-end PCs and is capable of running smoothly
on them. This is because the software core is lightweight and the underlying
network system is fast and optimized for such scenarios.

###  2.2 Fast Networking

AVIS Engine messaging system aims to deliver and transmit data in an efficient
manner, resulting in a fast and agile network. Section 3 explains the
networking of the simulator.

###  2.3 Global Positioning System (GPS)

The GPS unit of AVs is one of the most important components that plays an
important role in guiding and navigating the vehicle and can be used in
various applications such as navigation systems and path planning. The cost of
precise GPS receivers is the major drawback of GPS for these autonomous
applications [10].

###  2.4 Semantic Camera

A semantic segmentation that can be used for many purposes, such as problems
related to machine learning and deep learning. It can also be used to improve
photorealism (Figure 2).

Figure 2: Semantic view of the front camera

###  2.5 Depth camera

Concretely, there is no depth information in a normal camera, this is where
the depth camera comes to the rescue. The depth camera is used for depth
estimation and provides a better understanding of the environment. With this
feature, a third dimension is added to the camera information besides the 2-D
image, namely the depth (Figure 3).

Figure 3: Depth camera view in the simulator interface

###  2.6 Built-in filter for road features

In addition to the various cameras and sensors, the simulator has a built-in
road feature filter that highlights road features captured from the camera
(Figure 4).

Figure 4: Road feature filter view

###  2.7 Depth camera

Concretely, there is no depth information in a normal camera, this is where
the depth camera comes to the rescue. The depth camera is used for depth
estimation and provides a better understanding of the environment. With this
feature, a third dimension is added to the camera information besides the 2-D
image, namely the depth (Figure 3).

###  2.8 Built-in filter for road features

In addition to the various cameras and sensors, the simulator has a built-in
road feature filter that highlights road features captured from the camera
(Figure 4).

###  2.9 Light Detection and Ranging (LIDAR)

Detect and understand the environment in even greater detail. The sensor has
been simulated as a set of raycasts and the output is represented as a point
cloud and can be visualized in RViz [15] using the AVIS Engine ROS Bridge
(Figure 5). It is also customizable in 4 general specifications [24]
including:

  1. 1.

Ranging Specification

  2. 2.

Physical Specifications

  3. 3.

Laser Specification

  4. 4.

Optical Specification

Figure 5: Visualization of the LIDAR Sensor

###  2.10 Traffic System

Driving in an empty city with no other vehicles on the roads is not much of a
challenge. But when it comes to driving between other cars on the road, it
becomes even more difficult. AVIS Engine provides a customizable traffic
system that allows a more detailed study of AV (Figure 6).

###  2.11 Vehicle communication systems

Inter-vehicle communication is one of the most important components for AVs
and smart vehicles, and hopes to cover and improve issues around vehicle
safety. [18]. The AVIS Engine has a system that provides a vehicle
communication system for vehicle communication studies.

Figure 6: AVIS Engine Traffic System

###  2.12 ROS Bridge

ROS is a powerful software development kit for robotics [4, 23]. The AVIS
Engine features a compatibility layer that connects directly to ROS and
includes many built-in tools as well as community-developed tools and
packages. This allows for efficient and practical visualization, development
of control systems, simulations, and hardware implementations for which comple

##  3 Networking

Figure 7: Founding Layers and various parts of the simulator

The AVIS Engine is built on two basic building blocks: 1. AVIS Engine Server
software, 2. AVIS Engine Client (Figure 7). The server software is responsible
for the physics simulations, the 3D environment, the sound, the vehicle model,
the sensor simulations and all simulations related to the vehicle. The client,
on the other hand, is responsible for sending the client commands to the
server and receiving the response through an interface. The client can be
written in any programming language or platform that can establish a TCP
connection. The most important and extensive API developed for AVIS Engine is
written in Python programming language. Communication between these two basic
building blocks requires networking and is also the bottleneck of the whole
system.

Therefore, a fast networking system is crucial to achieve fast and stable
communication between the client and the server. Many issues must be
considered in its development. Dealing with different approaches to improve
and optimize the network such as compression algorithms.

We have experimented with different approaches to optimize the network layer
of the software, since a fast networking system is crucial. In early versions
of AVIS Engine, the software used Portable Network Graphics (PNG) format to
transmit images from the vehicle’s camera sensors. It was also low resolution
(256 pixels * 256 pixels) and slow overall. After performing techniques, the
simulator is now able to transmit 1080p images to the client without any speed
degradation. Here are some techniques that were used to optimize the network
layer:

###  3.1 Image Compression

The first and most effective step in the optimization process was the lossy
compression and conversion of the camera images to JPEG 222Joint Photographic
Experts Group format before they were sent over. This technique had a great
impact on increasing the speed of data transmission and reducing the delay. An
example of the reduced image of the camera sensor is described in the Table 2.

###  3.2 Color-Space Conversion

Color-space conversion is a widely used technique and is frequently used in
computer vision and artificial intelligence (see [7] for an example of an
artificial intelligence study). Also used for television systems and high-
definition television. This technique allowed us to slightly compress the
image camera data in the simulator, but it resulted in data loss. The
conversion to YCbCr is derived as follows and is defined by ITU-R BT.601 [22,
14, 27]:

  1. 1.

RGB to YCbCr Conversion:

$\begin{bmatrix}Y^{\prime}\\\ Cb\\\
Cr\end{bmatrix}=\begin{bmatrix}0.299&0.587&0.114\\\ -0.169&-0.331&0.499\\\
0.499&-0.418&-0.0813\end{bmatrix}\begin{bmatrix}R\\\ G\\\
B\end{bmatrix}+\begin{bmatrix}0\\\ 128\\\ 128\end{bmatrix}$ (1)

Which can be also written as:

$\begin{cases}Y=0.299R+0.587G+0.114B\\\ Cb=-0.169R+-0.332G+0.500B+128\\\
Cr=0.500R+-0.419G+-0.0813B+128\end{cases}$ (2)

  2. 2.

YCbCr to RGB Conversion:

$\begin{bmatrix}R\\\ G\\\ B\end{bmatrix}=\begin{bmatrix}1&0&1.402\\\
1&-0.334&-0.714\\\ 1&1.772&0\end{bmatrix}\begin{bmatrix}Y\\\ Cb-128\\\
Cr-128\end{bmatrix}$ (3)

Which can be written as:

$\begin{cases}R=Y+1.407\times(V-128)\\\
G=Y-0.345\times(Cb-128)-(0.716\times(Cr-128))\\\
B=Y+1.779\times(Cb-128)\end{cases}$ (4)

###  3.3 GZip Compression

This method is usually effective in compressing the data. However, in the case
of our simulator, it was unexpectedly ineffective and the change in data size
was small.

Table 2: Detailed information about how much data was compressed after each technique was performed Technique |  | Raw Camera Data
---
Uncompressed
(303459 Bytes)
Change in Data Size |  |
Image Compression (PNG to JPEG) | 32932 Bytes |  $\sim$90% Reduction |  |
Color-Space Conversion (RGB to YCbCr) | 100141 Bytes |  $\sim$33% Reduction |  |
GZip Compression | 12138 Bytes |  $\sim$4% Reduction |  |

##  4 A Performance Comparison

This section compares the performance of AVIS Engine with other available AV
simulation software. Note that all tests are performed on a laptop with the
specifications listed in the Table 3 . In particular, these tests compare the
AVIS engine to the CARLA simulator, a high-quality simulator for autonomous
vehicles. Both simulators provide a pedestrian system and a traffic system,
but to test performance, these features were not used.

Table 3: Test system specifications Configuration | Detail
---|---
CPU | Intel® Core™ i5-5257U CPU @ 2.70GHz
GPU | Intel® Iris Graphics 6100 1536 MB
RAM | 8 GB
Screen Resolution | 2560 x 1600
Simulator Resolution | 1920 x 1080

In the tests, no other cars or special objects are spawned in the scene, so
there is only simulator running with the main vehicle spawned and no client
connected.

Table 4 shows the results of the multiple tasks performed with each simulator.
The highest graphics setting in AVIS Engine is "Ultra" and the equivalent
setting in CARLA is "Epic". The lowest setting in AVIS Engine is called
"Fastest" and in CARLA it is called "Low". As the experiments show, the AVIS
Engine simulator is capable of running at high FPS even on a low-end laptop,
which can be useful when it comes to using efficient simulation software for
work or a study. With the AVIS Engine it is possible to achieve over 110 FPS
on a dedicated NVIDIA GeForce 1070 graphics card.

Table 4: Performance results from each simulator in different settings Simulator | FPS in Lowest Settings | FPS in Highest Settings
---|---|---
| Avg. | Top | Avg. | Top
AVIS Engine | 57 | 65 | 42 | 50
CARLA | 7 | 10 | 2 | 5

##  5 Conclusion

We have studied one state-of-the-art techniques that powers the development
autonomous vehicles which is simulation, and its necessity. We have also
described the AVIS Engine Simulator and its capabilities, features, and
capabilities, and analyzed the simulator’s good performance. In addition, we
have unpacked the ideas behind its fast networking, with the results clearly
showing how the techniques described in the paper, such as compression
algorithms and color-space conversions, have helped to reduce the data load
and thus enormously optimize the performance of the software.

## Acknowledgment

This work was supported by Bio-Inspired System Design Lab (BInSDeLa) at
Amirkabir University of Technology.

## References

  * [1] Matlab website (2022), https://www.mathworks.com/products/matlab.html
  * [2] Torcs - the open racing car simulator (2022), https://sourceforge.net/projects/torcs/
  * [3] Unity real-time development platform (2022), https://unity.com/
  * [4] Why ros? (2022), https://www.ros.org/blog/why-ros/
  * [5] Apollo: Apollo game engine based simulator - an autonomous vehicle simulator built with unity (2022), https://developer.apollo.auto/gamesim.html
  * [6] Best, A., Narang, S., Pasqualin, L., Barber, D., Manocha, D.: Autonovi-sim: Autonomous vehicle simulation platform with weather, sensing, and traffic control. In: 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW). pp. 1161–11618 (2018). https://doi.org/10.1109/CVPRW.2018.00152
  * [7] Billaut, V., de Rochemonteix, M., Thibault, M.: Colorunet: A convolutional classification approach to colorization. arXiv preprint arXiv:1811.03120 (2018)
  * [8] Burkhard, H.D., Duhaut, D., Fujita, M., Lima, P., Murphy, R., Rojas, R.: The road to robocup 2050. IEEE Robotics & Automation Magazine 9(2), 31–38 (2002)
  * [9] Cai, P., Lee, Y., Luo, Y., Hsu, D.: Summit: A simulator for urban driving in massive mixed traffic. In: 2020 IEEE International Conference on Robotics and Automation (ICRA). pp. 4023–4029 (2020). https://doi.org/10.1109/ICRA40945.2020.9197228
  * [10] Campbell, S., O’ Mahony, N., Krpalkova, L., Riordan, D., Walsh, J., Murphy, A., Ryan, C.: Sensor technology in autonomous vehicles : A review. pp. 1–4 (06 2018). https://doi.org/10.1109/ISSC.2018.8585340
  * [11] Dosovitskiy, A., Ros, G., Codevilla, F., Lopez, A., Koltun, V.: CARLA: An open urban driving simulator. In: Proceedings of the 1st Annual Conference on Robot Learning. pp. 1–16 (2017)
  * [12] Gerndt, R., Seifert, D., Baltes, J.H., Sadeghnejad, S., Behnke, S.: Humanoid robots in soccer: Robots versus humans in robocup 2050. IEEE Robotics & Automation Magazine 22(3), 147–154 (2015)
  * [13] Guerra, W., Tal, E., Murali, V., Ryou, G., Karaman, S.: Flightgoggles: Photorealistic sensor simulation for perception-driven robotics using photogrammetry and virtual reality. In: 2019 IEEE/RSJ International Conference on Intelligent Robots and Systems, IROS 2019, Macau, SAR, China, November 3-8, 2019. pp. 6941–6948. IEEE (2019). https://doi.org/10.1109/IROS40897.2019.8968116, https://doi.org/10.1109/IROS40897.2019.8968116
  * [14] Ibraheem, N., Hasan, M., Khan, R.Z.: Understanding color models: A review. ARPN Journal of Science and Technology 2 (01 2012)
  * [15] Kam, H.R., Lee, S.H., Park, T., Kim, C.H.: Rviz: a toolkit for real domain data visualization. Telecommunication Systems 60, 337–345 (2015)
  * [16] Kitano, H., Asada, M., Kuniyoshi, Y., Noda, I., Osawa, E., Matsubara, H.: Robocup: A challenge problem for ai. AI magazine 18(1), 73–73 (1997)
  * [17] Lima Azevedo, C., Oh, S., Marczuk, K., Soh, H., Basak, K., Toledo, T., Peh, L.S., Ben-Akiva, M., Deshmukh, N.: Simmobility short-term: An integrated microscopic mobility simulator (01 2017)
  * [18] Luo, J., Hubaux, J.P.: A survey of inter-vehicle communication (2004), http://infoscience.epfl.ch/record/28039
  * [19] Álvarez López, P., Behrisch, M., Bieker-Walz, L., Erdmann, J., Flötteröd, Y.P., Hilbrich, R., Lücken, L., Rummel, J., Wagner, P., Wießner, E.: Microscopic traffic simulation using sumo (11 2018). https://doi.org/10.1109/ITSC.2018.8569938
  * [20] Maria, A.: Introduction to modeling and simulation. In: Proceedings of the 29th Conference on Winter Simulation. p. 7–13. WSC ’97, IEEE Computer Society, USA (1997). https://doi.org/10.1145/268437.268440, https://doi.org/10.1145/268437.268440
  * [21] Nimalsiri, N.I., Mediwaththe, C.P., Ratnam, E.L., Shaw, M., Smith, D.B., Halgamuge, S.K.: A survey of algorithms for distributed charging control of electric vehicles in smart grid. IEEE Transactions on Intelligent Transportation Systems 21(11), 4497–4515 (2020). https://doi.org/10.1109/TITS.2019.2943620
  * [22] Podpora, M., Korba´s, G., Kawala-Janik, A.: Yuv vs rgb – choosing a color space for human-machine interaction. Annals of Computer Science and Information Systems Vol. 3 (09 2014). https://doi.org/10.15439/2014F206
  * [23] Quigley, M., Conley, K., Gerkey, B., Faust, J., Foote, T., Leibs, J., Wheeler, R., Ng, A.: Ros: an open-source robot operating system. vol. 3 (01 2009)
  * [24] Raj, T., Hashim, F.H., Huddin, A.B., Ibrahim, M.F., Hussain, A.: A survey on lidar scanning mechanisms. Electronics 9(5), 741 (2020)
  * [25] Sadeghnejad, S., Khadivar, F., Abdollahi, E., Moradi, H., Farahmand, F., Sadr Hosseini, S.M., Vossoughi, G.: A validation study of a virtual-based haptic system for endoscopic sinus surgery training. The International Journal of Medical Robotics and Computer Assisted Surgery 15(6), e2039 (2019)
  * [26] Shah, S., Dey, D., Lovett, C., Kapoor, A.: Airsim: High-fidelity visual and physical simulation for autonomous vehicles. In: Field and Service Robotics (2017), https://arxiv.org/abs/1705.05065
  * [27] Yang, Y., Yuhua, P., Zhaoguang, L.: A fast algorithm for ycbcr to rgb conversion. IEEE Transactions on Consumer Electronics 53(4), 1490–1493 (2007). https://doi.org/10.1109/TCE.2007.4429242