1 Introduction↩︎

Clustering is one of the most fundamental tasks in data science with wide-ranging applications across various domains, including unsupervised learning, text and image analysis, recommender systems, and bioinformatics. (Metric) \(k\)-Means is one of the most well-studied clustering objectives. Given a collection of \(n\) clients \(D\), a collection \(F\) of potential center locations in a metric space defined by distances \(d(.,.)\), and an integer \(k \in [n]\), the goal is to select \(k\) centers \(S \subseteq F\) such that \(\mathrm{conn}(S) := \sum_{j \in D} d^2(j,S)\) is minimized, where \(d(j,S)\) is the distance of point \(j\) to its nearest center in \(S\). A very important special case is Euclidean \(k\)-Means, where \(D\) are points in a Euclidean space, and \(F\) is implicitly given by any point in the Euclidean space.¹ By standard techniques [1], [2] and losing only a \((1+\varepsilon)\) factor in the approximation ratio for any constant \(\varepsilon> 0\), we can restrict \(F\) to a polynomial-time-computable subset of points. W.l.o.g., we will also assume that \(D\subseteq F\).

From a theoretical standpoint, the study of clustering problems through the lens of (polynomial-time) approximation algorithms has been extremely successful, yielding new algorithms and analysis techniques. In particular, \(k\)-Means and the closely related \(k\)-Median (whose objective function \(\sum_{j \in D} d(j, S)\) is the sum of non-squared distances) are two problems that have received the most attention. For \(k\)-Median, there have been gradual developments over the last three decades using various algorithmic techniques, including primal LP rounding [3], [4], primal-dual [5], greedy algorithm with dual fitting analysis [6]–[8], local search [9], and bipoint rounding [10]–[12]. The current best approximation ratio is \((2 + \varepsilon)\), achieved by [13]. It is known that one cannot improve the approximation ratio beyond \(1 + 2/e \approx 1.73\) [14].

For \(k\)-Means, the picture has been much more elusive. While new algorithmic techniques for \(k\)-Median often translate to \(k\)-Means, it takes nontrivial technical effort, and the approximation ratios usually end up being much larger: this is because only approximate versions of triangle inequalities are available for squared metrics, which present technical challenges for some \(k\)-Median algorithms, most notably the ones based on dual fitting.

For Metric \(k\)-Means, [15] gave a \(25\)-approximation using local search, which was improved to a \((9 + \varepsilon)\) approximation by [16] via primal-dual. Recently, the approximation factor was improved to \(5.83\) [17] via a novel dual fitting algorithm. For Euclidean \(k\)-Means, [18] gave a \((9 + \varepsilon)\)-approximation using local search, followed by improvements using primal-dual methods [16], [19], [20] and most recently to \(4 + \varepsilon\) [21] via dual fitting. The hardness for the metric case is \(1 + 8/e \approx 3.94\), while for Euclidean metrics the hardness is \(1.17\) [22], and under a stronger conjecture called the Johnson Coverage Hypothesis (JCH), the hardness can be strengthened to \(1.73\) [23]. We remark that we focus on the high-dimensional Euclidean case here, while better algorithms are known for the low-dimensional one [24]–[26].

In this work, we give a new dual fitting based algorithm that gives significantly better approximation ratios for both Euclidean \(k\)-Means and Metric \(k\)-Means. Our algorithm leads to a \(3+\ln 2+\varepsilon\;(\approx 3.694)\) approximation for the Euclidean case and a \(4.9 + \varepsilon\) approximation for the metric case. Our result for Euclidean metrics is the first algorithm that gives an approximation ratio below the \(3.94\)-hardness barrier for the metric case, showing that the optimal approximation ratios for these two versions will be different assuming \(\mathbf{P} \neq \mathbf{NP}\). Such a separation between general metrics and high-dimensional Euclidean metrics, even after extensive recent work [20], [23], [27], has not been achieved for \(k\)-Median and Facility Location.

Our main technical contribution is an improved Lagrangian Multiplier Preserving (LMP) approximation for the following (Uncapacitated) Facility Location problem (with squared distances and uniform opening costs), which we simply call Facility Location throughout the paper: the input consists of a set of clients \(D\), a collection of centers \(F\), a (possibly Euclidean) metric \(d\) on \(D \cup F\), and a facility opening cost \(f\), and the goal is to select \(S \subseteq F\) that minimizes \(f |S| + \mathrm{conn}(S)\). In particular, \(S \subseteq F\) is called an LMP \(c\)-approximation if \(c \cdot f|S| + \mathrm{conn}(S) \leq c \cdot {\rm OPT}_{\mathrm{LP}}(f)\), where \({\rm OPT}_{\mathrm{LP}}(f)\) denotes the optimal cost of the standard LP relaxation for Facility Location defined in 3. Our algorithms yield LMP \(4.9\) and \((3+\ln 2)\)-approximations for Facility Location in general and Euclidean metrics, respectively.

Our improved LMP approximation is obtained via a refined dual fitting algorithm that tightly accounts for the cost incurred, while still being able to facilitate an improved one-cluster based dual feasibility analysis. We also introduce a novel spectral analysis for \(k\)-Means, departing from the usual factor-revealing linear programs. The spectral analysis provides a more accurate method to analyze the squared objective in \(k\)-Means compared to the previous approaches. More specifically, upper bounding the approximation factor reduces to upper bounding the largest eigenvalue of a family of matrices. The latter task is highly non-trivial, but it leads to substantially lower approximation factors compared to known methods.

This is then combined with the two-step framework recently used in [13], [17], [21] to obtain our approximation algorithm. For the first step, building on our LMP approximation, we use the idea of log-adaptive algorithms to obtain a good approximate solution for \(k\)-Means that opens \(k + O(\log n / \varepsilon^3)\) centers. Let \(\text{OPT}_k\) be the optimal cost for an instance of \(k\)-Means.

The second step of the framework studies stable instances. A \(k\)-Means instance is \(\beta\)-stable if \(\text{OPT}_{k-1} \geq (1 + \beta)\text{OPT}_k\), where \(\text{OPT}_{k-1}\) denotes the optimal cost of \((k-1)\)-Means with the same clients and centers. We show that a slightly more careful analysis of the same algorithms as in [17], [21] yields the following results.

These results improve the \((5+\varepsilon)\)-approximation for general metrics and the \((4+\varepsilon)\)-approximation for Euclidean metrics [17], [21]. Similar to previous results [13], [17], [21], given [thr:mainBicriteria,thr:mainMetricStable,thr:mainEuclideanStable], we can easily prove our main results [thm:main95euclidean] and [thm:main95metric]; given a \(k\)-Means instance, if it is not \(\Omega(\varepsilon^4/\log n)\)-stable as a \(k'\)-Means instance for each \(k' = k, k-1, \dots, k - O(\log n /\varepsilon^3)\), then \(\text{OPT}_{k - O(\log n /\varepsilon^3)} \approx \text{OPT}_k\) and we can use 1 with \(k \leftarrow k - O(\log n /\varepsilon^3)\). Otherwise, the instance is stable for some \(k'\) in the above range, and we can use [thr:mainMetricStable,thr:mainEuclideanStable]. The formal proof appears in 11.

1.0.0.1 Concurrent work.

Concurrently and independently, Byrka, Guo, Hu, Li, Wan, and Wang [28] gave an iterative randomized rounding framework for \(k\)-clustering with \(\ell_p^p\)-type assignment costs. For Metric \(k\)-Means, their framework gives a \((5+\varepsilon)\)-approximation, improving the previous \(5.83\)-approximation. For Euclidean \(k\)-Means, they give a \((4+\varepsilon)\)-approximation and prove an \(11/3\)-LMP approximation. Their approach is LP-rounding based and is conceptually different from our dual-fitting/spectral analysis. Our results have better approximation guarantees for both settings, giving \((4.9+\varepsilon)\) for Metric \(k\)-Means and \((3+\ln 2+\varepsilon)\) for Euclidean \(k\)-Means.

2 Technical Overview↩︎

We start by revisiting the dual fitting algorithm for Facility Location of the two papers by Charikar et al. [17], [21] for Metric and Euclidean \(k\)-Means, respectively. The two algorithms are very similar, except for one difference that [21] lowers certain dual values which can lead to further improvements. In this overview, this difference is not important, but we will use the stronger algorithm from [21] which is still as good as the older algorithm [17] for the metric case as well.

We provide a description of their algorithm below.

2.0.0.1 Greedy Algorithm.

The algorithm maintains a set \(A\) of active clients, a set \(S \subseteq F\) of open facilities, and a partition of clients into three classes: active clients \(A\), indirectly connected clients \(IC\), and directly connected clients \(DC\). Intuitively, active clients are not yet connected to an open facility, and they are actively growing their dual variables. Indirectly connected clients are connected to an open facility, but they do not directly contribute (a strictly positive amount) to the opening cost of that facility. On the other hand, directly connected clients are connected to an open facility and they directly contribute to the opening cost of that facility. Initially, \(A=D\), \(IC=DC=S=\emptyset\), and \(\alpha_j=0\) for every client \(j\in D\). Fix a parameter \(\rho > 1\). We let \(\hat{f} = \Gamma f\) where \(\Gamma\) is some constant that depends on \(\rho\) to be specified later, and \(f\) is the (uniform) facility cost.

As long as \(A\neq\emptyset\), the algorithm increases uniformly the (dual) values \(\alpha_j\) for all clients \(j\in A\) until one of the following two events occurs.

1. A facility becomes tight. There exists an unopened facility \(i\in F\setminus S\) such that \[\label{eqn:facold} \sum_{j\in A \cup IC}\left[\alpha_j - \rho \cdot d^2(j,i)\right]^+ + \sum_{j\in DC}\left[\rho d^2(j,S) - \rho d^2 (j, i)\right]^+ = \hat{f}.\tag{1}\] In this case, the algorithm opens facility \(i\). Every client \(j\in DC\) with \(d^2(j,S)>d^2(j,i)\) is reconnected to \(i\). For every client \(j\in A\cup IC\), the algorithm acts as follows:

   • if \(\alpha_j \ge \rho \cdot d^2(j,i)\), then \(j\) is connected (or reconnected) to \(i\) and moved to \(DC\);

   • else, if \(\alpha_j \ge d^2(j,i)\), then \(j\) is connected (or reconnected) to \(i\), its value \(\alpha_j\) is lowered to \(d^2(j,i)\), and \(j\) is moved to \(IC\) if it was in \(A\).

2. An active client reaches its current service cost. For some client \(j\in A\), we have \[\alpha_j = d^2(j,S).\] In this case, the client \(j\) is moved to \(IC\) and connected to its closest facility in \(S\).

Intuitively, the quantities of type \([\cdot]^+\) above are bids that client \(j\) is willing to pay to open facility \(i\). Via these bids, \(j\) can contribute to the opening of multiple facilities. At the same time, the residual dual \(\alpha_j\) of \(j\) after subtracting these bids is sufficient to connect \(j\) to the closest open facility. The algorithm naturally defines a notion of time, such that at time \(\theta\), all the active clients have \(\alpha\)-value \(\theta\).

Then, they show the following two cost accounting lemmas (similar to [7]) for this algorithm (considering the values of \(\alpha\) and \(S\) at the end of the algorithm), that lead to an LMP \(\Gamma\)-approximation.

Note that the RHS of this inequality is the cost of the computed solution (with facility cost \(\hat{f}=\Gamma\,f\)). The proof of this lemma follows from the fact that a facility is only opened when it is fully paid for (by the bids), and the fact that for any client \(j\), after subtracting any bid to open facilities, the remaining part of \(\alpha_j\) can still pay for its connection cost \(d^2(j,S)\).

2 implies that \((\alpha / \Gamma)\) is dual-feasible for the Facility Location LP (see 7 ) with opening cost \(f\). Combining this with 1, one gets \[\Gamma f |S| + \mathrm{conn}(S) \leq \sum_{j \in D} \alpha_j \leq \Gamma \cdot {\rm OPT}_{LP}(f),\] concluding that \(S\) is an LMP \(\Gamma\)-approximate solution.

While the dual payment lemma more or less follows directly from the algorithm description, the key challenge in [17], [21] has been proving dual feasibility given the complexity of the algorithm due to the presence of both IC and DC clients and the non-monotonicity of the dual variables.

2.0.0.2 The barrier of previous works.

Our first observation is that there is a barrier for the dual feasibility analysis of the algorithm of [21] for Euclidean metrics — a very simple example shows that the LMP \(4\)-approximation analysis of 2 can be tight for some facility \(i\). As illustrated in 1, we will consider points on a line (i.e., \(1\)-dimensional Euclidean space). Let \(i\) be a facility and \(D^*=\{1,2\}\) located at \(x_i = 0\), \(x_1=-1\), \(x_2=1\). Then, \(d^2(1,i)+d^2(2,i)=2\). Let \(\hat{f}=20\). The algorithm simultaneously increases \(\alpha_1\) and \(\alpha_2\) until \([\alpha_1-\rho d^2(1,i)]^+ +[\alpha_2-\rho d^2(2,i)]^+ = [\alpha_1-2]^+ +[\alpha_2-2]^+=20\), i.e., when \(\alpha_1=\alpha_2=12\). At this time point, suppose that client \(1\) gets directly connected to the only open facility at \(x_{i'}=-3\) (which was opened via bids from client \(1\) and other clients not in \(D^*\)). This leads to \(d^2(1,S)=4\). From now on, the bid of client \(1\) to facility \(i\) becomes \(8-2=6\). We will keep increasing \(\alpha_2\) until one of the following holds:

1. \(\alpha_2\geq d^2(2,S)=(x_2-x_{i'})^2 = 16\). In this case, we have client 2 indirectly connected to \(i'\).

2. \(6+[\alpha_2-2]^+=20\), which means \(\alpha_2=16\). In this case, we have \(i\) opened and client 2 directly connected to \(i\).

Therefore, we have \(\sum_{j\in D^*} \alpha_j = 12+16=28\). On the other hand, we have \(\hat{f} + 4(d^2(1,i)+d^2(2,i)) = 20 + 4\cdot 2 = 28\). Hence, the dual feasibility is tight for this example, and we cannot hope for an LMP approximation ratio below \(4\) using this algorithm.

2.0.0.3 An issue with the payments.

Our new algorithm is inspired by the fact that while the dual feasibility can be tight, the dual payment 1 does not tightly account for the costs: \(\sum_{j} \alpha_j\) may actually be larger than the cost of the algorithm. In particular, we observe that for any client that is directly connected at the end, the algorithm always overpays the connection cost by a factor of \(\rho\). To see this, let \(DC^{end}\) denote the set of clients which are directly connected at the end. Let \(j \in DC ^{end}\), and suppose through the execution of the algorithm, \(j\) was connected to facilities \(i_1, i_2...\ldots i_t\). For simplicity, assume that \(j\) was directly connected to each of these facilities. Then the total bid (dual contribution, see 1 ) of \(j\) to these facilities is at most \[(\alpha_j - \rho d^2(j,i_1)) + \sum_{\ell=1}^{t-1} \rho d^2(j,i_{\ell})-\rho d^2(j, i_{\ell+1}) = \alpha_j - \rho d^2(j,i_t)~.\] This means that after accounting for facility costs, \(\alpha_j\) pays for \(\rho d^2(j, i_t)\), where \(\rho > 1\), leading to an overpayment of the connection cost by the factor \(\rho\). Unfortunately, one cannot set \(\rho = 1\) as it would not give even an LMP \(O(1)\)-approximation, as already observed in [17]. This means we could simply subtract \((\rho - 1 )d^2(j,i_t)\) in the dual feasibility analysis.

However, while this overpayment is not ideal, it seems challenging to account for this naturally in the analysis: dual feasibility analysis is highly local, which cannot rely on the sequence of connections a client has: in particular, it is not clear how to meaningfully incorporate the subtraction \((\rho - 1)d^2(j,i_t)\) of client \(j\), which is a global quantity, into the local one-cluster dual feasibility analysis. The previous example of 1 illustrates this difficulty as well; when \(\alpha_2 = 16\), the events triggering the actions of “having client 2 indirectly connected to \(i'\)” and “opening \(i\) and connecting client 1 and 2 to \(i\)” happen simultaneously, but depending on how the algorithm breaks the ties (which can be easily manipulated by the position of the external facility \(i'\)), the amount of the overpayment substantially changes. In the former case, it is \((\rho - 1) \cdot d^2(1, i') = 4\) just from client 1, and in the latter case, it is \((\rho - 1) \cdot (d^2(1, i) + d^2(2,i)) = 2\).

2.0.0.4 Our new algorithm with tight dual payments.

To resolve this issue, we propose a new algorithm (2) that allows us to get the best of both worlds. Our new algorithm still has an overpayment, but we subtract this overpayment at the end, while still being able to incorporate this subtraction in the dual feasibility analysis. In our new algorithm, each client reserves a part of its first direct connection cost, that it does not offer to other facilities for subsequent (direct) connections. The key is that clients only move to facilities where their connection cost is lower, so the first direct connection cost is an upper bound for every future direct connection.

Notice that until the first direct connection, an active or IC client \(j\) offers \(\alpha_j - \rho d^2(j,i)\), but subsequently, instead of offering the natural \(\rho d^2(j,S) - \rho d^2(j,i)\), it only offers \(d^2(j,S) - d^2(j,i)\). Thus if \(R_j = d(j,i_1)\) is the first direct connection distance, then \((\rho - 1)R_j^2\) is reserved by the client and never used for any subsequent bid to facilities.

Following the same analysis as before, consider any client \(j \in DC^{end}\). If it goes through a sequence of direct connections \(i_1, i_2 ...\ldots i_t\), then it offers \[\alpha_j - \rho d^2(j,i_1) + \sum_{\ell=1}^{t-1} d^2(j,i_{\ell}) - d^2(j, i_{\ell+1}) = \alpha_j - (\rho - 1) d^2(j,i_1) - d^2(j,i_t)~.\] Since the connection cost at the end is only \(d^2(j,i_t)\), this leads to an overpayment of \((\rho - 1)d^2(j,i_1)\), where \(i_1\) is the first facility \(j\) was directly connected to. But now we define \(R_j = d(j,i_1)\), and subtract \((\rho - 1)R_j^2\) in the dual feasibility analysis to make this accounting exact (we will define \(R_j = 0\) if \(j\) is never directly connected). This turns out to be helpful in the dual feasibility analysis because we know \(R_j = d(j, i_1) \geq d(j,i_2) \geq ...\ldots \geq d(j, i_t)\), as clients only move to closer facilities as the algorithm progresses. \(R_j\) is thus an upper bound on all these distances, which facilitates a refined dual feasibility analysis, as discussed in greater detail in the following subsections.

Going back to the example from 1 with \(\rho = 2\), we set \(R_1 = d(1,i') = 2\) when \(i'\) is open at time 12, and the new algorithm indirectly connects client \(2\) to \(i'\) at time 16 (which now happens regardless of the tiebreaking rule as the bid of client \(1\) toward \(i\) becomes smaller), leading to \(R_2 = 0\). Thus we have \(\sum_{j \in D^*} \alpha_j = \alpha_1 + \alpha_2 = 28 = \hat{f} + 4\bigl(d^2(1,i) + d^2(2,i)\bigr)\) just like before, but because of the subtraction of \((\rho - 1)R_1^2 = R_1^2\), the dual feasibility is no longer tight. In particular we have \(\alpha_1 - R_1^2 + \alpha_2 - \hat{f} = 4 = 2\bigl(d^2(1,i) + d^2(2,i)\bigr)\), and consequently, \(\sum_{j \in D^*}(\alpha_j - 2d^2(j,i) - R_j^2)^+ \leq \hat{f}\), which improves the dual feasibility factor to \(2\) instead of \(4\) for this example.

2.0.0.5 Analysis.

From the above discussion, our modified payment and dual feasibility lemmas can be stated as follows. Let \(\alpha\) and \(S\) be the respective values at the end of the algorithm. The following dual payment lemma almost directly follows from the algorithm description and the discussion above.

The crux of our analysis is now to prove a dual feasibility lemma for a small enough constant \(\Gamma\) (depending on \(\rho\)). In particular, we prove the following.

Intuitively, the above dual feasibility lemma shows that the variables \(\alpha'_j:=\frac{1}{\Gamma}(\alpha_j-(\rho-1)R^2_j)\) induce a dual-feasible solution for the Facility Location LP. An LMP \(\Gamma\)-approximation follows by the same argument as before.

Furthermore, note that we can interpret our algorithmic idea as having smaller bids for clients who are directly connected. The intuition behind our new algorithmic idea (for the Euclidean case) is that these directly connected clients are “bad” in the dual feasibility analysis of [21] only when all clients and facilities are “almost” on a line (i.e., one-dimensional Euclidean space) and its maximum direct connection cost \(R_j^2\) is close to the squared distances between it and other clients \(j'\) with large dual variable \(\alpha_{j'}\) (compared to \(d^2(j',i)\)). As we have this additional \(-(\rho-1)R_j^2\) term in the final dual variable \(\alpha'_j\), we can make the final dual variables of these “bad” directly connected clients much smaller.

However, to apply this intuition to the analysis (for the general Euclidean cases), we confront several technical challenges:

• We need to deal with high-dimensional Euclidean spaces.

• We need to transform these \(-(\rho-1)R_j^2\) terms into \(-d^2(j,i)\) terms (because, intuitively, the dual feasibility constraint compares \(\sum_{j\in D'} \alpha_j-(\rho-1)R_j^2\) with \(f+\sum_{j\in D'} \Gamma d^2(j,i)\) for every \(D'\subseteq D\)).

• Finally, in the Euclidean space, [21] uses the pairwise distance \(\frac{1}{|D'|} \sum_{j<\ell\in D'} d^2(j,\ell)\) as an intermediate quantity, which is upper bounded by \(\sum_{j\in D'} d^2(j,i)\) (see, e.g., 62). However, this inequality may not hold if the coefficients before each \(d^2(j,i)\) have an average of 1 but are “unbalanced”. After we transform the \(-(\rho-1)R_j^2\) terms into the \(-d^2(j,i)\) terms, we observe that the coefficients of the new terms have an average \(\Omega(1)\) in the cases [21]’s analysis is tight. However, the coefficients are “unbalanced”, and some of them can even be zero (due to \(R_j=0\) for the corresponding client). To prove a ratio strictly better than 4, the most significant challenge is to show that the pairwise distance terms combined with these new \(-d^2(j,i)\) terms can be upper bounded by \((1-\Omega(1))\cdot \sum_{j\in D'} d^2(j,i)\).

As an illustration on how we resolve these challenges, we will provide a framework of our spectral analysis in 2.1, and provide simple proofs that obtain improved ratios (which are slightly worse than our final results) on the worst cases for the analysis of [17], [21] in [sec:euclidean-warmup] [sec:metric-warmup].

We remark that slightly smaller values of \(\Gamma\) (hence slightly better approximation factors)² can be obtained, for both the Euclidean and the metric case, with different values of \(\rho\) and a substantially more technical and longer analysis. We decided to omit this for the sake of clarity (and given that the presented analysis is already very complex). We also remark that our approach cannot yield significantly better approximation factors. Our numerical simulations suggests that our spectral analysis framework has a barrier of 3.51 for Euclidean metrics and 4.67 for general metrics. See the end of 2.1 for a brief discussion of this.

2.1 Spectral analysis: the general framework proving the dual feasibility↩︎

In this part, we will introduce the general idea of our new spectral analysis of the dual feasibility, which is the key to improving the approximations for both Euclidean and Metric \(k\)-Means. For simplicity, we will focus on the special case where at each point in time each considered client is active or directly connected (without ever being indirectly connected). This case is already non-trivial, and allows us to illustrate the main ideas. To simplify the notation, we also assume that the final dual values of the mentioned clients are all distinct. We start with some general considerations. Then, in 2.2 and 2.3 we will detail the analysis (for the mentioned special case) in the Euclidean and metric case, respectively. The analysis for arbitrary instances is given in [sec:greedy] [sec:euclidean-lmp] for the Euclidean case, and in [sec:greedy] [sec:euclidean-lmp] [sec:metric-lmp] for the metric case.

2.1.0.1 The general idea of the spectral analysis.

We will consider a similar set of constraints on the dual variables \(\alpha_j\) as in the original analysis of [7] as the starting point of our spectral analysis. Similar constraints are also used in the analysis of [17], [21].

Fix a facility \(i\), and let \(D^* = \{j \mid \alpha_j - (\rho-1)R_j^2 \geq \Gamma d^2(j,i)\}\). We need to show that \[\begin{align} \label{eqn:overview-dual-feasibility-simpler} \sum_{j \in D^*} \alpha_j - (\rho-1)R_j^2 \leq \hat{f} + \Gamma \sum_{j\in D^*}d^2(j,i)~. \end{align}\tag{2}\] Let us assume by renumbering that \(j \in D^*\) is sorted in increasing order of \(\alpha_j\), and let us denote \(D^* = [s]\). Recall that we assume that no two \(\alpha_j\) are equal.

Then we consider the contributions to \(i\) at time \(\alpha_t - \varepsilon\) for some arbitrarily small \(\varepsilon> 0\) for each \(t \in [s]\). Observe that each client \(j < t\) is already directly connected to an open facility in our special case. Let \(S^t\) be the set of open facilities at time \(\alpha_t - \varepsilon\), and \(r_{jt} = d(j, S^t)\) be the distance at which \(j\) is connected at time \(\alpha_t - \varepsilon\). Then, we get the following inequality. \[\label{eqn:facility} \sum_{j < t} (r_{jt}^2 - d^2(j,i)) + \sum_{j \geq t} (\alpha_t - \rho d^2(j,i)) \leq \hat{f}.\tag{3}\] This inequality follows from the fact that clients \(j < t\) are already (directly) connected at the time \(\alpha_t - \varepsilon\), and hence they bid \([r_{jt}^2 - d^2(j,i)]^+\) to \(i\), and clients \(j \ge t\) are still active, hence they bid \([\alpha - \rho d^2(j,i)]^+\) to \(i\), where \(\alpha\) is their current \(\alpha\) value (which at time \(\alpha_t - \varepsilon\), is \(\alpha_t - \varepsilon\)). We relax the terms of type \([x]^+\) to \(x\), which can only make the LHS smaller. Taking \(\varepsilon\to 0\) and since the algorithm ensures that at all times, every facility receives a total bid of at most \(\hat{f}\), implying 3 .

Next, we eliminate \(r_{jt}^2\) via triangle inequality. Since we know that at time \(\alpha_t - \varepsilon\), \(t\) was not connected, we must have \[\begin{align} \label{eqn:tri-overview} \alpha_t \leq d^2(t,S) \leq (d(j,t) + {r_{jt}})^2 \end{align}\tag{4}\] See 11 for a formal proof. To get started with the spectral analysis, we will plug this lower bound on \(r_{jt}\) into 3 and obtain a new inequality that upper bounds the sum of \(\alpha_j-(\rho-1)R_j^2\) in terms of \(\hat{f}\) and the sum of \(d^2(j,i)\) (see [lem:dual-feasibility-simpler-metric] [lem:dual-feasibility-simpler-euclidean] for the formal statements we used for general and Euclidean metrics respectively).

We will use a slightly different approach to use the inequality \(\alpha_t\leq (d(j,t)+r_{jt})^2\) for the metric case and the Euclidean case.

• Consider the Euclidean space of dimension \(\mathsf{dim}\). We will assume w.l.o.g. that the facility \(i\) is located at the origin of the Euclidean space and use the coordinate of the clients, \((x_j)_{j\in [s]}\in \mathbb{R}^{\mathsf{dim}}\), to rewrite \(d(j,t)=\|x_j-x_t\|\) and \(d(j,i)=\|x_j\|\). Using 3 4 , some quadratic inequalities and the Cauchy-Schwarz inequality (to make every \(d(\cdot,\cdot)\) appear in a quadratic form in the inequalities), we can upper bound \(\sum_{j} \alpha_j - (\rho - 1)R_j^2\) by the sum of \(\hat{f}\) and a quadratic form of \(x=(x_j)_{j\in [s]} \in \mathbb{R}^{s \cdot \mathsf{dim}},\) denoted as \(\hat{Q}(x)\). Because we have \(\sum_{j} d^2(j,i) = \sum_{j} \|x_j\|^2 = \| x \|_2^2\), it suffices to upper bound \(\frac{\hat{Q}(x)}{\|x\|_2^2}\). As an \((s \cdot \mathsf{dim}) \times (s \cdot \mathsf{dim})\) matrix, \(\hat{Q}\) is a block-diagonal matrix with \(\mathsf{dim}\) many identical \(s \times s\) matrices on the diagonal, so an upper bound on the maximum eigenvalue of this \(s \times s\) matrix yields a bound on the approximation ratio.

• In the metric case, we use the metric triangle inequality \(d(j,t)\leq d(j,i) + d(t,i)\). Then, we will consider each \(d(j,i)\) as a variable \(a_j\), and use 3 4 and some quadratic inequalities to upper bound \(\sum_{j} \alpha_j - (\rho - 1)R_j^2\) by a sum of \(\hat{f}\) and a quadratic form of \((a_j)_{j\in [s]}\), denoted as \(\hat{Q}(a)\). Similarly, because we have \(\sum_{j} d^2(j,i) = \sum_{j} a_j^2=\|a\|^2\), it suffices to upper bound the maximum eigenvalue of \(\hat{Q}(a)\) (again a quadratic form on \(s\) variables) to prove the bound on our objective 2 .

Numerically, we found the following barriers for our approach:

• Euclidean case: After the Cauchy-Schwarz step (which makes the spectral analysis work), the underlying matrix \(\hat{Q}(x)\) has an eigenvalue \(\geq 3.51\) in the case where there are 300 clients: the first 260 clients are immediately directly connected when they become inactive; while the last 40 clients are never directly connected.

• Metric case: The underlying matrix \(\hat{Q}(x)\) has an eigenvalue \(\geq 4.67\) when the number of clients is large. More details can be found at the end of 2.3.

2.2 Spectral analysis for the Euclidean case↩︎

Let us fix \(\rho=2\). Recall that we are considering the special case where all the considered clients at any point of time are either active or directly connected. To further simplify the analysis, we present a weaker bound \(\Gamma\leq 3+\frac{27}{32}<3.85\) for this case.

Let \(x_1, x_2 \ldots x_s\) be the vectors denoting the locations of the \(s\) clients, and suppose that the facility \(i\) is located at the all-zeros location, so that \(d^2(j,i) = \|x_j\|^2\) for each \(j \in [s]\). For simplicity, we will use \(x=(x_j)_{j\in [s]}\) to denote the \(s \cdot \mathsf{dim}\)-dimensional vector by concatenating the \(s\) vectors. In the Euclidean case, we will rewrite 4 by \(\alpha_t \leq (r_{jt} + ||x_j - x_t||)^2\) for each \(t > j\). Thus we get \(\alpha_t \leq r_{jt}^2 + ||x_j - x_t||^2 + 2r_{jt}||x_j - x_t||\). Adding across all \(t > j\) for a fixed \(j\) and using \(r_{jt} \leq R_j\), we get \[\begin{align} \sum_{t > j} (\alpha_t - r_{jt}^2) \leq \sum_{t > j} 2R_j||x_j - x_t|| + \sum_{t > j}||x_j - x_t||^2 \end{align}\] Subtracting \(s(\rho - 1)R_j^2=sR_j^2\), \[\begin{align} \sum_{t > j} (\alpha_t -{r_{jt}^2}) - {sR_j^2} \leq \sum_{t > j} \bigl(2R_j||x_j - x_t|| + ||x_j - x_t||^2\bigr) - {sR_j^2} \end{align}\]Applying the quadratic inequality \(2ay - by^2 \leq \frac{a^2}{b}\), obtained from \(\frac{1}{b}(a-by)^2\geq 0\), we have \[\begin{align} \label{eqn:intro-euclidean-no-overbidding} \sum_{t > j} (\alpha_t - {r_{jt}^2}) - {sR_j^2} \leq \frac{\Bigl( \sum_{t > j}||x_j - x_t||\Bigr)^2}{{s}} + \sum_{t > j} ||x_j - x_t||^2 \end{align}\tag{5}\]

We sum 3 over all \(t\in [s]\) and the above 5 over all \(1\leq j<t\leq s\), and then divide by \(s\). We get \[\sum_{t \in [s]} \alpha_t - R_t^2\leq \hat{f} + Q(x)~,\] where \(Q(x)\) is a function of \(x_1, x_2 \ldots x_s\) alone: \[Q(x) =\frac{1}{s}\left( \sum_j\frac{\bigl( \sum_{t > j}||x_j - x_t||\bigr)^2}{{s}} + \sum_j\sum_{t > j} ||x_j - x_t||^2 + \sum_{t}\sum_{j < t} ||x_j||^2 + \sum_{t}\sum_{j \ge t} \rho ||x_j||^2\right)\,.\] The first two terms come from the lower bound plugged in for \(r_{jt}^2\) from above, and the last two terms come from the minus terms on the LHS of 3 .

Observe that the second term \(\frac{1}{s} \sum_{j} \sum_{t>j} \|x_j-x_t\|^2\) equals to the optimal single-cluster cost for \([s]\) (see, e.g., 62). Hence, it is no greater than \(\sum_{j} \|x_j\|^2\). Rearranging the third and the fourth term via \(\frac{1}{s} \sum_{t} \sum_{j} \rho \|x_j\|^2 - \frac{1}{s} \sum_{t} \sum_{j<t} (\rho-1) \|x_j\|^2\), we can rewrite them as \(\rho\sum_{j} \|x_j\|^2 - (\rho-1) \sum_{j} \frac{s-j}{s} \|x_j\|^2\). Using our parameter \(\rho=2\), we have the following upper bound for \(Q(x)\): \[\begin{align} Q(x) \leq \frac{1}{s^2} \sum_{j} \left(\sum\nolimits_{t>j} \|x_j-x_t\|\right)^2 - \sum_j \frac{s-j}{s} \|x_j\|^2 + 3\sum_{j} \|x_j\|^2. \end{align}\]

Now we apply Cauchy-Schwarz on the first term to obtain \(Q(x)\leq \hat{Q}(x) + 3\sum_{j} \|x_j\|^2\), where \[\hat{Q}(x) := \sum_j \frac{s - j}{s^2} \sum_{t > j} ||x_j - x_t||^2 - \sum_{j} \frac{s-j}{s} \|x_j\|^2~.\]

If we show that \(\hat{Q}(x) \leq \sum_j \lambda||x_j||^2\), then we obtain an LMP \((3+\lambda)\)-approximation. This now reduces to upper bounding the largest eigenvalue of the matrix given by the quadratic form \(\hat{Q}\). We next prove that \(\lambda\leq \frac{27}{32}<0.85\).

Consider two cases for \(j\): \(j\geq s/4\) and \(j<s/4\). If \(j\geq s/4\), we have \(\frac{s-j}{{s^2}} \leq \frac{3}{4s}\). The total contribution of \(j\) to \(\hat{Q}(x)\) is then upper bounded by \(\frac{3}{4s} \sum_{t>j} \|x_j-x_t\|^2\). If \(j<s/4\), we have (in the second inequality, we use \((a-b)^2 \leq 3a^2+1.5b^2\)) \[\begin{align} \frac{s-j}{s^2} \sum_{t>j} \|x_j-x_t\|^2 - \frac{s-j}{s} \|x_j\|^2 &< \frac{1}{s} \sum_{t>j} \|x_j-x_t\|^2 - \frac{3}{4} \|x_j\|^2 \\ &\leq \frac{3}{4s} \sum_{t>j} \|x_j-x_t\|^2 + \frac{1}{4s} \sum_{t>j} (3\|x_j\|^2+1.5\|x_t\|^2)- \frac{3}{4} \|x_j\|^2 \\ &\leq \frac{3}{4s} \sum_{t>j} \|x_j-x_t\|^2 + \frac{3}{8s} \sum_{t>j} \|x_t\|^2 \end{align}\] The intuition behind discussing these two different cases is as follows: for small \(j\), there is a large subtraction term \(-\frac{s-j}{s} \|x_j\|^2\); because of this, and the flexibility of squared triangle inequalities, we can let \(\|x_j\|^2\) pay more for the cost in \(\sum_{t>j} \|x_j-x_t\|^2\) (and let \(\|x_t\|^2\) pay less).

Summing over all \(j\), using the upper bound discussed above, we have \[\begin{align} \hat{Q}(x) &\leq \frac{3}{4s} \sum_{j<t} \|x_j-x_t\|^2 + \frac{3}{8s}\sum_{j<s/4} \sum_{t>j} \|x_t\|^2 \\ &\leq \frac{3}{4s} \sum_{j<t} \|x_j-x_t\|^2 + \frac{3}{32} \sum_{t} \|x_t\|^2 \\ &\leq \frac{3}{4} \sum_{j} \|x_j\|^2 + \frac{3}{32} \sum_{t} \|x_t\|^2 = \frac{27}{32} \sum_{j} \|x_t\|^2~, \end{align}\] where the last inequality follows from the fact that \(\frac{1}{s}\sum_{j<t} \|x_j-x_t\|^2\) is the optimal single-cluster cost for \([s]\) among all possible centers in \(\mathbb{R}^{\mathsf{dim}}\) while \(\sum_{j}\|x_j\|^2\) is the cost of opening the center at the origin for \([s]\).

In 5, we generalize and refine the above derivation to the general case with both \(IC\) and \(DC\) clients to obtain the claimed LMP \((3+\ln 2)\)-approximation for Facility Location in the Euclidean case. Our proof combines the above derivation with the techniques introduced in [21] to deal with clients in \(IC\), and it is substantially more complex.

2.3 Spectral analysis for the metric case↩︎

Now we present an overview of our spectral analysis for Metric \(k\)-Means in the special case when all clients are either directly connected or active (i.e., \(IC\) is empty at all times). We generalize it to all possible cases in 6.

Let \(a_j= d(j,i)\). According to 4 and the metric triangle inequality, we have \(\alpha_t \leq d^2(t,S_t)\leq (r_{jt}+d(j,t))^2 \leq (r_{jt}+d(j,i)+d(t,i))^2 = (r_{jt}+a_j+a_t)^2\). Thus, we can expand the square to get the upper bound \(\alpha_t - r^2_{jt} \leq 2{r_{jt}}(a_j + a_t) + (a_j + a_t)^2\).

Fix any \(j\). Adding across all \(t > j\), \[\begin{align} \sum_{t > j} (\alpha_t - r^2_{jt}) &\leq \sum_{t > j} 2{r_{jt}}(a_j + a_t) + (a_j + a_t)^2 \\ &\leq 2{R_j}\Bigl((s-j)a_j + \sum_{t > j}a_t\Bigr) + \sum_{t > j} (a_j + a_t)^2, \end{align}\] where we crucially use the fact that each \(r_{jt} \leq R_j\) by the property of our algorithm. Now we subtract \(s(\rho - 1)R_j^2\) for \(j\), and obtain \[\begin{align} \sum_{t > j} (\alpha_t - r_{jt}^2) - s(\rho - 1) R_j^2 &\leq 2{R_j}\Bigl((s-j)a_j + \sum_{t > j}a_t\Bigr) - s(\rho - 1)R_j^2 + \sum_{t > j} (a_j + a_t)^2 \\ &\leq \frac{1}{s(\rho - 1)}\Bigl((s-j)a_j + \sum_{t > j}a_t\Bigr)^2 + \sum_{t > j} (a_j + a_t)^2 \end{align}\] where the last inequality is obtained by observing that the quadratic terms are controlled by \(2ax - bx^2 \leq \frac{a^2}{b}\), which is equivalent to \(\frac{1}{b}(a-bx)^2 \geq 0\).

If we now plug this lower bound for \(r_{jt}\), and follow the same procedure of adding each inequality at time \(\alpha_t - \varepsilon\) for \(t \in [s]\) and divide by \(s\), we would get a bound of the form \(\sum_j (\alpha_j- (\rho - 1)R_j^2) \leq \hat{f} + \hat{Q}(a)\) where \(\hat{Q}(a)\) is the new quadratic form on the distances arising from this process. If the maximum eigenvalue of the quadratic form corresponding to \(\hat{Q}(a)\) is \(\Gamma\), then this is at most \(\hat{f} + \Gamma a^Ta\), giving us an LMP \(\Gamma\)-approximation.

Now we can write the quadratic form \(\hat{Q} (a)\) corresponding to our setup. \[\begin{align} \hat{Q} (a) = \frac{1}{s} \left( \sum_{j} \frac{1}{s(\rho - 1)}\Bigl((s-j)a_j + \sum_{t > j} a_t\Bigr)^2 + \sum_{j}\sum_{t > j} (a_j + a_t)^2 + \sum_{t} \sum_{j \geq t} \rho a_j^2 + \sum_{t} \sum_{j < t} a_j^2 \right) \end{align}\]

The first and second terms come from substituting the lower bounds for \(r_{jt}\) from the calculation above. The third and fourth terms come from the form of the bids of the clients at time \(\alpha_t - \varepsilon\), clients \(j\geq t\) are active and bid \(\alpha_t - \rho a_j^2\) (hence the coefficient of \(\rho\) for \(a_j^2\)), whereas clients \(j < t\) are directly connected and bid \(r_{jt}^2 - a_j^2\) (hence the coefficient of \(1\) for \(a_j^2\)).

First, let us show via elementary analysis that when \(\rho = 2\), we can bound \(\Gamma \leq 5\). The first term inside the parenthesis, \(\sum_j\frac{1}{s}\bigl((s-j)a_j + \sum_{t > j} a_t\bigr)^2\), is equal to \(\sum_j\frac{1}{s}\left(\sum_{t > j} (a_j + a_t)\right)^2\). Using the Cauchy-Schwarz inequality, this is at most \(\sum_j\frac{s-j}{s}\sum_{t > j}(a_j + a_t)^2\), which is at most \(\sum_j\frac{2(s-j)}{s}\sum_{t > j} (a_j^2 + a_t^2)\). Again, for the middle term, we get \(\sum_j\sum_{t > j} (a_j + a_t)^2 \leq \sum_j 2\sum_{t > j} (a_j^2 + a_t^2) = \sum_{j}2\sum_{t > j} (a_j^2 + a_t^2)\).

Thus with \(\rho = 2\), the expression is upper bounded by

\[\begin{align} \hat{Q}(a) \leq \frac{1}{s}\left (\sum_{j} \left(2 + \frac{2(s-j)}{s}\right) \sum_{t > j}(a_j^2 + a_t^2) + \sum_t \sum_{j \geq t} 2a_j^2 + \sum_t \sum_{j < t}a_j^2\right) \end{align}\]

Now let us fix an index \(j^*\) and collect coefficients of \(a_{j^{*}}^2\) for the term inside the parenthesis. The total coefficient is then

\[2(s-j^*)\left(1 + \frac{s-j^*}{s}\right) + 2\sum_{j < j^*}\left(1 + \frac{s-j}{s}\right) + 2j^* + (s - j^*) = 5s - j^* - 4 + \frac{j^*(j^* + 1)}{s}\]

Dividing by the outer \(s\) and using \(j^* \leq s\), the total coefficient is at most \(5\), which gives the desired inequality \(\hat{Q}(a) \leq 5a^Ta\), and hence \(\Gamma \leq 5\).

However, it is possible to do better. By running experiments for various values of \(\rho\), one can observe that the largest eigenvalue of the symmetric matrix describing this quadratic form is (numerically) minimized when \(\rho \approx 2.5\). With this value of \(\rho\), we numerically estimate the top eigenvalue to be \(4.67783\) as \(s \to \infty\). It can be shown that the matrix corresponding to the quadratic form \(\hat{Q}(a)\) is the sum of a diagonal and a rank one matrix, and hence it is not difficult to prove an analytic bound that reaches close to this constant for any \(s\). Thus using our new algorithm with the new spectral analysis, we get an approximation ratio of \(4.68\) for this special case when no client is indirectly connected.

Thus all that is left is to generalize this to the case when there are indirectly connected clients. Let \(DC_t\) and \(IC_t\) be the sets of directly and indirectly connected clients at time \(\alpha_t - \varepsilon\). Then the number of such valid sets³ is exponential in \(s\). Thus, we need a way to bound the general case directly.

For this purpose, we write down the equations for the most general \(IC\)/\(DC\) case, and try to solve them all at once, similarly to [17]. However, this analysis becomes more technically challenging as the matrices involved are no longer simple sums of diagonal and rank one matrices. Nevertheless, we show that the quadratic form in the general case can always be bounded by 4.9. Our proof uses the 34 which says that given a non-negative matrix \(M\), if we find a non-zero vector \(h\) such that \(\max_i \frac{(Mh)_i}{h_i}\) is at most \(\Gamma\), then the top eigenvalue of \(M\) is at most \(\Gamma\).

2.4 Turning LMP approximations to bicriteria approximations↩︎

To turn the LMP approximations to bi-criteria approximations (opening \(k+O(\log n)\) centers), we employ the framework of [13], [17], [21]. The only technical improvement is as follows. The most recent work [21] can only turn an LMP \(\Gamma\)-approximation to a \(\Gamma+\varepsilon\) bi-criteria approximations for any constant \(\Gamma\geq 4\). We remove this condition, and now this step works for any constant \(\Gamma>\rho\) (which is a clear lower bound for our LMP algorithm).

3 Preliminaries↩︎

For a client \(j\) and a set of facilities \(S\), let \(d(j,S) := \min_{i\in S} d(j,i).\) For any real number \(x\), let \([x]_+ := \max\{x,0\}\). For a positive integer \(n\), we write \([n] := \{1,2,\dots,n\}\), and more generally, for integers \(l \le r\), we write \([l:r] := \{l,l+1,\dots,r\}\). We use \(\mathbb{R}_{>0}\) and \(\mathbb{R}_{\ge 0}\) to denote the sets of positive and non-negative reals, respectively.

3.0.0.1 Facility Location with uniform opening costs and squared distances.

In the Facility Location with uniform opening costs and squared distances problem, we are given a set of clients \(D\), a set of facilities \(F\), and a facility opening cost \(f > 0\). The goal is to find a set \(S \subseteq F\) minimizing \[f|S| + \sum_{j\in D} d^2(j,S).\] We write \({\rm OPT}_{\mathrm{FL}}(f)\) for the optimum value of this problem. We will refer to this problem simply as Facility Location throughout the paper.

3.0.0.2 LP relaxation.

The standard LP relaxation of Facility Location is \[\label{eq:lp-primal} \begin{align} \min \quad & f\sum_{i\in F} y_i + \sum_{j\in D}\sum_{i\in F} d^2(j,i)x_{ij} \\ \text{s.t.}\quad & \sum_{i\in F} x_{ij} \ge 1 && \forall j\in D,\\ & x_{ij} \le y_i && \forall j\in D,\;i\in F,\\ & x_{ij} \ge 0,\;y_i \ge 0. && \end{align}\tag{6}\] Its dual can be written as \[\label{eq:lp-dual} \begin{align} \max \quad & \sum_{j\in D} \alpha_j \\ \text{s.t.}\quad & \sum_{j\in D} [\alpha_j - d^2(j,i)]^+ \le f && \forall i\in F,\\ & \alpha_j \ge 0 && \forall j\in D. \end{align}\tag{7}\] We write \({\rm OPT}_{\mathrm{LP}}(f)\) for the optimum value of the 7 .

3.0.0.3 Bounding distances.

Thanks to standard reductions (see, e.g., [17]), in the metric case we will assume that distances are integers in \([1,n^3/\varepsilon]\) while losing a factor \(1+O(\varepsilon)\) in the approximation.

In the Euclidean case, we will assume that non-zero distances are in \([1,n^2/\varepsilon^2]\) thanks to the following reduction in [21].

3.0.0.4 Triangle Inequalities for Squared Distances.

We will use the following simple inequalities taken from [17]. 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 Our Improved Greedy Algorithm↩︎

In this section, we present our improved greedy LMP approximate algorithm for the Facility Location problem with uniform facility cost \(f\). As in the previous works [17], [21], our algorithm is based on the well-known greedy algorithm by [7].

In our algorithm, we use parameters \(\rho,\Gamma>1\), which will be defined (differently for the metric and Euclidean cases) later. We also use \(\hat{f} = \Gamma \cdot f\), where \(f\) is the uniform facility opening cost. The following 2 is our new greedy algorithm achieving the improved LMP approximation. As standard in the area, we present the algorithm in a continuous version where there is a variable \(\theta\) (interpreted as time) that grows continuously from \(0\) to \(+\infty\). Initially, all clients are active. We grow the \(\alpha\)-values of active clients at the same speed as \(\theta\) grows, until one of two conditions is met: (1) a new facility becomes openable, or (2) a client’s \(\alpha\)-value reaches its connection cost to some open facility. It is easy to discretize this algorithm. Here, we say one facility is openable if the total bid of the clients to this facility reaches \(\hat{f}\), where the bid of client \(j\) to a facility \(i\) is defined as \([\alpha_j-\rho\cdot d^2(j,i)]^+\) when \(j\) is in \(A\) or \(IC\), or as \([d^2(j,S)-d^2(j,i)]^+\) when \(j\) is in \(DC\), at the time \(i\) is opened. In addition, for each client \(j \in D\), we explicitly track its first direct-connection cost, \(R^2_j\). We set \(R_j = d(j, i^*)\) where \(i^*\) is the facility that caused \(j\) to be moved to \(DC\) for the first time. For a client that remains only indirectly connected or active, we keep \(R_j = 0\) as it is initialized.

Because we freeze \(\alpha_j\) when a client \(j\) becomes directly connected, we can show that the \(\alpha\) values are non-increasing after connection. For some intermediate time \(\theta\), we will use \(\alpha_j(\theta)\) to denote the \(\alpha\) value of client \(j\in D\) at that time. Further, for each \(j\in D\), we define \(\hat{\alpha}_j\) to be its \(\alpha\) value at the first time when \(j\) is connected to some facility. Let \((\alpha^*_j)_{j\in D}\) be the final vector of \(\alpha\)-values.

Our goal is to show that \((\alpha^*_j - (\rho-1)R_j^2)_{j\in D}\) can characterize the total cost of 2 and can be scaled to a feasible dual solution to the 7 . First, we prove that the sum of \((\alpha^*_j - (\rho-1)R_j^2)_{j\in D}\) is an upper bound of the total cost. The proof relies on the following observation: for each client \(j\) that is finally directly connected, its \(\alpha\)-value can be written as the sum of its final connection cost \(d^2(j,S)\), \(\rho-1\) times its first direct-connection cost \(R_j^2\), and the bids of \(j\) to the open facilities.

Second, we prove that \(\big((\alpha^*_j - (\rho-1)R_j^2)/\Gamma\big)_{j\in D}\) can be scaled to a dual feasible solution for an appropriate choice of \(\rho, \Gamma\). Since the proof of this part is very long, we will first present a general framework for both the metric case and Euclidean case in ¿sec:subsec:dual-feasibility-framework?, where we establish a set of constraints for the variables \((\alpha^*_j)_{j\in D}\) and \((R_j)_{j\in D}\). Then, we will use the constraints to establish the following lemma for the metric case and the Euclidean case separately in [sec:metric-lmp] [sec:euclidean-lmp], using significantly different techniques.

Dual feasibility: the analysis framework↩︎

Next, we show the analysis framework to establish 10. We will use \([n]\) to denote the set \(\{1,2,\dots, n\}\) and \([l:r]\) to denote the set \(\{l, l+1,\dots, r\}\). From now on, we fix any facility \(i\in F\). Let \(D^*=\{j\in D:\alpha_j^*-(\rho-1)R^2_j > \Gamma\cdot d^2(j,i)\}\). Suppose that \(|D^*|=s\) and we index them by \([s]\) and order them in the increasing order of the time they become inactive.

That is, we have \(\hat{\alpha}_1\leq \hat{\alpha}_2 \leq \dots \leq \hat{\alpha}_{s}\). Then, our goal (which will be formally shown at the end of this subsection) becomes proving \[\begin{align} \label{eqn:new-goal-dualf} \sum_{j\in [s]} (\alpha^*_j-(\rho-1)R_j^2) \leq \hat{f} + \Gamma \cdot \sum_{j\in [s]} d^2(j,i)~. \end{align}\tag{9}\]

Consider the time \(\hat{\alpha}_t-\varepsilon\) for each client \(t\), where \(\varepsilon>0\) is sufficiently small so that neither condition in the algorithm is satisfied between time \(\hat{\alpha}_t-\varepsilon\) (inclusive) and \(\hat{\alpha}_t\) (exclusive). For simplicity, we define \(DC^t\) as the set of directly connected clients in \([t-1]\) (i.e., \(DC\cap [t-1]\)) at time \(\hat{\alpha}_t-\varepsilon\), and \(IC^t\) as the set of indirectly connected or active clients in \([t-1]\) (i.e., \((A\cup IC)\cap [t-1]\)) at time \(\hat{\alpha}_t-\varepsilon\). Moreover, we use \(S^t\) to denote the set \(S\) at time \(\hat{\alpha}_t-\varepsilon\). Because the sets \(IC, DC, A\) are disjoint and the set \(DC\) only grows, we have the following observation.

Suppose that \(r_{jt}\) (for \(t\in [s]\) and \(j\in DC^t\)) is the distance of client \(j\in DC^t\) to its closest opened facility at time \(\hat{\alpha}_t-\varepsilon\) (i.e., \(d(j,S^t)\)). One simple observation is that \(r_{jt}\leq R_j\) for each \(j\in DC^t\), as the set \(S\) only grows in the algorithm.

Next, we will establish two sets of constraints for the variables \((\alpha^*_j)_{j\in [s]}\), \((r_{jt})_{1\leq j<t\leq s}\), and \((R_j)_{j\in [s]}\). The first set of constraints is as follows. The proof of the first set of constraints relies on the fact that for each client \(t\in [s]\), its \(\alpha\)-value at time \(\hat{\alpha}_t-\varepsilon\) is strictly less than its distance to any open facility.

Our second set of constraints is as follows. The proof of the second set of constraints relies on the fact that for each client \(t\in [s]\), the total bid to the fixed facility \(i\) at time \(\hat{\alpha}_t-\varepsilon\) is strictly less than \(\hat{f}\), as otherwise \(t\) will be connected to \(i\) at time \(\hat{\alpha}_t-\varepsilon\) and become inactive, contradicting the definition of \(\hat{\alpha}_t\).

Finally, we only need to prove the following lemmas to complete the proof of 4. Because we will use the weighted version of [lem:dual-feasibility-simpler-metric] [lem:dual-feasibility-simpler-euclidean] in the log-adaptive algorithm and walking between two solutions later, we will directly claim the weighted generalization. To establish 10, we only need to consider the special case of them where \(w(j)=1\) for each \(j\in [s]\).

5 Establishing the LMP Approximations↩︎

In this section, we will prove [lem:dual-feasibility-simpler-euclidean]. For convenience, we restate the lemma below. Using the proof of [lem:dual-feasibility-simpler-euclidean], we will also prove a generalization of it, [lem:dual-feasibility-simpler-euclidean-special], which will be used in the later sections for the log-adaptive and robust analysis of the greedy algorithm for Euclidean metrics. Throughout this section, we will assume that \(\rho=2\). For simplicity, in the rest of this section, we will define \(\Phi=\sum_{j\in [s]} w(j) \cdot d^2(j,i)\).

Readers can assume \(w(j)=1\) for all \(j\in [s]\) in the first reading, which is enough for understanding the main ideas and establishing the LMP approximations.

5.1 Proof Roadmap↩︎

To combine the constraints in the lemma, we will create a set of non-negative parameters \(p_{tj\ell}\) for each \(t\in [s], j\in DC^t, \ell\in IC^t \cup \{t\}\) such that \(\sum_{\ell\in IC^t \cup\{t\}} p_{tj\ell}=1\) for each \(t\in [s]\) and each \(j\in DC^t\). ⁴

Using these parameters, we can lower bound the LHS of the third constraint [constr:euclidean-third] by adding the following quantity which is no greater than 0, according to the second constraint [constr:euclidean-second]: \[\begin{align} 0\geq \sum_{j\in DC^t} w(j)& \cdot \sum_{\ell\in IC^t\cup\{t\}} p_{tj\ell}(\alpha^*_\ell - (r_{jt}+d(j,\ell))^2) \\ &= \sum_{j\in DC^t} w(j) \cdot \left(-r_{jt}^2 + \sum_{\ell\in IC^t\cup\{t\}} p_{tj\ell} \left(\alpha^*_\ell - 2r_{jt}\cdot d(j,\ell) - d^2(j,\ell)\right)\right) \end{align}\] Therefore, we can lower bound the LHS of the third constraint [constr:euclidean-third] by the following quantity, where we use \(r_{jt}\leq R_j\) in the second inequality: \[\begin{align} &\text{LHS of constraint \ref{constr:euclidean-third}} \\ &~\geq \sum_{\ell\in [t:s]} w(\ell)\cdot \alpha^*_t + \sum_{\ell \in IC^t} w(\ell) \cdot \alpha^*_\ell + \sum_{j\in DC^t} w(j) \cdot \sum_{\ell\in IC^t\cup\{t\}} p_{tj\ell} \left(\alpha^*_\ell - 2r_{jt}\cdot d(j,\ell) - d^2(j,\ell)\right)\\ &~\geq \sum_{\ell\in [t:s]} w(\ell)\cdot \alpha^*_t + \sum_{\ell \in IC^t} w(\ell)\cdot \alpha^*_\ell + \sum_{j\in DC^t} w(j) \cdot \sum_{\ell\in IC^t\cup\{t\}} p_{tj\ell} \left(\alpha^*_\ell - 2R_j\cdot d(j,\ell) - d^2(j,\ell)\right). \end{align}\] Putting this inequality back to the third constraint [constr:euclidean-third], we have \[\begin{align} & \sum_{\ell\in [t:s]} w(\ell) \alpha^*_t + \sum_{\ell \in IC^t} w(\ell) \alpha^*_\ell + \sum_{j\in DC^t} w(j) \sum_{\ell\in IC^t\cup\{t\}} p_{tj\ell} \alpha^*_\ell \\ \leq& \hat{f} + \rho\sum_{\ell \in IC^t \cup [t:s]} w(\ell) d^2(\ell,i) + \sum_{j\in DC^t} w(j) d^2(j,i) + \sum_{j\in DC^t} w(j) p_{tj\ell} \sum_{\ell\in IC^t\cup \{t\}} 2R_j d(j,\ell) + d^2(j,\ell) \\ =& \hat{f} + 2\Phi - \sum_{j\in DC^t} w(j) d^2(j,i) + \sum_{j\in DC^t} w(j) p_{tj\ell} \sum_{\ell\in IC^t\cup \{t\}} 2R_j d(j,\ell) + d^2(j,\ell) \end{align}\]

Further, we create another set of non-negative parameters \(\beta_t\) for each \(t\in [s]\).⁵ The goal is to ensure that \(\sum_{t\in [s]} \beta_t \leq 1\) and the following inequality holds: \[\begin{align} \sum_{t\in [s]} \beta_t \cdot \left(\sum_{\ell \in [t:s]} w(\ell) \alpha^*_t + \sum_{\ell \in IC^t} w(\ell) \alpha^*_\ell + \sum_{j\in DC^t} w(j) \sum_{\ell\in IC^t\cup\{t\}} p_{tj\ell} \alpha^*_\ell\right) = \sum_{j \in [s]} w(j) \alpha^*_j ~. \end{align}\] In this way, we can have the following upper bound: \[\begin{align} &\sum_{j\in [s]} w(j)\left(\alpha^*_j - (\rho-1)R_j^2\right) \leq \hat{f} + 2 \Phi \notag \\ \label{eqn:lmp-dual-feasibility-goal-2} &~~ \underbrace{- \sum_{j\in [s]} w(j) R_j^2 + \sum_{t\in [s]} \beta_t \cdot \left(\sum_{j\in DC^t} w(j)p_{tj\ell} \left(\sum_{\ell\in IC^t\cup \{t\}} 2R_j d(j,\ell) + d^2(j,\ell)\right) - w(j) d^2(j,i)\right) }_{\text{the extra cost } \mathsf{EC}}. \end{align}\tag{14}\]

To obtain [lem:dual-feasibility-simpler-euclidean], it suffices to explicitly define the parameters \((p_{tj\ell})_{t\in [s], j\in DC^t, \ell\in [s]}\), \((\beta_{t})_{t\in IC^t\cup\{t\}}\) and prove the extra cost terms after \(\hat{f} + 2\Phi\) in the RHS of 14 is no greater than \((\Gamma-\rho)\cdot \Phi\), where \(\Gamma\) depends on whether the underlying metric is metric or Euclidean according to [lem:dual-feasibility-simpler-euclidean].

5.2 Defining the parameters↩︎

In this subsection, we will define the parameters \((p_{tj\ell})_{t\in [s], j\in DC^t, \ell\in [s]}\) and \((\beta_{t})_{t\in [s]}\) explicitly. These parameters are essentially the same as those in [21], but we write them in an explicit formula and prove some new properties here⁶. For the further establishment of the LMP approximation guarantee for this paper, we will also define \(P_{j\ell}\), which indicates how much \(d^2(j,\ell)\) will be charged in the extra cost term \(\mathsf{EC}\), and \(d^+_j\), which indicates how much the extra cost is reduced by \(d^2(j,i)\) (see 14 for a formal statement). After defining these parameters, we will prove 12 in next subsection.

For simplicity, we will use \(\mathsf{suf}(j)=\sum_{\ell \geq j} w(\ell)\) to denote the suffix sum of \(w(j)\) for each \(j\in [s]\). For each \(t\in [s]\) and each \(u\in [s]\), we consider \[\begin{align} \pi_t(u) = \prod_{j\in DC^t: j<u} \frac{\mathsf{suf}(j)}{\mathsf{suf}(j+1)}~, \end{align}\] where if \([u-1]\cap DC^t=\emptyset\), we simply let \(\pi_t(u)=1\). We define the following coefficients for each \(t\in [s]\) and each \(j\in [s]\): \[\begin{align} \label{eqn:euclidean-eta-tl} \eta_{t\ell} = \begin{cases} \mathsf{suf}(t) \cdot \pi_t(t) & \text{if } \ell=t~,\\ \hfil w(\ell) \cdot \pi_t(\ell) \hfil & \text{if } \ell\in IC^t~,\\ \hfil 0 \hfil & \text{otherwise}~, \end{cases} \end{align}\tag{15}\]

First, we present the definitions of the parameters \(p_{tj\ell}\) for each \(t\in [s]\), \(j\in DC^t\), and \(\ell \in IC^t \cup \{t\}\), so that \(\sum_{\ell \in IC^t \cup \{t\}} p_{tj\ell} =1\) for each \(t\in [s]\) and each \(j\in DC^t\). For simplicity, we define \(p_{tj\ell}=0\) if \(\ell \notin IC^t\cup\{t\}\). The definition of \(p_{tj\ell}\) will follow the rule below, \[\begin{align} \label{eqn:euclidean-p-tjl} p_{tj\ell} \propto \begin{cases} 0 & \text{if } \ell \le j \text{ or } \ell \notin IC^t\cup\{t\}~,\\ \eta_{t\ell} & \text{if } \ell>j \text{ and } \ell \in IC^t \cup \{t\}~. \end{cases} \end{align}\tag{16}\] Alternatively, we can define for each \(t\in [s]\) and each \(j\in DC^t\), the suffix sum of \(\eta_{t\ell}\) as follows: \[\begin{align} \label{eqn:euclidean-sigma} \sigma_t(j) = \sum_{\ell\geq j} \eta_{t\ell}~, \end{align}\tag{17}\] For any \(t\in[s]\), \(j\in DC^t\), and those \(\ell \in [j+1:t] \cap(IC^t \cup \{t\})\), we can rewrite \(p_{tj\ell } = \eta_{t\ell}/\sigma_t(j+1)\) (or, \(\eta_{t\ell}/\sigma_t(j)\) because \(\eta_{tj}=0\)). For simplicity, we assume that \(p_{tj\ell}=0\) if \(j\notin DC^t\). It is clear that \(\sum_{\ell\in IC^t \cup \{t\}} p_{tj\ell} =1\) for each \(t\in [s]\) and each \(j\in DC^t\). It is straightforward that we have \(p_{tj\ell}\geq 0\). Further, we have the following properties:

Next, we present the definitions of the parameters \(\beta_t\) for each \(t\in [s]\) by the following recurrence relation (note that \(\eta_{jj}>0\) for any \(j\in [s]\)): \[\begin{align} \label{eqn:Euclidean-beta-recurrence} \beta_j \cdot \eta_{jj} = w(j) - \sum_{t>j} \beta_t \cdot \eta_{tj}. \end{align}\tag{19}\] The definitions are in the same spirit as the definitions of the parameters \(\beta_t\) in the metric case (19). Note that because we have \(\eta_{tj}=0\) for any \(t<j\), we have

We will defer the proof that \(\beta_j\geq 0\) to 17 and the proof that \(\sum_{t\in [s]} \beta_t \leq 1\) to 7. Using the above definition and 13, we can show that

Finally, for simplicity, we construct the matrix \(P\), in which, for any \(j<\ell\), \(P_{j\ell}\) equals the coefficient (divided by \(2w(j)\)) of \(d(j,\ell)\) in the extra cost term \(\mathsf{EC}\) in 14 . We define \(P_{j\ell}=0\) for any \(j\geq \ell\). Intuitively, we can think of \(P_{j\ell}\) as how much we will charge \(d^2(j,\ell)\) in the extra cost term \(\mathsf{EC}\) (see 14 for a formal statement). Because we have \(p_{tj\ell}\neq 0\) only when \(j\in DC^t\), \(\ell\in IC^t\cup\{t\}\) and \(j<\ell\), we have following simple expression for \(P_{j\ell}\): \[\begin{align} P_{j\ell} = \sum_{t\in [s]} \beta_t \cdot p_{tj\ell}~. \end{align}\] We also define \(d^+_j = \sum_{\ell>j} P_{j\ell}\). Intuitively, this indicates how much the extra cost is reduced by \(d^2(j,i)\) (see 14 for a formal statement).

Further, we present a simple characterization of \(d^+_j\) in terms of \(\{\beta_t\}\), which will be used in the next subsection to upper bound the extra cost term \(\mathsf{EC}\) in 14 by \(\Phi\). We define \(t_j\geq j\) as the threshold such that for any \(t>t_j\), we have \(j\in DC^t\), and for any \(j<t\le t_j\), we have \(j\notin DC^t\).

5.3 Proof of [lem:dual-feasibility-simpler-euclidean]↩︎

5.3.1 Step 1: Bounding the extra cost term via \(P_{j\ell}\) and \(d^+_j\)↩︎

Using above definitions, we can rewrite the extra cost term in 14 :

5.3.2 Step 2: Bounding the Laplacian matrix of \(P\)↩︎

The following lemma upper bounds the Laplacian matrix of \(P\) by the identity matrix, which is essentially a key lemma in [21] to establish the LMP 4-approximation for Euclidean \(k\)-means. For the sake of completeness, we prove this lemma in 12.1.

As a corollary, we have the following properties of the matrix \(P\), which give bounds for the first term of 14.

Another direct corollary of this lemma is that the sum of \(P\) on any column \(\ell\) is upper bounded by \(w(\ell) (1-d^+_\ell)\).

5.3.3 Step 3: A key property on the conditional column sum of \(P\)↩︎

The third key property is an upper bound for the conditional column sum of \(P\) in terms of \(d^+_j\) and \(w(\ell)\), which will be used to upper bound the second and the third term in 14. Intuitively, consider the \(\ell\)-th column (corresponding to client \(\ell\)) in the matrix \(P\). The conditional column sum of \(P\) is the sum of all rows \(j\) in this column whose corresponding client \(j\) has \(d^+_j>\theta\). For a better understanding, \(d^+_j\) can be interpreted as a parameter indicating whether client \(j\) becomes directly connected early in the process of the greedy algorithm (according to [obv:dplus-vs46-sum-beta]). If \(d^+_j>d^+_{j'}\) for other client \(j'\), it means that client \(j\) becomes directly connected earlier than \(j'\). Hence, the conditional column sum of \(P\) can be interpreted as the sum of the coefficients of a set of clients \(j<\ell\) that become directly connected at some time point before \(\ell\) becomes directly connected in the process of the greedy algorithm. The reason we only need to take care of those clients \(j\) is because we only charge \(d^2(j,\ell)\) for those clients \(j\) in the upper bound of the extra cost term \(\mathsf{EC}\) in 14 (i.e., \(P_{j'\ell}=0\) for other clients \(j'\)).

To prove this lemma, we start with the following simple observation, which is a direct consequence of the definition of \(\pi_t(j)\).

Our first intermediate lemma is an upper bound for the conditional column sum of \(P\) restricted to a subset of clients \(DC^u\) for any \(u\in [s]\).

Next, we show the first deferred proof that \(\beta_u\geq 0\) for any \(u\in [s]\).

We show a lemma that could imply the second deferred proof that \(\sum_{t\in [s]} \beta_t \leq 1\). This lemma is more general in the sense that it also gives an upper bound for the tail sum of \(\beta_t\), which is useful for establishing 15, the focus of this part.

By choosing \(u=1\) for 22 in the proof, we have the following corollary as we have \(\pi_t(u)=1\) for \(t\in [s]\).

This further implies the following corollary because of [obv:dplus-vs46-sum-beta].

Finally, we are ready to present the proof of the lemma in this part.

5.3.4 Step 4(a): Proof for the Euclidean case↩︎

In this part, we will prove [lem:dual-feasibility-simpler-euclidean] for the Euclidean case. We will consider \(\Gamma = 3+\ln 2\) here. Without loss of generality, we assume the facility \(i\) is located at \(x_i=0^{\mathsf{dim}}\) and the corresponding \(s\) clients are located at \(x_1,\dots, x_s\) respectively. According to 14 14, we have \[\begin{align} &\sum_{j\in [s]} w(j) \left(\alpha^*_j - (\rho-1)R_j^2\right) \\ &\quad \leq \hat{f} + 2\Phi + \sum_{j<\ell} w(j) P_{j\ell} d^2(j,\ell) + \sum_{j\in [s]} w(j) d^+_j \sum_{\ell>j} P_{j\ell} d^2(j,\ell) - \sum_{j\in [s]} w(j) d^+_j d^2(j,i) \\ &\quad = \hat{f} + 2\Phi + \sum_{j<\ell} w(j) P_{j\ell} \|x_j-x_{\ell}\|^2 + \sum_{j\in [s]} w(j) d^+_j \sum_{\ell>j} P_{j\ell} \|x_j-x_{\ell}\|^2 - \sum_{j\in [s]} w(j) d^+_j \|x_j\|^2 \end{align}\] Because of 4, the first term after \(\hat{f} + 2\Phi\) can be upper bounded by \(\Phi\). Next, we will focus on upper bounding the remaining term (shown below) by \(\ln(2)\cdot \Phi\). \[\sum_{j\in [s]} w(j) d^+_j \cdot \sum_{\ell>j} P_{j\ell} \cdot \|x_j-x_\ell\|^2 - \sum_{j\in [s]} w(j) d^+_j \cdot \|x_j\|^2~.\]

Using 15, for any \(\theta\in [0,1]\), we have \[\begin{align} \sum_{j<\ell} w(j)&d^+_j P_{j\ell} \|x_j-x_\ell\|^2 - \sum_{j\in [s]} w(j) d^+_j \|x_j\|^2 \\ &\leq \theta \sum_{j<\ell} w(j)P_{j\ell}\cdot \|x_j-x_\ell\|^2 + \sum_{j<\ell} [d^+_j-\theta]^+ \cdot w(j)P_{j\ell} \|x_j-x_\ell\|^2 - \sum_{j\in [s]} w(j) d^+_j \|x_j\|^2 \\ &\leq \theta \cdot \Phi + \underbrace{\sum_{j<\ell} [d^+_i-\theta]^+ \cdot w(j)P_{j\ell} \|x_j-x_\ell\|^2 - \sum_{j\in [s]} w(j)d^+_j \|x_j\|^2}_{S_2} } \end{align}\] For the term \(S_2\), we derive the following upper bound: \[\begin{align} S_2&=\sum_{j\in [s]: d^+_j>\theta} \sum_{\ell>j} (d^+_j-\theta) w(j)P_{j\ell} \cdot \|x_j-x_\ell\|^2 - \sum_{j\in [s]} w(j) d^+_j \|x_j\|^2 \\ &\leq \sum_{j\in [s]: d^+_j>\theta} \sum_{\ell>j} (d^+_j-\theta) w(j)P_{j\ell} \cdot \left(\frac{\|x_j\|^2}{d^+_j-\theta} + \frac{\|x_\ell\|^2}{1-d^+_j +\theta}\right) - \sum_{j\in [s]} w(j)d^+_j \|x_j\|^2 } \\ &\leq \sum_{j\in [s]: d^+_j>\theta} w(j)\left(\sum_{\ell>j} P_{j\ell}\right) \|x_j\|^2 + \sum_{j\in [s]: d^+_j>\theta} \sum_{\ell>j} \frac{d^+_j-\theta}{1-d^+_j+\theta} \cdot w(j)P_{j\ell} \|x_\ell\|^2 - \sum_{j\in [s]} w(j)d^+_j \|x_j\|^2 \\ &\leq \sum_{j\in [s]: d^+_j>\theta} \sum_{\ell>j} \frac{d^+_j-\theta}{1-d^+_j+\theta} \cdot w(j)P_{j\ell} \|x_\ell\|^2\\ &= \sum_{\ell\in [s]} \|x_\ell\|^2 \cdot \sum_{j<\ell: d^+_j>\theta} \frac{d^+_j-\theta}{1-d^+_j+\theta} \cdot w(j)P_{j\ell}~. \end{align}\] Let \(g(y) = \frac{y-\theta}{1-y+\theta}\) for any \(y\in [\theta,1]\). This is easy to verify that \(g(\theta)=0\). Therefore, we have \(g(y) = \int_{\theta}^y g'(z) dz\). Hence, we can further upper bound \(S_2\) by \[\begin{align} S_2 &\leq \sum_{\ell \in [s]} \|x_\ell\|^2 \cdot \sum_{j<\ell: d^+_j>\theta} \int_{\theta}^{d^+_j} g'(z)\cdot w(j)P_{j\ell} dz \\ &= \sum_{\ell \in [s]} \|x_\ell\|^2 \cdot \int_{z=\theta}^1 g'(z) \cdot \sum_{j<\ell: d^+_j>z} w(j)P_{j\ell} dz~. \\ &\leq \sum_{\ell \in [s]} \|x_\ell\|^2 \cdot \int_{z=\theta}^1 g'(z) \cdot w(\ell)(1-z)dz } \\ &= \sum_{\ell \in [s]} w(\ell)\|x_\ell\|^2 \cdot \int_{z=\theta}^1 g'(z) (1-z)dz \end{align}\] Since \(g(z)=0\) when \(z=\theta\) and \(1-z=0\) when \(z=1\), we have \[\begin{align} \int_{z=\theta}^1 g'(z) \cdot (1-z)dz = \int_{z=\theta}^1 g(z)dz = \int_{z=\theta}^1 \frac{z-\theta}{1-z+\theta} dz = \theta-1-\ln(\theta)~. \end{align}\] We have \(S_2\leq (\theta-1-\ln(\theta)) \cdot \Phi\). Therefore, we have \[\begin{align} \sum_{j<\ell} w(j)d^+_j P_{j\ell} \|x_j-x_\ell\|^2 - \sum_{j\in [s]} w(j) d^+_j \|x_j\|^2 &\leq (\theta + \theta -1 - \ln(\theta))\cdot \Phi \\ &= (2\theta -1 - \ln(\theta)) \cdot \Phi~. \end{align}\] If we choose \(\theta = 1/2\), we get our desired bound that \[\sum_{j<\ell} w(j)d^+_j P_{j\ell} \|x_j-x_\ell\|^2 - \sum_{j\in [s]} w(j) d^+_j \|x_j\|^2 \leq \ln(2) \cdot \Phi~.\]

5.3.5 Step 4(b): Proof for the general metric case↩︎

In this part, we will prove [lem:dual-feasibility-simpler-euclidean] for the general metric case. We will consider \(\Gamma = 5\) here. For simplicity, we use \(a_j=d^2(j,i)\) for each \(j\in [s]\) in this part. According to 14 14, we have \[\begin{align} &\sum_{j\in [s]} w(j) \left(\alpha^*_j - (\rho-1)R_j^2\right) \\ &\quad \leq \hat{f} + 2\Phi + \sum_{j<\ell} w(j) P_{j\ell} d^2(j,\ell) + \sum_{j\in [s]} w(j) d^+_j \sum_{\ell>j} P_{j\ell} d^2(j,\ell) - \sum_{j\in [s]} w(j) d^+_j d^2(j,i) \end{align}\] First, because of 4, we can upper bound the first term after \(\hat{f}+2\Phi\) by \(2\Phi\). Next, we upper bound the remaining term, i.e., \[\begin{align} \label{eqn:remaining-term-general-metric} \sum_{j\in [s]} w(j) d^+_j \sum_{\ell>j} P_{j\ell} d^2(j,\ell) - \sum_{j\in [s]} w(j) d^+_j d^2(j,i) \end{align}\tag{23}\] by \(\Phi\) for the general metric case.

Using the definition of \(a_j\), we have \[\begin{align} \ref{eqn:remaining-term-general-metric} &\leq \sum_{j\in [s]} w(j) d^+_j \sum_{\ell>j} P_{j\ell} \cdot 2(a_j + a_\ell) - \sum_{j\in [s]} w(j) d^+_j a_j } \\ &= 2\sum_{j<\ell} w(j) d^+_j P_{j\ell} \cdot a_\ell + 2\sum_{j<\ell} w(j) d^+_j P_{j\ell} \cdot a_j - \sum_{j\in [s]} w(j) d^+_j a_j \\ &= 2\sum_{j<\ell} w(j) d^+_j P_{j\ell} \cdot a_\ell + \sum_{j\in [s]} w(j) (2(d^+_j)^2 - d^+_j) a_j\\ &= 2\sum_{j<\ell} w(j) \left(\int_{\theta=0}^{1} \mathbb{1}(d^+_j>\theta) d\theta\right) \cdot P_{j\ell} a_\ell + \sum_{j\in [s]} w(j) (2(d^+_j)^2 - d^+_j) a_j \\ &= 2\sum_{\ell\in [s]} \left(\int_{\theta=0}^{1} \sum_{j<\ell}\mathbb{1}(d^+_j>\theta) w(j) P_{j\ell} d\theta\right) \cdot a_\ell + \sum_{j\in [s]} w(j) (2(d^+_j)^2 - d^+_j) a_j \end{align}\] Note that for any \(\ell\in [s]\), the coefficient of \(a_\ell\) in the first term can be upper bounded as follows: \[\begin{align} \int_{\theta=0}^{1} \sum_{j<\ell}\mathbb{1}(d^+_j>\theta) w(j) P_{j\ell} d\theta &\leq \int_{\theta=0}^{d^+_\ell} \sum_{j<\ell} w(j) P_{j\ell} d\theta + \int_{\theta=d^+_\ell}^1 \sum_{j<\ell} \mathbb{1}(d^+_j>\theta) w(j) P_{j\ell} d\theta \\ &\leq d^+_\ell \cdot w(\ell) (1-d^+_\ell) + \int_{\theta=d^+_\ell}^1 w(\ell)(1-\theta) d\theta\ref{lem:sum-threshold-wPij-upperbound}} \\ &= w(\ell) \left(d^+_\ell (1-d^+_\ell) + \frac{(1-d^+_\ell)^2}{2}\right) \end{align}\] Therefore, we can further upper bound 23 by (the final inequality is because \(d^+_j\in [0,1]\) for each \(j\in [s]\)): \[\begin{align} \ref{eqn:remaining-term-general-metric} &\leq \sum_{j\in [s]} w(j) \cdot \left(2d^+_j(1-d^+_j) + (1-d^+_j)^2+2(d^+_j)^2 - d^+_j\right) \cdot a_j \\ &= \sum_{j\in [s]} w(j) \cdot \left(1-d^+_j+(d^+_j)^2\right)\cdot a_j \leq \Phi~. \end{align}\]

5.4 A generalization of [lem:dual-feasibility-simpler-euclidean] for the Euclidean case↩︎

In this subsection, we will prove a generalization of [lem:dual-feasibility-simpler-euclidean], where we only have some approximate version of the constraints. We will show that the same guarantee can be achieved up to a small additive factor. Later, we will use this generalization to obtain the bicriteria approximation guarantee for the Euclidean case. Because the proof is highly technical with dealing with the small additive factor and mostly follows the spirit of the earlier derivations, we defer it to 12.2.

6 An Improved LMP \(4.9\)-Approximation for the Metric Case↩︎

In this section, we will prove [lem:dual-feasibility-simpler-metric], restated below for convenience.

6.1 Unweighted↩︎

In this subsection, we prove [lem:dual-feasibility-simpler-metric] in the special case \(w(\ell)=1\) for all \(\ell\in[s]\). Define \(a_t \stackrel{\text{def}}{=}d(t,i)\) for all \(t\in[s].\) The final proof of this lemma for the unweighted case was obtained by GPT 5.2/5.4 Pro/Thinking, after several rounds of exchange with the authors implementing various theoretical ideas and experiments, and revised by the authors.

As briefly discussed in the overview in 2.1, the proof has four steps. First, from the inequalities at time \(\alpha_t-\varepsilon\) given in condition [constr:metric-third] of [lem:dual-feasibility-simpler-metric], we derive a family of constraints \(\eqref{eq:Pk-unweighted}\) by using the triangle inequality to lower bound the terms involving \(r_{jt}\). Second, we choose coefficients \((\beta)_{t\in[s]}\) so that, in the weighted combination of these constraints, the coefficient of each \(\alpha_t^*\) is exactly \(1\). Third, after collecting terms, we obtain a quadratic form in the distances to the facility \(i\). Finally, we bound the top eigenvalue of the corresponding matrix by the Collatz–Wielandt Formula, which yields the claimed constant.

6.1.0.1 Step 1: From \(\eqref{eq:Ik-unweighted}\) to \(\eqref{eq:Pk-unweighted}\).

For each \(t\in[s]\), given \(w(t)=1\) and \(a_t\stackrel{\text{def}}{=}d(t,i)\) for all \(t\in[s]\), rewriting and rearranging the condition [constr:metric-third] gives \[\begin{align} } &\sum_{j\in[t:s]} \alpha^*_t + \sum_{j \in IC^t} \alpha^*_j + \sum_{j\in DC^t} w(j) r_{jt}^2 \leq \hat{f} + \sum_{j \in IC^t \cup [t:s]} \rho d^2(j,i) + \sum_{j\in DC^t} w(j) d^2(j,i).\\ \label{eq:Ik-unweighted}\implies & \sum_{j\in DC^t}(r_{jt}^2-a_j^2) + \sum_{j\in IC^t}(\alpha_j^*-\rho a_j^2) + \sum_{j\ge t}(\alpha_t^*-\rho a_j^2) \le \hat{f}. \end{align}\tag{25}\] Fix \(t\in[s]\). For each \(j\in DC^t\), conditions [constr:metric-first]–[constr:metric-second] and the triangle inequality imply \[\alpha_t^* \le (r_{jt}+d(j,t))^2 \le (r_{jt}+a_j+a_t)^2 \le r_{jt}^2 + (a_j+a_t)^2 + 2(a_j+a_t)R_j.\] Rearranging and using \(2a_ja_t\le a_j^2+a_t^2\), we obtain \[r_{jt}^2-a_j^2 \ge \alpha_t^* - 2a_t^2 - 3a_j^2 - 2(a_j+a_t)R_j.\] Substituting into 25 gets us \[\begin{align} &\sum_{j\in DC^t}\alpha_t^* + \sum_{j\in IC^t}\alpha_j^* + \sum_{j\ge t}\alpha_t^* - \hat{f} \notag\\ &\qquad\le \sum_{j\in DC^t}\bigl(2a_t^2+3a_j^2+2(a_t+a_j)R_j\bigr) + \sum_{j\in IC^t}\rho a_j^2 + \sum_{j\ge t}\rho a_j^2. \label{eq:Pk-unweighted} \end{align}\tag{26}\]

6.1.0.2 Step 2: Choosing coefficients.

We now choose nonnegative coefficients \(\beta_1,\ldots,\beta_s\) and multiply each inequality \(\eqref{eq:Pk-unweighted}\) with \(\beta_t\) and sum over all \(t\in[s]\) such that the coefficient of each \(\alpha_t^*\) is exactly \(1\). This is very similar to the proof of how coefficients are derived in [17].

6.1.0.3 Step 3: The quadratic form

We next compute the quadratic form produced by the weighted sum of \(\eqref{eq:Pk-unweighted}_{t\in[s]}\).

6.1.0.4 Step 4: Matrix Form

It remains to bound \(Q_\rho\).

We now pick \(\rho=5/2\) and replace \(\eta_t\) by \(g_t\), since \(\eta_t\leq g_t\) (19). For every \(t\in[s]\) and every \(x\in\mathbb{R}_{\ge 0}^s\), 21 gives \[\label{eq:finalkentry} (Kx)_t \le \Bigl(2A_t + \frac{1}{2} g_t + \frac{2}{3} g_t^2\Bigr)x_t + \frac{2}{3} g_t U_t(x) + \frac{2}{3} \beta_t\sum_{u\in DC^t} U_u(x) + \frac{2}{3} \beta_t\sum_{u\in DC^t} g_u x_u.\tag{32}\]

The next lemma applies the Collatz–Wielandt Formula [29]

The remaining task is to prove an upper bound on \(F(x)\) as follows:

We can now finish the proof of the unweighted version of [lem:dual-feasibility-simpler-metric]. Using 20 with \(\rho=5/2\), we obtain \[\begin{align} \sum_{t\in[s]} \alpha_t^*-(\rho-1)R_t^2 - \hat{f} &\le \frac{5}{2} \sum_{t\in[s]} a_t^2 + Q_{5/2}(a) } \\ &\le \frac{5}{2} \sum_{t\in[s]} a_t^2 + \lambda_{\max}(K)\sum_{t\in[s]} a_t^2 } \\ &\le \frac{5}{2} \sum_{t\in[s]} a_t^2 + \left(\max_{x\in[1,\tau]}F(x)\right)\sum_{t\in[s]} a_t^2 } \\ &\le \frac{5}{2} \sum_{t\in[s]} a_t^2 + \frac{12}{5}\sum_{t\in[s]} a_t^2 } \\ &= 4.9\sum_{t\in[s]} a_t^2. \end{align}\] This proves [lem:dual-feasibility-simpler-metric] in the unweighted case, i.e., \(w(\ell)=1\) for all \(\ell\in[s]\).

6.2 Weighted↩︎

We now prove [lem:dual-feasibility-simpler-metric] in its full generality for any weight function \(w:[s]\to\mathbb{R}_{> 0}\).

To extend the argument to the weighted setting, the overall structure is the same as in the unweighted case: given the condition (3) of [lem:dual-feasibility-simpler-metric], we derive inequalities \(\eqref{eq:Pt-weighted-short}_{t\in[s]}\), choose coefficients so that the \(\alpha_t^*\) terms sum up to \(w(t)\) for each \(t\in[s]\), rewrite the resulting bound as a quadratic form, and then upper bound its top eigenvalue. The new issue is the bookkeeping. The cardinalities of clients are replaced by weighted masses, the normalization of the coefficients changes accordingly. After introducing the appropriate weighted analogues of the unweighted parameters and rescaling by \(\sqrt{w(t)}\), however, the same spectral argument goes through.

Since the proof is a weighted analogue of the unweighted subsection, we mention only the differences in the proof of the coefficients and in the quadratic form, and defer the full proof to 13.2.

Define

\[a_t \stackrel{\text{def}}{=}d(t,i), \qquad W(S)\stackrel{\text{def}}{=}\sum_{j\in S}w(j), \qquad W_t\stackrel{\text{def}}{=}W([t:s]), \qquad M_t\stackrel{\text{def}}{=}W(DC^t).\]

For each \(t\in[s]\), condition (3) gives \[\label{eq:Ik-weighted-short}\sum_{j\in DC^t} w(j)(r_{jt}^2-a_j^2) + \sum_{j\in IC^t} w(j)(\alpha_j^*-\rho a_j^2) + \sum_{j\ge t} w(j)(\alpha_t^*-\rho a_j^2) \le \hat{f}.\tag{41}\] Exactly as in the unweighted case, for every \(j\in DC^t\), \[r_{jt}^2-a_j^2 \ge \alpha_t^* - 2a_t^2 - 3a_j^2 - 2(a_j+a_t)R_j.\] Substituting into 41 yields \[\begin{align} &M_t\alpha_t^* + \sum_{j\in IC^t} w(j)\alpha_j^* + W_t\alpha_t^* - \hat{f} \notag\\ &\qquad\le 2M_ta_t^2 + \sum_{j\in DC^t} w(j)\bigl(3a_j^2+2(a_t+a_j)R_j\bigr) + \sum_{j\in IC^t}\rho w(j)a_j^2 + \sum_{j\ge t}\rho w(j)a_j^2. \label{eq:Pt-weighted-short} \end{align}\tag{42}\]

The proof is the same as that of 19, with cardinalities replaced by weighted masses. The only changes are the following two weighted analogues of Inequalities 29 , 30 \[(W_{t+1}+M_t)\sum_{u>t}\beta_u\le W_{t+1}, \qquad (W_t+M_t)g_t\le W_t.\]

We can now conclude exactly as in the unweighted case: \[\begin{align} &\sum_{t\in[s]} w(t)\bigl(\alpha_t^*-(\rho-1)R_t^2\bigr)-\hat{f} \\ &\qquad\le \frac{5}{2}\sum_{t\in[s]} w(t)a_t^2 + Q_{5/2}^w(a) } \\ &\qquad\le \frac{5}{2}\sum_{t\in[s]} w(t)a_t^2 + \lambda_{\max}(B)\sum_{t\in[s]} w(t)a_t^2 } \\ &\qquad\le 4.9\sum_{t\in[s]} w(t)a_t^2. } \end{align}\] This proves [lem:dual-feasibility-simpler-metric].

7 Log-Adaptive Algorithm for Weighted Facility Location↩︎

Let \(\varepsilon\in (0, 10^{-100})\) be a sufficiently small constant. We choose \(\delta = \varepsilon/2\). In weighted facility location, each client \(j \in D\) has a weight \(w(j) > 0\) and the objective is to minimize \(\sum_{j \in D} w(j) \cdot d^2(j, S) + |S|f\).

We present the modified greedy algorithm for weighted facility location that works in \(O(\log n / \varepsilon^3)\) phases.

The openability can be checked by solving a linear program.

4 is our log-adaptive algorithm. It works in \(O(\log n/\varepsilon^3)\) phases where, in the \(i\)th phase, every active client \(j\) has \(\alpha_j = (1+\varepsilon^3)^{i} = \theta\) 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\) becomes immediately put into \(DC\) after \(i\) is open. Unlike 2, we don’t lower the \(\alpha\)-values for clients in \(IC\) immediately after a facility is opened. Instead, we only lower the \(\alpha\)-values in the second stage of each phase.

For each client \(j \in D\), we explicitly track its first direct-connection cost, \(R_j\). We set \(R_j = \sqrt{1-\delta}d(j, i^*)\) where \(i^*\) is the facility that caused \(j\) to be moved to \(DC\) for the first time. For a client that remains only indirectly connected or active, we set \(R_j = 0\).

In the rest of this subsection, let \((\alpha^*_j)_{j\in D}\) be the final dual solution produced by the algorithm.

7.1 Dual Payment↩︎

7.2 Dual Feasibility↩︎

Next, we fix a facility \(i\) and show that the dual-feasibility for \(i\) is satisfied at the end of the algorithm. We say that facility \(i\) is frozen if there exists some open facility \(i_0\) such that \[d(i,i_0) \le 2\sqrt{\varepsilon\theta}~.\]

7.2.0.1 Modified \(\alpha\)-values.

For analysis purposes, we will introduce a new modification of the \(\alpha\)-values for each client \(j\) in terms of the fixed facility \(i\). We will denote this modified \(\alpha\)-value as \(\alpha^c_j\). This modified \(\alpha^c\)-value will be used to write the analog of the constraints in [lem:p-tjl,lem:beta-t].

The modifications happen only when we open a facility \(i_0\) that does not freeze \(i\), i.e., when \(d(i,i_0) > 2\sqrt{\varepsilon\theta}\). More specifically, at the stage with time \(\theta\), if we open such a facility \(i_0\) that does not freeze \(i\), we keep \(\alpha^c_j\) as the value before opening \(i_0\). That is, we don’t increase \(\alpha^c_j\) to \(\tau_j\) if \(j\in B(i_0,\sqrt{\varepsilon\theta})\) (step 1b) for such \(i_0\). In other steps (step 1b for a facility \(i_0\) that freezes \(i\) or step 2), we still modify \(\alpha^c_j\) in the same way as we modify \(\alpha_j\).

In the rest of this subsection, we will consider \(\alpha^c\) instead of \(\alpha^*\) in writing the constraints. To ensure that the \(\alpha^*\)-value can ensure dual feasibility, we will show that \(\alpha_j^c\) and \(\alpha_j^*\) are close to each other for every client \(j\).

Now, we assume that \[\begin{align} \label{eqn:Dstar} D^*= \{j\in D: \alpha^c_j - (\rho-1)R^2_j > \Gamma d^2(j,i)\}~. \end{align}\tag{45}\] Then, it suffices to show that \[\begin{align} \sum_{j \in D^*} w(j) \cdot (\alpha_j^c - (\rho-1)R^2_j) \le \hat{f} + (\Gamma+200\sqrt{\varepsilon}) \cdot \sum_{j \in D^*} w(j) \cdot d^2(j,i)~. \end{align}\] For simplicity, let us number the clients in \(D^*\) from \(1\) to \(s=|D^*|\) in the order they become inactive, breaking ties by order of increasing \(\alpha^c\)- value. We will use \(D^*\) and \([s]\) interchangeably.

7.2.0.2 Establishing the constraints.

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

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

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

For any client \(t\) and \(j\in DC^t\), define \[r_{jt}:= \sqrt{1-\delta}d(j,S^t).\] Next, we will prove 30 for the general metric setting and the Euclidean setting separately based appropriate choices of \(\rho\) and \(\Gamma\).

General metric↩︎

Now, we will show 30 for \(\rho = 2.5\) and \(\Gamma = 4.9\). To this end, let’s show that \((IC^t), (DC^t), (\alpha_j^c), (R_j), (r_{jt})\) satisfy all the constraints of [lem:dual-feasibility-simpler-metric], restated below for convenience.

By definition of \((IC^t)_t\) and \((DC^t)_t\), the two constraints on them hold. Let \(j\in DC^t\). Since \(R_j\) is the distance to the first facility that \(j\) directly connects to, and \(d(j, S^t)\) is non-increasing in \(t\), we have \(r_{jt}\le R_j\) for any \(t>j\). Now, let \(t\in [s]\), we have: \[\begin{align} \alpha_t^c &\le \alpha_t^*\\ &\le (1-\delta)d^2(t, S^t)\\ &\le (\sqrt{1-\delta}d(t, j) + \sqrt{1-\delta}d(j, S^t))^2\\ &\le (r_{jt} + d(j, t))^2. \end{align}\] [clm:3rdconstraint] gives the last constraint.

Hence, by [lem:dual-feasibility-simpler-metric], we have: \[\sum_{j\in [s]}w(j)\cdot (\alpha_j^c - (\rho-1)R_j^2)\le \hat{f} + \Gamma\cdot \sum_{j\in [s]}w(j)\cdot d^2(j,i).\] This completes the proof of 30 for \(\rho=2.5\) and \(\Gamma = 4.9\).

Euclidean metric↩︎

Now, we will show 30 for \(\rho=2\), \(\Gamma=3+\ln 2\), and \(d\) a Euclidean metric. Let \(a\) be the index of the last client that becomes inactive strictly before \(i\) is frozen.

First, let’s prove the following claim:

We will make several cases according to how \(i\) gets frozen.

7.2.0.7 Case 1: All clients in \(D^*\) become inactive strictly before \(i\) is frozen.

In this case, let’s show that \((IC^t), (DC^t), (\alpha_j^c), (R_j), (r_{jt})\) satisfy all the constraints of [lem:dual-feasibility-simpler-euclidean], restated below for convenience.

Since in this case \(a=s\), [clm:2nd-constraint-euclidean] implies constraint (2) is satisfied. [clm:3rdconstraint] implies constraint (3) is satisfied. Hence, we can apply [lem:dual-feasibility-simpler-euclidean], and get:

\[\sum_{j\in [s]}w(j)\cdot(\alpha_j^c-(\rho-1)R_j^2)\le \hat{f}+\Gamma\cdot \sum_{s\in [s]}w(j)\cdot d^2(j,i),\] which completes the proof in this case.

7.2.0.8 Case 2: Some clients in \(D^*\) become inactive when \(i\) is frozen.

In this case, we assume that \(i_0\) freezes \(i\). Let \(IC^F\) be the set of indirectly connected clients right before \(i_0\) is opened, and \[L:= \{j\in IC^F:d(j,i)>\sqrt{\varepsilon\theta}\}.\] Note that \(IC^F\) is different with \(IC^s\) in the way that \(IC^F\) does not include active clients right before \(i_0\) is opened, i.e., \(j\notin IC^F\) for any \(j\in [a+1:s]\) while \(IC^s=IC^F \cup [a+1:s-1]\).

We will further discuss two subcases to prove 30.

7.2.0.9 Subcase 2(a): \([a+1:s]\cap B(i_0,\sqrt{\varepsilon\theta})=\emptyset\) or \(\sum_{j\in L}w(j)>\varepsilon\cdot\sum_{j\in [a+1:s]}w(j)\).

We redefine \(DC^s\) and \(IC^s\) as the respectively directly connected and active or indirectly connected client at the beginning of the stage \(s\) becomes inactive. First, let’s bound the total increment of \(\alpha_j^c\) when \(i_0\) is opened.

Apply [clm:no95overbidding] at the beginning of the stage \(s\) becomes inactive, this lemma guarantees that \[\begin{align} w(s) (\alpha^c_s-\rho d^2(s,i)) + \sum_{\ell\in IC^s} w(\ell)(\alpha^c_\ell -\rho d^2(\ell,i)) + \sum_{j\in DC^s} w(j) (r^2_{js} - d^2(j,i)) \\ \leq \hat{f} + \varepsilon\sum_{j\in D^*} w(j)d^2(j,i). \end{align}\] Note that because we redefine \(IC^s, DC^s, S^s\) according to the sets \(IC,DC,S\) at the beginning of the stage \(s\) becomes inactive, it’s easy to see that [clm:2nd-constraint-euclidean] generalizes for \(t=s\). Thus, we can apply [lem:dual-feasibility-simpler-euclidean-special] (restated, see below), replacing \(\hat{f}\) with \(\hat{f} + \varepsilon\sum_{j\in D^*}w(j)d^2(j,i)\) and using \(L=\emptyset\) and \(\varepsilon'=0\) (the \(\varepsilon\)-parameter in the lemma). We get: \[\sum_{j\in [s]}w(j)\cdot (\alpha_j^c - (\rho-1)R_j^2)\le \hat{f} + (\Gamma + \varepsilon)\cdot\sum_{j\in [s]}w(j)\cdot d^2(j,i).\] This completes the proof in this case.

7.2.0.10 Subcase 2(b): \([a+1:s]\cap B(i_0,\sqrt{\varepsilon\theta})\neq \emptyset\) and \(\sum_{j\in L}w(j)\le \varepsilon\cdot\sum_{j\in [a+1:s]}w(j)\).

In this case, we define \(DC^s\) and \(IC^s\) as respectively the directly connected and the active or indirectly connected client right before \(s\) becomes inactive. Let’s show that \(a\), \(L\), \((IC^t), (DC^t), (\alpha_j^c), (R_j), (r_{jt})\) satisfy all the constraints of [lem:dual-feasibility-simpler-euclidean-special], restated below for convenience.

Just like for the other cases, constraints (a), (b), (c), (1), (2) and (3) hold. Remains to prove that constraint [constr:euclidean-special-second-approximate] is satisfied. Because \(L\subseteq IC^F\), \(IC^F\subseteq [a]\) and \(IC^F\subseteq IC^s\), we have \(L\subseteq [a]\cap IC^s\). Then, by proving the following claim, we establish constraint [constr:euclidean-special-second-approximate]. Note that because \(IC^F \setminus L \cup [a+1:s] = IC^s \setminus L\), the second and third bullet of the following claim together implies the second bullet of constraint [constr:euclidean-special-second-approximate]. The fourth bullet of the following claim implies the third bullet of constraint [constr:euclidean-special-second-approximate].

Thus, we can apply [lem:dual-feasibility-simpler-euclidean-special]. We get: \[\sum_{j\in [s]}w(j)\cdot(\alpha_j^c-(\rho-1)R_j^2)\le \hat{f} + (\Gamma + 200\sqrt{\varepsilon})\cdot\sum_{j\in [s]}w(j)\cdot d^2(j,i).\] This completes the proof.

8 Approximation with \(O(\log n / \varepsilon^3)\) Extra Centers↩︎

In this section, we present bi-criteria approximations with \(O(\log n / \varepsilon^3)\) extra centers, as claimed in 1 (restated in [thm:robust]). For simplicity, we will assume that \(\rho\in [2,3]\) in the rest of this section.

8.1 Pre-processing↩︎

For this section, we will consider the case where the distances between clients and facilities are lower bounded by \(1\), thanks to the following variant of 6:

We shall remark here that under this assumption, our log-adaptive 4 still runs in \(O(\log n/\varepsilon^3)\) phases.

8.2 Walking Between Two Solutions: Setup↩︎

This subsection is an adaptation of Subsection 5.2 of [21] to our new bidding strategy.

Given the description of the log-adaptive algorithm, we present our final bicriteria approximation MergeSolutions. It results in a solution that opens \(k+O(\log n)\) centers by maintaining two partial solutions that lead to opening at least \(k\) and at most \(k\) centers respectively, and gradually merging them into a single solution that opens exactly \(k+O(\log n)\) facilities, exploiting the \(O(\log n)\) adaptivity of the algorithm.

Before detailing the merging procedure 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\), some distance related to a free facility that will be introduced later, or some probability that a client is in \(IC\) or \(DC\).) 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, or that the dual constraint of some facilities is slightly violated. To address this, we generalize the definition of “openable” to “\(\eta\)-openable”, where the only difference is an \(O(n\eta)\) slack in the third and fourth conditions.

Note that, in the algorithm, each inactive client \(j\) is always either directly connected to some facility \(i^*\) or indirectly connected to some facility \(i^*\), where \(i^*\) is given by \(i^* = \mathop{\mathrm{arg\,min}}_{i\in S} d^2(j,i)\). In the rest of this section, we will maintain this invariant, but allow some flexibility to its connection type and \(\alpha\)-value. Specifically, we allow each client \(j\in D\) to have a set of copies \(D_j\), with each \(\tilde{j}\in D_j\) defined by (i) the connection type \(ct(\tilde{j})\), i.e., whether it is in \(A\), \(IC\) or \(DC\); (ii) the \(\alpha\)-value \(\tilde{\alpha}(\tilde{j})\); (iii) the probability mass \(\mu(\tilde{j})\). The probability mass should satisfy \(\sum_{\tilde{j}\in D_j} \mu(\tilde{j}) = 1\). In particular, for each client \(j\), \(D_j\) has two possibilities:

• \(D_j\) has only one copy that is of type \(A\) and with probability mass \(1\);

• \(D_j\) has some copies of type \(DC\) and at most one copy of type \(IC\).

In short, each \(D_j\) contains at most one copy that is of type \(IC\) or \(A\). We will say \(j\) is active if the first bullet holds, and inactive if the second bullet holds. In this way, we define the \(\alpha^I\)-value of a client as the \(\alpha\)-value of its active copy or the indirectly connected copy. If the client has no such copy, we say its \(\alpha^I\)-value is \(0\) (as this value will not be used). We remark this \(\alpha^I\)-value is different from the \(\alpha\)-values of the copies, and will only be used when we merge two solutions (8.3). Moreover, we will treat each copy \(\tilde{j}\) as a normal weighted client located at the same point as \(j\) in the metric space. When we compute the total bid, we average out over their different copies.

Further, we consider \(\mathcal{D}\) as the collection of (succinct) copies by dropping all the \(\alpha\)-values of directly connected copies from \(\cup_{j\in D} D_j\). This collection will be used to define the algorithm that merges two solutions. In this way, we can still deduce \(\alpha^I\) values for each client and this set can be independent of the \(\alpha\)-values of directly connected copies, which are not necessary in defining the algorithm.

Equivalently, for each inactive client \(j\), we consider its connection as a discrete probability distribution over (i) the connection type \(ct_j\), i.e., whether it is in \(IC\) or \(DC\); (ii) the \(\alpha\)-value of the client \(\alpha_j\). In particular, in the support of this distribution, there can be only one possible \(\alpha\)-value for indirect connects, i.e., there exists at most one value \(\alpha^I\) such that \(\Pr[(ct_j, \alpha_j) = (IC, \alpha^I)] > 0\). When we compute the total bid, we average out their contributions. Adapting to our new dual fitting approach, the expected contribution of client \(j\) to facility \(i\) now becomes: \[\begin{align} \mathbb{E}\left[(1-\delta)[d^2(j,S) - d^2(j, i)]^+ \mid ct_j=DC\right] \cdot \Pr[ct_j=DC]\\ + \mathbb{E}\left[[\alpha_j - \rho(1-\delta)d^2(j, i)]^+ \mid ct_j=IC\right] \cdot \Pr[ct_j=IC], \end{align}\] where the probability and the expectations are all taken over the discrete distribution. Letting \(\pi(j)=\Pr[ct_j=DC]\) and letting \(\alpha^I_j\) be the only possible \(\alpha\)-value for indirect connects of \(j\), we can rewrite the above as: \[\begin{align} \pi(j)\cdot (1-\delta)[d^2(j,S) - d^2(j, i)]^+ + (1-\pi(j))\cdot [\alpha^I_{j} - \rho(1-\delta)d^2(j, i)]^+ \end{align}\] For each active client \(j\), the connection distribution is a distribution with one element in the support having type \(ct_j=A\). Therefore, we can define \(\pi(j)=0\) and let \(\alpha^I_j\) be the \(\alpha\)-value of \(j\). In the rest of this section, we will use both points of view interchangeably.

Given the time \(\theta = (1+\varepsilon^3)^{p-1}\), the sequence \(\mathcal{H}=(H_1, \dots , H_p)\) of the sequences of facilities opened by each stage 1 of the algorithm that was executed so far, we easily can obtain the set of open facilities \(S\). Further given the collection of copies \(\mathcal{D}\) for the clients at the beginning of the phase, we can deduce the \(\alpha\)-value for each indirectly connected or active copy. Then, we can use the following process RetrieveNewCopies, which is the third step of stage 1 of our modified log-adaptive algorithm, to obtain the current collection of copies:

8.2.0.1 RetrieveNewCopies(\(\mathcal{D}, H_p\))

For each copy \(j\in \mathcal{D}\) that is active or indirectly connected:

• Let \(\alpha^I\) be its \(\alpha\)-value and \(j\) be the client it corresponds to.

• If \(\alpha^I \geq \rho(1-\delta)d^2(j, H_p)\), we modify the copy to be directly connected and drop its \(\alpha\)-value.

• Else if \(\alpha^I \geq (1-\delta)d^2(j, H_p)\), we modify the copy to be indirectly connected and keep its \(\alpha\)-value unchanged.

Given the output of the procedure, we can deduce all current copies for each client, and their corresponding connection types, probability mass and \(\alpha\)-values. Based on these variables, we have a new definition of \(\eta\)-openability tailored for our new dual fitting approach.

Similar as [13], [17], [21], we introduce the concept of free facilities. When we merge the two solutions, we introduce free facilities, whose opening costs are not necessarily paid. We will introduce at most \(O(1)\) such free facilities in each phase. However, we do not restrict the number of free copies a regular facility may have. Let \(S_{\mathsf{free}}\subseteq \mathcal{F}\) denote the multiset of free facilities that have been opened.

Moreover, we augment the metric space for the free facilities. For each free copy \(\tilde{i} \in S_{\mathsf{free}}\) of a regular facility \(i \in \mathcal{F}\), we define a parameter \(u(\tilde{i}) \ge 0\), which acts as the squared distance overhead of the free copy. Therefore, for any point \(x\) in the metric space (either a facility or a client), the distance squared is given by \(d^2(\tilde{i}, j) = u(\tilde{i}) + d^2(i, j)\).

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

Note that if \(d\) is Euclidean, we can have one additional dimension, which equals \(\sqrt{u(\tilde{i})}\) for any free facility \(\tilde{i}\) and equals \(0\) for any regular facility or client, to maintain the Euclidean property.

Finally, we will use the letters \(i\) and \(h\) to denote facilities. Specifically, we use \(i \in \mathcal{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 of the log-adaptive algorithm. 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 super-sequence 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.

We consider the following procedure CompleteCopies to return the collection of copies after opening a sequence of facilities \(H_p\) during the first stage of the phase. This procedure is identical to executing stage 2 of the log-adaptive algorithm on the collection of copies.

8.2.0.2 CompleteCopies\((\theta, \mathcal{D}, \mathcal{H}=(H_1,\dots, H_p))\)

Set \(S=\cup_{r\in [p]} H_r\) and \(\mathcal{D}'=\text{{\sc RetrieveNewCopies}}(\mathcal{D}, H_p)\). For each copy \(j\in \mathcal{D}'\):

• If \(j\) is of type \(IC\), set \(\alpha_j\leftarrow (1-\delta)d^2(j,S)\).

• If \(j\) is of type \(A\), set \(\alpha_j\leftarrow \min\{(1+\varepsilon^3)\theta, (1-\delta)d^2(j, S)\}\), and change its connection type to \(IC\) when \(\alpha_j=(1-\delta)d^2(j,S)\).

We remark that

• The log-adaptive algorithm generates \(0\)-valid sequences in the first stage of each phase. Specifically, it generates sequences \(\mathcal{H}= (H_1, H_2, \ldots, H_L)\) and a history of the collections of copies \(\mathcal{D}\mathcal{H}= (\mathcal{D}_0, \mathcal{D}_1, \ldots, \mathcal{D}_L)\), where each \(H_p\) is \(0\)-valid (with respect to \(\theta = (1+\varepsilon^3)^{p-1}\), \(\mathcal{D}_{p-1}\) and \(\mathcal{H}\)) and contains only regular facilities.

• The log-adaptive algorithm only maintains a single copy for each inactive client, i.e., each inactive client \(j\) is either directly connected to some facility \(i\) with probability \(1\) or indirectly connected to some facility \(i\) with probability \(1\). Therefore, each collection of copies \(\mathcal{D}_p\) maintained by the log-adaptive algorithm is perfectly valid (with respect to \(\theta=(1+\varepsilon^3)^{p-1}\), \(\mathcal{H}\) and \(\mathcal{D}_{p-1}\)).

• The set \(S\) of opened facilities is such that \(A = \emptyset\) at the final time \(\theta\).

For a general sequence of sequences \(\mathcal{H}= (H_1, H_2, \ldots, H_L)\) and a history of collection of copies \(\mathcal{D}\mathcal{H}= (\mathcal{D}_0, \mathcal{D}_1, \ldots, \mathcal{D}_L)\), let us say they consist of \(\eta\)-valid sequences and valid collections of copies if, for every \(p\in [L]\),

• \(H_p\) is an \(\eta\)-valid sequence with respect to \(\theta = (1+\varepsilon^3)^{p - 1}\), \(\mathcal{H}\), \(\mathcal{D}_{p-1}\); and

• \(\mathcal{D}_p\) is a valid collection of copies with respect to \(\theta = (1+\varepsilon^3)^{p-1}\), \(\mathcal{H}\), \(\mathcal{D}_{p-1}\).

It is easy to see from the definitions that the following claim holds.

Furthermore, we say that \(\mathcal{H}, \mathcal{D}\mathcal{H}\) form a solution with open facilities \(S =\cup_{r\in [L]} H_r\) if there are no active copies in \(\mathcal{D}_L\). Otherwise, we say that \(\mathcal{H}, \mathcal{D}\mathcal{H}\) form a partial solution.

The results in 7, namely 29 30, can thus be stated as: if \(\mathcal{H}, \mathcal{D}\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} w(j) \cdot d^2(j, S) \leq \frac{\Gamma}{1-\delta} \cdot \left( \text{opt}_{\text{LP}}(f) - f \cdot |S| \right)\,.\] We now allow an additive \(8n\eta\) slack in the third bullet (facility paying) and fourth bullet (dual feasibility) of the openability (Definition 3). In addition, in the second stage of each phase, we allow a directly connected copy \(j\) that satisfies \(\rho(1-\delta)d^2(j,S) \leq \alpha_j < \rho(1-\delta)d^2(j,S) + (1+\rho)\eta\) to split into two copies with arbitrary probability mass: one is identical, while the other is indirectly connected with \(\alpha\)-value \(\alpha_j=(1-\delta)d^2(j,S)\). All these modifications lead 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^3)\). To our help, we will use two basic routines defined as follows:

• \(\mathrm{\small CompleteSolution}_{f,u}(\mathcal{H}, \mathcal{D}\mathcal{H})\) takes a partial solution \(\mathcal{H}, \mathcal{D}\mathcal{H}\) as input and returns a solution by running the while-loop of the log-adaptive algorithm using \(f,u\) and the existing copies induced by the last collection in \(\mathcal{D}\mathcal{H}\).

  In particular, we also consider a slightly different input for this routine. If \(\mathcal{H}= (H_1, \ldots, H_{p})\) and \(\mathcal{D}\mathcal{H}= (\mathcal{D}_0, \ldots, \mathcal{D}_{p-1})\), where the first \(p-1\) sequences of \(\mathcal{H}\) form a partial solution with \(\mathcal{D}\mathcal{H}\), and \(H_p\) is \(\eta\)-valid. Then, we will first generate \(\mathcal{D}_p\) by CompleteCopies\((\theta=(1+\varepsilon^3)^{p-1}, \mathcal{D}_{p-1}, (H_1,\dots, H_p))\), append it to the end of \(\mathcal{D}\mathcal{H}\) (so that the resulting \(\mathcal{H},\mathcal{D}\mathcal{H}\) form a partial solution), and then start from phase \(p+1\) of the log-adaptive algorithm.

• \(\mathrm{\small CompleteSequence}_{f,u}(\mathcal{H}= (H_1, \ldots, H_{p-1}), H_p = \langle h_1, \ldots, h_\ell\rangle, \mathcal{D}\mathcal{H}=(\mathcal{D}_0,\dots,\mathcal{D}_{p-1}))\) where \(\mathcal{H}, \mathcal{D}\mathcal{H}\) form a partial solution, and \(H_p\) is an \(\eta\)-valid sequence except that it may not be maximal. \(\mathrm{\small CompleteSequence}\) then runs the while-loop of the log-adaptive algorithm starting with \(H_p\) and the existing copies induced by the last collection \(\mathcal{D}\mathcal{H}\) to return a maximal sequence of \(H_p\) (with respect to \(f, u\)).

We note that all facilities added by CompleteSolution and CompleteSequence are \(0\)-openable and regular. We also note that, given \(\mathcal{D}_{p-1}\) and \(\mathcal{H}\), CompleteSolution generates \(\mathcal{D}_{p}\) by CompleteCopies\((\theta=(1+\varepsilon^3)^{p-1}, \mathcal{D}_{p-1}, (H_1,\dots, H_p))\). Additionally, for intuition, observe that, with the above notation, the set of facilities opened by \(\mathrm{\small CompleteSolution}_{f,u}(\emptyset, (\mathcal{D}_0))\) is equivalent to that opened by the log-adaptive algorithm, where \(\mathcal{D}_0\) denotes the collection of copies in which each client has only one copy with \((ct=A,\mu=1,\alpha=1)\). (There is no free facility here.) To simplify notation, we only use \(\mathrm{\small CompleteSequence}\) when the partial solution \(\mathcal{H}\) and \(\mathcal{D}\mathcal{H}\) are 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.

8.3 MergeSolutions↩︎

This subsection is almost identical to Subsection 5.3 of [21]. The main difference is the proof of [claim:1ststep95etavalid] in 8.3.1.

We now present our procedure MergeSolutions that walks between two solutions. We will maintain two solutions \((\mathcal{H}, \mathcal{D}\mathcal{H}, f, u)\) and \((\mathcal{H}', \mathcal{D}\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\).

8.3.0.1 Initialization.

We binary search on \(f\) in the interval \([1/n^2, 4nM]\) and run the routine \(\mathrm{\small CompleteSolution}_{f, u}(\emptyset, \langle\mathcal{D}_0\rangle)\) (with zero \(u\)). When \(f = 1/n^2\), every client \(j \in D\) satisfies \(d(j, S) = 0\) if \(d(j, \mathcal{F}) = 0\) and \(d^2(j, S) \leq d^2(j, \mathcal{F}) + 2/n^2\) otherwise, as each client can open its closest facility by itself when \(\alpha_j = d^2(j, \mathcal{F}) + 2/n^2\). If this solution opens at most \(k\) centers, then it is already a \((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}, \mathcal{D}\mathcal{H}, f, u)\) and \((\mathcal{H}', \mathcal{D}\mathcal{H}', f', u)\) where \((\mathcal{H}, \mathcal{D}\mathcal{H}, f, u)\) opens less than \(k\) facilities and \((\mathcal{H}', \mathcal{D}\mathcal{H}', f', u)\) opens more than \(k\) facilities where \(f \leq f' \leq f + \eta\).

8.3.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}, \mathcal{D}\mathcal{H}, f, u)\) and \((\mathcal{H}', \mathcal{D}\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.

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

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

3. At most 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 S_{\mathsf{free}}\); otherwise if \(f=f'\) there is at most one free facility \(\tilde{i}\in S_{\mathsf{free}}\) so that \(u(\tilde{i}) \neq u'(\tilde{i})\) and \(|u(\tilde{i}) - u'(\tilde{i})| \leq \eta\). (We don’t consider \(\pi, \pi'\) defined in the next bullet as difference parameters.)

4. The probability that each client is directly connected to any facility is almost the same in both solutions, i.e., we have that for every client \(j\), \(|\pi(j) - \pi'(j)| \leq \eta/M\), where \(\pi, \pi'\) are the probabilities of direct connections for the collection of copies \(\mathcal{D}_{p-1}\) and \(\mathcal{D}'_{p-1}\), respectively.

5. \(\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

1. either find a solution of \(\eta\)-valid sequences that opens exactly \(k\) regular facilities and \(O(\log n)\) free facilities, or

2. 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^3)\), 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)\) free facilities. Combined with [thm:robust], this proves our main bicriteria guarantee. Indeed, we then have a solution consisting of open centers \(S\) with \(|S| = k + O(\log n / \varepsilon^3)\) such that, by [thm:robust], \[\begin{align} \sum_{j\in D} w(j) \cdot 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 non-zero distance is at least \(1\)). The main theorem 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 four subroutines, presented in the subsequent subsections respectively. Each of the first three subroutines achieves (1) or (2) without condition [condition-iv] above (which ends the whole procedure or lets us move on to the fourth subroutine), or sets up the stage for the next subroutine. The third subroutine always achieves (1) or (2) without condition [condition-iv]. The fourth subroutine, then achieves (2), for those not achieving (1), by ensuring condition [condition-iv].

8.3.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})\).

8.3.1.1 Outcome of this section of MergeSolutions.

Let \[\mathcal{H}'', \mathcal{D}\mathcal{H}'' = \mathrm{\small CompleteSolution}_{(f,u)}((H'_1, \ldots, H'_p), (\mathcal{D}_0, \ldots, \mathcal{D}_{p-1}))\,,\] which is a valid call by [claim:1ststep95etavalid]. If \(\mathcal{H}''\) opens \(k\) regular facilities, we return this solution. Otherwise, as \((\mathcal{H}, \mathcal{D}\mathcal{H}, f,u)\) and \((\mathcal{H}', \mathcal{D}\mathcal{H}', f', u')\) sandwich \(k\), we must have that either \((\mathcal{H}, \mathcal{D}\mathcal{H}, f, u)\) and \((\mathcal{H}'', \mathcal{D}\mathcal{H}'', f, u)\), or \((\mathcal{H}'', \mathcal{D}\mathcal{H}'', f, u)\) and \((\mathcal{H}', \mathcal{D}\mathcal{H}', f', u')\) sandwich \(k\).

On the one hand, if \((\mathcal{H}'', \mathcal{D}\mathcal{H}'', f, u)\) and \((\mathcal{H}', \mathcal{D}\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 exactly one parameter differs by at most \(\eta\). However, the new collection of copies \(\mathcal{D}''_p\) and \(\mathcal{D}'_p\) may not satisfy condition [condition-iv] regarding the bounds on the probabilities of direct connections (\(\pi\)). In this case, we proceed to 8.3.4 of phase \(p\).

On the other hand, if \((\mathcal{H}, \mathcal{D}\mathcal{H}, f,u)\) and \((\mathcal{H}'', \mathcal{D}\mathcal{H}'', f, u)\) sandwich \(k\), then these two solutions have the exact same parameters \(f\) and \(u\), have \(\mathcal{D}_{r}=\mathcal{D}''_{r}\) for every \(r\leq p-1\), 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 \(\mathcal{D}\mathcal{H}''\) to \(\mathcal{D}\mathcal{H}'\)) and proceed to the next step.

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

We are given two solutions \((\mathcal{H}, \mathcal{D}\mathcal{H}, f,u)\) and \((\mathcal{H}', \mathcal{D}\mathcal{H}', f, u)\) on the same parameters \(f,u\) that sandwich \(k\). It is guaranteed that \(\forall r\leq p-1, \mathcal{D}_{r}=\mathcal{D}'_{r}\) in the history of collections of copies. 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)\) and \(\mathcal{D}_1,\dots, \mathcal{D}_{p-1}\) are clear from the context, we omit them in the following when, e.g., using the definition of \(\eta\)-openability, and calling the procedures CompleteSolution and CompleteSequence.⁷

This subsection proceeds similarly to the approach in [17].

8.3.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}\), denoted \(\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 claim.

8.3.2.2 Outcome of this section of MergeSolutions.

We distinguish three cases.

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

8.3.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 \(\mathrm{\small CompleteSolution}(H_1, H_2, \ldots, H^{\text{before}}_p)\) and \(\mathrm{\small 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 kinds 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 \(\mathrm{\small CompleteSolution}(H_1, H_2, \ldots, H^{\text{before}}_p)\) is equivalent to \(\mathrm{\small 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 can thus proceed to 8.3.4 of phase \(p\), 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\).

8.3.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 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 \(\mathrm{\small CompleteSolution}(H_1, \ldots, H^{\text{before}}_p)\) and \(\mathrm{\small 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, \(\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 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\).

8.3.3 Removing Extra Facilities↩︎

This subsection proceeds almost identically to that of [17] (except that we will go to the next subroutine instead of finishing phase \(p\)). As in 8.3.2, the parameters \((f,u)\) and \(\mathcal{D}_1,\dots, \mathcal{D}_{p-1}\) are clear from the context, and we will omit them when using the definition of \(\eta\)-openability and calling the function CompleteSolution.

The input to this step consists of two solutions \(\mathrm{\small CompleteSolution}(H_1, \ldots, H^{\text{before}}_p)\) and \(\mathrm{\small 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 two potential 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:

8.3.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:3rdstep95outcome].

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

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

8.3.3.4 Case 2.

Otherwise, if the two solutions \(\mathrm{\small CompleteSolution}(H_1, H_2, \ldots, H^{\text{after}}_p)\) and
\(\mathrm{\small 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}_*\) was added at the end of \(H^{\text{after}}_p\) with \(u(\tilde{i}_*) = 0\). Note that the solution \(\mathrm{\small CompleteSolution}(H_1, H_2, \ldots, H^{\text{after}}_p)\) is equivalent to \(\mathrm{\small CompleteSolution}(H_1, H_2, \ldots, H^0_p)\) if we set \(u(\tilde{i}_*) = 10M\) (except for the presence of \(\tilde{i}_*\)). We can thus do a 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}_*\)). We can thus proceed to 8.3.4 of phase \(p\), since they have the first \(p\), instead of \(p-1\), \(\eta\)-valid sequences in common and exactly one difference parameter.

8.3.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 \(\mathrm{\small CompleteSolution}(H_1, H_2, \ldots, H^r_p)\) and \(\mathrm{\small 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}_p\). 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}_p\) (which clearly maintains that it 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 can proceed to 8.3.4 of phase \(p\), since they have the first \(p\), instead of \(p-1\), \(\eta\)-valid sequences in common and exactly one difference parameter.

8.3.4 Making the Probability of Direct Connects Almost Identical↩︎

The input of this step consists of solutions \((\mathcal{H}, \mathcal{D}\mathcal{H}, f, u)\) and \((\mathcal{H}, \mathcal{D}\mathcal{H}', f', u')\) such that \((f,u)\) and \((f',u')\) have at most one difference parameter. Moreover, the two solutions sandwich \(k\). However, \(\mathcal{D}_p\) and \(\mathcal{D}'_p\) are simply obtained from \(\mathcal{D}_{p-1}\) and \(\mathcal{D}'_{p-1}\) respectively by running CompleteCopies. In this step, we will walk between the collection of copies \(\mathcal{D}_{p}\) and \(\mathcal{D}'_p\) following 5 so that condition [condition-iv] holds for phase \(p\).

We follow the procedure below twice. Suppose that \(S=\cup_{r\in [p]} H_r\). Suppose that \(d, d'\) are the pairwise distances induced by \(u, u'\). One run is to transform \((\mathcal{D}_p, \mathcal{D}_{p-1}, d)\) towards \((\mathcal{D}'_p,\mathcal{D}'_{p-1}, d')\). The other run is to transform \((\mathcal{D}'_p, \mathcal{D}'_{p-1}, d')\) towards \((\mathcal{D}_p, \mathcal{D}_{p-1}, d)\).

8.3.4.1 Transforming connection types in \((\mathcal{D},\tilde{\mathcal{D}}, d)\) towards \((\mathcal{D}', \tilde{\mathcal{D}}', d')\)

Repeat the following.

• Find a copy \(j\) in \(\mathcal{D}\) and \(j'\) in \(\mathcal{D}'\) such that

  • \(j\) and \(j'\) correspond to the same client,

  • \(j\) and \(j'\) are both of type \(IC\) or \(A\) in \(\tilde{\mathcal{D}}\) and \(\tilde{\mathcal{D}}'\),

  • \(j\) is of type \(DC\) in \(\mathcal{D}\) while \(j'\) is of type \(IC\) or \(A\) in \(\mathcal{D}'\).

  If no such client exists, then stop the process.

• Update \(\mathcal{D}\) by changing \(j\) to \(IC\) type and update the \(\alpha\)-value \(\alpha_j\) to \((1-\delta)d^2(j,S)\).

After running the two procedures, we can obtain a sequence of solutions \(\tilde{\mathcal{D}}_0=\mathcal{D}_p, \tilde{\mathcal{D}}_1, \ldots, \tilde{\mathcal{D}}_m\) by the first run, while another sequence of solutions \(\tilde{\mathcal{D}}'_0=\mathcal{D}'_p, \tilde{\mathcal{D}}'_1, \ldots, \tilde{\mathcal{D}}'_{m'}\) is obtained by the second run. First, we show that all \(\tilde{\mathcal{D}}_r\) and \(\tilde{\mathcal{D}}'_{r}\) are valid collections of copies.

Second, we show that \(\tilde{\mathcal{D}}_m\) and \(\tilde{\mathcal{D}}'_{m'}\) are close enough so that condition [condition-iv] holds for them.

8.3.4.2 The outcome of this section of MergeSolutions.

Because the previous subroutines successfully aligned the opened facilities up to phase \(p\), we have \(H'_i = H_i\) for all \(i \le p\). We can thus obtain a sequence of solutions as follows: \[\begin{align} \mathcal{H}^{(0)}, \mathcal{D}\mathcal{H}^{(0)} &= \mathrm{\small CompleteSolution}_{f,u}((H_1,\dots, H_p), (\mathcal{D}_0, \dots, \mathcal{D}_{p-1}, \tilde{\mathcal{D}}_0))\\ \mathcal{H}^{(1)}, \mathcal{D}\mathcal{H}^{(1)} &= \mathrm{\small CompleteSolution}_{f,u}((H_1,\dots, H_p), (\mathcal{D}_0, \dots, \mathcal{D}_{p-1}, \tilde{\mathcal{D}}_1))\\ \hfill \vdots \hfill \\ \mathcal{H}^{(m)}, \mathcal{D}\mathcal{H}^{(m)} &= \mathrm{\small CompleteSolution}_{f,u}((H_1,\dots, H_p), (\mathcal{D}_0, \dots, \mathcal{D}_{p-1}, \tilde{\mathcal{D}}_m))\\ \mathcal{H}^{(m+1)}, \mathcal{D}\mathcal{H}^{(m+1)} &= \mathrm{\small CompleteSolution}_{f',u'}((H_1,\dots, H_p), (\mathcal{D}'_0, \dots, \mathcal{D}'_{p-1}, \tilde{\mathcal{D}}'_{m'}))\\ \mathcal{H}^{(m+2)}, \mathcal{D}\mathcal{H}^{(m+2)} &= \mathrm{\small CompleteSolution}_{f',u'}((H_1,\dots, H_p), (\mathcal{D}'_0, \dots, \mathcal{D}'_{p-1}, \tilde{\mathcal{D}}'_{m'-1}))\\ \hfill \vdots \hfill \\ \mathcal{H}^{(m+m'+1)}, \mathcal{D}\mathcal{H}^{(m+m'+1)} &= \mathrm{\small CompleteSolution}_{f',u'}((H_1,\dots, H_p), (\mathcal{D}'_0, \dots, \mathcal{D}'_{p-1}, \tilde{\mathcal{D}}'_0)) \end{align}\]

It is clear that the first one and the last one sandwich \(k\). Therefore, there must be two consecutive solutions in the above sequence that sandwich \(k\). We can discuss two cases then.

8.3.4.3 Case 1.

If \((\mathcal{H}^{(m)},\mathcal{D}\mathcal{H}^{(m)})\) and \((\mathcal{H}^{(m+1)}, \mathcal{D}\mathcal{H}^{(m+1)})\) sandwich \(k\), according to [clm:validD-duringwalk] and [clm:closeD-afterwalk], we can proceed to phase \(p+1\) of MergeSolutions with these two solutions. Since the connection types in \(\tilde{\mathcal{D}}_m\) and \(\tilde{\mathcal{D}}'_{m'}\) are perfectly aligned, the only difference between the two solutions is the original difference parameter (either \(f \neq f'\) or \(u \neq u'\)). They share the first \(p\) sequences in common, satisfy condition [condition-iv], and differ by exactly one difference parameter.

8.3.4.4 Case 2.

Suppose that \((\mathcal{H}^{(r)},\mathcal{D}\mathcal{H}^{(r)})\) and \((\mathcal{H}^{(r+1)},\mathcal{D}\mathcal{H}^{(r+1)})\) sandwich \(k\) for some \(r\neq m\). Without loss of generality, we assume that \(r < m\) (the other case is symmetric). Suppose \(j\) is the copy whose connection type differs in \(\tilde{\mathcal{D}}_r\) and \(\tilde{\mathcal{D}}_{r+1}\), where \(j\) is of type \(DC\) in \(\tilde{\mathcal{D}}_r\) while \(j\) is of type \(IC\) in \(\tilde{\mathcal{D}}_{r+1}\). We observe that condition [condition-iv] only fails on the client corresponding to copy \(j\).

We can further split this copy \(j\) by choosing \(\mu'\in [0,\mu(j)]\) such that one copy \(j^{(1)}\) has \(IC\) type with probability \(\mu'\) and the other copy \(j^{(2)}\) has \(DC\) type with probability \(\mu(j)-\mu'\). We denote the resulting collection of copies by \(\tilde{\mathcal{D}}^{\rm sp}(\mu')\). According to [clm:validD-duringwalk], because \(\tilde{\mathcal{D}}_r\) and \(\tilde{\mathcal{D}}_{r+1}\) are valid collections of copies, any interpolation \(\tilde{\mathcal{D}}^{\rm sp}(\mu')\) is a valid collection of copies. Therefore, we can obtain a solution \((\mathcal{H}^{\rm sp}(\mu'), \mathcal{D}\mathcal{H}^{\rm sp}(\mu'))\) by calling \(\mathrm{\small CompleteSolution}_{f,u}((H_1,\dots, H_p), (\mathcal{D}_0, \dots, \mathcal{D}_{p-1}, \tilde{\mathcal{D}}^{\rm sp}(\mu')))\).

By binary searching \(\mu'\), we can find two values \(\mu_1'\) and \(\mu_2'\) with \(|\mu_1'-\mu_2'|\leq \eta/M\) in polynomial time such that the two solutions \((\mathcal{H}^{\rm sp}(\mu'_1), \mathcal{D}\mathcal{H}^{\rm sp}(\mu'_1))\) and \((\mathcal{H}^{\rm sp}(\mu'_2), \mathcal{D}\mathcal{H}^{\rm sp}(\mu'_2))\) sandwich \(k\).

Because the overall probabilities of direct connections for client \(j\) in the two resultant solutions \((\mathcal{D}\mathcal{H}^{\rm sp}(\mu'_1))_p\) and \((\mathcal{D}\mathcal{H}^{\rm sp}(\mu'_2))_p\) differ by exactly \(|\mu'_1 - \mu'_2|\), and we bounded this difference by \(\eta/M\), these two solutions satisfy condition [condition-iv]. We can then proceed to phase \(p+1\) of MergeSolutions with these two solutions: they have the first \(p\) \(\eta\)-valid sequences in common, utilize identical parameters \((f, u)\), and the newly established parameter between them is exactly the connection probability mass \(\mu'\) (which is not viewed as a difference parameter).

8.4 Analysis of the Robust Algorithm↩︎

In this subsection, we prove [thm:robust], restated as follows, showing that a solution \((\mathcal{H}=(H_1,\dots, H_L), \mathcal{D}\mathcal{H}=(\mathcal{D}_0,\mathcal{D}_1,\dots, \mathcal{D}_L))\) consisting of \(\eta\)-valid sequences and valid collections of copies yields an almost LMP \(\Gamma\)-approximation. This subsection follows the same logic as Subsection 5.4 from [21], the main difference is in establishing the constraints for the Dual Payment, which are different for this new algorithm.

Our analysis closely follows that of LogAdaptiveAlgorithm in 7. However, one key difference is that in the current setting, each client \(j\in D\) may have multiple copies in its collection of copies \(D_j\). Older setting in 7 can be seen as a special case where each client has only one weighted copy. To close the gap between the two settings, for each copy \(\tilde{j} \in \mathcal{D}_L\) in the final collection of copies, it induces a corresponding weighted client with weight \(\mu(\tilde{j})\), which is equal to the probability mass of the copy \(\tilde{j}\) in \(\mathcal{D}_L\). For simplicity, in the rest of this subsection, we will call these induced clients by the copies simply as “clients” and call the clients in \(D\) as “original clients”. We will simply use \(j\) to denote such a client, and use \(oc(j)\) to denote its corresponding “original client”. Also, we will use the terminology clients and copies interchangeably later in the analysis.

We can keep track of how each copy splits in each \(\mathcal{D}_p\) for \(p=1,\dots,L\). In this way, for each copy in \(\mathcal{D}_L\), we can know which copy it belongs to in any \(\mathcal{D}_p\) for each \(p=1,\dots, L\). For simplicity, we will split each copy \(j\) in \(\mathcal{D}_p\) into multiple copies \(j_1, \dots, j_r\) in \(\mathcal{D}_L\) which all belong to it. These \(j_1, \dots, j_r\) will keep their own probability mass in \(\mathcal{D}_p\), but their connection types (\(A, DC, IC\)) and \(\tilde{\alpha}\)-values (if not of type \(DC\)) will be inherited from \(j\) in \(\mathcal{D}_p\). After processing this splitting, we can view \(\mathcal{D}_1,\dots, \mathcal{D}_{L}\) as the history of how each copy in \(\mathcal{D}_L\) evolves from phase \(1\) to phase \(L\). Further, the split process in the definition of “valid collection of copies” (5) can be regarded as turning some copies in each \(\mathcal{D}_p\) from \(DC\) to \(IC\).

One more issue to address is that in the collections of copies, we only maintain \(\alpha\)-values for active or indirectly connected copies. We may lose track of the \(\alpha\)-values of directly connected copies. To handle this issue, we define the \(\alpha\)-value for each copy using the following procedure. We remark that this procedure is only for analysis purposes. The sequence of sequences \(\mathcal{H}\) is enough for us to recover a solution. Hence, we do not need to pay attention to the time complexity of this procedure.

For RecoverAlpha, we have the following remarks:

• In the definition of the process, we will view free facilities and regular facilities equivalently (especially when we run RetrieveNewCopies), but in the later analysis of dual payment and [thm:robust], we still distinguish between them.

• The steps I) and II) in Stage 1 are well-defined since \(H_p\) is \(\eta\)-valid with respect to \(\theta, (H_1,\dots, H_{p-1}), \mathcal{D}_{p-1}\) (4).

• In Step IV) in Stage 1, we have two different ways to define \(R_j\) variables for newly directly connected clients. The reason we need different definitions is that \(h_t\) is only openable with respect to a super-sequence. If a client \(j\) is viewed as directly connected w.r.t. the super-sequence, we define \(R_j\) maximally, by enforcing \(R^2_j=\alpha_j/\rho\), because the \((\rho-1)\alpha_j/\rho\) part of \(\alpha_j\) will no longer be used for opening other facilities (see, the proof 33). Otherwise, we define \(R_j\) as usual, by \(R_j = \sqrt{1-\delta} d(j,h_t)\). This is a new feature for our analysis, compared with [21].

• We can show that RecoverAlpha produces the same solution as \((\mathcal{H}, \mathcal{D}\mathcal{H})\) by comparing RecoverAlpha with RetrieveNewCopies, CompleteCopies and the definition of valid collection of copies (5):

  • In Stage 1, step III) of RecoverAlpha opens facilities in \(H_p\) one by one; and then step IV) of RecoverAlpha updates the \(\alpha\)-values and connection types of all clients in exactly the same way as RetrieveNewCopies.

  • In Stage 2, steps i) and ii) of RecoverAlpha update the \(\alpha\)-values and connection types of all clients in exactly the same way as the second step of CompleteCopies. Step iii) of RecoverAlpha moves some directly connected clients to indirectly connected, which is exactly the same as the split process in 5 under our viewpoint. We shall remark that all copies moved from \(DC\) to \(IC\) are those in \(IC\cup A\) in \(\mathcal{D}_{p-1}\) because of our algorithm in 8.3.4.

Throughout the analysis, we let \(\alpha_j^*\) be the final \(\alpha\)-value of client \(j\in \mathcal{D}_L\). With this RecoverAlpha procedure, we can prove two key lemmas, 33 34, that are analogs of [lem:approx_guarantee_logadaptive,lem:dual-feasibility-logadaptive], respectively. The proofs of them are almost identical to that for LogAdaptiveAlgorithm. There are only three minor changes:

1. 3 of \(\eta\)-openability allows for a facility to be paid up to an \(8n\eta\)-additive difference and ensures every facility not to be overbid up to an \(8n\eta\)-additive difference. The first difference only appears in 33 below, while the second difference only appears in 34. These differences reflect in step I) and II) of stage 1.

2. 4 of an \(\eta\)-valid sequence allows \(h_t\) is \(\eta\)-openable with respect to a super sequence \(\mathcal{H}'\supseteq \mathcal{H}\oplus(\{h_1,\dots,h_{t-1}\})\). This difference appears in the proof of 33, and (the last part of) the proof of 34, an analog of 30. This difference reflects in step I), II) and IV) of stage 1. In particular, a significant difference is in the proof of \(r_{jt}\leq R_j\) (36).

3. 5 of valid collections of copies allows for some directly connected copies to be moved to indirectly connected in Stage 2 of each phase. This difference only appears in (the last part of) the proof of 33 and the (first part of) the proof of 34. This difference reflects in step iii) of stage 2.

The combination of these two lemmas will imply [thm:robust].

Note the weights satisfy \(\sum_{j\in \mathcal{D}_L} w(j) = n\) because each original client has total weight \(1\) over all its copies. In the rest of this subsection, we focus on proving 33 34. Because it is clear in the context we will only consider the copies as the clients, we will assume \(D=\mathcal{D}_L\) in the proof and omit the “original clients” for simplicity. We start by showing the dual payment guarantee with respect to \(\alpha^*\), and then we prove that \((\alpha^*_j-(\rho-1)R^2_j)_{j\in D}\) is a feasible solution to the dual after appropriate scaling.

8.4.1 Dual Payment↩︎

8.4.2 Dual Feasibility↩︎

We prove 34 for dual feasibility. Throughout this part, we will assume that \(\Gamma>\rho+100\sqrt{\varepsilon}\).

The following claim is an analog of [clm:no95overbidding].

Next, we fix an arbitrary regular facility \(i \in \mathcal{F}\) and show that the dual-feasibility for \(i\) is satisfied at the end of the robust algorithm. We say that facility \(i\) is frozen if there exists some open facility \(i_0\) such that \[d(i,i_0)\leq 2\cdot \sqrt{\varepsilon\theta}~.\]

8.4.2.1 Modified \(\alpha\)-values.

For analysis purposes, we introduce a modification of the \(\alpha\)-values for each client (copy) \(j \in \mathcal{D}_L\) in terms of the fixed facility \(i\). We denote this modified \(\alpha\)-value as \(\alpha^c_j\). This modified \(\alpha^c\)-value will be used to write the analog of the dual constraints.

The modifications happen only when we open a facility \(i_0\) that does not freeze \(i\), i.e., when \(d(i,i_0)>2\sqrt{\varepsilon\theta}\). More specifically, at the stage with time \(\theta\), if we open such a facility \(i_0\) that does not freeze \(i\), we keep \(\alpha^c_j\) as the value before opening \(i_0\). That is, we do not increase \(\alpha^c_j\) to \(\tau_j\) if \(j\in B(i_0,\sqrt{\varepsilon\theta})\) (step 1b) for such \(i_0\). In other steps (step 1b for a facility \(i_0\) that freezes \(i\) or step 2), we still modify \(\alpha^c_j\) in the same way we modify \(\alpha_j\).

In the rest of this subsection, we will consider \(\alpha^c\) instead of \(\alpha^*\) in writing the constraints. To ensure that the \(\alpha^*\)-value can ensure the dual feasibility, we first show that \(\alpha_j^c, \alpha_j^*\) are close to each other for every client \(j \in \mathcal{D}_L\).

The following lemma and its proof are identical to 31.

Now, we assume that \[\begin{align} \label{eqn:Dstar-robust} D^*= \{j\in \mathcal{D}_L: \alpha^c_j - (\rho-1)R_j^2 > \Gamma d^2(j,i)\}~. \end{align}\tag{52}\] Then, it suffices to show that \[\begin{align} \sum_{j \in D^*} w(j) \cdot (\alpha_j^c - (\rho-1)R_j^2) \le \hat{f} + 8n\eta + (\Gamma+200\sqrt{\varepsilon}) \cdot \sum_{j \in D^*} w(j) \cdot d^2(j,i)~. \end{align}\]

8.4.2.2 Establishing the constraints.

Up to this point, our arguments closely mirror the dual feasibility proof for the log-adaptive algorithm in 7. However, the remaining part of the proof needs some changes from 7, due to the fact that a facility may be \(\eta\)-openable with respect to a super-sequence instead of the sequence that the algorithm actually executes.

Let \(i_0\) be the facility that freezes \(i\), which opens in (stage 1 of) phase \(p\). Let \(\theta = (1+\varepsilon^3)^{p-1}\) be the active dual parameter at the beginning of phase \(p\). Let \(\mathcal{H}\) be the sequence of opened facilities up to \(i_0\), and let \(\mathcal{H'} \supseteq \mathcal{H}\) be the super-sequence with respect to which \(i_0\) is \(\eta\)-openable.

As in 7, we order the clients in \(D^*\) according to the order they are removed from \(A\) (in the algorithm), breaking ties in the increasing order of their \(\alpha^c\)-values. The only additional rule is for the clients who are in \(A\) right before \(i_0\) is opened in \(\mathcal{H}\); let us call this set of active clients \(A\). By [clm:robustfrozenstay], they are the last several clients in the above order of \(D^*\). Let \(A^F \subseteq A\) be the clients in \(D^*\) who are still active right before \(i_0\) is opened according to the collection of copies induced by \((\mathcal{D}_{p-1}, H'_p)\). Let \(A'=(D^* \cap A)\setminus A^F\) be the clients that become inactive during the extra prefix in \(\mathcal{H}'\). In our ordering for \(D^*\cap A\), we place \(A'\) first and \(A^F\) after, sorting within each group in the increasing order of their \(\alpha^c_j\) values.

Let \(a\) be the index of the last client in \(D^* \setminus A^F\), meaning the clients in \([a+1:s]\) are exactly \(A^F\). For each client \(t \in [a]\), we define \(S^t, DC^t,\) and \(IC^t\) to be the sets of opened facilities, directly connected clients in \([t-1]\), and indirectly connected (or active) clients in \([t-1]\), evaluated at the beginning of the exact stage \(t\) becomes inactive. For each client \(t\in [a+1:s]\), we define \(S^t, DC^t,\) and \(IC^t\) to be the sets of opened facilities, directly connected clients in \([t-1]\), and indirectly connected (or active) clients in \([t-1]\), evaluated right before \(t\) becomes inactive in \(\mathcal{H}'\). Similarly, let \(S^F, IC^F,\) and \(DC^F\) be the set of opened facilities, indirectly connected, and directly connected clients right before \(i_0\) is opened evaluated with respect to the super-sequence \(\mathcal{H}'\).

We remark that the definitions of \(S^F, IC^F, DC^F,\) and \(A^F\) here depend on the super-sequence \(\mathcal{H}'\), while in 7, they only depend on the sequence \(\mathcal{H}\) executed by the algorithm. To complete the proof, we will now apply the unified constraint lemmas ([lem:dual-feasibility-simpler-metric,lem:dual-feasibility-simpler-euclidean,lem:dual-feasibility-simpler-euclidean-special]) evaluated over this super-sequence state.

The two following claims are analogs of Claims [clm:notfrozen] and [clm:frozenstay]. The proof of [clm:robustnotfrozen] is identical to that of [clm:notfrozen]. For [clm:robustfrozenstay], we strengthen the argument to “client \(j\) becomes directly connected” upon “client \(j\) becomes active”.

The following claim is the analog of [clm:3rdconstraint].

For any client \(t \in [s]\) and \(j \in DC^t\), define \[r_{jt} := \sqrt{1-\delta}d(j,S^t).\] Next, we prove that \(r_{jt}\leq R_j\) for every \(t\in [s]\) and \(j\in DC^t\).

General metric↩︎

Now, we will show 34 for \(\rho = 2.5\) and \(\Gamma = 4.9\). To this end, we show that \((IC^t), (DC^t), (\alpha_j^c), (R_j), (r_{jt})\) satisfy all the constraints of [lem:dual-feasibility-simpler-metric] (with the right-hand side of the capacity constraint replaced by \(\hat{f} + 8n\eta\)).

The proof of this is strictly identical to the one in 7, with only two minor adaptations due to the super-sequence: First, the last constraint (the capacity bound) is now directly provided by [clm:robust3rdconstraint] with the updated right-hand side \(\hat{f} + 8n\eta\). Second, we will use 36 to argue that \(r_{jt}\leq R_j\) still holds.

The remaining bounds, including \(\alpha_t^c \le (r_{jt} + d(j,t))^2\), follow exactly the same triangle inequality logic as in the log-adaptive part because \(\alpha_t^c \le (1-\delta)d^2(t, S^t)\) remains valid for all \(t\) (for \(t \in [a+1:s]\), we have \(\alpha_t^c \le \tau_t \le (1-\delta)d^2(t, S^F) = (1-\delta)d^2(t, S^t)\) by the second bullet of \(\eta\)-openability).

Hence, applying [lem:dual-feasibility-simpler-metric] with the adjusted right-hand side, we get: \[\sum_{j\in [s]}w(j)\cdot (\alpha_j^c - (\rho-1)R_j^2)\le \hat{f} + 8n\eta + \Gamma\cdot \sum_{j\in [s]}w(j)\cdot d^2(j,i).\] This completes the proof of 34 for \(\rho=2.5\) and \(\Gamma = 4.9\).

Euclidean metric↩︎

Now, we will show 34 for \(\rho=2\), \(\Gamma=3+\ln 2\), and \(d\) a Euclidean metric. Let \(a\) be the index of the last client that becomes inactive strictly before \(i_0\) is opened in the super-sequence \(\mathcal{H}'\).

The following claim is the analog of [clm:2nd-constraint-euclidean].

We will make several cases according to how \(i\) gets frozen.

8.4.2.5 Case 1: All clients in \(D^*\) become inactive strictly before \(i\) is frozen.

In this case, we have \(a=s\). Let us show that \((IC^t), (DC^t), (\alpha_j^c), (R_j),\) and \((r_{jt})\) satisfy all the constraints of [lem:dual-feasibility-simpler-euclidean], where the right-hand side of the capacity bound is adjusted to \(\hat{f} + 8n\eta\).

First, constraint (1) holds because of 36. Second, since \(a=s\), [clm:robust2nd-constraint-euclidean] implies that constraint (2) is strictly satisfied for all \(t \in [s]\). Finally, [clm:robust3rdconstraint] directly establishes constraint (3) with the relaxed right-hand side \(\hat{f} + 8n\eta\).

Hence, we can seamlessly apply [lem:dual-feasibility-simpler-euclidean] (replacing \(\hat{f}\) with \(\hat{f} + 8n\eta\)) to obtain: \[\sum_{j\in [s]}w(j)\cdot(\alpha_j^c-(\rho-1)R_j^2)\le \hat{f} + 8n\eta + \Gamma\cdot \sum_{j\in [s]}w(j)\cdot d^2(j,i).\] This perfectly matches our target inequality, completing the proof for this case.

8.4.2.6 Case 2: Some clients in \(D^*\) become inactive when \(i\) is frozen.

In this case, we assume that \(i_0\) freezes \(i\) in the super-sequence \(\mathcal{H}'\). Let \(IC^F\) be the set of indirectly connected clients right before \(i_0\) is opened in \(\mathcal{H}'\), and let \(A^F\) be the set of active clients in \(D^*\) right before \(i_0\) is opened (meaning \(a < s\)). We define: \[L := \{j \in IC^F : d(j,i) > \sqrt{\varepsilon\theta}\}.\] Note that \(IC^F\) is disjoint from \(A^F\), meaning \(j \notin IC^F\) for any \(j \in [a+1 : s]\).

We will further discuss two subcases to prove 34.

8.4.2.7 Subcase 2(a): \([a+1:s]\cap B(i_0,\sqrt{\varepsilon\theta})=\emptyset\) or \(\sum_{j\in L}w(j)>\varepsilon\cdot\sum_{j\in [a+1:s]}w(j)\).

We redefine \(DC^s\) and \(IC^s\) as the respectively directly connected and active or indirectly connected client in \([s-1]\) at the beginning of the stage \(s\) becomes inactive. First, let us bound the total increment of \(\alpha_j^c\) when \(i_0\) is opened.

The following claim is the analog of [clm:increment-bound].

Applying [clm:robustnooverbidding] at the beginning of the stage that \(s\) becomes inactive, this lemma implies the following bound: \[\begin{align} w(s) (\alpha^c_s-\rho d^2(s,i)) + \sum_{\ell\in IC^s} w(\ell)(\alpha^c_\ell-\rho d^2(\ell,i)) + \sum_{j\in DC^s} w(j) (r^2_{js}-d^2(j,i)) \\ \leq (\hat{f} + 8n\eta) + \varepsilon\sum_{j\in D^*} w(j)d^2(j,i). \end{align}\]

Because we strictly redefined \(IC^s, DC^s,\) and \(S^s\) according to the sets \(IC, DC, S\) at the beginning of the stage that \(s\) becomes inactive, [clm:robust2nd-constraint-euclidean] perfectly generalizes for \(t=s\). Thus, we can cleanly apply [lem:dual-feasibility-simpler-euclidean-special], replacing \(\hat{f}\) with \((\hat{f} + 8n\eta) + \varepsilon\sum_{j\in D^*}w(j)d^2(j,i)\) and using \(L=\emptyset\) and \(\varepsilon'=0\) (the internal \(\varepsilon\)-parameter in the lemma). This directly gives: \[\sum_{j\in [s]}w(j)\cdot (\alpha_j^c - (\rho-1)R_j^2)\le (\hat{f} + 8n\eta) + (\Gamma + \varepsilon)\cdot\sum_{j\in [s]}w(j)\cdot d^2(j,i).\] This completes the proof for this subcase.

8.4.2.8 Subcase 2(b): \([a+1:s]\cap B(i_0,\sqrt{\varepsilon\theta})\neq \emptyset\) and \(\sum_{j\in L}w(j)\le \varepsilon\cdot\sum_{j\in [a+1:s]}w(j)\).

In this case, we again define \(DC^s = DC^F\) and \(IC^s = IC^F \cup A^F\) as the sets of directly connected and active/indirectly connected clients respectively, right before \(s\) becomes inactive (i.e., right before \(i_0\) is opened evaluated in the super-sequence \(\mathcal{H}'\)). Let \(S^s = S^F\).

We show that \(a\), \(L\), \((IC^t), (DC^t), (\alpha_j^c), (R_j),\) and \((r_{jt})\) strictly satisfy all the constraints of [lem:dual-feasibility-simpler-euclidean-special], with the capacity constraint relaxed by \(8n\eta\). Conditions (a), (b), (c), (1), (2), and (3) of the lemma hold by the same arguments established for Case 1 and the structural definitions of the sets. It remains to rigorously prove that condition (4) is satisfied. Because \(L \subseteq IC^F\), \(IC^F \subseteq [a]\), and \(IC^F \subseteq IC^s\), we have \(L \subseteq [a] \cap IC^s\). The following claim establishes the required algebraic bounds for condition (4).

The following claim is the analog of [clm:4th-constraint-euclidean].

Because [clm:robust4th-constraint-euclidean] perfectly satisfies condition (4) of [lem:dual-feasibility-simpler-euclidean-special] (since \(IC^F \setminus L \cup [a+1:s] = IC^s \setminus L\), the second and third bullets map to the second condition of the lemma, while the fourth bullet maps directly to the third condition), we can apply [lem:dual-feasibility-simpler-euclidean-special].

Replacing the capacity bound \(\hat{f}\) with \(\hat{f} + 8n\eta\) throughout the linear combination evaluated in the lemma, we obtain: \[\sum_{j\in [s]}w(j)\cdot(\alpha_j^c-(\rho-1)R_j^2)\le \hat{f} + 8n\eta + (\Gamma + 200\sqrt{\varepsilon})\cdot\sum_{j\in [s]}w(j)\cdot d^2(j,i).\] This identically satisfies our targeted dual feasibility bound required for the robust algorithm. This completes the proof of 34 for the Euclidean metric, and therefore the proof of [thm:robust].

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

In this section, we give a \(({4}+O(\sqrt{\varepsilon}))\)-approximation algorithm for \((\frac{\zeta}{\log n})\)-stable metric \(k\)-means instances, hence proving 2. We will reuse most of the notation and proofs in [17], and only mention the main differences.

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 \(({4}+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 result 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)})\).

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

9.0.0.2 Partition matroids.

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

9.0.0.3 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 [30] 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 [31], building upon the work of [32], 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.

9.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 5 the cost of a solution is a polynomially bounded integer and, at each improving step, the cost decreases by at least one.

9.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 \(26\) approximation for \(k\)-Means, see [15].

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}\).

9.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 \(26\)-approximation by 37. 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\).

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.

9.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 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\).

9.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^* = 26\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^*))\).

9.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 53 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. 39 bounds the cost of the clients in \(D^*_{\texttt{pure}}\) in \(S\) and 41 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.

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

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

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.

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

9.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 42 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 47 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 3d^2(p, c^*) + 6\text{avg}_{C^*, \text{OPT}} + 3\varepsilon\gamma_{C^*}\] (assuming the set \(\mathcal{B}\) of balls is valid and using 8 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)\).

9.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 algorithm that considers a valid ball set of balls \(\mathcal{B}\). In particular, at this point given that \(W\) is successful and \(\mathcal{B}\) 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 Lemma in [17] that takes as input the local search solution \(S\) and the set \(\Lambda\) of dummy centers.

9.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 46, 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.

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

In the following we slightly deviate from the approach in [17]. More specifically, the constant \(4\) below is a \(5\) in [17]. The rest of the analysis remains however very similar. 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)/{4}\), 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)/{4}\)), 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)/{4}\). The motivation for this reassignment is technical. More specifically, it will be used in Lemma 47 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)\).

9.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 9.3 is successful, the set of balls \(\mathcal{B}\) from 9.4 is valid, and \(\mathcal{Q}\) selected in this section satisfies the properties of 46. In [17] there is an analogue version of the next lemma, where (besides the mentioned difference in the constant \(4\) versus \(5\)) only the sum \(\text{cost}(M_O, \mu_O)+\text{cost}(M_D, \mu_D)\) is upper bounded. For reasons that will become clearer later, we need to bound the two quantities separately.

9.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\).

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

We next focus on a pair \((\widetilde{\mathcal{U}},\widetilde{\mu})\) that satisfies the properties in the claim of Lemma 48. 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.

9.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 48, \(\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 48). 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 9.3 is successful, the set of balls \(\mathcal{B}\) from 9.4 is valid, \(\mathcal{Q}\) selected in 9.6 satisfies the properties of 46, and \(\widetilde{\mathcal{U}},\widetilde{\mu}\) selected in this section satisfies the properties of 48. Again, there is an analogue of the next claim in [17], where however only the sum \(\text{cost}(M'_O, \mu'_O)+\text{cost}(M'_D, \mu'_D)\) is upper bounded. We need to upper bound the two quantities separately.

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 48. 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 48, 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{54}\\ \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{55} \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 49) \[\begin{align} \text{cost}(S^*) & \leq (1+6\sqrt{\varepsilon}) \left({(1-\frac{1}{e})}\text{cost}(M'_O, \mu'_O) +{\frac{1}{e}} \text{cost}(M'_D, \mu'_D) \right) + O(\sqrt{\varepsilon})\cdot \text{opt}\\ & \leq \sum_{p\in D} {4} \text{OPT}_p + O(\sqrt{\varepsilon}) \cdot \text{opt} \end{align}\] where the second inequality is by 47 and [claim:Mprimebounds]. So the proof of 49 in the next subsection is the final step in the proof of 2 (see also 9.11 where we put everything together to prove 2).

9.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}\).

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

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

9.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 48, 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:

9.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 38. 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\).

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

This directly implies 49.

9.11 Putting Everything Together – Proof of 2↩︎

We now turn to proving 2. We first describe the algorithm we analyse.

9.11.0.1 A \(({4}+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 9.1

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

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

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

   2. For each \(Q\in \mathcal{L}_{\text{exp}}\) produced as in 46 (9.6):

      1. For each \(\widetilde{\mathcal{U}}{\in \mathcal{L}_{cheap}}\), with the associated \(\widetilde{\mu}\), produced as in 48 (9.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 (9.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 \(({4}+O(\sqrt{\varepsilon}))\)-approximation with probability at least \((1-\varepsilon)(1-1/n)\). 2 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 \(({4}+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 49. The output of the \(s^*\)-\(D^2\)-Sample process (which is clearly a polynomial time procedure) is a subset of size \(s^*\). By 43, we construct \(\mathcal{L}_{\texttt{bal}}\) in polynomial time, and \(|\mathcal{L}_{\texttt{bal}}| \leq n^{\varepsilon^{-O(1)}}\). Similarly, by 46, we construct \(\mathcal{L}_{\text{exp}}\) in polynomial time and \(|\mathcal{L}_{\text{exp}}| \leq n^{\varepsilon^{-O(1)}}\). Finally, by 48, 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 \(({4}+O(\sqrt{\varepsilon}))\)-approximation follows by 49.

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

Next, 43 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 46 we have that with probability at least \(1-1/n\), the procedure from 46 produces a set of centers \(\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 47 applies. Condition on the event that we obtain a successful \((W, \mathcal{B}, \mathcal{Q})\), then, the procedure from 48 outputs a pair \((\widetilde{\mathcal{U}}, \widetilde{\mu})\) satisfying the properties in the claim of the lemma. 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 9.4 is valid, \(\mathcal{Q}\) selected in 9.6 satisfies the properties of 46, and \(\widetilde{\mathcal{U}},\widetilde{\mu}\) selected in this section satisfies the properties of 48 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.

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

In this section, we prove Theorem 3. Recall that our goal is to compute a \(2+O(\sqrt{\varepsilon})\) approximate solution for \(O(\zeta/\log n)\) stable instances. This improves on the \(4+O(\sqrt{\varepsilon})\) approximation in [21], that however works in the more general metric setting, under the only restriction that the client locations are valid facility locations (while our improvement only works for the Euclidean case). Our algorithm is the same as in [21], and our analysis is very similar too: in the following we will only focus on the differences in the analysis.

10.0.0.1 Notation and Preliminaries.

To simplify the formulas, we let \(\varepsilon\) be the minimum of \(\varepsilon\) and \(\zeta\) in the claim of Theorem 3 and further assume that \(\varepsilon<1/12\) is sufficiently small. With this notation, we present a \((2+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\).

Recall that, by known reductions, we can assume that we are given a set \(F\) of tentative centers, with \(D\subseteq F\), and we have to choose centers from that set. Recall also that non-zero distances are between \(1\) and \(n^2/\varepsilon^2\) according to Lemma 6. We also assume w.l.o.g. that the optimal cost is positive (hence at least \(1\)).

Given a cluster \(C^*\) with center \(c^*\in F\), in the analysis it is convenient to consider the actual center of mass \(\mu(C^*)\) of \(C^*\). Obviously \(\sum_{p\in C^*}d^2(\mu(C^*),p)\leq \sum_{p\in C^*}d^2(c^*,p)\). We will also need the following standard fact.

Given a set of centers \(A\subseteq F\) (not necessarily of cardinality \(k\)), we let \(A_p=d^2(p,A)\) denote the cost associated to \(p\) in the clustering induced by \(A\), and \(\text{cost}(A):=\sum_{p\in D}A_p\) be the total cost of \(A\). The cluster in \(A\) associated with \(c\in A\) is denoted by \(A(c)\). 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 under the restriction \(\text{OPT}\subseteq F\). Let also \(\text{opt}= \text{cost}(\text{OPT})\) be its cost.

10.0.0.2 Overview of our approach.

Our algorithm consists of a few steps which are discussed in the following subsections. The first step (see Section 10.1) is to compute a candidate solution \(S\) via standard local search. We will show that the points in almost all the clusters in \(\text{OPT}\) (excluding a logarithmic size set) pay in \(S\) roughly the same cost as in \(\text{OPT}\). We call the centers \(OPT_{pure}\) of such clusters pure. More precisely, we will show that all the points in \(OPT(c^*)\) with \(c^*\in \text{OPT}_{pure}\) can be reassigned to the same center \(c'\in S\) without increasing too much their cost.

For each remaining impure center \(c\), we identify a core set of clients in \(\text{OPT}(c)\), the leaders, which are sufficiently close to \(c\). We distinguish the impure clusters depending on whether the cost in \(S\) of the leaders is small (core-cheap) or large (core-expensive). The next step (see Section 10.2) is to identify one leader in most of the core-expensive clusters of \(\text{OPT}\). Let \(\Lambda\) be the set of such leaders, and \(\text{OPT}_{exp}\) be the core-expensive centers for which one leader was discovered. One can reassign all the clients \(\text{OPT}(c^*)\) with \(c^*\in \text{OPT}_{exp}\) to the corresponding leader while increasing their cost at most by a factor (roughly) \(2\). Let \(\text{OPT}_{cheap}\) be the remaining impure centers of \(\text{OPT}\). We show that we can reassign all the clients \(OPT(c^*)\) with \(c^*\in \text{OPT}_{cheap}\) to the same center in \(S\) while increasing their cost, again, by at most a factor (roughly) \(2\). The latter two mentioned factors \(2\) are replaced by \(4\) in the analysis in [21].

Consider next the solution that contains the leaders \(\Lambda\) plus the centers in \(S\) to which the points in the clusters of \(\text{OPT}_{pure}\cup \text{OPT}_{cheap}\) are reassigned (which are at most \(|\text{OPT}_{pure}|+|\text{OPT}_{cheap}|\) many). Notice that this is a feasible solution since \(|\Lambda|=|\text{OPT}_{exp}|\). Furthermore, it is (roughly) \(2\)-approximate by the previous discussion. Our final goal is to identify such a solution, or some good enough approximation of it. More precisely, we start with the solution \(\Lambda\cup S\) (which is infeasible), and remove a set \(S'_{disc}\subseteq S\) of \(|\Lambda|\) many carefully chosen centers. This removal is done in two steps. First (see Section 10.3), we remove a set \(\mathcal{Q}\subseteq S\) of centers which are expensive in \(S\). Then (see Section 10.4), we remove other \(\widetilde{\mathcal{U}}\) centers from \(S\), with \(|\widetilde{\mathcal{U}}|=|\Lambda|-|\mathcal{Q}|\).

10.1 Locally Optimal Solution and Pure Clusters↩︎

Let \(S\) be a locally optimal solution output by the following standard local search algorithm. Here, for technical reasons and differently from [33], we use a truncated version of local search: this way we do not need to assume that distances are integers.

10.1.0.1 \(localSearch()\)

\(S \gets\) Arbitrary solution for \(k\)-Means \(S \gets S'\) return \(S\)

We remark that the running time of the above procedure is polynomial. It is standard fact that local search achieves a \(26\) approximation for \(k\)-Means, see [15] (as we will see, the specific constant in not very relevant for us).

We next fix \(S\), and recall that we also fixed a reference optimal solution \(\text{OPT}\). We identify in \(\text{OPT}\) a set of pure clusters as follows. 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. By \(\text{OPT}_{pure}\subseteq \text{OPT}\) we denote the set of the centers of pure clusters (that we also call pure). The centers (resp, clusters) which are not pure are impure. The impure centers of \(\text{OPT}\) are denoted by \(\text{OPT}_{imp}\). We use a similar notation for the centers and clusters of \(S\). 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 \(26\)-approximation by 52. 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\). We will use the following Lemma from [17] (for technical reasons the original proof uses \(25\) as the approximation ratio of the local search solution, however their proof works exactly in the same way with \(26\)).

We next define a replacement cost \(r(p)\) for each client \(p\in D\). 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^*))\). For any impure cluster \(C^*\) with center \(c^*\), we let \(t(c^*)\) be the center of \(S\) closest to \({\mu(C^*)}\), and define \(r(p)=d^2(p,t(c^*))\) for each \(p\in C^*\).

The following lemma is the only one that requires that the solution \(S\) is obtained via 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. We will use the following Lemma from [17].

10.2 Leaders’ Sampling of Leader-Expensive Clusters of \(\text{OPT}\)↩︎

Let \(\text{OPT}_{imp}=\text{OPT}-\text{OPT}_{pure}\) be the impure centers of \(\text{OPT}\). We next partition \(\text{OPT}_{imp}\) according to the following definition. Consider any \(c^*\in \text{OPT}_{imp}\) and set \(C^*=\text{OPT}(c^*)\). We let \(\text{avg}(c^*)\) be the squared distance from \({\mu(C^*)}\) to the \(\varepsilon|C^*|\)-th closest client of \(C^*\) to \({\mu(C^*)}\). We further let \(\text{OPT}_{lead}(c^*)\) be the set of clients of \(C^*\) at a squared distance at most \(\text{avg}(c^*)\) from \({\mu(C^*)}\): we call such clients leaders of the cluster. So \(\text{avg}(c^*)\) is the maximum squared distance from a leader to \({\mu(C^*)}\). Also note that we have, \(|\text{OPT}_{lead}(c^*)|\geq \varepsilon|C^*|\). We thus have that \(\text{avg}(c^*)\) is roughly upper bounded by the average optimal cost of \(C^*\), namely \[\label{eqn:boundavg95eucl} \text{avg}(c^*)\leq {\frac{1}{(1-\varepsilon)|C^*|}\sum_{p\in C^*} d^2(p, \mu(C^*))}\leq \frac{1}{(1-\varepsilon)|C^*|}\sum_{p\in C^*} \text{OPT}_p.\tag{56}\] At the same time \(\text{OPT}_{lead}(c^*)\) contains a constant fraction of the clients in \(\text{OPT}(c^*)\). We say that an impure center \(c^*\in \text{OPT}_{imp}\) is leader-cheap if the total cost in \(S\) of the clients in \(\text{OPT}_{lead}(c^*)\) is less than \(\varepsilon^{5}\text{opt}/\log n\), i.e., \(\sum_{p\in \text{OPT}_{\text{avg}}(c^*)}S_p\leq \frac{\varepsilon^5}{\log n}\text{opt}\). The remaining impure centers are called leader-expensive. We name similarly the respective clusters. Let \(\text{OPT}_{lcheap}\) and \(\text{OPT}_{lexp}\) be the leader-cheap and leader-expensive centers, resp. The next lemma shows that leader-cheap clusters have a small replacement cost. There is a similar lemma in [21] where the factor \(2\) is replaced by \(4\).

Our next goal is to identify one leader for each leader-expensive center. We will actually miss each one of the latter centers with some small probability, however this is not a problem as we will see. The idea is to sample each client \(p\) with probability proportional to its cost \(S_p\) for a sufficiently large number of times: this way it is unlikely to not hit each set \(\text{OPT}_{lead}(c^*)\) of a leader-expensive center \(c^*\). Then one extracts a proper set of leaders by brute force. Here we critically exploit the fact that there is only a logarithmic number of impure centers (by 53), hence of leader-expensive centers. We will use the following result in [21] that works when \(S\) is \(O(1)\) approximate.

We next assume that the event in the above lemma happens, and focus on the corresponding value of \(\Lambda\). For \(c^*\in \text{OPT}_{exp}\), we let \(lead(c^*)\) be the leader in \(\text{OPT}(c^*)\cap \Lambda\). Let us define \(\text{OPT}_{exp}\) and \(\text{OPT}_{miss}\) as in the above lemma. Let also \(\text{OPT}_{cheap}=\text{OPT}_{lcheap}\cup \text{OPT}_{miss}\). We call expensive and cheap the centers in \(\text{OPT}_{exp}\) and \(\text{OPT}_{cheap}\), resp. Similarly for the respective clusters. The next lemma upper bounds the distance of the points in expensive clusters to the corresponding leaders. There is an analogous lemma in [21] where the \(2\) is replaced by a \(4\).

The next lemma shows that there exists a choice of \(|\text{OPT}_{exp}|\) many centers \(S_{disc}\subseteq S\), such that, replacing \(S_{disc}\) with \(\Lambda\) in \(S\) leads to a roughly \(2\) approximate solution. Let \(D_{exp}=\cup_{c^*\in \text{OPT}_{exp}}\text{OPT}(c^*)\) and \(D_{cheap}=D-D_{exp}\). The proof of the following lemma is almost identical to the one of an analogous lemma in [21] (where the \(2\) is replaced by a \(4\)).

10.3 Removal of Expensive Centers of \(S_{disc}\)↩︎

Let us start with the tentative solution \(S\cup \Lambda\) (which is infeasible since it contains more than \(k\) centers). Our goal is to remove a proper subset \(S'_{disc}\) of \(|\Lambda|\) many centers from \(S\), so that \(S\cup \Lambda -S'_{disc}\), which is a feasible solution, is roughly \(2\)-approximate. Ideally we would like to choose \(S'_{disc}=S_{disc}\), so that Lemma 58 directly gives what we need. We will rather find some \(S'_{disc}\) which is comparably good. For a subset \(A\subseteq S\), let \(S_{A}:=S\cup \Lambda-A\). As usual, we let \(S_A(c)\) be the cluster in \(S_A\) associated with the center \(c\in S_A\).

Consider any cluster \(C\) of \(S\cup \Lambda\) with center \(c\in S_{disc}\). For reasons that will become clearer later, if we remove \(c\) (which is then added to \(S'_{disc}\)), we would like to assign all the clients in \(C\) to the same center \(c'\) in \(S\cup \Lambda-\{c\}\). Similarly for the later steps. However, this might be too expensive when the cost of \(C\) in \(S\cup \Lambda\) is large. Our next goal is to identify a subset \(\mathcal{Q}\subseteq S_{disc}\) of centers such that, after removing \(\mathcal{Q}\) from \(S\cup \Lambda\), we enforce in some sense the above property. More precisely, we desire that the cost of the clusters with center in \(S_{disc}-\mathcal{Q}\) in the solution \(S_{\mathcal{Q}}:=S\cup \Lambda-\mathcal{Q}\) is very small compared with \(\text{opt}\). This is done via the following lemma in [33] that is based on the fact that \(\text{cost}(\Lambda\cup S-S_{disc})\leq O(1)\cdot\text{opt}\).

10.4 Removal of Cheap Centers from \(S\)↩︎

Our algorithm considers each choice of \(\mathcal{Q}\in \mathcal{L}_{exp}\). We next assume that there exists at least one such \(Q\) that satisfies the properties in Lemma 59 (event which happens with high probability), and focus on the execution of the algorithm corresponding to the latter \(\mathcal{Q}\). Consider the (possibly infeasible) solution \(S_{\mathcal{Q}}\). Our target reference solution is \(S_{S_{disc}}\), which is roughly \(2\)-approximate. Let \(\mathcal{U}:=S_{disc}-\mathcal{Q}\). Ideally we would like to remove \(\mathcal{U}\) from \(S_{\mathcal{Q}}\), hence getting the desired solution \(S_{\mathcal{Q}\cup \mathcal{U}}=S_{S_{disc}}\). We will rather remove a subset \(\mathcal{U}'\) of the same cardinality, whose removal leads to a solution of cost only marginally larger than the cost of \(S_{S_{disc}}\). In order to simplify our argument, for each candidate set of centers \(\mathcal{U}'\subseteq S-\mathcal{Q}\) with \(|\mathcal{U}'|=|\mathcal{U}|\), we will restrict our attention to suboptimal assignments \(\mu_{\mathcal{U}'}\) of points to centers in \(S_{\mathcal{Q}\cup \mathcal{U}'}\) such that, for each \(c\in \mathcal{U}'\), all the points \(p\in S_{\mathcal{Q}}(c)\) are assigned to the closest center \(t(c,\mathcal{U}')\in S_{\mathcal{Q}\cup \mathcal{U}'}\) to \(c\). The remaining points are assigned by \(\mu_{\mathcal{U}'}\) to the same center as in \(S_{\mathcal{Q}}\). We call consistent an assignment \(\mu_{\mathcal{U}'}\) of this type, and we denote the consistent cost of the corresponding solution by \[cost'(S_{\mathcal{Q}\cup \mathcal{U}'})=\sum_{p\in D}d^2(p,\mu_{\mathcal{U}'}(p)).\] Obviously \(cost'(S_{\mathcal{Q}\cup \mathcal{U}'})\) is an upper bound on the actual cost \(cost(S_{\mathcal{Q}\cup \mathcal{U}'})\) of the solution \(S_{\mathcal{Q}\cup \mathcal{U}'}\). The next lemma shows that, for \(\mathcal{U}'=\mathcal{U}\), the consistent cost is small enough. The first inequality in the following lemma follows directly from an analogous lemma in [21], while the second from 58.3.

Our next goal is to compute a set \(\widetilde{\mathcal{U}}\) of cardinality \(|\mathcal{U}|\) such that \[cost'(S_{\mathcal{Q}\cup \widetilde{\mathcal{U}}})\leq cost'(S_{\mathcal{Q}\cup \mathcal{U}})+O(\sqrt{\varepsilon})\text{opt}.\] This is done via the following lemma from [21].

10.5 Putting Everything Together – Proof of 3↩︎

We next summarize the overall algorithm:

10.5.0.1 \(apxStableEuclideankMeans()\)

\(\mathcal{L}_{apx}\leftarrow \emptyset\) \(S \gets localSearch()\) (see 10.1) Compute \(\mathcal{L}_{lead}\) according to Lemma 56 (see 10.2) Compute \(\mathcal{L}_{exp}\) according to Lemma 59 (see 10.3) Compute \({\cal L}_{cheap}\) according to Lemma 61 (see 10.4) \(\mathcal{L}_{apx}\leftarrow \mathcal{L}_{axp}\cup \{S\cup \Lambda - \mathcal{Q}-\widetilde{\mathcal{U}}\}\) return the cheapest solution \(APX\in \mathcal{L}_{apx}\)

We are now ready to prove Theorem 3.

Acknowledgement↩︎

The authors developed the new algorithm, found evidence that it might improve upon the previous work, and established the \((3+\ln 2+\varepsilon)\)-approximation for the Euclidean case without any help of AI.

For the metric case, the authors developed a \((5+\varepsilon)\)-analysis without any help of AI. The better \((4.9+\varepsilon)\)-analysis was obtained by GPT 5.2/5.4 Pro/Thinking after several rounds of communication between the authors and it. The writeup for the \((4.9+\varepsilon)\)-analysis has been verified and revised by the authors. The two analyses follow quite different routes within the spectral analysis framework and were found by different authors simultaneously.

Aditya Anand, Euiwoong Lee, and Amatya Sharma were supported in part by NSF grant 2236669. Ruiquan Gao was supported by NSF CCF-2112824 and Aviad Rubinstein’s David and Lucile Packard Fellowship. Fabrizio Grandoni was partially supported by the SNF Grants 200021-200731 and 200021-236706. Ernest van Wijland was supported in part by the French PEPR integrated projects EPIQ (ANR-22-PETQ-0007).

11 Proof of Theorems [thm:main95euclidean] and [thm:main95metric]↩︎

12 Missing Proofs in 5↩︎

12.1 Proof of [lem:stoc-key-lemma]↩︎

In this appendix, we prove [lem:stoc-key-lemma] for the sake of completeness. Our proof works for the more general case where \(x_i\in \mathbb{R}^{\mathsf{dim}}\) for some \(\mathsf{dim}\geq 1\).

12.1.1 Preliminaries on weighted \(k\)-means↩︎

We define the weighted versions of Euclidean \(k\)-Means that will be useful later in the proof.

One technical lemma we will use later is a characterization of the optimal cost of a cluster consisting of a set of clients among all possible choices of the center.

12.1.2 The proof↩︎

Let \(a\in [s-1]\) be the smallest integer such that \(a\in DC^t\) for some \(t\in [s]\). If no such \(a\) exists, we set \(a=s\) (note that otherwise \(a<s\)). We will prove by induction on \(a\).

12.1.2.1 Inductive hypothesis.

For any \(s'\in [s]\) and any client \(j\in [s]\), we let \(w_{s'}(j)=\sum_{t\in [s']} \beta_t\eta_{tj}\). The inductive hypothesis is then: for any \(s'\in [s]\), we have \[\begin{align} \label{eqn:stoc-inductive-hypothesis} \sum_{t\in [s']} \beta_t \sum_{j\in DC^t} w(j) \sum_{\ell\in IC^t\cup \{t\}} p_{tj\ell} \cdot \|x_j-x_\ell\|^2 \leq \text{cost}_{w_{s'}}([a:s])~. \end{align}\tag{57}\] Because of recurrence relation of \(\beta_t\) (19 ), we have \(w_s(j)=w(j)\) for each \(j\in [s]\). We also have \(P_{j\ell} = \sum_{t\in [s]: j\in DC^t} \beta_t p_{tj\ell}\). Therefore, the LHS of 57 equals \(\sum_{j\in [s]} w(j) \sum_{\ell\in [s]} P_{j\ell} \cdot \|x_j-x_\ell\|^2\). [lem:stoc-key-lemma] then follows from the inductive hypothesis for \(s'=s\): \[\begin{align} \sum_{j<\ell} P_{j\ell} \cdot \|x_j-x_\ell\|^2 \leq \text{cost}_w([a:s]) \leq \text{cost}_w([s]) \leq \sum_{j\in [s]} w(j) \|x_j\|^2~, \end{align}\] where the last inequality follows from the fact that \(\sum_{j\in [s]} w(j) \|x_j\|^2\) is the total weighted cost of opening the center at \(0^{\mathsf{dim}}\).

12.1.2.2 Base case.

The base case is when \(a=s\). In this case, because \(DC^t=\emptyset\) for any \(t\in [s]\), we have \(p_{tj\ell}=0\) for any \(t,j,\ell\in [s]\), \(j\in DC^t\). The proof of the base case is then straightforward because the LHS of 57 equals \(0\).

12.1.2.3 Inductive step.

We assume that the inductive hypothesis holds for any \(a'>a\) and we will prove that it also holds for \(a\). Let \(t_{a}\) the smallest index \(t\) such that \(a\in DC^t\). Note that we have \(a<t_a\leq s\) because \(DC^a\subseteq [a-1]\) and \(a\in DC^t\) for some \(t\in [s]\).

Consider the case \((\widetilde{DC}, \widetilde{IC}, w)\), where the weights are the same, but for each \(t\in [s]\), we have \(\widetilde{DC}^t = DC^t\setminus \{a\}\) and \(\widetilde{IC}^t = IC^t\cup \{a\}\). We will use \(\tilde{p}_{tj\ell}, \tilde{\eta}_{tj}, \tilde{\beta}_t, \tilde{a}\) to denote the corresponding \(p_{tj\ell}, \eta_{tj}, \beta_t, a\) in the case \((\widetilde{DC}, \widetilde{IC}, w)\). In this case, we have \(\tilde{a}>a\). Let \(\tilde{w}_{s'}(j) = \sum_{t\in [s']} \tilde{\beta}_t \tilde{\eta}_{tj}\) for any \(s'\in [s]\) and any client \(j\in [s]\). Hence, according to the inductive hypothesis, we have \[\begin{align} \sum_{t\in [s']} \tilde{\beta}_t \sum_{j\in \widetilde{DC}^t} w(j) \sum_{\ell\in \widetilde{IC}^t\cup \{t\}} \tilde{p}_{tj\ell} \cdot \|x_j-x_\ell\|^2 \leq \text{cost}_{\tilde{w}_{s'}}([\tilde{a}:s]) \end{align}\] Next, we compare the parameters in the case \((DC, IC, w)\) and the case \((\widetilde{DC}, \widetilde{IC}, w)\).

• For \(t<t_a\), we have \(\widetilde{DC}^t = DC^t\) and \(\widetilde{IC}^t = IC^t\). Therefore, we have \(\tilde{p}_{tj\ell} = p_{tj\ell}\) for any \(j,\ell\in [s]\).

• For \(t\geq t_a\), we have \(DC^t = \widetilde{DC}^t\cup \{a\}\) and \(\widetilde{IC}^t = IC^t\setminus \{a\}\). Therefore, we have the following properties for such \(t\).

  • We have \(\pi_t(j) = \frac{\mathsf{suf}(a)}{\mathsf{suf}(a+1)} \cdot \tilde{\pi}_{t}(j)\) for any \(j\in [a+1:s]\) and \(\pi_t(j)=\tilde{\pi}_t(j)=1\) for \(j\in [a]\).

  • We have \(\eta_{tj} = \frac{\mathsf{suf}(a)}{\mathsf{suf}(a+1)} \cdot \tilde{\eta}_{tj}\) for any \(j\in [a+1:s]\). We have \(\eta_{ta} = 0\), \(\tilde{\eta}_{ta}=1\). We have \(\eta_{tj}=\tilde{\eta}_{tj}=1\) for any \(j\in [a-1]\).

  • We have \(\sigma_t(j) = \frac{\mathsf{suf}(a)}{\mathsf{suf}(a+1)} \cdot \tilde{\sigma}_t(j)\) for any \(j\in [a+1:s]\) and \(\sigma_t(j)=\tilde{\sigma}_t(j)\) for \(j\in [a]\). In particular, we have \(\sigma_t(a) = \mathsf{suf}(a)\).

  • We have \(p_{tj\ell} = \tilde{p}_{tj\ell}\) for any \(j\in [a+1:s]\cap DC^t\) and \(\ell\in IC^t\cup \{t\}\). We have \(p_{ta\ell}=\eta_{t\ell}/\mathsf{suf}(a) = \tilde{\eta}_{t\ell}/\mathsf{suf}(a+1)\) and \(\tilde{p}_{ta\ell} = 0\) for any \(\ell\in [a+1:s]\cap IC^t\). We have \(p_{tj\ell}=\tilde{p}_{tj\ell}=0\) for all other cases.

Using the above properties, we have \(\beta_t = \frac{\mathsf{suf}(a+1)}{\mathsf{suf}(a)} \cdot \tilde{\beta}_t\) for \(t\geq t_a\) and \(\beta_t=\tilde{\beta}_t\) for \(a<t<t_a\). In addition, we have \(\beta_a=1-\sum_{j>a} \beta_j\). It can be verified by the recurrence relation of \(\beta_t\) (19 ). Therefore, we have \(\beta_t \eta_{tj} = \tilde{\beta}_t \tilde{\eta}_{tj}\) for any \(t\in [s]\), \(j\in [a+1:s]\). This implies that \(w_{s'}(j) = \tilde{w}_{s'}(j)\) for any \(s'\in [s]\) and any \(j\in [a+1:s]\).

12.1.2.4 case 1: \(s'<t_a\).

In this case, we have \(DC^t=\widetilde{DC}^t\) for any \(t\in [s']\), and we have \(\beta_t=\tilde{\beta}_t\) for any \(t\in [s']\) except for \(t=a\) (if \(s'\geq a\)). Because \(DC^a=\emptyset\) and 63, we have \[\begin{align} \sum_{t\in [s']} \beta_t \sum_{j\in DC^t} w(j) \sum_{\ell\in IC^t\cup \{t\}} p_{tj\ell} \|x_j-x_\ell\|^2 &= \sum_{t\in [s']} \tilde{\beta}_t \sum_{j\in \widetilde{DC}^t} w(j) \sum_{\ell\in \widetilde{IC}^t\cup \{t\}} \tilde{p}_{tj\ell} \|x_j-x_\ell\|^2\\ &\leq \text{cost}_{\tilde{w}_{s'}}([\tilde{a}:s]) = \text{cost}_{w_{s'}}([\tilde{a}:s]) \leq \text{cost}_{w_{s'}}([a:s]). \end{align}\]

12.1.2.5 case 2: \(s'\geq t_a\).

In this case, we have \(w_{s'}(a) = w(a)\) because \(\eta_{ta}=0\) for any \(t>s'\). Using the above properties, we have \[\begin{align} &\sum_{t\in [s']} \beta_t \sum_{j\in DC^t} w(j) \sum_{\ell\in IC^t\cup \{t\}} p_{tj\ell} \|x_j-x_\ell\|^2 \\ =\;& \sum_{t\in [t_a-1]} \beta_t \sum_{j\in DC^t} w(j) \sum_{\ell\in IC^t\cup \{t\}} p_{tj\ell} \|x_j-x_\ell\|^2 + \sum_{t\in [t_a:s']} \beta_t \sum_{j\in DC^t} w(j) \sum_{\ell\in IC^t\cup \{t\}} p_{tj\ell} \|x_j-x_\ell\|^2 \\ =\;& \sum_{t\in [t_a-1]} \tilde{\beta}_t \sum_{j\in \widetilde{DC}^t} w(j) \sum_{\ell\in \widetilde{IC}^t\cup \{t\}} \tilde{p}_{tj\ell} \|x_j-x_\ell\|^2 + \sum_{t\in [t_a:s']} \beta_t w(a) \sum_{\ell\in IC^t\cup \{t\}} p_{ta\ell} \|x_a-x_\ell\|^2 \\ & + \sum_{t\in [t_a:s']} \frac{\mathsf{suf}(a+1)}{\mathsf{suf}(a)}\tilde{\beta}_t \sum_{j\in \widetilde{DC}^t} w(j) \sum_{\ell\in \widetilde{IC}^t\cup \{t\}} \tilde{p}_{tj\ell} \|x_j-x_\ell\|^2 _{tja}=0}\\ =\;& \frac{w(a)}{\mathsf{suf}(a)} \sum_{t\in [t_a-1]} \tilde{\beta}_t \sum_{j\in \widetilde{DC}^t} w(j) \sum_{\ell\in \widetilde{IC}^t\cup \{t\}} \tilde{p}_{tj\ell} \|x_j-x_\ell\|^2 + \sum_{t\in [t_a:s']} \beta_t w(a) \sum_{\ell\in IC^t\cup \{t\}} p_{ta\ell} \|x_a-x_\ell\|^2 \\ & \frac{\mathsf{suf}(a+1)}{\mathsf{suf}(a)} \sum_{t\in [s']} \tilde{\beta}_t \sum_{j\in \widetilde{DC}^t} w(j) \sum_{\ell\in \widetilde{IC}^t\cup \{t\}} \tilde{p}_{tj\ell} \|x_j-x_\ell\|^2\\ \leq\;& \frac{w(a)}{\mathsf{suf}(a)} \cdot \text{cost}_{w_{t_a-1}}([\tilde{a}:s]) + \frac{\mathsf{suf}(a+1)}{\mathsf{suf}(a)} \cdot \text{cost}_{w_{s'}}([\tilde{a}:s]) + \underbrace{\sum_{t\in [t_a:s']} \beta_t w(a) \sum_{\ell\in IC^t\cup \{t\}} p_{ta\ell} \|x_a-x_\ell\|^2}_{S_3}and w_{\cdot}(j)=\tilde{w}_{\cdot}(j) for j\in [a+1:s]} \end{align}\] Next, we consider the new term \(S_3\). We have \[\begin{align} S_3 &= \sum_{t\in [t_a:s']} \beta_t w(a) \sum_{\ell\in IC^t\cup \{t\}} p_{ta\ell} \|x_a-x_\ell\|^2 \\ &= \frac{\mathsf{suf}(a+1) w(a)}{\mathsf{suf}(a)} \sum_{t\in [t_a:s']} \tilde{\beta}_t \sum_{\ell>a:\ell\in IC^t\cup \{t\}} \frac{\tilde{\eta}_{t\ell}}{\mathsf{suf}(a+1)} \|x_a-x_\ell\|^2~. \\ &= \frac{w(a)}{\mathsf{suf}(a)} \sum_{t\in [t_a:s']} \tilde{\beta}_t \sum_{\ell>a} \tilde{\eta}_{t\ell} \|x_a-x_\ell\|^2 _{t\ell}=0 for \ell\in DC^t} \\ &= \frac{w(a)}{\mathsf{suf}(a)} \sum_{\ell>a} \|x_a-x_\ell\|^2 \cdot \sum_{t\in [t_a:s']}\tilde{\beta}_t \tilde{\eta}_{t\ell} \\ &= \frac{w(a)}{\mathsf{suf}(a)} \sum_{\ell>a} \|x_a-x_\ell\|^2 \cdot \left(\sum_{t\in [s']}\tilde{\beta}_t \tilde{\eta}_{t\ell} - \sum_{t\in [t_a-1]}\tilde{\beta}_t \tilde{\eta}_{t\ell} \right) \\ &= \frac{w(a)}{\mathsf{suf}(a)} \sum_{\ell>a} \|x_a-x_\ell\|^2 \cdot \left(\tilde{w}_{s'}(\ell) - \tilde{w}_{t_a-1}(\ell) \right) \\ &= \frac{w_{s'}(a)}{\mathsf{suf}(a)} \sum_{\ell>a} w_{s'}(\ell) \|x_a-x_\ell\|^2 - \frac{w(a)}{\mathsf{suf}(a)} \sum_{\ell>a} w_{t_a-1}(\ell) \|x_a-x_\ell\|^2~. \end{align}\] According to 64 63, we have \[\begin{align} &\frac{w_{s'}(a)}{\mathsf{suf}(a)} \sum_{\ell>a} w_{s'}(\ell) \|x_a-x_\ell\|^2 + \frac{\mathsf{suf}(a+1)}{\mathsf{suf}(a)} \cdot \text{cost}_{w_{s'}}([\tilde{a}:s]) \\ &\quad \leq \frac{w_{s'}(a)}{\mathsf{suf}(a)} \sum_{\ell>a} w_{s'}(\ell) \|x_a-x_\ell\|^2 + \frac{\mathsf{suf}(a+1)}{\mathsf{suf}(a)} \cdot \text{cost}_{w_{s'}}([a+1:s]) = \text{cost}_{w_{s'}}([a:s])~. \end{align}\] Note that \(\sum_{t>a} w_{t_a-1}(t) \|x_a-x_t\|^2\) is the cost of the cluster \([a+1:s]\) with the center at \(x_a\). Therefore, we have \(\sum_{t>a} w_{t_a-1}(t) \|x_a-x_t\|^2 \geq \text{cost}_{w_{t_a-1}}([a+1:s])\). Therefore, putting \(S_3\) back to the earlier inequality, we have completed the inductive step: \[\begin{align} \sum_{t\in [s']} \beta_t \sum_{j\in DC^t} w(j) \sum_{\ell\in IC^t\cup \{t\}} p_{tj\ell} \|x_j-x_\ell\|^2 \leq \text{cost}_{w_s'}([a:s])~. \end{align}\]

12.2 Proof of [lem:dual-feasibility-simpler-euclidean-special]↩︎

In the rest of this subsection, we will prove the above lemma. We will use the same set of non-negative parameters \(p_{tj\ell}\) and \(\beta_t\) as in the proof of [lem:dual-feasibility-simpler-euclidean] (in terms of \((w(t))_{t\in [s]}\), \((DC^t)_{t\in [s]}\) and \((IC^t)_{t\in [s]}\)). One property we will use about the \(\beta\)-parameters is the following:

We also consider adding one term that is no greater than \(0\) to the LHS of the third constraint [constr:euclidean-special-third] for each \(t\in [a]\cup\{s\}\), using the same parameters \(p_{tj\ell}\) as in the proof of [lem:dual-feasibility-simpler-euclidean]. More specifically, we will consider the term on the LHS of 58 in the following lemma. This is only slightly different from the term we considered in the proof of [lem:dual-feasibility-simpler-euclidean], with a \((1+2\varepsilon)\) factor before each \((r_jt+d(j,\ell))^2\).

By adding the LHS of 58 to the LHS of the third constraint [constr:euclidean-special-third] for each \(t\in [a]\cup\{s\}\), we have the following inequality (for each \(t\in [a]\cup\{s\}\)): \[\begin{align} &\text{LHS of constraint \ref{constr:euclidean-special-third}} \\ &\geq \sum_{\ell\in[t:s]} w(\ell) \alpha^*_t + \sum_{\ell \in IC^t} w(\ell)\alpha^*_\ell + \sum_{j\in DC^t} w(j) \left(r_{jt}^2 + \sum_{\ell\in IC^t\cup\{t\}} p_{tj\ell}(\alpha^*_\ell - (1+2\varepsilon)(r_{jt}+d(j,\ell))^2)\right) \\ &= \sum_{\ell\in[t:s]} w(\ell) \alpha^*_t + \sum_{\ell \in IC^t} w(\ell)\alpha^*_\ell + \sum_{j\in DC^t} w(j) \sum_{\ell\in IC^t\cup \{t\}} p_{tj\ell}\alpha^*_\ell \\ & - \sum_{j\in DC^t} w(j) \left(2\varepsilon r_{jt}^2 + \sum_{\ell\in IC^t\cup\{t\}} p_{tj\ell} (1+2\varepsilon)(2r_{jt}d(j,\ell)+d^2(j,\ell))\right)p_{tj\ell} = 1} \\ &= \sum_{\ell\in IC^t\cup\{t\}} \eta_{t\ell} \alpha^*_\ell - \sum_{j\in DC^t} w(j) \left(2\varepsilon r_{jt}^2 + \sum_{\ell\in IC^t\cup\{t\}} p_{tj\ell} (1+2\varepsilon)(2r_{jt}d(j,\ell)+d^2(j,\ell))\right) } \\ &\geq \sum_{\ell\in IC^t\cup\{t\}} \eta_{t\ell} \alpha^*_\ell - \sum_{j\in DC^t} w(j) \left(2\varepsilon R_j^2 + \sum_{\ell\in IC^t\cup\{t\}} p_{tj\ell} (1+2\varepsilon)(2R_jd(j,\ell)+d^2(j,\ell))\right) \leq R_j} \end{align}\]

Despite not having the above inequality of \(t\in [a+1:s-1]\), because we have \(\beta_t=0\) for any \(t\in [a+1:s-1]\) (65), we have the following inequality for any \(t\in [s]\): \[\begin{align} \beta_t\cdot \left[\sum_{\ell\in IC^t\cup\{t\}} \eta_{t\ell} \alpha^*_\ell - \sum_{j\in DC^t} w(j) \left(2\varepsilon R_j^2 + \sum_{\ell\in IC^t\cup\{t\}} p_{tj\ell} (1+2\varepsilon)(2R_jd(j,\ell)+d^2(j,\ell))\right)\right] \\ \leq \beta_t \cdot \left(\hat{f} + \sum_{\ell \in IC^t \cup [t:s]} \rho w(\ell)d^2(\ell,i) + \sum_{j\in DC^t} w(j) d^2(j,i)\right)~. \end{align}\] Summing over all \(t\in [s]\), because of [obv:beta-eta-sum], \(\rho=2\), \(\Phi=\sum_{\ell} w(\ell)d^2(j,\ell)\), and \(\sum_{t\in [s]} \beta_t \leq 1\) (7), we have \[\begin{align} \sum_{\ell\in [s]} w(\ell)\alpha^*_\ell - \sum_{t\in [s]} \beta_t \cdot \sum_{j\in DC^t} w(j) \left(2\varepsilon R_j^2 + \sum_{\ell\in IC^t\cup\{t\}} p_{tj\ell} (1+2\varepsilon)(2R_jd(j,\ell)+d^2(j,\ell))\right) \\ \leq \hat{f} + 2\Phi - \sum_{t\in [s]} \beta_t \sum_{j\in DC^t} w(j) d^2(j,i) = \hat{f} + 2\Phi - \sum_{j\in [s]} w(j) d^+_j d^2(j,i)~. } \end{align}\] Therefore, we have \[\begin{align} &\sum_{j\in [s]} w(j) \cdot (\alpha^*_j - (\rho-1)R_j^2) - \hat{f} - 2\Phi\leq - \sum_{j\in [s]} w(j) R_j^2 - \sum_{j\in [s]} w(j) d^+_j d^2(j,i) \\ &+ \sum_{t\in [s]} \beta_t \sum_{j\in DC^t} w(j) \left(2\varepsilon R_j^2 + \sum_{\ell\in IC^t\cup\{t\}} p_{tj\ell} (1+2\varepsilon)(2R_jd(j,\ell)+d^2(j,\ell))\right) \\ &\leq \sum_{j\in [s]} w(j) \left(-(1-2\varepsilon) R_j^2 - d^+_j d^2(j,i) + \sum_{t\in [s]:j\in DC^t} \beta_t \sum_{\ell \in IC^t\cup\{t\}} p_{tj\ell} (1+2\varepsilon)(2R_jd(j,\ell)+d^2(j,\ell))\right)\beta_t =d^+_j \leq 1, by \ref{obv:dplus-vs46-sum-beta} \ref{cor:dplus-upperbound}} \\ &= \sum_{j\in [s]} w(j) \left(-(1-2\varepsilon) R_j^2 - d^+_j d^2(j,i) + \sum_{\ell>j} P_{j\ell} (1+2\varepsilon)(2R_jd(j,\ell)+d^2(j,\ell))\right) } \\ &\leq \sum_{j\in [s]} w(j) \left(- d^+_j d^2(j,i) + (1+2\varepsilon)\sum_{\ell>j} P_{j\ell} d^2(j,\ell) + \frac{(1+2\varepsilon)^2}{1-2\varepsilon} \cdot \left(\sum_{\ell>j} P_{j\ell} d(j,\ell)\right)^2 \right) \\ &\leq (1+2\varepsilon)\Phi + \sum_{j\in [s]} w(j) \cdot \left(-d^+_j d^2(j,i) + \frac{(1+2\varepsilon)^2}{1-2\varepsilon} \cdot \left(\sum_{\ell>j} P_{j\ell} d(j,\ell)\right)^2\right) } \\ &\leq (1+2\varepsilon)\Phi + \sum_{j\in [s]} w(j) \cdot \left(-d^+_j d^2(j,i) + (1+7\varepsilon) \cdot \sum_{\ell>j} d^+_j P_{j\ell} d^2(j,\ell)\right) \\ &\leq (1+9\varepsilon)\Phi + \sum_{j<\ell} w(j) d^+_j P_{j\ell} d^2(j,\ell) - \sum_{j\in [s]} w(j) d^+_j d^2(j,i) } \\ & \leq (1+\ln(2)+9\varepsilon)\Phi~. } \end{align}\]

13 Missing Proofs in 6↩︎

13.1 Proof of 23↩︎

Define \[Q(x)\stackrel{\text{def}}{=}288x^4 - 235\tau x^3 - 1071x - 792 + 1128\tau.\] Differentiating 33 gives \[\label{eq:FQ-derivatives} F'(x)=\frac{x}{2178}Q(x),\qquad Q'(x)=1152x^3-705\tau x^2-1071,\qquad Q''(x)=x(3456x-1410\tau).\tag{60}\] We shall use the rational bounds \[\label{eq:xi-rational-bounds} \frac{22894}{10000}<\tau<\frac{228943}{100000},\qquad \Bigl(\frac{22894}{10000}\Bigr)^3<12< \Bigl(\frac{228943}{100000}\Bigr)^3.\tag{61}\]

We split the interval into three pieces.

From 61 , \[3456-1410\tau > 3456 - 1410\cdot \frac{228943}{100000} = \frac{2279037}{10000} >0.\] Hence \(Q''(x)>0\) on \([1,\tau]\), so \(Q'\) is increasing. Therefore on \([1,\frac{4}{3}]\), \[Q'(x)\le Q'\Bigl(\frac{4}{3}\Bigr) = \frac{4979-3760\tau}{3} < \frac{4979-3760\cdot (22894/10000)}{3} = -\frac{453643}{375} <0.\] Thus \(Q\) is strictly decreasing on \([1,\frac{4}{3}]\).

At the endpoints, \[Q(1)=893\tau-1575 > 893\cdot \frac{22894}{10000} - 1575 = \frac{2347171}{5000} >0,\] while \[Q\Bigl(\frac{4}{3}\Bigr) = \frac{15416\tau-35364}{27} < \frac{15416\cdot (228943/100000)-35364}{27} = -\frac{876839}{337500} <0.\] Therefore there is a unique \(r\in (1,\frac{4}{3})\) such that \(Q(r)=0\). Since \(F'(x)\) has the sign of \(Q(x)\) by 60 , the function \(F\) increases on \([1,r]\) and decreases on \([r,\frac{4}{3}]\). Hence \[\max_{1\le x\le 4/3} F(x)=F(r).\]

We now compare \(F(r)\) with \(F(\frac{4}{3})\). By the mean value theorem, there exists \(\zeta\in (r,\frac{4}{3})\) such that \[Q\Bigl(\frac{4}{3}\Bigr)-Q(r)=Q'(\zeta)\Bigl(\frac{4}{3}-r\Bigr).\] Since \(Q(r)=0\) and \(Q'(\zeta)\le -453643/375\), while \[-Q\Bigl(\frac{4}{3}\Bigr) = \frac{35364-15416\tau}{27} < \frac{35364-15416\cdot (22894/10000)}{27} = \frac{44131}{16875},\] we get \[\frac{4}{3}-r = \frac{-Q(\frac{4}{3})}{-Q'(\zeta)} < \frac{44131/16875}{453643/375}.\] Because \(Q\) is decreasing on \([r,\frac{4}{3}]\) and vanishes at \(r\), one has \[0\le -Q(x)< \frac{44131}{16875}\qquad (r\le x\le 4/3).\] Therefore \[\begin{align} 0\le F(r)-F\Bigl(\frac{4}{3}\Bigr) &= \int_r^{4/3} -F'(x)\,dx = \int_r^{4/3}\frac{x(-Q(x))}{2178}\,dx \\ &\le \frac{4/3}{2178}\cdot \frac{44131}{16875}\cdot \Bigl(\frac{4}{3}-r\Bigr) \\ &< \frac{4/3}{2178}\cdot \frac{44131}{16875}\cdot \frac{44131/16875}{453643/375} \\ &= \frac{3895090322}{1125432995259375} <4\times 10^{-6}. \end{align}\] Next, \[F\Bigl(\frac{4}{3}\Bigr) = \frac{97760\tau}{264627} + \frac{822427}{529254} < \frac{97760\cdot (228943/100000)}{264627} + \frac{822427}{529254} = \frac{793785221}{330783750} <2.39971045.\] Hence \[\max_{1\le x\le 4/3} F(x) < 2.39971045 + 4\times 10^{-6} < \frac{12}{5}.\]

Since \(Q''>0\), the function \(Q\) is convex. We already know that \(Q(\frac{4}{3})<0\). Also, \[Q(2)=1674-752\tau < 1674-752\cdot \frac{22894}{10000} = -\frac{29768}{625} <0.\] A convex function that is negative at both endpoints of an interval is negative throughout that interval, so \(Q(x)<0\) on \([\frac{4}{3},2]\). By 60 , the function \(F\) is therefore strictly decreasing on \([\frac{4}{3},2]\), and hence \[\max_{4/3\le x\le 2}F(x)=F\Bigl(\frac{4}{3}\Bigr)<\frac{12}{5}.\]

Again \(Q''>0\), so \(Q'\) is increasing. Moreover, \[Q'(2)=8145-2820\tau > 8145-2820\cdot \frac{228943}{100000} = \frac{8444037}{5000} >0.\] Thus \(Q'(x)>0\) on \([2,\tau]\), so \(Q\) is increasing there. Since also \[Q(2)=1674-752\tau < 1674-752\cdot \frac{22894}{10000} = -\frac{29768}{625} <0,\] the function \(Q\) can cross zero at most once on \([2,\tau]\), and if it does, it crosses from negative to positive. Therefore \(F'(x)=xQ(x)/2178\) can change sign at most once, from negative to positive, so the maximum of \(F\) on \([2,\tau]\) is attained at one of the endpoints.

At \(x=2\), \[F(2)=\frac{752\tau+3297}{2178} < \frac{752\cdot (228943/100000)+3297}{2178} = \frac{31366571}{13612500} <\frac{12}{5}.\] At \(x=\tau\), using \(\tau^3=12\) and \(\tau^6=144\), we obtain \[F(\tau) = \frac{48\tau^6-47\tau^6-357\tau^3+(564\tau-396)\tau^2+4665}{2178} = \frac{\tau^6+207\tau^3-396\tau^2+4665}{2178} = \frac{221}{66}-\frac{2}{11}\tau^2.\] Hence, from the lower bound in 61 , \[F(\tau) < \frac{221}{66} - \frac{2}{11}\Bigl(\frac{22894}{10000}\Bigr)^2 = \frac{988148573}{412500000} <\frac{12}{5}.\] So \(F(x)<\frac{12}{5}\) throughout \([2,\tau]\) as well.

Combining the three steps proves the lemma.

13.2 Weighted: Full Proof↩︎

In this subsection, we prove [lem:dual-feasibility-simpler-metric]. Define \[a_t \stackrel{\text{def}}{=}d(t,i), \qquad W(S)\stackrel{\text{def}}{=}\sum_{j\in S} w(j) \quad (S\subseteq [s]), \qquad W_t \stackrel{\text{def}}{=}\sum_{u=t}^s w(u).\]

As briefly discussed in the overview in 2.1, the proof has four steps. First, from the inequalities at time \(\alpha_t-\varepsilon\) given in condition (3) of [lem:dual-feasibility-simpler-metric], we derive a family of constraints \(\eqref{eq:Pk-weighted}\) by using the triangle inequality to lower bound the terms involving \(r_{jt}\). Second, we choose coefficients \((\beta)_{t\in[s]}\) so that, in the weighted combination of these constraints, the coefficient of each \(\alpha_t^*\) is exactly \(w(t)\). Third, after collecting terms, we obtain a quadratic form in the distances to the facility \(i\). After rescaling by \(\sqrt{w(t)}\), the quadratic form is represented by an ordinary symmetric nonnegative matrix. Finally, we bound the top eigenvalue of the corresponding matrix by the Collatz–Wielandt Formula, which yields the claimed constant.

13.2.0.1 Step 1: From \((I_t^w)\) to \((P_t^w)\).

For each \(t\in[s]\), condition (3) gives \[\label{eq:Ik-weighted}\sum_{j\in DC^t} w(j)(r_{jt}^2-a_j^2) + \sum_{j\in IC^t} w(j)(\alpha_j^*-\rho a_j^2) + \sum_{j\ge t} w(j)(\alpha_t^*-\rho a_j^2) \le \hat{f}.\tag{62}\] Fix \(t\in[s]\). For each \(j\in DC^t\), conditions (1)–(2) and the triangle inequality imply \[\alpha_t^* \le (r_{jt}+d(j,t))^2 \le (r_{jt}+a_j+a_t)^2 \le r_{jt}^2 + (a_j+a_t)^2 + 2(a_j+a_t)R_j.\] Rearranging and using \(2a_ja_t\le a_j^2+a_t^2\), we obtain \[r_{jt}^2-a_j^2 \ge \alpha_t^* - 2a_t^2 - 3a_j^2 - 2(a_j+a_t)R_j.\] Substituting into 62 , we obtain \[\begin{align} &(W(DC^t)+W_t)\alpha_t^* + \sum_{j\in IC^t} w(j)\alpha_j^* - \hat{f} \notag\\ &\qquad\le 2W(DC^t)a_t^2 + \sum_{j\in DC^t} w(j)\bigl(3a_j^2+2(a_t+a_j)R_j\bigr) + \sum_{j\in IC^t}\rho w(j)a_j^2 + \sum_{j\ge t}\rho w(j)a_j^2. \label{eq:Pk-weighted} \end{align}\tag{63}\]

13.2.0.2 Step 2: Choice of the coefficients.

We now choose nonnegative coefficients \(\beta_1,\dots,\beta_s\) and multiply each inequality \(\eqref{eq:Pk-weighted}\) with \(\beta_t\) and sum over all \({t\in[s]}\) such that the coefficient of each \(\alpha_t^*\) is exactly \(w(t)\).

13.2.0.3 Step 3: The weighted quadratic form.

We now compute the quadratic form.

13.2.0.4 Step 4: Weighted matrix form.

It remains to bound \(Q_\rho^w(x)\) by a multiple of \[\sum_{t\in[s]} w(t)x_t^2.\] For \(x\in\mathbb{R}^s\), write \[\widetilde{x} \stackrel{\text{def}}{=}\bigl(\sqrt{w(1)}\,x_1,\dots,\sqrt{w(s)}\,x_s\bigr).\]

We now fix \(\rho=5/2\) and replace \(\eta_t\) by \(g_t\), since \(\eta_t\le g_t\) by 67. For every \(t\in[s]\) and every \(x\in\mathbb{R}_+^s\), 69 gives \[\label{eq:finalkentry-weighted} \frac{(B\widetilde{x})_t}{\sqrt{w(t)}} \le \Bigl(2A_t + \frac{1}{2} g_t + \frac{2}{3} g_t^2\Bigr)x_t + \frac{2}{3} g_t U_t(x) + \frac{2\beta_t}{3w(t)}\sum_{u\in DC^t} w(u)U_u(x) + \frac{2\beta_t}{3w(t)}\sum_{u\in DC^t} w(u)g_u x_u.\tag{67}\]

Since \(F(x)\) is identical to its unweighted counterpart, from 23, we know that \[\label{eq:functionupperbound-weighted} \max_{x\in[1,\tau]} F(x) < \frac{12}{5}.\tag{74}\]

We can now finish the proof of [lem:dual-feasibility-simpler-metric]. Using 68 with \(\rho=5/2\), we obtain \[\begin{align} &\sum_{t\in[s]} w(t)\bigl(\alpha_t^*-(\rho-1)R_t^2\bigr) - \hat{f} \\ &\qquad\le \frac{5}{2} \sum_{t\in[s]} w(t)a_t^2 + Q_{5/2}^w(a) } \\ &\qquad\le \frac{5}{2} \sum_{t\in[s]} w(t)a_t^2 + \lambda_{\max}(B)\sum_{t\in[s]} w(t)a_t^2 } \\ &\qquad\le \frac{5}{2} \sum_{t\in[s]} w(t)a_t^2 + \left(\max_{x\in[1,\tau]}F(x)\right)\sum_{t\in[s]} w(t)a_t^2 } \\ &\qquad\le \frac{5}{2} \sum_{t\in[s]} w(t)a_t^2 + \frac{12}{5}\sum_{t\in[s]} w(t)a_t^2 } \\ &\qquad= 4.9\sum_{t\in[s]} w(t)a_t^2. \end{align}\] This proves [lem:dual-feasibility-simpler-metric].

References↩︎

---

1. Sometimes this is called the continuous Euclidean \(k\)-Means problem in the literature.↩︎

2. For example, we know that \(\rho=1.96, \Gamma=3.692\) works for Euclidean metrics, and \(\rho\approx 2.430,\Gamma=4.865\) works for general metrics.↩︎

3. The only condition for validity is that \(DC_t \subseteq DC_{t+1}\) for all \(t\), since the algorithm moves clients only from \(IC\) and \(A\) to \(DC\), while never moving from \(DC\) to \(IC\).↩︎

4. These parameters are corresponding to the parameters \(\phi^{(t)}_{j,\ell}\) in [21]. We are using a slightly different notation here.↩︎

5. These parameters are corresponding to the parameters \(\beta_t\) in [21].↩︎

6. The definitions in [21] follow a recursion way. Properties here are either based on the explicit form or only implicitly proved in [21].↩︎

7. When we use CompleteSolution\((H_1,\dots, H_{p-1}, H''_p)\), we are actually referring to CompleteSolution\((\mathcal{H}''=(H_1,\dots, H_{p-1}, H''_p), \mathcal{D}\mathcal{H}''=(\mathcal{D}_0,\dots,\mathcal{D}_{p-1}))\). When we use CompleteSequence\((H''_p)\), we are actually referring to CompleteSequence\((\mathcal{H}''=(H_1,\dots, H_{p-1}), H''_{p}, \mathcal{D}\mathcal{H}''=(\mathcal{D}_0,\dots,\mathcal{D}_{p-1}))\).↩︎