vault

Closest Pair of Points

17 min read

Closest Pair of Points

Given a set P of n points in the Euclidean plane, the Closest Pair of Points problem asks for the two points p, q ∈ P that minimise the Euclidean distance ‖p − q‖₂. The naive O(n²) solution compares all pairs. The classical divide-and-conquer algorithm of Shamos and Hoey (Shamos-Hoey 1975) solves it in O(n log n) by sorting by x, recursing on each half, and then performing a clever strip merge along the dividing line — exploiting a tight geometric packing argument that says any point in the strip needs to be compared against at most a constant number (≤ 7) of subsequent strip points. The recurrence is T(n) = 2 T(n/2) + O(n) which by the Master Theorem resolves to Θ(n log n). The closest-pair problem is the canonical first example in any computational-geometry course (CLRS Ch. 33.4; de Berg et al. ch. 1) and a classical interview test of “do you understand the divide-and-conquer recurrence beyond Merge Sort and Quicksort?“

1. Intuition — Slice the Plane, Then Patch the Seam

If you have a thousand stars on a starfield and want to find the two closest, the obvious thing — measure every pair — is 1000 × 999 / 2 ≈ 500,000 measurements. Tedious but tractable for n = 1000; impossible for n = 10⁶.

Divide-and-conquer says: split the starfield down the middle by a vertical line. Find the closest pair in the left half; find the closest pair in the right half. Take the smaller distance, call it δ. The closest pair in the whole field is either this δ-pair or a cross-the-line pair where one point is left of the line and the other is right.

The crucial geometric observation: a cross-the-line pair that beats δ must have both endpoints within distance δ of the dividing line — otherwise the pair is already at least δ apart just from the horizontal gap. So we only need to look at points in a vertical strip of width 2δ centred on the dividing line, and among those, only at pairs whose y-coordinates are within δ of each other.

The remaining miracle is packing: in a δ × 2δ rectangle, you can fit at most a constant number of points that are themselves all ≥ δ apart from each other. This bounds the strip-merge work to O(n), which is what makes the recurrence T(n) = 2T(n/2) + O(n) (i.e., O(n log n)) instead of T(n) = 2T(n/2) + O(n²) (i.e., O(n² log n), which would be worse than the trivial brute force).

2. Tiny Worked Example

Take 8 points (chosen to spell out the algorithm clearly):

P = [ A=(2, 3), B=(3, 8), C=(5, 5), D=(6, 1),
      E=(8, 7), F=(9, 4), G=(11, 6), H=(12, 2) ]

Step 1: sort by x-coordinate. Already sorted: A, B, C, D, E, F, G, H.

Step 2: split into halves. L = {A, B, C, D} (x = 2..6), R = {E, F, G, H} (x = 8..12). Dividing line at x = 7 (midpoint).

Step 3: recurse on each half.

In L:

  • d(A, B) = √(1² + 5²) = √26 ≈ 5.10
  • d(A, C) = √(3² + 2²) = √13 ≈ 3.61
  • d(A, D) = √(4² + 2²) = √20 ≈ 4.47
  • d(B, C) = √(2² + 3²) = √13 ≈ 3.61
  • d(B, D) = √(3² + 7²) = √58 ≈ 7.62
  • d(C, D) = √(1² + 4²) = √17 ≈ 4.12

Closest in L: (A, C) or (B, C), both √13 ≈ 3.61. Call it δ_L = 3.61.

In R:

  • d(E, F) = √(1² + 3²) = √10 ≈ 3.16
  • d(E, G) = √(3² + 1²) = √10 ≈ 3.16
  • d(E, H) = √(4² + 5²) = √41 ≈ 6.40
  • d(F, G) = √(2² + 2²) = √8 ≈ 2.83
  • d(F, H) = √(3² + 2²) = √13 ≈ 3.61
  • d(G, H) = √(1² + 4²) = √17 ≈ 4.12

Closest in R: (F, G), distance √8 ≈ 2.83. Call it δ_R = 2.83.

Step 4: take δ = min(δ_L, δ_R) = 2.83. This is the best we know so far.

Step 5: build the strip. Keep only points with |x − 7| < 2.83. That’s 4.17 < x < 9.83:

  • D = (6, 1) ✓
  • E = (8, 7) ✓
  • F = (9, 4) ✓

(C = (5, 5) is at |5 − 7| = 2 < 2.83, also in the strip.)

Strip = {C=(5,5), D=(6,1), E=(8,7), F=(9,4)}.

Step 6: sort strip by y. D=(6,1) → F=(9,4) → C=(5,5) → E=(8,7) (y = 1, 4, 5, 7).

Step 7: for each strip point, compare with the next 7 points in y-order; if y-difference exceeds δ, stop early.

  • D=(6,1):
    • vs F=(9,4): |y diff| = 3 < δ; d = √(3² + 3²) = √18 ≈ 4.24 > 2.83, no improvement.
    • vs C=(5,5): |y diff| = 4 > δ. Stop.
  • F=(9,4):
    • vs C=(5,5): |y diff| = 1 < δ; d = √(4² + 1²) = √17 ≈ 4.12 > 2.83, no improvement.
    • vs E=(8,7): |y diff| = 3 > δ. Stop.
  • C=(5,5):
    • vs E=(8,7): |y diff| = 2 < δ; d = √(3² + 2²) = √13 ≈ 3.61 > 2.83, no improvement.

No cross-line pair beats 2.83. The closest pair is (F, G) with distance √8 ≈ 2.83.

Notice that we examined only 5 strip-pair distances, not the full 4 × 4 = 16 cross-product of strip points. That’s the packing argument paying off.

3. Pseudocode — The Classical Shamos-Hoey Recursion

closest_pair(P):
    if length(P) <= 3:
        return brute_force(P)                          # base case

    sort P by x-coordinate
    return closest_pair_recursive(P)

closest_pair_recursive(P):
    n := length(P)
    if n <= 3:
        return brute_force(P)

    mid     := n / 2
    mid_x   := P[mid].x
    L       := P[0 .. mid]
    R       := P[mid .. n]

    (p1L, p2L, dL) := closest_pair_recursive(L)
    (p1R, p2R, dR) := closest_pair_recursive(R)

    if dL < dR:
        (best_p, best_q, δ) := (p1L, p2L, dL)
    else:
        (best_p, best_q, δ) := (p1R, p2R, dR)

    # Build the strip of width 2·δ around mid_x
    strip := [point in P : |point.x − mid_x| < δ]
    sort strip by y-coordinate

    # Strip merge: each point compared with at most 7 forward neighbours
    for i := 0 to length(strip) − 1:
        for j := i + 1 to min(i + 7, length(strip) − 1):
            if strip[j].y − strip[i].y >= δ:
                break                                  # early termination
            d := distance(strip[i], strip[j])
            if d < δ:
                δ := d
                (best_p, best_q) := (strip[i], strip[j])

    return (best_p, best_q, δ)

The “magic 7” is justified rigorously in §5; informally, it’s the maximum number of points that can be packed in a δ × 2δ rectangle while remaining pairwise ≥ δ apart.

A common engineering optimisation: maintain Px (input sorted by x) and Py (input sorted by y) once, and split each into Lx, Rx, Ly, Ry in linear time per recursive call rather than re-sorting the strip each level. This drops the merge work from O(n log n) per level to O(n) per level, restoring the clean T(n) = 2T(n/2) + O(n) recurrence. (CLRS §33.4 presents this version.)

4. Python Implementation

import math
from typing import Sequence
 
Point = tuple[float, float]
 
 
def _dist_sq(p: Point, q: Point) -> float:
    """Squared Euclidean distance — avoid the sqrt in the inner loop for speed."""
    dx, dy = p[0] - q[0], p[1] - q[1]
    return dx * dx + dy * dy
 
 
def _brute_force(pts: Sequence[Point]) -> tuple[Point, Point, float]:
    """O(k^2) for k ≤ 3 base case."""
    n = len(pts)
    best = (pts[0], pts[1], _dist_sq(pts[0], pts[1]))
    for i in range(n):
        for j in range(i + 1, n):
            d2 = _dist_sq(pts[i], pts[j])
            if d2 < best[2]:
                best = (pts[i], pts[j], d2)
    return best
 
 
def closest_pair(points: list[Point]) -> tuple[Point, Point, float]:
    """
    Return (p, q, d) where (p, q) is the closest pair and d is the distance.
 
    O(n log n) time, O(n) space.
    """
    n = len(points)
    if n < 2:
        raise ValueError("need at least 2 points")
 
    # Sort once by x and once by y, in O(n log n) total
    Px = sorted(points, key=lambda p: (p[0], p[1]))
    Py = sorted(points, key=lambda p: (p[1], p[0]))
 
    p, q, d2 = _rec(Px, Py)
    return p, q, math.sqrt(d2)
 
 
def _rec(Px: list[Point], Py: list[Point]) -> tuple[Point, Point, float]:
    n = len(Px)
    if n <= 3:
        return _brute_force(Px)
 
    mid = n // 2
    midpoint = Px[mid]
    mid_x = midpoint[0]
 
    # Split Px into L and R in O(n) (already sorted by x).
    Lx = Px[:mid]
    Rx = Px[mid:]
 
    # Split Py into Ly (those whose x ≤ mid_x and are in left half) and Ry — O(n).
    # We use set membership for correctness when there are duplicates on the divider.
    left_set = set(id(p) for p in Lx)
    Ly = [p for p in Py if id(p) in left_set]
    Ry = [p for p in Py if id(p) not in left_set]
 
    p1, q1, d1 = _rec(Lx, Ly)
    p2, q2, d2 = _rec(Rx, Ry)
 
    if d1 < d2:
        best_p, best_q, best_d2 = p1, q1, d1
    else:
        best_p, best_q, best_d2 = p2, q2, d2
 
    delta = math.sqrt(best_d2)
 
    # Build strip of width 2·delta around mid_x, preserving y-sorted order.
    strip = [p for p in Py if abs(p[0] - mid_x) < delta]
 
    # Strip merge: for each point, compare to at most 7 forward neighbours.
    m = len(strip)
    for i in range(m):
        j = i + 1
        while j < m and (strip[j][1] - strip[i][1]) < delta:
            d2 = _dist_sq(strip[i], strip[j])
            if d2 < best_d2:
                best_d2 = d2
                best_p, best_q = strip[i], strip[j]
                delta = math.sqrt(best_d2)
            j += 1
            # The "j ≤ i + 7" bound is implicit in the y-cutoff (proof in §5).
 
    return best_p, best_q, best_d2

A few engineering notes:

  1. Squared distance everywhere internal. √ is expensive and unnecessary for comparisons — d²₁ < d²₂ iff d₁ < d₂. Take the square root only at the API boundary.
  2. id() for set membership. Two distinct Point tuples can be == (e.g., duplicate coordinates), but we want to split Py based on which physical point belongs to which half. id() (Python’s “address”) gives us identity rather than equality. (In a stricter implementation we’d index points and pass indices.)
  3. The while j < m and ... < delta inner loop is the strip merge. The “7” bound is implicit — the loop terminates when y diverges by delta, which (per the packing argument) bounds j − i by ≤ 7.

5. Complexity — Recurrence and Packing Proof

5.1 The Recurrence

T(n) = 2 · T(n/2) + O(n)

Two halves of size n/2, plus O(n) for the strip merge (and for splitting Py). By the Master Theorem Case 2 (a = 2, b = 2, f(n) = O(n), so n^log_b(a) = n matches), T(n) = Θ(n log n).

By hand: at recursion depth k, there are 2^k subproblems each of size n / 2^k, doing O(n / 2^k) strip-merge work each — total O(n) per level. There are log₂ n levels. Total work O(n log n). The initial sort is also O(n log n), so overall O(n log n) dominates.

5.2 Why the Strip Merge Is O(n) — The Packing Argument

This is the crucial step. Without bounding the inner-loop comparisons, the strip merge could be O(m²) where m = |strip|, blowing the recurrence up to T(n) = 2T(n/2) + O(n²), which the master theorem solves to Θ(n²) — no better than brute force.

Claim. For each point p in the y-sorted strip, the number of subsequent strip points q (with q.y ≥ p.y) such that q.y − p.y < δ is at most 7.

Proof. Consider p at position (x₀, y₀). The relevant region is the rectangle: {(x, y) : |x − mid_x| < δ, y₀ ≤ y < y₀ + δ} — width 2δ, height δ.

Inside this 2δ × δ rectangle, every pair of points from the left half is at least δ apart (since δ_L ≥ δ) and similarly for the right half. So we can pack at most:

  • 4 points from the left half (in the δ × δ left subrectangle, at the corners and possibly two in the middle, where any tighter packing would force two points within δ).
  • 4 points from the right half (symmetric).

Total: at most 8 points in the rectangle, including p itself. So at most 7 subsequent points to check.

The rigorous packing bound: in a δ × δ square, at most 4 points can be ≥ δ apart. Proof: divide the square into four (δ/2) × (δ/2) sub-squares. Each sub-square has diagonal δ/√2 < δ, so any two points in the same sub-square are < δ apart. Pigeonhole: 5 points in 4 sub-squares ⇒ two in the same sub-square ⇒ contradiction. So ≤ 4. The 2δ × δ rectangle is two δ × δ squares glued ⇒ ≤ 8 points total ⇒ ≤ 7 forward neighbours of p.

The literature is split on the constant

CLRS uses 7 (counting “points within δ in y”) and the analysis above gives 7. Some textbooks write 8 (including p itself) or 15 (using a coarser packing argument). All these constants are equivalent up to the bookkeeping convention; what matters is the bound is constant, which is what makes the strip-merge work O(m) ⊆ O(n). Shamos and Hoey’s 1975 paper does not focus on the exact constant — it focused on the existence of an O(n log n) algorithm. (CLRS §33.4)

5.3 Space

O(n) — for the sorted lists Px, Py, and the strip buffer at each recursion level (though levels share — at any moment, the live strip plus call stack is O(n) total).

6. Variants and Optimisations

6.1 Higher Dimensions

The same divide-and-conquer template works in d dimensions: split by hyperplane perpendicular to one axis, recurse, and check a “slab” of width 2δ on either side. The strip-merge constant grows exponentially in d (roughly c · 3^d), so the bound becomes T(n, d) = 2T(n/2, d) + O(n · 3^d) = O(n log n) for fixed d. For variable d, the curse of dimensionality kicks in. (Bentley 1976 has the multi-dimensional version.)

6.2 Randomised O(n) Algorithm

Khuller and Matias (1995, Information and Computation 118.1) show a randomised algorithm with expected O(n) time using floor-function-based hashing — for each point compute the bucket of side δ it lies in, hash into a grid, and check only neighbouring buckets. The randomisation handles the chicken-and-egg of “you need δ to bucket but you don’t know δ yet” by sampling. Beats Shamos-Hoey asymptotically, but the constant is large and the deterministic version is preferred in practice.

6.3 KD-Tree-Based Heuristic

Build a KD-Tree of all points, then for each point query its nearest neighbour. Worst case O(n²) but average-case O(n log n) and very fast for “well-distributed” data. Used in practice (scipy.spatial, sklearn) where worst-case bounds matter less than constant factors.

6.4 Sweep-Line Variant

Bentley and Ottmann’s sweep-line technique can also solve closest-pair in O(n log n) by maintaining a y-sorted active set as the sweep advances in x, again pruning by the δ-strip insight. Conceptually different but same asymptotic. Useful when you’re already doing other sweep-line geometry (e.g., computing line-segment intersections).

7. Diagram — The Strip Merge

flowchart TD
    P[All points sorted by x] --> S{Split at midpoint}
    S --> L[Left half<br/>recurse → δL]
    S --> R[Right half<br/>recurse → δR]
    L --> M[δ = min(δL, δR)]
    R --> M
    M --> ST[Build strip:<br/>points with |x − mid_x| < δ<br/>sorted by y]
    ST --> CHK[For each strip point,<br/>compare with ≤ 7 forward neighbours<br/>O(n) total]
    CHK --> ANS[Final answer:<br/>min(δ, best strip pair)]

    style ST stroke:#5af
    style CHK stroke:#5af

What this diagram shows. Top of the recursion: the points are split by the median x-coordinate. Each half is recursively closed-paired, yielding δ_L and δ_R. We take δ = min(δ_L, δ_R) as the current best. The remaining work is the strip merge (the two highlighted blue boxes): keep only points within distance δ (horizontally) of the dividing line; sort by y; for each strip point, compare against forward neighbours until the y-gap exceeds δ. The packing argument bounds the inner loop’s iterations by a constant (≤ 7), so the strip merge is O(n). Combining recursion (2T(n/2)) with linear merge (O(n)) gives T(n) = 2T(n/2) + O(n) = Θ(n log n).

8. Pitfalls

  1. Returning δ not the pair. Many tutorials show the algorithm returning just the distance. In interviews you usually want the pair, which means threading the (p, q) pair through every recursive call. Not hard but often forgotten.
  2. Resorting the strip by y at every level. This makes each level O(n log n) and bumps the recurrence to T(n) = 2T(n/2) + O(n log n) = O(n log² n) — still better than O(n²) but a notch worse than the optimal O(n log n). The fix: maintain Py (input sorted by y) once at the top, then split Py into Ly, Ry in linear time per recursion (using the dividing line’s mid_x to partition). This restores O(n) work per level.
  3. Off-by-one in the ”≤ 7” inner-loop bound. Don’t hardcode for j in range(i+1, min(i+8, m)). Instead use a while loop that terminates when strip[j].y − strip[i].y ≥ δ. The 7 bound is a consequence of the y-cutoff, not the loop’s stopping condition — using it as the stopping condition is brittle and obscures the algorithm.
  4. Comparing distances with < versus <=. Two pairs can tie; either is fine but be consistent. If the input has duplicate points (d = 0), the algorithm should handle that — the closest pair is the duplicate pair with d = 0.
  5. Floating-point distance comparisons. Squared distance avoids the sqrt and stays in integer arithmetic if input coordinates are integers, eliminating one source of round-off bugs. Always prefer d² < δ² to d < δ internally.
  6. Splitting by index vs by x-value. If the median x-coordinate is shared by many points (a vertical line of input), splitting strictly by index puts some “left” points to the right of mid_x. The strip-merge is robust to this (it tests |x − mid_x| < δ, not membership in a half) but the recursion balance can degrade. For uniformly random inputs it’s a non-issue.
  7. Forgetting the base case. if n ≤ 3 is the right base case — the brute force on 3 points is C(3, 2) = 3 distance checks, faster than another level of recursion. Going down to n = 1 would loop forever (no pairs to compare).
  8. Confusing this with k-d tree nearest-neighbour search. The closest-pair problem asks for the closest pair globally; nearest-neighbour search asks for “given a query point q, find the closest input point to q.” The algorithms are related but distinct.

9. Common Interview Problems

ProblemSourceNotes
Find Closest Pair (any platform)ClassicThe vanilla problem; expected O(n log n)
K Closest Points to OriginLeetCode 973Different problem (use Top K Elements / Quickselect); often confused with this
Closest Pair in 1DEasier variantSort, scan adjacent — O(n log n) for sort, O(n) for scan
K Closest Pairs of PointsLeetCode 658 (closely related)Heap-based, not D&C
Erect the Fence (Convex Hull)LeetCode 587Different geometry problem; D&C also applies
Minimum Spanning Tree on PointsApplicationClosest-pair gives a lower bound edge for MST; used in heuristics
Hierarchical ClusteringApplicationClosest-pair drives the merging step in single-linkage clustering

10. Open Questions

  • In practice on n = 10⁶ random points, how does the Shamos-Hoey D&C compare with the KD-tree heuristic? The asymptotics favour D&C; constant factors and cache behaviour favour KD-tree. (Empirical answer depends on input distribution.)
  • Khuller-Matias’s randomised O(n) algorithm beats Shamos-Hoey asymptotically but is rarely implemented. What’s the constant overhead in practice?
  • Is there an O(n log log n) or O(n α(n)) deterministic algorithm? The current best deterministic bound is O(n log n) (Fortune & Hopcroft 1979 showed this is tight in the algebraic decision tree model).
  • How does the ”≤ 7” bound generalise to higher dimensions? The proof above gives a 4 + 4 = 8 packing in 2D; 3D requires a 3D packing argument. Bentley 1976 has the multi-dimensional version.

11. See Also

Graph View

Backlinks