Closest pair of points
Also known as: Closest pair of points · closest pair problem · two nearest points · Shamos–Hoey algorithm · finding duplicate coordinates
Finding the two nearest of tens of thousands of points without trying billions of pairs: cut in half, solve the halves, and check a narrow strip at the cut, where each point has at most seven candidates. In real catalogues, earthquakes, airports, cities, the closest pair almost always turns out to be a duplicate, and that is why people look for it.
The USGS catalogue has 19,073 earthquakes of magnitude 5 and above for 2015–2025. Which two epicentres are closest to each other? Comparing all pairs means 181,880,128 distances. The algorithm below makes do with thirteen and a half thousand. Choose a set of points, earthquakes, airports or cities of the world, and a search method; the pair found is shown on the map and under the lens.
The first answer is disconcerting: 0 kilometres. Two earthquakes, a magnitude 6.1 on 22 April 2024 and a magnitude 5.0 on 2 June 2024, are recorded at the same point to a thousandth of a degree. The catalogue has three such coincidences. One of them is not two earthquakes but one recorded twice: in the USGS database the magnitude 5.3 event of 27 February 2021 at 18:59:25 has two different identifiers with the same time, place, depth and magnitude. For airports the closest pair is a record literally called “(Duplicate)YEG” and the real Edmonton airport. For cities, a city and its own district that the gazetteer put at one point.
This is what the closest pair looks like in almost any real catalogue: a duplicate. That is why people look for it: checking data, the first thing to look at is which records coincide or nearly coincide. Turn on “distinct points only”, and the answer becomes meaningful: two earthquakes on 24 May 2017, 68 metres apart, the second, stronger one 29 minutes after the first; two Krasnoyarsk airports 3.7 kilometres apart.
The idea
Sort the points by x and cut them with a vertical line into two equal halves. Find the closest pair in each half the same way; let the smaller of the two distances be d. It remains to check pairs that cross the cut, and such a pair can be closer than d only if both points lie in a strip 2d wide around the cut. Go through the points of the strip from bottom to top, and compare each only with the following ones while the difference in height is less than d: there are at most seven of them.
This is the classic algorithm of Michael Shamos and Dan Hoey. Two details make it fast. The first is the strip: the distance to the cut cuts off everything farther than d without computing a single distance. The second is sorting the strip by y: if it is sorted anew in every call, the combining work becomes n log n, and the whole thing n log² n. Instead, the recursion returns the points of its piece already ordered by y, and combining merges two ordered lists, as merge sort does, in linear time.
Step by step
Below are 174 towns of 10 thousand people and more around Almaty, between 40° and 46° north and 68° and 81° east, laid on a flat map. Start the search: the recursion cuts the points in half until a piece has two or three left, solves the pieces, and on the way back checks the strip at each cut. Thin lines are distances computed, the thick one the best pair so far.
The closest turn out to be Nizhnyaya Ala-Archa and Mayevka, 1.45 kilometres apart, and the search computed 168 distances out of a possible 15,051 for this. Watch the strips: while d is large, many towns fall into a strip, but as d decreases the strips narrow to a handful of points. Switch to “in the strip: every pair”, and you can see why the second detail is needed: without stopping by height, the strips at the upper cuts are compared each with each.
Why the strip has few candidates
Let d be the smallest of the distances inside the left half and inside the right half. If the points of the strip are ordered by y, then for each point p all points of the strip closer than d to it and above it are among the seven following p.
A point q above p and closer than d lies in a rectangle 2d wide and d high: horizontally within the strip, vertically between yp and yp + d. The cut divides the rectangle into two d × d squares; points of the left half lie in the left square, of the right half in the right one (points on the cut itself have the same x and lie in both, but each is counted in the square of its own half).
How many points of one half can be in its square? Any two of them are at least d apart, otherwise the half would contain a pair closer than d. Divide the square into four small squares with side d/2. Two points in one small square are at most its diagonal d/√2 < d apart, so each small square holds at most one point, and the square at most four. The rectangle holds at most eight points in all, counting p, and in y order they all come right after p (before any point above yp + d). The other seven are all there is, and nothing more needs checking.
The code does not count to seven: it goes up the strip while the difference in height is less than d, and stops. The lemma guarantees that there are few such steps from each point, and the height check cuts points off even earlier. The bound is generous: with more careful counting seven can be reduced, but that does not affect the asymptotics.
Code
The classic algorithm on the plane. Coordinates are whole metres, so squared distances are compared, without roots and without rounding errors. The recursion receives point numbers in increasing order of x and returns them in increasing order of y; ties in a coordinate are broken by the point number, so that both versions of the program take the same path. Besides the towns there are two inputs that look dangerous: all points on one vertical line, where everything falls into the strip, and a grid, where there are ties at every step.
one vertical line: 20000 points, closest 6, 15429 distances computed (0.77 per point) 100 x 100 grid: 10000 points, closest 10, 7712 distances computed (0.77 per point)
On the vertical line all points are in the strip, but the strip is checked only up to the first neighbour above d, and there is less than one distance per point.
The first scene is the same recursion in three dimensions. A point on the sphere is written as a unit vector (cos φ cos λ, cos φ sin λ, sin φ), and the straight distance through the Earth, the chord, grows together with the distance along the surface, so the nearest by chord are the nearest on the map. A cut by x gives a slab rather than a strip; the points of the slab are gone through by y, and a distance is computed only if the z coordinate is close too. The lemma carries over: in the slab each point has a bounded number of candidates, only the constant is larger.
What it costs
Sorting by x once is O(n log n). The recursion: two halves, and combining (the merge by y, selecting the strip and walking it) all in O(n): T(n) = 2T(n/2) + O(n), the second case of the master theorem, O(n log n). Memory is O(n) for the lists and log n recursion frames.
Nothing faster is possible in the worst case if the algorithm only compares and adds coordinates: the question “are all the numbers distinct” reduces to the closest pair: put the numbers on a line, and they are distinct only if the closest pair is at a non-zero distance. And for that question the algebraic decision tree model has a lower bound of Ω(n log n).
A measurement on random points, divide and conquer against trying all pairs and against the grid described below:
Brute force grows as n²: in Python two thousand points already take more than a second, and 32 thousand points would take about 250 times longer, some five minutes. Divide and conquer on the same 32 thousand takes about six tenths of a second.
A grid: O(n) on average
There is another way, without recursion. Suppose a distance d is known that is no less than the answer. Put the points into grid squares with side slightly more than d: then a point closer than d to a given one can be only in its square or in the eight neighbouring ones. If the squares are kept in a hash table, checking a point is a constant number of operations.
But d is not known in advance. Michael Rabin found a way out in 1976, and the variant in the listing, with a random order of points and rebuilding the grid, is covered in textbooks such as “Algorithm Design” by Jon Kleinberg and Éva Tardos. Take the points in random order and keep the grid for the current best d. If a new point gives a closer pair, rebuild the grid for the new d. A rebuild costs as many as the points already inserted, and it seems this might add up to n². But on average it does not.
If points are inserted in random order, rebuilding the grid takes O(n) insertions in total on average.
Look backwards. Among the first k points there is a closest pair. A rebuild at step k happens only if the k-th point is one of the two ends of this pair. With a random order any of the first k points is equally likely to be the last, so the probability of a rebuild at step k is at most 2/k (if there are several closest pairs, it is only smaller). A rebuild costs k insertions, and on average step k adds at most k · 2/k = 2. Over all steps, at most 2n.
The listing confirms the bound: on a hundred thousand random points the grid was rebuilt 20 times and moved 185,900 points in total, fewer than two per point. But “linear on average” does not mean “faster”: in Python the grid beats the recursion about twofold, while in Go, where a hash table with pair keys costs more than recursion over slices, divide and conquer is faster than the grid. The grid, on the other hand, generalises more easily: the same placement into squares answers the question “which pairs are closer than r”, the subject of the second problem.
Where it is used
- Checking data. Identical and nearly identical coordinates are the first sign of duplicates: one earthquake under two numbers, one airport as two records, a city and its district at one point. The closest pair answers “are there duplicates at all”, a grid placement “which ones”.
- Proximity control. Whether moving objects have come closer than allowed: at each moment this is a closest pair question, and that is how it is solved when there are many objects.
- Clustering. Single-linkage hierarchical clustering starts by merging the two nearest points, then the two nearest groups.
- Many questions about the same points. If you need not one pair but the nearest neighbour for many points or queries, you build a k-d tree or a Delaunay triangulation: the closest pair is always joined by an edge of the triangulation.
History
The closest pair problem is one of the first with which computational geometry began in the 1970s. Michael Shamos and Dan Hoey described an O(n log n) algorithm on the plane in 1975. In 1976 Jon Bentley and Shamos generalised it to any fixed dimension, with the same O(n log n), and in the same year Michael Rabin showed that with randomness and a grid it can be done in O(n) on average; this work became one of the first examples of randomness speeding up algorithms. In 1979 Steven Fortune and John Hopcroft obtained O(n log log n) without randomness, but with a floor operation that is absent from the model where the lower bound holds. In 1995 Samir Khuller and Yossi Matias proposed a very simple randomised algorithm: a sieve that first estimates the distance roughly and then refines it with a grid.
Problems
In the first problem you need the closest pair itself on the real cities of the world, and a test on which the popular “compare each point with a few neighbours by x” fails. In the second, count all pairs closer than a given distance: here a grid is handier.
The two closest towns Pairs closer than rData: Earthquakes: USGS catalogue, public domain Airports: OurAirports, public domain Places: GeoNames, CC BY 4.0