Near-Linear and Parameterized Approximations for Maximum Cliques in Disk Graphs

Clark, Colbourn, and Johnson [15] gave an elegant $O(n^{4.5})$-time algorithm for finding a maximum clique in a unit disk graph. The algorithm guesses in quadratic time the pair of most distant (diameter) disk centers in a maximum clique and considers the subgraph induced by all the disks whose centers are at most as distant from both candidate diameter centers. The induced graph is co-bipartite, hence computing a maximum clique is equivalent to finding a maximum independent set in the complement bipartite graph. This in turn reduces to finding a maximum matching. By reducing the necessary number of pairs of centers to search for and/or reducing the running time for finding a maximum independent set, the running time has been successively improved to $O(n^{3.5}\log n)$ by Breu [7], $O(n^3\log n)$ by Eppstein [18], $O(n^{2.5}\log n)$ by Espenant, Keil and Mondal [19] and, most recently, to $O(n^{7/3+o(1)})$ by Chan [11] as noted in Keil and Mondal [28]. In particular, the latter result is based on the clever divide and conquer approach by Espenant, Keil and Mondal [19], which considers only a linear number of center pairs, and on the observation that, for each such pair, one can use a small bi-clique cover of size¹¹1The standard $\tilde{O}(\cdot )$ notation hides polylog factors in $n$. $\tilde{O}(n^{4/3})$ of the resulting bipartite graph to compute a maximum matching in time almost linear to the size of cover. Very recently, Tkachenko and Wang [35] gave an algorithm that runs in time $O(n\log n+n{K}^{4/3+o(1)})$, where $K$ is the size of the maximum clique; this is an improvement to the above bound by Chan [11] when $K=o(n)$.

For general disk graphs, the complexity of computing the maximum clique has been a long-standing open problem [3, 19, 22]. Recently, Keil and Mondal [28] considered the case where there are $t$ distinct disk radii and gave an algorithm that runs in polynomial time for every constant $t$. In particular, the running time is $O(2^t{n}^{2t}(f(n)+n^2))$, where $f(n)$ is the time to compute a maximum matching in an $n$-vertex bipartite graph. A maximum matching in a bipartite graph can be computed in time $\tilde{O}(m)$, where $m$ is the number of edges in the graph, using the near-linear time max-flow algorithm in the recent breakthrough in [13] with standard reductions to bipartite matching, e.g., the $O(m\log m)$-time fractional flow rounding technique in [26]; alternatively, it can be computed in time $\tilde{O}(n^2)$ by a combinatorial algorithm [14]. Further, the problem admits a randomized EPTAS²²2EPTAS stands for efficient polynomial-time approximation scheme. with a running time of $2^{\tilde{O}(1/{\epsilon }^3)}{n}^{O(1)}$ and an exact sub-exponential time algorithm [5].

Finally, the problem has also been studied for intersection graphs of other classes of objects in the plane and has been shown to be NP-hard for rays [9], strings [29], ellipses and triangles [3], and for sets of both rectangles and unit disks [6].

The above recent exciting results and the absence of any non-trivial lower bounds for the problem leave some natural open questions.

For improving the currently best $O(n^{7/3+o(1)})$-time bound for unit disk graphs, any algorithm that follows the classical line of approach by Clark, Colbourn, and Johnson [15] would effectively need to reduce the number of candidate disk centers to $o(n)$ or compute a maximum matching in the resulting complement bipartite graph in $o(n^{4/3})$ time. Both of these options seem challenging. On the other hand, the algorithm by Tkachenko and Wang [35] does not bring an improvement for instances with $\mathrm{\Omega }(n)$-size cliques, and dense instances are in fact the hardest for the problem, see for example [5].

In this paper, we show that further improvement is possible, albeit by allowing for $(1-\epsilon )$-approximate solutions and by using randomization. Our first result is as follows:

1. (i)

   For unit disk graphs, there is an algorithm that computes with constant success probability a $(1-\epsilon )$-approximate maximum clique in expected time $\tilde{O}(n/{\epsilon }^2)$; this result is presented in Section 3.

In particular, we first random sample two disk centers inside a constant-size neighborhood of a properly defined grid-cell. We show that these centers can help us identify a small region that contains (the centers of) a $(1-\epsilon )$-approximate maximum clique of the neighborhood, with $\mathrm{\Omega }(\epsilon )$ probability. Moreover, the disks with centers in that small region induce a co-bipartite graph. To save time, our second step, which is given in Section 2, is to actually compute a $(1-\epsilon )$-approximate maximum clique in this co-bipartite graph in time that is almost linear in the number of disks of the graph.

Our second result is for the general disk graph setting. As mentioned above, the EPTAS in [5] runs in $2^{\tilde{O}(1/{\epsilon }^3)}\cdot O(n^4)$ time while, for disk graphs with $t$ distinct radii, the algorithm by Keil and Mondal [28] runs in $O^{\ast }(n^{2t})$ time³³3The standard notation $O^{\ast }(\cdot )$ hides polynomial factors.. This raises the question of whether such exponential dependencies w.r.t. $t$ (on $n$) and $\epsilon$ can be improved or even removed altogether. We show that this is possible by considering parameterized $(1-\epsilon )$-approximations w.r.t. $t$ and $\epsilon$ and, as in the case of unit disk graphs, by using randomization.

1. (ii)

   For disk graphs with $t$ distinct radii, there is an efficient parameterized approximation scheme (EPAS)⁴⁴4An efficient parameterized approximation scheme (EPAS), for some parameter $k$, is a $(1-\epsilon )$-approximation algorithm that runs in $f(k,\epsilon )n^{O(1)}$ expected time [21]. that computes with a constant success probability a $(1-\epsilon )$-approximate maximum clique in $\tilde{O}(f(t)\cdot {(1/\epsilon )}^{O(t)}\cdot n)$ time, for some (exponential) function $f(t)$; this result is presented in Section 4.

Our algorithm is based on the main idea behind the exact algorithm by Keil and Mondal. Suppose that for each of the radii one can guess the leftmost and rightmost disks in an optimal solution. Then, one can find two sets of disks $L$ and $R$ such that $L\cup R$ induces a co-bipartite graph and contains a maximum clique. The algorithm by Keil and Mondal [28] enumerates all possible choices to find the leftmost and rightmost disks for each radius. To speed up the running time for $(1-\epsilon )$-approximation, we bring in three more ideas. First, we can assume w.l.o.g that for disks of radius $r_i$ that appear in the optimal clique, at least $\mathrm{\Theta }(\epsilon k^{\ast }/t)$ of them appear in the max clique, where $k^{\ast }$ is the optimal clique size, which can be approximated by a factor of $1/5$ in $O(n\log n)$ time (see Section 4.2) – otherwise, we can skip disks of radius $r_i$ completely without losing more than $\epsilon$ fraction. Second, for each radius $r_i$ we can randomly sample two disks $a_i$ and $b_i$ – if we repeat this enough times, one such pair is close to, in terms of rank, the leftmost and rightmost disks of radius $r_i$ in the optimal solution. Last, as in the case of unit disks, we compute a $(1-\epsilon )$-approximate clique in the resulting co-bipartite graph using the algorithm from Section 2.

We note here that our result for unit disk graphs is not subsumed by the one for $t$ distinct radii since the latter gives a worse dependency on $\epsilon$ for $t=1$.

Because $|C^{\ast }|=n-|I^{\ast }|$, we have $n-c\le |I^{\ast }|$. We observe that $|I^{\ast }|\ge n/2$, as one of the sides of $\overline{G}$ consists of at least $n/2$ nodes forming an independent set. Then, $n-c\ge n-(1+\epsilon )|C^{\ast }|=n-(1+\epsilon )(n-|I^{\ast }|)=|I^{\ast }|-\epsilon \cdot n+\epsilon |I^{\ast }|\ge (1-\epsilon )|I^{\ast }|$. $◀$

We now focus on obtaining a $(1+\epsilon )$-approximation of $|C^{\ast }|$. For a bipartite graph, Kőnig’s theorem states that the cardinality of the minimum vertex cover is equal to the cardinality of the maximum matching. Thus, denoting by $M^{\ast }$ the maximum matching, it is enough to compute a $(1-\epsilon /2)$-approximation of $|M^{\ast }|$, denoted $m$. Indeed, we can assume $\epsilon \le 1$, and then $(1-\epsilon /2)|M^{\ast }|\le m\le |M^{\ast }|$ implies $|M^{\ast }|\le m/(1-\epsilon /2)\le |M^{\ast }|/(1-\epsilon /2)\le (1+\epsilon )|M^{\ast }|$, so $m/(1-\epsilon /2)$ can be returned as a $(1+\epsilon )$-approximation of $|C^{\ast }|$.

To design the algorithm, we need an efficient data structure for maintaining a dynamic set of points $S\subseteq {ℝ}^2$. Each $s\in S$ has its associated weight $w_s\in ℝ$, and we define a distance function $\delta :{ℝ}^2\times S\to ℝ$ as $\delta (p,s)=w_s+\mathrm{dist}(p,s)$. The query returns, for a given point $x\in {ℝ}^2$, the furthest point in $S$, that is, $\arg {\max }_{s\in S}\delta (x,s)$. This can be seen as a point location query in the furthest site additively weighted Voronoi diagram of $S$.

As explained by Agarwal, Efrat, and Sharir [2, Section 9.1], the problem reduces to maintaining the upper envelope of a set of $n$ totally defined continuous bivariate functions of constant description complexity in three dimensions, under insertions, deletions, and vertical ray shooting queries. The discussion of Agarwal, Efrat, and Sharir [2, Section 9.1] also shows that the complexity of these upper envelopes (and thus of the corresponding farthest-site Voronoi diagrams) is linear. Thus, the result follows from the data structure given in [27, Theorem 8.3] and its subsequent improvement [31]. $◀$

With the above data structure in hand, we can describe the algorithm for computing a $(1-\epsilon )$-approximate matching in $\overline{G}$. As explained above, returning $n$ minus its cardinality gives us the cardinality of a $(1-\epsilon )$-approximate clique in $G$.

Let the given co-bipartite graph be $G(X\cup Y)$. We will call $X$ the left side, and $Y$ the right side. We recall that the algorithm of Hopcroft and Karp [24] proceeds in phases while maintaining the current matching $M$ (initially empty). In the beginning of each phase, we run a breadth-first search to partition the nodes of the graph into layers. We start the search from the unmatched nodes on the left side of the graph that form layer $0$. On even layers of the search, we can traverse any unmatched edge. On odd layers of the search, we can only traverse a matched edge. The search terminates after the first layer containing an unmatched node on the right side of the graph (or when no further layers can be created, meaning that the current matching $M$ is a maximum matching). This gives us the length $k$ of a shortest augmenting path. Next, the algorithm finds a maximal set of node-disjoint augmenting paths of length $k$. This can be implemented by running a depth-first search from the unmatched nodes in the last layer. The depth-first search is only allowed to follow an edge leading to a yet unvisited node in the previous layer, and the paths in its tree must alternate between matched and unmatched edges, as in the breadth-first search.

The breadth-first search and the depth-first search run in time linear in the size of the graph, and the length of a shortest augmenting path increases by at least 2 in every phase. It is well-known that if the shortest augmenting path is of length at least $2/\epsilon$ then the current matching $M$ is a $(1-\epsilon )$-approximate maximum matching. Thus, it is enough to execute $1/\epsilon$ phases to obtain such a matching. However, the number of edges in our graph might be even quadratic in $n$, so to finish the proof we need to show how to implement a single phase in $O(n{\log }^{4}n)$ expected time.

We first explain how to implement the breadth-first search in expected $O(n{\log }^{4}n)$ time. To this end, we maintain an instance of the data structure from Lemma 2, initially containing the set $S$ of all points corresponding to the nodes on the right side of the graph. We set $w_s=-r$ for a point $s\in S$ corresponding to a disk of radius $r$. Then, on every even layer, we iterate over the nodes there, and for each such node $u$ corresponding to a disk $D(c,r)$, we need to iterate over the yet unvisited nodes $v$ such that $(u,v)$ is an edge. This is equivalent to obtaining (and removing) all points $s\in S$ such that $\delta (c,s)>r$, which in turn can be implemented by repeatedly retrieving $\arg {\max }_{s\in S}\delta (c,s)$, checking if $\delta (c,s)>r$, and if so, removing $s$ from $S$ and repeating. On every odd layer, we do not need the data structure, as each node has at most one incident matched edge. The overall number of operations on the data structure is $O(n)$, and the deletions dominate, so the total time for the breadth-first search is as claimed. For the depth-first search, we maintain a separate instance of the data structure from Lemma 2 for each even layer $\mathrm{\ell }$, initially containing the set $S_{\mathrm{\ell }}$ of all points corresponding to its nodes. Then, if we are currently at node $v$ on an odd layer $\mathrm{\ell }$, we need to iterate over the yet unvisited nodes $u$ from layer $\mathrm{\ell }-1$ such that $uv$ is an edge. Let the disk corresponding to $v$ be $D(c,r)$, Then, this is again equivalent to retrieving (and removing) all points $s\in S_{\mathrm{\ell }-1}$ such that $\delta (c,s)>r$, and it can be implemented as above. The overall number of operations on all data structures is $O(n)$, making the total time as required. $◀$

For the case of unit disk graphs, Lemma 3 can be improved. For this, we just need a data structure that plays the role of the one in Lemma 2 and that is optimized for unit disk graphs. The following lemma shows how to do this, adapting a method by Efrat, Itai, and Katz [17] to our setting:

Our data structure is very similar to the structure by Efrat, Itai, and Katz [17, Section 5.1], and we refer the reader to their paper for further details.

Let $𝒢$ be a regular grid with cells of edge length $1/2$. During the preprocessing phase, we locate for every disk center in $S$ the cell of $𝒢$ that contains it. Let $C$ be a nonemtpy grid cell, and let $S_C$ be the unit disks whose centers lie in $C$. Since the diameter of $C$ is less than $1$, the disks in $S_C$ have non-empty intersection (e.g., the center of $C$ must lie in all the disks of $S_C$). We compute the intersection $P_C$ of $S_C$. The intersection $P_C$ is a convex polygon-like structure whose boundary consists of arcs from $S_C$, where every disk from $S_C$ contributes at most one arc. The arcs appear along the upper and the lower boundary of $P_C$ in reverse order of the $x$-coordinates of the corresponding disk centers. Thus, the total complexity of $P_C$ is $O(|S_C|)$, and it can be computed in time $O(|S_C|\log |S_C|)$, using a divide-and-conquer strategy. As in the data structure by Efrat, Itai, and Katz, we compute for each $P_C$ two binary trees that represent the structure of the upper and the lower boundary of $P_C$, respectively, as in the well-known data structure by Overmars and van Leeuwen [33]. In these trees, the leaves correspond to the disks from $S_C$, in reverse $x$-order of their centers. The inner nodes represent the upper (or the lower) boundary of the intersection of the disks in the corresponding subtree. The boundary of each inner node is obtained by intersecting the boundaries of the two children. These boundary curves intersect in a single point. With each inner node, we store the intsersection point of the boundary curves of its children, as well as the part of the boundary curve for this node that does not appear on the parent boundary . The binary trees can be computed in $O(n\log n)$ time, and the total time for the preprocessing phase is thus $O(n\log n)$ (see the paper by Efrat, Itai, and Katz for the details).

A find-query can be answered as follows: given $q\in {ℝ}^2$, we first check if there is a non-empty grid cell $C$ that is far enough from $q$ such that $q$ must lie outside of $P_C$. If such a cell $C$ exists, we report an arbitrary disk from $P_C$ and are done (this takes $O(\log n)$ time, assuming that the nonempty grid cells are stored in a suitable data structure). Otherwise, there is a constant number of remaining candidate cells, and for each such candidate cell $C$, we determine whether $q$ lies outside of $P_C$. This can be done using the tree structures for $C$ to determine whether $q$ lies above or below the upper and the lower boundary of $P_C$, by descending into the binary tree, using the intersection points of the child boundaries as a guidance. If $q$ lies outside of $P_C$, we can use the result of the tree-searches to obtain a disk in $S_C$ that does not contain $q$, as desired. Again, this takes $O(\log n)$ time.

To delete a a disk $s$ from $S$, we find the cell $C$ that contains the center of $s$, and we delete $s$ from $S_C$. For this, we must update the tree structures that represent the intersection $P_C$. The details for the algorithm and its analysis are given by Efrat, Itai, and Katz. Briefly, the idea is to find the parts of the boundary that are occluded by the deleted disk, and to promote these parts up in the tree. The main observation is that every piece of the boundary can rise only upwards into the tree, and that the tree has only $O(\log n)$ levels, resulting in a total deletion time of $O(n\log n)$. $◀$

Using Lemma 4, we get the following improved version of Lemma 3 for unit disk graphs:

The proof is identical to the proof of Lemma 3, but using the data structure from Lemma 4 instead of Lemma 2. $◀$

Lemma 3 and Lemma 5 allow us to compute the cardinality of the sought $(1-\epsilon )$-approximate clique. Computing the corresponding subset of nodes is a bit more involved.

We provide the analysis for general disk graphs. The case of unit disk graphs is handled in the same way, using the faster data structure from Lemma 4.

Our strategy is to find a $(1+\epsilon )$-approximate vertex cover in $\overline{G}$ in $O((n/\epsilon ){\log }^{4}n)$ expected time, and we return its complement. By Proposition 1, this constitutes a $(1-\epsilon )$-approximate clique in $G$.

To obtain the approximate vertex-cover, we proceed as follows: recall that $X$ and $Y$ are the nodes on the left and on the right side of $G$. As in the proof of Lemma 3, we start with executing $\mathrm{\ell }=1/\epsilon$ phases of the Hopcroft-Karp algorithm. This takes expected time $O((n/\epsilon ){\log }^{4}n)$, and it produces a $(1-\epsilon )$-approximate maximum matching $M$. Next, we again run a breadth-first search to partition the nodes into layers. Now, however, we do not need to terminate the search after reaching the first layer containing an unmatched node of $Y$. The BFS gives us a partition of $X$ and $Y$ into subsets

corresponding to the consecutive layers, where $X_0,X_1,\mathrm{\cdots }\subseteq X$, and $Y_0,Y_1,\mathrm{\cdots }\subseteq Y$, as well as the remaining part $X^{\prime }=X\setminus {\bigcup }_i{X}_i$ and $Y^{\prime }=Y\setminus {\bigcup }_i{Y}_i$. Because we have run $\mathrm{\ell }$ phases of the Hopcroft-Karp algorithm, the layers $Y_0,Y_1,\mathrm{\dots },Y_{\mathrm{\ell }-1}$ do not contain any unmatched nodes, and hence ${\sum }_i|Y_i|=|M|$. We choose $i\in \{0,1,\mathrm{\dots },\mathrm{\ell }-1\}$ such that $|Y_i|\le |M|/\mathrm{\ell }$, and we define the following set $K$:

We will now establish that $K$ is a vertex cover for $G$, with size $(1+\epsilon )|M|$.

Consider an edge $e=xy$, where $x\in X$ and $y\in Y$. We need to show that $x\in K$ or $y\in K$. If $x\in K$, we are done. Thus, suppose that $x\notin K$, and we need to establish that $y\in K$. Since $x\notin K$, we have $x\in X_j$, for some $j\in \{0,1,\mathrm{\dots },i\}$. We consider two cases.

• $\mathrm{■}$

  Case 1: $e\mathbf{\in }M$: Because $x\in X_j$ is matched, we must have $j\ge 1$ and we have reached $x$ through an augmenting path ending with $e$. Thus, $y\in Y_{j-1}\subseteq K$.

• $\mathrm{■}$

  Case 2: $e\mathbf{\notin }M$: We can extend an augmenting path ending at $x$ with $e$ to obtain an augmenting path ending at $y$. Thus, $y\in Y_0\cup Y_1\cup \mathrm{\cdots }\cup Y_j\subseteq K$.

This shows that $K$ is a vertex cover.

We now move to analyzing the size of $K$. First, we argue that every node in $K$ is an endpoint of an edge in $M$. We prove this separately for the nodes in $X$ and in $Y$.

• $\mathrm{■}$

  Case 1: $x\mathbf{\notin }{X}_{\mathrm{𝟎}}\mathbf{\cup }{X}_{\mathrm{𝟏}}\mathbf{\cup }\mathbf{\cdots }\mathbf{\cup }{X}_i$: Such a node is matched, as it does not belong to $X_0$.

• $\mathrm{■}$

  Case 2: $y\mathbf{\in }{Y}_{\mathrm{𝟎}}\mathbf{\cup }{Y}_{\mathrm{𝟏}}\mathbf{\cup }\mathbf{\cdots }\mathbf{\cup }{Y}_i$: None of the sets $Y_0,Y_1,\mathrm{\dots },Y_{\mathrm{\ell }-1}$ contain an unmatched node, so such a node is matched.

Next, we consider the situation when an edge $xy\in M$ has two endpoints in $K$. Then, $x\notin X_0\cup X_1\cup \mathrm{\cdots }\cup X_i$ and $y\in Y_j$, for some $j\in \{0,1,\mathrm{\dots },i\}$. Now, it is possible to extend an augmenting path ending at $y$ by the matched edge $yx$ to obtain an augmenting path ending at $x$, so $x\in X_{j+1}$. Thus, $j=i$, and the considered matched edge connects $x\in X_{i+1}$ and $y\in Y_i$. By the choice of $i$, there are at most $|M|/\mathrm{\ell }$ such edges. This allows us to bound the number of nodes in $K$ by $|M|+|M|/\mathrm{\ell }=(1+\epsilon )|M|$.

Together, this shows that $K$ is an ($1+\epsilon )$-approximate vertex cover.

As in Lemma 3, running the breadth-first search takes $O(n{\log }^{4}n)$ expected time. After that, constructing $K$ can easily be done in additional linear time. $◀$

Let $p\in X$, and let $C\in 𝒞$ be the cell that contains $p$. Since every other disk center in $X$ has distance at most $2$ from $p$, it follows that $X$ must lie in $C^{+}$. Hence, we have $X\subseteq P_C$, as claimed. $\mathrm{⊲}$

Furthermore, a simple packing argument shows that every extension of a nonempty cell $C$ contains a clique of size linear in $|P_C|$.

The extension $C^{+}$ consists of $25$ grid cells with pairwise disjoint interiors, where each cell has diameter $2$. Hence, for any such cell $C^{\prime }\subset C^{+}$, the set $P_C\cap C^{\prime }$ must form a clique, since all these disks are pairwise intersecting. Recall that if a disk center lies on the boundary of a cell, it is assigned arbitrarily to one of the adjacent cells so that each point in $P_C$ is uniquely assigned to exactly one cell in $C^{+}$. Since there are $25$ such cells, there exists at least one cell in $C^{+}$ containing at least $|P_C|/25$ disk centers. It follows that $P_C$ contains a clique of size at least $|P_C|/25$, and therefore the maximum clique $X^{\ast }$ is at least of size $|P_C|/25$. $\mathrm{⊲}$

Now, let $C\in 𝒞$. We select two disk centers $p_1,p_2\in P_C$ independently and uniformly at random. Let $s$ be the line segment of length 2 that starts at $p_2$ and is directed toward $p_1$, so that $p_1$ lies either on $s$ or on its extension along the same direction. Let $x$ be the other endpoint of $s$, and let $L=D(x,2)\cap D(p_2,2)$ be the lens with axis $s$; see Figure 3 for an illustration. We determine the set $P_L=L\cap P_C$ of disks centers that lie in $L$, and we find an $(1-\epsilon /2)$-approximate maximum clique $X$ for $P_L$, using Theorem 6 from Section 2. The lens $L$ is partitioned by segment $s$ into two parts: the part to the left of $s$ and the part to the right of $s$. Without loss of generality, the points on $s$ belong to the left part. Any two points in the same part have a distance at most $2$. Thus, all disks with centers in the same part form a clique. In other words, the intersection graph of unit disks centered at $P_L$ is a co-bipartite graph. Consequently, as we show below, the clique $X$ is an $(1-\epsilon )$-approximate clique for $P_C$ with probability $\mathrm{\Omega }(\epsilon )$.

First, by Claim 8, we know that $|X^{\ast }|\ge |P_C|/25$. Thus, with probability at least $1/25$, we have that $p_1\in X^{\ast }$. Now, we sort the centers in $X^{\ast }$ in the increasing order of their distances from $p_1$, and let ${\pi }^{\ast }$ be the resulting ordering of $X^{\ast }$. Let $X^{+}$ be those centers from $X^{\ast }$ that have rank at least $(1-\epsilon /2)|X^{\ast }|$ in ${\pi }^{\ast }$, and let $X^{-}=X^{\ast }\setminus X^{+}$ (in other words, $X^{+}$ contains the $(\epsilon /2)|X^{\ast }|$ points in $X^{\ast }$ that are furthest away from $p_1$, and $X^{-}$ contains the $(1-\epsilon /2)|X^{\ast }|$ remaining points). Since

it follows that with probability at least $\epsilon /50$, the second sampled center $p_2$ lies in $X^{+}$.

Since $p_1$ and $p_2$ are sampled independently, we can conclude that with probability at least $\epsilon /1250$, it holds that $p_1\in X^{\ast }$ and that $p_2\in X^{+}$. Henceforth, we assume that indeed this is the case. Since $p_1$ and $p_2$ both lie in $X^{\ast }$, the distance $\mathrm{dist}(p_1,p_2)$ between $p_1$ and $p_2$ is at most $2$. From this, it follows that $D(p_1,\mathrm{dist}(p_1,p_2))\subseteq D(x,2)$ since by the choice of $x$, we have $\mathrm{dist}(x,p_1)+\mathrm{dist}(p_1,p_2)=2$, which implies that for every point $y\in D(p_1,\mathrm{dist}(p_1,p_2))$, the triangle inequality gives

as needed.

As $p_2\in X^{+}$, and since all the points in $X^{+}$ are at least as far away from $p_1$ as all the points in $X^{-}$, it follows that $X^{-}\subseteq D(p_1,\mathrm{dist}(p_1,p_2))\subseteq D(x,2)$. Moreover, since $X^{-}\subseteq X^{\ast }$ and $p_2\in X^{\ast }$, we also have $X^{-}\subseteq D(p_2,2)$. This implies that $X^{-}$ lies in the lens $L$ defined by $D(x,2)$ and $D(p_2,2)$. In particular, $P_L$ contains a clique of size $(1-\epsilon /2)|X^{\ast }|$, and since $X$ is a $(1-\epsilon /2)$-approximate clique for $P_L$, we have

This means that $X$ is an $(1-\epsilon )$-approximate maximum clique for $P_C$, as claimed. $◀$

We can now present the complete algorithm.

Let $C\in 𝒞$ be a nonempty grid cell. By Lemma 9, and using Theorem 6, one can compute, in $O((|P_C|/\epsilon )\log n)$ expected time, a clique for $P_C$ that, with probability $\mathrm{\Omega }(\epsilon )$, is an $(1-\epsilon )$-approximate clique for $P_C$. By repeating this process $O(1/\epsilon )$ times, and selecting the largest of the resulting cliques, one can compute in $O((|P_C|/{\epsilon }^2)\log n)$ expected time, a clique that is an $(1-\epsilon )$-approximate clique for $P_C$ with probability at least $1/2$.

We do this for $P_C$ of every nonempty cell $C\in 𝒞$, and we return the largest of the resulting cliques. (Note that all non-empty cells can be easily found in $O(n\log n)$ time [36]). By Claim 7, there exists a cell $C\in 𝒞$ such that $P_C$ contains a max-clique of $P$, so with probability at least $1/2$, this gives an $(1-\epsilon )$-approximate maximum clique for $P$. The total running time is proportional to

since every point from $P$ appears in $O(1)$ sets $P_C$. The theorem follows. $◀$

By a standard argument, for any value $\delta >0$, we can raise the probability of success to $1-\delta$ by independently repeating the algorithm of Theorem 10 $O(\log (1/\delta ))$ times and by taking the best solution. This increases the running time by an $O(\log (1/\delta ))$ factor.