Fully Scalable MPC Algorithms for Euclidean $k$-Center

A subset $S\subseteq P$ is called a $\tau$-independent set ($\tau$-IS) for $P$, if for every $x\ne y\in S$, $\mathrm{dist}(x,y)>\tau$, and we say $S\subseteq {ℝ}^d$ is a $\tau$-dominating set ($\tau$-DS) for $P$, if for every $x\in P$, $\mathrm{dist}(x,S)\le \tau$. A subset $S\subseteq P$ is a $(\tau ,\alpha )$-ruling set ($(\tau ,\alpha )$-RS) for $P$ if $S$ is both a $\tau$-IS and $\alpha$-DS for $P$. A $\tau$-MDS is a $\tau$-DS with the minimum size, denoted as ${\mathrm{MDS}}_{\tau }(P)$. A related well known notion is maximal independent set (MIS), where a $\tau$-MIS is $(\tau ,\tau )$-RS.

It is known (see, e.g., [35]) that any $\tau$-IS of $P$ for $\tau \ge 2\mathrm{OPT}$ must have at most $k$ points. Therefore, an MPC algorithm that computes a $(2\mathrm{OPT},\alpha )$-RS of $P$ would immediately yield $\alpha$-approximation for $k$-Center on $P$. On the other hand, for MDS, since the optimal solution for $k$-Center itself is a candidate for a $\tau$-MDS and has at most $k$ points for $\tau =\mathrm{OPT}$, a $(1+\epsilon )$-approximation to $\tau$-MDS would yield $(1+\epsilon )k$ centers. This relation helps to obtain the desired bi-criteria approximation. Compared with the setting of RS which could only leads to some $O(1)$-approximation, MDS operates on the entire ${ℝ}^d$. This is necessary for the $(1+\epsilon )$ ratio since centers need to be picked from ${ℝ}^d$ instead of only from $P$.

We obtain the following results for RS and MDS, in both low and high dimension. Combining with the above-mentioned reductions, these results readily imply Theorems 1.1, 1.2, and 1.3. All results run in $O({\log }_{s}n)$ rounds which is constant in the typical setup of $s=n^{\sigma }$ for constant . In low dimension, we obtain $(\tau ,(1+\epsilon )\tau )$-RS and a $(1+\epsilon )\tau$-DS whose size is at most $(1+\epsilon )$ times the $\tau$-MDS (which in a sense is a “bi-criteria” approximation). Both results use ${({\epsilon }^{-1}d)}^{O(d)}\cdot \mathrm{poly}(\log (n))$ local space, and $n\mathrm{poly}(d\log n)\cdot {({\epsilon }^{-1}d)}^{O(d)}$ total space. In high dimension, we obtain $(\tau ,({\epsilon }^{-1}\log n/\log \log n)\tau )$-RS, using (ideal) $\mathrm{poly}(d\log n)$ local space and $n^{1+\epsilon }\mathrm{poly}(d\log n)$ total space. These results are summarized in Table 1.

We remark that MIS, RS, and MDS are fundamental yet notoriously challenging problems in MPC. Existing studies on these problems are mostly under (general) graphs or general metric spaces, and they achieve worse bounds than ours, e.g., they need to use a super-constant number of rounds [25, 26, 27, 30, 40, 47], and/or are not fully scalable [10, 27, 31]. However, our results seem to suggest that these problems in Euclidean spaces behave very differently than in graphs/general metrics. On the one hand, we obtain fully-scalable algorithms in both low and high dimension, but on the other hand, our algorithms are only “approximations” to MIS and MDS; for instance, in low dimension, both our RS and MDS results have $(1+\epsilon )$ factor off in the dominating parameter to MIS and MDS. For RS/MIS and MDS without violating dominating parameter, we are only aware of a line of research in distributed computing for growth-bounded graphs, see Schneider and Wattenhofer [51], which indirectly lead to $O({\log }^{\ast }n)$-rounds fully scalable algorithms for MIS/MDS in ${ℝ}^d$ for constant $d$. It is still open to design fully scalable algorithms for MIS, even in 2D, in constant number of rounds. In fact, this problem is already challenging on a 2D input set with diameter $O(\tau )$. Nevertheless, our RS and MDS results suffice for approximations for $k$-Center.

In fact, on the rounded dataset $P$, we find a $(\tau ,\tau )$-RS which is as well $\tau$-MIS on $P$. A standard way of finding MIS in a graph is a greedy algorithm: start with the entire vertex set, and if the current vertex set is nonempty, then find any vertex $x$, add it to the output (MIS) set, remove all vertices that are adjacent to $x$ from the current vertex set, and repeat. In the geometric setting, we can use an improved version where in each iteration we identify a large number of vertices that are independent, instead of only one. Specifically, we partition the entire ${ℝ}^d$ into some $T$ groups ${𝒲}_1,\mathrm{\dots },{𝒲}_T$, such that each group consists of regions that are $\tau$ apart from each other. Furthermore, each region has a bounded diameter $\alpha \tau$ for some $\alpha \ge 1$.⁴⁴4The reader may notice the reminiscence with the widely-used notion of network decomposition in graphs [3, 50]. In that notion, the node set $V$ of the graph is partitioned into $T=O(\log n)$ groups $V_1$, …, $V_T$, such that in the subgraph induced by each $V_i$, each connected component has diameter at most $\alpha =O(\log n).$ For $d=2$, there is a simple way to partition with $T=O(1)$ and $\alpha =O(1)$, as illustrated in Figure 1. For general $d$, we give a partition that achieves $T=O(d)$ and $\alpha =O(d^{1.5})$. See Lemma 3.1 for the precise statement.

A key property of this decomposition is that the MIS computation in each region can be done independently within each group, as regions are $\tau$-separated. This yields an $O(T{\log }_{s}n)$-round MPC algorithm: iterating over the $T$ groups ${𝒲}_1,\mathrm{\dots },{𝒲}_T$, computing an MIS in each region in parallel, and removing points within $\tau$ from selected MIS points before proceeding to the next group. Since each region in any group is of diameter $\alpha \tau$ which is at most $d^{1.5}\tau$, there are at most $O{({\epsilon }^{-1}d)}^{O(d)}$ data points in each region, using the assumption of the rounded instance. Hence, as long as $s\ge \mathrm{\Omega }{({\epsilon }^{-1}d)}^{\mathrm{\Omega }(d)}$, the data points in each region can be stored in a single machine. Each such iteration can be implemented in $O({\log }_{s}n)$ MPC rounds, and thus the overall scheme which has $T$ iterations takes $O(T{\log }_{s}n)$ MPC rounds.

To reduce this to $O({\log }_{s}n)$ rounds, we exploit the locality of the above procedure. Instead of iterating sequentially over groups, each region $R\in {𝒲}_i$ for $1\le i\le T$ directly determines its final subset $R^{\prime }$ by identifying the influence from earlier groups that are within a bounded range $O(i\tau +(i-1)\alpha \tau )\le \mathrm{poly}(d)\cdot \tau$. Since an MIS selection only affects a $\tau$-radius per iteration, and each region has a diameter at most $\alpha \tau$, the cumulative affected area over $i$ iterations remains bounded. Leveraging the rounding property again, the number of relevant regions remains at most $O{({\epsilon }^{-1}d)}^{O(d)}$, allowing all necessary regions to be replicated locally. This ensures that each region can compute its MIS in parallel in $O({\log }_{s}n)$ rounds, provided $s\ge \mathrm{\Omega }{({\epsilon }^{-1}d)}^{\mathrm{\Omega }(d)}\mathrm{poly}(\log n)$.

We start with a weaker local space bound that requires a $2^{\mathrm{\Omega }(d^2)}$ dependence of $d$ in $s$, instead of the claimed $2^{\mathrm{\Omega }(d\log d)}$ bound. Our approach starts with a simple algorithm: partition ${ℝ}^d$ into hypercubes of side-length $\alpha \tau$, where $\alpha \ge 1$ is a parameter to be determined. Then send the data points in each hypercube to a single machine, and solve locally the exact MDS of each hypercube. Taking the union of these MDSs forms the final dominating set. Clearly, this returns a dominating set for the entire dataset, but it may not be of small size. Intuitively, the major gap in this algorithm is when the optimal solution uses points located at the boundary of the hypercubes to dominate the adjacent hypercubes, while the algorithm described here only uses a point to dominate points in a single hypercube.

To address this issue, we observe that bridging the gap between the algorithm’s solution and the optimal solution requires only bounding the number of points in the outer “$\tau$-extended region” of each hypercube $R$, denoted as $U_{\tau }^{\mathrm{\infty }}(R)$, that intersect with the optimal solution. Specifically, it suffices to show that this number is at most an $\epsilon$-fraction of the optimal solution size. To establish this bound, we apply an averaging argument: If we shift the entire hypercube partitioning by some multiple of $O(\tau )$, there must exist a shift that satisfies our requirement, provided that the hypercube side-length $\alpha \tau$ is sufficiently large. This may be visualized more easily in 1D: each $U_{\tau }^{\mathrm{\infty }}(R)$ is simply two intervals of length $\tau$ to the left and right of $R$ (which itself is also an interval, whose length is $\alpha \tau$). Then by shifting a multiple of $O(\tau )$, the two intervals of $U_{\tau }^{\mathrm{\infty }}$, after these shifts, form a partition of the entire $ℝ$. Unfortunately, this simple shifting of hypercubes only leads to $\alpha =2^{O(d)}$, which translates to a $2^{O(d^2)}$ dependence of $d$ in the local space $s$. The main reason for this $\alpha =2^{O(d)}$ bound is that the same point in the optimal solution may belong to up to $2^{O(d)}$ sets $U_{\tau }^{\mathrm{\infty }}(R)$ (for some hypercube $R$).

To further reduce $\alpha =\mathrm{poly}(d)$, which leads to the $d^{O(d)}$ local space bound, we need to employ a more sophisticated geometric hashing. We state this in Lemma 3.1 and ˜3.2. This hash maps each point in ${ℝ}^d$ into some bucket, and we would replace the hypercubes in the above-mentioned algorithm with such buckets. An important property of this bucketing is that any point in ${ℝ}^d$ (hence any point in the optimal solution) can intersect at most $\mathrm{poly}(d)$ number of sets $U_{\tau }^{\mathrm{\infty }}(R)$ over all buckets $R$, instead of $2^{O(d)}$ as in the simple hypercube partition. However, the use of this new bucketing also introduces additional issues. Specifically, since the buckets are of complicated structure, it is difficult to analyze the even more complex set $U_{\tau }^{\mathrm{\infty }}(R)$ for the averaging argument. To this end, we manage to show (in Lemma 3.1, third property) that the union of ${\bigcup }_R{U}_{\tau }^{\mathrm{\infty }}(R)$ is contained in the complement of a Cartesian power of (1D) intervals (e.g., ${([1,2]\cup [3,4])}^d$). This structure of Cartesian power is similar enough to the $U_{\tau }^{\mathrm{\infty }}$ annulus of hypercubes (which may be handled by projecting to each dimension), albeit taking the complement. This eventually enables us to use a modified averaging argument to finish the proof.

Hence, we need to introduce the second important modification to this one-round Luby’s algorithm in order to bypass the $O(\log n)$ factor. Specifically, before running the one-round Luby’s algorithm, we run a preprocessing step to map the data points to the buckets of a geometric hashing. The geometric hashing that we use is the consistent hashing [18] (see Lemma 3.4), and the concrete guarantee in our context is that, each hash bucket has diameter $\mathrm{\ell }:=O(\log n/\log \log n)\tau$, and for any subset $S\subseteq {ℝ}^d$ with diameter at most $O(\tau )$, the number of buckets that $S$ intersects is at most $\mathrm{\Lambda }:=\mathrm{poly}(\log n)$. We pick an arbitrary point in $P$ from each (non-empty) bucket, denoting the resultant set as $P^{\prime }$, and we run the one-round Luby’s on $P^{\prime }$. Clearly, this hashing step only additively increases the dominating parameter by $\mathrm{\ell }=O(\log n/\log \log n)\tau$ which we can afford. At a high level, the use of this rounding is to limit the size of the $O(\tau )$-neighborhood for every point, which is analogue to the degree of a graph. In a sense, what we prove is that one-round Luby’s on graphs with degree bound $\mathrm{poly}(\log n)$ yields an $O(\log n/\log \log n)$-ruling set.

Next, we explain in more detail why this hashing step helps to obtain $O(\log n/\log \log n)\tau$ dominating parameter. This requires us to do a formal analysis to one-round Luby’s algorithm (which did not seem to appear in the literature), and utilize the property of hashing that for every $x\in P^{\prime }$, $|B_{P^{\prime }}(x,\tau )|\le \mathrm{\Lambda }=\mathrm{poly}(\log n)$.

Let $R$ be the resultant set found by running one-round Luby’s on $P^{\prime }$. Fix a point $p\in P^{\prime }$, and we need to upper bound $\mathrm{dist}(x,R)$. To this end, we interpret the algorithm as defining an auxiliary sequence $S=(x_0:=p,x_1,\mathrm{\dots },x_T)$ (where $T$ is also random). Specifically, $S$ is formed in the following process: we start with $i=0$, whenever $x_i$ is not picked into $R$ we define $x_{i+1}$ as the point with the smallest $h$ value in $B_{P^{\prime }}(x,\tau )$, and we terminate the procedure otherwise. Clearly, $\mathrm{dist}(x,R)\le T\tau$, and it suffices to upper bound $T$ (which is a random stopping time). Indeed, a similar plan of re-assignment argument has been employed in [17], but the definition and analysis of $S$ is quite different since they focus on facility location.

Now, a crucial step is to bound the probability of the event that, given the algorithm picks a prefix $(x_0,\mathrm{\dots },x_i)$ of sequence $S$, the algorithm picks some new $x_{i+1}\in P^{\prime }$ instead of terminating. We call this extension probability. Ideally, if we can give this extension probability a universal upper bound , then for $t\ge 1$, the probability that $T>t$ is at most ${\gamma }^t$, which decreases exponentially with respect to $t$. However, this seemingly good bound may not immediately be sufficient, since one still needs to take a union bound over sequences of length $t$, because the sequence $S$ is random. A naïve bound is $n^t$ since each $x_i$ may be picked from any point in $P^{\prime }$, and this is positively large (i.e., ${\gamma }^t$ cannot “cancel it out”.). Moreover, an added difficulty is that the “ideal” case of having universal upper bound $\gamma$ for the extension probability may not be possible in the first place.

Our analysis resolves both difficulties. We show that the extension probability is roughly upper bounded by $|B_{P^{\prime }}(x_i,\tau )|/|{\bigcup }_{j=1}^i{B}_{P^{\prime }}(x_j,\tau )|$ (recalling that $(x_0,\mathrm{\dots },x_i)$ is given). For the union bound, since the hashing guarantees that $|B_{P^{\prime }}(x,\tau )|\le \mathrm{\Lambda }=\mathrm{poly}(\log n)$ for any $x\in P^{\prime }$, we can approximate the size $|B_{P^{\prime }}(x_i,\tau )|$ by rounding to the next power of $2$, denoted as ${\lceil |B_{P^{\prime }}(x_i,\tau )|\rceil }_2$, and this yields only $O(\log \log n)$ possibilities after rounding. For a (fixed) sequence $S^{\prime }:=(x_0,\mathrm{\dots },x_m)$, we define its configuration as the rounded value of $({\lceil |B_{P^{\prime }}(x_1,\tau )|\rceil }_2,\mathrm{\dots },{\lceil |B_{P^{\prime }}(x_m,\tau )|\rceil }_2)$. We can do a union bound with respect to the configuration of the sequences, and there are at most $O{(\log \log n)}^t$ number of configurations for length-$t$ sequences.

Finally, given a sequence $S^{\prime }=(x_0,\mathrm{\dots },x_t)$ with a given configuration (for some $t$), to upper bound the extension probability, we need to analyze ${\prod }_i\frac{|B_{P^{\prime }}(x_i,\tau )|}{|{\sum }_{j=1}^i{B}_{P^{\prime }}(x_j,\tau )|}$, and we show in a key lemma that this is upper bounded by $\exp (-\mathrm{\Omega }(t)\log (t/\log \mathrm{\Lambda }))$. Therefore, we can pick $t=\log n/\log \log n$, apply the union bound and conclude that $\mathrm{Pr}[T>t]\le {(\log \log n)}^t\cdot \exp (-\mathrm{\Omega }(t)\log (t/\log \mathrm{\Lambda }))\le 1/\mathrm{poly}(n)$.

We also provide a (sketch) of how to slightly modify our analysis to show the one-round Luby’s algorithm yields $O(\tau ,O(\log n)\tau )$-RS. As mentioned, this did not seem to appear in the literature. We also provide a tight instance for this one-round Luby’s algorithm.