Closest Pair of Points, in depth: the O(n log n) divide-and-conquer algorithm, why the strip check is linear, and how to implement it correctly

The brute-force baseline

Start with the baseline, because you will need it twice: as the base case of the recursion and as the oracle that tests your fast version. Compare squared distances rather than distances. Squaring is monotonic for non-negative numbers, so the closest pair under squared distance is the closest pair under distance, and with integer coordinates the squared distance is an exact integer with no rounding at all. Take a square root once, at the end, only if you need the actual distance.

from itertools import combinations
import math

def d2(p, q):
    return (p[0] - q[0]) ** 2 + (p[1] - q[1]) ** 2

def brute(points):
    """O(n^2) reference: (squared distance, (p, q))."""
    best = (math.inf, None)
    for p, q in combinations(points, 2):
        dd = d2(p, q)
        if dd < best[0]:
            best = (dd, (p, q))
    return best

For n = 10,000 the brute force makes about 50 million comparisons, which is slow but tolerable in pure Python. For n = 1,000,000 it is about 500 billion, which is not practical in any language. The goal is to touch each point only a logarithmic number of times.

Divide and conquer: the strip argument

Sort the points by x once. Split the sorted list at its middle index into a left half and a right half, separated by a vertical line at the middle point's x coordinate. Recursively find the closest pair inside each half; call the smaller of the two distances delta. The overall answer is either that pair or a pair with one point on each side of the line, so the whole difficulty is the cross pairs.

A cross pair closer than delta must have both points within delta of the dividing line horizontally, otherwise their horizontal separation alone already exceeds delta. So only points in a vertical strip of width 2 delta around the line can matter. That alone does not save us, because in the worst case every point is in the strip. The second observation does.

Sort the strip points by y. Take any strip point p and look only at strip points above it. A partner closer than delta lies inside a rectangle delta tall and 2 delta wide sitting on p. Cut that rectangle into eight squares of side delta over 2. Each square lies entirely on one side of the line, and two points in the same square would be at most delta over the square root of 2 apart, which is less than delta, contradicting the fact that each side has no pair closer than delta. So each square holds at most one point, the rectangle holds at most eight points including p, and p needs to be compared with at most the next 7 points in y order. Tighter constants are known, but 7 is a safe bound and is what most textbook implementations use.

In practice you also stop as soon as the next point is delta or more above p in y, which usually ends the inner loop after one or two comparisons. The strip step therefore costs O(n) comparisons, as long as you can get the strip in y order without sorting it from scratch.

A tested implementation

The implementation below returns the y-sorted list of its points from every recursive call, so the parent can merge two sorted lists in linear time, exactly as mergesort does. That merge is what makes the total O(n log n). Re-sorting the strip by y at every level instead gives O(n log squared n), which is still fast but not optimal, and is a common mistake in published code.

def closest_pair(points):
    """Return (squared distance, (p, q)) for len(points) >= 2."""
    px = sorted(points)                       # by x, ties by y

    def solve(lo, hi):                        # px[lo:hi] -> (d2, pair, by_y)
        if hi - lo <= 3:
            best = brute(px[lo:hi])
            return best[0], best[1], sorted(px[lo:hi], key=lambda p: p[1])
        mid = (lo + hi) // 2                  # split by index, not by x value
        mid_x = px[mid][0]
        dl, pl, yl = solve(lo, mid)
        dr, pr, yr = solve(mid, hi)
        best, pair = (dl, pl) if dl <= dr else (dr, pr)

        ys, i, j = [], 0, 0                   # merge step, as in mergesort
        while i < len(yl) and j < len(yr):
            if yl[i][1] <= yr[j][1]:
                ys.append(yl[i]); i += 1
            else:
                ys.append(yr[j]); j += 1
        ys += yl[i:] + yr[j:]

        strip = [p for p in ys if (p[0] - mid_x) ** 2 < best]
        for k, p in enumerate(strip):
            for q in strip[k + 1:k + 8]:      # at most 7 followers
                if (q[1] - p[1]) ** 2 >= best:
                    break                     # everything later is higher still
                dd = d2(p, q)
                if dd < best:
                    best, pair = dd, (p, q)
        return best, pair, ys

    best, pair, _ = solve(0, len(px))
    return best, pair

Three details are deliberate. Splitting by index rather than by an x threshold guarantees both halves shrink even when many points share an x coordinate; a threshold split can put every point on one side and recurse forever. Comparing squared values keeps integer inputs exact. And the strip uses strict less-than, so a strip point exactly delta away is skipped, which is correct because it cannot improve on delta.

Test it the way you would test any clever algorithm, against the brute force on many random small inputs with deliberately small coordinate ranges, so duplicates and collinear points occur often:

import random
for _ in range(3000):
    n = random.randint(2, 60)
    pts = [(random.randint(-20, 20), random.randint(-20, 20)) for _ in range(n)]
    assert closest_pair(pts)[0] == brute(pts)[0], pts

This exact harness passed 3,000 random cases against the implementation above when this article was written. Compare distances, not pairs, because ties can legitimately return different pairs.

Worked example: ten points, step by step

Run it on ten points: (2,3), (5,9), (8,1), (10,20), (12,10), (14,13), (20,18), (25,24), (31,33), (40,50). Sorted by x, the middle index is 5, so the left half is (2,3), (5,9), (8,1), (10,20), (12,10) and the right half starts at (14,13), putting the dividing line at x = 14.

The left half recurses again. Its own left part, (2,3) and (5,9), has squared distance 9 + 36 = 45. Its right part, (8,1), (10,20), (12,10), is a base case: the best there is (8,1) to (12,10) at 16 + 81 = 97. So the left half starts its strip check with delta squared 45, and the strip around x = 8 finds (2,3) to (8,1) at 36 + 4 = 40. The left half's answer is 40. The right half, by the same process, returns (14,13) to (20,18) at 36 + 25 = 61.

At the top level delta squared is min(40, 61) = 40, so delta is about 6.32 and the strip keeps points with x strictly between about 7.68 and 20.32: (8,1), (12,10), (14,13), (20,18) and (10,20). In y order they are (8,1), (12,10), (14,13), (20,18), (10,20). The scan goes:

1. From (8,1): the next point (12,10) is 9 higher, and 81 is at least 40, so stop at once.
2. From (12,10): (14,13) is 3 higher, 9 is below 40, so compute 4 + 9 = 13. That beats 40, so delta squared becomes 13. The next point (20,18) is 8 higher; 64 is at least 13, so stop.
3. From (14,13): (20,18) is 5 higher, 25 is at least 13, stop.
4. From (20,18): (10,20) is 2 higher, so compute 100 + 4 = 104, which is no improvement. The strip is exhausted.

The answer is (12,10) and (14,13), squared distance 13, distance about 3.606, a pair that neither recursive call could see because it straddles the line. The scan made two distance computations in total, which is the typical behavior: the y break does most of the pruning, and the 7-point bound only guarantees the worst case.

Cost, constants and alternatives

The recurrence is T(n) = 2T(n/2) + O(n): two half-size subproblems plus a linear merge and a linear strip scan, which solves to O(n log n), plus the O(n log n) initial sort. Memory is O(n) for the sorted copies and the per-level y lists. Recursion depth is about log base 2 of n, roughly 20 for a million points, so Python's default recursion limit of 1000 is not a concern for this algorithm.

The constant factors are worth knowing. The Python version allocates a new list at every merge, so for very large inputs a NumPy or compiled implementation will be one to two orders of magnitude faster. In a compiled language, avoid allocating per call: preallocate a scratch buffer for the merge, as an efficient mergesort does, and use 64-bit integers for squared distances. With coordinates up to 10^9 in magnitude, a coordinate difference can reach 2 times 10^9 and a sum of two squares can reach 8 times 10^18, which fits in a signed 64-bit integer with little room to spare; larger coordinates need 128-bit arithmetic or floats.

Failure modes

Most broken implementations fail in one of these ways:

• Duplicate points. Two identical points give distance 0. The algorithm handles it, but code that divides by the distance or uses it as a grid cell size will fail; check for zero explicitly if you do either.
• Splitting by x value. Choosing the side by comparing with a threshold sends all points sharing that x to one side, so a column of points recurses forever. Split by index.
• Floating-point comparisons. With float coordinates, distances computed in different orders can differ in the last bit, so a rare test failure may be a rounding tie rather than a logic bug. Compare against the oracle with a small tolerance, or use integers or exact rationals when correctness must be exact.
• Re-sorting the strip. Correct, but O(n log squared n). It matters only at large n, yet it is often the reason a 'fast' implementation underperforms.
• Wrong metric. Latitude and longitude are not planar coordinates; at high latitudes a degree of longitude is much shorter than a degree of latitude. Project to a local planar system first, or use a spatial index that supports great-circle distance.
• Mutating shared input. Sorting the caller's list in place changes their data. Sort a copy, as the reference does.

Trade-offs: when to use something else

Divide and conquer is the right default when you get one batch of points, need one answer and want a guaranteed bound. If points stream in sorted by x, a sweep line keeps an ordered set of the points within delta of the sweep position and checks each new point against its y neighbourhood; it has the same bound and processes points incrementally. If you have tens of millions of points and a good hash table, the randomized grid approach, whose idea goes back to Rabin's 1976 work on probabilistic algorithms, achieves expected linear time by bucketing points into cells of the current best distance and checking neighbouring cells only.

If the closest pair is one of many proximity questions you will ask, build an index instead. A k-d tree answers nearest-neighbour, radius and range queries after one O(n log n) build, and running a nearest-neighbour query from each point also yields the closest pair. The same sorted-order thinking appears in Graham scan, the other classic sort-then-sweep geometry algorithm.

Finally, higher dimensions: the divide-and-conquer idea generalizes, but the strip bound grows quickly with dimension, and in high dimensions tree indexes degrade too. For hundreds of dimensions, such as embedding vectors, approximate nearest-neighbour methods are the practical choice.

What to do next

1. Implement the brute force and the divide-and-conquer version from this page and run the random test harness until it passes.
2. Reproduce the ten-point trace by printing delta squared and the strip at each level.
3. Replace the merge with a re-sort, time both on a million random points, and observe the log factor.
4. Port the algorithm to a compiled language with a preallocated merge buffer and 64-bit squared distances.
5. Try inputs that break naive code: all points on one vertical line, many duplicates, and float coordinates with ties.
6. Implement the sweep-line or randomized grid variant and compare them on your real data distribution before choosing one.