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

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 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.

This value is called the cost of the clustering. Typically, the centres $c_{i}$ are selected to be the centroid (mean) of the cluster $C_{i}$. In other situations the centres must be from the data points themselves (i.e. $c_{i}\in C_{i}$) or from a given set ${\mathcal{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$ .

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 : 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 ).

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 . 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 ). 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 ), but none of them are known to have a bounded approximation ratio in the general setting. Arthur and Vassilvitskii show that Lloyd’s method with properly chosen initial centres will be an $O(\log k)$-approximation. Ostrovsky et al. 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 , Dastupta and Gupta , and Har-Peled and Sadri study convergence rate of Lloyd’s algorithm. In particular show that it can be super-polynomial. More recently, Vattani shows that it can take exponential time even in two dimensions. Arthur et al. proved that Lloyd’s algorithm has polynomial-time smoothed complexity. Kumar and Kannan , Ostrovsky et al. , and Awasthi et al. gave empirical and theoretical evidence for when and why the known heuristics work well in practice. For instance 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 . In fact, the $k$-means problem is NP-hard if $d$ is arbitrary even for $k=2$ . Also, if $k$ is arbitrary the problem is NP-hard even for $d=2$ . However, the $k$-means problem can be solved in polynomial time by the algorithm of when both $k$ and $d$ are constant.

$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 the sum of distances, instead of sum of square of distances as in $k$-means, i.e. minimize $\sum_{i=1}^{k}\sum_{x\in C_{i}}\delta(x,c_{i})$ where $\delta(x,c_{i})$ is the distance between $x$ and $c_{i}$. 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+\epsilon$ approximation by Arya et al. . The current best approximation uses different techniques and has an approximation ratio of $2.611+\epsilon$ . The local search $(9+\epsilon)$-approximation (for $k$-means) in can be seen as an extension of the analysis in 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. , building on the framework of Arora , gave the first PTAS. Kolliopoulos and Rao improved the time complexity and presented a PTAS for Euclidean $k$-median with time complexity $O(2^{O((1+\log\frac{1}{\epsilon})/\epsilon)^{d-1}}n\log n\log k)$. Such an approximation is also known as an efficient PTAS: the running time of the $(1+{\epsilon})$-approximation is of the form $f({\epsilon})\cdot{\rm 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 given opening costs $f_{i}$ for each $i\in{\mathcal{C}}$. The goal is to find a set of centers ${\mathcal{S}}\subseteq{\mathcal{C}}$ that minimizes $\sum_{x\in{\mathcal{X}}}\delta(x,{\mathcal{S}})+\sum_{i\in{\mathcal{S}}}f_{i}$. Currently the best approximation for uncapacitated facility location in general metrics is a 1.488-approximation . As with $k$-median, a PTAS is known for uncapacitated facility location in constant-dimensional Euclidean metrics , with the latter giving an efficient PTAS.

A precise description of the algorithm is given in Section 2. At a high level, we start with any set of $k$ centres ${\mathcal{S}}\subseteq{\mathcal{C}}$. Then, while there is some other set of $k$ centres ${\mathcal{S}}^{\prime}\subseteq{\mathcal{C}}$ with $|{\mathcal{S}}-{\mathcal{S}}^{\prime}|\leq\rho$ for some constant $\rho$ such that ${\mathcal{S}}^{\prime}$ is a cheaper solution than ${\mathcal{S}}$, we set ${\mathcal{S}}\leftarrow{\mathcal{S}}^{\prime}$. Repeat until ${\mathcal{S}}$ cannot be improved any further. Each iteration takes $|{\mathcal{C}}|^{O(\rho)}$ time, which is polynomial when $\rho$ is a constant. Such a solution is called a local optimum solution with respect to the $\rho$-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.

Let $\rho({\epsilon},d):=d^{O(d)}\cdot{\epsilon}^{O(-d/{\epsilon})}$. We will articulate the absolute constants suppressed by the $O(\cdot)$ notation later on in our analysis.

The local search algorithm that swaps up to $\rho({\epsilon},d)$ centres at a time is a $(1+\epsilon)$-approximation for $k$-means in metrics with doubling dimension $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, shows that local search can be used to get a $1+{\epsilon}$ approximation for $k$-median that uses up to $(1+{\epsilon})\cdot k$ centres.

As mentioned earlier, Awasthi et al. proved that $k$-means is APX-hard for $d=\Omega(\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+\epsilon)$-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<\epsilon<1$ and $d=\sigma\log n/\log\log n$ for sufficiently small absolute constant $\sigma$ we get a $(1+\epsilon)$-approximation for $k$-means that runs in time $O(2^{n^{\kappa}})$, for some constant $\kappa=\kappa(\sigma)<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 of $d=\Omega(\log n)$ to prove APX-hardness of $k$-means is almost tight unless ${\rm NP}\subseteq DTIME(2^{n^{\sigma}})$.

The notion of coresets and using them for finding faster algorithms for $k$-means has been studied extensively (e.g. 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\in{\mathcal{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)\cdot\delta(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 $\rho$ as in Theorem 1 (in fact, it can be slightly smaller) we consider the natural local search algorithm that returns a solution ${\mathcal{S}}\subseteq{\mathcal{C}}$ such that ${\rm cost}({\mathcal{S}})\leq{\rm cost}({\mathcal{S}}^{\prime})$ for all ${\mathcal{S}}^{\prime}\subseteq{\mathcal{C}}$ with both $|{\mathcal{S}}-{\mathcal{S}}^{\prime}|\leq\rho$ and $|{\mathcal{S}}^{\prime}-{\mathcal{S}}|\leq\rho$ (i.e. add and/or drop up to $\rho$ centres in ${\mathcal{C}}$).

The local search algorithm that adds and/or drops up to $\rho({\epsilon},d)$ centres at a time is a $(1+\epsilon)$-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 . 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 . 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 $\rho$ 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 .

The local search algorithm that adds and/or drops up to $\rho({\epsilon},d)$ centres at a time (provided the resulting solution has at most $k$ facilities) is a a $(1+\epsilon)$-approximation for generalized $k$-median in metrics with doubling dimension $d$.

Note: Shortly after we announced our result , Cohen-Addad, Klein, and Mathieu 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.

2 Proof Outline

The general framework for analysis of local search algorithms for $k$-median and $k$-means in is as follows. Let ${\mathcal{S}}$ and ${\mathcal{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, if ${\mathcal{B}}({\mathcal{S}})$ is the set of solutions obtained by performing swaps from $Q$, the main task is to show that $\sum_{{\mathcal{S}}^{\prime}\in{\mathcal{B}}({\mathcal{S}})}({\rm cost}({\mathcal{S}}^{\prime})-{\rm cost}({\mathcal{S}}))\leq\alpha\cdot{\rm cost}({\mathcal{O}})-{\rm cost}({\mathcal{S}})$ for some constant $\alpha$. Given that $0\leq{\rm cost}({\mathcal{S}}^{\prime})-{\rm cost}({\mathcal{S}})$ for all ${\mathcal{S}}^{\prime}\in{\mathcal{B}}({\mathcal{S}})$ (because ${\mathcal{S}}$ is a local optimum), ${\rm cost}({\mathcal{S}})\leq\alpha\cdot{\rm cost}({\mathcal{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.

Let us define ${\mathcal{S}}$ and ${\mathcal{O}}$ as before. First, we do a filtering over ${\mathcal{S}}$ and ${\mathcal{O}}$ to obtain subsets $\overline{{\mathcal{S}}}\subseteq{\mathcal{S}}$ and $\overline{{\mathcal{O}}}\subseteq{\mathcal{O}}$ such that every centre in ${\mathcal{S}}-\overline{{\mathcal{S}}}$ (in ${\mathcal{O}}-\overline{{\mathcal{O}}}$) is “close” to a centre in $\overline{{\mathcal{S}}}$ (in $\overline{{\mathcal{O}}}$) while these filtered centres are far apart. We define a “net” around each centre $i\in\overline{{\mathcal{S}}}$ which captures a collection of other filtered centres in $\overline{{\mathcal{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 ${\mathcal{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}-c_{j}$. 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({\epsilon})\cdot(c^{*}_{j}+c_{j})$ and we will generate the $c^{*}_{j}-c_{j}$ 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 ${\mathcal{S}}$ than it is to $i$.

Notation and Preliminaries

We usually refer to a potential centre in ${\mathcal{C}}$ by a simple index $i$ and a point in ${\mathcal{X}}$ by a simple index $j$ (or slight variants like $i^{*}$ or $\overline{i^{\prime}}$). 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 5.1. So, for any set $S\subseteq{\mathcal{C}}$ and any $j\in{\mathcal{X}}$, let $\delta(j,S)=\min_{i\in S}\delta(j,i)$. We also define ${\rm cost}(S)=\sum_{j\in{\mathcal{X}}}\delta(j,S)^{2}$.

Our goal in (discrete) $k$-means is to find a set of centres $S\subseteq{\mathcal{C}}$ of size $k$ to minimize ${\rm cost}(S)$. Note that once we fix the set of centres we can find a partitioning of ${\mathcal{X}}$ that realizes $\sum_{j\in{\mathcal{X}}}\delta(j,S)^{2}$ by assigning each $j\in{\mathcal{X}}$ to the nearest centre in $S$, breaking ties arbitrarily.

The simple $\rho$-swap local search heuristic shown in Algorithm 1 is essentially the same one considered in .

Recall that we defined $\rho({\epsilon},d)=d^{O(d)}\cdot{\epsilon}^{O(d/{\epsilon})}$, where the constants will be specified later and consider the local search algorithm with $\rho=\rho({\epsilon},d)$ swaps. By a standard argument (as in ) one can show that replacing the condition of the while loop with ${\rm cost}(({\mathcal{S}}-Q)\cup P)\leq(1-\frac{{\epsilon}}{k})\cdot{\rm cost}({\mathcal{S}})$, the algorithm terminates in polynomial time. Furthermore, if $\alpha$ is such that any locally optimum solution returned by Algorithm 1 has cost at most $\alpha\cdot{\rm cost}({\mathcal{O}})$ where ${\mathcal{O}}$ denotes a global optimum solution, then any ${\mathcal{S}}$ such that ${\rm cost}(({\mathcal{S}}-Q)\cup P)<(1-\frac{{\epsilon}}{k})\cdot{\rm cost}({\mathcal{S}})$ for any possible swap $P,Q$ satisfies ${\rm cost}({\mathcal{S}})\leq\frac{\alpha}{1-{\epsilon}}{\rm cost}({\mathcal{O}})$. This follows by arguments in and the fact that our local search analysis uses at most $k$ “test swaps”.

For ease of exposition, we ignore this factor $1+\epsilon$ loss, and consider the solution ${\mathcal{S}}$ returned by Algorithm 1. Recall that we use ${\mathcal{O}}$ to denote the global optimum solution. For $j\in{\mathcal{X}}$, let $c^{*}_{j}=\delta(j,{\mathcal{O}})^{2}$ and $c_{j}=\delta(j,{\mathcal{S}})^{2}$, so ${\rm cost}({\mathcal{O}})=\sum_{j\in{\mathcal{X}}}c^{*}_{j}$ and ${\rm cost}({\mathcal{S}})=\sum_{j\in{\mathcal{X}}}c_{j}$. We also denote the centre in ${\mathcal{O}}$ nearest to $j$ by $\sigma^{*}(j)$ and the centre in ${\mathcal{S}}$ nearest to $j$ by $\sigma(j)$. Define $\phi:{\mathcal{O}}\cup{\mathcal{S}}\rightarrow{\mathcal{O}}\cup{\mathcal{S}}$ to be the function that assigns $i^{*}\in{\mathcal{O}}$ to its nearest centre in ${\mathcal{S}}$ and assigns $i\in{\mathcal{S}}$ to its nearest centre in ${\mathcal{O}}$. For any two sets $S,T\subseteq{\mathcal{O}}\cup{\mathcal{S}}$, we let $S\triangle T=(S\cup T)-(S\cap T)$.

We assume ${\mathcal{O}}\cap{\mathcal{S}}=\emptyset$. This is without loss of generality because we could duplicate each location in ${\mathcal{C}}$ and say ${\mathcal{O}}$ uses the originals and ${\mathcal{S}}$ the duplicates. It is easy to check that ${\mathcal{S}}$ would still be a locally optimum solution in this instance. We can also assume that these are the only possible colocated facilities, so $\delta(i,i^{\prime})>0$ for distinct $i,i^{\prime}\in{\mathcal{O}}$ or distinct $i,i^{\prime}\in{\mathcal{S}}$. Finally, we will assume ${\epsilon}$ is sufficiently small (independent of all other parameters, including $d$) so that all of our bounds hold.

To prove this, we will construct a set of test swaps that yield various inequalities which, when combined, provide the desired bound on ${\rm cost}({\mathcal{S}})$. That is, we will partition ${\mathcal{O}}\cup{\mathcal{S}}$ into sets where $|P\cap{\mathcal{O}}|=|P\cap{\mathcal{S}}|\leq\rho({\epsilon},d)$ for each part $P$. For each such set $P$, $0\leq{\rm cost}({\mathcal{S}}\triangle P)-{\rm cost}({\mathcal{S}})$ because ${\mathcal{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 ${\rm cost}({\mathcal{S}})\leq(1+O({\epsilon}))\cdot{\rm cost}({\mathcal{O}})$. For example, for a point $j\in{\mathcal{X}}$ if $\sigma^{*}(j)\in P$ then the change in $j$’s assignment cost is at most $c^{*}_{j}-c_{j}$ because we could assign $j$ from $\sigma(j)$ to $\sigma^{*}(j)$. The problem is that points $j$ with $\sigma(j)\in P$ but $\sigma^{*}(j)\not\in 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 ${\mathcal{O}}\cup{\mathcal{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.

For $i^{*}\in{\mathcal{O}}$ let $D_{i^{*}}:=\delta(i^{*},{\mathcal{S}})=\delta(i^{*},\phi(i^{*}))$. For $i\in{\mathcal{S}}$ let $D_{i}:=\delta(i,{\mathcal{O}})=\delta(i,\phi(i))$.

The first thing is to sparsify ${\mathcal{O}}$ and ${\mathcal{S}}$ using a simple filtering step. Algorithm 2 filters ${\mathcal{O}}$ to a set that is appropriately sparse for our analysis.

Think of $\eta(i^{*})$ as a proxy for $i^{*}$ that is very close to $i^{*}$. Using a similar process, we filter ${\mathcal{S}}$ to get $\overline{{\mathcal{S}}}$ and proxy centres $\eta(i)\in\overline{{\mathcal{S}}}$ for each $i\in{\mathcal{S}}$. The idea is that the set of centres left in $\overline{{\mathcal{O}}}$ and $\overline{{\mathcal{S}}}$ are somewhat far apart yet any point that was assigned to a centre in ${\mathcal{O}}-\overline{{\mathcal{O}}}$ (or in ${\mathcal{S}}-\overline{{\mathcal{S}}}$) can be “cheaply” reassigned to a proxy.

For each $i\in{\mathcal{O}}\cup{\mathcal{S}}$ we have $\delta(i,\eta(i))\leq{\epsilon}\cdot D_{i}$. For any distinct $i,i^{\prime}\in\overline{{\mathcal{O}}}\cup\overline{{\mathcal{S}}}$ we have $\delta(i,i^{\prime})\geq{\epsilon}\cdot\max\{D_{i},D_{i^{\prime}}\}$.

Proof. That $\delta(i,\eta(i))\leq{\epsilon}\cdot D_{i}$ follows immediately by construction. If $i\in\overline{{\mathcal{O}}},i^{\prime}\in\overline{{\mathcal{S}}}$ or vice-versa, then in fact $\delta(i,i^{\prime})\geq\max\{D_{i},D_{i^{\prime}}\}$ simply by definition of $D_{i},D_{i^{\prime}}$.

Now suppose $i,i^{\prime}\in\overline{{\mathcal{O}}}$ and that $i^{\prime}$ was considered after $i$ in Algorithm 2 (so $D_{i^{\prime}}\geq D_{i}$). The fact that $i^{\prime}$ was added to $\overline{{\mathcal{O}}}$ even though $i$ was already in $\overline{{\mathcal{O}}}$ means $\delta(i,i^{\prime})\geq{\epsilon}\cdot D_{i^{\prime}}$. The same argument works if $i,i^{\prime}\in\overline{{\mathcal{S}}}$.

Next we define mappings similar to $\phi,\sigma,\sigma^{*}$ except they only concern centres that were not filtered out.

${\overline{\phi}}:\overline{{\mathcal{O}}}\cup\overline{{\mathcal{S}}}\rightarrow\overline{{\mathcal{O}}}\cup\overline{{\mathcal{S}}}$ maps each $i\in\overline{{\mathcal{O}}}$ to its nearest location in $\overline{{\mathcal{S}}}$ and vice versa.

${\overline{\sigma}}^{*}:{\mathcal{X}}\rightarrow\overline{{\mathcal{O}}}$ defined by ${\overline{\sigma}}^{*}(j)=\eta(\sigma^{*}(j))$.

${\overline{\sigma}}:{\mathcal{X}}\rightarrow\overline{{\mathcal{S}}}$ defined by ${\overline{\sigma}}(j)=\eta(\sigma(j))$.

Finally, for each $i\in{\overline{\phi}}(\overline{{\mathcal{O}}})$, let ${\rm cent}(i)$ be the centre in ${\overline{\phi}}^{-1}(i)$ that is closest to $i$, breaking ties arbitrarily.

Note ${\overline{\sigma}}(j)$ may not necessarily be the centre in $\overline{{\mathcal{S}}}$ that is closest to $j$. Also note that if one considers a bipartite graph with parts $\overline{{\mathcal{O}}}$ and $\overline{{\mathcal{S}}}$, then ${\overline{\phi}}$ maps centres from one side to the other.

For each $i^{\prime}\in\overline{{\mathcal{O}}}\cup\overline{{\mathcal{S}}}$, $D_{i^{\prime}}\leq\delta(i^{\prime},{\overline{\phi}}(i^{\prime}))\leq(1+{\epsilon})\cdot D_{i^{\prime}}$.

Proof. Suppose $i^{\prime}\in\overline{{\mathcal{O}}}$, the proof is essentially the same for $i^{\prime}\in\overline{{\mathcal{S}}}$. On one hand, we know

because $\overline{{\mathcal{S}}}\subseteq{\mathcal{S}}$. On the other hand,

Conclude by observing $D_{\phi(i^{\prime})}\leq\delta(i^{\prime},\phi(i^{\prime}))=D_{i^{\prime}}$.

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.

${\mathcal{T}}:=\{({\rm cent}(i),i):i\in{\overline{\phi}}(\overline{{\mathcal{O}}}){\rm~{}and~{}}{\epsilon}\cdot\delta({\rm cent}(i),i)\leq D_{i}\}$

${\mathcal{N}}:=\{(i^{*},i)\in\overline{{\mathcal{O}}}\times\overline{{\mathcal{S}}}:\delta(i,i^{*})\leq{\epsilon}^{-1}\cdot D_{i}{\rm~{}and~{}}D_{i^{*}}\geq{\epsilon}\cdot D_{i}\}$

For each $i\in\overline{{\mathcal{S}}}$, the set $\{i^{*}:(i^{*},i)\in{\mathcal{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 ${\mathcal{T}}$ are not separated by the partition. Our requirement for ${\mathcal{N}}$ is not quite as strong. The partition is constructed randomly and it will be sufficient to have each pair in ${\mathcal{N}}$ being separated by the partition with probability at most ${\epsilon}$.

The following says that if at least one centre of each pair in ${\mathcal{T}}$ is open after a swap, then every centre in $\overline{{\mathcal{O}}}\cup\overline{{\mathcal{S}}}$ is somewhat close to some open centre. The bound is a bit big, but it will be multiplied by $O({\epsilon})$ whenever it is used in the local search analysis.

Let $A\subseteq\overline{{\mathcal{O}}}\cup\overline{{\mathcal{S}}}$ be such that $A\cap\{{\rm cent}(i),i\}\neq\emptyset$ for each $({\rm cent}(i),i)\in{\mathcal{T}}$. Then $\delta(i^{\prime},A)\leq 5\cdot D_{i^{\prime}}$ for any $i^{\prime}\in{\mathcal{O}}\cup{\mathcal{S}}$.

Proof. We first prove the statement for $i^{\prime}\in{\mathcal{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 $i_{0}:=i^{\prime}$, $i_{1}:=\eta(i_{0})$, and $i_{2}:={\overline{\phi}}(i_{1})$. We build the rest inductively, noting that we guarantee $i_{a}\in{\overline{\phi}}(\overline{{\mathcal{O}}})$ for even indices $a\geq 2$ (so ${\rm cent}(i_{a})$ is defined).

Inductively, for even $a\geq 2$ we do the following. If $i_{a}\in A$ then we stop. Otherwise, if $({\rm cent}(i_{a}),i_{a})\in{\mathcal{T}}$ then by assumption it must be that ${\rm cent}(i_{a})\in A$ so we let $i_{a+1}:={\rm cent}(i_{a})$ and stop. Finally, if $({\rm cent}(i_{a}),i_{a})\not\in{\mathcal{T}}$ then we set $i_{a+1}:={\overline{\phi}}(i_{a})$ and $i_{a+2}:={\overline{\phi}}(i_{a+1})$ and iterate with $a^{\prime}=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\geq 2$ such that the walk did not stop at $i_{a}$. If $i_{a+1}={\overline{\phi}}(i_{a})$ then

If $i_{a+1}\neq{\overline{\phi}}(i_{a})$ then it must be $i_{a}={\rm cent}(i_{a-1})$. In this case, it must be $i_{a}={\overline{\phi}}(i_{a-1})$, so because ${\rm cent}(i_{a})$ is the closest centre in ${\overline{\phi}}^{-1}(i_{a})$ to $i_{a}$ (by definition), we have

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.

For every even $a\geq 2$ such that the walk did not end at $i_{a}$ or $i_{a+1}$, we have $\delta(i_{a},i_{a+1})\leq 2{\epsilon}\cdot\delta(i_{a-1},i_{a})$.

Proof. Because the walk did not end at $i_{a}$ or $i_{a+1}$, $({\rm cent}(i_{a}),i_{a})\not\in{\mathcal{T}}$ meaning $D_{i_{a}}<{\epsilon}\cdot\delta({\rm cent}(i_{a}),i_{a})\leq{\epsilon}\cdot\delta(i_{a-1},i_{a})$. Using this and Lemma 2 in the first bound below, we see

Let $m$ be the index of the last centre in the walk. Thus, $\delta(i_{a},i_{a+1})\leq\delta(i_{a-1},i_{a})$ for all $2\leq a\leq m-1$. From this and Claim 1 we have

The second last step uses Lemma 2 and the last step uses the fact that $D_{\eta(i^{\prime})}\leq D_{i^{\prime}}$ (either $i^{\prime}=\eta(i^{\prime})$ or else $i^{\prime}$ was filtered out by $\eta(i^{\prime})$, in which case it has a larger $D$-value) and the assumption that ${\epsilon}$ is small enough.

Now suppose $i^{\prime}\in{\mathcal{S}}$. We bound $\delta(i^{\prime},A)$ mostly using what we have done already. That is, we have

Note $D_{\phi(i^{\prime})}\leq\delta(i^{\prime},\phi(i^{\prime}))=D_{i^{\prime}}$ again by Lemma 2. So,

2 Good Partitioning of 𝒪∪𝒮𝒪𝒮{\mathcal{O}}\cup{\mathcal{S}} and Proof of Theorem 4

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

There is a randomized algorithm that samples a partitioning $\pi$ of ${\mathcal{O}}\cup{\mathcal{S}}$ such that:

For each part $P\in\pi$, $|P\cap{\mathcal{O}}|=|P\cap{\mathcal{S}}|\leq\rho({\epsilon},d)$.

For each part $P\in\pi$, ${\mathcal{S}}\triangle P$ includes at least one centre from every pair in ${\mathcal{T}}$.

For each $(i^{*},i)\in{\mathcal{N}}$, $\Pr[i,i^{*}{\rm~{}lie~{}in~{}different~{}parts~{}of~{}}\pi]\leq{\epsilon}$.

The following gives a way to handle the fact that the triangle inequality does not hold with squares of the distances.

For any real numbers $x,y$ we have $(x+y)^{2}\leq 2(x^{2}+y^{2})$.

Proof. $(x+y)^{2}\leq(x+y)^{2}+(x-y)^{2}=2x^{2}+2y^{2}.$

For each point $j\in{\mathcal{X}}$, $D_{{\overline{\sigma}}(j)}\leq D_{\sigma(j)}\leq\delta(j,\sigma(j))+\delta(j,\sigma^{*}(j))$. Similarly, $D_{{\overline{\sigma}}^{*}(j)}\leq D_{\sigma^{*}(j)}\leq\delta(j,\sigma(j))+\delta(j,\sigma^{*}(j))$.

Proof. As usual, we only prove the first statement since the second is nearly identical. If ${\overline{\sigma}}(j)=\sigma(j)$ then $D_{{\overline{\sigma}}(j)}=D_{\sigma(j)}$ is trivially true. Otherwise, ${\overline{\sigma}}(j)=\eta(\sigma(j))$ was already in $\overline{{\mathcal{S}}}$ when $\sigma(j)$ was considered by the filtering algorithm meaning $D_{{\overline{\sigma}}(j)}\leq D_{\sigma(j)}$.

Proof of Theorem 4. Let $\pi$ be a partition sampled by the algorithm from Theorem 5. For each point $j\in{\mathcal{X}}$ and each part $P$ of $\pi$, let $\Delta^{P}_{j}:=\delta(j,{\mathcal{S}}\triangle P)^{2}-\delta(j,{\mathcal{S}})^{2}$ denote the change in assignment cost for the point after swapping in the centers in $P\cap{\mathcal{O}}$ and swapping out $P\cap{\mathcal{S}}$. Local optimality of ${\mathcal{S}}$ and $|P\cap{\mathcal{S}}|=|P\cap{\mathcal{O}}|\leq\rho({\epsilon},d)$ means $0\leq\sum_{j}\Delta^{P}_{j}$ for any part $P$.

Classify each point $j\in{\mathcal{X}}$ in one of the following ways:

Lucky: $\sigma(j)$ and ${\overline{\sigma}}(j)$ do not lie in the same part of $\pi$.

Long: $j$ is not lucky but $\delta({\overline{\sigma}}(j),{\overline{\sigma}}^{*}(j))>{\epsilon}^{-1}\cdot D_{{\overline{\sigma}}(j)}$.

Bad: $j$ is not lucky or long and $({\overline{\sigma}}^{*}(j),{\overline{\sigma}}(j))\in{\mathcal{N}}$ yet ${\overline{\sigma}}(j),{\overline{\sigma}}^{*}(j)$ lie in different parts of $\pi$.

Good: $j$ is neither lucky, long, nor bad.

We now place an upper bound on $\sum_{P\in\pi}\Delta^{P}_{j}$ for each point $j\in{\mathcal{X}}$. Note that each centre in ${\mathcal{S}}$ is swapped out exactly once over all swaps $P$ and each centre in ${\mathcal{O}}$ is swapped in exactly once. With this in mind, consider the following cases for a point $j\in{\mathcal{X}}$. In the coming arguments, we let $\delta_{j}:=\delta(j,\sigma(j))$ and $\delta^{*}_{j}:=\delta(j,\sigma^{*}(j))$ for brevity. Note $c_{j}=\delta_{j}^{2}$ and $c^{*}_{j}=\delta^{*2}_{j}$.

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 ${\mathcal{S}}\triangle P$ by first moving it to either $\sigma(j)$ or $\sigma^{*}(j)$ and then moving it a distance of $O({\epsilon})\cdot(\delta_{j}+\delta^{*}_{j})$ to reach an open facility. Considering that we reassigned $j$ from $\sigma(j)$, the reassignment cost will be

Case: $j$ is lucky For the part $P\in\pi$ with $\sigma^{*}(j)\in P$, we have $\Delta^{P}_{j}\leq c^{*}_{j}-c_{j}$ as we could move $j$ from $\sigma(j)$ to $\sigma^{*}(j)$. If $\sigma(j)$ is swapped out in a different swap $P^{\prime}$, we move $j$ to ${\overline{\sigma}}(j)$ (which remains open because $j$ is lucky) and bound $\Delta^{P^{\prime}}_{j}$ by:

again using the assumption that ${\epsilon}$ is sufficiently small. For every other swap $P^{\prime\prime}$, we have that $\sigma(j)$ remains open after the swap so $\Delta^{P^{\prime\prime}}_{j}\leq 0$ as we could just leave $j$ at $\sigma(j$). In total, we have

Case: $j$ is long Again, for $P\in\pi$ with $\sigma^{*}(j)\in P$ we get $\Delta^{P}_{j}\leq c^{*}_{j}-c_{j}$. If $\sigma(j)$ is swapped out in a different swap $P^{\prime}$, then we bound $\Delta^{P^{\prime}}_{j}$ by moving $j$ from $\sigma(j)$ to the open centre nearest to $\sigma(j)$. Note that ${\mathcal{S}}\triangle P^{\prime}$ contains at least one centre from every pair in ${\mathcal{T}}$, so we bound this distance using Lemma 3. This case is depicted in Figure 3. We have

Using this, we bound $\Delta^{P^{\prime}}_{j}$ as follows.

In every other swap we could leave $j$ at $\sigma(j)$. Thus, for a long point $j$ we have

Case: $j$ is bad We only move $j$ in the swap $P$ when $\sigma(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 $c_{j}$. This will eventually be compensated by the fact that $j$ is bad with probability at most ${\epsilon}$ over the random choice of $\pi$. For now, we just provide the reassignment cost bound for bad $j$.

Case: $j$ is good This breaks into two subcases. We know $\delta({\overline{\sigma}}^{*}(j),{\overline{\sigma}}(j))\leq{\epsilon}^{-1}\cdot D_{{\overline{\sigma}}(j)}$ because $j$ is not long. In one subcase, $D_{{\overline{\sigma}}^{*}(j)}\geq{\epsilon}\cdot D_{{\overline{\sigma}}(j)}$ so $({\overline{\sigma}}^{*}(j),{\overline{\sigma}}(j))\in{\mathcal{N}}$. Since $j$ is not bad and not lucky, we have $\sigma(j),{\overline{\sigma}}(j),{\overline{\sigma}}^{*}(j)\in P$ for some common part $P\in\pi$. In the other subcase, $D_{{\overline{\sigma}}^{*}(j)}<{\epsilon}\cdot D_{{\overline{\sigma}}(j)}$. Note in this case we still have $\sigma(j),{\overline{\sigma}}(j)\in P$ for some common part $P$ because $j$ is not lucky.

Subcase: $D_{{\overline{\sigma}}^{*}(j)}\geq{\epsilon}\cdot D_{{\overline{\sigma}}(j)}$ The only time we move $j$ is when $\sigma(j)$ is closed. As observed in the previous paragraph, this happens in the same swap when ${\overline{\sigma}}^{*}(j)$ is opened, so send $j$ to ${\overline{\sigma}}^{*}(j)$. This is illustrated in Figure 4.

Subcase: $D_{{\overline{\sigma}}^{*}(j)}<{\epsilon}\cdot D_{{\overline{\sigma}}(j)}$ Again, the only time we move $j$ is when $\sigma(j)$ is closed. We reassign $j$ by first moving it to $\sigma^{*}(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 $\delta(\sigma^{*}(j),{\overline{\sigma}}^{*}(j))\leq{\epsilon}D_{\sigma^{*}(j)}$ by our filtering. This implies that $D_{\sigma^{*}(j)}\leq D_{{\overline{\sigma}}^{*}(j)}+{\epsilon}D_{\sigma^{*}(j)}$ which in turn implies $D_{\sigma^{*}(j)}\leq\frac{1}{1-{\epsilon}}D_{{\overline{\sigma}}^{*}(j)}$; thus $D_{\sigma^{*}(j)}\leq(1+2{\epsilon})D_{{\overline{\sigma}}^{*}(j)}$ which in turn is bounded by ${\epsilon}(1+{\epsilon})D_{{\overline{\sigma}}(j)}$ in this subcase (by the assumption of subcase). Thus,

Considering both subcases we can say that for any good point $j$ that

Aggregating these bounds and remembering $0\leq\sum_{j\in{\mathcal{X}}}\Delta^{P}_{j}$ for each $P\in\pi$ because ${\mathcal{S}}$ is a locally optimum solution, we have

The last step is to average this inequality over the random choice of $\pi$. Note that any point $j\in{\mathcal{X}}$ is bad with probability at most ${\epsilon}$ by the guarantee in Theorem 5 and the definition of bad. Thus, we see

3 Running Time Analysis

Recall that we can go to the discrete case by finding a set ${\mathcal{C}}$ of size $O(n\epsilon^{-d}\log(1/\epsilon)$. The analysis of Arya et al. shows that the number of local search steps is at most $\log(cost(S_{0})/cost({\mathcal{O}}))/\log\frac{1}{1-{\epsilon}/k}$, where $S_{0}$ 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 $\rho=\rho({\epsilon},d)$, a crude upper bound on the time complexity of each step is $O((n{\epsilon}^{-d}\log(1/{\epsilon}))^{\rho})$.

We can speed up this algorithm by using the idea of coresets. First observe that our local search algorithm extends to the weighted setting where each point $j\in{\mathcal{X}}$ has a weight $w(j)$ and the cost of a clustering with centres $C$ is $\sum_{j\in{\mathcal{X}}}w(j)\cdot\delta^{2}(j,C)$.

Earlier works on coresets for $k$-means imply the existence of a $(k,{\epsilon})$-coresets of small size. For example, show existence of $(k,{\epsilon})$-coresets of size $O(k^{3}/{\epsilon}^{d+1})$. Applying the result of , we get a $(k,{\epsilon})$-centroid set of size $O(\log(1/{\epsilon})k^{3}/{\epsilon}^{2d+1})$ over a $(k,{\epsilon})$-coreset of size $O(k^{3}/{\epsilon}^{d+1})$. Thus running our local search $\rho$-swap algorithm takes $O((k/{\epsilon})^{\zeta})$ time per iteration where $\zeta=d^{O(d)}\cdot{\epsilon}^{-O(d/{\epsilon})}$.

The Partitioning Scheme: Proof of Theorem 5

We treat $G_{-\infty}$ differently in our partitioning algorithm. It is important to note that no pair in ${\mathcal{T}}$ or ${\mathcal{N}}$ has precisely one point in $G_{-\infty}$, as the following shows.

For each pair of centres $(i^{*},i)\in{\mathcal{T}}\cup{\mathcal{N}}$, $|\{i,i^{*}\}\cap G_{-\infty}|\neq 1$.

Proof. Consider some colocated pair $(i^{*},i)\in{\mathcal{S}}\times{\mathcal{O}}$. As $D_{i}=D_{i^{*}}=0$ and because no other centre in ${\mathcal{S}}\cup{\mathcal{O}}$ is colocated with $i$ and $i^{*}$, then $i\in\overline{{\mathcal{S}}}$ and $i^{*}\in\overline{{\mathcal{O}}}$. Thus, ${\overline{\phi}}(i^{*})=i$ and it is the unique closest facility in $\overline{{\mathcal{O}}}$ to $i$ so $i^{*}={\rm cent}(i)$. This shows every pair $(i^{*},i)\in{\mathcal{T}}$ with either $D_{i}=0$ or $D_{i^{*}}=0$ must have both $D_{i}=D_{i^{*}}=0$ so $i^{*},i\in G_{-\infty}$.

Next consider some $(i^{*},i)\in{\mathcal{N}}$. We know ${\epsilon}\cdot D_{i}\leq D_{i^{*}}$ by definition of ${\mathcal{N}}$. Thus, if $i^{*}\in G_{-\infty}$ then $i\in G_{-\infty}$ as well. Conversely, suppose $i\in G_{-\infty}$. Since $(i^{*},i)\in{\mathcal{N}}$ then $\delta(i,i^{*})\leq{\epsilon}^{-1}\cdot D_{i}=0$. So, $D_{i^{*}}=0$ meaning $i^{*}\in G_{-\infty}$ as well.

2 Properties of the Partitioning Scheme

We will show the parts formed in the partitioning scheme so far have constant size (depending only on ${\epsilon}$ and $d$) and also that each pair in ${\mathcal{N}}$ is cut with low probability. Figure 6 illustrates some key concepts in these proofs.

The last bound follows for sufficiently small $\epsilon$ and because $b$ is the smallest integer at least $4/{\epsilon}$.

For each $(i^{*},i)\in{\mathcal{N}}$ with $i,i^{*}\not\in G_{-\infty}$ we have $\Pr[i,i^{*}{\rm~{}lie~{}in~{}different~{}parts}]\leq{\epsilon}$.

The probability that $i$ and $i^{*}$ are cut by the random offset along one of the $d$-dimensions is then at most $\frac{{\epsilon}}{4d}$. Taking the union bound over all dimensions, the probability that $i$ and $i^{*}$ lie in different cells is at most $\frac{{\epsilon}}{4}$.

Finally, we bound the probability that $i^{*}$ will be moved to a different part when fixing ${\mathcal{T}}$. If $(i^{*},{\overline{\phi}}(i^{*}))\not\in{\mathcal{T}}$ or if $i^{*}\not\in\mathcal{I}$ then this will not happen. So, we will bound $\Pr[i^{*}\in\mathcal{I}]$ if $(i^{*},{\overline{\phi}}(i^{*}))\in{\mathcal{T}}$. That is, we bound the probability that $i^{\prime},i^{*}$ lie in different parts where $i^{\prime}$ is such that $i^{*}={\rm cent}(i^{\prime})$ and $(i^{*},i^{\prime})\in{\mathcal{T}}$. Note that $D_{i^{\prime}}\leq\delta(i^{\prime},i^{*})\leq(1+{\epsilon})\cdot D_{i^{*}}\leq{\epsilon}^{-1}\cdot D_{i^{*}}$ by Lemma 2 and $D_{i^{\prime}}\geq{\epsilon}\cdot\delta(i^{\prime},i^{*})={\epsilon}\cdot D_{i^{*}}$ by definition of ${\mathcal{T}}$. So $i^{\prime}$ and $i^{*}$ lie in different bands with probability at most $\frac{{\epsilon}}{4}$ by the same argument as with $(i,i^{*})$. Similarly, conditioned on $i^{\prime},i^{*}$ lying in the same band, the probability they lie in different cells is at most $\frac{{\epsilon}}{4}$ again by the same arguments as with $(i,i^{*})$.