Better Guarantees for k-Means and Euclidean k-Median by Primal-Dual Algorithms

Introduction

A commonly used heuristic for $k$-means is Lloyd’s algorithm , which is based on iterative improvement. However, despite its ubiquity in practice, Lloyd’s algorithm has, in general, no worst-case guarantee and may not even converge in polynomial time . To overcome some of these limitations, Arthur and Vassilvitskii proposed a randomized initialization procedure for Lloyd’s algorithm, called $k$-means$++$, that leads to a $\Theta(\log k)$ expected approximation guarantee in the worst case. Under additional assumptions about the clusterability of the input dataset, Ostrovsky et al. showed that this adaptation of Lloyd’s algorithm gives a PTAS for $k$-means clustering. However, under no such assumptions, the best approximation algorithm in the general case has for some time remained a $(9+\epsilon)$-approximation algorithm based on local search, due to Kanungo et al. . Their analysis also shows that no natural local search algorithm performing a fixed number of swaps can improve upon this ratio. This leads to a barrier for these techniques that are rather far away from the best-known inapproximability result which only says that it is NP-hard to approximate this problem to within a factor better than $1.0013$ .

In summary, while $k$-means is perhaps the most widely used clustering problem in computer science, the only constant-factor approximation algorithm for the general case is based on simple local search heuristics that, for inherent reasons, give guarantees that are rather far from known hardness results. This is in contrast to many other well-studied clustering problems, such as facility location and $k$-median. Over the past several decades, a toolbox of core algorithmic techniques such as dual fitting, primal-dual and LP-rounding, has been refined and applied to these problems leading to improved approximation guarantees . In particular, the current best approximation guarantees for both facility location (a 1.488-approximation due to Li ) and $k$-median (a 2.675-approximation due to Byrka et al. ) are LP-based and give significantly better results than previous local search algorithms . However, such LP-based techniques have not yet been able to attain similar improvements for $k$-means. One reason for this is that they have relied heavily on the triangle inequality, which does not hold in the case of $k$-means.

For any $\epsilon>0$, there is a $(\rho_{\textsf{\tiny mean}}+\epsilon)$-approximation algorithm for the $k$-means problem, where $\rho_{\textsf{\tiny mean}}\approx 6.357$. Moreover, the integrality gap of the standard LP is at most $\rho_{\textsf{\tiny mean}}$.

Using Lagrangian relaxation, we can then consider the resulting discrete problem using the standard linear programming formulation for facility location. This general approach was pioneered in this context by Jain and Vazirani who gave primal-dual algorithms for the $k$-median problem. In their paper, they first present a Lagrangian Multiplier Preserving ($\mathsf{LMP}$) $3$-approximation algorithm for the facility location problem. Then they run binary search over the opening cost of the facilities and use the aforementioned algorithm to get two solutions: one that opens more than $k$ facilities and one that opens fewer than $k$, such that the opening cost of facilities in these solutions are close to each other. These solutions are then combined to obtain a solution that opens exactly $k$ facilities. This step results in losing another factor $2$ in the approximation guarantee, which results in a $6$-approximation algorithm for $k$-median. The factor $6$ was later improved by Jain, Mahdian, and Saberi who obtained a $4$-approximation algorithm for $k$-median by developing an $\mathsf{LMP}$ $2$-approximation algorithm for facility location.

One can see that the same approach gives a much larger constant factor for the $k$-means problem since one cannot anymore rely on the triangle inequality. We use two main ideas to overcome this obstacle: (1) we exploit the geometric structure of $k$-means to obtain an improved $\mathsf{LMP}$-approximation, and (2) we develop a new primal-dual algorithm that opens exactly $k$ facilities while losing only an arbitrarily small factor.

For our first contribution, we modify the primal-dual algorithm of Jain and Vazirani into a parameterized version which allows us to regulate the “aggressiveness” of the opening strategy of facilities. By using properties of Euclidean metrics we show that this leads to improved $\mathsf{LMP}$ approximation algorithms for $k$-means.

By the virtue of , these results already imply upper bounds on the integrality gaps of the standard LP relaxations, albeit with an exponential time rounding algorithm. Our second and more technical contribution is a new primal-dual algorithm that accomplishes the same task in polynomial time. In other words, we are able to turn an $\mathsf{LMP}$ approximation algorithm into an algorithm that opens at most $k$ facilities without deteriorating the approximation guarantee. We believe that this contribution is of independent interest. Indeed, all recent progress on the approximation of $k$-median beyond long-standing local search approaches has involved reducing the factor $2$ that is lost by Jain and Vazirani when two solutions are combined to open exactly $k$ facilities (i.e. in the rounding of a so-called bipoint solution) . Here, we show that it is possible to reduce this loss all the way to $(1+\epsilon)$ by fundamentally changing the way in which dual solutions are constructed and maintained.

Instead of finding two solutions by binary search as in the framework of , we consider a sequence of solutions such that the opening costs and also the dual values of any two consecutive solutions are close in $L^{\infty}$-norm. We show that this latter property allows us to combine two appropriate, consecutive solutions in the sequence into a single solution that opens exactly $k$ facilities while losing only a factor of $1+\epsilon$ (rather than 2) in the approximation guarantee. Unfortunately, the dual solutions produced by the standard primal-dual algorithm approach are unstable, in the sense that a small change in opening price may result in drastic changes in the value of the dual variables. Thus, we introduce a new primal-dual procedure which instead iteratively transforms a dual solution for one price into a dual solution for another price. By carefully constraining the way in which the dual variables are altered, we show that we can obtain a sequence of “close” solutions that can be combined as desired.

We believe that this technique may be applicable in other settings, as well. An especially interesting open question is whether it is possible combine stronger $\mathsf{LMP}$ approximation algorithms, such as the one by Jain, Mahdian, Saberi , with our lossless rounding to obtain an improved $(2+\epsilon)$-approximation algorithm for $k$-median.

For any $\epsilon>0$, there is a $(\rho_{\textsf{\tiny med}}+\epsilon)$-approximation algorithm for the Euclidean $k$-median problem, where $\rho_{\textsf{\tiny med}}\approx 2.633$. Moreover, the integrality gap of the standard LP is at most $\rho_{\textsf{\tiny med}}$.

In the second extension, we consider a variant of the $k$-means problem in which each $c(j,S)$ corresponds to the squared distance in an arbitrary (possibly non-Euclidean) metric on $\cD\cup\cF$. For this problem, the best-known approximation algorithm is a 16-approximation due to Gupta and Tangwongsan . In this paper, we obtain the following improvement:

For any $\epsilon>0$, there is a $(9+\epsilon)$-approximation algorithm for the $k$-means problem in general metrics. Moreover, the integrality gap of the standard LP is at most $9$.

We remark that the same hardness reduction as used for $k$-median immediately yields a much stronger hardness result for the above generalization than what is known for the standard $k$-means problem: it is hard to approximate the $k$-means problem in general metrics within a factor $1+8/e-\epsilon\approx 3.94$ for any $\epsilon>0$.

In Section 2 we review the standard LP formulation that we use, as well as its Lagrangian relaxation. We then in Section 3 show how to exploit the geometric structure of $k$-means and Euclidean $k$-median to give improved $\mathsf{LMP}$ guarantees. In Section 4 we show the main ideas behind our new rounding approach by giving an algorithm that runs in quasi-polynomial time. These results are then generalized in Sections 5, 6, and 7 to obtain an algorithm that runs in polynomial time.

The standard LP relaxation and its Lagrangian relaxation

Here and in the remainder of the paper, we shall consider the discrete $k$-means problem, where we are given a discrete set $\cF$ of facilities (corresponding to candidate centers).As discussed in the introduction, it is well-known that a $\rho$-approximation algorithm for this case can be turned into a $(\rho+\epsilon)$-approximation algorithm for the standard $k$-means problem, for any constant $\epsilon>0$ (see e.g., ). Henceforth, we will simply refer to the discrete $k$-means problem as the $k$-means problem.

Given an instance $(\cD,\cF,d,k)$ of the $k$-means problem or the $k$-median problem, let $c(j,i)$ denote the connection cost of client $j$ if connected to facility $i$. That is, $c(j,i)=d(j,i)$ in the case of $k$-median and $c(j,i)=d(j,i)^{2}$ in the case of $k$-means. Let $n=|\cD|$ and $m=|\cF|$.

The standard linear programming (LP) relaxation of these problems has two sets of variables: a variable $y_{i}$ for each facility $i\in\cF$ and a variable $x_{ij}$ for each facility-client pair $i\in\cF,j\in\cD$. The intuition of these variables is that $y_{i}$ should indicate whether facility $i$ is opened and $x_{ij}$ should indicate whether client $j$ is connected to facility $i$. The standard LP relaxation can now be formulated as follows.

A main difficulty for approximating the $k$-median and the $k$-means problems is the hard constraint that at most $k$ facilities can be selected, i.e., constraint (2.3) in the above relaxation. A popular way of overcoming this difficulty, pioneered in this context by Jain and Vazirani , is to consider the Lagrangian relaxation where we multiply the constraint (2.3) times a Lagrange multiplier $\lambda$ and move it to the objective. This results, for every $\lambda\geq 0$, in the following relaxation and its dual that we denote by LP$(\lambda)$ and DUAL$(\lambda)$, respectively.

LP$(\lambda)$ $\displaystyle\min$ $\displaystyle\sum_{i\in\mathcal{F},j\in\mathcal{D}}x_{ij}\cdot c(j,i)+\lambda\cdot\left(\sum_{i\in\mathcal{F}}y_{i}-k\right)$ s.t. (2.1), (2.2), and (2.4). DUAL$(\lambda)$ $\displaystyle\max$ $\displaystyle\sum_{j\in\mathcal{D}}\alpha_{j}-\lambda\cdot k$ s.t. $\displaystyle\sum_{j\in\mathcal{D}}[\alpha_{j}-c(j,i)]^{+}$ $\displaystyle\leq\lambda\qquad\forall i\in\mathcal{F}$ (2.5) $\displaystyle\alpha$ $\displaystyle\geq 0.$

If we disregard the constant term $\lambda\cdot k$ in the objective functions, LP$(\lambda)$ and DUAL$(\lambda)$ become the standard LP formulation and its dual for the facility location problem where the opening cost of each facility equals $\lambda$ and the connection costs are defined by $c(\cdot,\cdot)$. Recall that the facility location problem (with uniform opening costs) is defined as the problem of selecting a set $S\subseteq\cF$ of facilities to open so as to minimize the opening cost $|S|\lambda$ plus the connection cost $\sum_{j\in\mathcal{D}}c(j,S)$. Jain and Vazirani introduced the following method for addressing the $k$-median problem motivated by simple economics. On the one hand, if $\lambda$ is selected to be very small, i.e., it is cheap to open facilities, then a good algorithm for the facility location problem will open many facilities. On the other hand, if $\lambda$ is selected to be very large, then a good algorithm for the facility location problem will open few facilities. Ideally, we want to use this intuition to find an opening price that leads to the opening of exactly $k$ facilities and thus a solution to the original, constrained problem.

To make this intuition work, we need the notion of Lagrangian Multiplier Preserving ($\mathsf{LMP}$) approximations: we say that a $\rho$-approximation algorithm is $\mathsf{LMP}$ for the facility location problem with opening costs $\lambda$ if it returns a solution $S\subseteq\mathcal{F}$ satisfying

Exploiting Euclidean metrics via primal-dual algorithms

In this section we show how to exploit the structure of Euclidean metrics to achieve better approximation guarantees. Our $\mathsf{LMP}$ approximation algorithm builds upon the primal-dual algorithm for the facility location problem by Jain and Vazirani . We refer to their algorithm as the JV algorithm. The main modification to their algorithm is that we allow for a more “aggressive” opening strategy of facilities. The amount of aggressiveness is measured by the parameter $\delta$: we devise an algorithm JV$(\delta)$ for each parameter $\delta\geq 0$, where a smaller $\delta$ results in a more aggressive opening strategy. We first describe JV$(\delta)$ and we then optimize $\delta$ for the considered objectives to obtain the claimed approximation guarantees.

We remark that the result in (non-constructively) upper bounds the integrality gap of the standard LP relaxation of $k$-median in terms of the $\mathsf{LMP}$ approximation guarantee of JV. This readily generalizes to the $k$-means problem and JV$(\delta)$. Consequently, our guarantees presented here upper bound the integrality gaps as the theorems state in the introduction.

As alluded to above, the algorithm is a modification of JV, and Remark 3.2 below highlights the difference. The algorithm consists of two phases: the dual-growth phase and the pruning phase.

In this stage, we construct a feasible dual solution $\alpha$ to DUAL$(\lambda)$. Initially, we set $\alpha={\bm{0}}$ and let $A=\cD$ denote the set of active clients (which is all clients at first). We then repeat the following until there are no active clients, i.e., $A=\emptyset$: increase the dual-variables $\{\alpha_{j}\}_{j\in A}$ corresponding to the active clients at a uniform rate until one of the following events occur (if several events happen at the same time, break ties arbitrarily):

A dual constraint $\sum_{j\in\cD}[\alpha_{j}-c(j,i)]^{+}\leq\lambda$ becomes tight for a facility $i\in\cF$. In this case we say that facility $i$ is tight or temporarily opened. We update $A$ by removing the active clients with a tight edge to $i$, that is, a client $j\in A$ is removed if $\alpha_{j}-c(j,i)\geq 0$. For future reference, we say that facility $i$ is the witness of these removed clients.

An active client $j\in A$ gets a tight edge, i.e., $\alpha_{j}-c(j,i)=0$, to some already tight facility $i$. In this case, we remove $j$ from $A$ and let $i$ be its witness.

This completes the description of the dual-growth phase. Before proceeding to the pruning phase, let us remark that the constructed $\alpha$ is indeed a feasible solution to DUAL$(\lambda)$ by design. It is clear that $\alpha$ is non-negative. Now consider a facility $i\in\mathcal{F}$ and its corresponding dual constraint $\sum_{j\in\mathcal{D}}[\alpha_{j}-c(j,i)]^{+}\leq\lambda$. On the one hand, the constraint is clearly satisfied if it never becomes tight during the dual-growth phase. On other hand, if it becomes tight, then all clients with a tight edge to it are removed from the active set of clients by Event 1. Moreover, if any client gets a tight edge to $i$ in subsequent iterations it gets immediately removed from the set of active clients by Event 2. Therefore the left-hand-side of the constraint will never increase (nor decrease) after it becomes tight so the constraint remains satisfied. Having proved that $\alpha$ is a feasible solution to DUAL$(\lambda)$, let us now describe the pruning phase.

After the dual-growth phase (too) many facilities are temporarily opened. The pruning phase will select a subset of these facilities to open. In order to formally describe this process, we need the following notation. For a client $j$, let $N(j)=\{i\in\mathcal{F}:\alpha_{j}-c(j,i)>0\}$ denote the facilities to which client $j$ contributes to the opening cost. Similarly, for $i\in\mathcal{F}$, let $N(i)=\{j\in\mathcal{D}:\alpha_{j}-c(j,i)>0\}$ denote the clients with a positive contribution toward $i$’s opening cost. For a temporarily opened facility $i$, let

and by convention let $t_{i}=0$ if $N(i)=\emptyset$ (this convention will be useful in future sections and will only be used when the opening cost $\lambda$ of facilities are set to ). Note that, if $N(i)\neq\emptyset$, then $t_{i}$ equals the “time” that facility $i$ was temporarily opened in the dual-growth phase. A crucial property of $t_{i}$ that follows from the construction of $\alpha$ is the following.

For a client $j$ and its witness $i$, $\alpha_{j}\geq t_{i}$. Moreover, for any $j^{\prime}\in N(i)$ we have $t_{i}\geq\alpha_{j^{\prime}}$.

For the pruning phase, it will be convenient to define the client-facility graph $G$ and the conflict graph $H$. The vertex set of $G$ consist of all the clients and all the facilities $i$ such that $\sum_{j\in\cD}[\alpha_{j}-c(j,i)]^{+}=\lambda$ (i.e., the tight or temporarily open facilities). There is an edge between facility $i$ and client $j$ if $i\in N(j)$. The conflict graph $H$ is defined based on the client-facility graph $G$ and $t$ as follows:

The vertex set consists of all facilities in $G$.

There is an edge between two facilities $i$ and $i^{\prime}$ if some client $j$ is adjacent to both of them in $G$ and $c(i,i^{\prime})\leq\delta\min(t_{i},t_{i^{\prime}})$.

The pruning phase now finds a (inclusion-wise) maximal independent set $\mathsf{IS}$ of $H$ and opens those facilities; clients are connected to the closest facility in $\mathsf{IS}$.

The difference between the original algorithm JV and our modified JV$(\delta)$ is the additional condition $c(i,i^{\prime})\leq\delta\min(t_{i},t_{i^{\prime}})$ in the definition of the conflict graph. Notice that if we select a smaller $\delta$ we will have fewer edges in $H$. Therefore a maximal independent set will likely grow in size, which results in a more “aggressive” opening strategy. Adjusting $\delta$ will allow us to achieve better $\mathsf{LMP}$ approximation guarantees.

2 Analysis of JV(δ)𝛿(\delta) for the considered objectives

In the following subsections, we optimize $\delta$ and analyze the guarantees obtained by the algorithm JV$(\delta)$ for the objective functions: k-means objective in general metrics, standard $k$-means objective (in Euclidean metrics), and k-median objective in Euclidean metrics. The first analysis is very similar to the original JV analysis and also serves as a motivation for the possible improvements in Euclidean metrics.

We consider the case when $c(j,i)=d(j,i)^{2}$ and $d$ forms a general metric. We let $\delta=\infty$ so JV$(\delta)$ becomes simply the JV algorithm. We prove the following.

Let $d$ be any metric on $\cD\cup\cF$ and suppose that $c(j,i)=d(j,i)^{2}$ for every $i\in\cF$ and $j\in\cD$. Then, for any $\lambda\geq 0$, Algorithm JV$(\infty)$ constructs a solution $\alpha$ to DUAL$(\lambda)$ and returns a set $\mathsf{IS}$ of opened facilities such that

Consider any client $j\in\cD$. We shall prove that

The statement then follows by summing up over all clients and noting that any facility $i\in\mathsf{IS}$ was temporarily opened and thus we have $\sum_{j\in\mathcal{D}}[\alpha_{j}-c(j,i)]^{+}=\lambda$.

To prove (3.1), we first note that $|\mathsf{IS}\cap N(j)|\leq 1$. Indeed, consider $i\neq i^{\prime}\in N(j)$. Then $(j,i)$ and $(j,i^{\prime})$ are edges in the client-facility graph $G$ and as $\delta=\infty$, $i$ and $i^{\prime}$ are adjacent in the conflict graph $H$. Hence, the temporarily opened facilities in $N(j)$ form a clique in $H$ and at most one of them can be selected in the maximal independent set $\mathsf{IS}$. We complete the analysis by considering the two cases $|\mathsf{IS}\cap N(j)|=1$ and $|\mathsf{IS}\cap N(j)|=0$.

Let $i^{*}$ be the unique facility in $\mathsf{IS}\cap N(j)$. Then

Notice the amount of slack in the above analysis (specifically, the first inequality). In the Euclidean case, we exploit this slack for a more aggressive opening and to improve the approximation guarantee.

Let $i_{1}$ be $j$’s witness. First, if $i_{1}\in\mathsf{IS}$ then by the same arguments as above we have the desired inequality; specifically, since $j$ has a tight edge to $i_{1}$ but $i_{1}\not\in N(j)$ we must have $\alpha_{j}=c(j,i_{1})$. Now consider the more interesting case when $i_{1}\not\in\mathsf{IS}$. As $\mathsf{IS}$ is a maximal independent set in $H$, there must be a facility $i_{2}\in\mathsf{IS}$ that is adjacent to $i_{1}$ in $H$. By definition of $H$, there is a client $j_{1}$ such that $(j_{1},i_{1})$ and $(j_{1},i_{2})$ are edges in the client-facility graph $G$, i.e., $j_{1}\in N(i_{1})\cap N(i_{2})$. By the definition of witness and $N(\cdot)$, we have

and by the description of the algorithm (see Claim 3.1 in Section 3) we have $\alpha_{j}\geq t_{i_{1}}\geq\alpha_{j_{1}}$. Hence, using the triangle inequality and that $(a+b+c)^{2}\leq 3(a^{2}+b^{2}+c^{2})$,

As $\sum_{i\in N(j)\cap\mathsf{IS}}(\alpha_{j}-c(j,i))=0$, this completes the proof of this case and thus the theorem.

2.2 k𝑘k-Means objective in Euclidean metrics

We start with some intuition that illustrates our approach. From the standard analysis of JV (and our analysis of $k$-means in general metrics), it is clear that the bottleneck for the approximation guarantee comes from the connection-cost analysis of clients that need to do a “$3$-hop” as illustrated in the left part of Figure 1: client $j$ is connected to open facility $i_{2}$ and the squared-distance is bounded by the path $j-i_{1}-j_{1}-i_{2}$. Moreover, this analysis is tight when considering \textsf{\small JV}=\textsf{\small JV(\infty)}. Our strategy will now be as follows: Select $\delta$ to be a constant smaller than $4$. This means that in the configurations of Figure 1, we will also open $i_{2}$ if the distance between $i_{1}$ and $i_{2}$ is close to $2$. Therefore, if we do not open $i_{2}$, the distance between $i_{1}$ and $i_{2}$ is less than $2$ (as in the right part of Figure 1) which allows us to get an approximation guarantee better than $9$. However, this might result in a client contributing to the opening cost of many facilities in $\mathsf{IS}$. Nonetheless, by using the properties of Euclidean metrics, we show that even in this case, we are able to achieve an $\mathsf{LMP}$ approximation guarantee with ratio better than $9$.

Specifically, define $\delta_{\textsf{\tiny mean}}$ to be the constant larger than $2$ that minimizes

Let $d$ be a Euclidean metric on $\cD\cup\cF$ and suppose that $c(j,i)=d(j,i)^{2}$ for every $i\in\cF$ and $j\in\cD$. Then, for any $\lambda\geq 0$, Algorithm JV$(\delta_{\textsf{\tiny mean}})$ constructs a solution $\alpha$ to DUAL$(\lambda)$ and returns a set $\mathsf{IS}$ of opened facilities such that

To simplify notation, we use $\delta$ instead of $\delta_{\textsf{\tiny mean}}$ throughout the proof. Consider any client $j\in\cD$. We shall prove that

Similarly to the proof of Theorem 3.3, the statement then follows by summing up over all clients. A difference compared to the standard analysis of JV is that in our algorithm we may open several facilities in $N(j)$, i.e., client $j$ may contribute to the opening of several facilities. We divide our analysis into the three cases $|N(j)\cap\mathsf{IS}|=1$, $|N(j)\cap\mathsf{IS}|>1$, and $|N(j)\cap\mathsf{IS}|=0$. For brevity, let $S$ denote $N(j)\cap\mathsf{IS}$ and $s=|S|$.

If we let $i^{*}$ be the unique facility in $S$,

In this case, there are multiple facilities in $\mathsf{IS}$ that $j$ is contributing to. We need to show that $\alpha_{j}-\sum_{i\in S}(\alpha_{j}-c(j,i))\geq\frac{1}{\rho_{\textsf{\tiny mean}}}c(j,\mathsf{IS})$.

The sum $\sum_{i\in S}c(j,i)$ is the sum of square distances from $j$ to facilities in $S$ which is at least the sum of square distances of these facilities from their centroid $\mu$, i.e., $\sum_{i\in S}c(j,i)\geq\sum_{i\in S}c(i,\mu)$. Moreover, by the identity, $\sum_{i\in S}c(i,\mu)=\frac{1}{2s}\sum_{i\in S}\sum_{i^{\prime}\in S}c(i,i^{\prime})$, we get

As there is no edge between any pair of distinct facilities $i$ and $i^{\prime}$ in $S\subseteq\mathsf{IS}$, we must have

where the last inequality follows because $j$ is contributing to both $i$ and $i^{\prime}$ and hence $\min(t_{i},t_{i^{\prime}})\geq\alpha_{j}$. By the above,

Now, since $\delta\geq 2$ the above upper bound is a non-increasing function of $s$. Therefore, since $s\geq 2$ we always have

We also know that $\alpha_{j}>c(j,i)$ for any $i\in S$. Therefore, $\alpha_{j}>c(j,\mathsf{IS})$ and, since $\delta\geq 2$:

We conclude the analysis of this case by rearranging the above inequality and recalling that $\rho_{\textsf{\tiny mean}}\geq\frac{1}{\delta/2-1}$.

Here, we claim that there exists a tight facility $i$ such that

To see that such a facility $i$ exists, consider the witness $w(j)$ of $j$. By Claim 3.1, we have $\alpha_{j}\geq t_{w(j)}$ and since $j$ has a tight edge to its witness $w(j)$, $\alpha_{j}\geq c(j,w(j)=d(j,w(j))^{2}$; or, equivalently, $\sqrt{\alpha_{j}}\geq\sqrt{t_{w(j)}}$ and $\sqrt{\alpha_{j}}\geq d(j,w(j))$ which implies that there is a tight facility, namely $w(j)$, satisfying (3.4).

Since $\mathsf{IS}$ is a maximal independent set of $H$, either $i\in\mathsf{IS}$, in which case $d(j,\mathsf{IS})\leq d(j,i)$, or there is an $i^{\prime}\in\mathsf{IS}$ such that the edge $(i^{\prime},i)$ is in $H$, in which case

where the second inequality follows from $d(i,i^{\prime})^{2}=c(i,i^{\prime})\leq\delta\min(t_{i},t_{i^{\prime}})$ by the definition of $H$. In any case, we have by (3.4)

Squaring both sides and recalling that $\rho_{\textsf{\tiny mean}}\geq(1+\sqrt{\delta})^{2}$ completes the last case and the proof of the theorem.

2.3 k𝑘k-Median objective in Euclidean metrics

We use a very similar approach as the one for $k$-means (in Euclidean metrics) to address the $k$-median objective in Euclidean metrics. In this section, we have $c(j,i)=d(j,i)$, i.e., the distances are not squared. Define,

We have $\delta_{\textsf{\tiny med}}\approx 1.633$ and $\rho_{\textsf{\tiny med}}\approx 2.633$.

Let $d$ be a Euclidean metric on $\cD\cup\cF$ and suppose that $c(j,i)=d(j,i)$ for every $i\in\cF$ and $j\in\cD$. Then, for any $\lambda\geq 0$, Algorithm JV$(\delta_{\textsf{\tiny med}})$ constructs a solution $\alpha$ to DUAL$(\lambda)$ and returns a set $\mathsf{IS}$ of opened facilities such that

To simplify notation, we use $\delta$ instead of $\delta_{\textsf{\tiny med}}$ throughout the proof. Similar to the proof of the previous theorem, we proceed by considering a single client $j$ and prove

Let $S$ denote $N(j)\cap\mathsf{IS}$ and $s=|S|$. We again proceed by case distinction on $s$. We first bound the number of cases.

Using the centroid property of squared distances in Euclidean space,

where the last inequality follows from the fact that each pair of facilities $i,i^{\prime}\in S\subseteq\mathsf{IS}$ are not adjacent in $H$ so $d(i,i^{\prime})>\delta\min(t_{i},t_{i^{\prime}})$ and $\min(t_{i},t_{i^{\prime}})\geq\alpha_{j}$ since $i,i^{\prime}\in S\subseteq N(j)$. Since the left-hand-side is upper bounded by $s\alpha_{j}^{2}$, we get $s>\frac{\delta^{2}(s-1)}{2}$. Therefore $s<\frac{\delta^{2}}{\delta^{2}-2}=4$. ∎

We now proceed by considering the cases $s=0,1,2,3$.

Consider the witness $i_{1}$ of $j$. We have $\alpha_{j}\geq t_{i_{1}}$ by Claim 3.1 and also $\alpha_{j}\geq c(j,i_{1})=d(j,i_{1})$. Since $\mathsf{IS}$ is a maximal independent set of $H$, either $i_{1}\in\mathsf{IS}$, in which case $c(j,\mathsf{IS})=d(j,\mathsf{IS})\leq d(j,i_{1})\leq\alpha_{j}$, or there is an $i_{2}\in\mathsf{IS}$ such that the edge $(i_{1},i_{2})$ is in $H$, in which case

In any case, we have $c(j,\mathsf{IS})/\rho_{\textsf{\tiny med}}\leq\alpha_{j}$ as required.

If we let $i^{*}$ be the unique facility in $S$,

where we used the triangle inequality and that $c(i_{1}^{*},i_{2}^{*})>\delta\min(t_{i^{*}_{1}},t_{i^{*}_{2}})\geq\delta\alpha_{j}$ since both $i_{1}^{*}$ and $i_{2}^{*}$ are in $S$ and hence $i_{1}^{*}$ and $i_{2}^{*}$ are not adjacent in $H$. Rearranging the above inequality noting that $\alpha_{j}\geq c(j,\mathsf{IS})$, we have

and the case follows because $\rho_{\textsf{\tiny med}}\geq 1/(\delta-1)$.

using the triangle inequality. Rearranging the above inequality noting that $\alpha_{j}\geq c(j,\mathsf{IS})$, we have

and the lemma follows because $\rho_{\textsf{\tiny med}}\geq 1/(\frac{3\delta}{2}-2)$.

Quasi-polynomial time algorithm

In this section, we present a quasi-polynomial time approach that turns the $\mathsf{LMP}$ approximation algorithms presented in the previous section into approximation algorithms for the original problems ($k$-means and $k$-median), i.e., into algorithms that find solutions satisfying the strict constraint that at most $k$ facilities are opened. This is achieved by only deteriorating the approximation guarantee by an arbitrarily small factor regulated by $\epsilon$. We also introduce several of the ideas used in the polynomial time approach. Although the results obtained in this section are weaker (quasi-polynomial instead of polynomial), we believe that the easier quasi-polynomial algorithm serves as a good starting point before reading the more complex polynomial time algorithm. From now on, we concentrate on the $k$-means problem and we let $\rho=\rho_{\textsf{\tiny mean}}$ denote the approximation guarantee and $\delta=\delta_{\textsf{\tiny mean}}$ denote the parameter to our algorithm, where $\rho_{\textsf{\tiny mean}}$ and $\delta_{\textsf{\tiny mean}}$ are defined as in Section 3.2.2 (it will be clear that the techniques presented here are easily applicable to the other considered objectives, as well). Throughout this section we fix $\epsilon>0$ to be a small constant, and we assume for notational convenience and without loss of generality that $n\gg 1/\epsilon$. We shall also assume that the distances satisfy the following:

By losing a factor $(1+100/n^{2})$ in the approximation guarantee, we can assume that the squared-distance between any client $j$ and any facility $i$ satisfies: $1\leq d(i,j)^{2}\leq n^{6}$, where $n=|\cD|$.

The proof follows by standard discretization techniques and is presented in Appendix C.

Our algorithm will produce a $(\rho+O(\epsilon))$-approximate solution. In the algorithm, we consider separately the two phases of the primal-dual algorithm from Section 3.2.2. Suppose that the first phase produces a set of values $\alpha=\{\alpha_{j}\}_{j\in\mathcal{D}}$ satisfying the following definition:

A feasible solution $\alpha$ of DUAL$(\lambda)$ is good if for every $j\in\mathcal{D}$ there exists a tight facility $i$ such that $(1+\sqrt{\delta}+\epsilon)\sqrt{\smash[b]{\alpha_{j}}}\geq d(j,i)+\sqrt{\smash[b]{\delta t_{i}}}$.

Recall that for a dual solution $\alpha$, $t_{i}$ is defined to be the largest $\alpha$-value out of all clients that are contributing to a facility $i$: $t_{i}=\max_{j\in N(i)}\alpha_{j}$ where $N(i)=\{j\in\cD:\alpha_{j}-d(i,j)^{2}>0\}$.

Two solutions $\alpha$ and $\alpha^{\prime}$ are close if $|\alpha^{\prime}_{j}-\alpha_{j}|\leq\frac{1}{n^{2}}$ for all $j\in\cD$.

We first describe our procedure for generating a close sequence of good solutions. Select the following parameters

We also use the notion of buckets that partition the real line:

We say that $B(v)$ is the index of the bucket containing $v$.

The buckets will be used to partition the $\alpha$-values of the clients. As, in every constructed solution $\alpha$, each client will have a tight edge to a facility, Lemma 4.1 implies that $\alpha_{j}$ will always be at least $1$. Therefore, the definition gives the property that the $\alpha$-values of any two clients $j$ and $j^{\prime}$ in the same bucket differ by at most a factor of $1+\epsilon$. In other words, the buckets will be used to classify the clients according to similar $\alpha$-values.

Note that this implies that each solution in our sequence is good. Indeed, consider a dual solution $\alpha$ that satisfies Invariant 1. Then, for any client $j$, we have some $i$ ($=w(j)$) such that $\sqrt{\alpha_{j}}\geq d(i,j)$ (since $j$ has a tight edge to $w(j)$) and $\sqrt{(1+\epsilon)\delta\alpha_{j}}\geq\sqrt{\delta t_{i}}$ where we used that $B(\alpha_{j})\geq B(t_{i})$ implies $(1+\epsilon)\alpha_{j}\geq t_{i}$. Hence,

and so $\alpha$ is good (here, for the first inequality we have used that $\sqrt{1+\epsilon}\leq 1+\epsilon/2$ and $\sqrt{\delta}\leq 2$). We observe that our initial solution $\alpha^{(0)}$ has $t_{i}=0$ for all $i\in\cF$, and so Invariant 1 holds trivially. In our following analysis, we will show that each call to Sweep preserves Invariant 1.

We now formally describe the procedure QuasiSweep that, given the last previously generated solution $\alpha^{\mbox{\tiny{in}}}$ in our sequences produces the solution $\alpha^{\mbox{\tiny{out}}}$ returned next.

We initialize the algorithm by setting $\alpha_{j}=\alpha^{\mbox{\tiny{in}}}_{j}$ for each $j\in\cD$ and by increasing the opening prices of each facility from $\lambda$ to $\lambda+\epsilon_{z}$. At this point, no facility is tight and therefore the solution $\alpha$ is not a good solution of DUAL$(\lambda+\epsilon_{z})$. We now describe how to modify $\alpha$ to obtain a solution $\alpha^{\mbox{\tiny{out}}}$ satisfying Invariant 1 (and hence into a good solution). The algorithm will maintain a current set $A$ of active clients and a current threshold $\theta$. Initially, $A=\emptyset$, and $\theta=0$. We slowly increase $\theta$ and whenever $\theta=\alpha_{j}$ for some client $j$, we add $j$ to $A$. While $j\in A$, we increase $\alpha_{j}$ at the same rate as $\theta$. We remove a client $j$ from $A$, whenever the following occurs:

$j$ has a tight edge to some tight facility $i$ with $B(\alpha_{j})\geq B(t_{i})$. In this case, we say that $i$ is the witness of $j$.

Note that if a client $j$ satisfies this condition when it is added to $A$, then we remove $j$ from $A$ immediately after it is added. In this case, $\alpha_{j}$ is not increased.

Increasing the $\alpha$-values for clients in $A$, may cause the contributions to some facility $i$ to exceed the opening cost $\lambda+\epsilon_{z}$. To prevent this from happening, we also decrease every value $\alpha_{j}$ with $B(\alpha_{j})>B(\theta)$ at a rate of $|A|$ times the rate that $\theta$ is increasing. Observe that while there exists any such $j\in N(i)$, the total contribution of the clients toward opening this $i$ cannot increase, and so $i$ cannot become tight. It follows that once any facility $i$ becomes tight, $B(\alpha_{j})\leq B(\theta)$ for every $j\in N(i)$ and so $i$ is presently a witness for all clients $j\in N(i)\cap A$. At this moment all such clients in $N(i)\cap A$ will be removed from $A$ and their $\alpha$-values will not subsequently be changed. Thus, $i$ remains tight until the end of QuasiSweep. Moreover, observe any other client $j^{\prime}$ that is added to $A$ later will immediately be removed from $A$ as soon as it has a tight edge to $i$. Thus, neither $t_{i}$ nor the total contribution to $i$ change throughout the remainder of QuasiSweep. In particular, $i$ remains a witness for all such clients $j$ for the remainder of QuasiSweep.

We stop increasing $\theta$ once every client $j$ has been added and removed from $A$. The procedure QuasiSweep then terminates and outputs $\alpha^{\mbox{\tiny{out}}}=\alpha$. As we have just argued, the contributions to any tight facility can never increase, and every client that is removed from $j$ will have a witness through the rest of QuasiSweep (in particular, in $\alpha^{\mbox{\tiny{out}}}$). Thus, $\alpha^{\mbox{\tiny{out}}}$ is a feasible solution of DUAL$(\lambda+\epsilon_{z})$ in which every client $j$ has a witness $w(j)$, i.e., $j$ has a tight edge to the tight facility $w(j)$ and $B(t_{w(j)})\leq B(\alpha_{j})$. Hence, the output of Sweep always satisfies Invariant 1.

This completes the description of QuasiSweep. For a small example of its execution see Figure 2. We now show that the produced sequence of solutions is close and to analyze the running time.

1.2 Closeness and running time

We begin by showing that QuasiSweep produces a close sequence of solutions.

For each client $j\in\cD$, we have $|\alpha^{\mbox{\tiny{in}}}_{j}-\alpha^{\mbox{\tiny{out}}}_{j}|\leq 1/n^{2}$.

We first note that the largest $\alpha$-value at any time is at most $(\lambda+\epsilon_{z})+n^{6}\leq L\epsilon_{z}+n^{6}=4n^{7}+n^{6}\leq 5n^{7}$. This follows from the feasibility of $\alpha$ because, by Lemma 4.1, no squared-distance is larger than $n^{6}$ and the opening cost of any facility is at most $\lambda+\epsilon_{z}\leq L\epsilon_{z}$. Hence, $B(\alpha_{j})\leq 1+\lfloor\log_{1+\epsilon}(5n^{7})\rfloor\leq 10\log_{1+\epsilon}(n)$ for any client $j$ and dual solution $\alpha$.

Any $\alpha_{j}$ can increase by at most $\epsilon_{z}n^{3b}$ while $B(\theta)\leq b$.

The proof is by induction on $b=0,1,\dots,10\log_{1+\epsilon}(n)$.

This case is trivially true because there are no clients $j$ with $B(\alpha_{j})=0$, and so no clients can have been added to $A$ while $B(\theta)=0$. Indeed, any client $j$ had a tight edge to some facility in $\alpha^{\mbox{\tiny{in}}}$, which by Lemma 4.1 implies $\alpha^{\mbox{\tiny{in}}}_{j}\geq 1$, and a client’s $\alpha$-value can decrease only while some smaller $\alpha$-value is increasing.

Now, we suppose some $\alpha_{j}$ is increasing while $B(\theta)\leq b$. Note that we then must have $\alpha_{j}=\theta$. Let $i$ be the witness of $j$ in $\alpha^{\mbox{\tiny{in}}}$, and let $N^{\mbox{\tiny{in}}}(i)$ be the set of clients contributing to $i$ in $\alpha^{\mbox{\tiny{in}}}$. We further suppose that $\alpha_{j}$ is increased by at least $\epsilon_{z}$ while $B(\theta)\leq b$; otherwise the claim follows immediately, since $\epsilon_{z}\leq n^{3b}\epsilon_{z}$ for all $b\geq 0$.

First, suppose that $\alpha_{j}<\alpha^{\mbox{\tiny{in}}}_{j}$ and so $\alpha_{j}$ was previously decreased by QuasiSweep. Moreover, since $\alpha_{j}$ has increased by at least $\epsilon_{z}$ while $B(\theta)\leq b$, we must have previously decreased $\alpha_{j}$ while $B(\alpha_{j})\leq b$. In particular, at the last moment $\alpha_{j}$ was decreased, we must have had $B(\alpha_{j})\leq b$, and since $\alpha_{j}$ was decreasing at this moment, we also had $B(\theta)<B(\alpha_{j})$. Then, $\alpha_{j}$ was decreased only while $B(\theta)<b$. Moreover, during this time, $j$’s $\alpha$-value was decreased at a rate of $|A|$, and so was decreased (in total) at most $n$ times the largest amount that any other $\alpha_{j^{\prime}}$ was increased. By the induction hypothesis, any $\alpha_{j^{\prime}}$ was increased at most $\epsilon_{z}\cdot n^{3b-3}$ while $B(\theta)<b$, and so $\alpha_{j}$ was decreased at most $\epsilon_{z}\cdot n^{3b-2}$. Thus, after $\alpha_{j}$ increases by at most $\epsilon_{z}\cdot n^{3b-2}$ we will again have $\alpha_{j}=\alpha^{\mbox{\tiny{in}}}_{j}$.

Now, we consider how much $\alpha_{j}$ may increase while $\alpha_{j}\geq\alpha^{\mbox{\tiny{in}}}_{j}$ (and still $B(\theta)\leq b$). For each $j^{\prime}\in N^{\mbox{\tiny{in}}}(i)$ we must have initially had $B(\alpha^{\mbox{\tiny{in}}}_{j^{\prime}})\leq B(\alpha^{\mbox{\tiny{in}}}_{j})$ since $i$ is a witness for $j$ in $\alpha^{\mbox{\tiny{in}}}$. Additionally, by our assumptions in this case, $B(\alpha^{\mbox{\tiny{in}}}_{j})\leq B(\alpha_{j})\leq b$. Thus, the $\alpha$-value of any $j^{\prime}\in N^{\mbox{\tiny{in}}}(i)$ was decreased by QuasiSweep only while $B(\theta)\leq b-1$ and so by the same argument as above the $\alpha$-value of any $j^{\prime}\in N^{\mbox{\tiny{in}}}(i)$ can decrease at most $\epsilon_{z}n^{3b-2}$ throughout QuasiSweep. Thus, the total contribution to $i$ from all $j^{\prime}\neq j$ can decrease at most $(n-1)\cdot\epsilon_{z}\cdot n^{3b-2}$. After increasing $\alpha_{j}$ at most $(n-1)\cdot\epsilon_{z}\cdot n^{3b-2}+\epsilon_{z}$, $i$ will again be tight. Moreover, at this moment, any client $j^{\prime}$ contributing to $i$ was either already added to $A$ (and potentially also removed) in which case $\alpha_{j^{\prime}}\leq\theta=\alpha_{j}$ or it was not already added to $A$ in which case $\alpha_{j^{\prime}}\leq\alpha^{\mbox{\tiny{in}}}_{j^{\prime}}$ since $\alpha_{j^{\prime}}$ has not been increased yet. In either case, $B(\alpha_{j^{\prime}})\leq B(\alpha_{j})$ so $i$ is a witness for $j$, and $j$ will be removed from $A$.

Altogether, the total amount $\alpha_{j}$ can increase while $B(\theta)\leq b$ is then the sum of these two increases, which is $\epsilon_{z}\cdot n^{3b-2}+(n-1)\cdot\epsilon_{z}\cdot n^{3b-2}+\epsilon_{z}\leq\epsilon_{z}\cdot n^{3b}$, as required. ∎

The claim immediately bounds the increase $\alpha^{\mbox{\tiny{out}}}_{j}-\alpha^{\mbox{\tiny{in}}}_{j}$ by $\tfrac{1}{n^{3}}\leq\tfrac{1}{n^{2}}$ as required (recall that $\epsilon_{z}=n^{-3-10\log_{1+\epsilon}n}$)). Moreover, as shown in the proof of the claim above, the $\alpha$-value of every client decreases by no more than $n$ times the maximum increase in the $\alpha$-value of any client. Then, the desired bound $\tfrac{1}{n^{2}}$ on $\alpha^{\mbox{\tiny{in}}}_{j}-\alpha^{\mbox{\tiny{out}}}_{j}$ follows as well. ∎

For the sake of clarity, we have presented the QuasiSweep procedure in a continuous fashion. We show in Appendix A how to implement QuasiSweep as a discrete algorithm running in polynomial time. We conclude the analysis of this section by noting that, as Sweep is repeated $L=n^{O(\epsilon^{-1}\log n)}$ times, the total running time for producing the sequence $\alpha^{(0)},\alpha^{(1)},\dots,\alpha^{(L)}$ is $n^{O(\epsilon^{-1}\log n)}$.

2 Finding a solution of size k𝑘k

2.2 Analysis

The total running time is $n^{O(\epsilon^{-1}\log n)}$ since the number of steps $L$ (and the number of dual solutions in our sequence) is $n^{O(\epsilon^{-1}\log n)}$ and each step runs in polynomial time since it involves the construction of at most $O(|\cF|)$ conflict graphs and maximal independent sets.

For any $j\in\cD_{>0}$, $d(j,\mathsf{IS})^{2}\leq\rho\cdot\left(\alpha_{j}-\sum_{i\in S_{j}}\beta_{ij}\right)$.

Consider some $j\in\cD_{>0}$ and first suppose that $|S_{j}|=1$. Then, if we let $S_{j}=\{i\}$, $\alpha_{j}=\beta_{ij}+d(j,i)^{2}\geq\beta_{ij}+d(j,\mathsf{IS})^{2}$ just as in “Case $s=1$” of Theorem 3.4. Next, suppose that $|S_{j}|=s>1$. In other words, $j$ is contributing to multiple facilities in $\mathsf{IS}$. By construction we have $\alpha_{j}\leq\min(t_{i},t_{i^{\prime}})$ for any two facilities $i,i^{\prime}\in S_{j}$. Thus, $\alpha_{j}-\sum_{i\in S_{j}}\beta_{ij}\geq\frac{1}{\rho}d(j,\mathsf{IS})^{2}$ by the exact same arguments as in “Case $s>1$” of Theorem 3.4. ∎

Next, we bound the total service cost of all those clients that do not contribute to any facility in $\mathsf{IS}$. The proof is very similar to “Case $s=0$” in the proof of Theorem 3.4.