Local Search Yields a PTAS for k-Means in Doubling Metrics
Zachary Friggstad†

Mohsen Rezapour

∗

Mohammad R. Salavatipour‡

arXiv:1603.08976v2 [cs.DS] 9 Jan 2017

Department of Computing Science
University of Alberta

Abstract
The most well known and ubiquitous clustering problem encountered in nearly every branch
of science is undoubtedly k-means: given a set of data points and a parameter k, select k centres
and partition the data points into k clusters around these centres so that the sum of squares of
distances of the points to their cluster centre is minimized. Typically these data points lie in
Euclidean space Rd for some d ≥ 2.
k-means and the first algorithms for it were introduced in the 1950’s. Over the last six
decades, hundreds of papers have studied this problem and different algorithms have been proposed for it. The most commonly used algorithm in practice is known as Lloyd-Forgy, which is
also referred to as “the” k-means algorithm, and various extensions of it often work very well
in practice. However, they may produce solutions whose cost is arbitrarily large compared to
the optimum solution. Kanungo et al. [2004] analyzed a very simple local search heuristic to get
a polynomial-time algorithm with approximation ratio 9 +  for any fixed  > 0 for k-means in
Euclidean space.
Finding an algorithm with a better worst-case approximation guarantee has remained one
of the biggest open questions in this area, in particular whether one can get a true PTAS for
fixed dimension Euclidean space. We settle this problem by showing that a simple local search
algorithm provides a PTAS for k-means for Rd for any fixed d.
More precisely, for any error parameter  > 0, the local search algorithm that considers
swaps of up to ρ = dO(d) · −O(d/) centres at a time will produce a solution using exactly k
centres whose cost is at most a (1 + )-factor greater than the optimum solution. Our analysis
extends very easily to the more general settings where we want to minimize the sum of q’th
powers of the distances between data points and their cluster centres (instead of sum of squares
of distances as in k-means) for any fixed q ≥ 1 and where the metric may not be Euclidean but
still has fixed doubling dimension.
Finally, our techniques also extend to other classic clustering problems. We provide the first
demonstration that local search yields a PTAS for uncapacitated facility location and the
generalization of k-median to the setting with non-uniform opening costs in doubling metrics.

1

Introduction

With advances in obtaining and storing data, one of the emerging challenges of our age is data analysis. It is hard to find a scientific research project which does not involve some form of methodology
∗

A preliminary version of this paper appeared in Proceedings of 57th FOCS 2016 [27]
This research was undertaken, in part, thanks to funding from the Canada Research Chairs program and an
NSERC Discovery Grant.
‡
Supported by NSERC.
†

1

to process, understand, and summarize data. A large portion of data analysis is concerned with
predicting patterns in data after being trained with some training data set (machine learning).
Two problems often encounterd in data analysis are classification and clustering. Classification
(which is an instance of supervised learning) is the task of predicting the label of a new data point
after being trained with a set of labeled data points (called a training set). Basically, after given
a training set of correctly labeled data points the program has to identify the label of a given new
(unlabeled) data point. Clustering (which is an instance of unsupervised learning) is the task of
grouping a given set of objects or data points into clusters/groups such that the data points that
are more similar fall into the same cluster while data points (objects) that do not seem similar are
in different clusters. Some of the main purposes of clustering are to understand the underlying
structure and relation between objects and find a compact representation of data points.
Clustering and different methods to achieve it have been studied since the 1950’s in different
branches of science: biology, statistics, medical sciences, computer science, social sciences, engineering, physics, and more. Depending on the notion of what defines a cluster, different models of
clustering have been proposed and studied by researchers. Perhaps the most widely used clustering
model is the k-means clustering: Given a set X of n data points in d-dimensional Euclidean space
Rd , and an integer k, find a set of k points c1 , . . . , ck ∈ Rd to act as as centres that minimize the
sum of squared distances of each data point to its nearest centre. In other words, we would like to
partition X into k cluster sets, {C1 , . . . , Ck } and find a centre ci for each Ci to minimize
k X
X

||x − ci ||22 .

i=1 x∈Ci

Here, ||x − ci ||2 is the standard Euclidean distance in Rd between points x and ci .
This value is called the cost of the clustering. Typically, the centres ci are selected to be the
centroid (mean) of the cluster Ci . In other situations the centres must be from the data points
themselves (i.e. ci ∈ Ci ) or from a given set C. This latter version is referred to as discrete k-means
clustering. Although in most application of k-means the data points are in some Euclidean space,
the discrete variant can be defined in general metrics. The k-means clustering problem is known
to be an NP-hard problem even for k = 2 or when d = 2 [1, 45, 23, 50].
Clustering, in particular the k-means clustering problem as the most popular model for it, has
found numerous applications in very different areas. The following is a (short) list of applications of
clustering that have been addressed by Jain [35]: image segmentation, information access, grouping customers into different types for efficient marketing, grouping delivery services for workforce
management and planning, and grouping genome data in biology. For instance, clustering is used
to identify groups of genes with related expression patterns in a key step of the analysis of gene
functions and cellular processes. It is also extensively used to group patients based on their genetic,
pathological, and cellular features which is proven useful in analyzing human genetics diseases (e.g.,
see [22, 33]).
The most widely used algorithm for k-means (which is also sometimes referred to as “the”
k-means algorithm) is a simple heuristic introduced by Lloyd in 1957 [44]. This algorithm starts
from an initial partition of the points into k clusters and it repeats the following two steps as long
as it improves the quality of the clustering: pick the centroids of the clusters as centres, and then
re-compute a new clustering by assigning each point to the nearest centre. Although this algorithm
works well in practice it is known that the ratio of the cost of the solution computed by this
algorithm vs the optimum solution cost (known as the “approximation ratio”) can be arbitrarily
large (see [36]). Various modifications and extensions of this algorithm have been produced and
studied, e.g. ISODATA, FORGY, Fuzzy C-means, k-means++, filtering using kd-trees (see [35]),
2

but none of them are known to have a bounded approximation ratio in the general setting. Arthur
and Vassilvitskii [6] show that Lloyd’s method with properly chosen initial centres will be an
O(log k)-approximation. Ostrovsky et al. [48] show that under some assumptions about the data
points the approximation ratio is bounded by a constant. The problem of finding an efficient
algorithm for k-means with a proven theoretical bound on the cost of the solution returned is
probably one of the most well studied problems in the whole field of clustering with hundreds of
research papers devoted to this.
Arthur and Vassilvitskii [5], Dastupta and Gupta [20], and Har-Peled and Sadri [32] study
convergence rate of Lloyd’s algorithm. In particular [5] show that it can be super-polynomial.
More recently, Vattani [51] shows that it can take exponential time even in two dimensions. Arthur
et al. [4] proved that Lloyd’s algorithm has polynomial-time smoothed complexity. Kumar and
Kannan [38], Ostrovsky et al. [48], and Awasthi et al. [9] gave empirical and theoretical evidence
for when and why the known heuristics work well in practice. For instance [9] show that when the
size of an optimum (k − 1)-means is sufficiently larger than the cost of k-means then one can get
a near optimum solution to k-means using a variant of Lloyd’s algorithm.
The k-means problem is known to be NP-hard [1, 23, 45]. In fact, the k-means problem is
NP-hard if d is arbitrary even for k = 2 [1, 23]. Also, if k is arbitrary the problem is NP-hard
even for d = 2 [45, 50]. However, the k-means problem can be solved in polynomial time by the
algorithm of [34] when both k and d are constant.
A polynomial-time approximation scheme (PTAS) is an algorithm that accepts an additional
parameter  > 0. It finds solutions whose cost is at most 1 +  times the optimum solution cost
and runs in polynomial time when  can be regarded as a fixed constant (i.e. O(nf () ) for some
2
function f ). Matoušek [46] gave a PTAS for fixed k, d with running time O(n(log n)k −2k d ). Since
then several other PTASs have been proposed for variant settings of parameters but all need k to
be constant [47, 11, 21, 31, 39, 40, 25, 30]. Bādoiu et al. [11] gave a PTAS for fixed k and any
O(1)
d with time O(2(k/) poly(d)n logk n). De la Vega et al. [21] proposed a (1 + )-approximation
3 8
O(1)
with running time O(2(k / )(ln(k/)) ln k dn logk n). This was improved to O(2(k/) dn) by Kumar
et al. [39, 40]. By building coresets of size O(k−d log n), Har-Peled and Mazumdar [31] presented
a (1 + )-approximation with time O(n + k k+1 −(2d+1)k logk+1 n logk 1 ). This was slightly improved
by Har-Peled and Kushal [30]. Feldman et al. [25] developed another (1 + )-approximation with
running time Õ(nkd+d·poly(k/)+2Õ(k/) ), where the Õ(.) notation hides polylog factors. Recently,
Bandyapadhyay and Varadarajan [12] presented a pseudo-approximation for k-means in fixeddimensional Euclidean space: their algorithm finds a solution whose cost is at most 1 +  times of
the optimum but might use up to (1 + ) · k clusters.
The result of Matoušek [46] also shows that one can select a set C of “candidate” centers in Rd
from which the k centres should be chosen from with a loss of at most (1 + ) and this set can be
computed in time O(n−d log(1/)). This reduces the k-means problem to the discrete setting where
along with X we have a set C of candidate centres and we have to select k centres from C. Kanungo
et al. [36] proved that a simple local search heuristic yields an algorithm with approximation ratio
9 +  for Rd . This remains the best known approximation algorithm with polynomial running time
for Rd . They also present an algorithm which is a hybrid of local search and Lloyd algorithm and
give empirical evidence that this works much better in practice (better than Lloyd’s algorithm) and
has proven constant factor ratio. For general metrics, Gupta and Tangwongsan [29] proved that
local search is a (25 + )-approximation. It was an open problem for a long time whether k-means
is APX-hard or not in Euclidean metrics. This was recently answered positively by Awasthi et
al. [10] where they showed that the problem is APX-hard in Rd but the dimension d used in the
proof is Ω(log n). Blomer et al. [13] have a nice survey of theoretical analysis of different k-means

3

algorithms.
k-median
Another very well studied problem that is also closely related to k-means is k-median. The only
difference is that the goal (objective function) in k-median is to minimize
P sum of distances,
Pk the
instead of sum of square of distances as in k-means, i.e. minimize
i=1
x∈Ci δ(x, ci ) where
δ(x, ci ) is the distance between x and ci . This problem occurs in operations research settings.
There are constant factor approximation algorithms for k-median in general metrics. The
simple local search (which swaps in and out a constant number of centres in each iteration) is
known to give a 3 +  approximation by Arya et al. [7, 8]. The current best approximation uses
different techniques and has an approximation ratio of 2.611 +  [43, 14]. The local search (9 + )approximation (for k-means) in [36] can be seen as an extension of the analysis in [7] for k-median.
One reason that analysis of k-means is more difficult is that the squares of distances do not
necessarily satisfy the triangle inequality. For instances of k-median on Euclidean metrics, Arora
et al. [3], building on the framework of Arora [2], gave the first PTAS. Kolliopoulos and Rao [37]
improved the time complexity and presented a PTAS for Euclidean k-median with time complexity
1
d−1
O(2O((1+log  )/) n log n log k). Such an approximation is also known as an efficient PTAS: the
running time of the (1 + )-approximation is of the form f () · poly(n) (in fixed-dimension metrics).
Uncapacitated facility location
The uncapacitated facility location problem is the same as k-median except instead of a
cardinality constraint bounding the number of open facilities, we are instead
P costs
P given opening
fi for each i ∈ C. The goal is to find a set of centers S ⊆ C that minimizes x∈X δ(x, S) + i∈S fi .
Currently the best approximation for uncapacitated facility location in general metrics is
a 1.488-approximation [42]. As with k-median, a PTAS is known for uncapacitated facility
location in constant-dimensional Euclidean metrics [3, 37], with the latter giving an efficient
PTAS.

1.1

Our result and technique

Although a PTAS for k-median in fixed dimension Euclidean space has been known for almost two
decades, getting a PTAS for k-means in fixed dimension Euclidean space has remained an open
problem. We provide a PTAS for this setting. We focus on the discrete case where we have to select
a set of k centres from a given set C. Our main result is to show that a simple local search heuristic
that swaps up to ρ = dO(d) · −O(d/) centres at a time and assigns each point to the nearest centre
is a PTAS for k-means in Rd .
A precise description of the algorithm is given in Section 2. At a high level, we start with any
set of k centres S ⊆ C. Then, while there is some other set of k centres S 0 ⊆ C with |S − S 0 | ≤ ρ
for some constant ρ such that S 0 is a cheaper solution than S, we set S ← S 0 . Repeat until S
cannot be improved any further. Each iteration takes |C|O(ρ) time, which is polynomial when ρ is
a constant. Such a solution is called a local optimum solution with respect to the ρ-swap heuristic.
We still have to ensure that the algorithm only iterates a polynomial number of times; a standard
modification discussed in Section 2 ensures this.
Recall that the doubling dimension of a metric space is the smallest τ such that any ball of
radius 2r around a point can be covered by at most 2τ balls of radius r. If the doubling dimension
can be regarded as a constant then we call the metric a doubling metric (as in [49]). The Euclidean
metric over Rd has doubling dimension O(d). Our analysis implies a PTAS for more general settings
where the data points are in a metric space with constant doubling dimension (described below)
4

and when the objective function is to minimize the sum of q’th power of the distances for some
fixed q ≥ 1.
Let ρ(, d) := dO(d) · O(−d/) . We will articulate the absolute constants suppressed by the O(·)
notation later on in our analysis.
Theorem 1 The local search algorithm that swaps up to ρ(, d) centres at a time is a (1 + )approximation for k-means in metrics with doubling dimension d.
Now consider the generalization where the objective function measures the sum of q’th power of
the distances between points and their assigned cluster centre, we call this `qq -norm k-clustering.
Here, we are given the points X in a metric space δ(·, ·) along with a set C of potential centres. We
are to select k centres from C and partition the points into k cluster sets C1 , . . . , Ck with each Ci
having a centre ci so as to minimize
k X
X
δ(x, ci )q .
i=1 x∈Ci

Note that the case of q = 2 is the k-means problem and q = 1 is k-median. We note that our
analysis extends to provide a PTAS for this setting when q is fixed. That is, Theorem 1 holds for
`qq -norm k-clustering for constants q ≥ 1, except that we require that the local search procedure
q
be the ρ0 -swap heuristic where ρ0 = dO(d) · (2q /)O(2 ·d/) .
Note that even for the case of k-median, this is the first PTAS for metrics with constant doubling
dimension. Also, while a PTAS was known for k-median for constant-dimensional Euclidean
metrics, determining if local search provided such a PTAS was an open problem. For example, [18]
shows that local search can be used to get a 1 +  approximation for k-median that uses up to
(1 + ) · k centres.
For the ease of exposition, we initially present the proof restricted to k-means in Rd . We explain
in Section 5.1 how to extend the analysis to doubling metrics and describe in Section 5.2 how this
can be easily extended to prove Theorem 1 when the objective is `qq -norm k-clustering.
As mentioned earlier, Awasthi et al. [10] proved that k-means is APX-hard for d = Ω(log n) and
they left the approximability of k-means for lower dimensions as an open problem. A consequence of
our algorithm is that one can get a (1 + )-approximation for k-means that runs in sub-exponential
time for values of d up to O(log n/ log log n). More specifically, for any given 0 <  < 1 and
d = σ log n/ log log n for sufficiently small absolute constant σ we get a (1 + )-approximation for kκ
means that runs in time O(2n ), for some constant κ = κ(σ) < 1; for d = O(log log n/ log log log n)
we get a quasi-polytime approximation scheme (QPTAS). Therefore, our result in a sense shows
that the requirement of [10] of d = Ω(log n) to prove APX-hardness of k-means is almost tight
σ
unless NP ⊆ DT IM E(2n ).
The notion of coresets and using them for finding faster algorithms for k-means has been studied
extensively (e.g. [31, 30, 15, 25, 26] and references there). A coreset is a small subset of data points
(possibly with weights associated to them) such that running the clustering algorithm on them
(instead of the whole data set) generates a clustering of the whole data set with approximately
good cost. In order to do this, one can go to the discrete case. To use coresets, we need to be able
to solve (discrete) k-means in the more general setting where each centre i ∈ C has an associated
weight w(i), and the cost of assigning a point j to i (if i is selected to be a centre) is w(i) · δ(i, j)2 .
Our local search algorithm works for this weighted setting as well. We will show how to use these
ideas to improve the running time of the local search algorithm.
Finally, we observe that our techniques easily extend to show a natural local-search heuristic
for uncapacitated facility location is a PTAS. In particular, for the same constant ρ as in
Theorem 1 (in fact, it can be slightly smaller) we consider the natural local search algorithm that
5

returns a solution S ⊆ C such that cost(S) ≤ cost(S 0 ) for all S 0 ⊆ C with both |S − S 0 | ≤ ρ and
|S 0 − S| ≤ ρ (i.e. add and/or drop up to ρ centres in C).
Theorem 2 The local search algorithm that adds and/or drops up to ρ(, d) centres at a time is a
(1+)-approximation for uncapacitated facility location in metrics with doubling dimension
d.
This seems to be the first explicit record of a PTAS for uncapacitated facility location
in doubling metrics, but one was known in constant-dimensional Euclidean metrics [3, 37]. Prior to
this work, local search was only known to provide a PTAS for uncapacitated facility location
in constant-dimensional Euclidean metrics if all opening costs are the same [19]. The basic idea
why this works is our analysis in Theorem 1 considers test swaps that, overall, swap in each local
optimum centre exactly once and swap out each global optimum centre exactly once (i.e. the swaps
come from a full partitioning of the local and global optimum). The partitioning we use to obtain
these swaps only uses the assumption that precisely k centres are open in a feasible solution at one
point, and this step can be safely ignored in the case of uncapacitated facility location.
Lastly, we consider the common generalization of k-median and uncapacitated facility
location where all centres have costs and at most k centres may be chosen (sometimes called
the generalized k-median or k-uncapacitated facility location problem). We consider
the local search operation that tries to add and/or drop up to ρ centres at a time as long as the
candidate solution being tested includes at most k centres. We get a PTAS in this case as well.
To the best of our knowledge, the previous best approximation was a 5-approximation in general
metrics, also obtained by local search [29].
Theorem 3 The local search algorithm that adds and/or drops up to ρ(, d) centres at a time (provided the resulting solution has at most k facilities) is a a (1 + )-approximation for generalized
k-median in metrics with doubling dimension d.
The results of Theorems 2 and 3 extend to the setting where we are considering the `qq -norm
for distances instead of the `1 -norm.
Note: Shortly after we announced our result [28], Cohen-Addad, Klein, and Mathieu [17, 16]
announced similar results for k-means on Euclidean and minor-free metrics using the local search
method. Specifically, they prove that the same local search algorithm yields a PTAS for k-means
on Euclidean and minor-free metrics. These results are obtained independently.

1.2

Proof Outline

The general framework for analysis of local search algorithms for k-median and k-means in [8, 36,
29] is as follows. Let S and O be a local optimum and a global optimum solution, respectively.
They carefully identify a set Q of potential swaps between local and global optimum. In each
such swap, the cost of assigning a data point x to the nearest centre after a swap is bounded
with respect to the local and global cost assignment. In other words,
Pif B(S) is the set of solutions
obtained by performing swaps from Q, the main task is to show that S 0 ∈B(S) (cost(S 0 )−cost(S)) ≤
α · cost(O) − cost(S) for some constant α. Given that 0 ≤ cost(S 0 ) − cost(S) for all S 0 ∈ B(S)
(because S is a local optimum), cost(S) ≤ α · cost(O).
Our analysis has the same structure but has many more ingredients and several intermediate
steps to get us what we want. Note that the following only describes steps used in the analysis of
the local search algorithm; we do not perform any of the steps described below in the algorithm
itself.

6

Let us define S and O as before. First, we do a filtering over S and O to obtain subsets S ⊆ S
and O ⊆ O such that every centre in S − S (in O − O) is “close” to a centre in S (in O) while
these filtered centres are far apart. We define a “net” around each centre i ∈ S which captures a
collection of other filtered centres in O that are relatively close to i. The idea of the net is that
if we choose to close i (in a test swap) then the data points that were to be assigned to i will be
assigned to a nearby centre in the net of i. Since the metric is a constant-dimensional Euclidean
metric (or, more generally, a doubling metric), we can choose these nets to have constant size.
For each j assigned to i in S, if the centre i∗ that j is assigned to in the optimum solution lies
somewhat close to i then we can reassign j to a facility in the net around i that is close to i∗ . In
this case, the reassignment cost for j will be close to c∗j − cj . Otherwise, if i∗ lies far from i then
we can reassign j to a facility near i in the net around i and the reassignment cost will only be
O() · (c∗j + cj ) and we will generate the c∗j − cj term for the local search analysis when i∗ is opened
in another different swap.
Of course, there are some complications in that we need to do something else with j if the net
around i is not open. This will happen infrequently, but we need a somewhat reasonable bound
when it does happen. Also, for reasons that will become apparent in the analysis we do something
different in the case that i∗ is somewhat close to i but i∗ is much closer to a different facility in S
than it is to i.
One main part of our proof is to show that there exists a suitable randomized partitioning of
S ∪ O such that each part has small size (which will be a function of only  and d and, ultimately,
determines the size of the swaps in our local search procedure), and for any pair (i, i∗ ) ∈ S ×O where
i∗ lies in the net of i we have Pr[i, i∗ lie in the same part] ≥ 1 − . This randomized partitioning is
the only part of the proof that we rely on properties of doubling metrics (or Rd ). For those small
portion of facilities that their net is “cut” by our partitioning, we show that when we close those
centres in our test swaps then the reassignment cost of a point j that was assigned to them is only
O(1) times more than the sum of their assignment costs in O and S. Given that this only happens
with small probability (due to our random partition scheme), this term is negligible in the final
analysis. So we get an overall bound of 1 + O() on the ratio of cost of S over O.
Outline of the paper: We start with some basic definitions and notation in Section 2. In
Sections 3 and 4 we show that the local search algorithm with an appropriate number of swaps
provides a PTAS for k-means in Rd . In Section 5 we show how this can be extended to prove
Theorem 1 and to the setting where we measure the `qq -norm of the solution for any constant q ≥ 1.
Finally, before concluding, in Section 6 we show that our analysis easily extends to prove Theorems
2 and 3.

2

Notation and Preliminaries

Recall that in the k-means problem we are given a set X of n points in Rd and an integer k ≥ 1;
we have to find k centres c1 , . . . , ck ∈ Rd so as to minimize the sum of squares of distances of each
point to the nearest centre. As mentioned earlier, by using the result of [46], at a loss of (1 + )
factor we can assume we have a set C of “candidate” centres from which the centres can be chosen
from. This set can be computed in time O(n−d log(1/)) and |C| = O(n−d log(1/)). Therefore,
we can reduce the problem to the discrete case.
Formally, suppose we are given a set C of points (as possible cluster centres) along with X and
we have to select the k centres from C. Furthermore, we assume the points are given in a metric
space (V, δ) (not necessarily Rd ). For any two points p, q ∈qV , δ(p, q) denotes the distance between
Pd
2
them: for the case of the metric being Rd , then δ(p, q) =
`=1 |p` − q` | .
7

Algorithm 1 ρ-Swap Local Search
Let S be an arbitrary set of k centres from C
while ∃ sets P ⊆ C − S, Q ⊆ S with |P | = |Q| ≤ ρ s.t. cost((S − Q) ∪ P ) < cost(S) do
S ← (S − Q) ∪ P
return S
We usually refer to a potential centre in C by a simple index i and a point in X by a simple
index j (or slight variants like i∗ or i0 ). This is to emphasize that we do not need to talk about
specific coordinates of points in Euclidean space. In fact, only once in our proof do we rely on the
particular embedding of the points in Euclidean space. This argument will also be replaced by a
more general argument when discussing doubling metrics in Section
P 5.1. So, for any set S ⊆ C and
any j ∈ X , let δ(j, S) = mini∈S δ(j, i). We also define cost(S) = j∈X δ(j, S)2 .
Our goal in (discrete) k-means is to find a set of centres S ⊆ C of size k to minimize
cost(S).
P
Note that once we fix the set of centres we can find a partitioning of X that realizes j∈X δ(j, S)2
by assigning each j ∈ X to the nearest centre in S, breaking ties arbitrarily.
For ease of exposition we focus on k-means on Rd and then show how the analysis can be
extended
to work for metrics with constant doubling dimension and when we want to minimize
P
q
j∈X δ(j, S) for a fixed q ≥ 1 instead of just q = 2.
The simple ρ-swap local search heuristic shown in Algorithm 1 is essentially the same one
considered in [36].
Recall that we defined ρ(, d) = dO(d) · O(d/) , where the constants will be specified later and
consider the local search algorithm with ρ = ρ(, d) swaps. By a standard argument (as in [7, 36])
one can show that replacing the condition of the while loop with cost((S −Q)∪P ) ≤ (1− k )·cost(S),
the algorithm terminates in polynomial time. Furthermore, if α is such that any locally optimum
solution returned by Algorithm 1 has cost at most α · cost(O) where O denotes a global optimum
solution, then any S such that cost((S − Q) ∪ P ) < (1 − k ) · cost(S) for any possible swap P, Q
α
satisfies cost(S) ≤ 1−
cost(O). This follows by arguments in [7, 36] and the fact that our local
search analysis uses at most k “test swaps”.
For ease of exposition, we ignore this factor 1 +  loss, and consider the solution S returned
by Algorithm 1. Recall that we use O to denote
the global optimum
P
P solution. For j ∈ X , let
∗
2
2
∗
cj = δ(j, O) and cj = δ(j, S) , so cost(O) = j∈X cj and cost(S) = j∈X cj . We also denote the
centre in O nearest to j by σ ∗ (j) and the centre in S nearest to j by σ(j). Define φ : O ∪ S → O ∪ S
to be the function that assigns i∗ ∈ O to its nearest centre in S and assigns i ∈ S to its nearest
centre in O. For any two sets S, T ⊆ O ∪ S, we let S4T = (S ∪ T ) − (S ∩ T ).
We assume O ∩ S = ∅. This is without loss of generality because we could duplicate each
location in C and say O uses the originals and S the duplicates. It is easy to check that S would
still be a locally optimum solution in this instance. We can also assume that these are the only
possible colocated facilities, so δ(i, i0 ) > 0 for distinct i, i0 ∈ O or distinct i, i0 ∈ S. Finally, we will
assume  is sufficiently small (independent of all other parameters, including d) so that all of our
bounds hold.

3

Local Search Analysis for Rd

In this section we focus on Rd (for fixed d ≥ 2) and define ρ(, d) = 32 · (2d)8d · −36·d/ . Our goal in
this section is to prove that the ρ-swap local search with ρ = ρ(, d) is a PTAS for k-means in Rd .

8

Algorithm 2 Filtering O
O←∅
for each i∗ ∈ O in nondecreasing order of Di∗ do
if ∃ i∗ ∈ O such that δ(i∗ , i∗ ) ≤  · Di∗ then
η(i∗ ) ← i∗
else
η(i∗ ) ← i∗
O ← O ∪ {i∗ }
return O
Theorem 4 Let S be a locally-optimum solution with respect to the ρ(, d)-swap local search heuristic when the points lie in Rd . Then cost(S) ≤ (1 + O()) · cost(O).
To prove this, we will construct a set of test swaps that yield various inequalities which, when
combined, provide the desired bound on cost(S). That is, we will partition O ∪ S into sets where
|P ∩ O| = |P ∩ S| ≤ ρ(, d) for each part P . For each such set P , 0 ≤ cost(S4P ) − cost(S) because
S is a locally optimum solution. We will provide an explicit upper bound on this cost change that
will reveal enough information to easily conclude cost(S) ≤ (1 + O()) · cost(O). For example, for
a point j ∈ X if σ ∗ (j) ∈ P then the change in j’s assignment cost is at most c∗j − cj because we
could assign j from σ(j) to σ ∗ (j). The problem is that points j with σ(j) ∈ P but σ ∗ (j) 6∈ P must
go somewhere else; most of our effort is ensuring that the test swaps are carefully chosen so such
reassignment cost increases are very small.
First we need to describe the partition of O ∪ S. This is a fairly elaborate scheme that involves
several steps. As mentioned earlier, the actual algorithm for k-means is the simple local search
we described and the algorithms we describe below to get this partitioning scheme are only for the
purpose of proof and analysis of the local search algorithm.
Definition 1 For i∗ ∈ O let Di∗ := δ(i∗ , S) = δ(i∗ , φ(i∗ )). For i ∈ S let Di := δ(i, O) = δ(i, φ(i)).
The first thing is to sparsify O and S using a simple filtering step. Algorithm 2 filters O to a set
that is appropriately sparse for our analysis.
Think of η(i∗ ) as a proxy for i∗ that is very close to i∗ . Using a similar process, we filter S to
get S and proxy centres η(i) ∈ S for each i ∈ S. The idea is that the set of centres left in O and S
are somewhat far apart yet any point that was assigned to a centre in O − O (or in S − S) can be
“cheaply” reassigned to a proxy.
Lemma 1 For each i ∈ O ∪ S we have δ(i, η(i)) ≤  · Di . For any distinct i, i0 ∈ O ∪ S we have
δ(i, i0 ) ≥  · max{Di , Di0 }.
Proof. That δ(i, η(i)) ≤  · Di follows immediately by construction. If i ∈ O, i0 ∈ S or vice-versa,
then in fact δ(i, i0 ) ≥ max{Di , Di0 } simply by definition of Di , Di0 .
Now suppose i, i0 ∈ O and that i0 was considered after i in Algorithm 2 (so Di0 ≥ Di ). The fact
that i0 was added to O even though i was already in O means δ(i, i0 ) ≥  · Di0 . The same argument
works if i, i0 ∈ S.
Next we define mappings similar to φ, σ, σ ∗ except they only concern centres that were not
filtered out.

9

i∗
4

Di∗
1

i∗
1
Di∗
1

Di

i∗ = cent(i1 )
3

i2 = φ(i∗ )
4

i1 = φ(i∗ ) = φ̄(i∗ )
2
2

i∗ = η(i∗ )
2
1

i3 = φ̄(i∗ ) = η(i2 )
4

Figure 1: The square centres lie in O and the circle centres lie in S. The black ones survived the
filtering (i.e. lie in O∪S) and the grey ones were filtered out. Note that the grey square i∗1 on the left
was too close to i∗2 so it was not added to O and η(i∗1 ) = i∗2 . Also note that φ(i∗4 ) 6= φ(i∗4 ) because the
facility i2 that defined Di∗4 := δ(i∗4 , i2 ) was not added to S. Still, δ(i∗4 , φ(i∗4 )) = δ(i∗4 , i3 ) ≤ (1 + )Di∗4 .
Definition 2
• φ : O ∪ S → O ∪ S maps each i ∈ O to its nearest location in S and vice versa.
• σ ∗ : X → O defined by σ ∗ (j) = η(σ ∗ (j)).
• σ : X → S defined by σ(j) = η(σ(j)).
Finally, for each i ∈ φ(O), let cent(i) be the centre in φ
arbitrarily.

−1

(i) that is closest to i, breaking ties

Note σ(j) may not necessarily be the centre in S that is closest to j. Also note that if one considers
a bipartite graph with parts O and S, then φ maps centres from one side to the other.
Lemma 2 For each i0 ∈ O ∪ S, Di0 ≤ δ(i0 , φ(i0 )) ≤ (1 + ) · Di0 .
Proof. Suppose i0 ∈ O, the proof is essentially the same for i0 ∈ S. On one hand, we know
Di0 = δ(i0 , φ(i0 )) ≤ δ(i0 , φ(i0 ))
because S ⊆ S. On the other hand,
δ(i0 , φ(i0 )) = δ(i0 , η(φ(i0 ))) ≤ δ(i0 , φ(i0 )) + δ(φ(i0 ), η(φ(i0 ))) ≤ Di0 +  · Dφ(i0 ) .
Conclude by observing Dφ(i0 ) ≤ δ(i0 , φ(i0 )) = Di0 .
Figure 1 depicts many of the concepts covered above. Finally, the last definition in this section
identifies pairs of centres that we would like to have in the same part of the partition we construct.
Definition 3
• T := {(cent(i), i) : i ∈ φ(O) and  · δ(cent(i), i) ≤ Di }
• N := {(i∗ , i) ∈ O × S : δ(i, i∗ ) ≤ −1 · Di and Di∗ ≥  · Di }
For each i ∈ S, the set {i∗ : (i∗ , i) ∈ N } is the “net” for centre i that was discussed in the proof
outline in Section 1.2.
Ultimately we will require that pairs in T are not separated by the partition. Our requirement
for N is not quite as strong. The partition is constructed randomly and it will be sufficient to have
each pair in N being separated by the partition with probability at most .
10

3.1

Every Centre is Close to Some Pair in T

The following says that if at least one centre of each pair in T is open after a swap, then every
centre in O ∪ S is somewhat close to some open centre. The bound is a bit big, but it will be
multiplied by O() whenever it is used in the local search analysis.
Lemma 3 Let A ⊆ O ∪ S be such that A ∩ {cent(i), i} 6= ∅ for each (cent(i), i) ∈ T . Then
δ(i0 , A) ≤ 5 · Di0 for any i0 ∈ O ∪ S.
Proof. We first prove the statement for i0 ∈ O, the other case is similar but requires one additional
step so we will discuss it below. Consider the following sequence of centres. Initially, set i0 := i0 ,
i1 := η(i0 ), and i2 := φ(i1 ). We build the rest inductively, noting that we guarantee ia ∈ φ(O) for
even indices a ≥ 2 (so cent(ia ) is defined).
Inductively, for even a ≥ 2 we do the following. If ia ∈ A then we stop. Otherwise, if
(cent(ia ), ia ) ∈ T then by assumption it must be that cent(ia ) ∈ A so we let ia+1 := cent(ia )
and stop. Finally, if (cent(ia ), ia ) 6∈ T then we set ia+1 := φ(ia ) and ia+2 := φ(ia+1 ) and iterate
with a0 = a + 2. This walk is depicted in Figure 2.
We will soon show the walk terminates. For now, we observe that, apart from the first step,
the steps decrease in length geometrically. In particular, consider some a ≥ 2 such that the walk
did not stop at ia . If ia+1 = φ(ia ) then
δ(ia , ia+1 ) = δ(ia , φ(ia )) ≤ δ(ia , ia−1 ).
If ia+1 6= φ(ia ) then it must be ia = cent(ia−1 ). In this case, it must be ia = φ(ia−1 ), so because
−1
cent(ia ) is the closest centre in φ (ia ) to ia (by definition), we have
δ(ia , ia+1 ) = δ(ia , cent(ia )) ≤ δ(ia , ia−1 ).
We prove that the lengths of the edges traversed decrease geometrically with every other step.
This, along with the fact that the lengths of the steps of the walk are nonincreasing (except,
perhaps, the first two steps), will show that this process eventually terminates and also bounds the
cost of the path.
Claim 1 For every even a ≥ 2 such that the walk did not end at ia or ia+1 , we have δ(ia , ia+1 ) ≤
2 · δ(ia−1 , ia ).
Proof. Because the walk did not end at ia or ia+1 , (cent(ia ), ia ) 6∈ T meaning Dia <  ·
δ(cent(ia ), ia ) ≤  · δ(ia−1 , ia ). Using this and Lemma 2 in the first bound below, we see
δ(ia , ia+1 ) ≤ (1 + ) · Dia < (1 + ) · δ(ia−1 , ia ) ≤ 2 · δ(ia−1 , ia )

11

i3 = φ̄(i2 )
i0 = i0

i5 = φ̄(i4 )
i4 = φ̄(i3 )

≤ Di0
≤ (1 + )Di0

i2 = φ̄(i1 )

i1 = η(i0 )

Figure 2: Illustration of the walk i0 , i1 , . . . in the proof of Lemma 3. Apart from a = 1, the step
(ia , ia+1 ) is shorter than the step (ia−1 , ia ). Furthermore, every other step decreases in length
geometrically.
Let m be the index of the last centre in the walk. Thus, δ(ia , ia+1 ) ≤ δ(ia−1 , ia ) for all 2 ≤ a ≤
m − 1. From this and Claim 1 we have
δ(i0 , A) ≤

m−1
X

δ(ij , ij+1 )

j=0

≤ δ(i0 , i1 ) +

X

2δ(ia−1 , ia ) + δ(im−1 , im )

2≤a≤m
a even

≤ δ(i0 , i1 ) + 2δ(i1 , i2 )

X
(2)a + δ(im−1 , im )
a≥0

2
· δ(i1 , i2 ) + δ(i1 , i2 )
1 − 2
3 − 2
= δ(i0 , i1 ) +
· δ(i1 , i2 )
1 − 2
(3 − 2)(1 + )
· Dη(i0 )
≤  · D i0 +
1 − 2
≤ (3 + 7) · Di0 .
≤ δ(i0 , i1 ) +

The second last step uses Lemma 2 and the last step uses the fact that Dη(i0 ) ≤ Di0 (either i0 = η(i0 )
or else i0 was filtered out by η(i0 ), in which case it has a larger D-value) and the assumption that
 is small enough.
Now suppose i0 ∈ S. We bound δ(i0 , A) mostly using what we have done already. That is, we
have
δ(i0 , A) ≤ δ(i0 , φ(i0 )) + δ(φ(i0 ), A) ≤ Di0 + (3 + 7)Dφ(i0 ) .
Note Dφ(i0 ) ≤ δ(i0 , φ(i0 )) = Di0 again by Lemma 2. So,
δ(i0 , A) ≤ (4 + 7)Di0 ≤ 5Di0 .

3.2

Good Partitioning of O ∪ S and Proof of Theorem 4

The main tool used in our analysis is the existence of the following randomized partitioning scheme.

12

Theorem 5 There is a randomized algorithm that samples a partitioning π of O ∪ S such that:
• For each part P ∈ π, |P ∩ O| = |P ∩ S| ≤ ρ(, d).
• For each part P ∈ π, S4P includes at least one centre from every pair in T .
• For each (i∗ , i) ∈ N , Pr[i, i∗ lie in different parts of π] ≤ .
We prove this theorem in Section 4. For now, we will complete the analysis of the local search
algorithm using this partitioning scheme. Note that in the following we do not use the geometry of
the metric (i.e. all arguments hold for general metrics); it is only in the proof of Theorem 5 that
we use properties of Rd .
The following gives a way to handle the fact that the triangle inequality does not hold with
squares of the distances.
Lemma 4 For any real numbers x, y we have (x + y)2 ≤ 2(x2 + y 2 ).
Proof. (x + y)2 ≤ (x + y)2 + (x − y)2 = 2x2 + 2y 2 .
Lemma 5 For each point j ∈ X , Dσ(j) ≤ Dσ(j) ≤ δ(j, σ(j)) + δ(j, σ ∗ (j)). Similarly, Dσ∗ (j) ≤
Dσ∗ (j) ≤ δ(j, σ(j)) + δ(j, σ ∗ (j)).
Proof. As usual, we only prove the first statement since the second is nearly identical. If σ(j) =
σ(j) then Dσ(j) = Dσ(j) is trivially true. Otherwise, σ(j) = η(σ(j)) was already in S when σ(j)
was considered by the filtering algorithm meaning Dσ(j) ≤ Dσ(j) .
For the other inequality, note
Dσ(j) = δ(σ(j), φ(σ(j))) ≤ δ(σ(j), σ ∗ (j)) ≤ δ(j, σ(j)) + δ(j, σ ∗ (j)).
Proof of Theorem 4.
Let π be a partition sampled by the algorithm from Theorem 5. For each point j ∈ X and each
part P of π, let ∆Pj := δ(j, S4P )2 − δ(j, S)2 denote the change in assignment cost for the point
after swapping in the centers in P ∩ P
O and swapping out P ∩ S. Local optimality of S and
|P ∩ S| = |P ∩ O| ≤ ρ(, d) means 0 ≤ j ∆Pj for any part P .
Classify each point j ∈ X in one of the following ways:
• Lucky: σ(j) and σ(j) do not lie in the same part of π.
• Long: j is not lucky but δ(σ(j), σ ∗ (j)) > −1 · Dσ(j) .
• Bad: j is not lucky or long and (σ ∗ (j), σ(j)) ∈ N yet σ(j), σ ∗ (j) lie in different parts of π.
• Good: j is neither lucky, long, nor bad.
P
We now place an upper bound on P ∈π ∆Pj for each point j ∈ X . Note that each centre in
S is swapped out exactly once over all swaps P and each centre in O is swapped in exactly once.
With this in mind, consider the following cases for a point j ∈ X . In the coming arguments, we let
δj := δ(j, σ(j)) and δj∗ := δ(j, σ ∗ (j)) for brevity. Note cj = δj2 and c∗j = δj∗2 .
In all cases for j except when j is bad, the main idea is that we can bound the distance from
j to some point in S4P by first moving it to either σ(j) or σ ∗ (j) and then moving it a distance

13

of O() · (δj + δj∗ ) to reach an open facility. Considering that we reassigned j from σ(j), the
reassignment cost will be
(δj + O() · (δj + δj∗ ))2 − cj = O() · (cj + c∗j )
or
(δj∗ + O() · (δj + δj∗ ))2 − cj = (1 + O()) · c∗j − (1 − O()) · cj .
Case: j is lucky
For the part P ∈ π with σ ∗ (j) ∈ P , we have ∆Pj ≤ c∗j − cj as we could move j from σ(j) to σ ∗ (j).
If σ(j) is swapped out in a different swap P 0 , we move j to σ(j) (which remains open because j is
0
lucky) and bound ∆Pj by:
∆Pj

0

≤ δ(j, σ(j))2 − δj2
≤ (δj + δ(σ(j), σ(j)))2 − cj
≤ (δj +  · Dσ(j) )2 − cj
2

= 2 · δj · Dσ(j) +  ·

(Lemma 1)

2
Dσ(j)
2

≤ 2 · δj · (δj + δj∗ ) +  · (δj + δj∗ )2
≤ 2 · (δj +

δj∗ )2

2

+  · (δj +

(Lemma 5)

δj∗ )2

≤ 4 · (c∗j + cj ) + 22 · (c∗j + cj )
≤ 6 ·

(c∗j

(Lemma 4)

+ cj ),

again using the assumption that  is sufficiently small. For every other swap P 00 , we have that σ(j)
00
remains open after the swap so ∆Pj ≤ 0 as we could just leave j at σ(j). In total, we have
X

∆Pj ≤ c∗j − cj + 6 · (c∗j + cj ).

P ∈π

Case: j is long
Again, for P ∈ π with σ ∗ (j) ∈ P we get ∆Pj ≤ c∗j − cj . If σ(j) is swapped out in a different swap
0
P 0 , then we bound ∆Pj by moving j from σ(j) to the open centre nearest to σ(j). Note that S4P 0
contains at least one centre from every pair in T , so we bound this distance using Lemma 3. This
case is depicted in Figure 3. We have
Dσ(j) ≤ δ(σ(j), σ(j)) + δ(σ(j), φ(σ(j)))
≤ Dσ(j) + Dσ(j)

(Lemma 1)
∗

∗

≤ (δ(j, σ(j)) + δ(j, σ (j)) + δ(σ(j), σ (j))
≤ (δj +
≤ (δj +
≤ (δj +

δj∗ )
δj∗ )
δj∗ )

(since j is long)
∗

+ (δ(σ(j), σ(j)) + δ(σ(j), j) + δ(j, σ (j)) + δ(σ ∗ (j), σ ∗ (j)))
+ (Dσ(j) + δj + δj∗ + Dσ∗ (j) )
+ ((δj +

= 2(1 + )(δj +

δj∗ )

+ δj +

δj∗ ).

14

δj∗

+ (δj +

(Lemma 1)
δj∗ ))

(Lemma 5)

0

Using this, we bound ∆Pj as follows.
∆Pj

0

≤ (δj + δ(σ(j), S4P ))2 − cj
≤ (δj + 5Dσ(j) )2 − cj
≤
≤
≤
≤

(Lemma 3)

(δj + 10(1 + ) · (δj + δj∗ ))2 − cj
21 · δj (δj + δj∗ ) + 1012 · (δj + δj∗ )2
22 · (δj + δj∗ )2
44 · (cj + c∗j ).

bound for Dσ(j) above

i = σ̄(j)
σ(j)

Di


Di

σ̄ ∗ (j)
ī ∈ S4P

j

Figure 3: Illustrating the case when j is long. It is moved from σ(j) to the nearest open centre at
an additional distance of O(Dσ(j) ). This is negligible compared to δj + δj∗ .
In every other swap we could leave j at σ(j). Thus, for a long point j we have
X
∆Pj ≤ c∗j − cj + 44 · (c∗j + cj ).
P ∈π

Case: j is bad
We only move j in the swap P when σ(j) is closed. In this case, we assign j in the same way as if it
was long. Our bound is weaker here and introduces significant positive dependence on cj . This will
eventually be compensated by the fact that j is bad with probability at most  over the random
choice of π. For now, we just provide the reassignment cost bound for bad j.
∆Pj

≤ (δj + δ(σ(j), S4P ))2 − cj
≤ (δj + 5Dσ(j) )2 − cj
5δj∗ )2

≤ (6δj +
≤ 71 ·

(c∗j

− cj

(Lemma 4)

+ cj ).

So for bad points we have
X

(Lemma 3)

∆Pj ≤ 71 · (c∗j + cj ).

P ∈π

15

Case: j is good
This breaks into two subcases. We know δ(σ ∗ (j), σ(j)) ≤ −1 · Dσ(j) because j is not long. In
one subcase, Dσ∗ (j) ≥  · Dσ(j) so (σ ∗ (j), σ(j)) ∈ N . Since j is not bad and not lucky, we have
σ(j), σ(j), σ ∗ (j) ∈ P for some common part P ∈ π. In the other subcase, Dσ∗ (j) <  · Dσ(j) . Note
in this case we still have σ(j), σ(j) ∈ P for some common part P because j is not lucky.
Subcase: Dσ∗ (j) ≥  · Dσ(j)
The only time we move j is when σ(j) is closed. As observed in the previous paragraph, this
happens in the same swap when σ ∗ (j) is opened, so send j to σ ∗ (j). This is illustrated in Figure 4.

σ ∗ (j)

σ¯∗ (j) ∈ S4P
i = σ̄(j)
σ(j)
Di

Di


j

Figure 4: Illustrating the subcase when Dσ∗ (j)) ≥  · Dσ(j) . In this case, the only additional distance
j travels after first being reassigned to σ ∗ (j) is to the nearby proxy σ ∗ (j)

∆Pj

≤ (δj∗ + δ(σ ∗ (j), σ ∗ (j)))2 − cj
≤ (δj∗ +  · Dσ∗ (j) )2 − cj
≤
=
≤
≤

(Lemma 1)

(δj∗ +  · (δj + δj∗ ))2 − cj
c∗j + 2 · (δj + δj∗ )2 + 2 · δj∗ (δj + δj∗ ) − cj
c∗j + 22 · (cj + c∗j ) + 2 · (δj + δj∗ )2 − cj
c∗j − cj + 5 · (c∗j + cj ).

(Lemma 4)
(Lemma 4)

Subcase: Dσ∗ (j) <  · Dσ(j)
Again, the only time we move j is when σ(j) is closed. We reassign j by first moving it to σ ∗ (j)
and then using Lemma 3 to further bound the cost. Figure 5 depicts this reassignment.
We bound the cost change for this reassignment as follows. Recall that δ(σ ∗ (j), σ ∗ (j)) ≤ Dσ∗ (j)
by our filtering. This implies that Dσ∗ (j) ≤ Dσ∗ (j) + Dσ∗ (j) which in turn implies Dσ∗ (j) ≤
1
∗
∗
1− Dσ (j) ; thus Dσ ∗ (j) ≤ (1 + 2)Dσ (j) which in turn is bounded by (1 + )Dσ(j) in this subcase

16

ī ∈ S4P
ī∗ = σ̄ ∗ (j)
j

i = σ̄(j)
Di


σ(j)
Di

Figure 5: Illustrating the subcase when Dσ∗ (j) <  · Dσ(j) . The only additional distance j moves
after being reassigned to σ ∗ (j) is small, as it travels first to the nearby proxy σ ∗ (j) and then at
distance at most 5Dσ∗ (j) ≤ 5 · Dσ(j) by Lemma 3.
(by the assumption of subcase). Thus,
∆Pj

≤ (δj∗ + δ(σ ∗ (j), S4P ))2 − cj
≤ (δj∗ + 5Dσ∗ (j) ) − cj
≤ (δj∗ + 5(1 + 2) · Dσ(j) )2 − cj
≤ (δj∗ + 5(1 + 2) · (δj + δj∗ ))2 − cj
≤
≤
≤

c∗j
c∗j
c∗j

+ 11 ·

δj∗ (δj

+

δj∗ )

− cj + 12 · (δj +

2

+ 26 (δj +

(Lemma 5)
δj∗ )2

− cj

δj∗ )2

− cj + 24 · (c∗j + cj ).

Considering both subcases we can say that for any good point j that
X
∆Pj ≤ c∗j − cj + 24 · (c∗j + cj ).
P ∈π

P
Aggregating these bounds and remembering 0 ≤ j∈X ∆Pj for each P ∈ π because S is a locally
optimum solution, we have
XX
X 
 X
0≤
∆Pj ≤
(1 + 44)c∗j − (1 − 44)cj +
71(c∗j + cj ).
j∈X P ∈π

j∈X
j not bad

j∈X
j bad

The last step is to average this inequality over the random choice of π. Note that any point
j ∈ X is bad with probability at most  by the guarantee in Theorem 5 and the definition of bad.

17

Thus, we see

0 ≤ Eπ 


XX

∆Pj 

j∈X P ∈π

≤

X



Pr[j not bad] · (1 + 44)c∗j − (1 − 44)cj + Pr[j bad] · 71(c∗j + cj )

j∈X

≤

X

(1 + 44)c∗j − (1 − ) · (1 − 44)cj + 71 · (c∗j + cj )

j∈X

≤

X

(1 + 115)c∗j − (1 − 116)cj .

j∈X

Rearranging shows
cost(S) ≤

3.3

1 + 115
· cost(O) ≤ (1 + O()) · cost(O).
1 − 116

Running Time Analysis

Recall that we can go to the discrete case by finding a set C of size O(n−d log(1/). The analysis of
1
Arya et al. [7, 8] shows that the number of local search steps is at most log(cost(S0 )/cost(O))/ log 1−/k
,
where S0 is the initial solution. This is polynomial in the total bit complexity of the input (i.e.
the input size). So we focus on bounding the time complexity of each local search step. Since
the number of swaps in each step is bounded by ρ = ρ(, d), a crude upper bound on the time
complexity of each step is O((n−d log(1/))ρ ).
We can speed up this algorithm by using the idea of coresets. First observe that our local search
algorithm extends to the weighted
Psetting where each point j ∈ X has a weight w(j) and the cost
of a clustering with centres C is j∈X w(j) · δ 2 (j, C).
Let X ⊆ Rd be a weighted set P
of points with weight function w : X → R. For any set C ⊆ Rd
of size k we use µC (X ) to denote j∈X w(x)δ(j, C)2 . The cost of the optimum k-means solution
for X is µO (X , k) = minC⊆Rd :|C|=k µC (X ).
Definition 4 Let X ⊆ Rd and S ⊆ Rd be two weighted point sets. For any k, , we say S is a
(k, )-coreset if for every set of k centres C ⊆ Rd :
(1 − )µC (X ) ≤ µC (S) ≤ (1 + )µC (X ).
A set Z ⊆ Rd is a (k, )-centroid set for X if there is a set C ⊆ Z of size k such that µC (X ) ≤
(1 + )µO (X , k).
Earlier works on coresets for k-means [30, 25, 26, 41] imply the existence of a (k, )-coresets
of small size. For example, [30] show existence of (k, )-coresets of size O(k 3 /d+1 ). Applying the
result of [46], we get a (k, )-centroid set of size O(log(1/)k 3 /2d+1 ) over a (k, )-coreset of size
O(k 3 /d+1 ). Thus running our local search ρ-swap algorithm takes O((k/)ζ ) time per iteration
where ζ = dO(d) · −O(d/) .

18

4

The Partitioning Scheme: Proof of Theorem 5

Here is the overview of our partitioning scheme. First we bucket the real values by defining bucket
a to have real values r where −a ≤ r < −(a+1) . Fix b = Θ(1/) and then choose a random off-set
a0 and consider b consecutive buckets a0 + ` · b, . . . , a0 + (` + 1) · b − 1 to form “bands”. We get
a random partition of points i ∈ O ∪ S into bands B1 , B2 , . . . such that all centres i in the same
band B` have the property that their Di values are in buckets a0 + ` · b, . . . , a0 + (` + 1) · b − 1
0
0
(so −(a +`·b) ≤ Di < −(a +(`+1)·b−1) ). Given that for each pair (i, i∗ ) ∈ N the Di , Di∗ values are
within a factor 1/ of each other, the probability that they fall into different bands is small (like
/4). Next we impose a random grid partition on each band B` to cut it into cells with cell width
roughly Θ(d/`·b ); again the randomness comes from choosing a random off-set for our grid. The
randomness helps to bound the probability of any pair (i, i∗ ) ∈ N being cut into two different cells
to be small again. This will ensure the 3rd property of the theorem holds. Furthermore, since
0
0
each i ∈ B` has Di ≥ −(a +`·b) and by the filtering we did (Lemma 1), balls of radius 2 · −(a +`·b)
around them must be disjoint; hence using a simple volume/packing argument we can bound the
number of centres from each B` that fall into the same grid cell by a function of  and d. This
helps of establish the first property. We have to do some clean-up to ensure that for each pair
(i, cent(i)) ∈ T they belong to the same part (this will imply property 2 of the Theorem) and that
at the end for each part P , |P ∩ S| = |P ∩ O|. These are a bit technical and are described in details
below.
Start by geometrically grouping the centres in O ∪ S. For a ∈ Z, let


1
1
Ga := i ∈ O ∪ S : a ≤ Di < a+1 .


Finally, let G−∞ := {i ∈ O ∪ S : Di = 0}. Note that each i ∈ O ∪ S appears in exactly one set
among {Ga : a ∈ Z} ∪ {G−∞ }.
We treat G−∞ differently in our partitioning algorithm. It is important to note that no pair in
T or N has precisely one point in G−∞ , as the following shows.
Lemma 6 For each pair of centres (i∗ , i) ∈ T ∪ N , |{i, i∗ } ∩ G−∞ | =
6 1.
Proof. Consider some colocated pair (i∗ , i) ∈ S × O. As Di = Di∗ = 0 and because no other centre
in S ∪ O is colocated with i and i∗ , then i ∈ S and i∗ ∈ O. Thus, φ(i∗ ) = i and it is the unique
closest facility in O to i so i∗ = cent(i). This shows every pair (i∗ , i) ∈ T with either Di = 0 or
Di∗ = 0 must have both Di = Di∗ = 0 so i∗ , i ∈ G−∞ .
Next consider some (i∗ , i) ∈ N . We know  · Di ≤ Di∗ by definition of N . Thus, if i∗ ∈ G−∞
then i ∈ G−∞ as well. Conversely, suppose i ∈ G−∞ . Since (i∗ , i) ∈ N then δ(i, i∗ ) ≤ −1 · Di = 0.
So, Di∗ = 0 meaning i∗ ∈ G−∞ as well.

4.1

Partitioning {Ga : a ∈ Z}

Step 1: Forming Bands
Fix b to be the smallest integer that is at least 4/. We first partition {Ga : a ∈ Z} into bands.
Sample an integer shift a0 uniformly at random in the range {0, 1, . . . , b − 1}. For each ` ∈ Z, form
the band B` as follows:
[
B` :=
Ga0 +j+`·b .
0≤j≤b−1

19

Step 2: Cutting Out Cells
Focus on a band B` . Let W` := 4d · −(`+2)·b−1 (recall that d is the dimension of the Euclidean
metric). This will be the “width” of the cells we create. It is worth mentioning that this is the
first time we are going to use the properties of Euclidean metrics as all our arguments so far were
treating δ(·, ·) as a general metric.
We consider a random grid with cells having width W` in each dimension. Choose an offset
β ∈ Rd uniformly at random from the cube [0, W` ]d . For emphasis, we let pij refer to the j’th
component of i’s point in Euclidean space (this is the only time we will refer specifically to the
coordinates for a point). So centre i has Euclidean coordinates (pi1 , pi2 , . . . , pid ).
For any tuple of integers a ∈ Zd define the cell Ca` as follows.
Ca` = {i ∈ B` : βj + W` · aj ≤ pij < βj + W` · (aj + 1) for all 1 ≤ j ≤ d}.
These are all points in the band B` that lie in the half-open cube with side length W` and lowest
corner β + W` · a. Each i ∈ B` lies in precisely one cell as these half-open cubes tile Rd .
Step 3: Fixing T
Let I ⊆ O be the centres of the form cent(i) for some (cent(i), i) ∈ T where i, cent(i) 6∈ G−∞ and i
and cent(i) lie in different cells after step 2. We will simply move each i ∈ I to the cell containing
φ(i). More precisely, for ` ∈ Z and a ∈ Z we define the part Pa` as
Pa` = (Ca` − I) ∪ {i ∈ I : φ(i) ∈ Ca` }.
We will show that the parts Pa` essentially satisfy most of the desired properties stated about
the random partitioning scheme promised in Theorem 5. It is easy to see that property 2 holds for
each part Pa` since for each pair (cent(i), i) ∈ T both of them belong to the same part. All that
remains is to ensure they are balanced (i.e. have include the same number of centres in O ∪ S)
and to incorporate G−∞ and the centres in O ∪ S that were filtered out. These are relatively easy
cleanup steps.

4.2

Properties of the Partitioning Scheme

We will show the parts formed in the partitioning scheme so far have constant size (depending only
on  and d) and also that each pair in N is cut with low probability. Figure 6 illustrates some key
concepts in these proofs.
Lemma 7 For each ` ∈ Z and a ∈ Zd we have |Pa` | ≤ 2 · (2d)2d · −9d/ .
Proof. Note that each centre i ∈ Ca` witnesses the inclusion of at most one additional centre of I
into Pa` (in particular, cent(i)). So, |Pa` | ≤ 2 · |Ca` | meaning it suffices to show |Ca` | ≤ (2d)2d · −9d/ .
For each i ∈ Ca` we have Di ≥ −`·b . By Lemma 1, the balls of radius 2 · −`·b around the centres
in Ca` must be disjoint subsets of Rd .
 √ d
2π
· Rd , so the total
The volume of a ball of radius R in Rd can be lower bounded by
d
volume of all of these balls of radius


2

· −`·b is at least

|Ca` |

√
d
π·
· √
· −d·`·b .
2·d

On the other hand, these balls have their centres located in a common cube with side length
W` and, thus, all such balls are contained in a slightly larger cube with side length W` + Di ≤
20

Wl

Wl

Di
i
Figure 6: The random grid & balls of radii ·D
2 (solid) and  (dashed) around each i ∈ B` . Note
·Di
the balls of radius 2 are disjoint, which is used to bound |Ca` | by a volume argument in Lemma
7. For each pair (i, i∗ ) ∈ T ∪ N in the band B` , one centre must lie in the ball of radius Di around
the other. The proof of Lemma 8 essentially shows that each such ball crosses a grid line with very
small probability, so such a pair is probably not cut by the random grid.

W` + −(`+2)·b+1 . The total volume of all balls is at most (W` + −(`+2)·b+1 )d ≤ (5d−(2b+1) )d · −d·`·b .
Combining this with the lower bound on the total volume, we see
!d
√
50 · d2 · −(2b+2)
`
√
|Ca | ≤
≤ (2d)2d · −9d/ .
π
The last bound follows for sufficiently small  and because b is the smallest integer at least 4/.
Lemma 8 For each (i∗ , i) ∈ N with i, i∗ 6∈ G−∞ we have Pr[i, i∗ lie in different parts] ≤ .
Proof. We first bound the probability they lie in different bands. Note that Di∗ ≤ δ(i, i∗ ) ≤ Di /
and also Di ≤ Di∗ because (i∗ , i) ∈ N . Thus, if we say i ∈ Ga and i∗ ∈ Ga∗ then |a − a∗ | ≤ 1.
The probability that Ga and Ga∗ are separated when forming the bands B` is at most 1/b ≤ /4.
Next, conditioned on the event that i, i∗ lie in the same band B` , we bound the probability they
lie in different cells. Say i, i∗ ∈ B` and note

· W` .
δ(i, i∗ ) ≤ −1 · Di ≤ −(`+2)·b ≤
4d
The probability that i and i∗ are cut by the random offset along one of the d-dimensions is then at

most 4d
. Taking the union bound over all dimensions, the probability that i and i∗ lie in different
cells is at most 4 .
21

Finally, we bound the probability that i∗ will be moved to a different part when fixing T . If
(i∗ , φ(i∗ )) 6∈ T or if i∗ 6∈ I then this will not happen. So, we will bound Pr[i∗ ∈ I] if (i∗ , φ(i∗ )) ∈ T .
That is, we bound the probability that i0 , i∗ lie in different parts where i0 is such that i∗ = cent(i0 )
and (i∗ , i0 ) ∈ T . Note that Di0 ≤ δ(i0 , i∗ ) ≤ (1 + ) · Di∗ ≤ −1 · Di∗ by Lemma 2 and Di0 ≥
 · δ(i0 , i∗ ) =  · Di∗ by definition of T . So i0 and i∗ lie in different bands with probability at most

∗
0 ∗
4 by the same argument as with (i, i ). Similarly, conditioned on i , i lying in the same band, the
probability they lie in different cells is at most 4 again by the same arguments as with (i, i∗ ).
Note that if i, i∗ lie in different parts then they were cut by the random band or by the random
box, or else i0 , i∗ were cut by the random band or the random box (the latter may not apply if i∗
is not involved in a pair in T ). Thus, by the union bound we have
Pr[i, i∗ are in different parts] ≤ Pr[i, i∗ lie in different B` ] +
Pr[i, i∗ lie in different Ca` | i, i∗ lie in the same B` ] +
Pr[i∗ , i0 lie in different B` ] +
Pr[i∗ , i0 lie in different Ca` | i∗ , i0 lie in the same B` ]




+ + + =
≤
4 4 4 4

4.3

Balancing the Parts

We have partitioned {Ga : a ∈ Z} into parts Pa` for various ` ∈ Z and a ∈ Zd . We extend this to
a partition of all of O ∪ S into constant-size parts that are balanced between O and S to complete
the proof of Theorem 5. Let P = {Pa` : ` ∈ Z, a ∈ Zd and Pa` 6= ∅} be the collection of nonempty
parts formed so far.
The proof of Lemma 6 shows that G−∞ partitions naturally into colocated centres. So let P−∞
denote the partition of G−∞ into these pairs. Finally let P 0 denote the partition of (O −O)∪(S −S)
into singleton sets.
We summarize important properties of P ∪ P−∞ ∪ P 0 .
• P ∪ P−∞ ∪ P 0 is itself a partitioning of O ∪ S into parts with size at most 2 · (2d)2d · −9d/
(Lemma 7).
• Every (cent(i), i) ∈ T pair has both endpoints in the same part. This is true for pairs with
endpoints in P by Step 3 in the partitioning of {Ga : a ∈ Z}. If i ∈ G−∞ this is trivially true
as colocated pairs (i∗ , i) ∈ O ∪ S must have i∗ = cent(i).
• Over the randomized formation of P, each (i∗ , i) ∈ N has endpoints in different parts with
probability at most . For i, i∗ lying in P this follows from Lemma 8. For i, i∗ lying in
P−∞ this follows because they must then be colocated pairs so they always form a part by
themselves.
From now on, we simply let P = P ∪ P−∞ ∪ P 0 . We show how to combine parts of P into
constant-size parts that are also balanced between O and S. Since merging parts does not destroy
the property of two centres lying together, we will still have that pairs of centres in T appear
together and any pair of centres in N lie in different parts with probability at most .
For any subset A ⊆ O ∪ S, let µ(A) = |A ∩ O| − |A ∩ S| denote the imbalance of A.
Lemma
P 9 Let Y ≥ 1 be an integer and A a collection of disjoint, nonempty subsets of O ∪ S such3
that A∈A
Pµ(A) = 0. If |A| ≤ Y for each A then there is some nonempty B ⊆ A where |B| ≤ 2Y
such that B∈B µ(B) = 0.
22

Proof. If µ(A) = 0 for some A ∈ A then simply let B = {A}. Otherwise, note |µ(A)| ≤ |A| ≤ Y for
each A ∈ A and partition A into sets Ax := {A ∈ A : µ(A) = x} for x ∈ {1, . . . , Y }∪{−1, . . . , −Y }.
We have
Y
Y
Y
Y
X
X
X
X
|A−x | ≤
|A−x | · x =
|Ax | · x ≤ Y ·
|Ax |,
x=1

x=1

x=1

x=1

PY

P
so if x=1 |Ax | ≤ Y 2 then |A| ≤ 2Y 3 and we can take B = A. Similarly if Yx=1 |A−x | ≤ Y 2 then
we can take B = A.
P
P
Finally, we are left with the case that Yx=1 |Ax | and Yx=1 |A−x | both exceed Y 2 . By the
pigeonhole principle, there are values 1 ≤ x, y ≤ Y such that |Ax | ≥ Y and |A−y | ≥ Y . In this
case, we take B to be any y sets from Ax plus any x sets from A−y . Then |B| = x + y ≤ 2Y and B
is balanced, which is what we needed to show.
To complete the partitioning we iteratively apply Lemma 9 to P with Y = 2 · (2d)2d · −9d/
where we note |P | ≤ Y by Lemma 7. Each iteration, we find some nonempty Q ⊆ P such that
P
Q∈Q µ(Q) = 0. Remove Q from P and repeat until all sets from P have been removed. Each
3
balanced part obtained is the union of at most 2 · 2 · (2d)2d · −9d/ different parts in P, meaning
4
each part has size at most 2 · 2 · (2d)2d · −9d/ = 32 · (2d)8d · −36·d/ .

Proof of Theorem 1 and Extension to `qq -norm k-clustering

5

In this section we show how to extended the analysis presented in the previous sections for Rd to
prove Theorem 1. For doubling metrics we use ρ = ρ(, d) = 32 · d16d · −256d/ and consider the
ρ-swap local search. First we argue why we get a PTAS when the metric δ(·, ·) is a doubling metric.

5.1

Extending to Doubling Metrics

Let (V, δ) be a doubling metric space with doubling dimension d (so each ball of radius 2r can be
covered with 2d balls of radius r in V ). Note that the only places in the analysis where we used
the properties of metric being Rd was in the proof of Theorem 5 and in particular in Step 2 where
we cut out cells from bands and then later in the proof of Lemma 7. The rest of the proof remains
unchanged.
Given a metric (V, δ), the aspect ratio, denoted by ∆, is the ratio of the largest distance to the
maxu,v∈V δ(u,v)
smallest non-zero distance in the metric: ∆ = minu,v∈V,u6
. Talwar [49] gave a hierarchical
=v δ(u,v)
decomposition of a doubling metric using an algorithm similar to one by Fakcharoenphol et al. [24].
It assumes minu,v∈V,u6=v δ(u, v) = 1, which can be accomplished by scaling the distances. So, ∆ is
the maximum distance between two points in the metric.
Theorem 6 [49] Suppose minu,v∈V :u6=v δ(u, v) = 1. There is a randomized hierarchical decomposition of V , which is a sequence of partitions P0 , P1 , . . . , Ph , where Pi−1 is a refinement of Pi ,
Ph = {V }, and P0 = {{v}}v∈V . The decomposition has the following properties:
1. P0 corresponds to the leaves and Ph corresponds to the root of the split-tree T , and the height
of T is h = ϕ + 2, where ϕ = log ∆ and ∆ is the aspect ratio of metric.
2. For each level i and each S ∈ Pi , S has diameter at most 2i+1 .
3. The branching factor b of T is at most 12d .

23

4. For any u, v ∈ V , the probability that they are in different sets corresponding to nodes in level
i of T is at most 5d · 2−i · δ(u, v).
The proof of randomized partitioning scheme for the case of doubling metrics follows similarly
as in Section 4. We create bands B` as before. For each band B` we cut out the cells (that define
1
the parts) in the following way. First, observe that for all i ∈ B` we have 1`b ≤ Di < (`+2)b−1
.
1
∗
∗
Additionally, by Lemma 1 for all i ∈ S ∩ B` , i ∈ O ∩ B` we have `b−1 ≤ δ(i, i ). So, if we scale all
distances between points in B` so the minimum distance is 1, for all i ∈ B` we have 1 ≤  · Di < 12b
and for any two distinct i, i∗ ∈ B` we have 1 ≤ δ(i, i∗ ).
Consider the hierarchical decomposition of the points in B` using Theorem 6 in this scaled
metric and consider the clusters defined by the sets at level λ = log d + (2b + 3) log 1 . These will
define our cells for B` , let’s call them C1` , . . . , Cq` .1 The only two things that we need to do is: a)
bound the size of each cell Cr` and b) to show that Lemma 8 still holds. Bounding |Cr` | will be
equivalent to Lemma 7. We then we apply Step 3 to fix T as before, this will at most double the
size of each cell. The rest of the proof of Theorem 5 remains unchanged as it does not rely on the
geometry of Rd .
To bound |Ci` | we use property 3 of Theorem 6 and the fact that the bottom level P0 consists
of singleton sets:
|Ci` | ≤ 12d·λ
1

= 2(log 12)d·[(2b+3) log  +log d]
≤ −64d/ · d4d ,
where we used the definition of b in the last inequality. Using similar argument as in the case of
Rd , we can bound the size of each part by 4Y 4 where Y ≤ 2|Ci` |, which implies an upper bound of
32 · d16d · −256d/ .
As for proving equivalent version of Lemma 8, it is enough to show that Pr[i, i∗ lie in different Cr` ]
given that they are in the same band B` is at most /4 since the rest of that proof is the same. For
such pair i, i∗ , note that δ(i, i∗ ) ≤ −1 · Di ≤ −2b−1 . Using property 4 of Theorem 6:
δ(i, i∗ )
2λ
−2b−1
≤ 5d −2b−3
d
2
= 5

Pr[i, i∗ lie in different Cr` | i, i∗ lie in the same B` ] ≤ 5d

≤ /4.
The rest of the proof of Theorem 5 remains the same as in Subsection 4.3.

5.2

Extending to `qq -norm k-clustering

It is fairly straightfoward to see that all the calculations in the proof of Theorem 4 where we bound
∆Pj will hold with similar bounds if the objective function is sum of distances (instead of squares of
1

It is easy to see that we really don’t need the full hierarchical decomposition as in Theorem 6. One can simply
run one round of the Decompose algorithm of [49] that carves out clusters with diameter size at most 2λ and show
that the relevant properties of Theorem 6 hold. For ease of presentation we simply invoke Theorem 6.

24

Algorithm 3 ρ-Swap Local Search for uncapacitated facility location
Let S ← C
while ∃ set ∅ ( S 0 ⊆ C, |S − S 0 |, |S 0 − S| ≤ ρ s.t. costufl (S 0 ) < costufl (S) do
S ← S0
return S
distances), that is `1 -clustering (i.e. k-median). For the case of `qq (with q > 2) it is also possible
to verify that whenever we have a term of the form O()(c∗j + cj ) when upper bounding ∆Pj in the
analysis of k-means (e.g. when j is lucky, or long, or good) then the equivalent term for `qq will be
O(2q )(δjq + δj∗ q ). For the case when j is bad we still get an upper bound of ∆Pj ≤ O(2q (δjq + δj∗ q )).
One can use the crude bound of (x +  · y)q ≤ xq + 2q  · max{x, y} for all x, y ≥ 0 to help bound
terms like (dj +  · Dσ(j) )q . We skip the straightforward (but tedious) details.
This means we get a (1 + O())-approximation for k-clustering when the objective function is
measured according to the `qq -norm for any fixed q ≥ 1 by considering swaps of size ρ = dO(d) ·
q
(2q /)−d·2 / in the local search algorithm.

6

Extending to the uncapacitated facility location and generalized
k-median problems

Recall that the setting of uncapacitated facility location is similar to k-median except we
are given opening costs for each i ∈ C instead of a cardinality bound k. A feasible solution is any
nonempty S ⊆ C and its cost is
X
X
costufl (S) =
δ(j, S) +
fi .
j∈X

i∈S

The standard local search algorithm is the following. Let ρ be defined as before.
We will show every locally optimum solution has cost at most (1 + O()) · OP T using at
0
most
2 · |C|
 test swaps. Thus, replacing the cost condition in Algorithm 1 with costufl (S ) ≤


1 − 2·|C|
· costufl (S) is a polynomial-time variant that is also a PTAS.
Let S be a locally optimum solution and O a global optimum. We may assume, by duplicating
points if necessary (for the analysis only), that S ∩ O = ∅, thus |S| + |O| ≤ 2 · |C| (where |C| refers
to the size of the original set of facilities before duplication). Our analysis proceeds in a manner
that is nearly identical to our approach for k-means.
In particular, we prove the following. It is identical to Theorem 5 in every way except the first
point does not require P ∩ O and P ∩ S to have equal size.
Theorem 7 There is a randomized algorithm that samples a partitioning π of O ∪ S such that:
• For each part P ∈ π, |P ∩ O|, |P ∩ S| ≤ ρ(, d).
• For each part P ∈ π, S4P includes at least one centre from every pair in T .
• For each (i∗ , i) ∈ N , Pr[i, i∗ lie in different parts of π] ≤ .
The proof is identical to Theorem 5, except we do not to the “balancing” step in Section 4.3. Indeed
this was the only step in the proof of Theorem 5 that required |S| = |O|.
We now complete the proof of Theorem 2.
Proof. Let cj = δ(j, S) and c∗j = δ(j, O) for each j ∈ X . Sample a partition π of S ∪ O as per
25

Algorithm 4 ρ-Swap Local Search for uncapacitated facility location
Let S be any k centres from S
while ∃ set S 0 ⊆ C, with |S − S 0 |, |S 0 − S| ≤ ρ and 1 ≤ |S 0 | ≤ k s.t. costgkm (S 0 ) < costgkm (S)
do
S ← S0
return S
Theorem 7. For each part P ∈ π and each j ∈ X we let ∆Pj denote δ(j, S4P ) − δ(j, S). Using the
same bounds considered in the proof of Theorem 4 we have
"
#
X p
Eπ
∆j ≤ (1 + O()) · c∗j − (1 − O())cj .
P ∈π

As each i ∈ S is closed exactly once and each i∗ ∈ O is opened exactly once, the fact that S is
locally optimal and fact that each P ∈ π has |P ∩ S|, |P ∩ O| ≤ ρ shows
"
#
X
0 ≤ Eπ
costufl (S4P ) − costufl (S)
P ∈π


≤ 


X

!

(1 + O()) · c∗j − (1 − O())cj  +

X
i∗ ∈O

j∈X

fi∗ −

X

fi

.

i∈S

Rearranging shows costufl (S) ≤ (1 + O())costufl (O).
Finally we conclude by analyzing a local search algorithm for generalized k-median. Here
we are given both a cardinality bound k and opening costs for each i ∈ C. The goal is to open some
S ⊆ C with 1 ≤ |S| ≤ k to minimize
X
X
costgkm (S) =
δ(j, S) +
fi .
j∈X

i∈S

Note this is the same as costufl , we use this slightly different notation to avoid confusion as we are
discussing a different problem now.
We can prove Theorem 3 very easily.
Proof. Let S and O denote a local optimum and a global optimum (respectively). By adding
artificial points i with fi = 0 that are extremely far away from X , we may assume S = O. Then
we simply use our original partitioning result, Theorem 5, to define test swaps. Noting that each
i ∈ S is closed exactly once and each i∗ ∈ O is opened exactly once over the swaps induced by a
partition π, we proceed in the same way as above in the proof of Theorem 2 and see
costgkm (S) ≤ (1 + O()) · costgkm (O).
Using almost identical arguments as we did in Section 5.2 for extension of k-means to `qq -norm
in doubling dimensions, one can easily extend the results of Theorem 2 and 3 to the case where we
consider the sum of q’th power of distances (`qq -norm).

26

7

Conclusion

We have presented a PTAS for k-means in constant-dimensional Euclidean metrics and, more
generally, in doubling metrics and when the objective is to minimizing the `qq -norm of distances
between points and their nearest centres for any constant q ≥ 1. This is also the first approximation
for k-median in doubling metrics and the first demonstration that local search yields a true PTAS
for k-median even in the Euclidean plane. Our approach extends to showing local search yields a
PTAS for uncapacitated facility location and generalized k-median in doubling metrics.
The running time of a single step of the local search algorithm is O(k ρ ) where ρ = dO(d) ·−O(d/)
for the case of k-means when the metric has doubling dimension d. We have not tried to optimize
the constants in the O(·) notations in ρ. The dependence on d cannot be improved much under
the Exponential Time Hypothesis (ETH). For example, if the running time of a single iteration
1−δ
was only O(k exp(d )·f () ) for some constant δ then we would have a sub-exponential time (1 + )approximation when d = O(log n). Recall that k-means is APX-hard when d = Θ(log n) [10], so
this would refute the ETH.
It may still be possible to obtain an EPTAS for any constant dimension d. That is, there could
be a PTAS with running time of the form O(g()·nexp(d) ) for some function g() (perhaps depending
also on d). Finally, what is the fastest PTAS we can obtain in the special case of the Euclidean
plane (i.e. d = 2)? It would be interesting to see if there is an EPTAS whose running time is linear
or near linear in n for any fixed constant .

References
[1] Daniel Aloise, Amit Deshpande, Pierre Hansen, and Preyas Popat. NP-hardness of Euclidean
sum-of-squares clustering. Mach. Learn., 75(2):245–248, 2009.
[2] Sanjeev Arora. Polynomial time approximation schemes for Euclidean traveling salesman and
other geometric problems. J. ACM, 45(5):753–782, 1998.
[3] 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
Theory of Computing (STOC ’98), pages 106–113. ACM, 1998.
[4] David Arthur, Bodo Manthey, and Heiko Röglin. Smoothed analysis of the K-means method.
J. ACM, 58(5):19:1–19:31, 2011.
[5] David Arthur and Sergei Vassilvitskii. How slow is the K-means method? In Proceedings of the
Twenty-second Annual Symposium on Computational Geometry (SoCG ’06), pages 144–153.
ACM, 2006.
[6] 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
’07), pages 1027–1035. SIAM, 2007.
[7] Vijay Arya, Naveen Garg, Rohit Khandekar, Adam Meyerson, Kamesh Munagala, and
Vinayaka Pandit. Local search heuristic for K-median and facility location problems. In
Proceedings of the Thirty-third Annual ACM Symposium on Theory of Computing (STOC
’01), pages 21–29. ACM, 2001.

27

[8] 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.
[9] Pranjal Awasthi, Avrim Blum, and Or Sheffet. Stability yields a PTAS for K-median and
K-means clustering. In Proceedings of the 2010 IEEE 51st Annual Symposium on Foundations
of Computer Science (FOCS ’10), pages 309–318. IEEE Computer Society, 2010.
[10] Pranjal Awasthi, Moses Charikar, Ravishankar Krishnaswamy, and Ali Kemal Sinop. The
Hardness of Approximation of Euclidean K-Means. In Proceedings of 31st International Symposium on Computational Geometry (SoCG ’15), Leibniz International Proceedings in Informatics (LIPIcs), pages 754–767, 2015.
[11] Mihai Bādoiu, Sariel Har-Peled, and Piotr Indyk. Approximate clustering via Coresets. In
Proceedings of the Thiry-fourth Annual ACM Symposium on Theory of Computing (STOC
’02), pages 250–257. ACM, 2002.
[12] Sayan Bandyapadhyay and Kasturi Varadarajan. On variants of K-means clustering. In
Proceedings of the 32nd International Symposium on Computational Geometry (SoCG ’16),
Leibniz International Proceedings in Informatics (LIPIcs), 2016.
[13] Johannes Blömer, Christiane Lammersen, Melanie Schmidt, and Christian Sohler. Theoretical
analysis of the K-means algorithm - A survey. CoRR, abs/1602.08254, 2016.
[14] 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 ’15), pages 737–756. SIAM, 2015.
[15] Ke Chen. On Coresets for K-median and K-means clustering in metric and Euclidean spaces
and their applications. SIAM J. Comput., 39:923–947, 2009.
[16] 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.
[17] Vincent Cohen-Addad, Philip N. Klein, and Claire Mathieu. The power of local search for
clustering. CoRR, abs/1603.09535, 2016.
[18] Vincent Cohen-Addad and Claire Mathieu. Effectiveness of local search for geometric optimization. In Proceedings of the 31st International Symposium on Computational Geometry
(SoCG ’15), Leibniz International Proceedings in Informatics (LIPIcs), 2015.
[19] Vincent Cohen-Addad and Claire Mathieu. Effectiveness of local search for geometric optimization. In Proceedings of 31st Annual Symposium on Computational Geometry (SOCG ’15),
pages 329–344, 2015.
[20] Sanjoy Dasgupta. How fast is K-means? In Proceedings of the 16th Annual Conference on
Learning Theory and 7th Kernel Workshop (COLT/Kernel ’03), pages 735–735, 2003.

28

[21] W. Fernandez de la Vega, Marek Karpinski, Claire Kenyon, and Yuval Rabani. Approximation
schemes for clustering problems. In Proceedings of the Thirty-fifth Annual ACM Symposium
on Theory of Computing (STOC ’03), pages 50–58. ACM, 2003.
[22] Patrik D’haeseleer et al. How does gene expression clustering work? Nature biotechnology,
23(12):1499–1502, 2005.
[23] P. Drineas, A. Frieze, R. Kannan, S. Vempala, and V. Vinay. Clustering large graphs via the
singular value decomposition. Mach. Learn., 56(1-3):9–33, 2004.
[24] Jittat Fakcharoenphol, Satish Rao, and Kunal Talwar. A tight bound on approximating arbitrary metrics by tree metrics. In Proceedings of the Thirty-fifth Annual ACM Symposium on
Theory of Computing (STOC ’03), pages 448–455. ACM, 2003.
[25] Dan Feldman, Morteza Monemizadeh, and Christian Sohler. A PTAS for K-means clustering
based on weak Coresets. In Proceedings of the Twenty-third Annual Symposium on Computational Geometry (SoCG ’07), SoCG ’07, pages 11–18. ACM, 2007.
[26] Dan Feldman, Melanie Schmidt, and Christian Sohler. Turning big data into tiny data:
Constant-size Coresets for K-means, PCA and projective clustering. In Proceedings of the
Twenty-Fourth Annual ACM-SIAM Symposium on Discrete Algorithms (SODA ’13), pages
1434–1453. SIAM, 2013.
[27] 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.
[28] Zachary Friggstad, Mohsen Rezapour, and Mohammad R. Salavatipour. Local search yields a
PTAS for k-means in doubling metrics. CoRR, abs/1603.08976, 2016.
[29] A. Gupta and T. Tangwongsan. Simpler analyses of local search algorithms for facility location.
CoRR, abs/0809.2554, 2008.
[30] Sariel Har-Peled and Akash Kushal. Smaller Coresets for K-median and K-means clustering. In
Proceedings of the Twenty-first Annual Symposium on Computational Geometry (SoCG ’05),
pages 126–134. ACM, 2005.
[31] Sariel Har-Peled and Soham Mazumdar. On Coresets for K-means and K-median clustering.
In Proceedings of the Thirty-sixth Annual ACM Symposium on Theory of Computing (STOC
’04), pages 291–300. ACM, 2004.
[32] Sariel Har-Peled and Bardia Sadri. How fast is the K-means method? Algorithmica, 41(3):185–
202, 2005.
[33] Matan Hofree, John P Shen, Hannah Carter, Andrew Gross, and Trey Ideker. Network-based
stratification of tumor mutations. Nature methods, 10(11):1108–1115, 2013.
[34] Mary Inaba, Naoki Katoh, and Hiroshi Imai. Applications of weighted Voronoi diagrams and
randomization to variance-based K-clustering. In Proceedings of the tenth annual Symposium
on Computational Geometry (SoCG ’94), pages 332–339. ACM, 1994.

29

[35] Anil K. Jain. Data clustering: 50 years beyond K-means. Pattern Recogn. Lett., 31(8):651–666,
2010.
[36] 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. Theory Appl., 28(2-3):89–112, 2004.
[37] Stavros G. Kolliopoulos and Satish Rao. A nearly linear-time approximation scheme for the
Euclidean Kappa-median problem. In Proceedings of the 7th Annual European Symposium on
Algorithms (ESA ’99), pages 378–389. Springer-Verlag, 1999.
[38] Amit Kumar and Ravindran Kannan. Clustering with spectral norm and the K-means algorithm. In Proceedings of the 2010 IEEE 51st Annual Symposium on Foundations of Computer
Science (FOCS ’10), pages 299–308. IEEE Computer Society, 2010.
[39] Amit Kumar, Yogish Sabharwal, and Sandeep Sen. A simple linear time (1+ ) -approximation
algorithm for K-means clustering in any dimensions. In Proceedings of the 45th Annual IEEE
Symposium on Foundations of Computer Science (FOCS ’04), pages 454–462. IEEE Computer
Society, 2004.
[40] Amit Kumar, Yogish Sabharwal, and Sandeep Sen. Linear-time approximation schemes for
clustering problems in any dimensions. J. ACM, 57(2):5:1–5:32, 2010.
[41] Michael Langberg and Leonard J. Schulman. Universal -approximators for integrals. In Proceedings of the Twenty-first Annual ACM-SIAM Symposium on Discrete Algorithms (SODA
’10), pages 598–607. SIAM, 2010.
[42] Shi Li. A 1.488-approximation for the uncapacitated facility location problem. In Proceedings of
the 38th Annual International Colloquium on Automata, Languages and Programming (ICALP
’11), pages 45–58, 2011.
[43] Shi Li and Ola Svensson. Approximating K-median via pseudo-approximation. In Proceedings
of the Forty-fifth Annual ACM Symposium on Theory of Computing (STOC ’13), pages 901–
910. ACM, 2013.
[44] S. Lloyd. Least squares quantization in pcm. IEEE Trans. Inf. Theor., 28(2):129–137, 2006.
[45] Meena Mahajan, Prajakta Nimbhorkar, and Kasturi Varadarajan. The planar K-means problem is NP-hard. In Proceedings of the 3rd International Workshop on Algorithms and Computation (WALCOM ’09), pages 274–285. Springer-Verlag, 2009.
[46] Jirı Matoušek. On approximate geometric k-clustering. Discrete & Computational Geometry,
24(1):61–84, 2000.
[47] Rafail Ostrovsky and Yuval Rabani. Polynomial-Time Approximation Schemes for geometric
min-sum median clustering. J. ACM, 49(2):139–156, 2002.
[48] 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:1–28:22, 2013.
[49] Kunal Talwar. Bypassing the embedding: Algorithms for low dimensional metrics. In Proceedings of the Thirty-sixth Annual ACM Symposium on Theory of Computing (STOC ’04),
pages 281–290. ACM, 2004.
30

[50] Andrea Vattani. The hardness of K-means clustering in the plane. Manuscript, 2009.
[51] Andrea Vattani. K-means requires exponentially many iterations even in the plane. Discrete
Comput. Geom., 45(4):596–616, 2011.

31

