Simple Randomized Algorithms for Closest Pair Problems

Mordecai Golin^{∗} Rajeev Raman^{†} Christian Schwarz^{†}

Michiel Smid^{†}

December 30, 1992

Abstract

We present a conceptually simple, randomized incremental algorithm for finding the closest pair in a set of n points in D-dimensional space, where D ≥ 2 is a fixed constant. Using dynamic perfect hashing, the algorithm runs in O(n) expected time. In addition to being quick on the average, this algorithm is reliable: we show that it runs in O(n log n/ log log n) time with high probability.

1 Introduction

The closest pair problem is to find a closest pair in a given set of points. The problem has a long history in computational geometry and has been extensively studied. It is well known, for example, that finding the closest pair in a set of n points requires Ω(n log n) time in the algebraic decision tree model of computation [1] and that there are optimal algorithms which match this lower bound. It is also well known that if the model of computation is changed then the lower bound no longer holds. This was first shown by Rabin [9] who described an algorithm that combines the use of the floor function with randomization to achieve an O(n) expected running time. The expectation is taken over choices made by the algorithm and not over possible inputs. Recently, Khuller and Matias [6] have described a radically different algorithm that also uses the floor function and randomization to achieve an O(n) expected running time.

In this paper, we present yet another algorithm that combines the use of the floor function and randomization to solve the problem in O(n) expected time. Our algorithm is conceptually simpler than the ones in [9] and [6]. It also differs from them in that it is a randomized incremental algorithm. The other two algorithms are inherently static. The algorithm that we present here assumes that the input points are fed to it in a sequence which is a random permutation p_{1}, p_{2}, …, p_{n} of the n points. The i-th stage of the algorithm will find the closest pair of the set S_{i} := {p_{1}, p_{2}, …, p_{i}} in O(1) expected time leading to an overall O(n) expected running time. To prove this running time, we use the technique of backwards analysis, due to Seidel. See [12, 13].

∗INRIA-Rocquencourt, 78153 Le Chesnay Cedex, France, and Hongkong UST. This author was supported by the ESPRIT Basic Research Actions Program, under contract No. 7141 (project ALCOM II)

and by NSF grant CCR-8918152. Part of this work was done while this author was visiting the Max-Planck-Institut für Informatik.

†Max-Planck-Institut für Informatik, W-6600 Saarbrücken, Germany. These authors were supported by the ESPRIT Basic Research Actions Program, under contract No. 7141 (project ALCOM II).

In Section 2, we present the algorithm and its analysis. We give an O(n log n) expected time tree based algorithm and show how to use dynamic perfect hashing [3] to improve its expected running time to O(n).

In Section 3, we consider the problem of returning not only the closest pair, but all the k closest pairs, where k is an integer between 1 and \binom{n}{2}. For this problem, there are deterministic algorithms with running time O(n log n + k), which is optimal in the algebraic decision tree model of computation. (See [7, 10].) Combined with randomization and the floor function, we get a simple algorithm with expected running time O(k n), which is better than the algorithms in [7, 10] if k = o(log n).

In Section 4, we show that our closest pair algorithms are reliable: For example, the linear expected time algorithm actually runs in O(n log n/ log log n) time with high probability.

All algorithms of this paper assume that the floor function can be computed at unit cost. In Section 5, however, we give a variant of the closest pair algorithm presented in Section 2 that has O(n log n) expected running time, even with high probability, without using this non-algebraic function. This algorithm fits in the algebraic decision tree model of computation, extended with the power of randomization.

2 The closest pair algorithm

To keep our exposition simple, we give the algorithm for the two-dimensional case and the euclidean metric. The extension to arbitrary, but fixed, dimension D ≥ 2 and arbitrary L_{t}-metric, 1 ≤ t ≤ ∞, is straightforward. Let d(p, q) denote the distance between the points p = (p^{(1)}, p^{(2)}) and q = (q^{(1)}, q^{(2)}), i.e.,

d(p, q) = \sqrt{(p^{(1)} − q^{(1)})^{2} + (p^{(2)} − q^{(2)})^{2}}.

Let S = {p_{1}, p_{2}, …, p_{n}} be a set of points. The closest pair distance in S is

δ(S) := min{d(p, q) : p, q ∈ S, p ≠ q}.

The closest pair problem is to find a pair of points p, q ∈ S such that d(p, q) = δ(S).

Our algorithm will be based upon the following simple observation. Let S_{i} := {p_{1}, p_{2}, …, p_{i}} be the set containing the first i points of S. Then δ(S_{i+1}) < δ(S_{i}) if and only if there is some point p ∈ S_{i} such that d(p, p_{i+1}) < δ(S_{i}).

Suppose a square grid with mesh size δ(S_{i}) is laid over the plane^{1} and each point of S_{i} is stored in the grid box in which it appears. Let b be the grid box in which the new point p_{i+1} is located. Then every point in S_{i} that is within distance δ(S_{i}) of p_{i+1} must be located in one of the 9 grid boxes that are adjacent to b. (We consider the box b as being adjacent to itself.) We call these 9 boxes the neighbors of b.

We note that each grid box can only contain at most four points from S_{i}. This is because if a grid box contained more than four points then some pair of them would be less than δ(S_{i}) apart, contradicting the definition of δ(S_{i}).

1To exactly specify the grid we will always assume that (0, 0) is one of its lattice points.

The above observations lead to the following generic algorithm for finding the closest pair in S. The points will be fed to the algorithm in random order p_{1}, p_{2}, …, p_{n}, i.e., each of the n! possible orders is equally likely. The algorithm starts by calculating δ(S_{2}) = d(p_{1}, p_{2}) and inserting the points of S_{2} into the grid with mesh size δ(S_{2}). It then proceeds incrementally, always keeping set S_{i} stored in a grid with mesh size δ(S_{i}). When fed point p_{i+1} it finds the at most 36 points in the 9 grid boxes neighboring the box in which p_{i+1} is located and computes d_{i+1}, the minimum distance between p_{i+1} and these at most 36 points. If there are no points in these boxes then d_{i+1} = ∞. From the discussion above we know that δ(S_{i+1}) = min (d_{i+1}, δ(S_{i})) .

If d_{i+1} ≥ δ(S_{i}) then δ(S_{i+1}) = δ(S_{i}) and the algorithm inserts p_{i+1} into the current grid. Otherwise, δ(S_{i+1}) = d_{i+1} < δ(S_{i}) and the algorithm discards the old grid, creates a new one with mesh size δ(S_{i+1}), and inserts the points of S_{i+1} into this grid.

The algorithm thus calculates δ(S_{i}) for i = 2, 3, …, n, in this order. Then it outputs the value δ(S_{n}) = δ(S). An example of our algorithm is shown in Figure 1.

Remark: Instead of one closest pair, we can also incrementally maintain a list of all point pairs in S_{i} with distance δ(S_{i}), for 2 ≤ i ≤ n, without additional cost. Consider the i-th stage of the algorithm, when we add p_{i+1} to S_{i}. Since the algorithm checks all points within distance δ(S_{i}) of p_{i+1} anyway, we find all points q in S_{i} such that d(q, p_{i+1}) = δ(S_{i}). The pairs (q, p_{i+1}) with δ(q, p_{i+1}) = δ(S_{i+1}) — there can be only O(1) of them — are exactly the new closest pairs after the i-th stage of the algorithm. If δ(S_{i+1}) = δ(S_{i}), we add these pairs to the current closest pair list. If δ(S_{i+1}) < δ(S_{i}), we discard the old list and make a new list that contains just these pairs.

To actually implement the above algorithm we will need the following: let P be a point set, d a positive real number, p a point, G a grid, and b the name of a box in a grid. We define the following operations.

• Build(P, d) : Return a grid G with mesh size d that contains the points in P.

• Insert(G, p) : Insert point p into grid G.

• Report(G, b) : Return all points in grid box b.

Pseudocode for the closest pair algorithm using these operations is presented in Figure 2.

How do we actually implement these grid operations? One easy way is to use balanced binary search trees: Consider a grid G with mesh size d, let p = (p^{(1)}, p^{(2)}) be a point in the plane, and denote the box containing p in G by b_{p}. The integer pair (⌊p^{(1)}/d⌋, ⌊p^{(2)}/d⌋) is called the index of b_{p}. The point set is stored as follows: We determine the indices of the non-empty boxes and store them in lexicographical order in a balanced binary search tree. Moreover, with each box in this tree, we store a list of all points that are contained in this box. To insert a point into the grid, we use the floor function to compute the index of the grid box that contains this point. Then we search in the tree for this box. If it is storet�����������