1 Introduction↩︎

In a generic clustering problem we are given a collection of points together with a dissimilarity measure (i.e. distance function) between pairs of points. The high-level goal is to partition the points into a “small” number of clusters so that similar points are clustered together while dissimilar ones are clustered separately. One of the most fundamental and best-studied clustering problems is \(k\)-Means. Here we are given a collection \(D\) of \(n\) points (or clients) and a collection \(F\) of centers (or facilities), as well as a integer \(k>0\). We are also given metric distances \(d:(D\cup F) \times (D\cup F)\rightarrow \mathbb{R}_{\geq 0}\). Our goal is to select a set \(S\) of \(k\) centers (the open centers) so as to minimize the sum of the squared distances from each client to the closest open center², i.e., \[\sum_{j\in D}d^2(j,S)=\sum_{j\in D}\min_{i\in S}d^2(j,i).\] Observe that each feasible solution \(S\) naturally induces \(k\) clusters where the cluster associated with \(i\in S\) is given by the clients that are closer to \(i\) than to any other center in \(S\) (breaking ties arbitrarily). We next use \(\text{OPT}_k\) to denote a reference optimal solution and \(\text{opt}_k\) to denote its cost. When \(k\) is clear from the context, we simply use \(\text{OPT}\) and \(\text{opt}\).

\(k\)-Means is NP-hard and well-studied in terms of approximation algorithms. It is impossible to approximate it (in polynomial time) better than a factor \(1+8/e \approx 3.94\) [1]. The current-best (polynomial-time) approximation factor for this problem is \(9+\varepsilon\) for any constant \(\varepsilon>0\) by Ahmadian, Norouzi-Fard, Svensson, and Ward [2]. Their algorithm is based on a primal-dual Lagrangian Multiplier Preserving (LMP) \(9\)-approximation algorithm for facility location with squared metric connection costs (similar in spirit to a classical primal-dual \(3\)-approximation algorithm by Jain and Vazirani [3] for metric connection costs), plus a careful way to combine so-called bi-point solutions that introduces a \(1+\varepsilon\) factor only in the approximation (rather than a larger constant factor as in prior work). This improved on an earlier \(25\)-approximation by Gupta and Tangwongsan [4] based on local search: we will need this result later. Interestingly a \((1+8/e+\varepsilon)\)-approximation can be obtained in FPT time [5].

Very often in practice one considers the Euclidean special case of \(k\)-Means (which in the literature is often just called \(k\)-Means), where the clients are \(n\) points in the (\(d\)-dimensional) Euclidean space, with the respective distances. In this case one is allowed to select any point in the Euclidean space as a center. However known reductions [6]–[8] allow one to focus on a set of candidate centers \(F\) of size \(n^{O_\varepsilon(1)}\) while introducing a factor \(1+\varepsilon\) in the approximation for any constant \(\varepsilon>0\). Therefore one can (essentially) see the Euclidean case as a special case of the metric one. The best-known approximation factor for the Euclidean case is \(5.92\) by Cohen-Addad, Esfandiari, Mirrokni, Narayanan [9]. This improves on a sequence of increasingly better approximations for the problem: a \((9+\varepsilon)\)-approximation by Kanungo, Mount, Netanyahu, Piatko, Silverman, and Wu [10], a \(6.36\)-approximation by Ahmadian et al. [2], and a \(6.13\)-approximation by Grandoni, Ostrovsky, Rabani, Schulman, and Venkat [11].

Our main result is a \(5.83\)-approximation for (metric) \(k\)-Means. More precisely, we get the following:

We remark that the above result not only substantially improves the current best \((9+\varepsilon)\)-approximation for the general metric case [2], but even slightly the current best \(5.92\)-approximation for the Euclidean case [9]. We overview our approach in Section 2.

1.1 Related Work↩︎

\(k\)-Means belongs to the family of \(k\)-clustering problems where the target number of clusters \(k\) is fixed. Other famous examples are \(k\)-Center and \(k\)-Median. In \(k\)-Center one wishes to select a set \(S\) of \(k\) centers so as to minimize the maximum distance from any client to \(S\), i.e., the objective function to minimize is \(\max_{j\in D}d(j,S)\). This problem admits a simple greedy \(2\)-approximation which is best possible unless \(P=NP\) [12], [13]. \(k\)-Median is defined like \(k\)-Means, except that here one wishes to minimize the sum of the distances rather than squared distances, i.e., the objective function is \(\sum_{j\in D}d(j,S)\). \(k\)-Median is very close to \(k\)-Means in terms of results and techniques. \(k\)-Median is hard to approximate below a factor \(1 + 2/e\) [1]. For general metric distances, the first constant approximation was achieved by Charikar, Guga, Tardos and Shmoys [14]. After a very long sequence of improvements [15]–[21], the current-best \((2+\varepsilon)\)-approximation for this problem was very recently achieved by Cohen-Addad, Grandoni, Lee, Schwiegelshohn, and Svensson [22]. This is also the best result for the Euclidean case, improving on an earlier \(2.406\)-approximation [9].

2 Overview of our Approach↩︎

At a very high level, our approach looks similar in spirit to the one leading to a recent \((2+\varepsilon)\)-approximation for the related \(k\)-Median problem by Cohen-Addad, Grandoni, Lee, Schwiegelshohn, and Svensson [22]. In more detail, we exploit a combination of two different algorithms. The first one is a bicriteria approximation algorithm with the desired approximation factor, which opens \(O(\log n/\varepsilon^2)\) more centers than the \(k\) allowed ones.

The second algorithm is a radically different \((5+\varepsilon)\)-approximation algorithm that only works for instances which are stable in the following sense. We say that an instance of \(k\)-Means is \(\beta\)-stable if \(\text{opt}_{k-1}\geq (1+\beta)\text{opt}_k\), where \(\text{opt}_{k-1}\) is the optimal cost for the same instance but with the target number of centers being \(k-1\).

It is relatively easy to derive Theorem 1 from the above two theorems.

It remains to describe how the above two theorems are obtained. The algorithm from Theorem 3 is rather complex and its analysis is highly non-trivial (see Section 6). However, it is a relatively easy adaptation of the \((2+\varepsilon)\)-approximation algorithm in [22] for a similar notion of stable \(k\)-Median instances. The main difference is that we have to carefully use an approximate form of triangle inequality (see Lemmas 2 and 3). This also explains why we get a higher approximation factor \(5+\varepsilon\).

The main contribution of this paper is our proof of Theorem 2. In more detail, we consider the related facility location problem (with uniform opening costs). Recall that here, instead of a bound \(k\) on the number of centers, we are given a uniform facility cost \(f\). We are now allowed to open an arbitrary number of centers/facilities \(S\), and the objective function is the total cost of the open facilities plus the squared distance from each client to the closest (open) facility in \(S\), i.e., \[\text{cost}_{FL}(S):=f|S|+\sum_{j\in D}d^2(j,S)\,.\] We say that an algorithm for the above problem is LMP \(\Gamma\)-approximate if it produces a feasible solution \(S\) such that \[\Gamma\cdot f|S|+\sum_{j\in D}d^2(j,S)\leq \Gamma\cdot \text{opt}_{LP}(f). \label{eq:LMP}\tag{1}\] where \(\text{opt}_{LP}(f)\leq \text{opt}\) is the optimal cost of a standard LP relaxation for the problem (see Section 3). In other words, the solution cost at most \(\Gamma \cdot \text{opt}\) even if we increase the facility cost of \(S\) by a factor \(\Gamma\).

Recall that the current-best \((9+\varepsilon)\)-approximation for (Metric) \(k\)-Means by Ahmadian et al. [2] builds upon a primal-dual LMP \(9\)-approximation for facility location (with squared metric connection costs), which is inspired by a classical primal-dual \(3\)-approximation by Jain and Vazirani [3] for the case of metric connection costs. Our main contribution is a greedy LMP \((3+2\sqrt{2})\)-approximation for facility location with squared metric connection costs which is a non-trivial (but relatively simple) variant of the classical greedy LMP \(2\)-approximation (JMMSV) for facility location with metric connection costs by Jain, Mahdian, Markakis, Saberi, and Vazirani [16]. It is instructive to recall how JMMSV works in order to see the differences.

2.0.0.1 The JMMSV Greedy Algorithm.

The algorithm has a variable \(\alpha_j\) per client \(j\) (initialized to \(0\)), a set \(S\) of open facilities (initialized to \(\emptyset\)), and a set \(A\) of active clients (initialized to \(D\) ). At each point of time, the bid \(bid(j,i)\) of client \(j\) towards facility \(i\) is \([\alpha_j-d(j,i)]^+\) if \(j\) is active and \([d(j,S)-d(j,i)]^+\) otherwise³. The variables \(\alpha_j\) of active clients are increased uniformly until one of the following events happens:

1. For some client \(j\) and \(i\in S\), \(\alpha_j\geq d(j,i)\). In that case \(j\) is removed from \(A\) and we say that \(j\) is connected to \(i\);

2. For some (not open) facility \(i\notin S\), one has \(\sum_{j\in D}bid(j,i)=f\).⁴ In that case we open \(i\), i.e., we add \(i\) to \(S\). Furthermore, each \(j\in A\) with \(\alpha_j\geq d(j,i)\) is removed from \(A\) (and connected to \(j\)). Also, each inactive client \(j\in D-A\) with \(bid(j,i)>0\) is reconnected to \(i\).

The intuition behind the bids is as follows. On the one hand, the active clients \(j\) (which are not yet connected to an open facility) are willing to pay the difference \(\alpha_j-d(j,i)\) (if positive) to open \(i\), while reserving \(d(j,i)\) for their own connection cost to \(i\). On the other hand, the inactive clients \(j\in D-A\) which are already connected to the closest facility in \(S\), are willing to offer \(d(j,S)-d(j,i)\) (if positive) towards the opening of facility \(i\): if then \(i\) is actually opened, \(j\) reconnects to \(i\) which is closer (while altogether still spending \(\alpha_j\) in total). This also motivates the term greedy. Notice that at any point of time active clients \(j\) satisfy \(\alpha_j<{d}(j,S)\).

2.0.0.2 Counterexample for the Naïve Extension.

A naïve idea is to adapt the above algorithm to the case of squared distances by simply replacing \(d(\cdot,\cdot)\) with \(d^2(\cdot,\cdot)\). It is instructive to see why this attempt fails miserably (which also explains why prior work [2], [9] use variants of the Jain-Vazirani algorithm [3] instead, despite the worse approximation factor). Consider the example in Figure 1, where we have two facilities \(i\) and \(i'\), and two clients \(1\) and \(2\), all on a line, with distances specified in the figure. The facility opening cost is \(f\). Consider the following (hypothetical but valid) execution of the algorithm.

• Time \(f/2+0.25-\varepsilon\): facility \(i'\) is open by client \(2\) and other clients not in the figure.

  • At that time, client \(1\) and \(2\)’s bid to \(i\) was \(f/2 - \varepsilon\) each (so \(i\) was very close to be opened), but after the opening of \(i'\), client \(2\) lowers its bid to \(i\) to \((\sqrt{f/2}-0.5)^2 - 0.5^2= f/2-\sqrt{f/2}\).

• Time \(f/2+\sqrt{f/2} + 0.25\): facility \(i\) is open as the bid from client \(1\) becomes \(f/2+\sqrt{f/2} + 0.25 - 0.25 = f/2 + \sqrt{f/2}\).

• To prove an LMP \(\Gamma\)-approximation, the dual constraint corresponding to a solution where \(1\) and \(2\) are connected to \(i\) (paying \(f\) for opening and \(2 \cdot 0.5^2 = 0.5\) for connecting) requires that \(\alpha_1+\alpha_2 \leq f + 0.5\Gamma\). Our \(\alpha_1+\alpha_2\) is \(f + \sqrt{f/2} + 0.5-\varepsilon\), which is only an LMP \(\Omega(\sqrt{f})\)-approximation.

In this example, client \(2\) contributes to the opening of some center far away from it, with a contribution much less than its \(\alpha\) value. Or equivalently, the connection cost \(d^2(2,i')\) is strictly but only slightly smaller than \(\alpha_2\). This is the key (and essentially the only) reason why this naïve variant of JMMSV does not work; one can see that any other assumption between \(\alpha_2\) and \(d^2(2, i')\) in the above example (realized by moving \(i'\) and clients not in the figure) would lead to a good LMP approximation as follows.

• If \(d^2(2,i')=\alpha_2\), then \(2\)’s bid to \(i\) is \(\alpha_2 - 0.25\) even after the opening of \(i'\), so the fact that \(i\) was open at time \(\alpha_1-\varepsilon\) (for infinitesimally small \(\varepsilon> 0\)) implies \(\alpha_1+\alpha_2 \leq f + 0.5\).

• If \(d^2(2,i')\leq (1-\delta)\alpha_2\) for some constant \(\delta > 0\), we could simply use an approximate version of the triangle inequality to prove a good LMP approximation. For example, if \(\delta=1/2\), we would have \(\alpha_1 \leq d^2(1,i') \leq 2d^2(2,i') + 4(d^2(1,i)+d^2(2,i)) \leq \alpha_2 + 2\). Since \(\alpha_2 \leq f/2+0.25\) in this example, this would give us \(\alpha_1+\alpha_2 = f + O(1)\).

2.0.0.3 Our Algorithm.

This inspired us to construct a different variant of the JMMSV algorithm in the following spirit: if a client can contribute to the opening of a facility, its \(\alpha\) value must be much higher than its connection cost to the facility (e.g., at least twice as high as it). In more detail, our LMP algorithm uses the following different bidding strategy. Let \(\gamma>1\) be a parameter to be fixed later. The bid of active clients is \([\alpha_j-\gamma d^2(j,i)]^+\). In particular notice that \(j\) cannot contribute to the opening of \(i\) if \(\alpha_j<\gamma d^2(j,i)\) (while one might naturally expect it for \(\alpha_j>d^2(j,i)\)). However, whenever \(\alpha_j\geq d^2(j,i)\) for some already open facility \(i\in S\) and active clients \(j\), we still connect \(j\) to \(i\) and make \(j\) inactive. In the above case, if \(\gamma d^2(j,i)>\alpha_j\), we say that \(j\) is indirectly connected to \(i\) since \(j\) did not contribute to the opening of \(i\) with a positive bid. The inactive clients \(j\) that are connected to an open facility \(i\) whose opening they contributed to, are called directly connected. The bids for directly and indirectly connected clients in the later steps are different: they are \([\alpha_j-\gamma d^2(j,i)]^+\) for the indirectly connected ones and \([\gamma d^2(j,S)-\gamma d^2(j,i)]^+\) for the directly connected ones. When an indirectly connected client contributes to the opening of a facility, it becomes directly connected. For the above process, we are able to show that the final solution \(S\) satisfies: \[\sum_{j\in D}\alpha_j\geq f|S|+\sum_{j\in D}d^2(j,S),\] i.e., the cost of the solution is upper bounded by the sum of the variables \(\alpha_j\). Furthermore, we can establish the following inequality showing that \(\alpha_j\)’s, after scaling, form a dual-feasible solution. \[\sum_{j\in D} \left[\alpha_j - \left(\gamma+2+\frac{2}{\gamma-1}\right)d^2(j,i)\right]^+ \leq f,\quad \forall i\in F.\] The above equations together imply that the overall algorithm is an LMP \((\gamma+2+\frac{2}{\gamma-1})\)-approximation for the considered facility location problem. Fixing \(\gamma=1+\sqrt{2}\) gives the claim.

2.0.0.4 Opening \(k\) Centers.

With the above LMP approximation for facility location at hand, we can now derive the result in Theorem 2 following an approach again close to [22]. A standard way to derive a \(\rho_{kMeans}\)-approximation algorithm for \(k\)-Means from an LMP \(\rho_{FL}\)-approximation for facility location with squared metric connection costs is to perform a binary search over the uniform facility cost \(f\) so as to obtain a solution with \(k_1<k\) open facilities for some facility cost \(f_1\), and a solution with \(k_2>k\) open facilities for some facility cost \(f_2<f_1\), where the difference between \(f_1\) and \(f_2\) is very tiny. This is called a bi-point solution. Then the two solutions are combined together to obtain a solution opening \(k\) facilities. This combination is typically costly, hence leading to a factor \(\rho_{kMeans}\) substantially larger than \(\rho_{FL}\)⁵. We instead carefully modify the mentioned LMP \((3+2\sqrt{2})\)-approximation for facility location (while increasing by \(\varepsilon\) the approximation factor) so that the variables \(\alpha_j\) are increased in a logarithmic number of rounds only (with a small multiplicative increase at each round). Exploiting this fact and the walking-between-solutions framework in [22], we are able to construct a bi-point solution where \(k_2\leq k+O(\log n/\varepsilon^2)\). This leads to Theorem 2. Again, the details are rather technical but they are similar in spirit to [22].

3 Preliminaries↩︎

3.0.0.1 LP Relaxations and Basic Assumptions.

Given \((D, F, d)\), the standard LP relaxations for \(k\)-Means and facility location with squared distances and uniform opening cost \(f\) are as follows.

Let \(\text{opt}_{LP}(k)\) and \(\text{opt}_{LP}(f)\) be the optimal values for \(LP_{km}\) and \(LP_{FL}(f)\) respectively. Recall that \([a]^+ := \max(a, 0)\). The dual for \(LP_{FL}(f)\) is: \[\begin{align} \max \quad & \sum_{j\in D}\alpha_j & (DP_{FL}(f))\\ s.t. \quad & \sum_{j\in D} [\alpha_j - d^2(i,j)]^+ \leq f & \forall i\in F \\ & \alpha \geq 0. \end{align}\]

Thanks to standard reductions, we can assume that distances are integers in \([1,n^3/\varepsilon]\) while loosing a factor \(1+O(\varepsilon)\) in the approximation.

3.0.0.2 Triangle Inequalities for Squared Distances.

Given three distances \(d,d_1,d_2\) with \(d\leq d_1+d_2\), the following lemma provides a useful upper bound on \(d^2\) as a function of \(d_1\) and \(d_2\). The next lemma provides a similar upper bound in an analogous setting where \(d\leq d_1+d_2+d_3\).

4 Greedy Algorithm↩︎

We formally present the new greedy algorithm for facility location with squared distances outlined in 2. Let \(\gamma=1+\sqrt{2}\) and \(\Gamma = \gamma+2+\frac{2}{\gamma-1}\approx 5.828\) be the desired approximation ratio. Given an instance \((D, F, d, f)\) for facility location, let \(\hat{f}=\Gamma f\). The following is our new greedy algorithm⁶ achieving an LMP \(\Gamma\)-approximation (see Algorithm 2).

We say that the clients in \(A\) are active, that the clients in \(IC\) are indirectly connected, and that the clients in \(DC\) are directly connected. The following standard lemma proves that these \(\alpha\) values are enough to pay the connection costs and the (augmented) opening costs.

Let \((\alpha_j)_{j\in D}\) be the final vector of \(\alpha\) values. We following lemma is the core of our analysis: it implies that \(\alpha / \Gamma\) is a feasible solution for \(DP_{FL}(f)\). Combined with 4, it shows that our solution \(S\) satisfies \[\sum_{j \in D} d^2(j, S) + \Gamma |S| f \leq \sum_{j \in D} \alpha_j \leq \Gamma \cdot \text{opt}_{LP}(f) \leq \Gamma \cdot \text{opt}_{FL}(f),\] which finishes the proof of the LMP \(\Gamma\)-approximation.

5 \(5.83\)-Approximation with \(O(\log n/\varepsilon^2)\) Extra Centers↩︎

Let \(\gamma=1+\sqrt{2}\) and \(\Gamma = \gamma+2+\frac{2}{\gamma-1}\approx 5.828\). In this section, we present a polynomial-time algorithm that opens \(k + O(\log n/ \varepsilon^2)\) many centers and returns a solution of cost at most \((\Gamma + \varepsilon)\text{opt}\), proving 2.

This section is close to Section 3 of [22] that used the classical LMP 2-approximation [16] for a \((2+\varepsilon)\)-approximation for \(k\)-Median. Here, we start from the LMP \(\Gamma\)-approximation (with squared distances) presented in 4 instead. 5.1 presents a modified algorithm LogAdaptiveAlgorithm that runs in \(O(\log n / \varepsilon^2)\) adaptive steps (that we call phases) and 5.2 analyzes the algorithm. Finally, building on additional setups in 5.3, 5.4 presents our final algorithm MergeSolutions achieving the guarantee of 2.

5.1 Log Adaptive Algorithm↩︎

We present the modified greedy algorithm that works in \(O(\log n / \varepsilon^2)\) phases. At high level, instead of considering the continuous time where \(\alpha\) values of the active clients are increased at unit rate, the idea is to consider discrete time steps (that we call phases), where the \(i\)th phase \((i \geq 0)\) corresponds to time \((1+\varepsilon^2)^i\). Then, since we assumed that the minimum nonzero distance is \(1\) and the maximum distance is \(\poly(n)\), the number of phases is \(O(\log(n)/\varepsilon^2)\).

We will ensure that the dual feasibility (5) will be satisfied without any slack, which could be violated if \(\alpha\) values of the active clients discretely jump from one value to a strictly larger value. In order to deal with this issue, [22] uses two ideas, namely (1) scaling down the metric by a factor of \((1-\delta)\) with \(\delta = 3\varepsilon\) and (2) allowing clients extremely close to a facility \(i\) to continuously increase their \(\alpha\) values in order to open \(i\). We use a similar definition adapted to our new greedy algorithm for squared distances and its analysis. A key ingredient is the following technical definition of facilities that can be opened.

The openability can be easily checked by solving a linear program.

3 is our log-adaptive algorithm. It works in \(O(\log n/\varepsilon^2)\) phases where, in the \(i\)th phase, every active client \(j\) has \(\alpha_j = (1+\varepsilon^2)^{i}\) and the algorithm opens openable facilities in an arbitrary order by increasing \(\alpha_j\)’s to \(\tau_j\)’s. However, when \(i\) is open, the only scenario \(j\) with \(\tau_j > \alpha_j\) is when \(j \in B(i, \sqrt{\varepsilon\theta})\), so \(j\) is immediately put into \(DC\) after \(i\) is open.

5.2 Analysis↩︎

In this subsection, we show that 3 achieves an LMP \({\frac{\Gamma}{1-\delta}=}\Gamma + O(\varepsilon)\) approximation. This is a direct consequence of the following two lemmas, proved in 5.2.1 and 5.2.2, resp. Let \((\alpha^*_j)_{j \in D}\) be the final \(\alpha\) values.

In more detail, Lemma 7 shows that the dual values \(\frac{\alpha^*_j}{1-\delta}\) can pay the cost of the solution with facility opening cost \(\frac{\Gamma}{1-\delta}f\). Lemma 8 implies that \(\alpha^*/\Gamma\) is a feasible solution for the dual of the facility location LP with facility cost \(f\). We conclude as desired that \[\frac{\Gamma}{1-\delta}\text{opt}_{\text{LP}}(f)\geq \sum_{j\in D}d^2(j,S)+\frac{\Gamma}{1-\delta}|S|f.\]

5.2.1 Approximation Guarantee↩︎

5.2.2 Dual Feasibility↩︎

Here we prove 8. First, we prove the following lemma, showing intuitively that no facility is overpaid during the course of our algorithm.

Let \((\alpha_j^*)_j\) be the final vector of \(\alpha\)-values. We say that a facility \(i\) is frozen (by \(i_0\)) if there exists some open facility \(i_0\) such that \(d(i,i_0)\leq (1+\sqrt{\Gamma\varepsilon}+3\varepsilon)\sqrt{\theta/\Gamma}\). Notice that a facility \(i\) can be frozen by \(i\) itself. We next focus on a specific facility \(i\), with the goal of proving Lemma 8. Let \(D^*=\{j:\alpha_j^*>\Gamma d^2(i,j)\}\). Similarly to the proof of 5, we define the set \(D^*\) for \(i\) and assume w.l.o.g. that \(D^*=[s]\) and that \(\alpha^*_j\leq \alpha^*_{j+1}\) for \(j=1,\ldots,s-1\).

Next, for each client \(j\in [s]\):

• If \(j\) becomes inactive strictly before \(i\) is frozen, we define \(DC^j\) as the set of clients in \([j-1]\) that are directly connected at the beginning of the stage \(j\) becomes inactive, \(IC^j = [j-1]\setminus DC^j\), and \(S^j\) as the set of opened facilities at the beginning of the stage \(j\) becomes inactive.

• Else, we define \(DC^j\) as the set of clients in \([j-1]\) that are directly connected right before \(j\) becomes inactive, \(IC^j = [j-1] \setminus DC^j\), and \(S^j\) as the set of opened facilities right before \(j\) becomes inactive.

The last two claims together with the fourth bullet point of 1 imply the following fact:

We now are ready to prove Lemma 8. The following proof is almost identical to the final part of the proof of Lemma 5, the main difference being that we use [fct:bidbound] here.

5.3 Walking Between Two Solutions: Setup↩︎

Given the description of the LogAdaptiveAlgorithm, we present our final algorithm MergeSolutions, resulting in a solution that opens \(k+O(\log n{/\varepsilon^2})\) centers by maintaining two partial solutions that lead to opening at least \(k\) and at most \(k\) centers respectively and gradually merge them into a single solution that opens exactly \(k+O(\log n{/\varepsilon^2})\) facilities, exploiting the \(O(\log n{/\varepsilon^2})\) adaptivity of LogAdaptiveAlgorithm.

Before detailing the merging procedure, which we call MergeSolutions, it is helpful to observe that we may generalize the analysis to make it more “robust”. We introduce this more robust analysis in this subsection, where we also define the notation that will be used later to describe MergeSolutions.

Starting from the two solutions obtained in the standard way whose facility costs differ by at most \(\eta := 2^{-n}\), the two solutions we maintain will have an exponentially small difference \(\eta\) in one of their parameters. (This parameter will be either the facility cost \(f\) or some distance related to a free facility that will be introduced later.) For this reason, when we merge the two solutions, it may occur that some facilities in one solution are almost completely paid for, rather than fully paid for. To address this, we generalize the definition of "openable" to "\(\eta\)-openable," where the only difference is in the third condition: we now require each facility that is opened to be paid \(\hat{f} - \eta \gamma n\) instead of \(\hat{f}\) (as in Definition 1). We also note that being openable only depends on the time \(\theta=(1+\varepsilon^2)^p\) and the sequence \(\mathcal{H}=(H_1, \dots , H_p)\) of the sets of facilities open by each stage 1 of the algorithm that was executed so far. Indeed, letting \(S=\cup_{i\leq p}H_i\), we can obtain \(A\), \(IC\) and \(DC\) as follows: \(A=\{j\in D: \theta < (1-\delta)d^2(j, S)\}\), \(IC=\{j\in D\setminus A:\theta<(1-\delta)\gamma d^2(j, S)\}\) and \(DC=D\setminus (IC\cup A)\). Furthermore, the \(\alpha\)-values of the indirectly connected clients only depend on \(\mathcal{H}\). This allows us to simplify notation and we have the following generalization of Definition 1.

We also introduce the concept of free facilities. During the process of merging the two solutions, we introduce free facilities, whose opening costs are not necessarily paid. While we will introduce at most \(O(\log n{/\varepsilon^2})\) such free facilities, we do not restrict the number of free copies a regular facility may have. Let \(S_{\mathsf{free}}\subseteq S\) denote the multiset of free facilities that have been opened.

Moreover, for each free copy \(\tilde{i} \in S_{\mathsf{free}}\) of \(i \in F\), we define a parameter \(u(\tilde{i})\), which represents the distance between the freely opened facility \(\tilde{i}\) and its original copy \(i\). So for any point \(x\) in the metric space (either a facility or a client), the distance is given by \(d(\tilde{i}, x) = \sqrt{u(\tilde{i}) + d^2(i, x)}\). We remark that the space remains metric after this extension. In particular, for any two points \(x_1,x_2\) other than \(\tilde{i}\), one has \(d(x_1,x_2)\leq d(\tilde{i},x_1)+d(\tilde{i},x_2)\).

We further introduce the following terminology: as mentioned, we refer to the facilities in \(S_{\mathsf{free}}\) as free facilities, the facilities in \(F\) and \(S_{\mathsf{reg}}= S - S_{\mathsf{free}}\) as regular facilities, and when we simply mention "facilities," we are referring to their union.

Finally, we will use the letters \(i\) and \(h\) to denote facilities. Specifically, we use \(i \in F\) to denote a regular facility and use \(\tilde{i}\) to denote a free copy of \(i\). When referring to a facility that can be either free or regular, we use the letter \(h\).

We now introduce the definition of valid sequences, which aims to capture the sequence of facilities that are opened during the while-loop of stage 1 in LogAdaptiveAlgorithm. Let us say \(\mathcal{H}'=(H_1',\dots,H_p')\) is a super-sequence of \(\mathcal{H}=(H_1,\dots,H_p)\) and denote it \(\mathcal{H}'\supseteq \mathcal{H}\) if they have a common prefix of length \(p-1\), i.e. for all \(i<p\), \(H_i'=H_i\), and \(H_p'\) is a superset of \(H_p\). Given two sequences \(\mathcal{H}_1=(H_1,\dots,H_r)\) and \(\mathcal{H}_2=(H_{r+1}, \dots,H_l)\), we define \(\mathcal{H}_1\oplus\mathcal{H}_2=(H_1,\dots,H_l)\) as the concatenation of the two sequences.

The reader might wonder how it is possible to efficiently check whether a sequence is \(\eta\)-valid. We will not need to do that since during the execution of our algorithm a convenient \(\mathcal{H}'\) will be always available when needed.

We remark that LogAdaptiveAlgorithm generates \(0\)-valid sequences in each phase. Specifically, it generates sequences \(\mathcal{H}= (H_1, H_2, \ldots, H_L)\), where each \(H_p\) is \(0\)-valid (with respect to \(\theta = (1+\varepsilon^2)^{p-1}\) and \((H_1,\ldots,H_{p-1})\)) and contains only regular facilities. The set \(S\) of opened facilities is such that \(A = \emptyset\) at the final time \(\theta\).

For a general sequence of sets of facilities \(\mathcal{H}= (H_1, H_2, \ldots, H_L)\), let us say it consists of \(\eta\)-valid sequences if, for every \(p \in [L]\), \(H_p\) is an \(\eta\)-valid sequence with respect to \(\theta = (1+\varepsilon^2)^{p - 1}\) and \((H_1,\ldots,H_{p-1})\). Furthermore, we say that \(\mathcal{H}\) is a solution if the set \(S=\cup_{q\leq p}H_q\) of opened facilities ensures that \(A = \emptyset\) at the end, meaning that at the time \(\theta = (1+\varepsilon^2)^L\) of the \((L+1)\)st phase, the set of clients \(\{ j \mid \theta \geq (1-\delta) d^2(j,S) \}\) equals all clients in \(D\). Otherwise, we say that \(\mathcal{H}\) is a partial solution.

The results in 5.2, namely 7 and 8, can thus be stated as if \(\mathcal{H}\) is a solution of \(0\)-valid sequences consisting of regular facilities then the set \(S\) of opened facilities satisfies \[\sum_{j\in D} d^2(j, S) \leq \frac{\Gamma}{1-\delta} \cdot \left( \text{opt}_{\text{LP}}(f) - f \cdot |S| \right)\,.\] We now allow to open facilities that are paid up to an amount \(\hat{f}- \eta\gamma n\) instead of \(\hat{f}\), and for free facilities that can be opened without paying them. This leads to the following generalization.

In the next subsection we will merge two solutions into one with \(|S_{\mathsf{reg}}| =k\) and \(|S_{\mathsf{free}}| =O(\log n/\varepsilon^2)\). To our help, we will have two basic routines \(\mathrm{\small CompleteSolution}\) and \(\mathrm{\small CompleteSequence}\):

• \(\mathrm{\small CompleteSolution}_{f,u}(\mathcal{H})\) takes a partial solution \(\mathcal{H}\) as input and returns a solution by running LogAdaptiveAlgorithmfrom the while-loop of stage \(1\), starting with the open facilities in \(\mathcal{H}\) and using the opening cost \(f\) and the distances \(u\) for free facilities already opened.

• \(\mathrm{\small CompleteSequence}_{f,u}(\mathcal{H}= (H_1, \ldots, H_{p-1}), H_p = \langle h_1, \ldots, h_\ell\rangle )\) where \(\mathcal{H}\) is partial solution and \(H_p\) is an \(\eta\)-valid sequence except that it may not be maximal. \(\mathrm{\small CompleteSequence}\) then completes \(H_p\) into a maximal sequence by adding the open facilities obtained by running the while-loop of stage \(1\) in LogAdaptiveAlgorithm starting with \(H_p\) (with respect to opening cost \(f\) and distances \(u\) for the free facilities already opened).

We note that all facilities added by CompleteSolution and CompleteSequence are \(0\)-openable and regular. Additionally, for intuition, observe that, with the above notation, the set of facilities opened by \(\mathrm{\small CompleteSolution}_{f,u}(\emptyset)\) is equivalent to that opened by LogAdaptiveAlgorithm. (There is no free facility so \(u\) is meaningless here.) To simplify notation, we only use \(\mathrm{\small CompleteSequence}\) when the partial solution \(\mathcal{H}\) is clear from the context, so we simply write \(\mathrm{\small CompleteSequence}_{(f,u)}(\langle h_1, \ldots, h_\ell \rangle)\). Furthermore, we omit \((f,u)\) when the parameters are clear from the context.

5.4 MergeSolutions↩︎

We now present our procedure MergeSolutions that walks between two solutions. We will maintain two solutions \((\mathcal{H}, f, u)\) and \((\mathcal{H}', f', u')\) whose numbers of open regular facilities sandwich \(k\). (Formally, let us say \(a\) and \(b\) sandwich \(k\) if \(a<k<b\) or \(b<k<a\).) We will gradually make them closer until one of them opens exactly \(k\) regular facilities. We start from the two initial solutions obtained using the standard method. Note that the minimum nonzero distance is \(1\) and the maximum squared distance is \(M= \poly(n)\), which implies that \(\text{opt}\leq nM\).

5.4.0.1 Initialization.

Binary search on \(f\) in the interval \([1/n^2, 4nM]\) and run \(\mathrm{\small CompleteSolution}_{f, u}(\emptyset)\) (with empty \(u\)). When \(f = 1/n^2\), every client \(j \in D\) satisfies \(d(j, S) = 0\) if \(d(j, F) = 0\) and \(d^2(j, S) \leq d^2(j, F) + 2/n^2\) otherwise, as each client can open its closest facility by itself when \(\alpha _j = d^2(j, F) + 2/n^2\). If this solution opens at most \(k\) centers, then it is already an \((1+2/n)\)-approximation. On the other hand, when \(f = 4nM\), the algorithm will open exactly one facility, since the first facility is open when \(\theta \geq 2M\), so every client becomes immediately inactive when it is open. Therefore, the binary search yields two solutions \((\mathcal{H}, f, u)\) and \((\mathcal{H}', f', u)\) where \((\mathcal{H}, f, u)\) opens less than \(k\) facilities and \((\mathcal{H}', f', u)\) opens more than \(k\) facilities where \({f'} \leq {f} \leq {f'} + \eta\).

5.4.0.2 General Step.

MergeSolutions works on phase \(p = 1, 2, \dots\) iteratively, and ensures that the two solutions are identical up to the \(p\)-th phase. Formally, MergeSolutions takes two solutions \((\mathcal{H}, f, u)\) and \((\mathcal{H}', f', u')\) as well as the current phase \(p\) as input. The following promises on them will be satisfied whenever we call MergeSolutions.

• Both consist of \(\eta\)-valid sequences.

• The numbers of open regular facilities by the two solutions sandwich \(k\).

• Exactly one parameter, which we call the difference parameter, differs by at most \(\eta\) between the two solutions, i.e, if \(f\neq f'\) then \(|f-f'| \leq \eta\) and \(u(\tilde{i}) = u'(\tilde{i})\) for every \(\tilde{i} \in F_{\mathsf{free}}\); otherwise if \(f=f'\) there is exactly one free facility \(\tilde{i}\in F_{\mathsf{free}}\) so that \(u(\tilde{i}) \neq u'(\tilde{i})\) and \(|u(\tilde{i}) - u'(\tilde{i})| \leq \eta\).

• \(\mathcal{H}= (H_1, H_2, \ldots, H_L)\) and \(\mathcal{H}' = (H'_1, H'_2, \ldots, H'_{L'})\) have a common prefix of length \(p-1\), i.e., \(H_i = H'_i\) for \(i=1,2, \ldots, p-1\). Furthermore, the remaining sequences \(H_p, \ldots, H_L\) and \(H'_p,\ldots, H'_{L'}\) only contain regular facilities and are \(0\)-valid with respect to the parameters \(f, u\) and \(f', u'\), respectively. (This implies that \(\mathcal{H}\) and \(\mathcal{H}'\) open the same set of free facilities.)

Given the promises, the goal of MergeSolutions is to

• either find a solution of \(\eta\)-valid sequences that opens exactly \(k\) regular facilities and \(O(\log n{/\varepsilon^2})\) free facilities, or

• produce two new solutions that satisfy the above promises for phase \(p + 1\).

Since the maximum squared distance between any two points is \(M\leq \poly(n)\) and \(f \leq 4nM\), before \(\theta\) reaches \(6M\), at least one facility will open and every client will become inactive. Therefore, for some threshold \(p^* = O(\log n / \varepsilon^2)\), once \(\mathcal{H}\) and \(\mathcal{H}'\) agree on the first \(p^*\) phases, they open exactly the same set of facilities which violates the promise that the numbers of open regular facilities sandwich \(k\). Therefore, case (1) must happen for some \(p \leq p^*\), which yields a solution of \(\eta\)-valid sequences that opens exactly \(k\) regular facilities and \(O(\log n{/\varepsilon^2})\) free facilities. Combined with 4, this proves 2. Indeed, we then have a solution consisting of open centers \(S\) with \(|S| = k + O(\log(n)/\varepsilon^2)\) and \(k\) regular facilities. Then, by Theorem 4, \[\begin{align} \sum_{j\in D} d^2(j, S) &\leq \frac{\Gamma}{1-\delta} \cdot \left( \text{opt}_{\text{LP}}(f) - (f-\eta n) \cdot k \right) \leq (\Gamma+O(\varepsilon)) \text{opt}\,, \end{align}\] where, for the last inequality, we used that \(\eta = 2^{-n}\) and \(\text{opt}\geq 1\) (since any distance is at least \(1\)). 2 then follows by scaling \(\varepsilon\) appropriately.

In the remainder of this section, we present MergeSolutions for two solutions and a given phase \(p\). It consists of three subroutines, presented in Section 5.4.1 5.4.2, and 5.4.3 respectively. Each subroutine achieves (1) or (2) above (which ends the whole procedure or lets us move on to the next phase), or sets up the stage for the next subroutine. (The final third subroutine always achieves (1) or (2).)

5.4.1 Making Parameters Identical↩︎

By renaming, we assume that if \(f\neq f'\) is the difference parameter then \(f' < f\), and if the difference parameter is \(u(\tilde{i}) \neq u'(\tilde{i})\) then \(u(\tilde{i}) < u'(\tilde{i})\).

5.4.1.1 Outcome of this section of MergeSolutions.

Let \[\mathcal{H}'' = \mathrm{\small CompleteSolution}_{(f,u)}(H'_1, \ldots, H'_p)\,,\] which is a valid call by the above claim. If \(\mathcal{H}''\) opens \(k\) regular facilities, we return this solution. Otherwise, as \((\mathcal{H}, f,u)\) and \((\mathcal{H}', f', u')\) sandwich \(k\), we must have that either \((\mathcal{H}, f, u)\) and \((\mathcal{H}'', f, u)\), or \((\mathcal{H}'', f, u)\) and \((\mathcal{H}', f', u')\) sandwich \(k\). On the one hand, if \((\mathcal{H}'', f, u)\) and \((\mathcal{H}', f',u')\) sandwich \(k\), then we made progress in that we now have two solutions with a common prefix of length \(p\) (instead of \(p-1\)) and only one parameter differs by at most \(\eta\). In this case, we are done with phase \(p\) and can continue MergeSolutions with these two solutions and phase \(p+1\).

On the other hand, if \((\mathcal{H}, f,u)\) and \((\mathcal{H}'', f, u)\) sandwich \(k\), then these two solutions have the same parameters and their first \(\eta\)-valid sequences that differ are \(H_p\) and \(H'_p\). In this case, we rename \(\mathcal{H}''\) to \(\mathcal{H}'\) and proceed to the next step.

5.4.2 Maximize Common Prefix of \(H_p\) and \(H'_p\)↩︎

We are given two solutions \((\mathcal{H}, f,u)\) and \((\mathcal{H}', f, u)\) on the same parameters \(f,u\) that sandwich \(k\). Moreover, if we let \(\mathcal{H}= (H_1, \ldots, H_L)\) and \(\mathcal{H}' = (H'_1, \ldots, H'_{L'})\), then \(H_i = H'_i\) for \(i= 1, \ldots, p-1\) and \(H_p \neq H'_p\). In addition \(H_p\) and \(H'_p\) only contain regular facilities.

In this step of MergeSolutions, we increase the common prefix of \(H_p\) and \(H'_p\) by “slowly” modifying \(H'_p\). As \(H_p\) only contains regular facilities, we let \(H_p = \langle i_1, i_2, \ldots, i_\ell \rangle\). Moreover, as the parameters \((f,u)\) are clear from the context, we omit them in the following when, e.g., calling the procedures CompleteSolution and CompleteSequence.

5.4.2.1 Transforming \(H'_p\) step-by-step so that \(H_p\) becomes a prefix.

Repeat the following until all of \(H_p\) becomes a prefix of \(H'_p\).

• Let \(q\) be the largest index so that \(H'_p = \langle i_1, i_2, \ldots, i_q, h'_{q+1}, \ldots, h'_{\ell_1}\rangle\), i.e., \(q\) equals the length of the common prefix \(i_1, i_2, \ldots, i_q\) of \(H_p\) and \(H'_p\).

• Update \(H'_p\) by inserting a free copy of \(i_{q+1}\), \(\tilde{i}_{q+1}\), at the end of \(H'_p\) with \(u(\tilde{i}_{q+1}) = 0\). That is, we get the sequence \[\begin{gather} H'_p = \langle i_1, i_2, \ldots, i_q, h'_{q+1}, \ldots, h'_{\ell_1}, \tilde{i}_{q+1} \rangle,. \end{gather}\]

• For \(t = q+1 , \ldots, \ell_1\),

  • Update \(H'_p\) by first removing \(h'_t\) and then completing the sequence by a call to CompleteSequence. More explicitly, if \[\begin{align} H'_p = &\langle i_1, \ldots, i_q, h'_t, h'_{t+1} \ldots, h'_{\ell_1}, \tilde{i}_{q+1}, h'_{\ell_1+1}, \ldots, h'_{\ell_2}\rangle\\ \intertext{then we delete h'_t to obtain} H''_p = &\langle i_1, \ldots, i_q, h'_{t+1} \ldots, h'_{\ell_1}, \tilde{i}_{q+1}, h'_{\ell_1+1}, \ldots, h'_{\ell_2}\rangle\\ \intertext{and finally update H'_p = \mathrm{\small CompleteSequence}(H''_p) which will be of the form} &\langle i_1, \ldots, i_q, h'_{t+1}, \ldots, h'_{\ell_1}, \tilde{i}_{q+1}, h'_{\ell_1+1}, \ldots, h'_{\ell_2}, \ldots, h'_{\ell_3}\rangle\,, \end{align}\] where \(h'_{\ell_2+1},\ldots, h'_{\ell_3}\) are the facilities added by CompleteSequence(which could be empty).

• After removing \(h'_{q+1}, \ldots, h'_{\ell_1}\), we have \[\begin{align} H'_p &= \langle i_1, i_2, \ldots, i_q, \tilde{i}_{q+1}, h'_{x}, \ldots, h'_{y}\rangle\,.\\ \intertext{We update H'_p by replacing \tilde{i}_{q+1} with its regular copy i_{q+1}, i.e., we obtain} H'_p &= \langle i_1, i_2, \ldots, i_q, {i}_{q+1}, h'_{x}, \ldots, h'_{y}\rangle\,. \end{align}\]

We have the following claims.

5.4.2.2 Outcome of this section of MergeSolutions.

We distinguish three cases.

5.4.2.3 Case 1.

If at some point \(H'_p\) is such that \(\mathrm{\small CompleteSolution}(H'_1, H'_2, \ldots, H'_p)\) opens exactly \(k\) facilities, then we return that solution.

5.4.2.4 Case 2.

Otherwise, suppose there is an update to \(H'_p\) so that if we let \(H^{\text{before}}_p\) and \(H^{\text{after}}_p\) denote \(H'_p\) before and after this update, respectively, then CompleteSolution\((H_1, H_2, \ldots, H^{\text{before}}_p)\) and CompleteSolution\((H_1, H_2, \ldots, H^{\text{after}}_p)\) sandwich \(k\). There are three kinds of updates to \(H'_p\). However, note that the update where \(\tilde{i}_{q+1}\) is replaced by its regular copy \(i_{q+1}\) cannot satisfy the conditions of this case. Indeed, as \(u(\tilde{i}_{q+1}) = 0\), we have \[\mathrm{\small CompleteSolution}(H_1, \ldots, H^{\text{before}}_p) = \mathrm{\small CompleteSolution}(H_1, \ldots, H^{\text{after}}_p)\] in this case, which contradicts that they sandwich \(k\).

We proceed to analyze the two other kind of updates. If the update to \(H'_p\) is to insert a free facility \(\tilde{i}_{q+1}\) with \(u(\tilde{i}_{q+1}) = 0\). Then \[\begin{align} H^{\text{before}}_p&= \langle i_1, i_2, \ldots, i_q, h'_{q+1}, \ldots, h'_{\ell_1} \rangle\\ H^{\text{after}}_p&= \langle i_1, i_2, \ldots, i_q, h'_{q+1}, \ldots, h'_{\ell_1}, \tilde{i}_{q+1} \rangle \end{align}\] Note that the solution CompleteSolution\((H_1, H_2, \ldots, H^{\text{before}}_p)\) is equivalent to CompleteSolution\((H_1, H_2, \ldots, H^{\text{after}}_p)\) if we set \(u(\tilde{i}_{q+1}) = 10M\) (except for the free facility \(\tilde{i}_{q+1}\)). We can thus do binary search on \(u(\tilde{i}_{q+1})\) to find two values \(u'(\tilde{i}_{q+1})\) and \(u(\tilde{i}_{q+1})\) with \(|u'(\tilde{i}_{q+1}) - u(\tilde{i}_{q+1})| \leq \eta\) so that the solutions \(\mathrm{\small CompleteSolution}_{(f,u')}(H_1, \ldots, H^{\text{after}}_p)\) and \(\mathrm{\small CompleteSolution}_{(f,u)}(H_1, \ldots, H^{\text{after}}_p)\) sandwich \(k\) (where \(u\) and \(u'\) are identical except for \(\tilde{i}_{q+1}\)). We have thus finished phase \(p\), and we can restart MergeSolutions with the two solutions and phase \(p+1\), since they have the first \(p\), instead of \(p-1\), \(\eta\)-valid sequences in common and exactly one difference parameter. Their maximality follows from that of \(H^{\text{before}}_p\).

Finally, consider when the update to \(H'_p\) is of the second type, i.e., \(u(\tilde{i}_{q+1}) = 0\) and \[\begin{align} H^{\text{before}}_p= &\langle i_1, \ldots, i_q, h'_t, h'_{t+1} \ldots, h'_{\ell_1}, \tilde{i}_{q+1}, h'_{\ell_1+1}, \ldots, h'_{\ell_2}\rangle\\ H^{\text{after}}_p= &\langle i_1, \ldots, i_q, h'_{t+1}, \ldots, h'_{\ell_1}, \tilde{i}_{q+1}, h'_{\ell_1+1}, \ldots, h'_{\ell_2}, \ldots, h'_{\ell_3}\rangle \end{align}\] In this case, \(\mathrm{\small CompleteSolution}(H_1, \ldots, H^{\text{before}}_p)\) and \(\mathrm{\small CompleteSolution}(H_1, \ldots, H^{\text{after}}_p)\) are two solutions with the same parameters \(f, u\), that sandwich \(k\) and we proceed to the next step with the following properties of \(H^{\text{before}}_p\) and \(H^{\text{after}}_p\):

• Both sequences contain one free facility, \(\tilde{i}_{q+1}\).

• \(H^{\text{before}}_p\) and \(H^{\text{after}}_p\) are identical except that (1) \(H^{\text{before}}_p\) contains a regular facility \(h'_t\) not present in \(H^{\text{after}}_p\), and (2) \(H^{\text{after}}_p\) contains a (potentially empty) suffix of regular facilities not present in \(H^{\text{before}}_p\).

5.4.2.5 Case 3.

None of the previous cases apply. Then, completing the solution with the initial \(H'_p\) and the final \(H'_p\) cannot sandwich \(k\). Since then (assuming the first case does not apply) there must be an update to \(H'_p\) that results in two solutions that sandwich \(k\), since the the completion of \(H_p\) and the completion of the initial \(H'_p\) sandwich \(k\). The above implies that if we let \(H^{\text{before}}_p= H_p\) and \(H^{\text{after}}_p\) equal the final \(H'_p\) then CompleteSolution\((H_1, \ldots, H^{\text{before}}_p)\) and CompleteSolution\((H_1, \ldots, H^{\text{after}}_p)\) sandwich \(k\). Further note that since \(H^{\text{after}}_p\) equals \(H'_p\) after the final updates \(H^{\text{before}}_p\) appears as a prefix of \(H^{\text{after}}_p\) and both sequences only contain regular facilities. To summarize, in this case, CompleteSolution\((H_1, \ldots, H^{\text{before}}_p)\) and \(\mathrm{\small CompleteSolution}(H_1, \ldots, H^{\text{after}}_p)\) are two solutions with the same parameters \(f, u\), that sandwich \(k\) and we proceed to the next step with the following properties of \(H^{\text{before}}_p\) and \(H^{\text{after}}_p\):

• Both sequences contain no free facilities.

• \(H^{\text{before}}_p\) and \(H^{\text{after}}_p\) are identical except that \(H^{\text{after}}_p\) contains a suffix of regular facilities not present in \(H^{\text{before}}_p\).

5.4.3 Removing Extra Facilities↩︎

The input to this step consists of two solutions \(\mathrm{\small CompleteSolution}(H_1, \ldots, H^{\text{before}}_p)\) and CompleteSolution\((H_1,\ldots, H^{\text{after}}_p)\) that have the same parameters \(f, u\) and sandwich \(k\). Moreover, \(H^{\text{before}}_p\) and \(H^{\text{after}}_p\) both contain at most one free facility and they are identical except for potential two differences:

• \(H^{\text{before}}_p\) may contain a facility not present in \(H^{\text{after}}_p\);

• \(H^{\text{after}}_p\) contains a suffix of regular facilities (which may be empty) that is not present in \(H^{\text{before}}_p\).

To simplify notation in this section, we let \(i_*\) denote the regular facility in \(H^{\text{before}}_p\) that is not present in \(H^{\text{after}}_p\) if it exists (so \(i_*\) corresponds to \(h'_t\) in the previous section). We further let \(h_1, \ldots, h_\ell\) denote the suffix of \(H^{\text{after}}_p\) not present in \(H^{\text{before}}_p\).

We now transform \(H^{\text{after}}_p\) into \(H^{\text{before}}_p\) one change at a time:

5.4.3.1 Generation of \(\eta\)-valid sequences from \(H^{\text{after}}_p\) to \(H^{\text{before}}_p\).

• If \(H^{\text{before}}_p\) contains the regular facility \(i_*\) not present in \(H^{\text{after}}_p\), then let \(\tilde{i}_*\) be a free copy of \(i_*\) with \(u(\tilde{i}_*) =0\) and obtain \(H^{0}_p\) from \(H^{\text{after}}_p\) by adding \(\tilde{i}_*\) at the end. Otherwise, we let \(H^0_p\) equal \(H^{\text{after}}_p\).

• For \(r= 1, \ldots, \ell\), obtain \(H^r_p\) from \(H^{r-1}_p\) by removing \(h_r\).

We remark that \(H^\ell_p\) is identical to \(H^{\text{before}}_p\) since we removed the whole suffix \(h_1, \ldots h_\ell\) present in \(H^{\text{after}}_p\), except for potentially the facility \(i_*\). That is the only difference in \(H^{\text{before}}_p\) and \(H^\ell_p\) may be the placement of \(\tilde{i}_*\) and that \(\tilde{i}_*\) is a free copy of the regular facility \(i_*\) present in \(H^{\text{before}}_p\). However, since \(u(\tilde{i}_*) = 0\), we have that the two sequences \(H^{\text{before}}_p\) and \(H^\ell_p\) have an identical impact on CompleteSolution, which is a fact we use in the proof of Claim [claim:3rdstep95outcome].

5.4.3.2 Outcome of this section of MergeSolutions.

We distinguish three cases based on the first bullet point that is true in the above claim.

5.4.3.3 Case 1.

If one of the sequences \(H^r_p\) is such that \(\mathrm{\small CompleteSolution}(H_1, H_2, \ldots, H^r_p)\) opens \(k\) regular facilities, then we return that solution.

5.4.3.4 Case 2.

Otherwise, if the solutions \(\mathrm{\small CompleteSolution}(H_1, H_2, \ldots, H^{\text{after}}_p)\) and CompleteSolution\((H_1, H_2, \ldots, H^{0}_p)\) sandwich \(k\), then \(H^0_p\) is equal to \(H^{\text{after}}_p\) except that \(\tilde{i}_p\) was added in the end of \(H^{\text{after}}_p\) with \(u(\tilde{i}_p) = 0\). Note that the solution CompleteSolution\((H_1, H_2, \ldots, H^{\text{after}}_p)\) is equivalent to CompleteSolution\((H_1, H_2, \ldots, H^0_p)\) if we set \(u(\tilde{i}_{q+1}) = 10M\) (except the presence of \(\tilde{i}_*\)). We can thus do binary search on \(u(\tilde{i}_*)\) to find two values \(u'(\tilde{i}_*)\) and \(u(\tilde{i}_*)\) with \(|u'(\tilde{i}_*) - u(\tilde{i}_*)| \leq \eta\) so that the two solutions \(\mathrm{\small CompleteSolution}_{(f,u')}(H_1, \ldots, H^0_p)\) and \(\mathrm{\small CompleteSolution}_{(f,u)}(H_1, \ldots, H^0_p)\) sandwich \(k\) (where \(u\) and \(u'\) are identical except for \(\tilde{i}_{q+1}\)). We have thus finished phase \(p\), and we can restart MergeSolutions with the two solutions and phase \(p+1\), since they have the first \(p\), instead of \(p-1\), \(\eta\)-valid sequences in common and exactly one difference parameter.

5.4.3.5 Case 3.

Finally, consider the case when none of the above cases apply and (thus by the above claim) there is an \(r\in \{0, 1, \ldots, \ell-1\}\) so that CompleteSolution\((H_1, H_2, \ldots, H^r_p)\) and CompleteSolution\((H_1, H_2, \ldots, H^{r+1}_p)\) sandwich \(k\). The only difference between \(H^r_p\) and \(H^{r+1}_p\) is that \(H^r_p\) contains \(h_{r+1}\) which is not present in \(H^{r+1}\). Let \(H_p\) be the sequence obtained by adding a free copy \(\tilde{i}_{r+1}\) of \(h_{r+1}\) at the end of \(H^{r+1}\) (which clearly maintains that is \(\eta\)-valid), then \[\begin{align} \mathrm{\small CompleteSolution}(H_1, H_2, \ldots, H_p) &= \mathrm{\small CompleteSolution}(H_1, H_2, \ldots, H^r_p)if u(\tilde{i}_{r+1}) = 0\\ \intertext{and} \mathrm{\small CompleteSolution}(H_1, H_2, \ldots, H_p) &= \mathrm{\small CompleteSolution}(H_1, H_2, \ldots, H^{r+1}_p)if u(\tilde{i}_{r+1}) = 10M. \end{align}\] Hence, doing a binary search on \(u(\tilde{i}_{r+1})\) as in the previous case finishes phase \(p\). We can restart MergeSolutions with the two solutions and phase \(p+1\), since they have the first \(p\), instead of \(p-1\), \(\eta\)-valid sequences in common and exactly one difference parameter.

5.5 Analysis of the Robust Algorithm↩︎

In this section, we prove 4, showing that a solution \(\mathcal{H}\) of \(\eta\)-valid sequences yields an almost \(\Gamma\)-LMP approximation. The proof is almost identical to that for 3. There are only two minor changes:

1. 2 of \(\eta\)-openability allows for a facility to be paid up to an \(n\gamma\eta\)-additive difference. This difference only appears in 9 below, which is an analog of 7.

2. 3 of an \(\eta\)-valid sequence allows \(h_t\) is \(\eta\)-openable with respect to a super sequence \(\mathcal{H}'\supseteq \mathcal{H}\oplus(\langle h_1,\dots,h_{t-1}\rangle)\). This difference only appears in (the last part of) the proof of 10, an analog of 8.

These two lemmas immediately imply 4.

In order to use mostly the same analysis, one syntactic difference between the setting of 5.2 and the current setting that we have to reconcile is the the former analyzes the outcome of 3 while the latter analyzes a solution \(\mathcal{H}=(H_1,\dots,H_L)\) of \(\eta\)-valid sequences. To use the same terminology from 5.2, let us consider the execution of the solution \(\mathcal{H}\) where for \(p=1,\dots,L\), we run the \(p\)-th phase according to the sequence \(H_p\); in Stage 1, each facility in the sequence becomes open one-by-one with the corresponding \((\tau_j)_j\) values, and in Stage 2 at the end of the phase, \(\theta\leftarrow (1+\varepsilon^2)\theta\).

Then, an atomic step can still be defined in the same way as in 3; it corresponds to opening a facility in Stage 1, possibly a free facility, (and updating \(\alpha\), \(S\), \(A\), \(IC\), \(DC\) accordingly) or increasing \(\alpha\)-values in Stage 2 simultaneously (and updating \(IC\) accordingly). By “at any point in the algorithm”, we mean any point in the algorithm’s execution that is not in the middle of an atomic step.

Throughout the analysis, we let \(\alpha_j^*\) be the final \(\alpha\)-value of client \(j\in D\), and we let \(S^*\) be the final set of opened facilities. We start by analyzing the approximation guarantee of the algorithm with respect to \(\alpha^*\), and then we prove that \(\alpha^*/\Gamma\) is a feasible solution to the dual.

5.5.1 Approximation Guarantee↩︎

The following lemma is the analogue of Lemma 7.

5.5.2 Dual Feasibility↩︎

We next prove the following lemma, which is analogous to 8.

The following three claims and their proofs are identical to [clm:nooverbidding], [clm:notfrozen] and [clm:frozenstay], without any modification and using a similar notation.

As usual, fix an arbitrary facility \(i\in F\). The setup in this paragraph is also identical to 5.2.2. Let \(D^*=\{j\in D:\alpha_j^*>\Gamma d^2(i, j)\}\) be those clients that contribute to the left-hand side. We also reuse the definitions of freezing and \(IC^k, DC^k, A^k\) for for \(k\in D^*\) that becomes inactive strictly before \(i\) becomes frozen.

The remaining part of the proof needs some changes from 5.2.2, due to the fact that a facility may be openable with respect to a super-sequence instead of a sequence that the algorithm executes. Let \(i_0\) be the facility that freezes \(i\), which opens in (stage 1) of phase \(p\). Let \(\theta = (1+\varepsilon^2)^p\) be the \(\theta\) value in phase \(p\). Let \(\mathcal{H'} \supseteq \mathcal{H}\) be the super-sequence with respect to which \(i_0\) is openable. As in 5.2, we order the clients in \(D^*\) according to the time they are removed from \(A\) (in the algorithm) and break ties according to \(\alpha^*\)-values in increasing order. The only additional rule is for the clients who are in \(A\) right before \(i_0\) is open; let us call them \(A\). They are, by [clm:robustfrozenstay], the last group of clients in \(D^*\). Let \(A' \subseteq A\) be the clients who are active right before \(i_0\) is open according to the supersequence \(\mathcal{H'}\). Then for the order within \(A\), we place \(A \setminus A'\) first and \(A'\) after, and within each group we sort by the \(\alpha^*_j\) values. Similarly, let \(IC'\) and \(DC'\) be the set of indirectly and directly connected clients right before \(i_0\) is open in the supersequence \(\mathcal{H}'\), and \(I'=A\setminus A'\). For \(k\in D^*\) that becomes inactive when \(i\) becomes frozen, we define \(DC^k\) by \(DC'\cap [k-1]\), and \(IC^k\) by \((IC'\cup A')\cap [k-1]\). We prove the following analog of [fct:bidbound].

We are now ready to prove 10.

We conclude with the proof of 4.

6 \((5+O(\sqrt{\varepsilon}))\)-Approximation for \((\frac{\zeta}{\log n})\)-Stable Instances↩︎

In this section, we give a \((5+O(\sqrt{\varepsilon}))\)-approximation algorithm for \((\frac{\zeta}{\log n})\)-stable \(k\)-means instances.

To simplify notation, we let \(\varepsilon\) be the minimum of \(\varepsilon\) and \(\zeta\) in the above theorem and further assume that \(\varepsilon<1/12\) is sufficiently small. With this notation, we present a \((5+O(\sqrt{\varepsilon})\)-approximation algorithm for \(\varepsilon/\log(n)\)-stable instances. This implies the above theorem as any \(\zeta/\log(n)\)-stable instance is also \(\zeta'/\log(n)\)-stable with \(\zeta' \leq \zeta\).

The algorithm and the proof of the above theorem follow several steps where we iteratively simplify the instance by adding more structure. Finally, we obtain enough structure to solve the problem by maximizing a submodular function subject to a partition matroid constraint. We have not optimized constants in favor of simplifying the description, and throughout this section we let \(O_\varepsilon(\log n))\) denote \(O(\log(n)/\varepsilon^{O(1)})\). An overview of this section and the algorithm is as follows:

6.0.0.1 Overview of a \((5+O(\sqrt{\varepsilon})\)-Approximation for \(k\)-means on \((\varepsilon/\log n)\)- Stable Instances \(k, D, F, d\).

1. We start, in 6.1, by running a standard Local Search on \((k, D, F, d)\) to obtain a set of centers \(S\). It is well-known that \(S\) is a constant-factor approximation. Moreover, for a stable instance, we show that \(S\) “correctly” implicitly identifies all but \(O_\varepsilon(\log n)\) centers of \(\text{OPT}\). Specifically, this ensures that all but \(O_\varepsilon(\log n)\) clusters have the same cost as in \(\text{OPT}\) up to a \((1+O(\sqrt{\varepsilon}))\) factor.

2. In 6.3, we use a technique inspired from the classic \(k\)-means++ algorithm and its so called \(D^2\)-sampling procedure to sample clients proportional to their cost in solution \(S\), and we show that if we take \(O_\varepsilon(\log n)\) samples \(W\) then with good probability we “hit” all but an insignificant fraction of clusters of high cost in \(S\).

3. In 6.4, we then, guess a subset of the sampled clients \(W\) so that each hit cluster has exactly one client, which we refer to as the leader of that cluster. In addition, for each client \(\ell\) that is a leader, we guess the (approximate) distance from \(\ell\) to the optimal center in the optimal cluster it hits. This defines a set \(\mathcal{B}\) of balls with the guarantee that the optimal solution opens one center in each of these balls. Moreover, if we let \(\hat{\text{OPT}}\) be the optimal centers corresponding to these balls, we prove in 6.5 that there is a swap of centers \(M_O = S-S_0 \cup \hat{\text{OPT}}\) that removes \(S_0\) and adds the optimal centers \(\hat{\text{OPT}}\) so that the cost of the solution \(M_O\) is good and in particular a \((4+O(\sqrt{\varepsilon}))\)-approximation.

4. The remaining steps are then focused on approximating this swap, i.e., intuitively to that of finding \(S_0\) and \(\hat{\text{OPT}}\):

   1. We first find the set \(\mathcal{Q}\) of expensive centers of \(S_0\) so that the clusters of centers \(S_0-\mathcal{Q}\) have total cost \(O(\varepsilon\text{opt})\) in 6.6. This allows us to think of the clusters of centers in \(S_0 - \mathcal{Q}\) as contracted, i.e., as single weighted points, and in particular, that all clients of these clusters are assigned to the same center when their center in \(S_0\) is removed. We then, in 6.7, use this structure to partition \(S_0 - \mathcal{Q}\) into two sets \(\mathcal{U}\) and \(\mathcal{R}\) where \(\mathcal{U}\) are those centers that are reassigned to other centers in \(S - S_0\) and \(\mathcal{R}\) are those centers that are assigned to \(\hat{\text{OPT}}\).

   2. The set \(\mathcal{U}\) is then approximated with a set \(\widetilde{\mathcal{U}}\) (with close to identical properties) in 6.8.

   3. At this stage, we have thus “guessed” all centers \(\mathcal{Q}\cup \widetilde{\mathcal{U}}\) to remove from \(S\) except for those in \(\mathcal{R}\). This is achieved in 6.10 where we reduce the problem of both finding \(\mathcal{R}\) and an approximate version of \(\hat{\text{OPT}}\) to that of submodular function optimization subject to a partition matroid (since we should open one center per ball in \(\mathcal{B}\)).

We remark that in each of Steps 3, 4a, 4b, and 4c, we will make guesses by enumerating polynomially \(n^{\varepsilon^{-O(1)}}\) many potential choices. In the following, we describe the algorithm and its analysis by presenting each of these steps. When presenting the next step, we analyze the branch of the algorithm that took the successful guesses up to that point. The algorithm outputs the solution of the smallest cost among the \(n^{\varepsilon^{-O(1)}}\) constructed solutions obtained by enumerating all possible guesses (see 6.11 for a more formal treatment where we put everything together). In particular, the final solution will have cost at most the solution that is output in the analyzed branch of “correct” guesses.

We end this section by introducing basic definitions and concepts, followed by the definition of matroids and submodular functions. We also state the known result, an \((1-1/e)\)-approximation algorithm for maximizing monotone submodular functions subject to a matroid constraint, that we use in the last step (6.10).

6.0.0.2 Definitions of leaders.

We use the following definitions throughout this section. Given a set of centers \(A\subseteq F\), we let \(\text{cost}(A):=\sum_{p\in D}d^2(p,A)\) denote the cost of the corresponding clustering (where each client \(p\) is assigned to the closest center in \(A\) and pays the corresponding distance). Sometimes we will consider a (possibly suboptimal) assignment \(\mu:D\rightarrow A\) of clients to centers in \(A\), and let \(\text{cost}(A,\mu):=\sum_{p\in D}d^2(p,\mu(p))\) be the corresponding cost. Let \(\text{OPT}\) be a fixed optimum solution for the \(k\)-means instance, and let \(\text{opt}= \text{cost}(\text{OPT})\) be its cost. For a client \(p\in D\), let \(\text{OPT}_p:=d^2(p,\text{OPT})\) be the cost paid by \(p\) in \(\text{OPT}\).

Let \(C^*\) be an \(\text{OPT}\) cluster with center \(c^*\in \text{OPT}\). We let \(\text{avg}_{C^*, \text{OPT}}\) be the squared distance from \(c^*\) to the \(\varepsilon|C^*|\)th closest client of \(C^*\) to \(c^*\). We further let \(C^*_\text{avg}\) be the set of clients of \(C^*\) at a squared distance at most \(\text{avg}_{C^*, \text{OPT}}\) from \(c^*\), we refer to these clients as leaders for \(C^*\). So \(\text{avg}_{C^*, \text{OPT}}\) is the maximum squared distance from a leader to its center \(c^*\). Also note that we have, \(|C^*_{avg}|\geq {\varepsilon}|C^*|\). We thus have that \(\text{avg}_{C^*, \text{OPT}}\) is upper bounded by \(\frac{1}{1-\varepsilon}\) times the average cost \(\frac{1}{|C^*|} \sum_{p\in C^*} d^2(p, c^*)\), and at the same time \(C^*_{avg}\) contains a constant fraction of the clients.

6.0.0.3 Partition matroids.

In 6.10, we formulate a submodular function that we then maximize over a partition matroid using known results. We define those concepts and also state the result that we will use. First, recall the definition of a matroid and partition matroid. A matroid is a tuple \((E, \mathcal{I})\) defined on a ground set \(E\) with a family of independent sets \(\mathcal{I}\) that satisfy: (1) if \(A \in \mathcal{I}\) and \(A' \subseteq A\), then \(A' \in \mathcal{I}\); and (2) if \(A, B \in \mathcal{I}\) and \(|A| < |B|\), then there exists an element \(e \in B - A\) such that \(A \cup \{e\} \in \mathcal{I}\). In the special case of a partition matroid, the ground set \(E\) is partitioned into disjoint subsets \(E_1, E_2, \dots, E_k\), and there are non-negative integers \(r_1, r_2, \dots, r_k\) such that:

• A subset \(I \subseteq E\) is independent if and only if \(|I \cap E_i| \leq r_i\) for each \(i = 1, 2, \dots, k\).

In other words, each subset \(E_i\) has a capacity limit \(r_i\), and an independent set \(I\) cannot contain more than \(r_i\) elements from \(E_i\). To work with matroids (to get a running time that is polynomial in the size of the ground set \(E\)), we often use independence queries, which allow us to determine if a given subset \(A \subseteq E\) is independent (i.e., whether \(A \in \mathcal{I}\)). Specifically, an independence query on a subset \(A\) returns true if \(A\) is independent and false otherwise.

6.0.0.4 Submodular functions.

Having defined matroids, we proceed to define non-negative monotone submodular functions. For a finite ground set \(E\), a set function \(f: 2^E \to \mathbb{R}\) is submodular if it satisfies the diminishing returns property: for every pair of sets \(A \subseteq B \subseteq E\) and any element \(e \in E - B\), it holds that \[f(A \cup \{e\}) - f(A) \geq f(B \cup \{e\}) - f(B).\] Intuitively, this means that adding an element \(e\) to a smaller set \(A\) provides at least as much additional value as adding \(e\) to a larger set \(B\). In addition, \(f\) is called monotone if for any pair of sets \(A \subseteq B \subseteq E\), we have \(f(A) \leq f(B)\). In other words, adding elements to a set does not decrease the function value. Finally, \(f\) is non-negative if \(f(A) \geq 0\) for all \(A \subseteq E\).

A celebrated result [23] gives a polynomial-time \((1-1/e)\)-approximation algorithm for the problem of maximizing a non-negative monotone submodular function over a matroid constraint. While the original result was a randomized algorithm, a recent striking result [24], building upon the work of [25], gives a deterministic algorithm with the same guarantee. We remark that polynomial time here refers to a running time that is polynomial in the size of the ground set \(E\). To summarize, they show the following theorem.

6.1 Properties of a Locally Optimal Solution↩︎

Let \(S\) be a locally optimal solution output by the following standard local search algorithm. We remark that its running time is polynomial since by 1 the cost of a solution is a polynomially bounded integer and, at each improving step, the cost decreases by at least one.

6.1.0.1 Classic Local Search Procedure for \(k\)-Means.

1. \(S \gets\) Arbitrary solution for \(k\)-Means.

2. While there exists a solution \(S'\) such that \(|S \Delta S'|\le 2\) and \(\text{cost}(S') <\text{cost}(S)\):  

   1. \(S \gets S'\).

3. Output \(S\).

Local search achieves a \(25\) approximation for \(k\)-Means, see [4].

In the remaining part of the algorithm and its analysis we fix \(S\), and recall that we fixed the reference optimal solution \(\text{OPT}\). We distinguish different types of clusters of \(\text{OPT}\) and analyze their properties in the next two sections. We let \(S_p:=d^2(p,S)\) be the cost paid by \(p\) in \(S\) and recall that \(\text{OPT}_p\) denotes the cost \(p\in D\) pays in \(\text{OPT}\).

6.2 Pure Clusters↩︎

We say that a cluster \(C^*\) of \(\text{OPT}\) is pure (w.r.t. \(S\)) if there exists a cluster \(C'\) of \(S\) such that \(|C' \Delta C^*| \le {3\varepsilon} \min(|C^*|,|C'|)\). In which case we also say that \(C'\) is pure w.r.t. \(\text{OPT}\). We say that \(C'\) and \(C^*\) are associated. The following lemma follows from the fact that the input is \(\beta\)-stable for \(\beta := \varepsilon/\log n\) and that the solution \(S\) is a \(25\)-approximation by 11. Given a client or center \(a\) and a distance \(r\), we let \(B(a,r)\) be the set of clients and centers at distance at most \(r\) from \(a\). This lemma is originally proved in [26], we include a proof for completeness.

Let \(C^*\) be a pure cluster of \(\text{OPT}\) and \(C'\) be the associated cluster of \(S\). Let \(c^*\) be the center of \(C^*\) and \({c'}\) be the center of \({C'}\). We define \(t(c^*) := {c'}\). For any client \(p \in C^*\), let \(r(p) := d^2(p, {c'}) = d^2(p, t(c^*))\). We have the following lemma which is the only lemma that requires that the solution \(S\) is obtained via the local search (the proof of the previous lemma only used that it was a constant-factor approximation). In words, it says that \(S\) approximates the connection cost of clients belonging to pure clusters of the optimal solution almost perfectly.

6.3 Cheap Clusters and \(D\)-Sample Process to Hit Expensive Clusters↩︎

In the previous section, we showed that the cost of pure clusters in \(\text{OPT}\) is close to to their optimal cost in our local search solution \(S\). Our goal in this section is to “hit” those non-pure clusters of \(\text{OPT}\) that have a high cost in \(S\). We say that a non-pure cluster \(C^*\) of \(\text{OPT}\) with center \(c^*\) is basic-cheap if the total cost in \(S\) of the clients in \(C^*_\text{avg}\) is less than \(\varepsilon^{5}\text{opt}/\log n\) or if there exists a center of \(S\) at squared distance at most \((1+\varepsilon)\text{avg}_{C^*,\text{OPT}} + \varepsilon\gamma_{C^*}\) from \(c^*\), where \(\gamma_{C^*} := \frac{1}{|C^*|} \sum_{p\in C^*} (\text{OPT}_p + S_p)\); in the latter case we also say that \(C^*\) is covered by \(S\).

The following procedure, that aims to “hit” all clusters of \(\text{OPT}\) that are non-pure and non-basic-cheap, is inspired by the classic \(k\)-means++ algorithm and its so called \(D^2\)-sampling procedure to sample clients proportional to their cost in solution \(S\).

6.3.0.1 \(s\)-\(D^2\)-Sample Process.

1. Let \(W\) be the set obtained by taking \(s\) independent samples of clients where, in each sample, client \(p\) is sampled with probability \(\frac{d^2(p, S)}{\text{cost}(S)}\).

2. Output \(W\).

Our algorithm obtains a set \(W\) of clients by sampling \(s^*\) using the \(s^*\)-\(D^2\)-sample process, where \[s^* = 25\cdot \frac{\log n}{\varepsilon^{5}}\ln\left( \frac{{81}}{\varepsilon^{2}}\right)\,.\] We have the following lemma that says that we hit all non-pure and non-basic-cheap clusters except for a set \(V\) of clusters that have an insignificant cost in \(S\).

For a non-pure cluster \(C^*\) of \(\text{OPT}\) with center \(c^*\in \text{OPT}\), we let \(t(c^*)\) be the closest center in \(S\) to \(c^*\), and for any \(p\in C^*\), we let the replacement cost of \(p\) be \(r(p):=d^2(p,t(c^*))\).

6.3.0.2 Successful \(W\).

We say that the set of clients \(W\) obtained from the \(s^*\)-\(D^2\)-sample process is a successful sample if 21 holds. We define the following two types of \(\text{OPT}\) clusters. We say that a non-pure cluster \(C^*\) of \(\text{OPT}\) with center \(c^*\) is cheap if it is basic-cheap or is in \(V\); a non-pure cluster is expensive otherwise. In other words, a cluster \(C^*\) of the optimal solution is expensive if it is a non-pure and non-basic-cheap cluster such that \(C^*_{avg} \cap W \neq \emptyset\).

Assuming that \(W\) is successful, we have a good bound on the cost of the clients assigned to cheap centers in the optimal solution.

Given \(S\) and \(W\), the set of clusters of the optimal solution \(\text{OPT}\) is thus partitioned into pure clusters, cheap clusters, and expensive clusters. We let \(\hat{\text{OPT}}\) be the centers of \(\text{OPT}\) that are expensive. We further define the set \(D^*_{\texttt{exp}}\subseteq D\) to be the subset of clients that belong to the expensive clusters in \(\text{OPT}\). We have thus partitioned \(D\) into sets \(D^*_{\texttt{pure}}, D^*_{\texttt{cheap}}\), and \(D^*_{\texttt{exp}}\) depending on the type of optimal cluster they belong to. 13 bounds the cost of the clients in \(D^*_{\texttt{pure}}\) in \(S\) and 15 bounds the cost of the clients in \(D^*_{\texttt{cheap}}\). For future reference, we summarize these two lemmas (by weakening the upper bound for clients in \(D^*_{\texttt{pure}}\)).

It follows that the connection cost in solution \(S\) of clients in \(D- D^*_{\texttt{exp}}\) is roughly within a factor \({4}\) of their cost in the optimal solution. The remaining part of the algorithm is thus devoted to modifying \(S\) to obtain a small connection cost of the clients in \(D^*_{\texttt{exp}}\) without significantly increasing the cost of the other clients.

6.4 Identifying Balls of Expensive Clusters↩︎

By the definition of expensive clusters, each expensive cluster \(C^*\) of \(\text{OPT}\) satisfies \(C^*\cap W \neq \emptyset\) where \(W\) is the sample obtained by running the sampling procedure of the last subsection. Here, our goal is to guess a subset \(\mathcal{B}\) of balls so that for each expensive cluster \(C^*\) there is a ball \(B(\ell, {\sqrt{\rho}}) \in \mathcal{B}\) so that \(\ell\), which we refer to as a leader, is in \(C^*_{avg}\) and \(\rho\) is a very close approximation to \(\text{avg}_{C^*, \text{OPT}}\). This guarantees that the center \(c^* \in \hat{\text{OPT}}\) of cluster \(C^*\) is in \(B(\ell, {\sqrt{\rho}})\) and that no other center in \(B(\ell, {\sqrt{\rho}})\) is "too" far from \(c^*\) since the radius of the ball is \(\approx {\sqrt{\text{avg}_{C^*, \text{OPT}}}}\). The balls \(\mathcal{B}\) will then be used when we approximate the centers in \(\hat{\text{OPT}}\) via submodular function maximization subject to a partition matroid constraint (the balls will correspond to the partitions of the matroid). As the algorithm does not know the set of expensive clusters, it enumerates all subsets of \(W\), and for each guessed point, it enumerates a constant number of radii.

6.4.0.1 Ball Procedure.

1. \(\mathcal{L}_\texttt{bal} \gets \emptyset\).

2. For each subset \(\{p_1, p_2, \ldots, p_q\} \subseteq W\):

   1. For all integers \(i_1, i_2, \ldots, i_q\) such that \[\varepsilon^3 {S_{p_j}} \le (1+\varepsilon^3)^{i_j} \le {S_{p_j}}/\varepsilon^3 \qquad for j=1, 2,\ldots, q\] add \(\mathcal{B}= \left\{B(p_1, {\sqrt{(1+\varepsilon^3)^{i_1}}}), B(p_2, {\sqrt{(1+\varepsilon^3)^{i_2}}}), \ldots, B(p_q, {\sqrt{(1+\varepsilon^3)^{i_q}}}\right\}\) to \(\mathcal{L}_{\texttt{bal}}\).

We show that the ball procedure runs in polynomial time and that one of the subsets \(\mathcal{B}\) of balls that the procedure constructs is a correct guess: every expensive cluster \(C^*\) has a leader, and \(\text{avg}_{C^*, \text{OPT}}\) is guessed approximately correct. We say that a set of balls \(\mathcal{B}\) is a valid set of balls if it satisfies the following: For each expensive cluster \(C^*\) of \(\text{OPT}\), we have a ball \(B(\ell, {\sqrt{\rho}}) \in \mathcal{B}\) such that \(\ell\in C^*_{avg}\) is a leader and \[\text{avg}_{C^*,\text{OPT}} \leq \rho \leq \text{avg}_{C^*,\text{OPT}} + \varepsilon\gamma_{C^*},\] where we recall that \(\gamma_{C^*} = \frac{1}{|C^*|} \sum_{p \in C^*} (\text{OPT}_p + S_p)\).

We thus have that the ball procedure produces a family \(\mathcal{L}_{\texttt{bal}}\) such that at least one set \(\mathcal{B}\in \mathcal{L}_{\texttt{bal}}\) satisfies the conditions of the lemma. While our algorithm will try all possible sets in \(\mathcal{L}_{\texttt{bal}}\) (since it does not know which one is valid), we do our analysis by considering the run of the algorithm when it selects this set \(\mathcal{B}\) of valid balls.

6.5 Dummy Centers and Mixed Solutions \(M_O\) and \(M_D\)↩︎

Next, the algorithm creates, for each ball \(B(\ell, {\sqrt{\rho}}) \in \mathcal{B}\) with center \(\ell\) and radius \({\sqrt{\rho}}\), a dummy center \(\delta\). The center \(\delta\) is at distance \({\sqrt{\rho}}\) from \(\ell\) and and at distance \({\sqrt{\rho}} + d(\ell, p)\) from any other input point \(p\).

Let \(\Lambda\) be the set of dummy centers, so \(|\Lambda| = |\mathcal{B}|\). Recall that \(D^*_{\texttt{exp}}\) is the set of clients that are in an expensive cluster of \(\text{OPT}\), i.e., served by a center in \(\hat{\text{OPT}}\) in the optimal solution.

6.5.0.1 Definition of the mixed solutions \(M_O\) and \(M_D\).

Next, let us analyze the cost of the swap \((\hat{\text{OPT}}, S_0)\). Namely, of the solution \(M_O = S - S_0 \cup \hat{\text{OPT}}\), where \(M\) in \(M_O\) stands for “mixed solution” and the subscript \(O\) stands for that the mixed solution is obtained by swapping in some elements from \(\text{OPT}\) (and removing \(S_0\)). We will also analyze the “dummy” version of \(M_O\) where instead of swapping in \(\hat{\text{OPT}}\), we swap in the dummy centers \({\Lambda}\). We denote that solution by \(M_D = S - S_0 \cup \Lambda\), where subscript \(D\) stands for “dummy”.

Notice that \(M_O\) contains \(S- S_0\) and \(\hat{\text{OPT}}\). Therefore, if we sum up the above bounds for all clients we get \[\begin{align} \text{cost}(M_O) & \leq \sum_{p\in D^*_{\texttt{exp}}} \text{OPT}_p + \sum_{p\in D{-} D^*_{\texttt{exp}}} r(p)\,, \end{align}\] which by 16 is at most \(\sum_{p\in D^*_{\texttt{exp}}} \text{OPT}_p + \sum_{p\in D{-}D^*_{\texttt{exp}}} {4} \text{OPT}_p + O({\sqrt{\varepsilon}}) \cdot \text{opt}\) (assuming that that the sample \(W\) was successful). By swapping \(S_0\) with \(\hat{\text{OPT}}\) we thus get a \(1\)-approximation on the clients \(D^*_{\texttt{exp}}\) and a \({4}\)-approximation of the remaining clients (plus a small error term). Our goal in the next sections will thus be to approximate this swap. The next few steps will be aimed at simplifying (or guessing parts of) \(S_0\). The last step will then approximate \(\hat{\text{OPT}}\) and the remaining part of \(S_0\) by submodular function maximization subject to a partition matroid constraint. For that part, it will be helpful to have a bound on the dummy solution as well; specifically, a modified version of it (see 21 and [claim:Mprimebounds]). For intuition of those statements, let us here say that one can show \[\begin{align} \text{cost}(M_D) & \leq \sum_{p\in D^*_{\texttt{exp}}} {9}\,\text{OPT}_p + \sum_{p\in D- D^*_{\texttt{exp}}} r(p) + O({\sqrt{\varepsilon}}) \cdot \text{opt}\,\, \end{align}\] by observing that only the clients \(D^*_{\texttt{exp}}\) have a different connection cost in \(M_D\) than in \(M_O\). Moreover, each such client \(p\in D^*_{\texttt{exp}}\) that was previously assigned to a center \(c^* \in \hat{\text{OPT}}\), and belongs to cluster \(C^*\) of \(\text{OPT}\) with leader \(\ell\), has an associated dummy center \(\delta\) at squared distance \[{d^2(p,\delta)\leq 3d^2(p, c^*)+3d^2(c^*,\ell)+3d^2(\ell,\delta)\leq} {3}d^2(p, c^*) + {6}\text{avg}_{C^*, \text{OPT}} + {3}\varepsilon\gamma_{C^*}\] (assuming the set \(\mathcal{B}\) of balls is valid and using 3 with \(\gamma=3\)). Simplifications then give the stated inequality as \(|C^*| \text{avg}_{C^*, \text{OPT}} \leq (1+2\varepsilon) \sum_{p\in C^*} \text{OPT}_p\) and \(|C^*| \varepsilon\gamma_{C^*} = \varepsilon\sum_{p\in C^*} ( \text{OPT}_p + S_p)\).

6.6 Removal of Expensive Centers of \(S_0\)↩︎

Thus, it remains to show that our algorithm can yield a good approximation to the gain that the swap \((\hat{\text{OPT}}, S_0)\) would provide. Of course, it will not necessarily be optimum, but we can show it is enough to conclude the proof of our theorem. We focus on the iteration of the procedure that produces valid ball guesses. In particular, at this point given that \(W\) is successful and the ball guesses are valid, the number of balls in our ball guess is equal to \(|S_0|=|\hat{\text{OPT}}|{=|\Lambda|}\leq \frac{\log n}{\varepsilon^3}\). Our algorithm then makes use of the following procedure that takes as input the local search solution \(S\) and the set \(\Lambda\) of dummy centers.

6.6.0.1 \(expRem(S,\Lambda)\)

\(\mathcal{L}_{\text{exp}}\leftarrow \emptyset\) \(\mathcal{Q}\leftarrow \emptyset\) With probability 1/2 do the following: Consider the clustering induced by \(\Lambda\cup S-\mathcal{Q}\) and, for each \(c\in S-\mathcal{Q}\), let \(\text{cost}(c)\) be the cost of the cluster associated with \(c\) in the considered solution Sample one \(c\in S-\mathcal{Q}\) with probability \(\frac{\text{cost}(c)}{\sum_{c'\in S-\mathcal{Q}}\text{cost}(c')}\) Set \(\mathcal{Q}\leftarrow \mathcal{Q}\cup \{c\}\) \(\mathcal{L}_{\text{exp}}\leftarrow \mathcal{L}_{\text{exp}}\cup \{\mathcal{Q}\}\) return \(\mathcal{L}_{\text{exp}}\)

The goal of the above procedure is to guess a subset \(\mathcal{Q}\) of all “costly” centers in \(S_0\) (the set \(\mathcal{L}_{\text{exp}}\) contains all guesses). Specifically, if we consider the solution \((S - \mathcal{Q}) \cup \Lambda\), then we wish that the cost of the clusters corresponding to the remaining centers in \(S_0 - \mathcal{Q}\) is at most \(\varepsilon\cdot {\text{opt}}\). Formally, the clusters corresponding to \(S_0 - \mathcal{Q}\) in \(S - \mathcal{Q}\cup \Lambda\) consists of all the clients whose closest center in \(S - \mathcal{Q}\cup \Lambda\) is from \(S_0 - \mathcal{Q}\). If we let \(D' \subseteq D\) be those clients then the total cost in solution \(S- \mathcal{Q}\cup \Lambda\) of the clusters with centers in \(S_0 - \mathcal{Q}\) is defined as \[\sum_{p\in D'} d^2(p, S- \mathcal{Q}\cup \Lambda) = \sum_{p\in D'} d^2(p, S_0 - \mathcal{Q})\,,\] where the equality holds because of the definition \(D'\).

6.7 Consistently Assigning Clients in Mixed Solutions and the sets \(\mathcal{U}\) and \(\mathcal{R}\)↩︎

While our algorithm tries all possible \(\mathcal{Q}\in \mathcal{L}_{\text{exp}}\), we focus on the execution path of a set \(\mathcal{Q}\) of centers satisfying the properties of 20, i.e.,

1. \(\mathcal{Q}\subseteq S_0\), and

2. the cost in the solution \(S - \mathcal{Q}\cup \Lambda\) of the clusters with centers in \(S_0 - \mathcal{Q}\) is at most \(\varepsilon\cdot \text{opt}\).

We define the solution \(S_{\mathcal{Q}} := S- \mathcal{Q}\cup \Lambda\) and for a center \(c\in S_{\mathcal{Q}}\) we let \(S_{\mathcal{Q}}(c)\) be the subset of clients closest to \(c\) (i.e., assigned to \(c\)) in the solution \(S_{\mathcal{Q}}\). The second property above allows us to consider each cluster \(S_{\mathcal{Q}}(c)\), \(c\in S_0 - \mathcal{Q}\), as contracted by paying a small extra cost of \(\varepsilon\cdot \text{opt}\). In particular, this allows us to give a “consistent” assignment \(\mu_O\) of clients in the mixed solutions. Here, consistent means that all clients in \(S_Q(c)\), for \(c\in S_0 - \mathcal{Q}\), are assigned to the same center in the mixed solutions where all of \(S_0\) is removed.

6.7.0.1 Definition of \(\mu_O\) and \(\mu_D\).

We modify how he clients are assigned in the solution \(M_O\) to obtain the assignment \(\mu_O\). For every \(c \in S_0 - \mathcal{Q}\), all the clients in \(S_{\mathcal{Q}}(c)\) are reassigned to the same center of \(M_O\) as follows. Let \(c_1\) be the cluster of \(\hat{\text{OPT}}\) that is the closest to \(c\) and let \(c_2\) be the center of \(M_O-\hat{\text{OPT}}= S - S_0\) that is the closest to \(c\). The clients of \(S_{\mathcal{Q}}(c)\) are all assigned to either \(c_1\) or \(c_2\):

• If \(d^2(c,c_1) \le d^2(c,c_2)/5\), assign \(S_{\mathcal{Q}}(c)\) to \(c_1\); \(S_{\mathcal{Q}}(c)\) is called a type-1 cluster.

• Otherwise (\(d^2(c,c_1) > d^2(c,c_2)/5\)), assign \(S_{\mathcal{Q}}(c)\) to \(c_2\); \(S_{\mathcal{Q}}(c)\) is called a type-2 cluster.

Let \(T_1\) be the set of clients in a type-1 cluster and \(T_2\) the set of clients in a type-2 cluster. The remaining clients in \(D^*_{\texttt{exp}}- (T_1 \cup T_2)\) are assigned to their closest centers in \(\hat{\text{OPT}}\), and the remaining clients in \((D- D^*_{\texttt{exp}})- (T_1 \cup T_2)\) are assigned to their closest centers in \(S- S_0\). Notice that, even neglecting the extra cost due to the consistent assignment, the above assignment of clients in \(T_2\) is suboptimal when \(d(c,c_2)>d(c,c_1)>d(c,c_2)/5\). The motivation for this reassignment is technical. More specifically, it will be used in Lemma 21 to have a better upper bound on the cost of clients in \(C^*\cap T_1\) for a non-expensive cluster \(C^*\) (case b2). This leads to a larger upper bound on the cost of clients in \(C^*\cap T_2\) for an expensive cluster \(C^*\) (case a3), which is however tolerable.

We further obtain an assignment \(\mu_D\) of clients in \(M_D\) by modifying \(\mu_O\) as follows: each client \(p\) with \(\mu_O(p) \in \hat{\text{OPT}}\) is assigned to the associated dummy center by \(\mu_D\). So for each client \(p\) we have \(\mu_D(p) = \mu_O(p)\) unless \(\mu_O(p)\in \hat{\text{OPT}}\) in which case \(\mu_D(p)\) equals the dummy center associated with \(\mu_O(p)\).

6.7.0.2 Sets \(\mathcal{R}\) and \(\mathcal{U}\).

By the definition of type-1 and type-2 cluster, we can split the set of centers of \(S_0 - \mathcal{Q}\) into two groups \(\mathcal{R}\) and \(\mathcal{U}\). Let \(\mathcal{R}\) be the set of centers of \(S_0- \mathcal{Q}\) whose set of clients is completely assigned to a center of \(\hat{\text{OPT}}\) in the assignment \(\mu_O\), and let \(\mathcal{U}= S_0 - \mathcal{Q}- \mathcal{R}\) be the remaining ones that are assigned to centers in \(S- S_0\).

We end this section by analyzing the costs of the assignments \(\mu_O\) and \(\mu_D\). Recall that we use the notation \(\text{cost}(M_O,\mu_O)\) to denote the cost of \(M_O\) equipped with assignment \(\mu_O\) and, similarly \(\text{cost}(M_D, \mu_D)\) for the cost of \(M_D\) equipped with assignment \(\mu_D\). We say that the selection of \((W, \mathcal{B}, \mathcal{Q})\) is successful, if the sample \(W\) selected in 6.3 is successful, the set of balls \(\mathcal{B}\) from 6.4 is valid, and \(\mathcal{Q}\) selected in this section satisfies the properties of 20.

6.8 Removal of Cheap Centers in \(S_0\)↩︎

At this stage we have defined a clustering \(S_{\mathcal{Q}} = S - \mathcal{Q}\cup \Lambda\), where \(\Lambda\) is the set of dummy centers. Recall that the mixed solution \(M_O\) equals \(S-S_0 \cup \hat{\text{OPT}}\). So \(S_\mathcal{Q}- \Lambda\) represents progress compared to \(S\) in that we have already removed a subset \(\mathcal{Q}\) of \(S_0\). Furthermore, with respect to the assignment of clients \(\mu_O\), we have that the remaining centers \(S_0 - Q\), which we want to remove, are divided into two sets \(\mathcal{R}\) and \(\mathcal{U}\). The clusters with centers in \(\mathcal{R}\) are completely assigned (by \(\mu_O\)) to centers in \(\hat{\text{OPT}}\), and every cluster with center in \(\mathcal{U}\) is completely assigned to a center in \(S - S_0\). The task of this section is to make further progress by guessing the centers in \(\mathcal{U}\) and their reassignment to centers in \(S- S_0\). This turns out to be a difficult task, and we will instead “approximately” guess a set \(\widetilde{\mathcal{U}}\) and reassignment \(\widetilde{\mu}\) of the clients of those clusters. Specifically, we will output a solution \(S_{\mathcal{Q}\cup \widetilde{\mathcal{U}}} = S_{\mathcal{Q}} - \widetilde{\mathcal{U}}\) (obtained by removing \(\widetilde{\mathcal{U}}\) from \(S_{\mathcal{Q}}\)) and an assignment of clients \(\widetilde{\mu}\) to centers in \(S_\mathcal{Q}-\widetilde{\mathcal{U}}\). This is done through the following lemma. Recall that for \(c\in S_{\mathcal{Q}}\), \(S_{\mathcal{Q}}(c)\) is the set of clients assigned to \(c\) in \(S_\mathcal{Q}\), i.e., closest to \(c\) among all centers in \(S_\mathcal{Q}\).

In words, the above properties of \(\widetilde{\mu}\) say that we maintain the assignment of clients whose closest center remains the same, and the clients associated to the removed centers \(\widetilde{\mathcal{U}}\) are reassigned to centers not in \(\mathcal{R}\) at no higher cost (up to \(O({\sqrt{\varepsilon}} \cdot \text{opt})\)) than the cost of the reassignment of \(\mathcal{U}\) in the modified solution \(M_O\).

Let \(\mathcal{X}\) be the centers in \(S-S_0\) to which the points \(p\in S_{\mathcal{Q}}(c)\) with \(c\in \mathcal{U}\) are assigned according to \(\mu_O\). Notice that \(|\mathcal{X}|\leq |\mathcal{U}|\) and \(\mathcal{X}\cap (\mathcal{Q}\cup \mathcal{U}\cup \mathcal{R})=\emptyset\). We assume that the values of \(|\mathcal{U}|\), \(|\mathcal{R}|\), and \(|\mathcal{X}|\) are known. This can be achieved by trying, for each such value, all the integers between \(0\) and \(|S_0|\leq \frac{\log n}{\varepsilon^3}\). Let \(\ell:=|\mathcal{U}|+|\mathcal{R}|+|\mathcal{X}|\leq 2|\mathcal{U}|+|\mathcal{R}|\leq 2|S_0 - \mathcal{Q}|\). Then we run the recursive procedure \(cheapRem()\) described in the box with input \((\mathcal{U}',\mathcal{R}',\mathcal{X}',\mathcal{N},\widetilde{\mathcal{U}})=(\emptyset,\emptyset,\emptyset,\emptyset,\emptyset)\). We assume that \(cheapRem()\) has access to the quantities \(S\), \(\mathcal{Q}\), \(|\mathcal{U}|\), \(|\mathcal{R}|\), \(|\mathcal{X}|\) and \(\ell\), as well as to a global variable \(\mathcal{L}_{cheap}\), which is initialized to \(\emptyset\). The intuition for the parameters is as follows. Intuitively \(\mathcal{U}'\), \(\mathcal{R}'\) and \(\mathcal{X}'\) are subsets of \(\mathcal{U}\), \(\mathcal{R}\) and \(\mathcal{X}\), resp., that we have already identified, while \(\widetilde{\mathcal{U}}\) is the current value of the set \(\widetilde{\mathcal{U}}\) under construction. Intuitively, \(\mathcal{N}\) are the centers to which the points \(p\in S_{\mathcal{Q}}(c)\) with \(c\in \widetilde{\mathcal{U}}\) are reassigned, namely \(\widetilde{\mu}(p)\in \mathcal{N}\). At the end of the root call, \(\mathcal{L}_{cheap}\) contains the desired collection of sets \(\widetilde{\mathcal{U}}\) as in 22. Whenever some \(\widetilde{\mathcal{U}}\) is added to \(\mathcal{L}_{cheap}\), we define a corresponding \(\widetilde{\mu}\) as follows: For each \(p\in S_{\mathcal{Q}}(c)\), \(\widetilde{\mu}(p)=c\) if \(c\in S_{\mathcal{Q}}-\widetilde{\mathcal{U}}\). Otherwise, i.e., if \(c\in \widetilde{\mathcal{U}}\), \(\widetilde{\mu}(p)=next(c)\in \mathcal{N}\), where \(next(c)\) is defined in the recursive call when \(c\) is added to \(\widetilde{\mathcal{U}}\).

6.8.0.1 \(cheapRem(\mathcal{U}',\mathcal{R}',\mathcal{X}',\mathcal{N},\widetilde{\mathcal{U}})\)

Add \(\widetilde{\mathcal{U}}\) to \(\mathcal{L}_{cheap}\) and, for each \(c\in \widetilde{\mathcal{U}}\), set \(next(c)\) to the closest center to \(c\) in \(S-\mathcal{Q}-\widetilde{\mathcal{U}}-\widetilde{\mathcal{R}}\) halt Add \(\widetilde{\mathcal{U}}\) to \(\mathcal{L}_{cheap}\) and halt Let \(c\in S-\mathcal{Q}-\mathcal{R}'-\mathcal{X}'-\mathcal{N}-\widetilde{\mathcal{U}}\) minimize \(R(c):=\sum_{p\in S_{\mathcal{Q}}(c)}d^2(p,next(c))\), where \(next(c)\) is the closest center to \(c\) in \(S-\mathcal{Q}-\mathcal{U}'-\mathcal{R}'-\widetilde{\mathcal{U}}-\{c\}\). \(cheapRem(\mathcal{U}',\mathcal{R}',\mathcal{X}'\cup \{c\},\mathcal{N},\widetilde{\mathcal{U}})\) \(cheapRem(\mathcal{U}',\mathcal{R}'\cup \{c\},\mathcal{X}',\mathcal{N},\widetilde{\mathcal{U}})\) \(cheapRem(\mathcal{U}',\mathcal{R}',\mathcal{X}',\mathcal{N}\cup \{next(c)\},\widetilde{\mathcal{U}}\cup \{c\})\) \(cheapRem(\mathcal{U}',\mathcal{R}'\cup \{next(c)\},\mathcal{X}',\mathcal{N},\widetilde{\mathcal{U}})\) \(cheapRem(\mathcal{U}'\cup \{next(c)\},\mathcal{R}',\mathcal{X}',\mathcal{N},\widetilde{\mathcal{U}})\) \(cheapRem(\mathcal{U}',\mathcal{R}',\mathcal{X}',\mathcal{N}\cup \{next(c)\},\widetilde{\mathcal{U}}\cup \{c\})\)

6.9 Mixed Solutions after the Removal of \(\widetilde{\mathcal{U}}\)↩︎

We let \((S_{\mathcal{Q}\cup \widetilde{\mathcal{U}}}, \widetilde{\mu})\) be the output of a successful run of the \(cheapRem\) procedure, i.e., where \(\widetilde{\mathcal{U}}\in \mathcal{L}_{cheap}\) and its associated assignment \(\widetilde{\mu}\) satisfies the properties of 22. Recall that \(S_{\mathcal{Q}\cup \widetilde{\mathcal{U}}}\) consists of the centers \(S' = S - \mathcal{Q}- \widetilde{\mathcal{U}}\) and the dummy centers \(\Lambda\). Further recall that \(M_O\) consists of centers \(S- S_0 = S - \mathcal{Q}- \mathcal{U}- \mathcal{R}\) and \(\hat{\text{OPT}}\). \(M_D\) is the same set of centers except that \(\hat{\text{OPT}}\) is replaced by \(\Lambda\). Notice that some centers of \(M_O\) and \(M_D\) are now removed if \(\widetilde{\mathcal{U}}- \mathcal{U}\neq \emptyset\). To take care of this, we modify these mixed solutions to obtain \(M'_O\) and \(M'_D\).

We first define \(M'_O\); the definition of \(M'_D\) is then very similar. The centers of \(M'_O\) are \(S- \mathcal{Q}- \widetilde{\mathcal{U}} - \mathcal{R}\) and \(\hat{\text{OPT}}\). Hence, the difference between \(M_O\) and \(M'_O\) is that in \(M_O\) we remove \(\mathcal{U}\) from \(S\) and in \(M'_O\) we remove \(\widetilde{\mathcal{U}}\). In other words, \(M'_O = M_O \cup \mathcal{U}- \widetilde{\mathcal{U}}\). Similarly, we let \(M'_D = M_D \cup \mathcal{U}- \widetilde{\mathcal{U}}\). We update the assignment \(\mu_O\) of \(M_O\) to an assignment \(\mu'_O\) of \(M'_O\) (and the assignment \(\mu_D\) to \(\mu'_D\)). Recall that \(\mu_O\) assigns each client \(p\) to its closest center in \(\hat{\text{OPT}}\) (if \(p\in D^*_{\texttt{exp}}\)) or its closest center in \(S- S_0\) (if \(p\in D- D^*_{\texttt{exp}}\)) except for those clients that belong to a cluster \(S_{\mathcal{Q}}(c)\) with \(c\in \mathcal{U}\cup \mathcal{R}\). Indeed, the clients in \(S_\mathcal{Q}(c)\) are all assigned to the center of \(\hat{\text{OPT}}\) that is closest to \(c\) if \(c\in \mathcal{R}\), and if \(c\in \mathcal{U}\) they are all assigned to the center in \(S - S_0\) that is closest to \(c\). We also recall that \(\mu_D\) is the same as \(\mu_O\) except when \(\mu_O(p) \in \hat{\text{OPT}}\) in which case \(\mu_O(p)\) is replaced by its corresponding dummy center.

6.9.0.1 Definitions of \(\mu_O'\) and \(\mu_D'\).

For clients \(p\in S_{\mathcal{Q}}(c)\) with \(c \not\in \mathcal{U}\cup \widetilde{\mathcal{U}}\cup \Lambda\), we define \(\mu_O'(p) = \mu_O(p)\). For \(p\in S_{\mathcal{Q}}(c)\) with \(c\in \Lambda\), we let \(\mu'_O(p)\) be the closest center in \(\hat{\text{OPT}}\). Finally, for \(p\in S_{\mathcal{Q}}(c)\) with \(c\in \mathcal{U}\cup \widetilde{\mathcal{U}}\), we define \(\mu'_O(p) = \widetilde{\mu}(p)\).

The assignment \(\mu'_D\) is obtained in the same way from \(\mu_D\) with the difference that we use \(\Lambda\) instead of \(\hat{\text{OPT}}\): For clients \(p\in S_{\mathcal{Q}}(c)\) with \(c \not\in \mathcal{U}\cup \widetilde{\mathcal{U}}\cup \Lambda\), we define \(\mu_D'(p) = \mu_D(p)\). For \(p\in S_{\mathcal{Q}}(c)\) with \(c\in \Lambda\), we let \(\mu'_D(p)\) be the closest center in \(\Lambda\). Finally, for \(p\in S_{{\mathcal{Q}}}(c)\) with \(c\in \mathcal{U}\cup \widetilde{\mathcal{U}}\), we define \(\mu'_D(p) = \widetilde{\mu}(p)\).

This completes the definition of \(\mu'_O\) and \(\mu'_D\). We remark that we have \(\mu'_O(p) = \mu_O(p)\) and \(\mu'_D(p) = \mu_D(p)\) for all \(p\in S_\mathcal{Q}(c)\) with \(c\not \in \mathcal{U}\cup \widetilde{\mathcal{U}}\cup \Lambda\).

We continue by arguing that \(\mu'_O\) is well-defined, i.e., that \(\mu'_O(p) \in M_O'\) for every \(p\in D\) (the proof for \(\mu'_D\) is the same). This is immediate for a client \(p\in S_\mathcal{Q}(c)\) with \(c \in \Lambda\). For a client \(p\in S_{{\mathcal{Q}}}(c)\) with \(c\in \cal U \cup \widetilde{\mathcal{U}}\), it holds because, by 22, \(\widetilde{\mu}(p)\) equals \(c\in M'_O\) if \(c \in \mathcal{U}- \widetilde{\mathcal{U}}\) and otherwise \(\widetilde{\mu}(p) \in S- \mathcal{Q}- \widetilde{\mathcal{U}}- \mathcal{R}\subseteq M'_O\). It remains to verify that \(\mu_O(p) \in M_O'\) for a client \(p\in S_\mathcal{Q}(c)\) with \(c\not \in \mathcal{U}\cup \widetilde{\mathcal{U}}\cup \Lambda\). If \(c\in \mathcal{R}\), we have \(\mu_O(p) \in \hat{\text{OPT}}\subseteq M'_O\). Similarly, if \(c \not \in \mathcal{U}\cup \widetilde{\mathcal{U}}\cup \Lambda \cup \mathcal{R}\) and \(p\in D^*_{\texttt{exp}}\) we have \(\mu_O(p) \in \hat{\text{OPT}}\subseteq M'_O\). In the remaining case when \(c \not \in \mathcal{U}\cup \widetilde{\mathcal{U}}\cup \Lambda \cup \mathcal{R}\) and \(p\in D- D^*_{\texttt{exp}}\), \(\mu_O(p)\) equals \(p\)’s closest center in \(S- S_0\). As \(p\in S_\mathcal{Q}(c)\), we thus have \(\mu_O(p) = c \in S - (S_0 \cup \widetilde{\mathcal{U}}) \subseteq M'_O\).

The following upper bound on \(\text{cost}(M'_O, \mu'_O) + \text{cost}(M'_D, \mu'_D)\) is a fairly immediate consequence of the assumption that the cheap-removal process is successful (Property 3c of Lemma 22). Similarly to before, we say that the selection of \((W, \mathcal{B}, \mathcal{Q}, \widetilde{\mathcal{U}},\widetilde{\mu})\) is successful, if the sample \(W\) selected in 6.3 is successful, the set of balls \(\mathcal{B}\) from 6.4 is valid, \(\mathcal{Q}\) selected in 6.6 satisfies the properties of 20, and \(\widetilde{\mathcal{U}},\widetilde{\mu}\) selected in this section satisfies the properties of 22.

To better understand \(\text{cost}(M'_O, \mu'_O)\), let us consider the cost \(d(p, \mu'_O(p))\) of a single client \(p\):

• If \(p \in S_{\mathcal{Q}}(c)\) with \(c\not\in \widetilde{\mathcal{U}}\cup \mathcal{R}\) then \(\widetilde{\mu}(p) = c\) by Property 3a of 22. Moreover, \(\mu'_O(p)\) assigns \(p\) either to a center in \(\hat{\text{OPT}}\) or to one in \(S - \mathcal{Q}- \widetilde{\mathcal{U}}- \cal R\) and we have \(d^2(p, c) \leq d^2(p, S - \mathcal{Q}- \widetilde{\mathcal{U}}- \cal R)\). Hence, in either case, we have \[d^2(p, \mu'_O(p)) \geq d^2(p, \{\widetilde{\mu}(p)\} \cup \hat{\text{OPT}}) = d^2(p, \{\widetilde{\mu}(p)\} \cup \hat{\text{OPT}}\cup \Lambda) \,.\]

• If \(p\in S_{\mathcal{Q}}(c)\) with \(c\in \mathcal{R}\) then \(\mu'_O(p) = \mu_O(p) = c^*\in \hat{\text{OPT}}\) where \(c^*\) is the center of \(\hat{\text{OPT}}\) that is closest to \(c\). So \[d^2(p, \mu_O'(p)) = d^2(p, c^*)\,.\]

• If \(p \in S_{\mathcal{Q}}(c)\) with \(c\in \widetilde{\mathcal{U}}\) then \(\mu'_O(p) = \widetilde{\mu}(p)\) and so \[d^2(p, \mu_O'(p)) \geq d^2(p, \{\widetilde{\mu}(p)\} \cup \hat{\text{OPT}}\cup \Lambda).\]

Similarly, we can analyze \(d^2(p, \mu_D(p))\) by replacing \(\hat{\text{OPT}}\) with the set \(\Lambda\) of dummy centers (and \(c^*\) by its associated dummy center). Summarizing, we have \[\begin{align} \text{cost}(M'_O, \mu'_O) &\geq \sum_{c\in S_{\mathcal{Q}} - \mathcal{R}}\sum_{p\in S_{\mathcal{Q}}(c)} d^2\left(p, \{c\} \cup \hat{\text{OPT}}\cup \Lambda\right) + \sum_{c\in \mathcal{R}} \min_{c' \in \hat{\text{OPT}}} \sum_{p\in S_{\mathcal{Q}}(c)} d^2\left(p, c' \right)\\ \intertext{and} \text{cost}(M'_D, \mu'_D) &\geq \sum_{c\in S_{\mathcal{Q}} - \mathcal{R}}\sum_{p\in S_{\mathcal{Q}}(c)} d^2\left(p, \{c\} \cup \Lambda\right) + \sum_{c\in \mathcal{R}} \min_{c'\in \Lambda}\sum_{p\in S_{\mathcal{Q}}(c)} d^2\left(p, c'\right) \end{align}\]

Furthermore, by Properties 3a and 3b of 22, we have that, for every \(c\in \mathcal{R}\), the set \(S_{\mathcal{Q}}(c)\) of clients assigned to \(c\) in \(S_{\mathcal{Q}}\) equals the set \(\widetilde{\mu}^{-1}(c)\) of clients assigned by \(\widetilde{\mu}\). So, if we let \(\{C_c\}_{c\in S_{\mathcal{Q}\cup \widetilde{\mathcal{U}}}}\) be the partitioning of the set \(D\) of clients according to \(\widetilde{\mu}\), i.e., \(C_c = \{p\in D\mid \widetilde{\mu}(p) =c\}\), we can thus rewrite the above bounds on the cost as \[\begin{align} \text{cost}(M'_O, \mu'_O) &\geq \sum_{c\not \in \mathcal{R}}\sum_{p\in C_c} d^2\left(p, \{c\} \cup \hat{\text{OPT}}\cup \Lambda\right) + \sum_{c\in \mathcal{R}}\min_{c' \in \hat{\text{OPT}}}\sum_{p\in C_c} d^2\left(p, c'\right) \tag{30}\\ \intertext{and} \text{cost}(M'_D, \mu'_D) &\geq \sum_{c \not\in \mathcal{R}}\sum_{p\in C_c} d^2\left(p, \{c\} \cup \Lambda\right) + \sum_{c\in \mathcal{R}}\min_{c' \in \Lambda}\sum_{p\in C_c} d^2\left(p, c'\right)\,.\tag{31} \end{align}\] In the next section, we use these bounds to give a polynomial-time algorithm that outputs a solution \(S^*\) with \(k\) centers so that (see 23) \[\begin{align} \text{cost}(S^*) \leq \frac{(1+{6\sqrt{\varepsilon}})}{2} \left( \text{cost}(M'_O, \mu'_O) + \text{cost}(M'_D, \mu'_D) \right) + O({\sqrt{\varepsilon}}\cdot \text{opt}) \leq \sum_{p\in D} {5} \text{OPT}_p + O({\sqrt{\varepsilon}} \cdot \text{opt}) \end{align}\] where the second inequality is by [claim:Mprimebounds]. So the proof of 23 in the next subsection is the final step in the proof of [thm:mainadditivecenters] (see also 6.11 where we put everything together to prove [thm:mainadditivecenters]).

6.10 Finding \(S^*\) via Submodular Optimization↩︎

In this section, we give a polynomial-time algorithm for finding the solution \(S^*\) by reducing the problem to maximizing a submodular function subject to a partition matroid constraint. Specifically, we prove the following lemma:

Recall that at this point the algorithm calculated the local search solution \(S\), sampled \(W\), guessed \(\mathcal{B}\), \(\mathcal{Q}\), and \(\widetilde{\mathcal{U}}\) with the assignment \(\widetilde{\mu}\) of the clustering \(S_{\mathcal{Q}\cup \widetilde{\mathcal{U}}}\). We assume that all these choices were successful guesses, i.e., that \((W, \mathcal{B}, \mathcal{Q}, \widetilde{\mathcal{U}}, \widetilde{\mu})\) is successful. Further recall the notation that we partition the clients into the clusters \(\{C_c\}_{c\in S_{\mathcal{Q}\cup \widetilde{\mathcal{U}}}}\) where \(C_c = \{ p\in D: \widetilde{\mu}(p) =c\}\).

Our goal is to give a polynomial-time algorithm that finds approximations of the sets \(\hat{\text{OPT}}\) and \(\mathcal{R}\) that only have slightly worse cost. We start by defining the feasible set of candidates for \(\hat{\text{OPT}}\) as a partition matroid. Recall that \(\hat{\text{OPT}}\) contains one center from each ball in \(\mathcal{B}\).

6.10.0.1 Definition of partition matroid \(\mathcal{M}_\mathcal{B}\).

We define the partition matroid that captures the constraint that we wish to open a center in each ball in \(\mathcal{B}\). Let \(F_\mathcal{B}\) be the (multi) subset of facility/center locations containing \(B\cap F\) for each ball \(B\in \mathcal{B}\). If a center \(c\) is in multiple balls in \(\mathcal{B}\), then \(F_\mathcal{B}\) contains one distinct copy of \(c\) for each ball. We let \(F_{\mathcal{B}}(B) \subseteq F_\mathcal{B}\) denote the centers associated with \(B\in \mathcal{B}\). These sets satisfy the following two properties:

• The sets \(F_{\mathcal{B}}(B)\) partition \(F_\mathcal{B}\).

• For \(B \in \mathcal{B}\), \(F_{\mathcal{B}}(B)\) contains (a copy of) every center in \(F\cap B\).

The first property holds since we took a unique copy of each center for each ball, and they are thus disjoint: \(F_\mathcal{B}(B) \cap F_\mathcal{B}(B') = \emptyset\) for distinct \(B,B' \in \mathcal{B}\). Indeed, while it is not better for the cost to open multiple copies of a center, we make the copies to ensure the above two properties. This allows us to define the partition matroid \(\mathcal{M}_\mathcal{B}= (F_\mathcal{B}, \mathcal{I})\) where \[\mathcal{I}= \{X \subseteq F_\mathcal{B}: |X\cap F_\mathcal{B}(B)| \leq 1for everyB\in \mathcal{B}\}\,.\] Moreover, since there is exactly one ball in \(\mathcal{B}\) for each center in \(\hat{\text{OPT}}\), we have \(\hat{\text{OPT}}\in \mathcal{I}\) (where we slightly abuse notation as we should take the copy of center \(c^*\in \hat{\text{OPT}}\) that belongs to its associated ball).

6.10.0.2 Core, concentrated and hit clusters.

Our goal is now to define a submodular function \(f\) so that we obtain a good approximation to \(\hat{\text{OPT}}\) by maximizing \(f\) over the matroid constraint \(\mathcal{M}_\mathcal{B}\). In particular, the domain of \(f\) is every subset of \(F_\mathcal{B}\). However, we need some additional steps before defining \(f\). In particular, we introduce the concept of the core of a cluster \(C_c\) and the notions of concentrated and hit clusters, which allow us to simplify the structure of centers in \(\mathcal{R}\). For a cluster \(C_c\), we define the core of \(C_c\) as \[\begin{gather} C^{\text{core}}_c = \left\{p \in C_c: d^2(p, c) \leq \varepsilon\cdot \frac{\text{cost}(S)}{|\mathcal{R}|\cdot |C_c|}\right\}\,. \end{gather}\] We further say that cluster \(C_c\) is concentrated if \[\begin{gather} |C^{\text{core}}_c| \geq (1-\varepsilon) |C_c| \end{gather}\] and it is hit by a set \(X\) of centers if there is a point \(p \in C^{\text{core}}_c\) such that \[\begin{gather} d^2(p, X) < d^2(p,c)\,. \end{gather}\] For shortness we will sometimes say that a center is concentrated (resp., hit), if the corresponding cluster is so.

6.10.0.3 Guessing the centers of \(\mathcal{R}\) that are not concentrated.

We further simplify the task of finding the set of centers \(\mathcal{R}\) by guessing the centers of \(\mathcal{R}\) that are not concentrated. Specifically, partition the set \(\mathcal{R}\) into \(\mathcal{R}_0\) and \(\mathcal{R}_1\), where \(\mathcal{R}_1\) contains those centers of \(\mathcal{R}\) that are concentrated and \(\mathcal{R}_0\) contains those that are not. We can correctly guess \(\mathcal{R}_0\) in polynomial time since the following simple claim shows that it is a subset of the centers whose cluster costs at least \(\varepsilon^2 \cdot \text{opt}/|\mathcal{R}|\), of which there are only \(O(|\mathcal{R}|/\varepsilon^2)\) many, and so guessing \(\mathcal{R}_0\) can thus be done in time \(2^{O(|\mathcal{R}|/\varepsilon^2)}\).

Now let \(\mathcal{P}\) be the potential centers that can be in \(\mathcal{R}\). Specifically, we let \(\mathcal{P}\) be the set that contains a center \(c\) in \(S -\mathcal{Q}- \widetilde{\mathcal{U}}\) if \(C_c\) equals the set of clients in \(S_{\mathcal{Q}}(c)\), where we recall that \(S_{\mathcal{Q}}(c)\) is the set of clients that are closest to \(c\) in the clustering \(S_{\mathcal{Q}}\). Notice that the algorithm has all the information \(S, \mathcal{B}, \mathcal{Q}, \widetilde{\mathcal{U}}\) and \(\widetilde{\mu}\) to calculate \(\mathcal{P}\). Furthermore, by 22, we have that no clients from \(\widetilde{\mathcal{U}}\) were reassigned by \(\widetilde{\mu}\) to a center in \(\mathcal{R}\). So for \(c\in \mathcal{R}\) we have \(C_c = S_{\mathcal{Q}}(c)\). Hence, \(\cal R \subseteq \mathcal{P}\) and we can guess \(\mathcal{R}_0\) as follows:

6.10.0.4 Guessing \(\mathcal{R}_0\)

1. Let \(\mathcal{C}\) be the clusters of \(\mathcal{P}\) whose cost is at least \(\varepsilon^2 \cdot \text{opt}/|\mathcal{R}|\).

2. Output each subset of \(\mathcal{C}\).

By the above claim and the definition of \(\mathcal{P}\), one of the outputs is \(\mathcal{R}_0\). Furthermore, the above guessing procedure outputs a family of polynomial many subsets. Indeed, we have \(|\mathcal{R}| \leq \log(n)/\varepsilon^3\) by 12. Moreover, the total cost of the clusters in \(\mathcal{P}\) is \(O({\text{opt}})\). To see this notice that the clusters corresponding to \(\mathcal{P}\) is a subset of the clusters of \(S_{\mathcal{Q}}\) and the cost of \(S_{\mathcal{Q}}\) is at most \(\text{cost}(S - S_0 \cup \Lambda)\) (since \(S_{\mathcal{Q}}\supseteq S- S_0 \cup \Lambda\)), which has cost at most \(O(\text{opt})\) by [lemma:S0properties]. This implies that \(|\mathcal{C}| = O(\log(n)/\varepsilon^5)\), so the total number of subsets is \(n^{\varepsilon^{-O(1)}}\). The algorithm proceeds by trying all possible subsets in the output, and we analyze the algorithm when it takes the correct guess of \(\mathcal{R}_0\).

6.10.0.5 Definition of the submodular function \(f\).

We first define another function \(g\) on the same domain as \(f\), i.e., on all subsets \(X\) of \(F_\mathcal{B}\). We will then define \(f\) by \(f(X) = g(\emptyset) - g(X)\). Let the closed cost of a cluster \(C_c\) be defined as \[\text{closedcost}_c(X) := \begin{cases} \sum_{p\in C_c} d^2(p, \{c\} \cup X \cup \Lambda) & if C_c is hit by X, \\ \min_{c'\in X \cup \Lambda} \sum_{p\in C^{\text{core}}_c} d^2(p, c') + \sum_{p\in C_c - C^{\text{core}}_c} d^2(p, \{c\} \cup X \cup \Lambda) & otherwise. \end{cases}\] Further, let \(\mathcal{P}_1\) be the potential centers for \(\mathcal{R}_1\): it contains each center \(c\in \mathcal{P}- \mathcal{R}_0\) so that \(c\) is concentrated, i.e., \(|C^{\text{core}}_c| \geq (1-\varepsilon) |C_c|\). We remark that the algorithm can calculate this set \(\mathcal{P}_1\) as it only depends on \(\mathcal{P}\), the guessed set \(\mathcal{R}_0\), the value \(\text{cost}(S)\), and the clusters \(C_c\) defined by \(\widetilde{\mu}\). Moreover, by definition, we have \(\mathcal{R}_1 \subseteq \mathcal{P}_1\). For a subset \(X \subseteq F_\mathcal{B}\), we then then define \(g(X)\) to be the minimum value of \[\begin{gather} \sum_{c \not \in \mathcal{R}_0 \cup \mathcal{R}'_1}\sum_{p\in C_c} d^2\left(p, \{c\} \cup X \cup \Lambda\right) + \sum_{c\in \mathcal{R}_0} \sum_{p\in C_c} d^2(p, X \cup \Lambda) + \sum_{c\in \mathcal{R}'_1} \text{closedcost}_c(X) \end{gather}\] over all subsets \(\mathcal{R}'_1 \subseteq \mathcal{P}_1\) with \(|\mathcal{R}_1'| = |\mathcal{R}_1| = |\mathcal{R}| - |\mathcal{R}_0|\). In words, over the best \(\mathcal{R}_1'\), \(g(X)\) is the cost of the solution obtained by removing the centers \(\mathcal{R}_0 \cup \mathcal{R}_1'\) and assigning clients as follows:

• If \(p\in C_c\) for a remaining center \(c \not\in \mathcal{R}_0 \cup \mathcal{R}_1'\) or \(p \not \in C^{\text{core}}_c\) with \(c\in \mathcal{R}'_1\), \(p\) is assigned to its closest center in \(\{c\} \cup X \cup \Lambda\).

• If \(p \in C_c\) for a removed center \(c\in \mathcal{R}_0\), \(p\) is assigned to its closest center in \(X \cup \Lambda\).

• If \(p \in C^{\text{core}}_c\) for a removed center \(c\in \mathcal{R}'_1\), \(p\) is assigned to its closest center in \(\{c\} \cup X \cup \Lambda\) if \(X\) hits \(C_c\), and otherwise all clients in \(C^{\text{core}}_c\) are assigned to the same center \(c'\in X \cup \Lambda\) that minimizes the cost.

We remark that this assignment is infeasible in the sense that it may assign clients to removed centers in \(\mathcal{R}_1\). Nevertheless, we relate the values of \(g(\emptyset)\) and \(g(\hat{\text{OPT}})\) to \(\text{cost}(M'_D, \mu'_D)\) and \(\text{cost}(M'_O, \mu'_O)\), respectively; and, we show that given an \(X \subseteq F_\mathcal{B}\) that is independent in \(M_\mathcal{B}\), i.e., \(X\in \mathcal{I}\), we can in polynomial-time output \(k\) centers \(S^*\) whose cost is at most \((1+{6\sqrt{\varepsilon}})g(X) + {O(\sqrt{\varepsilon}}) \text{opt}\). Finally, the definition of \(g\) allows us to prove that \(f\) (defined by \(f(X) = g(\emptyset) - g(X)\)) is a non-negative monotone submodular function.

We give the proof of the lemma in the next subsection. We explain here how it implies 23.

6.10.1 Proof of 24↩︎

6.10.1.1 Proof of Property 1.

Given \(X \subseteq F_{\mathcal{B}}\), we argue that we can evaluate \(g(X)\) in polynomial time. For a center \(c \in \mathcal{P}_1\), define \[\text{cost-inc}_c(X) = \text{closedcost}_c(X) - \sum_{p\in C_c} d^2(p, \{c\} \cup X \cup \Lambda)\,.\] In other words, using that \(d^2(p,c ) \leq d^2(p, X)\) for \(p\in C^{\text{core}}_c\) if \(C_c\) is not hit by \(X\), \[\begin{align} \text{cost-inc}_c(X) = \begin{cases} 0 & if C_c is hit by X, \\ \min_{c'\in X \cup \Lambda} \sum_{p\in C^{\text{core}}_c} (d^2(p, c') - d^2(p,\{c\} \cup \Lambda)) & otherwise. \end{cases} \label{eq:increase} \end{align}\tag{32}\] We remark that \(\text{cost-inc}_c(X)\geq 0\) by definition. Let \(D_0\) denote the clients belonging to clusters \(C_c\) with \(c\in \mathcal{R}_0\). With this notation, \(g(X)\) is the minimum value \[\begin{align} \sum_{p\in D- D_0} d^2(p, \{\widetilde{\mu}(p)\} \cup X \cup \Lambda) + \sum_{p\in D_0} d^2(p, X \cup \Lambda) + \sum_{c \in \mathcal{R}'_1} \text{cost-inc}_c(X) \label{eq:nicedefofg} \end{align}\tag{33}\] over all subsets \(\mathcal{R}'_1 \subseteq \mathcal{P}_1\) with \(|\mathcal{R}'_1| = |\mathcal{R}_1|\). We further have that the value of \(\text{cost-inc}_c(X)\) for \(c\in \mathcal{R}'_1\) is independent of other centers in \(\mathcal{R}'_1\). We can thus obtain the best choice of \(\mathcal{R}'_1\) by evaluating \(\text{cost-inc}_c(X)\) for every \(c\in \mathcal{P}_1\) and choose the \(|\mathcal{R}_1|\) ones of smallest value. Notice that the algorithm knows \(|\mathcal{R}_1|\) since it equals \(|S_0| - |\mathcal{Q}| - |\widetilde{\mathcal{U}}| - |\mathcal{R}_0| =|\mathcal{B}| - |\mathcal{Q}| - |\widetilde{\mathcal{U}}| - |\mathcal{R}_0|\). After we have obtained \(\mathcal{R}'_1\) in polynomial time, we can evaluate \(g(X)\) in polynomial time by simply calculating the sum.

6.10.1.2 Proof of Property 2.

In the proof of this property it will be convenient to use definition 33 of \(g\). We start by verifying that \(f\) is monotone. I.e., that \(f(Y) \geq f(X)\) when \(X \subseteq Y\), which is equivalent to verifying that \(g(Y) \leq g(X)\). Let \(g(X)\) equal \[\sum_{ p\in D- D_0} d^2(p, \{\widetilde{\mu}(p)\} \cup X \cup \Lambda) + \sum_{p\in D_0} d^2(p, X \cup \Lambda)+ \sum_{c \in \mathcal{R}'_1} \text{cost-inc}_c(X)\] for some set \(\mathcal{R}'_1\). Then \(g(Y)\) is at most \[\sum_{p\in D- D_0} d^2(p, \{\widetilde{\mu}(p)\} \cup Y \cup \Lambda)+ \sum_{p\in D_0} d^2(p, {Y} \cup \Lambda) + \sum_{c \in \mathcal{R}'_1} \text{cost-inc}_c(Y)\] As trivially \(d^2(p, \{\widetilde{\mu}(p)\} \cup Y \cup \Lambda) \leq d^2(p, \{\widetilde{\mu}(p)\}\cup X \cup \Lambda)\), \(d^2(p, Y \cup \Lambda) \leq d^2(p, X \cup \Lambda)\) and \(\text{cost-inc}_c(Y) \leq \text{cost-inc}_c(X)\) (see 32 ), we have \(g(Y) \leq g(X)\) as required. Moreover, as \(f(\emptyset) = 0\) by definition, non-negativity follows.

We proceed to verify that \(f\) is submodular, i.e., that for every \(X \subseteq Y \subseteq F_B\) and \(c_0 \in F_B - Y\), \(f(X \cup \{c_0\}) - f(X) \geq f(Y \cup \{c_0\}) - f(Y)\), or equivalently \[\begin{align} g(X \cup \{c_0\}) - g(X) \leq g(Y \cup \{c_0\}) - g(Y)\,. \label{eq:submodularineq} \end{align}\tag{34}\] Let \(\mathcal{R}_1^X, \mathcal{R}_1^Y\), and \({\mathcal{R}}_1^{Y+c_0}\) be subsets of \(\mathcal{P}_1\) of cardinality \(|\mathcal{R}_1|\) such that \[\begin{align} g(X) &= \sum_{p\in D- D_0} d^2(p, \{\widetilde{\mu}(p)\} \cup X \cup \Lambda) + \sum_{p\in D_0} d^2(p, X \cup \Lambda) + \sum_{c \in \mathcal{R}^X_1} \text{cost-inc}_c(X)\\ g(Y) &= \sum_{p\in D- D_0} d^2(p, \{\widetilde{\mu}(p)\} \cup Y \cup \Lambda)+ \sum_{p\in D_0} d^2(p, Y \cup \Lambda) + \sum_{c \in \mathcal{R}^Y_1} \text{cost-inc}_c(Y)\\ g(Y \cup \{c_0\}) & = \sum_{p\in D- D_0} d^2(p, \{\widetilde{\mu}(p)\} \cup (Y \cup \{c_0\}) \cup \Lambda)+ \sum_{p\in D_0} d^2(p, (Y \cup \{c_0\}) \cup \Lambda) \\ & + \sum_{c \in {\mathcal{R}}^{Y+c_0}_1} \text{cost-inc}_c(Y\cup \{c_0\}) \end{align}\] We shall define \({\mathcal{R}}^{X+c_0}_1\) to be a subset of \(\mathcal{P}_1\) of cardinality \(|\mathcal{R}_1|\) so that the upper bound \[\begin{align} g(X\cup \{c_0\}) & \leq \sum_{p\in D- D_0} d^2(p, \{\widetilde{\mu}(p)\} \cup (X\cup \{c_0\}) \cup \Lambda)+ \sum_{p\in D_0} d^2(p, (X \cup \{c_0\}) \cup \Lambda) \\ & + \sum_{c \in {\mathcal{R}}^{X+c_0}_1} \text{cost-inc}_c(X \cup \{c_0\}) \end{align}\] allows us to verify Inequality 34 . First notice that no matter the definition of \({\mathcal{R}}^{X+c_0}_1\), we have that \[\begin{align} & \sum_{p\in D- D_0} (d^2(p, \{\widetilde{\mu}(p)\} \cup (X\cup \{c_0\}) \cup \Lambda) - d^2(p, \{\widetilde{\mu}(p)\} \cup X \cup \Lambda) )\nonumber \\ & +\sum_{p\in D_0} (d^2(p, (X\cup \{c_0\}) \cup \Lambda) - d^2(p, X \cup \Lambda) ) \label{eq:boringsubmodular1} \end{align}\tag{35}\] is upper bounded by \[\begin{align} & \sum_{p\in D- D_0} (d^2(p, \{\widetilde{\mu}(p)\} \cup (Y \cup \{c_0\}) \cup \Lambda) - d^2(p, \{\widetilde{\mu}(p)\} \cup Y \cup \Lambda)) \nonumber \\ & + \sum_{p\in D_0}(d^2(p, (Y \cup \{c_0\}) \cup \Lambda) - d^2(p, Y \cup \Lambda))\,. \label{eq:boringsubmodular2} \end{align}\tag{36}\] Indeed, if we let \(D' \subseteq D- D_0\) be the subset of clients for which \(d^2(p, c_0) < d^2(p, \{\widetilde{\mu}(p)\} \cup Y \cup \Lambda)\) and \(D'_0 \subseteq D_0\) be the subset of clients for which \(d^2(p, c_0) < d^2(p, Y \cup \Lambda)\), then 36 equals \[\sum_{p\in D'} (d^2(p, c_0) - d^2(p, \{\widetilde{\mu}(p)\} \cup Y \cup \Lambda))+ \sum_{p\in D'_0} (d^2(p, c_0) - d^2(p, Y \cup \Lambda))\] and we have \[\eqref{eq:boringsubmodular1} \leq \sum_{p\in D'} (d^2(p, c_0) - d^2(p, \{\widetilde{\mu}(p)\} \cup X \cup \Lambda)) + \sum_{p\in D'_0} (d^2(p, c_0) - d^2(p, X \cup \Lambda))\,.\] It follows that 35 is at most 36 because \[\sum_{p\in D'} ( d^2(p, \{\widetilde{\mu}(p)\} \cup Y \cup \Lambda) - d^2(p, \{\widetilde{\mu}(p)\} \cup X \cup \Lambda)) + \sum_{p\in D'_0} ( d^2(p, Y \cup \Lambda) - d^2(p, X \cup \Lambda)) \leq 0\,.\]

To prove that \(f\) is submodular, it is thus sufficient to define \({\mathcal{R}}^{X+c_0}_1\) so that \[\begin{align} \sum_{c\in {\mathcal{R}}^{X+c_0}_1} \text{cost-inc}_c(X \cup \{c_0\}) - \sum_{c\in {\mathcal{R}}^X_1} \text{cost-inc}_c(X)\leq \sum_{c\in {\mathcal{R}}^{Y+c_0}_1} \text{cost-inc}_c(Y \cup \{c_0\}) - \sum_{c\in {\mathcal{R}}^Y_1} \text{cost-inc}_c(Y)\,. \label{eq:mainsubmodular} \end{align}\tag{37}\] To this end, let \(\Delta_{\text{new}}\) contain the centers in \({\mathcal{R}}^{Y+c_0}_1 - {\mathcal{R}}^Y_1\) and those centers \(c \in {\mathcal{R}}^{Y+c_0}_1 \cap \mathcal{R}^Y_1\) with \(\text{cost-inc}_c(Y \cup \{c_0\}) < \text{cost-inc}_c(Y)\). Notice that we may assume that, for each \(c\in \Delta_{\text{new}}\), \[\text{cost-inc}_c(Y \cup \{c_0\})= \text{cost-inc}_c(\{c_0\}) = \begin{cases} 0 & if C_c is hit by c_0, \\ \sum_{p\in C^{\text{core}}_c} (d^2(p, c_0) - d^2(p,\{c\} \cup \Lambda)) & otherwise. \end{cases}\] because the addition of \(c_0\) caused a decrease in \(\text{cost-inc}(\cdot)\) for these clusters. Further, let \(\Delta_{\text{old}}^Y\) contain \(\Delta_{\text{new}} \cap \mathcal{R}^{Y}_1\) and the centers \(\mathcal{R}^{Y}_1 - \mathcal{R}^{Y+c_0}_1\). So \(\mathcal{R}^{Y+c_0}_1 = (\mathcal{R}^{Y}_1 - \Delta_{\text{old}}^Y) \cup \Delta_{\text{new}}\) and \(|\Delta_{\text{new}}| = |\Delta_{\text{old}}^Y|\). We now define \(\mathcal{R}^{X+c_0}_1 = (\mathcal{R}^{X}_1 - \Delta_{\text{old}}^X) \cup \Delta_{\text{new}}\) where \(\Delta^X_{\text{old}}\) is obtained as follows:

• Initialize \(\Delta_{\text{old}}^X =\mathcal{R}^X_1\cap \Delta_{\text{new}}\).

• While \(|\Delta_{\text{old}}^X|< |\Delta_{\text{new}}|\), add a center according to the following priorities:

  • If there is a center \(c\in \mathcal{R}^X_1 - \Delta_{\text{old}}^X\) so that \(c\not \in \mathcal{R}^{Y}_1\), add \(c\) to \(\Delta_{\text{old}}^X\).

  • Else add a center \(c\in \mathcal{R}^X_1- \Delta_{\text{old}}^X\) with \(c \in \Delta_{\text{old}}^Y\).

We claim that the above is well-defined, i.e., that if \(|\Delta_{\text{old}}^X|< |\Delta_{\text{new}}|\) and there is no center \(c\in \mathcal{R}^X_1 - \Delta_{\text{old}}^X\) so that \(c\not \in \mathcal{R}^{Y}_1\), then there must be a center \(c\in \mathcal{R}^X_1- \Delta_{\text{old}}^X\) with \(c \in \Delta_{\text{old}}^Y\). This is because \(|\Delta_{\text{old}}^X| < |\Delta_{\text{new}}| = |\Delta_{\text{old}}^Y|\) and all remaining centers \(\mathcal{R}_1^X - \Delta_{\text{old}}^X\) are in \(\mathcal{R}^Y_1\) in this case, so at least one of them must be in \(\Delta_{\text{old}}^Y\) (recall that \(|\mathcal{R}_1^X| = |\mathcal{R}_1^Y|\)).

With this notation we have \[\begin{align} \notag \sum_{c\in {\mathcal{R}}^{Y+c_0}_1} \text{cost-inc}_c(Y \cup \{c_0\}) - \sum_{c\in {\mathcal{R}}^Y_1} \text{cost-inc}_c(Y ) & = \sum_{c\in \Delta_{\text{new}}} \text{cost-inc}_c(Y \cup \{c_0\}) - \sum_{c\in \Delta_{\text{old}}^Y}\text{cost-inc}_{c}(Y) \\ &= \sum_{c\in \Delta_{\text{new}}} \text{cost-inc}_c(\{c_0\}) - \sum_{c\in \Delta_{\text{old}}^Y}\text{cost-inc}_{c}(Y)\,. \label{eq:firstsubmodular} \end{align}\tag{38}\] Similarly, \[\begin{align} \sum_{c\in {\mathcal{R}}^{X+c_0}_1} \text{cost-inc}_c(X \cup \{c_0\}) - \sum_{c\in {\mathcal{R}}^X_1} \text{cost-inc}_c(X ) & \leq \sum_{c\in \Delta_{\text{new}}} \text{cost-inc}_c(X \cup \{c_0\}) - \sum_{c\in \Delta_{\text{old}}^X} \text{cost-inc}_{c}(X) \notag\\ & = \sum_{c\in \Delta_{\text{new}}} \text{cost-inc}_c(\{c_0\}) - \sum_{c\in \Delta_{\text{old}}^X} \text{cost-inc}_{c}(X)\,, \label{eq:secondsubmoular} \end{align}\tag{39}\] where we used that \(\text{cost-inc}_c(X \cup \{c_0\}) \leq \text{cost-inc}_c(X)\) for any center \(c\).

To analyze this let, \(\pi\) be a bijection from \(\Delta_{\text{old}}^Y\) to \(\Delta_{\text{old}}^X\) so that \(\pi(c) =c\) for all \(c\in \Delta_{\text{old}}^Y \cap \Delta_{\text{old}}^X\).

By the above claim, \[\begin{align} \sum_{c\in \Delta_{\text{old}}^Y}\text{cost-inc}_{c}(Y) - \sum_{c\in \Delta_{\text{old}}^X} \text{cost-inc}_{c}(X) = \sum_{c\in \Delta_{\text{old}}^Y}(\text{cost-inc}_{c}(Y) - \text{cost-inc}_{\pi(c)}(X)) \leq 0 \end{align}\] and thus by the above arguments (see 38 and 39 ) \[\sum_{c\in {\mathcal{R}}^{X+c_0}_1} \text{cost-inc}_c(X \cup \{c_0\}) - \sum_{c\in {\mathcal{R}}^X_1} \text{cost-inc}_c(X)\leq \sum_{c\in {\mathcal{R}}^{Y+c_0}_1} \text{cost-inc}_c(Y \cup \{c_0\}) - \sum_{c\in {\mathcal{R}}^Y_1} \text{cost-inc}_c(Y)\,.\] We have thus verified 37 , which concludes the proof that \(f\) is submodular.

6.10.1.3 Proof of Property 3.

We have \[\begin{align} g(\emptyset) &\leq \sum_{c \not \in \mathcal{R}}\sum_{p\in C_c} d^2\left(p, \{c\} \cup \Lambda\right) + \sum_{p\in D_0} d^2(p, \Lambda) + \sum_{c\in \mathcal{R}_1} \text{closedcost}_c(\emptyset) \\ &\leq \sum_{c \not \in \mathcal{R}_1}\sum_{p\in C_c} d^2\left(p, \{c\} \cup \Lambda\right)+ \sum_{p\in D_0} d^2(p, \Lambda) + \sum_{c\in \mathcal{R}_1} \min_{c'\in \Lambda} \sum_{p\in C_c} d^2(p, c') \\ & \leq \sum_{c \not \in \mathcal{R}}\sum_{p\in C_c} d^2\left(p, \{c\} \cup \Lambda\right) + \sum_{c\in \mathcal{R}} \min_{c'\in \Lambda} \sum_{p\in C_c} d^2(p, c')\,\overset{\eqref{eq:MprimeDbound}}{\leq}\text{cost}(M'_D, \mu'_D). \end{align}\]

6.10.1.4 Proof of Property 4.

Similarly to the third property, we have \[\begin{align} g(\hat{\text{OPT}}) &\leq \sum_{c \not \in \mathcal{R}}\sum_{p\in C_c} d^2\left(p, \{c\} \cup \hat{\text{OPT}}\cup \Lambda\right) + \sum_{p\in D_0}d^2(p, \hat{\text{OPT}}\cup \Lambda)+ \sum_{c\in \mathcal{R}_1} \text{closedcost}_c(\hat{\text{OPT}}) \\ &\leq \sum_{c \not \in \mathcal{R}}\sum_{p\in C_c} d^2\left(p, \{c\} \cup \hat{\text{OPT}}\cup \Lambda\right)+ \sum_{p\in D_0} d^2(p, \hat{\text{OPT}}\cup \Lambda) + \sum_{c\in \mathcal{R}_1} \min_{c'\in \hat{\text{OPT}}} \sum_{p\in C_c} d^2(p, c') \\ & \leq \sum_{c \not \in \mathcal{R}}\sum_{p\in C_c} d^2\left(p, \{c\} \cup \hat{\text{OPT}}\cup \Lambda\right) + \sum_{c\in \mathcal{R}} \min_{c'\in \hat{\text{OPT}}} \sum_{p\in C_c} d^2(p, c')\overset{\eqref{eq:MprimeObound}}{\leq}\,\text{cost}(M'_O, \mu'_O)\,. \end{align}\]

6.10.1.5 Proof of Property 5.

Let \(X \subseteq F_\mathcal{B}\) be independent in \(M_\mathcal{B}\), i.e., \(X\in \mathcal{I}\). We give a polynomial-time algorithm that outputs \(k\) centers \(S^*\) whose cost is at most \((1+{6\sqrt{\varepsilon}})g(X) +{O(\sqrt{\varepsilon})} \text{opt}\). Let \(\mathcal{R}'_1\) be a subset of \(\mathcal{P}_1\) of cardinality \(|\mathcal{R}_1|\) so that \[\begin{gather} g(X) = \sum_{c \not \in \mathcal{R}_0 \cup \mathcal{R}'_1}\sum_{p\in C_c} d^2\left(p, \{c\} \cup X \cup \Lambda\right)+ \sum_{p \in D_0} d^2(p, X \cup \Lambda) + \sum_{c\in \mathcal{R}'_1} \text{closedcost}_c(X) \end{gather}\] By the arguments in the proof of the first property, we can calculate \(\mathcal{R}'_1\) in polynomial time. Now we obtain \(S^*\) from \(S_{\mathcal{Q}\cup \widetilde{\mathcal{U}}}\) as follows:

• Initialize \(S^* = S_{\mathcal{Q}\cup \widetilde{\mathcal{U}}} = S \cup \Lambda- \mathcal{Q}- \widetilde{\mathcal{U}}\).

• Remove centers \(\mathcal{R}_0\) and \(\mathcal{R}'_1\) from \(S^*\) to obtain a set of cardinality \(k\).

• Finally, for each ball in \(\mathcal{B}\), if \(X\) contains a center \(c\) in that ball, replace the associated dummy center by \(c\); otherwise, replace the dummy center with an arbitrary center in the ball.

By the definition of the dummy centers \(\Lambda\), the cost of \(S^*\) is at most \[\sum_{c \not \in \mathcal{R}_0 \cup \mathcal{R}'_1}\sum_{p\in C_c} d^2\left(p, \{c\} \cup X \cup \Lambda\right)+ \sum_{p \in D_0} d^2(p, X \cup \Lambda) + \sum_{c\in \mathcal{R}'_1} \sum_{p\in C_c} d^2(p, X \cup \Lambda)\,.\] To relate this to \(g(X)\) we need to compare \(\text{closedcost}_c(X)\) to \(\sum_{p\in C_c} d^2(p, X\cup \Lambda)\) for every \(c\in \mathcal{R}'_1\). Recall that \[\text{closedcost}_c(X) = \begin{cases} \sum_{p\in C_c} d^2(p, \{c\} \cup X \cup \Lambda) & if C_c is hit by X, \\ \min_{c'\in X \cup \Lambda} \sum_{p\in C^{\text{core}}_c} d^2(p, c') + \sum_{p\in C_c - C^{\text{core}}_c} d^2(p, \{c\} \cup X \cup \Lambda) & otherwise. \end{cases}\] Consider a cluster \(C_c\) with \(c\in \mathcal{R}'_1\) and partition \(C_c- C^{\text{core}}_c\) into sets \(A\) and \(B\) where \(A\) contain those clients \(p\) so that \(d^2(p, \{c\} \cup X \cup \Lambda) =d^2(p, c)\) and \(B\) contain the remaining clients \(p\) with \(d^2(p, \{c\} \cup X \cup \Lambda) = d^2(p, X \cup \Lambda)\). Suppose first that there is a client \(p_0 \in C^{\text{core}}_c\) so that \(d^2(p_0, X \cup \Lambda) < d^2(p_0, c)\), i.e., \(C_c\) is hit by \(X \cup \Lambda\). By the definition of \(C^{\text{core}}_c\), we have \(d^2(p_0, c) \leq \varepsilon\cdot \frac{\text{cost}(S)}{|\mathcal{R}| \cdot |C_c|}\leq {25} \varepsilon\cdot \frac{\text{opt}}{|\mathcal{R}| \cdot |C_c|}\) where we used that \(\text{cost}(S) \leq {25} \cdot \text{opt}\) by 11. So, by 3 with \(\gamma=1+\sqrt{\varepsilon}\), \[\begin{align} \sum_{p\in C_c} d^2(p, X \cup \Lambda)& \leq \sum_{p\in C_c - B} ({1+\sqrt{\varepsilon}})d^2(p, c) + {\frac{3}{\sqrt{\varepsilon}}}(d^2(c, p_0) + d^2(p_0, X\cup \Lambda))) + \sum_{p\in B} d^2(p, X \cup \Lambda)\\ &\leq \sum_{p\in A} ({1+\sqrt{\varepsilon}})d^2(p, c) + \sum_{p\in B} d^2(p, X \cup \Lambda) + {O(\sqrt{\varepsilon})} \cdot \frac{\text{opt}}{|\mathcal{R}|}\,, \end{align}\] where we used \(d^2(p_0, X\cup \Lambda) \leq d^2(p_0,c) \leq {25}\varepsilon\cdot \frac{\text{opt}}{|\mathcal{R}| \cdot |C_c|}\) and \(d^2(p, c) \leq {25}\varepsilon\cdot \frac{\text{opt}}{|\mathcal{R}| \cdot |C_c|}\) for any \(p\in C^{\text{core}}_c\) for the last inequality. Moreover, by definition of \(A\) and \(B\), we have \(d^2(p, c) = d^2(p, \{c\} \cup X \cup \Lambda)\) for any \(p\in A\) and \(d^2(p, X \cup \Lambda) = d^2(p, \{c\} \cup X \cup \Lambda)\) for any \(p\in B\). Hence, in this case, \[\begin{align} \sum_{p\in C_c} d^2(p, X \cup \Lambda) & \leq \sum_{p\in C_c - C^{\text{core}}_c} ({1+\sqrt{\varepsilon}})d^2(p, \{c\} \cup X \cup Y) + {O(\sqrt{\varepsilon})} \cdot \frac{\text{opt}}{|\mathcal{R}|} \\ & \leq ({1+\sqrt{\varepsilon}})\text{closedcost}_c(X) + {O(\sqrt{\varepsilon})} \cdot \frac{\text{opt}}{|\mathcal{R}|}\,. \end{align}\]

Let us now consider a cluster \(C_c\) with \(c\in \mathcal{R}'_1\) that is not hit by \(X \cup \Lambda\), particularly not by \(X\). Then \[\begin{align} \text{closedcost}_c(X) &= \min_{c'\in X \cup \Lambda} \sum_{p\in C^{\text{core}}_c} d^2(p, c') + \sum_{p\in C_c - C^{\text{core}}_c} d^2(p, \{c\} \cup X \cup \Lambda) \\ & \geq \sum_{p\in C^{\text{core}}_c \cup B} d^2(p, X \cup \Lambda) + \sum_{q\in A} d^2(q, c) \end{align}\] Using again 3 with \(\gamma=1+\sqrt{\varepsilon}\), we further have \[\begin{align} & \sum_{p\in C_c} d^2(p, X \cup \Lambda) = \sum_{p\in C^{\text{core}}_c \cup B} d^2(p, X \cup \Lambda) + \sum_{q\in A} d^2(q, X \cup \Lambda)\\ \leq & \sum_{p\in C^{\text{core}}_c \cup B} d^2(p, X \cup \Lambda) + \sum_{q\in A} \frac{1}{|C^{\text{core}}_c|}\sum_{p\in C^{\text{core}}_c} ({(1+\sqrt{\varepsilon})}d^2(q, c) + {\frac{3}{\sqrt{\varepsilon}}}(d^2(c, p) + d^2(p, X \cup \Lambda))). \end{align}\] Now using that \(d^2(c,p) \leq {25} \varepsilon\cdot \frac{\text{opt}}{|\mathcal{R}| \cdot |C_c|}\) for \(p\in C^{\text{core}}_c\) we obtain the upper bound \[\begin{align} \sum_{p\in C_c} d^2(p, X \cup \Lambda) & \leq \left(1 + {\frac{3}{\sqrt{\varepsilon}}}\frac{|A|}{|C^{\text{core}}_c|}\right) \sum_{p\in C^{\text{core}}_c\cup B} d(p, X \cup \Lambda) + \sum_{q\in A} ({1+\sqrt{\varepsilon}})d^2(q,c) + {O(\sqrt{\varepsilon})} \frac{\text{opt}}{|\mathcal{R}|} \\ &\leq \left(1 + \frac{{3\sqrt{\varepsilon}}}{(1-\varepsilon)}\right) \sum_{p\in C^{\text{core}}_c \cup B} d^2(p, X \cup \Lambda) + \sum_{q\in A} ({1+\sqrt{\varepsilon}})d^2(q,c) + {O(\sqrt{\varepsilon})} \frac{\text{opt}}{|\mathcal{R}|} \\ & \leq (1+ {6\sqrt{\varepsilon}}) \text{closedcost}_c(X) + {O(\sqrt{\varepsilon})} \frac{\text{opt}}{|\mathcal{R}|} \end{align}\] where for the second inequality we used that every cluster in \(\mathcal{R}_1' \subseteq \mathcal{P}_1\) is concentrated. In more detail \(|C^{\text{core}}_c|\geq (1-\varepsilon)|C_c|\), hence \(|A|\leq |C_c-C^{\text{core}}_c|\leq \varepsilon|C_c|\).

Hence, for any \(c\in \mathcal{R}_1'\) we proved that \(\sum_{p\in C_c} d^2(p,X \cup \Lambda) \leq (1+{6\sqrt{\varepsilon}}) \text{closedcost}_c(X) + {O(\sqrt{\varepsilon})} \frac{\text{opt}}{|\mathcal{R}|}\). We thus get that the cost of \(S^*\) is at most

\[\begin{align} & \sum_{c \not \in \mathcal{R}_0 \cup \mathcal{R}'_1}\sum_{p\in C_c} d^2\left(p, \{c\} \cup X \cup \Lambda\right)+ \sum_{p \in D_0} d^2(p, X \cup \Lambda) + \sum_{c\in \mathcal{R}'_1}\left( (1+{6\sqrt{\varepsilon}}) \text{closedcost}_c(X) + {O(\sqrt{\varepsilon})} \frac{\text{opt}}{|\mathcal{R}|}\right) \\ \leq & (1+{6\sqrt{\varepsilon}}) \left(\sum_{c \not \in \mathcal{R}_0 \cup \mathcal{R}'_1}\sum_{p\in C_c} d^2\left(p, \{c\} \cup X \cup \Lambda\right)+ \sum_{p \in D_0} d^2(p, X \cup \Lambda) + \sum_{c\in \mathcal{R}'_1} \text{closedcost}_c(X) \right) + {O(\sqrt{\varepsilon})} \text{opt} \end{align}\] which equals \((1+{6\sqrt{\varepsilon}}) g(X) + {O(\sqrt{\varepsilon})} \text{opt}\).

6.11 Putting Everything Together – Proof of [thm:mainadditivecenters]↩︎

We now turn to proving [thm:mainadditivecenters]. We first describe the algorithm we analyse.

6.11.0.1 A \(({5+O(\sqrt{\varepsilon})})\)-Approximation for \(k\)-Means on \(\varepsilon/\log n\)- Stable Instances \(k, D, F, d\).

1. \(S \gets\) Local Search on \((k, D, F, d)\), see 6.1

2. Let \(W\) be the set produced by the \(s^*\)-\(D^2\)-Sample process (6.3) on \(S\).

3. For each \(\mathcal{B}{\in \mathcal{L}_{\texttt{bal}}}\) produced as in 17 on \(W\) (6.4)

   1. Compute the set of dummy centers \(\Lambda\).

   2. For each \(Q\in \mathcal{L}_{\text{exp}}\) produced by \(expRem()\) as in 20 (6.6):

      1. For each \(\widetilde{\mathcal{U}}{\in \mathcal{L}_{cheap}}\), with the associated \(\widetilde{\mu}\), produced by \(cheapRem()\) as in 22 (6.8):
         solve the Submodular Instance associated with \((W,\mathcal{B},\mathcal{Q},\widetilde{\mathcal{U}},\widetilde{\mu})\) and compute the resulting solution \(S_{(W,\mathcal{B},\mathcal{Q},\widetilde{\mathcal{U}},\widetilde{\mu})}\) for the \(k\)-Means problem (6.10)

4. Output the solution \(S_{(W,\mathcal{B},\mathcal{Q},\widetilde{\mathcal{U}},\widetilde{\mu})}\) of smallest cost.

Our approach consists of running the above algorithm \(\log n\) times and take the minimum cost solution output. In the remaining, we argue that the above algorithm yields a \({(5+O(\sqrt{\varepsilon}))}\)-approximation with probability at least \((1-\varepsilon)(1-1/n)\). [thm:mainadditivecenters] then follows as the probability that all \(\log n\) executions of the above algorithm would fail is at most \(\left( 1- (1-\varepsilon)(1-1/n)\right)^{\log n} < 1/n\), and so with high probability we output a \({(5+O(\sqrt{\varepsilon}))}\)-approximation.

We first discuss the algorithm’s running time. The local search algorithm runs in polynomial time. The submodular optimization step also runs in polynomial time by 23. The output of the \(s^*\)-\(D^2\)-Sample process (which is clearly a polynomial time procedure) is a subset of size \(s^*\). By 17, we construct \(\mathcal{L}_{\texttt{bal}}\) in polynomial time the and \(|\mathcal{L}_{\texttt{bal}}| \leq n^{\varepsilon^{-O(1)}}\). Similarly, by 20, we construct \(\mathcal{L}_{\text{exp}}\) in polynomial time and \(|\mathcal{L}_{\text{exp}}| \leq n^{\varepsilon^{-O(1)}}\). Finally, by 22, we construct in polynomial-time \(n^{\varepsilon^{-O(1)}}\) candidate pairs \(\widetilde{\mathcal{U}},\widetilde{\mu}\). It follows that we solve \(n^{\varepsilon^{-O(1)}}\) submodular function optimization problems, each in polynomial time, and thus the total running time is polynomial.

We turn to proving the approximation guarantee. We aim to show that there exists a \((W,\mathcal{B}, \mathcal{Q}, \widetilde{\mathcal{U}}, \widetilde{\mu})\) that is successful with probability at least \((1-\varepsilon/2)\). Assuming this, the \(({5+O(\sqrt{\varepsilon})})\)-approximation follows by combining [claim:Mprimebounds] and 23.

Our algorithm runs the local search algorithm and thus provides us with a \(25\) approximation \(S\) that satisfies 13. Our algorithm next applies the \(s^*\)-\(D^2\)-Sample process and by 14 finds a successful \(W\) with probability at least \(1-\varepsilon\). Condition on having a successful \(W\), we then have that 15 holds, and so does 16 by combining 13 and 15.

Next, 17 implies that one set of balls \(\mathcal{B}\) is a valid set of balls. This set of balls induces a set of dummy centers \(\Lambda\). From there, by 20 we have that with probability at least \(1-1/n\), the \(expRem()\) procedure produces a set of center \(\mathcal{Q}\) such that

1. \(\mathcal{Q}\subseteq S_0\), and

2. The total cost in solution \(S {-} Q \cup \Lambda\) of the clusters in \(S_0 {-} \mathcal{Q}\) is at most \(\varepsilon\cdot \text{opt}\).

and so we have obtained a successful \((W,\mathcal{B},\mathcal{Q})\) and 21 applies. Condition on the event that we obtain a successful \((W, \mathcal{B}, \mathcal{Q})\), then, applying 22 we have that the \(cheapRem()\) procedure outputs a pair \((\widetilde{\mathcal{U}}, \widetilde{\mu})\) satisfying the properties of 22. We thus have a \((W, \mathcal{B}, \mathcal{Q}, \widetilde{\mathcal{U}},\widetilde{\mu})\) such that \(W\) is successful, the set of balls \(\mathcal{B}\) from 6.4 is valid, \(\mathcal{Q}\) selected in 6.6 satisfies the properties of 20, and \(\widetilde{\mathcal{U}},\widetilde{\mu}\) selected in this section satisfies the properties of 22 and so \((W, \mathcal{B}, \mathcal{Q}, \widetilde{\mathcal{U}},\widetilde{\mu})\) is successful. Finally, note that the only probabilistic steps were the sampling of \(W\) and of \(\mathcal{Q}\), and both steps are successful with probability at least \((1-\varepsilon) \cdot (1-1/n)\), as required.

References↩︎

---

1. R. Gao was supported by a Stanford Graduate Fellowship. F. Grandoni was partially supported by the SNF Grants 200021-200731 and 200021-236706. E. Lee was partially supported by the NSF Grant 2236669. E. van. Wijland. was partially supported by the French PEPR integrated projects EPIQ (ANR-22-PETQ-0007).↩︎

2. It can equivalently be phrased as the problem of partitioning the given points (\(D\)) into \(k\) disjoint clusters and finding the best center (amongst \(F\)) for each cluster, to minimize the \(k\)-Means objective.↩︎

3. We let \([a]^+:=\max\{0,a\}\).↩︎

4. To be precise, this version of JMMSV establishes an equivalent notion of approximation, but not exactly the LMP 2-approximation defined in 1 ; the latter is achieved when we run the algorithm with scaled opening cost \(\hat{f} = 2f\), which our actual algorithms do with different scaling factors. In this overview, let us ignore this subtlety and conflate these two notions of approximations while just using \(f\).↩︎

5. An exception is the approach in [2] which loses only a factor \((1+\varepsilon)\). However their approach is tailored to the specific LMP algorithms considered in that work and hard to adapt to others, as the authors openly admit.↩︎

6. As standard for primal-dual-based algorithms, we describe the algorithm in a more intuitive continuous version. It is easy to convert it into an equivalent discrete algorithm.↩︎

7. We remark that we could select \(\zeta\) to be arbitrarily small and get a better guarantee in the statement of 23. We have chosen this value to simplify the calculations as improving the guarantee of 23 does not improve the overall result.↩︎