On the Local Structure of Stable Clustering Instances
Vincent Cohen-Addad∗
University of Copenhagen

arXiv:1701.08423v3 [] 10 Aug 2017

Chris Schwiegelshohn

†

Sapienza University of Rome

Abstract
We study the classic k-median and k-means clustering objectives in the beyond-worst-case
scenario. We consider three well-studied notions of structured data that aim at characterizing
real-world inputs:
• Distribution Stability (introduced by Awasthi, Blum, and Sheffet, FOCS 2010)
• Spectral Separability (introduced by Kumar and Kannan, FOCS 2010)
• Perturbation Resilience (introduced by Bilu and Linial, ICS 2010)
We prove structural results showing that inputs satisfying at least one of the conditions are
inherently “local”. Namely, for any such input, any local optimum is close both in term of
structure and in term of objective value to the global optima.
As a corollary we obtain that the widely-used Local Search algorithm has strong performance
guarantees for both the tasks of recovering the underlying optimal clustering and obtaining a
clustering of small cost. This is a significant step toward understanding the success of local
search heuristics in clustering applications.

∗
†

vcohenad@gmail.com
chris.schwiegelshohn@tu-dortmund.de

1

Introduction

Clustering is a fundamental, routinely-used approach to extract information from datasets. Given
a dataset and the most important features of the data, a clustering is a partition of the data such
that data elements in the same part have common features. The problem of computing a clustering
has received a considerable amount of attention in both practice and theory.
The variety of contexts in which clustering problems arise makes the problem of computing a
“good” clustering hard to define formally. From a theoretician’s perspective, clustering problems
are often modeled by an objective function we wish to optimize (e.g., the famous k-median or
k-means objective functions). This modeling step is both needed and crucial since it provides a
framework to quantitatively compare algorithms. Unfortunately, the most popular objectives for
clustering, like the k-median and k-means objectives, are hard to approximate, even when restricted
to Euclidean spaces.
This view is generally not shared by practitioners. Indeed, clustering is often used as a preprocessing step to simplify and speed up subsequent analysis, even if this analysis admits polynomial
time algorithms. If the clustering itself is of independent interest, there are many heuristics with
good running times and results on real-world inputs.
This induces a gap between theory and practice. On the one hand, the algorithms that are
efficient in practice cannot be proven to achieve good approximation to the k-median and k-means
objectives in the worst-case. Since approximation ratios are one of the main methods to evaluate
algorithms, theory predicts that determining a good clustering is a difficult task. On the other
hand, the best theoretical algorithms turn out to be noncompetitive in applications because they
are designed to handle “unrealistically” hard instances with little importance for practitioners. To
bridge the gap between theory and practice, it is necessary to go beyond the worst-case analysis by,
for example, characterizing and focusing on inputs that arise in practice.

1.1

Real-world Inputs

Several approaches have been proposed to bridge the gap between theory and practice. For example,
researchers have considered the average-case scenario (e.g., [26]) where the running time of an
algorithm is analyzed with respect to some probability distribution over the set of all inputs.
Smooth analysis (e.g., [90]) is another celebrated approach that analyzes the running time of an
algorithm with respect to worst-case inputs subject to small random perturbations.
Another successful approach, the one we take in this paper, consists in focusing on structured
inputs. In a seminal paper, Ostrovsky, Rabani, Schulman, and Swamy [85] introduced the idea
that inputs that come from practice induce a ground-truth or a meaningful clustering. They argued
that an input I contains a meaningful clustering into k clusters if the optimal k-median cost of
a clustering using k centers, say OPTk (I), is much smaller than the optimal cost of a clustering
using k − 1 centers OPTk−1 (I). This is also motivated by the elbow method 1 (see Section 7 for
more details) used by practitioners to define the number of clusters. More formally, an instance I
of k-median or k-means satisfies the α-ORSS property if OPTk (I)/OPTk−1 (I) ≤ α.
α-ORSS inputs exhibit interesting properties. The popular k-means++ algorithm (also known
as the D2 -sampling technique) achieves an O(1)-approximation for these inputs2 . The condition is
also robust with respect to noisy perturbations of the data set. ORSS-stability also implies several
1

The elbow-method consists in running an (approximation) algorithm for an incrementally increasing number of
clusters until the cost drops significantly.
2
For worst-case inputs, the k-means++ achieves an O(log k)-approximation ratio [9, 31, 66, 85].

1

other conditions aiming to capture well-clusterable instances. Thus, the inputs satisfying the ORSS
property arguably share some properties with the real-world inputs. In this paper, we also provide
experimental results supporting this claim, see Appendix C.
These results have opened new research directions and raised several questions. For example:
• Is it possible to obtain similar results for more general classes of inputs?
• How does the parameter α impact the approximation guarantee and running time?
• Is it possible to prove good performance guarantees for other popular heuristics?
• How close to the “ground-truth” clustering are the approximate clusterings?
We now review the most relevant work in connection to the above open questions, see Sections 2
for other related work.
Distribution Stability (Def. 4.1) Awasthi, Blum and Sheffet [12] have tackled the first two
questions by introducing the notion of distribution stable instances. Distribution stable instances
are a generalization of the ORSS instances (in other words, any instance satisfying the ORSS
property is distribution stable). They also introduced a new algorithm tailored for distribution
stable instances that achieves a (1 + ε)-approximation for α-ORSS inputs (and more generally αdistribution stable instances) in time nO(1/εα) . This was the first algorithm whose approximation
guarantee was independent from the parameter α for α-ORSS inputs.
Spectral Separability (Def. 6.1) Kumar and Kannan [74] tackled the first and third questions
by introducing the proximity condition3 . This condition also generalizes the ORSS condition. It
is motivated by the goal of learning a distribution mixture in a d-dimensional Euclidean space.
Quoting [74], the message of their paper can loosely be stated as:
If the projection of any data point onto the line joining its cluster center to any other
cluster center is γk times standard deviations closer to its own center than the other
center, then we can cluster correctly in polynomial time.
In addition, they have made a significant step toward understanding the success of the classic
k-means by showing that it achieves a 1 + O(1/γ)-approximation for instances that satisfy the
proximity condition.
Perturbation Resilience (Def. 5.1) In a seminal work, Bilu and Linial [29] introduced a new
condition to capture real-world instances. They argue that the optimal solution of a real-world
instance is often much better than any other solution and so, a slight perturbation of the instance
does not lead to a different optimal solution. Perturbation-resilient instances have been studied in
various contexts (see e.g., [13, 16, 20, 21, 27, 76]). For clustering problems, an instance is said to
be α-perturbation resilient if an adversary can change the distances between pairs of elements by a
factor at most α and the optimal solution remains the same. Recently, Angelidakis, Makarychev,
and Makarychev [80] have given a polynomial-time algorithm for solving 2-perturbation-resilient instances4 . Balcan and Liang [21] have tackled
√ the third question by showing that a classic algorithm
for hierarchical clustering can solve 1 + 2-perturbation-resilient instances. This very interesting
3

In this paper, we work with a slightly more general condition called spectral separability but the motivations
behind the two conditions are similar.
4
We note that it is NP-hard to recover the optimal clustering of a < 2-perturbation-resilient instance [27].

2

result leaves open the question as whether classic algorithms for (“flat”) clustering could also be
proven to be efficient for perturbation-resilient instances.
Main Open Questions Previous work has made important steps toward bridging the gap between theory and practice for clustering problems. However, we still do not have a complete
understanding of the properties of “well-structured” inputs, nor do we know why the algorithms
used in practice perform so well. Some of the most important open questions are the following:
• Do the different definitions of well-structured input have common properties?
• Do heuristics used in practice have strong approximation ratios for well-structured inputs?
• Do heuristics used in practice recover the “ground-truth” clustering on well-structured inputs?

1.2

Our Results: A unified approach via Local Search

We make a significant step toward answering the above open questions. We show that the classic
Local Search heuristic (see Algorithm 1), that has found widespread application in practice (see
Section 2), achieves good approximation guarantees for distribution-stable, spectrally-separable,
and perturbation-resilient instances (see Theorems 4.2, 5.2, 6.2).
More concretely, we show that Local Search is a polynomial-time approximation scheme (PTAS)
for both distribution-stable and spectrally-separable5 instances. In the case of distribution stability,
we also answer the above open question by showing that most of the structure of the optimal
underlying clustering is recovered by the algorithm. Furthermore, our results hold even when only
a δ fraction (for any constant δ > 0) of the points of each optimal cluster satisfies the β-distributionstability property.
For γ-perturbation-resilient instances, we show that if γ > 3 then any solution is the optimal
solution if it cannot be improved by adding or removing 2γ centers. We also show that the analysis
is essentially tight.
These results show that well-structured inputs have the property that the local optima are close
both qualitatively (in terms of structure) and quantitatively (in terms of objective value) to the
global “ground-truth” optimum.
These results make a significant step toward explaining the success of Local Search approaches
for solving clustering problems in practice.

Algorithm 1 Local Search(ε) for k-Median and k-Means
1:
2:
3:
4:
5:
6:
7:
8:

5

Input: A, F, cost, k
Parameter: ε
S ← Arbitrary subset of F of cardinality at most k.
while ∃ S 0 s.t. |S 0 | ≤ k and |S − S 0 | + |S 0 − S| ≤ 2/ε and cost(S 0 ) ≤ (1 − ε/n) cost(S)
do
S ← S0
end while
Output: S

Assuming a standard preprocessing step consisting of a projection onto a subspace of lower dimension.

3

1.3

Organization of the Paper

Section 2 provides a more detailed review of previous work on worst-case approximation algorithms
and Local Search. Further comments on stability conditions not covered in the introduction can
be found in Section 7 at the end of the paper. Section 3 introduces preliminaries and notation.
Section 4 is dedicated to distribution-stable instances, Section 5 to perturbation-resilient instances,
and Section 6 to spectrally-separated instances. All the missing proofs can be found in the appendix.

2

Related Work

Worst-Case Hardness The problems we study are NP-hard: k-median and k-means are already
NP-hard in the Euclidean plane (see Meggido and Supowit [83], Mahajan et al. [79], and Dasgupta
and Freud [43]). In terms of hardness of approximation, both problems are APX-hard, even in the
Euclidean setting when both k and d are part of the input (see Guha and Khuller [56], Jain et
al. [64], Guruswami et al. [59], and Awasthi et al. [14]). On the positive side, constant factor
approximations are known in metric space for both k-median and k-means (see [3, 33, 77, 65, 84]).
For Euclidean spaces we have a PTAS for both problems, either assuming d fixed and k arbitrary [7,
37, 52, 62, 63, 72], or assuming k fixed and d arbitrary [48, 75].
Local Search Local Search is an all-purpose heuristic that may be applied to any problem,
see Aarts and Lenstra [1] for a general introduction. For clustering, there exists a large body
of bicriteria approximations for k-median and k-means [23, 34, 38, 73]. Arya et al. [11] showed
that Local Search with a neighborhood size of 1/ε gives a 3 + 2ε approximation to k-median, see
also [58]. Kanungo et al. [70] proved an approximation ratio of 9 + ε for k-means clustering by
Local Search, which was until very recently [3] the best known algorithm with a polynomial running
time in metric and Euclidean spaces.6 Recently, Local Search with an appropriate neighborhood
size was shown to be a PTAS for k-means and k-median in certain restricted metrics including
constant dimensional Euclidean space [37, 52]. Due to its simplicity, Local Search is also a popular
subroutine for clustering tasks in various more specialized computational models [24, 30, 57]. For
more theoretical clustering papers using Local Search, we refer to [39, 45, 53, 60, 95].
Local Search is also often used for clustering in more applied areas of computer science (e.g., [92,
54, 4, 61]). Indeed, the use of Local Search with a neighborhood of size 1 for clustering was first
proposed by Tüzün and Burke [93], see also Ghosh [55] for a more efficient version of the same
approach. Due the ease by which it may be implemented, Local Search has become one of the
most commonly used heuristics for clustering and facility location, see Ardjmand [5]. Nevertheless,
high running times is one of the biggest drawbacks of Local Search compared to other approaches,
though a number of papers have engineered it to become surprisingly competitive, see Frahling and
Sohler [51], Kanungo et al. [69], and Sun [91].

3

Definitions and Notations

The problem The problem we consider in this work is the following slightly more general version
of the k-means and k-median problems.
6
They combined Local Search with techniques from Matousek [81] for k-means clustering in Euclidean spaces. The
running time of the algorithm as stated incurs an additional factor of ε−d due to the use of Matousek’s approximate
centroid set. Using standard techniques (see e.g. Section B of this paper), a fully polynomial running time in n, d,
and k is also possible without sacrificing approximation guarantees.

4

Definition 3.1 (k-Clustering). Let A be a set of clients, F a set of centers, both lying in a metric
space (X , dist), cost a function A × F → R+ , and k a non-negative integer. The k-clustering
problem asks for a subset S of F , of cardinality at most k, that minimizes
X
cost(S) =
min cost(x, c).
x∈A

c∈S

The clustering of A induced by S is the partition of A into subsets C = {C1 , . . . Ck } such that
Ci = {x ∈ A | ci = argmin cost(x, c)} (breaking ties arbitrarily).
c∈S

The well known k-median and k-means problems correspond to the special cases cost(a, c) =
dist(a, c) and cost(a, c) = dist(a, c)2 respectively. Throughout the rest of this paper, let OPT denote
the value of an optimal solution. To give slightly simpler proofs for β-distribution-stable and αperturbation-resilient instances, we will assume that cost(a, b) = dist(a, b). If cost(a, b) = dist(a, b)p ,
then α depends exponentially on the p for perturbation resilience. For distribution stability, we still
have a PTAS by introducing a dependency in 1/εO(p) in the neighborhood size of the algorithm.
The analysis is unchanged save for various applications of the following lemma at different steps of
the proof.
Lemma 3.2. Let p ≥ 0 and 1/2 > ε > 0. For any a, b, c ∈ A ∪ F , we have cost(a, b) ≤ (1 +
ε)p cost(a, c) + cost(c, b)(1 + 1/ε)p .

4

Distribution Stability

We work with the notion of β, δ-distribution stability which generalizes β-distribution stability.
This extends our result to datasets that exhibit a slightly weaker structure than the β-distribution
stability. Namely, the β, δ-distribution stability only requires that for each cluster of the optimal
solution, most of the points satisfy the β-distribution stability condition.
Definition 4.1 ((β, δ)-Distribution Stability). Let (A, F, cost, k) be an instance of k-clustering
where A ∪ F lie in a metric space and let S ∗ = {c∗1 , . . . , c∗k } ⊆ F be a set of centers and C ∗ =
{C1∗ , . . . , Ck∗ } be the clustering induced by S ∗ . Further, let β > 0 and 0 ≤ δ ≤ 1. Then the
pair (A, F, cost, k), (C ∗ , S ∗ ) is a (β, δ)-distribution stable instance if, for any i, there exists a set
∆i ⊆ Ci∗ such that |∆i | ≥ (1 − δ)|Ci∗ | and for any x ∈ ∆i , for any j 6= i,
cost(x, c∗j ) ≥ β

OPT
,
|Cj∗ |

where cost(x, c∗j ) is the cost of assigning x to c∗j .
For any instance (A, F, cost, k) that is (β, δ)-distribution stable, we refer to (C ∗ , S ∗ ) as a (β, δ)clustering of the instance. We show the following theorem for the k-median problem. For the
k-clustering problem with parameter p, the constant η becomes a function of p.
Theorem 4.2. Let p > 0, β > 0, and ε < min(1 − δ, 1/3). For a (β, δ)-stable instance with
(β, δ) clustering (C ∗ , S ∗ ) and an absolute constant η, the cost of the solution output by Local
Search(4ε−3 β −1 + O(ε−2 β −1 )) (Algorithm 1) is at most (1 + ηε)cost(C ∗ ).
Moreover, let L = {L1 , . . . , Lk } denote the clusters of the solution output by Local Search(4ε−3 β −1 +
O(ε−2 β −1 )). If δ = 0 (i.e.: the instance is simply β-distribution-stable), there exists a bijection
φ : L 7→ C ∗ such that for at least m = k − O(ε−3 β −1 ) clusters L01 , . . . , L0m ⊆ L, the following two
statements hold.
5

IRεi
2

IRεi

c∗i

β·

OPT
|Ci |

c∗j
p ∈ ∆j

L(i)

Ci∗

Cj∗
2

Figure 1: Example of aScluster Ci∗ 6∈ Z ∗ . An important fraction of the points in IRεi are served by
L(i) and few points in j6=i ∆j are served by L(i).

2

• At least a (1 − ε) fraction of IRεi ∩ Ci∗ are served by a unique center L(i) in solution L.
S
2
• The total number of clients p ∈ j6=i Cj∗ served by L(i) in L is at most ε|IRεi ∩ Ci∗ |.
We first give a high-level description of the analysis. Assume for simplicity that all the optimal
clusters cost less than an ε3 fraction of the total cost of the optimal solution. Combining this
assumption with the β-distribution-stability property, one can show that the centers and points
close to the center are far away from each other. Thus, guided by the objective function, the local
search algorithm identifies most of these centers. In addition, we can show that for most of these
good centers the corresponding cluster in the local solution is very similar to the optimal cluster
(see Figure 1). In total, only very few clusters (a function of ε and β) of the optimal solution are not
present in the local solution. We conclude our proof by using local optimality. Our proof includes
a few ingredients from [12] such as the notion of inner-ring (we work with a slightly more general
definition) and distinguishes between cheap and expensive clusters. Nevertheless our analysis is
slightly stronger as we consider a significantly weaker stability condition and can not only analyze
the cost of the solution of the algorithm, but also the structure of its clusters.
Throughout this section, we consider a set of centers S ∗ = {c∗1 , . . . , c∗k } whose induced clustering
is C ∗ = {C1∗ , . . . , Ck∗ } and such that the instance is (β, δ)-stable with respect (C ∗ , S ∗ ). We denote
P P
by clusters the parts of a partition C ∗ = {C1∗ , . . . , Ck∗ }. Let cost(C ∗ ) = ki=1 x∈C ∗ cost(x, c∗i ).
i
Moreover, for any cluster Ci∗ , for any client x ∈ Ci∗ , denote by gx the cost of client x in solution
C ∗ : gx = cost(x, c∗i ) = dist(x, c∗i ) since we consider the k-median problem. Let L denote the
output of LocalSearch(β −1 ε−3 ) and
P lx the cost induced by client x in solution L, namely lx =
min`∈L cost(x, `), and cost(L) = x∈A lx . The following definition is a generalization of the innerring definition of [12].

6

Definition 4.3. For any ε0 , we define the inner ring of cluster i, IRεi 0 , as the set of x ∈ A ∪ F
such that dist(x, c∗i ) ≤ ε0 βOPT/|Ci∗ |.
P
We say that cluster i is cheap if x∈C ∗ gx ≤ ε3 βOPT, and expensive otherwise. We aim at
i
proving the following structural lemma.
Lemma 4.4. There exists a set of clusters Z ∗ ⊆ C ∗ of size at most 2ε−3 β −1 + O(ε−2 β −1 ) such
that for any cluster Ci∗ ∈ C ∗ − Z ∗ , we have the following properties
1. Ci∗ is cheap.
2

2. At least a (1 − ε) fraction of IRεi ∩ Ci∗ are served by a unique center L(i) in solution L.
S
2
3. The total number of clients p ∈ j6=i ∆j served by L(i) in L is at most ε|IRεi ∩ Ci∗ |.
See Fig 1 for a typical cluster of C ∗ − Z ∗ . We start with the following lemma which generalizes
Fact 4.1 in [12].
Lemma 4.5. Let Ci∗ be a cheap cluster. For any ε0 , we have |IRεi 0 ∩ Ci∗ | > (1 − ε3 /ε0 )|Ci∗ |.
We then prove that the inner rings of cheap clusters are disjoint for δ +
Lemma 4.6. Let δ +

ε3
ε0

ε3
ε0

< 1 and ε0 < 13 .

< 1 and ε0 < 13 . If Ci∗ 6= Cj∗ are cheap clusters, then IRεi 0 ∩ IRεj 0 = ∅.

For each cheap cluster Ci∗ , let L(i) denote a center of L that belongs to IRεi if there exists
exactly such center and remain undefined otherwise. By Lemma 4.6, L(i) 6= L(j) for i 6= j.
Lemma 4.7. Let ε < 31 . Let C ∗ − Z1 denote the set of clusters Ci∗ that are cheap, such that L(i) is
2
2
defined and such that at least (1 − ε)|IRεi ∩ Ci∗ | clients of IRεi ∩ Ci∗ are served in L by L(i). Then
|Z1 | ≤ (2ε−3 + 11.25 · ε−2 + 22.5 · ε−1 )β −1 .
Proof. There are five different types of clusters in C ∗ :
1. k1 expensive clusters
2. k2 cheap clusters with no center of L belonging to IRεi
3. k3 cheap clusters with at least two centers of L belonging to IRεi
2

2

4. k4 cheap clusters with L(i) being defined and less than (1 − ε)|IRεi ∩ Ci∗ | clients of IRεi ∩ Ci∗
are served in L by L(i)
2

2

5. k5 cheap clusters with L(i) being defined and at least (1 − ε)|IRεi ∩ Ci∗ | clients of IRεi ∩ Ci∗
are served in L by L(i)
The definition of cheap clusters immediately yields k1 ≤ ε−3 β −1 .
Since L and C ∗ both have k clusters and the inner rings of cheap clusters are disjoint (Lemma 4.6),
we have c1 k1 + c3 k3 + k4 + k5 = k1 + k2 + k3 + k4 + k5 = |Z1 | + k5 = k with c1 ≥ 0 and c3 ≥ 2
resulting in k3 ≤ (c3 − 1)k3 = (1 − c1 )k1 + k2 ≤ k1 + k2 .
Before bounding k2 and k4 , we discuss the impact of a cheap cluster Ci∗ with at least a p fraction
2
of the clients of IRεi ∩ Ci∗ being served in L by some centers that are not in IRεi . By the triangular
inequality, the cost for any client x of this p fraction is at least (ε − ε2 )βcost(C ∗ )/|Ci∗ |. Then the

7

2

total cost of all clients of this p fraction in L is at least p|IRεi ∩ Ci∗ |(1 − ε)εβcost(C ∗ )/|Ci∗ |. By
2
Lemma 4.5, substituting |IRεi ∩ Ci∗ | yields for this total cost
2

p|IRεi ∩ Ci∗ |(1 − ε)εβ

cost(C ∗ )
cost(C ∗ )
≥ p(1 − ε)2 |Ci∗ |εβ
= p(1 − ε)2 εβcost(C ∗ ).
∗
|Ci |
|Ci∗ |

To determine k2 , we must use p = 1 while we have p > ε for k4 . Therefore, the total costs of all
clients of the k2 and the k4 clusters in L are at least k2 (1−ε)2 εβcost(C ∗ ) and k4 (1−ε)2 ε2 βcost(C ∗ ),
respectively.
Now, since cost(L) ≤ 5OPT ≤ 5cost(C ∗ ), we have (k2 + k4 ε)εβ ≤ 5/(1 − ε)2 ≤ 45/4.
Therefore, we have |Z1 | = k1 +k2 +k3 +k4 ≤ 2k1 +2k2 +k4 ≤ (2ε−3 +11.25·ε−2 +22.5·ε−1 )β −1 .
We continue with the following lemma, whose proof relies on similar arguments.
Lemma 4.8. There exists a set Z2 ⊆ C ∗ − Z1 of S
size at most 11.25ε−1 β −1 such that for any cluster
∗
∗
Cj ∈ C − Z2 , the total number of clients x ∈ i6=j ∆i , that are served by L(j) in L is at most
2

ε|IRεi ∩ Ci∗ |.
Therefore, the proof of Lemma 4.4 follows from combining Lemmas
S 4.7 and 4.8.
We now turn to the analysis of the cost of L. Let C(Z ∗ ) = C ∗ ∈Z ∗ Ci∗ . For any cluster
i

2

Ci∗ ∈ C ∗ − Z ∗ , let L(i) be the unique center of L that serves at least (1 − ε)|IRεi ∩ Ci∗ | > (1 − ε)2 |Ci |
S
2
b to be the set
clients of IRεi ∩ Ci∗ , see Lemmas 4.4 and 4.5. Let Lb = C ∗ ∈C ∗ −Z ∗ L(i) and define A
i
b Finally, let A(L(i)) be the set of clients
of clients that are served in solution L by centers of L.
b
that are served by L(i) in solution L. Observe that the A(L(i)) partition A.
Lemma 4.9. We have
X

−ε · cost(L)/n +

X

lx ≤

∗)
b
x∈A−C(Z

gx +

∗)
b
x∈A−C(Z

2ε
· (cost(C ∗ ) + cost(L)).
(1 − ε)2

Proof. Consider the following mixed solution M = Lb ∪ {c∗i | Ci∗ ∈ Z ∗ }. We start by bounding the
b the center that serves it in L belongs to M. Thus its cost in M
cost of M. For any client x ∈ A,
is at most lx . Now, for any client x ∈ C(Z ∗ ), the center that serves it in C ∗ is in M, so its cost in
M is at most gx .
b ∪ C(Z ∗ )). Consider such a client x and
Finally, we evaluate the cost of the clients in A − (A
let Ci∗ be the cluster it belongs to in solution C ∗ . Since Ci∗ ∈ C ∗ − Z ∗ , L(i) is defined and we have
L(i) ∈ Lb ⊆ M. Hence, the cost of x in M is at most cost(x, L(i)). Observe that by the triangular
inequality, cost(x, L(i)) ≤ cost(x, c∗i ) + cost(c∗i , L(i)) = gx + cost(c∗i , L(i)).
2
Now consider a client x0 ∈ IRεi ∩Ci∗ ∩A(L(i)). By the triangular inequality, we have cost(c∗i , L(i)) ≤
cost(c∗i , x0 ) + cost(x0 , L(i)) = gx0 + lx0 . Hence,
cost(c∗i , L(i)) ≤

1
2
|IRεi

∩

Ci∗

X

∩ A(L(i))|

(gx0 + lx0 ).

2
x0 ∈IRεi ∩Ci∗ ∩A(L(i))

b to L(i) induces a cost of at most
It follows that assigning the clients of Ci∗ ∩ (A − A)
X
b
x∈Ci∗ ∩(A−A)

gx +

b
|Ci∗ ∩ (A − A)|

X

2

|IRεi ∩ Ci∗ ∩ A(L(i))|
8

2
x0 ∈IRεi ∩Ci∗ ∩A(L(i))

(gx0 + lx0 ).

2

2

2

b ≤
Due to Lemma 4.4, we have |IRεi ∩Ci∗ ∩A(L(i))| ≥ (1−ε)·|IRεi ∩Ci∗ | and |(IRεi ∩Ci∗ )∩(A− A)|
ε2
ε2
ε2
ε2
∗
∗
∗
∗
∗
b
ε · |IRi ∩ Ci |. Further, |(Ci − IRi ) ∩ (A − A)| ≤ |(Ci − IRi )| = |Ci | − |IRi ∩ Ci |. Combining
these three bounds, we have
b
|Ci∗ ∩ (A − A)|
2

|IRεi ∩ Ci∗ ∩ A(L(i))|

2
b + |(C ∗ ∩ IRε2 ) ∩ (A − A)|
b
|(Ci∗ − IRεi ) ∩ (A − A)|
i
i

=

2

|IRεi ∩ Ci∗ ∩ A(L(i))|
2

|Ci∗ | − (1 − ε)|IRεi ∩ Ci∗ |

≤

2

|Ci∗ |

=

2

(1 − ε) · |IRεi ∩ Ci∗ |
(1 − ε) · |IRεi ∩ Ci∗ |
|Ci∗ |
2ε
2ε − ε2
<
,
−
1
≤
∗
(1 − ε)2 · |Ci |
(1 − ε)2
(1 − ε)2

≤

−1
(1)

where the inequality in (1) follows from Lemma 4.5.
Summing over all clusters Ci∗ ∈ C ∗ − Z ∗ , we obtain that the cost in M for the clients in
b ∩ C ∗ is less than
(A − A)
i
X
2ε
· (cost(C ∗ ) + cost(L)).
gx +
(1 − ε)2
∗ ))
b
c∈A−(A∪C(Z

By Lemmas 4.7 and 4.8, we have |M − L| + |L − M| = 2 · |Z ∗ | ≤ (4ε−3 + O(ε−2 ))β −1 . By
selecting the neighborhood size of Local Search (Algorithm 1) to be greater than this value, we
have (1 − ε/n) · cost(L) ≤ cost(M). Therefore, combining the above observations, we have
X
X
X
2ε
ε
lx +
gx +
gx +
· (cost(C ∗ ) + cost(L)).
(1 − ) · cost(L) ≤
2
n
(1
−
ε)
∗
∗)
b
x∈A−C(Z

x∈C(Z )

∗ ))
b
x∈A−(A∪C(Z

By simple transformations, we then obtain
X
ε
− · cost(L) +
lx ≤
n
∗)
b
x∈A−(A)∪C(Z

X

gx +

∗)
b
x∈A−(A)∪C(Z

2ε
· (cost(C ∗ ) + cost(L)).
(1 − ε)2

b − C(Z ∗ ). For any cluster
We now turn to evaluate the cost for the clients that are in A
∗
∗
∗
∈ C − C(Z ) and for any x ∈ Ci − A(L(i)) define Reassign(x) to be the cost of x with
respect to the center in L(i). Note that there exists only one center of L in IRiε for any cluster
Ci∗ ∈ C ∗ − C(Z ∗ ). Before going deeper in the analysis, we need the following lemma.

Ci∗

Lemma 4.10. For any Ci∗ ∈ C ∗ − C(Z ∗ ), we have
X
X
Reassign(x) ≤
x∈Ci∗ −A(L(i))

x∈Ci∗ −A(L(i))

gx +

X
2ε
(lx + gx ).
2
(1 − ε)
∗
x∈Ci

We now partition the clients of cluster Ci∗ ∈ C ∗ − Z ∗ . For any i, let Bi be the set of clients of
that are served in solution L by a center L(j) for some j 6= i and Cj∗ ∈ C ∗ − Z ∗ . Moreover, let
S
b − S Dj .
Di = (A(L(i)) ∩ ( j6=i Bj )). Finally, define Ei = (Ci∗ ∩ A)
j6=i

Ci∗

Lemma 4.11. Let Ci∗ be a cluster in C ∗ − Z ∗ . Define the solution Mi = L − {L(i)} ∪ {c∗i } and
denote by mix the cost of client x in solution Mi . Then
X
X
X
X
X
X
ε
mix ≤
lx +
gx +
Reassign(x) +
lx +
(
gx + lx ).
(1 − ε)
x∈A

x∈A−
(A(L(i))∪Ei )

x∈Ei

x∈Di

x∈A(L(i))−
(Ei ∪Di )

9

x∈Ei

We can thus prove the following lemma, which concludes the proof.
Lemma 4.12. We have
X

−ε · cost(L) +

lx ≤

∗)
b
x∈A−C(Z

X

gx +

∗)
b
x∈A−C(Z

3ε
· (cost(L) + cost(C ∗ )).
(1 − ε)2

The proof of Theorem 4.2 follows from (1) summing the equations from Lemmas 4.9 and 4.12
and (2) Lemma 4.4. The comparison of the structure of the local solution to the structure of C ∗ is
an immediate corollary of Lemma 4.4.

5

Perturbation Resilience

We first give the definition of α-perturbation-resilient instances.
Definition 5.1. Let I = (A, F, cost, k) be an instance for the k-clustering problem. For α ≥ 1, I
is α-perturbation-resilient if there exists a unique optimal set of centers C ∗ = {c∗1 , . . . , c∗k } and for
any instance I 0 = (A, F, cost0 , k, p), such that
∀ a, b ∈ P, cost(a, b) ≤ cost0 (a, b) ≤ αcost(a, b),
the unique optimal set of centers is C ∗ = {c∗1 , . . . , c∗k }.
For ease of exposition, we assume that cost(a, b) = dist(a, b) (i.e., we work with the k-median
problem). Given solution S0 , we say that S0 is 1/ε-locally optimal if any solution S1 such that
|S0 − S1 | + |S1 − S0 | ≤ 2/ε has at least cost(S0 ).
Theorem 5.2. Let α > 3. For any instance of the k-median problem that is α-perturbationresilient, any 2(α − 3)−1 -locally optimal solution is the optimal set of centers {c∗1 , . . . , c∗k }.
Moreover, define lc to be the cost for client c in solution L and gc to be its cost in the optimal
solution C ∗ . Finally, for any sets of centers S and S0 ⊂ S, define NS (S0 ) to be the set of clients
served by a center of S0 in solution S, i.e.: NS (S0 ) = {x | ∃s ∈ S0 , dist(x, s) = mins0 ∈S dist(x, s0 )}.
The proof of Theorem 5.2 relies on the following theorem of particular interest.
Theorem 5.3 (Local-Approximation Theorem.). Let L be a 1/ε-locally optimal solution and C ∗
be any solution. Define S = L ∩ C ∗ and L̃ = L − S and C˜∗ = C ∗ − S. Then
X
X
X
X
lc +
lc ≤
gc + (3 + 2ε)
gc .
c∈NC ∗ (C˜∗ )−NL (L̃)

c∈NL (L̃)

c∈NC ∗ (C˜∗ )−NL (L̃)

c∈NL (L̃)

We first show how Theorem 5.3 allows us to prove Theorem 5.2.
Proof of Theorem 5.2. Given an instance (A, F, dist, k), we define the following instance I 0 =
(A, F, dist0 , k), where dist0 (a, b) is a distance function defined over A ∪ F that we detail below.
For each client c ∈ NL (L̃) ∪ NC ∗ (C˜∗ ), let `c be the center of L that serves it in L, for any point
p 6= `c , we define dist0 (c, p) = αdist(c, p) and dist0 (c, `c ) = dist(c, `c ). For the other clients we set
dist0 = dist. Observe that by local optimality, the clustering induced by L is {c∗1 , . . . , c∗k } if and
only if L = C ∗ . Therefore, the cost of C ∗ in instance I 0 is equal to
X
X
X
α
gc +
min(αgc , lc ) +
gc .
c∈NL (L̃)

c∈NC ∗ (C˜∗ )−NL (L̃)

c6∈NC ∗ (C˜∗ )∪NL (L̃)

10

On the other hand, the cost of L in I 0 is the same as in I. By Theorem 5.3
X

X

lc +

c∈NC ∗ (C˜∗ )−NL (L̃)

X

lc + ≤

gc + (3 +

c∈NC ∗ (C˜∗ )−NL (L̃)

c∈NL (L̃)

2(α − 3)
)
2

X

gc

c∈NL (L̃)

and by definition of S we have, for each element c ∈
/ NC ∗ (C˜∗ ) ∪ NL (L̃), lc = gc .
0
Thus the cost of L in I is at most
(3 +

2(α − 3)
)
2

X

X

gc +

c∈NC ∗ (C˜∗ )−NL (L̃)

c∈NL (L̃)

X

gc +

gc

c6∈NC ∗ (C˜∗ )∪NL (L̃)

Now, observe that for the clients in NC ∗ (C˜∗ ) − NL (L̃) = NC ∗ (C˜∗ ) ∩ NL (S), we have lc ≥ gc .
Therefore, we have that the cost of L is at most the cost of C ∗ in I 0 and so by definition of
α-perturbation-resilience, we have that the clustering {c∗1 , . . . , c∗k } is the unique optimal solution in
I 0 . Therefore L = C ∗ and the Theorem follows.
We now turn to the proof of Theorem 5.3.
Consider the following bipartite graph Γ = (L̃ ∪ C˜∗ , E) where E is defined as follows. For any
center f ∈ C˜∗ , we have (f, `) ∈ E where ` is the center of L̃ that is the closest to f . Denote NΓ (`)
the neighbors of the point corresponding to center ` in Γ.
For each edge (f, `) ∈ E, for any client c ∈ NC ∗ (f ) − NL (`), we define Reassignc as the cost of
reassigning client c to `. We derive the following lemma.
Lemma 5.4. For any client c, Reassignc ≤ lc + 2gc .
Proof. By definition we have Reassignc = dist(c, `). By the triangle inequality dist(c, `) ≤ dist(c, f )+
dist(f, `). Since f serves c in C ∗ we have dist(c, f ) = gc , hence dist(c, `) ≤ gc + dist(f, `). We now
bound dist(f, `). Consider the center `0 that serves c in solution L. By the triangle inequality we
have dist(f, `0 ) ≤ dist(f, c) + dist(c, `0 ) = gc + lc . Finally, since ` is the closest center of f in L, we
have dist(f, `) ≤ dist(f, `0 ) ≤ gc + lc and the lemma follows.
We partition the centers of L̃ as follows. Let L̃0 be the set of centers of L̃ that have degree 0
in Γ. Let L̃≤ε−1 be the set of centers of L̃ that have degree at least one and at most 1/ε in Γ. Let
L̃>ε−1 be the set of centers of L̃ that have degree greater than 1/ε in Γ.
We now partition the centers of L̃ and C˜∗ using the neighborhoods of the vertices of L̃ in Γ. We
start by iteratively constructing two set of pairs S≤ε−1 and S>ε−1 . For each center ` ∈ L̃≤ε−1 ∪L̃>ε−1 ,
we pick a set A` of |NΓ (`)| − 1 centers of L̃0 and define a pair ({`} ∪ A` , NΓ (`)). We then remove
A` from L̃0 and repeat. Let S≤ε−1 be the pairs that contain a center of L̃≤ε−1 and let S>ε−1 be the
remaining pairs.
The following lemma follows from the definition of the pairs.
˜∗

Lemma 5.5. Let (RL̃ , RC ) be a pair in S≤ε−1 ∪ S>ε−1 . If ` ∈ RL̃ , then for any f such that
˜∗
(f, `) ∈ E, f ∈ RC .
˜∗

Lemma 5.6. For any pair (RL̃ , RC ) ∈ S≤ε−1 we have that
X
c∈NC

˜∗
∗ (RC )

lc ≤

X
c∈NC

˜∗
∗ (RC )

11

gc + 2

X
NL (RL̃ )

gc .

˜∗

Proof. Consider the mixed solution M = L − RL̃ ∪ RC . For each point c, let mc denote the cost
of c in solution M . We have the following upper bounds

˜∗

if c ∈ NC ∗ (RC ).
gc
∗
mc ≤ Reassignc if c ∈ NL (RL̃ ) − NC ∗ (RC˜ ) and by Lemma 5.5.


lc
Otherwise.
Now, observe that the solution M differs from L by at most 2/ε centers. Thus, by 1/ε-local
optimality we have cost(L) ≤ cost(M ). Summing over all clients and simplifying, we obtain
X
X
X
X
lc +
lc ≤
gc +
Reassignc .
c∈NC ∗ (RC ∗ )
˜

c∈NL (RL̃ )−NC ∗ (RC ∗ )

c∈NC ∗ (RC ∗ )

˜

c∈NL (RL̃ )−NC ∗ (RC ∗ )

˜

˜

The lemma follows by combining with Lemma 5.4.
We now analyze the cost of the clients served by a center of L that has degree greater than ε−1
in Γ. The argument is very similar.
˜∗

Lemma 5.7. For any pair (RL̃ , RC ) ∈ S>ε−1 we have that
X
X
X
lc ≤
gc + 2(1 + ε)
gc .
c∈NC ∗ (RC ∗ )
˜

c∈NC ∗ (RC ∗ )
˜

NL (RL̃ )

ˆ For
Proof. Consider the center `ˆ ∈ RL̃ that has in-degree greater than ε−1 . Let L̂ = RL̃ − {`}.
∗
˜
each ` ∈ L̂, we associate a center f (`) in RC in such a way that each f (`) 6= f (`0 ), for ` 6= `0 . Note
˜∗
˜∗
that this is possible since |L̂| = |RC | − 1. Let f˜ be the center of RC that is not associated with
any center of L̂.
Now, for each center ` of L̂ we consider the mixed solution M ` = L − {`} ∪ {f (`)}. For each
client c, we bound its cost m`c in solution M ` . We have


if c ∈ NC ∗ (f (`)).
gc
`
mc = Reassignc if c ∈ NL (`) − NC ∗ (f (`)) and by Lemma 5.5.


lc
Otherwise.
Summing over all center ` ∈ L̂, we have by ε−1 -local optimality
X
X X
lc +
lc ≤
˜
c∈NC ∗ (RC ∗ )−NC ∗ (f˜)

`∈RL̃ c∈NL (`)

X

gc +

˜
c∈NC ∗ (RC ∗ )−NC ∗ (f˜)

X

X

Reassignc .

(2)

`∈RL̃ c∈NL (`)

We now complete the proof of the lemma by analyzing the cost of the clients in NC ∗ (f˜). We
∗
consider the
the center `∗
P center ` ∈ L̂ that minimizes the reassignment cost of its clients. (`Namely,
∗ ,f˜)
such that c∈NL (`∗ ) Reassignc is minimized. We then consider the solution M
= L−{`∗ }∪{f˜}.
(`∗ ,f˜)

∗

˜

For each client c, we bound its cost mc
in solution M (` ,f ) . We have


if c ∈ NC ∗ (f˜).
gc
(`∗ ,f˜)
mc
≤ Reassignc if c ∈ NL (`∗ ) − NC ∗ (f˜) and by Lemma 5.5.


lc
Otherwise.
12

Thus, summing over all clients c, we have by local optimality
X
X
X
lc ≤
gc +
lc +
c∈NC ∗ (f˜)

c∈NL (`∗ )−NC ∗ (f (`∗ ))

X

Reassignc .

(3)

c∈NL (`∗ )−NC ∗ (f (`∗ ))

c∈NC ∗ (f˜)

By Lemma 5.4, combining Equations 2 and 3 and averaging over all centers of L̂ we have
X
X
X
lc ≤
gc + 2(1 + ε)
gc .
c∈NC ∗ (RC ∗ )

c∈NC ∗ (RC ∗ )

˜

˜

NL (RL̃ )

We now turn to the proof of Theorem 5.3.
Proof of Theorem 5.3. Observe first that for any c ∈ NL (L̃) − NC ∗ (C˜∗ ), we have lc ≤ gc . This
follows from the fact that the center that serves c in C ∗ is in S and so in L and thus, we have
lc ≤ gc . Therefore
X
X
lc ≤
gc .
(4)
c∈NL (L̃)−NC ∗ (C˜∗ )

c∈NL (L̃)−NC ∗ (C˜∗ )

We now sum the equations of Lemmas 5.6 and 5.7 over all pairs and obtain


X
X
X 
X
X

lc ≤
gc + (2 + 2ε)
gc 

(RL̃ ,RC ∗ ) c∈NC ∗ (RC ∗ )∪NL (RL̃ )
˜

(RL̃ ,RC ∗ )

˜

X

˜

lc ≤

c∈NC ∗ (C˜∗ )∪NL (L̃)

c∈NC ∗ (RC ∗ )∪NL (RL̃ )
˜

X

gc + (2 + 2ε)

c∈NC ∗ (C˜∗ )∪NL (L̃)

NL (RL̃ )

X

gc .

c∈NL (L̃)

Therefore,
X
c∈NC ∗ (C˜∗ )−NL (L̃)

lc +

X
NL (L̃)

X

lc ≤

c∈NC ∗ (C˜∗ )−NL (L̃)

X

gc + (3 + 2ε)

gc .

c∈NL (L̃)

Additionally, we show that the analysis is tight (up to a (1 + ε) factor):
Proposition 5.8. For any ε > 0, there exists an infinite family of 3 − ε-perturbation-resilient
instances such that for any constant ε > 0, there exists a locally optimal solution that has cost at
least 3OPT.
Proof. Consider a tripartite graph with nodes O, C, and L, where O is the set of optimal centers, L
is the set of centers of a locally optimal solution, and C is the set of clients. We have |O| = |L| = k
and |C| = k 2 . We specify the distances as follows. First, assume some arbitrary but fixed ordering
on the elements of O, L, and C. Then dist(Oi )(Ci,j ) = 1 + ε/3 and dist(Li )(Cj,i ) = 3 for any
i, j ∈ [k]. All other distances are induced by the shortest path metric along the edges of the graph,
i.e. dist(Oi )(Cj,` ) = 7 + ε/3 and dist(Li )(Cj,` ) = 5 + 2ε/3 for j, ` 6= i. We first note that O is
indeed the optimal solution with a cost of k 2 · (1 + ε/3). Multiplying the distances dist(Oi )(Ci,j ) by
a factor of (3 − ε) for all i ∈ [k] and j mod k = i, still ensures that O is an optimal solution with a
13

cost of k 2 · (1 + ε/3) · (3 − ε) = k 2 · 3(1 − ε2 ), which shows that the instance is (3 − ε)-perturbation
resilient.
What remains to be shown is that L is locally optimal. Assume that we swap out s centers.
Due to symmetry, we can consider the solution {Oi |i ∈ [s]} ∪ {Li |i ∈ [k] − [s]}. Each of centers
{Oi |i ∈ [s]} serve k clients with a cost of k · s · (1 + ε/3). The remaining clients are served by
{Li |i ∈ [k] − [s]}, as 5 + 2ε/3 < 7 + ε/3. The cost amounts to s · (k − s) · 5 + 2ε/3 for the clients
that get reassigned and (k − s)2 · 3 for the remaining clients. Combining these three figures gives
2
us a cost of k 2 · 3 + ksε − s2 · (2 + 2ε/3) > k 2 · 3 + ksε + s2 · 3. For k > 3s
ε , this is greater than k 3,
the cost of L.

6

Spectral Separability

In this section we will study the spectral separability condition for the Euclidean k-means problem.
Definition 6.1 (Spectral Separation [74]7 ). Let (A, Rd , || · ||2 , k) be an input for k-means clustering in Euclidean space and let {C1∗ , . . . Ck∗ } denote an optimal clustering of A with centers
S = {c∗1 , . . . c∗k }. Denote by C an n × d matrix such that the row Ci = argmin ||Ai − c∗j ||2 . Denote
c∗j ∈S

by || · ||2 the spectral norm of a matrix. Then
(i, j) the following condition holds:


{C1∗ , . . . Ck∗ }

is γ-spectrally separated, if for any pair


1
1 
||c∗i − c∗j || ≥ γ ·  p ∗ + q
||A − C||2 .
|Ci |
|Cj∗ |
Nowadays, a standard preprocessing step in Euclidean k-means clustering is to project onto the
subspace spanned by the rank k-approximation. Indeed, this is the first step of the algorithm by
Kumar and Kannan [74] (see Algorithm 2).
Algorithm 2 k-means with spectral initialization [74]
Project points onto the best rank k subspace
Compute a clustering C with constant approximation factor on the projection
3: Initialize centroids of each cluster of C as centers in the original space
4: Run Lloyd’s k-means until convergence
1:
2:

In general, projecting onto the best rank k subspace and computing a constant approximation
on the projection results in a constant approximation in the original space. Kumar and Kannan [74]
and later Awasthi and Sheffet [15] gave tighter bounds if the spectral separation is large enough. Our
algorithm omits steps 3 and 4. Instead, we project onto slightly more dimensions and subsequently
use Local Search as the constant factor approximation in step 2. To utilize Local Search, we further
require a candidate set of solutions, which is described in Section B. For pseudocode, we refer to
Algorithm 3. Our main result is to show that, given spectral separability, this algorithm is PTAS
for k-means (Theorem 6.2).
Theorem 6.2. Let (A, Rd , || · ||2 , k) be an instance of Euclidean k-means √
clustering with optimal
clustering C = {C1∗ , . . . Ck∗ } and centers S = {c∗1 , . . . c∗k }. If C is more than 3 k-spectrally separated,
then Algorithm 3 is a polynomial time approximation scheme.
7

The proximity condition of Kumar and Kannan [74] implies the spectral separation condition.

14

Algorithm 3 SpectralLS
1:
2:
3:
4:
5:

Project points A onto the best rank k/ε subspace
Embed points into a random subspace of dimension O(ε−2 log n)
Compute candidate centers (Corollary B.3)
Local Search(Θ(ε−4 ))
Output clustering

We first recall the basic notions and definitions for Euclidean k-means. Let A ∈ Rn×d be a set of
points in d-dimensional Euclidean space, where the row Ai contains the coordinates of the ith point.
The singular value decomposition is defined as A = U ΣV T , where U ∈ Rn×d and V ∈ Rd×d are
orthogonal and Σ ∈ Rd×d is a diagonal matrix containing the singular values where per convention
the singular values are given in descending order, i.e. Σ1,1 = σ1 q
≥ Σ2,2 = σ2 ≥ . . . Σd,d = σd .
Pd
2
Denote the Euclidean norm of a d-dimensional vector x by ||x|| =
i=1 xi . The spectral norm
qP
d
2
and Frobenius norm are defined as ||A||2 = σ1 and ||A||F =
i=1 σi , respectively.
The best rank k approximation
min ||A − X||F is given via Ak = Uk ΣV T = U Σk V T =
rank(X)=k

U ΣVkT ,

VkT

where Uk , Σk and
consist of the first k columns of U , Σ and V T , respectively, and
are zero otherwise. The best rank k approximation also minimizes the spectral norm, that is
||A − Ak ||2 = σk+1 is minimal among all matrices of rank k. The following fact is well known
throughout k-means literature and will be used frequently throughout this section.
1 P
Fact 6.3. Let A be a set of points in Euclidean space and denote by c(A) = |A|
x∈A x the centroid
of A. Then the 1-means cost of any candidate center c can be decomposed via
X
X
||x − c||2 =
||x − c(A)||2 + |A| · ||c(A) − c||2
x∈A

x∈A

and
X
x∈A

||x − c(A)||2 =

1 XX
||x − y||2 .
2 · |A|
x∈A y∈A

Note that the centroid is the optimal 1-means
. . Ck } of
P center
P of A. For a clustering C = {C1 , . P
A with centers S = {c1 , . . . ck }, the cost is then ki=1 p∈Ci ||p − ci ||2 . Further, if ci = |C1i | p∈Ci p,
we can rewrite the objective function in matrix form by associating the
 ith point with the ith row of
 √ 1 ∗ if Ai ∈ C ∗
j
|Cj |
n×k
some matrix A and using the cluster matrix X ∈ R
with Xi,j =
to denote
0
else
membership. Note that X T X = I, i.e. X is an orthogonal projection and that ||A−XX T A||2F is the
cost of the optimal k-means clustering. k-means is therefore a constrained rank k-approximation
problem.
We first restate the separation condition.
Definition 6.4 (Spectral Separation). Let A be a set of points and let {C1 , . . . Ck } be a clustering
of A with centers {c1 , . . . ck }. Denote by C an n×d matrix such that Ci = argmin ||Ai −cj ||2 . Then
j∈{1,...,k}

{C1 , . . . Ck } is γ spectrally separated, if for any pair of centers ci and cj the following condition
holds:
!
1
1
+p
||A − C||2 .
||ci − cj || ≥ γ · p
|Cj |
|Ci |
15

The following crucial lemma relates spectral separation and distribution stability.
Lemma 6.5. For a point set A, let C = {C1 , . . . , Ck } be an√optimal clustering with centers S =
{c1 , . . . , ck } associated clustering matrix X that is at least γ · k spectrally separated, where γ > 3.
For ε > 0, let Am be the best rank m = k/ε approximation of A. Then there exists a clustering
K = {C10 , . . . C20 } and a set of centers Sk , such that
1. the cost of clustering Am with centers Sk via the assignment of K is less than ||Am −
XX T Am ||2F and
2. (K, Sk ) is Ω((γ − 3)2 · ε)-distribution stable.
We note that this lemma would also allow us to use the PTAS of Awasthi et al. [12]. Before
giving the proof, we outline how Lemma 6.5 helps us prove Theorem 6.2. We first notice that if the
rank of A is of order k, then elementary bounds on matrix norm show that spectral separability
implies distribution stability. We aim to combine this observation with the following theorem due
to Cohen et al. [36]. Informally, it states that for every rank k approximation, (an in particular for
every constrained rank k approximation such as k-means clustering), projecting to the best rank
k/ε subspace is cost-preserving.
Theorem 6.6 (Theorem 7 of [36]). For any A ∈ Rn×d , let A0 be the rank dk/εe-approximation of
A. Then there exists some positive number c such that for any rank k orthogonal projection P ,
||A − P A||2F ≤ ||A0 − P A0 ||2F + c ≤ (1 + ε)||A − P A||2F .
The combination of the low rank case and this theorem is not trivial as points may be closer to a
wrong center after projecting, see also Figure 2. Lemma 6.5 determines the existence of a clustering
whose cost for the projected points Am is at most the cost of C ∗ . Moreover, this clustering has
constant distribution stability as well which, combined with the results from Section B, allows us to
use Local Search. Given that we can find a clustering with cost at most (1 + ε) · ||Am − XX T Am ||2F ,
Theorem 6.6 implies that we will have a (1 + ε)2 -approximation overall.
To prove the lemma, we will require the following steps:
• A lower bound on the distance of the projected centers ||ci Vm VmT − cj Vm VmT || ≈ ||ci − cj ||.
∗ = {c V V T , . . . , c∗ V V T } of A
• Find a clustering K with centers Sm
m with cost less than
1 m m
k m m
T
2
||Am − XX Am ||F .

• Show that in a well-defined sense, K and C ∗ agree on a large fraction of points.
• For any point x ∈ Ki , show that the distance of x to any center not associated with Ki is
large.
We first require a technical statement.
Lemma 6.7. For a point set A, let C = {C1 , . . . Ck } be a clustering with associated clustering
matrix X and let A0 and A00 be optimal low rank approximations where without loss of generality
k ≤ rank(A0 ) < rank(A00 ). Then for each cluster Ci
s


1 X
A00j − A0j
|Ci |
j∈Ci

≤
2

16

k
· ||A − XX T A||2 .
|Ci |

p

cm
i

ci
∆i
γ(∆i + ∆j )

cj
∆j
cm
j

Figure 2: Despite the centroids of each cluster being close after computing the best rank m approximation, the projection of a point p to the line connecting the centroid of cluster Ci and Cj
can change after computing the best
q rank m approximation. In this case ||p − cj || < ||p − ci || and
m
||p − cm
i || < ||p − cj ||. (Here ∆i =

k
|Ci | ||A

− XX T A||2 .)

17

P
Proof. By Fact 6.3 |Ci | · || |C1i | j∈Ci (A00i − A0i )||22 is, for a set of point indexes Ci , the cost of moving
the centroid of the cluster computed on A00 to the centroid of the cluster computed on A0 . For a
clustering matrix X, ||XX T A0 − XX T A0 ||2F is the sum of squared distances of moving the centroids
computed on the point set A00 to the centroids computed on A0 . We then have
2

1 X 00
(Aj − A0j )
|Ci |·
|Ci |
j∈Ci

2
≤ ||XX T A00 −XX T A0 ||2F ≤ ||X||2F ·||A00 −A0 ||22 ≤ k·σk+1
≤ k·||A−XX T A||22 .
2

Proof of Lemma 6.5. For any point p associated with some row of A, let pm = pVm VmT be the
corresponding row in Am . Similarly, for some cluster Ci , denote the center in A by ci and the
k
k
center in Am by cm
i . Extend these notion analogously for projections p and ci to the span of the
best rank k approximation Ak .
We have for any m ≥ k i 6= j
m
m
m
||cm
i − cj || ≥ ||ci − cj || − ||ci − ci || − ||cj − cj ||
!
√
1
1
≥ γ· √ +p
k||A − XX T A||2
Ci
|Cj |
1 √
1 √
k||A − XX T A||2 − p
k||A − XX T A||2
−√
Ci
|Cj |
!
√
1
1
= (γ − 1) · √ + p
k||A − XX T A||2 ,
Ci
|Cj |

(5)

where the second inequality follows
√ from Lemma 6.7.
k
√
||A − XX T A||2 . We will now construct our target clustering
In the following, let ∆i =
|Ci |

K. Note that we require this clustering (and its properties) only for the analysis. We distinguish
between the following three cases.
m
Case 1: p ∈ Ci and cm
i = argmin||p − cj ||:
j∈{1,...,k}

m
These points remain assigned to cm
i . The distance between pm and a different center cj is at
γ−1
1
m
m
least 2 ||ci − cj || ≥ 2 ε(∆i + ∆j ) due to Equation 5.
m
k
k
k
Case 2: p ∈ Ci , cm
i 6= argmin||p − cj ||, and ci 6= argmin||p − cj ||:
j∈{1,...,k}

j∈{1,...,k}

These points will get reassigned to their closest center.
1
m
m
The distance between pm and a different center cm
j is at least 2 ||ci − cj || ≥
due to Equation 5.

γ−1
2 ε(∆i

+ ∆j )

m
m
k
k
k
Case 3: p ∈ Ci , cm
i 6= argmin||p − cj ||, and ci = argmin||p − cj ||:
j∈{1,...,k}

pm

j∈{1,...,k}

cm
i

We assign
to
at the cost of a slightly weaker movement bound on the distance between
pm and cm
.
Due
to
orthogonality
of V , we have for m > k, (Vm − Vk )T Vk = VkT (Vm − Vk ) = 0.
j
Hence Vm VmT Vk = Vm VkT Vk + Vm (Vm − Vk )T Vk = Vk VkT Vk + (Vm − Vk )VkT Vk = Vk VkT Vk = Vk .
Then pk = pVk VkT = pVm VmT Vk VkT = pm Vk VkT .
18

Further, ||pk − ckj || ≥ 12 ||ckj − cki || ≥
between pm and a different center cm
j

γ−1
2 (∆i

+ ∆j ) due to Equation 5. Then the distance

m
k
m
k
||pm − cm
j || ≥ ||p − cj || − ||cj − cj || =

≥ ||pk − ckj || − ∆j ≥

q
k
||pm − pk ||2 + ||pk − ckj ||2 − ||cm
j − cj ||

γ−3
(∆i + ∆j ),
2

where the equality follows from orthogonality and the second to last inequality follows from
Lemma 6.7.
m
Now, given the centers {cm
1 , . . . ck }, we obtain a center matrix MK where the ith row of MK
is the center according to the assignment of above. Since both clusterings use the same centers
but K improves locally on the assignments, we have ||Am − MK ||2F ≤ ||Am − XX T Am ||2F , which
proves the first statement of the lemma. Additionally, due to the fact that Am − XX T Am has rank
m = k/ε, we have

||Am − MK ||2F ≤ ||Am − XX T Am ||2F ≤ m · ||Am − XX T Am ||22 ≤ k/ε · ||A − XX T A||2F

(6)

To ensure stability, we will show that for each element of K there exists an element of C, such
that both clusters agree on a large fraction of points. This can be proven by using techniques from
Awasthi and Sheffet [15] (Theorem 3.1) and Kumar and Kannan [74] (Theorem 5.4), which we
repeat for completeness.
Lemma 6.8. Let K = {C10 , . . . Ck0 } and C = {C1 , . . . Ck } be defined as above. Then there exists a
bijection b : C → K such that for any i ∈ {i, . . . , k}




32
32
1−
|Ci | ≤ b(|Ci |) ≤ 1 +
|Ci |.
(γ − 1)2
(γ − 1)2
Proof. Denote by Ti→j the set of points from Ci such that ||cki − pk || > ||ckj − pk ||. We first
2
note that ||Ak − XX T A||2F ≤ 2k · ||Ak − XX T A||22 ≤ 2k · ||A − Ak ||2 + ||A − XX T A||2 ≤ 8k ·
||A − XX T A||22 ≤ 8 · |Ci | · ∆2i for any i ∈ {1, . . . , k}. The distance ||pk − cki || ≥ 12 ||cki − ckj || ≥


√
γ−1
1
1
√
√
+
k||A − XX T A||22 . Assigning these points to cki , we can bound the total
·
2
C
i

|Cj |

number of points added to and subtracted from cluster Cj by observing

γ−1 2
|Ti→j | ≤
|Ti→j | ·
· (∆i + ∆j )2 ≤ ||Ak − XX T A||2F ≤ 8 · |Cj | · ∆2j
2
i6=j
i6=j


X
X
γ−1 2
2
∆j
|Tj→i | ≤
|Tj→i | ·
· (∆i + ∆j )2 ≤ ||Ak − XX T A||2F ≤ 8 · |Cj | · ∆2j .
2
∆2j

X

X

i6=j

j6=i




Therefore, the cluster sizes are up to some multiplicative factor of 1 ±

19

32
(γ−1)2



identical.

We now have for each point pm ∈ Ci0 a minimum cost of
!2
!
√
1
1
γ−3
T
m
m 2
· k · ||A − XX A||2
||p − cj || ≥
+p
· p
2
|Cj |
|Ci |


2

v
v
u
u
1
1
 γ − 3 u 
 √

u
T


≥ 
+ t
· t
 · k · ||A − XX A||2 
32
32
2
1 + (γ−1)2 · |Ci0 |
1 + (γ−1)2 · |Cj0 |
≥

4 · (γ − 3)2 ||Am − MK ||2F
·ε
81
|Cj0 |

where the first inequality holds due to Case 3, the second inequality holds due to Lemma 6.8 and
the last inequality follows from γ > 3 and Equation 6. This ensures that the distribution stability
condition is satisfied.
Proof of Theorem 6.2. Given the optimal clustering C ∗ of A with clustering matrix X, Lemma 6.5
guarantees the existence of a clustering K with center matrix MK such that ||Am − MK ||2F ≤
||Am − XX T Am || and that C has constant distribution stability. If ||Am − MK ||2F is not a constant
factor approximation, we are already done, as Local Search is guaranteed to find a constant factor
approximation. Otherwise due to Corollary B.3 (Section B in the appendix), there exists a discretization (Am , F, ||·||2 , k) of (Am , Rd , ||·||2 , k) such that the clustering C of the first instance has at
most (1+ε) times the cost of C in the second instance and such that C has constant distribution stability. By Theorem 4.2, Local Search with appropriate (but constant) neighborhood size will find a
clustering C 0 with cost at most (1+ε) times the cost of K in (Am , F, ||·||2 , k). Let Y be the clustering
matrix of C 0 . We then have ||Am −Y Y T Am ||2F +||A−Am ||2F ≤ (1+ε)2 ||Am −MK ||2F +||A−Am ||2F ≤
(1 + ε)2 ||Am − XX T Am ||2F + ||A − Am ||2F ≤ (1 + ε)3 ||A − XX T A||2F due to Theorem 6.6. Rescaling
ε completes the proof.
Remark. Any (1 + ε)-approximation will not in general agree with a target clustering. To see this
consider two clusters: (1) with mean on the origin and (2) with mean δ on the the first axis and
0 on all other coordinates. We generate points via a multivariate Gaussian distribution with an
identity covariance matrix centered on the mean of each cluster. If we generate enough points, the
instance will have constant spectral separability. However, if δ is small and the dimension large
enough, an optimal 1-clustering will approximate the k-means objective.

7

A Brief Survey on Stability Conditions

There are two general aims that shape the definitions of stability conditions. First, we want the
objective function to be appropriate. For instance, if the data is generated by mixture of Gaussians,
the k-means objective will be more appropriate than the k-median objective. Secondly, we assume
that there exists some ground truth, i.e. a correct assignment of points into clusters. Our objective
is to recover this ground truth as well as possible. These aims are not mutually exclusive. For
instance, an ideal objective function will allow us to recover the ground truth. We refer to Figure 3
for a visual overview of stability conditions and their relationships.

7.1

Cost-Based Separation

Given that an algorithm optimized with respect to some objective function, it is natural to define
a stability condition as a property the optimum clustering is required to have.
20

ORSS-Stability [85] Assume that we want to cluster a data set with respect to the k-means
objective, but have not decided on the number of clusters. A simple way of determining the
”correct” value of k is to run a k-means algorithm for k{1, 2, . . . m} until the objective value
decreases only marginally (using m centers). At this point, we set k = m − 1. The reasoning behind
this method, commonly known as the elbow-method is that we do not gain much information by
using m instead of m − 1 clusters, so we should favor the simpler model. Contrariwise, this implies
that we did gain information going from m − 2 to m − 1 and, in particular, that the m − 2-means
cost was considerably larger than the m − 1-means cost.
Ostrovsky et al. [85] considered whether such discrepancies in the cost also allow us to solve the
k-means problem more efficiently, see also Schulman [88] for an earlier condition for two clusters
and the irreducibility condition by Kumar et al. [75]. Specifically, they assumed that the optimal
k-means clustering has only an ε2 -fraction of the cost of the optimal (k − 1)-means clustering. For
such cost separated instances, the popular D2 -sampling technique has an improved performance
compared to the worst-case O(log k)-approximation ratio [9, 31, 66, 85]. Awasthi et al. [12] showed
that if an instance is cost-stable, it also admits a PTAS. In fact, they also showed that the weaker
condition β-stability is sufficient. β-stability states that the cost of assigning a point of cluster Ci to
another cluster Cj costs at least β times the total cost divided by the size of cluster Ci . Despite its
focus on the properties of the optimum, β-stability has many connections to target-clustering (see
below). Nowadays, the cost-stable property is one of the strongest stability conditions, implying
both distribution stability and spectral separability (see below). It is nevertheless the arguably
most intuitive stability condition.
Perturbation Resilience The other main optimum-based stability condition is perturbation
resilience. It was originally considered for the weighted max-cut problem by Bilu et al. [29, 28].
There, the optimum max cut is said to be α-perturbation resilient, if it remains the optimum even
if we multiply any edge weight up to a factor of α > 1. This notion naturally extends to metric
clustering problems, where, given a n×n distance matrix, the optimum clustering is α-perturbation
resilient if it remains optimal if we multiply entries by a factor α. Perturbation resilience has some
similarity to smoothed analysis (see Arthur et al. [8, 10] for work on k-means). Both smoothed
analysis and perturbation stability aim to study a smaller, more interesting part of the instance
space as opposed to worst case analysis that covers the entire space. Perturbation resilience assumes
that the optimum clustering stands out among any alternative clustering and measures the degree
by which it stands out via α. Smooth analysis is motivated by considering a problem after applying
a random perturbation, which for example accounts for measurement errors.
Perturbation resilience is unique among the considered stability conditions in that we aim to
recover the optimum solution, as opposed to finding a good (1+ε) approximation. Awasthi et al. [13]
showed that 3-perturbation resilience is sufficient√to find the optimum k-median clustering, which
was further improved by Balcan and Liang to 1+ 2 [21] 8 and finally to 2 by Angelidakis et al. [80].
Ben-David and Reyzin [27] showed that recovering the optimal clustering is NP-hard if the instance
is less than 2-perturbation resilient. Balcan et al. [20] gave an algorithm that optimally solves
symmetric and asymmetric k-center on 2-perturbation resilient instances. Recently, Angelidakis
et al. gave an algorithm that determines the optimum cluster for almost all used center-based
clustering if the instance is 2-perturbation resilient [80].
8

These results also holds for a slightly more general condition called the center proximity condition.

21

Cost Separation
Ostrovsky, Rabani,
Schulman, Swamy [85]
Jaiswal, Garg [66]

Perturbation Resilience
Bilu, Daniely,
Linial, Saks [28, 29]
Awasthi, Blum, Sheffet [13]
Balcan, Liang [21]

Distribution Stability
Awasthi, Blum, Sheffet [12]

Center Proximity
Awasthi, Blum, Sheffet [13]
Balcan, Liang [21]

Approximation Stability
Balcan, Blum, Gupta [17, 18]

Spectral Separation
Kumar, Kannan [74]
Awasthi, Sheffet [15]

Figure 3: An overview over all definitions of well-clusterability. Arrows correspond to implication.
For example, if an instance is cost-separated then it is distribution-stable; therefore the algorithm
by Awasthi, Blum and Sheffet [12] also works for cost-separated instances. The three highlighted
stability definitions in the middle of the figure are considered in this paper.

7.2

Target-Based Stability

The notion of finding a target clustering is more prevalent in machine learning than minimizing an
objective function. Though optimizing an objective value plays an important part in this line of
research, our ultimate goal is to find a clustering C that is close to the target clustering C ∗ . The
distance between two clusterings is the fraction of points where C and C ∗ disagree when considering
an optimal matching of clusters in C to clusters in C ∗ .
When the points are generated from some (unknown) mixture model, we are also given an
implicit target clustering. As a result, much work has focused on finding such clusterings using
probabilistic assumptions, see, for instance, [2, 6, 25, 32, 40, 41, 42, 44, 68, 82, 94]. We would like
to highlight two conditions that make no probabilistic assumptions and have a particular emphasis
on the k-means and k-median objective functions.
Approximation Stability The first assumption is that finding the target clustering is related
to optimizing the k-means objective function. In the simplest case, the target clustering coincides
with the optimum k-means clustering, but this a strong assumption that Balcan et al. [17, 18] avoid.
Instead they consider instances where any clustering with cost within a factor c of the optimum has
a distance at most ε to the target clustering, a condition they call (c, ε)-approximation stability.
Balcan et al. [17, 18] then showed that this condition is sufficient to both bypass worst-case lower
bounds for the approximation factor, and to find a clustering with distance O(ε) from the target
clustering. The condition was extended to account for the presence of noisy data by Balcan et
al. [22]. This approach was improved for other min-sum clustering objectives such as correlation
clustering by Balcan and Braverman [19]. For constant c, (c, ε) approximation stability also implies
the β-stability condition of Awasthi et al. [12] with constant β, if the target clusters are greater

22

than εn.
Spectral Separability Another condition that relates target clustering recovery via the k-means
objective was introduced by Kumar and Kannan [74]. In order to give an intuitive explanation,
consider a mixture model consisting of k centers. If the mixture is in a low-dimensional space, and
assuming that we have, for instance, approximation stability with respect to the k-means objective,
we could simply use the algorithm by Balcan et al. [18]. If the mixture has many additional
dimensions, the previous conditions have scaling issues, as the k-means cost may increase with each
dimension, even if many of the additional dimensions mostly contain noise. The notion behind the
spectral separability condition is that if the means of the mixture are well-separated in the subspace
containing their centers, it should be possible to determine the mixture even with the added noise.
Slightly more formally, Kumar and Kannan state that a point satisfies a proximity condition
if the projection of a point onto the line connecting its cluster center to another cluster center
is Ω(k) standard deviations closer to its own center than to the other. The standard deviations
are scaled with respect to the spectral norm of the matrix in which the ith row is the difference
vector between the ith point and its cluster mean. Given that all but an ε-fraction of points satisfy
the proximity condition, Kumar and Kannan [74] gave an algorithm that computes a clustering
with distance O(ε) to the target. They also show that their condition is (much) weaker than the
cost-stability condition by Ostrovsky et al. [85] and discuss some implications
of cost-stability on
√
approximation factors. Awasthi and Sheffet [15] later showed that Ω( k) standard deviations are
sufficient to recover most of the results by Kumar and Kannan.

8

Acknowledgments

The authors thank their dedicated advisor for this project: Claire Mathieu. Without her, this
collaboration would not have been possible.
The second author acknowledges the support by Deutsche Forschungsgemeinschaft within the
Collaborative Research Center SFB 876, project A2, and the Google Focused Award on Web
Algorithmics for Large-scale Data Analysis.

Appendix
A

(β, δ)-Stability

Lemma 4.5. Let Ci∗ be a cheap cluster. For any ε0 , we have |IRεi 0 ∩ Ci∗ | > (1 − ε3 /ε0 )|Ci∗ |.
Proof. Observe that each client that is not in IRiε0 is at a distance larger than ε0 βcost(C ∗ )/|Ci∗ |
from c∗i . Since Ci∗ is cheap, the total cost of the clients in Ci∗ = (IRεi 0 ∩ Ci∗ ) ∪ (Ci∗ − IRεi 0 ) is at most
ε3 βcost(C ∗ ) and in particular, the total cost of the clients in Ci∗ −IRεi 0 does not exceed ε3 βcost(C ∗ ).
Therefore, the total number of such clients is at most ε3 βcost(C ∗ )/(ε0 βcost(C ∗ )/|Ci∗ |) = ε3 |Ci∗ |/ε0 .
Lemma 4.6. Let δ +

ε3
ε0

< 1. If Ci∗ 6= Cj∗ are cheap clusters, then IRεi 0 ∩ IRεj 0 = ∅.

Proof. Assume that the claim is not true and consider a client x ∈ IRεi 0 ∩ IRεj 0 . Without loss of
generality assume |Ci∗ | ≥ |Cj∗ |. By the triangular inequality, we have cost(c∗j , c∗i ) ≤ cost(c∗j , x) +
cost(x, c∗i ) ≤ ε0 βcost(C ∗ )/|Cj∗ | + ε0 βcost(C ∗ )/|Ci∗ | ≤ 2ε0 βcost(C ∗ )/|Cj∗ |. Since the instance is
(β, δ)-distribution stable with respect to (C ∗ , S ∗ ) and due to Lemma 4.5, we have |∆i |+|IRεi 0 ∩Ci∗ | >
23

(1−δ)|Ci∗ |+(1−ε3 /ε0 )|Ci∗ | = (2−δ−ε3 /ε0 )|Ci∗ |. For δ+ε3 /ε0 < 1, there exists a client x0 ∈ IRεi 0 ∩∆i .
Thus, we have cost(x0 , c∗j ) ≤ cost(x0 , c∗i )+cost(c∗j , c∗i ) ≤ 3ε0 βcost(C ∗ )/|Cj∗ | < βcost(C ∗ )/|Cj∗ |. Since
x0 is in ∆i , we have cost(x0 , c∗j ) ≥ βcost(C ∗ )/|Cj∗ | resulting in a contradiction.
Lemma 4.8. There exists a set Z2 ⊆ C ∗ − Z1 of S
size at most 11.25ε−1 β −1 such that for any cluster
∗
∗
Cj ∈ C − Z2 , the total number of clients x ∈ i6=j ∆i , that are served by L(j) in L, is at most
2

ε|IRεi ∩ Ci∗ |.
Proof. Consider a cheap cluster Cj∗ ∈ C ∗ − Z1 such that the total number of clients x ∈ ∆i for
2

i 6= j, that are served by L(j) in L, is greater than ε|IRεj ∩ Cj∗ |. By the triangular inequality and
the definition of (β, δ)-stability, the total cost for each x ∈ ∆i with i 6= j served by L(j) is at least
2
(1 − ε)βcost(C ∗ )/|Cj∗ |. Since there are at least ε|IRεj ∩ Cj∗ | such clients, their total cost is at least
2

ε|IRεj ∩ Cj∗ |(1 − ε)βcost(C ∗ )/|Cj∗ |. By Lemma 4.5, this total cost is at least
2

ε|IRεj ∩ Cj∗ |(1 − ε)β

∗
cost(C ∗ )
2
∗ cost(C )
≥
ε(1
−
ε)
|C
|β
.
j
|Cj∗ |
|Cj∗ |

Recall that by [11], L is a 5-approximation and so there exist at most 11.25·ε−1 β −1 such clusters.
Lemma 4.10. Let Ci∗ be a cluster in C ∗ − Z ∗ . Define the solution Mi = L − {L(i)} ∪ {c∗i } and
denote by mix the cost of client x in solution Mi . Then
X
X
X
X
X
X
ε
Reassign(x) +
lx +
gx +
mix ≤
lx +
gx + lx ).
(
(1 − ε)
x∈A

x∈A−
(A(L(i))∪Ei )

x∈Di

x∈Ei

x∈Ei

x∈A(L(i))−
(Ei ∪Di )

Proof. Consider a client x ∈ Ci∗ − A(L(i)). By the triangular inequality, we have Reassign(x) =
cost(x, L(i)) ≤ cost(x, c∗i ) + cost(c∗i , L(i)) = gx + cost(c∗i , L(i)). Then,
X
X
Reassign(x) ≤
gx + |Ci∗ − A(L(i))| · cost(c∗i , L(i)).
x∈Ci∗ −A(L(i))

x∈Ci∗ −A(L(i))

Now consider the clients in Ci∗ ∩ A(L(i)). By the triangular inequality, we have cost(c∗i , L(i)) ≤
cost(c∗i , x0 ) + cost(x0 , L(i)) ≤ gx + lx . Therefore,
cost(c∗i , L(i)) ≤

|Ci∗

1
∩ A(L(i))|

X

(gx + lx ).

x∈Ci∗ ∩A(L(i))

|Ci∗ −A(L(i))|
ε2
∗
∗
|Ci∗ ∩A(L(i))| . Due to Lemma 4.5, we have |IRi ∩ Ci | ≥ (1 − ε)|Ci | and due to
2
2
Lemma 4.4, we have |IRεi ∩Ci∗ ∩A(L(i))| ≥ (1−ε)|IRεi ∩Ci∗ |. Therefore |Ci∗ ∩A(L(i))| ≥ (1−ε)2 |Ci∗ |
|C ∗ −A(L(i))|
2ε
and |Ci∗ − A(L(i))| ≤ (1 − (1 − ε)2 )|Ci∗ | ≤ 2ε|Ci∗ |, yielding |Ci∗ ∩A(L(i))| ≤ (1−ε)
2.
i

We now bound

Combining, we obtain
X

Reassign(x) ≤

x∈Ci∗ −A(L(i))

X

gx +

x∈Ci∗ −A(L(i))

≤

X

gx +

x∈Ci∗ −A(L(i))

24

|Ci∗ − A(L(i))|
|Ci∗ ∩ A(L(i))|
2ε
(1 − ε)2

X

(gx + lx )

x∈Ci∗ ∩A(L(i))

X
x∈Ci∗ ∩A(L(i))

(gx + lx ).

Lemma 4.11. Let Ci∗ be a cluster in C ∗ − Z ∗ . Define the solution Mi = L − {L(i)} ∪ {c∗i } and
denote by mic the cost of client c in solution Mi . Then
X
x∈A

mix ≤

X

lx +

x∈A−
(A(L(i))∪C̃i∗ )

X

X

gx +

x∈Di

x∈C̃i∗

X

Reassign(x) +

lx +

x∈A(L(i))−
(C̃i∗ ∪Di )

X
ε
(
gx + lx ).
(1 − ε)
∗
x∈C̃i

Proof. For any client x ∈ A − A(L(i)), the center that serves it in L belongs to Mi . Thus its cost
is at most lx . Moreover, observe that any client x ∈ Ei ⊆ Ci∗ can now be served by c∗i , and so its
cost is at most gx . For each client x ∈ Di , we bound its cost by Reassign(x) since all the centers of
L except for L(i) are in Mi and x ∈ Bj∗ ⊆ Cj∗ ∈ C ∗ − C(Z ∗ ).
Now, we bound the cost of a client x ∈ A(L(i)) − (Ei ∪ Di ) ⊆ A(L(i)). The closest center in
i
M for a client x0 ∈ A(L(i)) is not farther than c∗i . By the triangular inequality, the cost of such
client x0 is at most cost(x0 , c∗i ) ≤ cost(x0 , L(i)) + cost(L(i), c∗i ) = lx0 + cost(L(i), c∗i ), and so
X
X
mix ≤ |A(L(i)) − (Ei ∪ Di )| · cost(L(i), c∗i ) +
lx .
(7)
x∈A(L(i))−
(Ei ∪Di )

x∈A(L(i))−
(Ei ∪Di )

Now, observe that, for any client x ∈ |A(L(i)) ∩ Ei |, by the triangular inequality, we have
cost(L(i), c∗i ) ≤ cost(L(i), x) + cost(x, c∗i ) = lx + gx . Therefore,
1
|A(L(i)) ∩ Ei |

cost(L(i), c∗i ) ≤

X

(lx + gx ).

(8)

x∈A(L(i))∩Ei

Combining Equations 7 and 8, we have
X

X

mix ≤

x∈A(L(i))−
(Ei ∪Di )

lx +

x∈A(L(i))−
(Ei ∪Di )

X

≤

lx +

|A(L(i)) − (Ei ∪ Di )|
|A(L(i)) ∩ Ei |

X

(lx + gx )

x∈A(L(i))∩Ei

|A(L(i)) − Ei | X
(lx + gx ).
|A(L(i)) ∩ Ei |

(9)

x∈Ei

x∈A(L(i))
−(Ei ∪Di )

We now remark that since Ei is in C ∗ − Z ∗ , we have by Lemmas 4.7 and 4.8, |A(L(i)) − Ei | ≤
2
2
ε · |IRiε ∩ Ci∗ | and (1 − ε) · |IRiε ∩ Ci∗ | ≤ |A(L(i)) ∩ Ei |. Thus, combining with Equation 9 yields
the lemma.
Lemma 4.12. We have
X

−ε · cost(L) +

X

lx ≤

∗)
b
x∈A−C(Z

gx +

∗)
b
x∈A−C(Z

3ε
· (cost(L) + cost(C ∗ )).
(1 − ε)2

Proof. We consider a cluster Ci∗ in C ∗ − Z ∗ and the solution Mi = L − {L(i)} ∪ {c∗i }. Observe that
Mi and L only differ by L(i) and c∗i . Therefore, by local optimality we have (1 − nε ) · cost(Li ) ≤
cost(Mi ). Then Lemma 4.11 yields
ε
(1 − ) · cost(Li ) ≤
n

X
x∈A−
(A(L(i))∪Ei )

lx +

X
x∈Ei

gx +

X

Reassign(x) +

x∈Di

25

X
x∈A(L(i))−
(Ei ∪Di )

lx +

X
ε
·
(gx + lx )
(1 − ε)
x∈E

and so, simplifying
−

X
X
X
X
X
ε
ε
· cost(Li ) +
lx +
lx ≤
gx +
Reassign(x) +
·
(gx + lx )
n
(1 − ε)
x∈Ei

x∈Di

x∈Ei

x∈Di

x∈Ei

We now apply this analysis to each cluster Ci∗ ∈ C ∗ − Z ∗ . Summing over all clusters Ci∗ , we obtain,


|C ∗ −Z ∗ |
X
X
X
ε

− · cost(L)+
lx +
lx  ≤
n
i=1
x∈Ei
x∈Di


∗
∗
|C −Z |
X
X
X
ε

· (cost(L) + cost(C ∗ ))
gx +
Reassign(c) +
(1 − ε)
i=1

x∈Ei

x∈Di

By Lemma 4.10 and the definition of Ei ,
ε
− · cost(L) +
n
|C ∗ −Z ∗ |

≤

Therefore, −

B



X

X

i=1

b
x∈Ci∗ ∩A

ε
· cost(L) +
n

gx +

X
∗)
b
x∈A−C(Z

ε
2ε
+
1 − ε (1 − ε)2

lx ≤

X



gx +

∗)
b
x∈A−C(Z

|C ∗ −Z ∗ |

X

X

i=1

b
x∈Ci∗ ∩A

lx

· (cost(L) + cost(C ∗ )).

3ε
· (cost(L) + cost(C ∗ )).
(1 − ε)2

Euclidean Distribution Stability

In this section we show how to reduce the Euclidean problem to the discrete version. Our analysis
is focused on the k-means problem, however we note that the discretization works for all values of
cost = distp , where the dependency on p grows exponentially. For constant p, we obtain polynomial
sized candidate solution sets in polynomial time. For k-means itself, we could alternatively combine
Matousek’s approximate centroid set [81] with the Johnson Lindenstrauss lemma and avoid the
following construction; however this would only work for optimal distribution stable clusterings
and the proof Theorem 6.2 requires it to hold for non-optimal clusterings as well.
First, we describe a discretization procedure. It will be important to us that the candidate
solution preserves (1) the cost of any given set of centers and (2) distribution stability.
For a set of points P , a set of points Nε is an ε-net of P if for every point x ∈ P there exists some
point y ∈ Nε with ||x − y|| ≤ ε. It is well known that for unit Euclidean ball of dimension d, there
exists an ε-net of cardinality (1 + 2/ε)d , see for instance Pisier [87], though in this case the proof
is non-constructive. Constructive methods yield slightly worse, but asymptotically similar bounds
of the form ε−O(d) , see for instance Chazelle [35] for an extensive overview on how to construct
such nets. Note that having constructed an ε-net for the unit sphere, we also have an ε · r-net
for any sphere with radius r. The following lemma shows that a sufficiently small ε-net preserves
distribution stability. Again for ease of exposition, we only give the proof for p = 1, and assuming
we can construct an appropriate ε-net, but similar results also hold for (k, p) clustering as long as
p is constant.

26

Lemma B.1. Let A be a set of n points in d-dimensional Euclidean space and let β, ε > 0 with
min(β, ε) > 2η > 0 be constants. Suppose there exists a clustering C = {C1 , . . . , Ck } with centers
S = {c1 , . . . ck } such that
P P
1. cost(C, S) = ki=1 x∈Ci ||x − ci || is a constant approximation to the optimum clustering and
2. C is β-distribution stable.
Then there exists a discretization D of the solution space such that there exists a subset S 0 =
{c01 , . . . c0k } ⊂ D of size k with
Pk P
0
1.
i=1
x∈Ci ||x − ci || ≤ (1 + ε) · cost(C, S) and
2. C with centers S 0 is β/2-distribution stable.
The discretization consists of O(n · log n · η d+2 ) many points.
Proof. Let OPT being the cost of an optimal k-median clustering. Define an exponential sequence
to the base of (1 + η) starting at (η · OPT
n ) and ending at (n · OPT). The sequence contains
t = log1+η (n2 /η) ∈ O(η −1 log n) many elements for 1/η < n. For each point p ∈ A, define B(p, `i )
as the d-dimensional ball centered at p with radius (1+η)i ·η· OPT
n . We cover the ball B(p, `i ) with an
η/8 · `i net denoted by Nη/8 (p, `i ). As the set of candidate centers, we let D = ∪p∈A ∪ti=0 Nη/8 (p, `i ).
Clearly, |D| ∈ O(n · log n · (1 + 16/η)d+2 ).
Now for each ci ∈ S, set c0i = argmin ||q − ci ||. We will show that S 0 = {c01 , . . . c0k } satisfies the
q∈D

two conditions of the lemma.
0
For (1), we first consider the points p with ||p − ci || ≤ ε/8 · OPT
n . Then there exists a ci such
OPT
OPT
0
that ||p − ci || ≤ (η/8 + ε/8) n ≤ ε/4 n and summing up over all such points, we have a total
contribution to the objective value of at most ε/4 · OPT.
Now consider the remaining points. Since the cost(C, S) is a constant approximation, the center
i+1 ·η· OPT for some i ∈ {0, . . . t}. Then
ci of each point p satisfies (1+η)i ·η· OPT
n ≤ ||ci −p|| ≤ (1+η)
n
there exists some point q ∈ Nη/8 (p, `i+1 ) with ||q−ci || ≤ η/8·(1+η)i+1 ·η· OPT
n ≤ η/8·(1+η)||p−ci || ≤
η/4||p − ci ||. We then have ||p − c0i || ≤ (1 + η/4)||p − ci ||. Summing up over both cases, we have a
total cost of at most ε/4 · OPT + (1 + η/4) · cost(C, S 0 ) ≤ (1 + ε/2) · cost(C, S 0 ).
OPT
OPT
To show (2), let us consider some point p ∈
/ Cj with ||p−cj || > β · OPT
|Cj | . Since β · |Cj | ≥ 2η · n ,
i+1 · OPT .
there exists a point q and an i ∈ {0, . . . t} such that β/8·(1+η)i · OPT
n ≤ ||ci −q|| ≤ β/8·(1+η)
n
0 satisfies ||p − c0 || ≥
.
Similarly
to
above,
the
point
c
Then ||c0j − cj || ≤ β · (1 + η)i+1 · OPT
j
j
n
OPT
OPT
OPT
||p − cj || − ||cj − c0j || ≥ β · OPT
|Cj | − β/8(1 + η) · n ≥ (1 − 1/4)β · |Cj | > β/2 · |Cj | .

To reduce the dependency on the dimension, we combine this statement with the seminal
theorem originally due to Johnson and Lindenstrauss [67].
Lemma B.2 (Johnson-Lindenstrauss lemma). For any set of n points N in d-dimensional Euclidean space and any 0 < ε < 1/2, there exists a distribution F over linear maps f : `d2 → `m
2 with
m ∈ O(ε−2 log n) such that
2
Pf ∼F [∀x, y ∈ N, (1 − ε)||x − y|| ≤ ||f (x) − f (y)|| ≤ (1 + ε)||x − y||] ≥ .
3
It is easy to see that Johnson-Lindenstrauss type embeddings preserve the Euclidean k-means
cost of any clustering, as the cost of any clustering can be written in terms of pairwise distances
(see also Fact 6.3 in Section 6). Since the distribution over linear maps F can be chosen obliviously
27

with respect to the points, this extends to distribution stability of a set of k candidate centers as
well.
Combining Lemmas B.2 and B.1 gives us the following corollary.
Corollary B.3. Let A be a set of points in d-dimensional Euclidean space with a clustering C =
{C1 , . . . Ck } and centers S = {c1 , . . . ck } such that C is β-perturbation stable. Then there exists a
−1
(A, F, || · ||2 , k)-clustering instance with clients A, npoly(ε ) centers F and a subset S 0 ⊂ F ∪ A of
k centers such that C and S 0 is O(β) stable and the cost of clustering A with S 0 is at most (1 + ε)
times the cost of clustering A with S.
Remark. This procedure can be adapted to work for general powers of cost functions. For Lemma B.1,
we simply rescale η. The Johnson-Lindenstrauss lemma can also be applied in these settings, at
a slightly worse target dimension of O((p + 1)2 log((p + 1)/ε)ε−3 log n), see Kerber and Raghvendra [71].

C

Experimental Results

In this section, we discuss the empirical applicability of stability as a model to capture real-world
−3 −1
data. Theorem 4.2 states that local search with neighborhood of size nΩ(ε β ) returns a solution
of cost at most (1 + ε)OPT. Thus, we ask the following question.
For which values of β are the random and real instances β-distribution-stable?
We focus on the k-means objective and we consider real-world and random instances with
ground truth clustering and study under which conditions the value of the solution induced by the
ground truth clustering is close to the value of the optimal clustering with respect to the k-means
objective. Our aim is to determine (a range of) values of β for which various data sets satisfy
distribution stability.

Setup
The machines used for the experiments have a processor Intel(R) Core(TM) i73770 CPU, 3.40GHz
with four cores and a total virtual memory of 8GB running on an Ubuntu 12.04.5 LTS operating
system. We implemented the Algorithms in C++ and Python. The C++ compiler is g++ 4.6.3.
Our experiments always used Local Search with a neighborhood of size 1. At each step, the neighborhood of the current solution was explored in parallel: 8 threads were created by a Python script
and each of them correspond to a C++ subprocess that explores a 1/8 fraction of the space of
the neighboring solutions. The best neighboring solution found by the 8 threads was taken for the
next step. For Lloyd’s algorithm we use the C++ implementation by Kanungo et al. [70] available
online.
To determine the stability parameter β, we also required a lower bound on the cost. This was
done via a linear relaxation describe in Algorithm 4. The LP for the linear program was generated
via a Python script and solved using the solver CPLEX. The average ratio between our upper bound
given via Local Search and lower bounds given via Algorithm 4 is 1.15 and the variance for the
value of the optimal fractional solution that is less than 0.5% of the value of the optimal solution.
Therefore, our estimate of β is quite accurate.

28

Algorithm 4 Linear relaxation for the k-means problem.
Input: A set of clients A, a set of candidates centers F , a number of centers k, a distance function
dist.
XX
min
xa,b · dist(a, b)2
a∈A b∈F

subject to,
P
P b∈F yb
∀a ∈ A,
b∈F xa,b
∀a ∈ A, ∀b ∈ F,
yb
∀a ∈ A, ∀b ∈ F,
xa,b

C.1

≤ k
= 1
≥ xa,b
≥ 0

Real Data

In this section, we focus on four classic real-world datasets with ground truth clustering: abalone,
digits, iris, and movement libras. abalone, iris, and movement libras have been used in
various works (see [46, 47, 49, 50, 89] for example) and are available online at the UCI Machine
learning repository [78].
The abalone dataset consists of 8 physical characteristics of all the individuals of a population
of abalones. Each abalone corresponds to a point in a 8-dimensional Euclidean space. The ground
truth clustering consists in partitioning the points according to the age of the abalones.
The digits dataset consists of 8px-by-8px images of handwritten digits from the standard
machine learning library scikit-learn [86]. Each image is associated to a point in a 64-dimensional
Euclidean space where each pixel corresponds to a coordinate. The ground truth clustering consists
in partitioning the points according to the number depicted in their corresponding images.
The iris dataset consists of the sepal and petal lengths and widths of all the individuals of a
population of iris plant containing 3 different types of iris plant. Each plant is associated to a point
in 4-dimensional Euclidean space. The ground truth clustering consists in partitioning the points
according to the type of iris plant of the corresponding individual.
The Movement libras dataset consists of a set of instances of 15 hand movements in LIBRAS9 .
Each instance is a curve that is mapped in a representation with 90 numeric values representing
the coordinates of the movements. The ground truth clustering consists in partitioning the points
according to the type of the movement they correspond to.
Properties
Number of points
Number of clusters
Value of ground truth clustering
Value of fractional relaxation
Value of Algorithm 1
% of pts correct. class. by Alg. 1
β-stability

Abalone
636
28
169.19
4.47
4.53
17
1.27e-06

Digits
1000
10
938817.0
855567.0
855567.0
76.2
0.0676

Iris
150
3
96.1
83.96
83.96
90
0.2185

Movement libras
360
15
780.96
366.34
369.65
39
0.0065

Table 1: Properties of the real-world instances with ground truth clustering. The neighborhood
size for Algorithm 1 is 1.
Table 1 shows the properties of the four instances.
For the Abalone and Movement libras instances, the values of an optimal solution is much
9

LIBRAS is the official Brazilian sign language

29

smaller than the value of the ground truth clustering. Therefore the k-means objective function
might not be ideal as a recovery mechanism. Since Local Search optimizes with respect to the
k-means objective, the clustering output by Local Search is far from the ground truth clustering
for those instances: the percentage of points correctly classified by Algorithm 1 is at most 17% for
the Abalone instance and at most 39% for the Movement libras instance. For the Digits and
Iris instances the value of the ground truth clustering is at most 1.15 times the optimal value.
In those cases, the number of points correctly classified is much higher: 90% for the Iris instance
and 76.2% for the Digits instance.
The experiments also show that the β-distribution-stability condition is satisfied for β > 0.06
for the Digits, Iris and Movement libras instances. This shows that the β-distribution-stability
condition captures the structure of some famous real-world instances for which the k-means objective is meaningful for finding the optimal clusters. We thus make the following observations.
Observation C.1. If the value of the ground truth clustering is close to the value of the optimal
solution, then one can expect the instance satisfy the β-distribution stability property for some
constant β.
The experiments show that Algorithm 1 with neighborhood size 1 (s = 1) is very efficient for
all those instances since it returns a solution whose value is within 2% of the optimal solution for
the Abalone instance and a within 0.002% for the other instances. Note that the running time
of Algorithm 1 with s = 1 is Õ(k · n/ε) (using a set of O(n) candidate centers) and less than 15
minutes for all the instances. We make the following observation.
Observation C.2. If the value of the ground truth clustering is close to the value of the optimal
solution, then one can expect both clusterings to agree on a large fraction of points.
Finally, observe that for those instances the value of an optimal solution to the fractional
relaxation of the linear program is very close to the optimal value of an optimal integral solution
(since the cost of the integral solution is smaller than the cost returned by Algorithm 1). This
suggests that the fractional relaxation (Algorithm 4) might have a small integrality gap for realworld instances.
Open Problem: We believe that it would be interesting to study the integrality gap of the classic
LP relaxation for the k-median and k-means problems under the stability assumption (for example
β-distribution stability).

C.2

Data generated from a mixture of k Gaussians

The synthetic data was generated via a Python script using numpy. The instances consist of 1000
points generated from a mixture of k Gaussians with the same variance σ lying in d-dimensional
space, where d ∈ {5, 10, 50} and k ∈ {5, 50, 100}. We generate 100 instances for all possible combinations of the parameters. The means of the k Gaussians are chosen uniformly and independently
at random in Qd ∩ (0, 1)d . The ground truth clustering is the family of sets of points generated
by the same Gaussian. We compare the value of the ground truth clustering to the optimal value
clustering.
The results are presented in Figures 4 and 5. We observe that when the variance σ is large, the
ratio between the average value of the ground truth clustering and the average value of the optimal
clustering becomes more important. Indeed, the ground truth clusters start to overlap, allowing to
improve the objective value by defining slightly different clusters. Therefore, the use of the k-means

30

or k-median objectives for modeling the recovery problem is not suitable anymore. In these cases,
since Local Search optimizes the solution with respect to the current cost, the clustering output by
local search is very different from the ground truth clustering. We thus identify instances for which
the k-means objective is meaningful and so, Local Search is a relevant heuristic. This motivates
the following defintion.
Definition C.3. We say that a variance σ̂ is relevant if, for the k-means instances generated
with variance σ̂ the ratio between the average value of the ground truth clustering and the optimal
clustering is less than 1.05.
We summarize in Table 2 the relevant variances observed.
hhh
hh

hhh
hh

Values of k
hhhh

Number of dimensions
2
10
50

5

50

100

hhh
h
h

< 0.05
< 15
< 1000000.0

< 0.002
<1
< 100

< 0.0005
< 0.5
<7

Table 2: Relevant variances for k ∈ {5, 50, 100} and d ∈ {2, 10, 50}.
We consider the β-distribution-stability condition and ask whether the instances generated
from a relevant variance satisfy this condition for constant values of β. We remark that β can take
arbitrarily small values.
We thus identify relevant variances (see Table 2) for each pair k, d, such that optimizing the
k-means objective in a d-dimensional instances generated from a relevant variance corresponds to
finding the underlying clusters.
On stability conditions. We now study the β-distribution-stability condition for random instances generated from a mixture of k Gaussians. The results are depicted in Figures 7 and 6.
We observe that for random instances that are not generated from a relevant variance, the
instances are β-distribution-stable for very small values of β (e.g., β < 1e − 07). We also make the
following observation.
Observation C.4. Instances generated using relevant variances satisfy the β-distribution-stability
condition for β > 0.001.
We remark that the number of dimensions is constant here and that having more dimensions
might incur slightly different values for β. It would be interesting to study this dependency in a
new study.

References
[1] Emile Aarts and Jan K. Lenstra, editors. Local Search in Combinatorial Optimization. John
Wiley & Sons, Inc., New York, NY, USA, 1st edition, 1997.
[2] Dimitris Achlioptas and Frank McSherry. On spectral learning of mixtures of distributions. In
Learning Theory, 18th Annual Conference on Learning Theory, COLT 2005, Bertinoro, Italy,
June 27-30, 2005, Proceedings, pages 458–469, 2005.

31

[3] Sara Ahmadian, Ashkan Norouzi-Fard, Ola Svensson, and Justin Ward. Better guarantees for
k-means and euclidean k-median by primal-dual algorithms. CoRR, abs/1612.07925, 2016.
[4] Ehsan Ardjmand, Namkyu Park, Gary Weckman, and Mohammad Reza Amin-Naseri. The
discrete unconscious search and its application to uncapacitated facility location problem.
Computers & Industrial Engineering, 73:32 – 40, 2014.
[5] Ehsan Ardjmand, Namkyu Park, Gary R. Weckman, and Mohammad Reza Amin-Naseri.
The discrete unconscious search and its application to uncapacitated facility location problem.
Computers & Industrial Engineering, 73:32–40, 2014.
[6] Sanjeev Arora and Ravi Kannan. Learning mixtures of arbitrary gaussians. In Proceedings
on 33rd Annual ACM Symposium on Theory of Computing, July 6-8, 2001, Heraklion, Crete,
Greece, pages 247–257, 2001.
[7] Sanjeev Arora, Prabhakar Raghavan, and Satish Rao. Approximation schemes for Euclidean
k -medians and related problems. In Proceedings of the Thirtieth Annual ACM Symposium on
the Theory of Computing, Dallas, Texas, USA, May 23-26, 1998, pages 106–113, 1998.
[8] David Arthur, Bodo Manthey, and Heiko Röglin. Smoothed analysis of the k-means method.
J. ACM, 58(5):19, 2011.
[9] David Arthur and Sergei Vassilvitskii. k-means++: the advantages of careful seeding. In
Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA
2007, New Orleans, Louisiana, USA, January 7-9, 2007, pages 1027–1035, 2007.
[10] David Arthur and Sergei Vassilvitskii. Worst-case and smoothed analysis of the ICP algorithm,
with an application to the k-means method. SIAM J. Comput., 39(2):766–782, 2009.
[11] Vijay Arya, Naveen Garg, Rohit Khandekar, Adam Meyerson, Kamesh Munagala, and
Vinayaka Pandit. Local search heuristics for k-median and facility location problems. SIAM
J. Comput., 33(3):544–562, 2004.
[12] Pranjal Awasthi, Avrim Blum, and Or Sheffet. Stability yields a PTAS for k-median and
k-means clustering. In 51th Annual IEEE Symposium on Foundations of Computer Science,
FOCS 2010, October 23-26, 2010, Las Vegas, Nevada, USA, pages 309–318, 2010.
[13] Pranjal Awasthi, Avrim Blum, and Or Sheffet. Center-based clustering under perturbation
stability. Inf. Process. Lett., 112(1-2):49–54, 2012.
[14] Pranjal Awasthi, Moses Charikar, Ravishankar Krishnaswamy, and Ali Kemal Sinop. The
hardness of approximation of Euclidean k-means. In 31st International Symposium on Computational Geometry, SoCG 2015, June 22-25, 2015, Eindhoven, The Netherlands, pages 754–
767, 2015.
[15] Pranjal Awasthi and Or Sheffet. Improved spectral-norm bounds for clustering. In Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques - 15th
International Workshop, APPROX 2012, and 16th International Workshop, RANDOM 2012,
Cambridge, MA, USA, August 15-17, 2012. Proceedings, pages 37–49, 2012.
[16] Ainesh Bakshi and Nadiia Chepurko. Polynomial time algorithm for 2-stable clustering instances. CoRR, abs/1607.07431, 2016.
32

[17] Maria-Florina Balcan, Avrim Blum, and Anupam Gupta. Approximate clustering without the
approximation. In Proceedings of the Twentieth Annual ACM-SIAM Symposium on Discrete
Algorithms, SODA 2009, New York, NY, USA, January 4-6, 2009, pages 1068–1077, 2009.
[18] Maria-Florina Balcan, Avrim Blum, and Anupam Gupta. Clustering under approximation
stability. J. ACM, 60(2):8, 2013.
[19] Maria-Florina Balcan and Mark Braverman. Finding low error clusterings. In COLT 2009
- The 22nd Conference on Learning Theory, Montreal, Quebec, Canada, June 18-21, 2009,
2009.
[20] Maria-Florina Balcan, Nika Haghtalab, and Colin White. k-center clustering under perturbation resilience. In 43rd International Colloquium on Automata, Languages, and Programming,
ICALP 2016, July 11-15, 2016, Rome, Italy, pages 68:1–68:14, 2016.
[21] Maria-Florina Balcan and Yingyu Liang. Clustering under perturbation resilience. SIAM J.
Comput., 45(1):102–155, 2016.
[22] Maria-Florina Balcan, Heiko Röglin, and Shang-Hua Teng. Agnostic clustering. In Algorithmic
Learning Theory, 20th International Conference, ALT 2009, Porto, Portugal, October 3-5,
2009. Proceedings, pages 384–398, 2009.
[23] Sayan Bandyapadhyay and Kasturi R. Varadarajan. On variants of k-means clustering. CoRR,
abs/1512.02985, 2015.
[24] MohammadHossein Bateni, Aditya Bhaskara, Silvio Lattanzi, and Vahab S. Mirrokni. Distributed balanced clustering via mapping coresets. In Advances in Neural Information Processing Systems 27: Annual Conference on Neural Information Processing Systems 2014, December
8-13 2014, Montreal, Quebec, Canada, pages 2591–2599, 2014.
[25] Mikhail Belkin and Kaushik Sinha. Polynomial learning of distribution families. In 51th
Annual IEEE Symposium on Foundations of Computer Science, FOCS 2010, October 23-26,
2010, Las Vegas, Nevada, USA, pages 103–112, 2010.
[26] Shai Ben-David, Benny Chor, Oded Goldreich, and Michel Luby. On the theory of average
case complexity. Journal of Computer and system Sciences, 44(2):193–219, 1992.
[27] Shalev Ben-David and Lev Reyzin. Data stability in clustering: A closer look. Theor. Comput.
Sci., 558:51–61, 2014.
[28] Yonatan Bilu, Amit Daniely, Nati Linial, and Michael E. Saks. On the practically interesting
instances of MAXCUT. In 30th International Symposium on Theoretical Aspects of Computer
Science, STACS 2013, February 27 - March 2, 2013, Kiel, Germany, pages 526–537, 2013.
[29] Yonatan Bilu and Nathan Linial. Are stable instances easy? Combinatorics, Probability &
Computing, 21(5):643–660, 2012.
[30] Guy E. Blelloch and Kanat Tangwongsan. Parallel approximation algorithms for facilitylocation problems. In SPAA 2010: Proceedings of the 22nd Annual ACM Symposium on
Parallelism in Algorithms and Architectures, Thira, Santorini, Greece, June 13-15, 2010, pages
315–324, 2010.

33

[31] Vladimir Braverman, Adam Meyerson, Rafail Ostrovsky, Alan Roytman, Michael Shindler,
and Brian Tagiku. Streaming k-means on well-clusterable data. In Proceedings of the TwentySecond Annual ACM-SIAM Symposium on Discrete Algorithms, SODA 2011, San Francisco,
California, USA, January 23-25, 2011, pages 26–40, 2011.
[32] S. Charles Brubaker and Santosh Vempala. Isotropic PCA and affine-invariant clustering. In
49th Annual IEEE Symposium on Foundations of Computer Science, FOCS 2008, October
25-28, 2008, Philadelphia, PA, USA, pages 551–560, 2008.
[33] Jaroslaw Byrka, Thomas Pensyl, Bartosz Rybicki, Aravind Srinivasan, and Khoa Trinh. An improved approximation for k -median, and positive correlation in budgeted optimization. In Proceedings of the Twenty-Sixth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA
2015, San Diego, CA, USA, January 4-6, 2015, pages 737–756, 2015.
[34] Moses Charikar and Sudipto Guha. Improved combinatorial algorithms for facility location
problems. SIAM J. Comput., 34(4):803–824, 2005.
[35] Bernard Chazelle. The discrepancy method - randomness and complexity. Cambridge University Press, 2001.
[36] Michael B. Cohen, Sam Elder, Cameron Musco, Christopher Musco, and Madalina Persu.
Dimensionality reduction for k-means clustering and low rank approximation. In Proceedings
of the Forty-Seventh Annual ACM on Symposium on Theory of Computing, STOC 2015,
Portland, OR, USA, June 14-17, 2015, pages 163–172, 2015.
[37] Vincent Cohen-Addad, Philip N. Klein, and Claire Mathieu. Local search yields approximation schemes for k-means and k-median in euclidean and minor-free metrics. In IEEE 57th
Annual Symposium on Foundations of Computer Science, FOCS 2016, 9-11 October 2016,
Hyatt Regency, New Brunswick, New Jersey, USA, pages 353–364, 2016.
[38] Vincent Cohen-Addad and Claire Mathieu. Effectiveness of local search for geometric optimization. In 31st International Symposium on Computational Geometry, SoCG 2015, June
22-25, 2015, Eindhoven, The Netherlands, pages 329–343, 2015.
[39] David Cohen-Steiner, Pierre Alliez, and Mathieu Desbrun. Variational shape approximation.
ACM Trans. Graph., 23(3):905–914, 2004.
[40] Amin Coja-Oghlan. Graph partitioning via adaptive spectral techniques. Combinatorics,
Probability & Computing, 19(2):227–284, 2010.
[41] Anirban Dasgupta, John E. Hopcroft, Ravi Kannan, and Pradipta Prometheus Mitra. Spectral
clustering with limited independence. In Proceedings of the Eighteenth Annual ACM-SIAM
Symposium on Discrete Algorithms, SODA 2007, New Orleans, Louisiana, USA, January 7-9,
2007, pages 1036–1045, 2007.
[42] Sanjoy Dasgupta. Learning mixtures of gaussians. In 40th Annual Symposium on Foundations
of Computer Science, FOCS ’99, 17-18 October, 1999, New York, NY, USA, pages 634–644,
1999.
[43] Sanjoy Dasgupta and Yoav Freund. Random projection trees for vector quantization. IEEE
Trans. Information Theory, 55(7):3229–3242, 2009.

34

[44] Sanjoy Dasgupta and Leonard J. Schulman. A probabilistic analysis of EM for mixtures of
separated, spherical gaussians. Journal of Machine Learning Research, 8:203–226, 2007.
[45] Inderjit S. Dhillon, Yuqiang Guan, and Jacob Kogan. Iterative clustering of high dimensional
text data augmented by local search. In Proceedings of the 2002 IEEE International Conference
on Data Mining (ICDM 2002), 9-12 December 2002, Maebashi City, Japan, pages 131–138,
2002.
[46] Daniel B Dias, Renata CB Madeo, Thiago Rocha, Helton H Bı́scaro, and Sarajane M Peres.
Hand movement recognition for brazilian sign language: a study using distance-based neural
networks. In Neural Networks, 2009. IJCNN 2009. International Joint Conference on, pages
697–704. IEEE, 2009.
[47] Richard O Duda, Peter E Hart, et al. Pattern classification and scene analysis, volume 3.
Wiley New York, 1973.
[48] Dan Feldman and Michael Langberg. A unified framework for approximating and clustering
data. In Proceedings of the 43rd ACM Symposium on Theory of Computing, STOC 2011, San
Jose, CA, USA, 6-8 June 2011, pages 569–578, 2011.
[49] Bernd Fischer, Johann Schumann, Wray Buntine, and Alexander G Gray. Automatic derivation of statistical algorithms: The em family and beyond. In Advances in Neural Information
Processing Systems, pages 673–680, 2002.
[50] Ronald A Fisher. The use of multiple measurements in taxonomic problems. Annals of
eugenics, 7(2):179–188, 1936.
[51] Gereon Frahling and Christian Sohler. A fast k-means implementation using coresets. Int. J.
Comput. Geometry Appl., 18(6):605–625, 2008.
[52] Zachary Friggstad, Mohsen Rezapour, and Mohammad R. Salavatipour. Local search yields
a PTAS for k-means in doubling metrics. In IEEE 57th Annual Symposium on Foundations
of Computer Science, FOCS 2016, 9-11 October 2016, Hyatt Regency, New Brunswick, New
Jersey, USA, pages 365–374, 2016.
[53] Zachary Friggstad and Yifeng Zhang. Tight analysis of a multiple-swap heurstic for budgeted
red-blue median. In 43rd International Colloquium on Automata, Languages, and Programming, ICALP 2016, July 11-15, 2016, Rome, Italy, pages 75:1–75:13, 2016.
[54] Diptesh Ghosh. Neighborhood search heuristics for the uncapacitated facility location problem.
European Journal of Operational Research, 150(1):150 – 162, 2003. O.R. Applied to Health
Services.
[55] Diptesh Ghosh. Neighborhood search heuristics for the uncapacitated facility location problem.
European Journal of Operational Research, 150(1):150–162, 2003.
[56] Sudipto Guha and Samir Khuller. Greedy strikes back: Improved facility location algorithms.
J. Algorithms, 31(1):228–248, 1999.
[57] Sudipto Guha, Adam Meyerson, Nina Mishra, Rajeev Motwani, and Liadan O’Callaghan.
Clustering data streams: Theory and practice. IEEE Trans. Knowl. Data Eng., 15(3):515–
528, 2003.
35

[58] Anupam Gupta and Kanat Tangwongsan. Simpler analyses of local search algorithms for
facility location. CoRR, abs/0809.2554, 2008.
[59] Venkatesan Guruswami and Piotr Indyk. Embeddings and non-approximability of geometric problems. In Proceedings of the Fourteenth Annual ACM-SIAM Symposium on Discrete
Algorithms, January 12-14, 2003, Baltimore, Maryland, USA., pages 537–538, 2003.
[60] Pierre Hansen and Nenad Mladenovic. J-means : a new local search heuristic for minimum
sum of squares clustering. Pattern Recognition, 34(2):405–413, 2001.
[61] Pierre Hansen and Nenad Mladenović. Variable neighborhood search: Principles and applications. European journal of operational research, 130(3):449–467, 2001.
[62] Sariel Har-Peled and Akash Kushal. Smaller coresets for k-median and k-means clustering.
Discrete & Computational Geometry, 37(1):3–19, 2007.
[63] 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.
[64] Kamal Jain, Mohammad Mahdian, and Amin Saberi. A new greedy approach for facility
location problems. In Proceedings on 34th Annual ACM Symposium on Theory of Computing,
May 19-21, 2002, Montréal, Québec, Canada, pages 731–740, 2002.
[65] Kamal Jain and Vijay 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.
[66] Ragesh Jaiswal and Nitin Garg. Analysis of k-means++ for separable data. In Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques - 15th
International Workshop, APPROX 2012, and 16th International Workshop, RANDOM 2012,
Cambridge, MA, USA, August 15-17, 2012. Proceedings, pages 591–602, 2012.
[67] W. B. Johnson and J. Lindenstrauss. Extensions of Lipschitz mapping into Hilbert space.
In Conf. in modern analysis and probability, volume 26 of Contemporary Mathematics, pages
189–206. American Mathematical Society, 1984.
[68] Ravindran Kannan, Hadi Salmasian, and Santosh Vempala. The spectral method for general
mixture models. SIAM J. Comput., 38(3):1141–1156, 2008.
[69] Tapas Kanungo, David M. Mount, Nathan S. Netanyahu, Christine D. Piatko, Ruth Silverman,
and Angela Y. Wu. An efficient k-means clustering algorithm: Analysis and implementation.
IEEE Trans. Pattern Anal. Mach. Intell., 24(7):881–892, 2002.
[70] Tapas Kanungo, David M. Mount, Nathan S. Netanyahu, Christine D. Piatko, Ruth Silverman,
and Angela Y. Wu. A local search approximation algorithm for k-means clustering. Comput.
Geom., 28(2-3):89–112, 2004.
[71] Michael Kerber and Sharath Raghvendra. Approximation and streaming algorithms for projective clustering via random projections. In Proceedings of the 27th Canadian Conference
on Computational Geometry, CCCG 2015, Kingston, Ontario, Canada, August 10-12, 2015,
2015.
36

[72] Stavros G. Kolliopoulos and Satish Rao. A nearly linear-time approximation scheme for the
euclidean k-median problem. SIAM J. Comput., 37(3):757–782, June 2007.
[73] Madhukar R. Korupolu, C. Greg Plaxton, and Rajmohan Rajaraman. Analysis of a local
search heuristic for facility location problems. J. Algorithms, 37(1):146–188, 2000.
[74] Amit Kumar and Ravindran Kannan. Clustering with spectral norm and the k-means algorithm. In 51th Annual IEEE Symposium on Foundations of Computer Science, FOCS 2010,
October 23-26, 2010, Las Vegas, Nevada, USA, pages 299–308, 2010.
[75] Amit Kumar, Yogish Sabharwal, and Sandeep Sen. Linear-time approximation schemes for
clustering problems in any dimensions. J. ACM, 57(2), 2010.
[76] Shrinu Kushagra, Samira Samadi, and Shai Ben-David. Finding meaningful cluster structure
amidst background noise. In Algorithmic Learning Theory - 27th International Conference,
ALT 2016, Bari, Italy, October 19-21, 2016, Proceedings, pages 339–354, 2016.
[77] Shi Li and Ola Svensson. Approximating k-median via pseudo-approximation. In Symposium
on Theory of Computing Conference, STOC’13, Palo Alto, CA, USA, June 1-4, 2013, pages
901–910, 2013.
[78] Moshe Lichman. UCI machine learning repository, 2013.
[79] Meena Mahajan, Prajakta Nimbhorkar, and Kasturi R. Varadarajan. The planar k-means
problem is NP-hard. Theor. Comput. Sci., 442:13–21, 2012.
[80] Konstantin Makarychev and Yury Makarychev.
resilient problems.

Algorithms for stable and perturbation-

[81] Jirı́ Matousek. On approximate geometric k-clustering. Discrete & Computational Geometry,
24(1):61–84, 2000.
[82] Frank McSherry. Spectral partitioning of random graphs. In 42nd Annual Symposium on
Foundations of Computer Science, FOCS 2001, 14-17 October 2001, Las Vegas, Nevada, USA,
pages 529–537, 2001.
[83] Nimrod Megiddo and Kenneth J. Supowit. On the complexity of some common geometric
location problems. SIAM J. Comput., 13(1):182–196, 1984.
[84] Ramgopal R. Mettu and C. Greg Plaxton. The online median problem. SIAM J. Comput.,
32(3):816–832, 2003.
[85] Rafail Ostrovsky, Yuval Rabani, Leonard J. Schulman, and Chaitanya Swamy. The effectiveness of Lloyd-type methods for the k-means problem. J. ACM, 59(6):28, 2012.
[86] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion,
Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, Jake
Vanderplas, Alexandre Passos, David Cournapeau, Matthieu Brucher, Matthieu Perrot, and
Édouard Duchesnay. Scikit-learn: Machine learning in python. J. Mach. Learn. Res., 12:2825–
2830, November 2011.
[87] Gilles Pisier. The volume of convex bodies and Banach space geometry. Cambridge Tracts in
Mathematics. 94, 1999.
37

[88] Leonard J. Schulman. Clustering for edge-cost minimization (extended abstract). In Proceedings of the Thirty-Second Annual ACM Symposium on Theory of Computing, May 21-23,
2000, Portland, OR, USA, pages 547–555, 2000.
[89] Edward Snelson, Carl Edward Rasmussen, and Zoubin Ghahramani. Warped gaussian processes. Advances in neural information processing systems, 16:337–344, 2004.
[90] Daniel A. Spielman and Shang-Hua Teng. Smoothed analysis of algorithms: Why the simplex
algorithm usually takes polynomial time. J. ACM, 51(3):385–463, 2004.
[91] Minghe Sun. Solving the uncapacitated facility location problem using tabu search. Computers
& OR, 33:2563–2589, 2006.
[92] Dilek Tuzun and Laura I Burke. A two-phase tabu search approach to the location routing
problem. European journal of operational research, 116(1):87–99, 1999.
[93] Dilek Tüzün and Laura I. Burke. A two-phase tabu search approach to the location routing
problem. European Journal of Operational Research, 116(1):87–99, 1999.
[94] Santosh Vempala and Grant Wang. A spectral algorithm for learning mixture models. J.
Comput. Syst. Sci., 68(4):841–860, 2004.
[95] Yi Yang, Min Shao, Sencun Zhu, Bhuvan Urgaonkar, and Guohong Cao. Towards event source
unobservability with minimum network traffic in sensor networks. In Proceedings of the First
ACM Conference on Wireless Network Security, WISEC 2008, Alexandria, VA, USA, March
31 - April 02, 2008, pages 77–88, 2008.

38

102

2.2
2.0
1.6

ratio

beta

1.8
101

1.4
1.2
1.0
10-2

10-1

covariance

0
1010
-3

101

100

(a) k = 5, d = 2.

10-2

10-1
covariance

100

101

(b) k = 50, d = 2.

ratio

102

101

0
1010
-4

10-3

10-1
10-2
covariance

100

101

(c) k = 100, d = 2.

Figure 4: The ratio of the average k-means cost induced by the means over the average optimal
cost vs the variance for 2-dimensional instances generated from a mixture of k Gaussians (k ∈
{5, 50, 100}). We observe that the k-means objective becomes “relevant” (i.e., is less than 1.05
times the optimal value) for finding the clustering induced by Gaussians when the variance is less
than 0.1 for k = 5, less than 0.02 when k = 50, and less than 0.0005 when k = 100.

39

1.8

1.15

1.6

1.10

1.4

ratio

ratio

1.20

1.05

1.2

1.00
0.95
10-2

1.0
10-1

101
100
covariance

102

103

10-2

(a) k = 5, d = 10.

10-1

100
covariance

101

102

(b) k = 50, d = 10.

1.6
1.5

beta

1.4
1.3
1.2
1.1
1.0
10-2

10-1

covariance

100

101

(c) k = 100, d = 10.

Figure 5: The ratio of the average k-means cost induced by the means over the average optimal
cost vs the variance for 10-dimensional instances generated from a mixture of k Gaussians (k ∈
{5, 50, 100}). We observe that the k-means objective becomes “relevant” (i.e., is less than 1.05
times the optimal value) for finding the clustering induced by Gaussians when the variance is less
than 0.1 for k = 5, less than 0.02 when k = 50, and less than 0.0005 when k = 100.

40

100

101

10-1

100

10-2

10-1

10-3

beta

beta

102

10-2

10-4

10-3

10-5

-4
1010
-2

10-1

101
100
covariance

102

-6
1010
-2

103

(a) k = 5, d = 10.

10-1

100
covariance

101

102

(b) k = 50, d = 10.

100

beta

10-1
10-2
10-3
10-4
-5
1010
-2

10-1

covariance

100

101

(c) k = 100, d = 10.

Figure 6: The average minimum value of β for which the instance is β-distribution-stable vs the
variance for 10-dimensional instances generated from a mixture of k Gaussians (k ∈ {5, 50, 100}).
We observe that for relevant variances, the value of β is greater than 0.001.

41

100

10-2
10-3

10-2

beta

beta

10-1

10-3

10-4
10-5
10-6

-4
1010
-2

10-1

covariance

-7
1010
-3

101

100

(a) The average minimum value of β for which the
instances is β-distribution-stable vs the variance for 2dimensional instances generated from a mixture of 5
Gaussians.

10-2

10-1
covariance

100

101

(b) The average of the minimum value of β for which
the instances is β-distribution-stable vs the variance
for 2-dimensional instances generated from a mixture
of 50 Gaussians.

10-1
10-2
beta

10-3
10-4
10-5
10-6
-7
1010
-4

10-3

10-1
10-2
covariance

100

101

(c) The average minimum value of β for which the
instance is β-distribution-stable vs the variance for 2dimensional instances generated from a mixture of 100
Gaussians.

Figure 7: The average minimum value of β for which the instance is β-distribution-stable vs the
variance for 2-dimensional instances generated from a mixture of k Gaussians (k ∈ {5, 50, 100}).
We observe that for relevant variances, the value of β is greater than 0.001.

42

