Closest pair of points in Python

Why this matters

A point is a pair of numbers, such as (3, 4), that marks a place on a flat map. Programs use points all the time. A delivery app might store the location of every driver. A game might store the position of every enemy. A map app might need to know which two cafes are closest to each other.

The simple way to find the closest pair is to measure the distance between every two points and keep the smallest. That is correct, but it gets slow. With 1,000 points it does about half a million distance checks. With 1,000,000 points it does about half a trillion. This lesson shows a method called divide and conquer that does the same job with far fewer checks. Divide and conquer means you split a problem into smaller copies of itself, solve each copy, and then combine the answers. The method in this lesson checks only a small band of points around the middle, not every pair.

The previous lesson in this module, Binary search on the answer, used a similar idea of cutting a problem in half. The next lesson, Master theorem, shows how to read the running time of a method like this one from a formula.

The idea from the first principle

Start with the problem stated in plain words. You have a list of points. Each point has an x value (how far right it is) and a y value (how far up it is). The distance between two points is the straight-line distance between them. You get it with Pythagoras: take the difference in x, take the difference in y, square both, add them, and take the square root. For (0, 0) and (3, 4) the differences are 3 and 4, the squares are 9 and 16, the sum is 25, and the square root is 5.

A brute force method is the one that tries every option. For points, that means checking every pair. With 4 points there are 6 pairs. With 6 points there are 15 pairs. In general, n points have n times (n minus 1) divided by 2 pairs. For 1,000 points that is 499,500 pairs.

Now the idea that makes divide and conquer work. Sort the points by x. Draw an imaginary vertical line through the middle point of that sorted list. The points are now split into a left group and a right group, and each group has about half the points.

The closest pair is in one of three places:

1. Both points are in the left group.
2. Both points are in the right group.
3. One point is in the left group and the other is in the right group.

The first two cases are smaller copies of the same problem. You solve each one with the same method. The third case is the hard one, because the two points cross the dividing line. Call the smaller of the two answers from the halves d. Any crossing pair that beats d must have both points within d of the dividing line. Points farther away than d cannot form a better pair, so you can ignore them. The points within d of the line form a narrow vertical band. Call that band the strip.

The strip is the key to the speed. Inside the strip you still need to check pairs, but there are very few useful pairs to check. The next section shows why.

Why the strip only needs a few checks

Sort the strip by y. Take one point in the strip. Which other points could be closer than d to it? Only points whose y value is within d of its y value, because a bigger difference in y already means a bigger distance. So for each point, you look at the points that come after it in y order, and you stop as soon as the y difference reaches d.

Why does that stop quickly? Picture the strip as a tall, narrow box. It is 2d wide, because it reaches d to the left of the line and d to the right. It is d tall, because you only care about points within d of the current point. Cut that box into 8 small squares, each d / 2 on a side. Two points in the same small square would be at most d / 2 times the square root of 2 apart, which is less than d. But two points on the same side of the dividing line are always at least d apart, because d is the smallest distance found inside each half. So each small square holds at most one point, and the whole box holds at most 8 points.

That means for any point in the strip, at most 7 points after it in y order can be within d of it. So the strip check does a fixed, small amount of work per point.

The method in steps

Here is the method written out in order:

1. Make a copy of the points and sort the copy by x.
2. If the group has 3 or fewer points, check every pair. Then sort the group by y and return the smallest distance.
3. Otherwise, find the middle position in the x-sorted group. Remember the x value there. Call it mid_x.
4. Solve the left half and the right half with this same method. Each call returns two things: the closest distance in its half, and its points sorted by y.
5. Merge the two y-sorted halves into one list sorted by y. Merging is the step from merge sort: look at the front of each list and take the smaller one, until both are empty.
6. Build the strip: keep the points whose x is less than d away from mid_x. Because the merged list is in y order, the strip is already sorted by y.
7. For each point in the strip, check the next points in y order until the y gap reaches d. Update d when you find a smaller distance.
8. Return d.

Python has no built-in point type, so the code uses a dataclass. A dataclass is a short class that stores a few named values. The line @dataclass(frozen=True) makes the values read-only. Point(3, 4) makes a point, and p.x reads its x value.

The standard library module heapq has a function called merge. It joins sorted lists into one sorted list, which is the merge step from merge sort. It produces its results one at a time, so list(...) collects them. key=lambda p: p.y tells merge and sorted to compare points by their y value.

The Java version of this method sorts one shared array in place. The Python version below returns a new list at each level. It is easier to follow, and it uses more memory, because each level builds new lists from its points.

import heapq
import math
from dataclasses import dataclass
 
 
@dataclass(frozen=True)
class Point:
    x: int
    y: int
 
 
distance_checks = 0
 
 
def distance(a: Point, b: Point) -> float:
    global distance_checks
    distance_checks += 1
    dx = a.x - b.x
    dy = a.y - b.y
    return math.sqrt(dx * dx + dy * dy)
 
 
def brute_force(pts: list[Point]) -> float:
    best = math.inf
    for i in range(len(pts)):
        for j in range(i + 1, len(pts)):
            best = min(best, distance(pts[i], pts[j]))
    return best
 
 
def solve(pts_by_x: list[Point]) -> tuple[float, list[Point]]:
    # Returns the closest distance and the same points sorted by y.
    n = len(pts_by_x)
    if n <= 3:
        best = brute_force(pts_by_x)
        return best, sorted(pts_by_x, key=lambda p: p.y)
 
    mid = n // 2
    mid_x = pts_by_x[mid].x
    best_left, left = solve(pts_by_x[:mid])
    best_right, right = solve(pts_by_x[mid:])
    best = min(best_left, best_right)
 
    # Both halves are sorted by y. heapq.merge joins them into one list sorted by y.
    merged = list(heapq.merge(left, right, key=lambda p: p.y))
 
    strip = [p for p in merged if abs(p.x - mid_x) < best]
    for a in range(len(strip)):
        for b in range(a + 1, len(strip)):
            if strip[b].y - strip[a].y >= best:
                break
            best = min(best, distance(strip[a], strip[b]))
    return best, merged
 
 
def closest_fast(points: list[Point]) -> float:
    pts = sorted(points, key=lambda p: p.x)  # sorted returns a new list, so the caller's list is not changed
    best, _ = solve(pts)
    return best

Several details in this code matter.

The n <= 3 case is the base case. A base case is the smallest input the function answers directly, without calling itself. With 3 points, checking every pair is only 3 checks, so brute force is the cheapest choice there.

The slices pts_by_x[:mid] and pts_by_x[mid:] make two new lists, one for each half. A slice is a range of positions in a list. The first slice holds the positions from 0 up to mid, not including mid. The second holds the rest. Each half is then solved on its own.

The _ in best, _ = solve(pts) means "ignore this value". The function returns two values, and this caller needs only the first.

The global line inside distance and run says that distance_checks is the counter defined at the top of the file. Without it, Python would treat the assignment as a new local variable, and the counter outside the function would never change.

The line if strip[b].y - strip[a].y >= best: break makes the strip check fast. Once the y gap is at least best, every later point in the strip is even farther away in y, so the loop can stop.

math.inf is positive infinity. Any real distance is smaller than it, so it works as a starting value for best.

The method takes a list of points and sorts a copy. The caller's list is left alone. A function that sorts its input should say so, or sort a copy.

Worked example 1: a pair that crosses the line

Start with four points: A (0, 0), B (4, 9), C (5, 11), and D (9, 0). Look at the points by eye first. The two closest points look like B and C, because they are near each other in the top middle. Check that with the method.

Step 1 sorts by x. The order is already A, B, C, D, since 0, 4, 5, 9 are in increasing order.

Step 2 handles the base case. There are 4 points, so the method does not stop yet. It takes mid = 2, so the left half is A and B, and the right half is C and D. The value mid_x is the x of the point at position 2, which is C, so mid_x is 5.

Step 3 solves the left half. A and B are two points, so the base case checks their one pair. The distance from (0, 0) to (4, 9) is the square root of 16 plus 81, which is the square root of 97, about 9.849. The left half is sorted by y: A (y = 0), then B (y = 9).

Step 4 solves the right half. The distance from C (5, 11) to D (9, 0) is the square root of 16 plus 121, which is the square root of 137, about 11.705. The right half is sorted by y: D (y = 0), then C (y = 11).

So d is 9.849 after the two halves.

Step 5 merges the halves by y. The merged order is A (y = 0), D (y = 0), B (y = 9), C (y = 11).

Step 6 builds the strip. Check each point's x distance from 5 against d, which is 9.849. A is 5 away, which is less than 9.849, so A is in the strip. D is 4 away, so D is in. B is 1 away, so B is in. C is 0 away, so C is in. All four points are in the strip, in y order: A, D, B, C.

Step 7 checks pairs in the strip.

• Start with A. Compare it with D. The y gap is 0, which is less than d. The distance is 9, so d becomes 9.
• Still with A, compare it with B. The y gap is 9, which is not less than 9, so the loop stops for A.
• Next is D. Compare it with B. The y gap is 9, which again stops the loop.
• Next is B. Compare it with C. The y gap is 2, which is less than 9. The distance from (4, 9) to (5, 11) is the square root of 1 plus 4, which is the square root of 5, about 2.236. So d becomes 2.236.
• C is the last point, so there is nothing to compare.

The answer is 2.236, the pair B and C. The halves alone would have said 9.849, which is wrong. The strip check found the pair that crossed the line.

Run this example with the code above:

points = [Point(0, 0), Point(4, 9), Point(5, 11), Point(9, 0)]
print(closest_fast(points))

Output:

2.23606797749979

The method made 4 distance calculations in total: one in the left base case, one in the right base case, and two in the strip. The brute force method would have made 6. For 4 points that is a small saving, but the next example shows why the saving grows.

Worked example 2: a strip with several candidates

Now use six points: (0, 0), (1, 9), (4, 5), (6, 0), (7, 8), and (9, 2). Sorted by x they are already in that order.

With 6 points, mid is 3. The left half is the first three points, (0, 0), (1, 9), and (4, 5). The right half is (6, 0), (7, 8), and (9, 2). mid_x is 6, the x of the point at position 3.

The left base case checks three pairs. The distances are about 9.055, 6.403, and 5. So the left answer is 5, and the left half sorted by y is (0, 0), then (4, 5), then (1, 9).

The right base case checks three pairs. The distances are about 8.062, 3.606, and 6.325. So the right answer is 3.606, the square root of 13, from the pair (6, 0) and (9, 2). The right half sorted by y is (6, 0), then (9, 2), then (7, 8).

So d is 3.606 after the halves.

The merge step produces the list (0, 0), (6, 0), (9, 2), (4, 5), (7, 8), (1, 9), all sorted by y.

The strip keeps points whose x is less than 3.606 away from 6. The x distances are 6 for (0, 0), 0 for (6, 0), 3 for (9, 2), 2 for (4, 5), 1 for (7, 8), and 5 for (1, 9). So the strip is (6, 0), (9, 2), (4, 5), (7, 8). Notice that (9, 2) and (4, 5) are on opposite sides of the line, so this strip check looks at a pair that crosses it.

Now check pairs in the strip, in y order.

• From (6, 0): compare with (9, 2). The y gap is 2, which is less than 3.606. The distance is the square root of 13, which is not less than 3.606, so d stays the same. Compare with (4, 5). The y gap is 5, which is not less than 3.606, so the loop stops.
• From (9, 2): compare with (4, 5). The y gap is 3, which is less than 3.606. The distance is the square root of 34, about 5.831, which does not help. Compare with (7, 8). The y gap is 6, so the loop stops.
• From (4, 5): compare with (7, 8). The y gap is 3. The distance is the square root of 18, about 4.243, which does not help.
• From (7, 8): there are no points after it.

The strip check ran three distance calculations and found nothing better. The answer is 3.606, from the right half. The total for the whole run is 9 distance calculations: 3 in the left base case, 3 in the right base case, and 3 in the strip. Brute force would make 15.

points = [Point(0, 0), Point(1, 9), Point(4, 5), Point(6, 0), Point(7, 8), Point(9, 2)]
print(closest_fast(points))

Output:

3.605551275463989

This example shows the most common result. The strip check often runs, and often finds nothing better, but it has to run to be sure. The check is cheap because the loop stops early.

Worked example 3: 1,000 random points

The third example is too large for hand work, so the code below runs both methods on 1,000 random points. The points come from random.Random(42), a random number generator with a fixed seed. A seed is a starting number. The same seed gives the same sequence of numbers, so every run makes the same points and prints the same numbers. The points differ from the ones the Java version of this lesson makes, because Python and Java use different generators.

The distance function counts its own calls in distance_checks, so the program can compare the two methods by how many distances each one calculates. The code below uses the run function and a main block. Add import random to the top of the file, and put this code below the code from the method section.

def run(label: str, points: list[Point]) -> None:
    global distance_checks
    distance_checks = 0
    fast = closest_fast(points)
    fast_checks = distance_checks
 
    distance_checks = 0
    slow = brute_force(points)
    slow_checks = distance_checks
 
    print(f"{label}: n={len(points)} fast={fast:.4f} brute={slow:.4f} "
          f"fast_checks={fast_checks} brute_checks={slow_checks}")
 
 
if __name__ == "__main__":
    rng = random.Random(42)
    big = [Point(rng.randrange(10000), rng.randrange(10000)) for _ in range(1000)]
    run("random1000", big)

The output of the program is:

random1000: n=1000 fast=10.0000 brute=10.0000 fast_checks=1131 brute_checks=499500

Both methods agree on the answer, 10.0000. The divide and conquer method made 1,131 distance calculations. Brute force made 499,500. That is about 440 times fewer. The gap grows with the input size, which is the whole point of the method.

The count stays low because each point in a strip is compared only until the y gap reaches the current best distance, and then the loop stops.

Why the running time is O(n log n)

The method splits the points into two halves, solves both, and then does work proportional to the number of points. Write the total work for n points as T(n). It is the work for the two halves, plus the merge and the strip work for all n points:

T(n) = 2 * T(n / 2) + n

The n at the end is the merge and the strip, which look at each point a fixed number of times. The notation O(n log n) is a short way to say how that total grows. Read O(...) as an upper limit: the work never grows faster than the expression inside. The term log n is the number of times you can halve n before you reach 1. For 1,000 points that is about 10, since 2 to the power 10 is 1,024. Each of those roughly 10 levels of splitting does about n units of merge and strip work, so the total is about n times log n.

The next lesson, Master theorem, shows how to read a formula like this one and get the n log n answer without drawing the whole tree.

The sorting at the start also costs O(n log n). So the whole method takes O(n log n) time, which is much better than the O(n squared) time of the brute force method, where n squared is n times n.

Common mistakes

Forgetting the strip. If the code takes the smaller of the two half answers and stops, it misses pairs that cross the line. The first worked example shows this. If solve returns only min(best_left, best_right), the first worked example prints 9.848857801796104, which is about 9.849. The true answer is 2.236. A quick test that uses four points on both sides of a line catches this kind of bug.

Sorting the strip by x or by nothing. The early break only works when the strip is sorted by y. If the strip is in the wrong order, the loop can stop before it reaches a close point, and the method returns a wrong answer that looks plausible.

Sorting inside every level. Some people sort the strip from scratch at each level of recursion. The answer stays correct, but the work grows to about O(n log squared n). The merge step avoids this by reusing the y order that the recursive calls already built.

Overflow does not happen, but / does. The Java version of this method can overflow an int when it squares a large difference, and the result is a wrong number. Python's int has no size limit, so 50000 * 50000 is an exact 2500000000, and math.sqrt of it is 50000.0. The mistake that replaces it is the division operator. / always gives a float, and a float cannot be used as a list index.

pts = [Point(0, 0), Point(1, 9), Point(4, 5), Point(6, 0)]
mid = len(pts) / 2
print(pts[mid])

The last line of the error that Python prints is:

TypeError: list indices must be integers or slices, not float

Use // for every position in a list, such as len(pts) // 2.

Sorting objects without telling Python how to compare them. sort() compares objects with <. A dataclass does not define < unless you ask for it, so this program fails:

pts = [Point(3, 1), Point(1, 2)]
pts.sort()

The last line of the error is:

TypeError: '<' not supported between instances of 'Point' and 'Point'

The fix in this lesson is to pass a key function, such as sorted(points, key=lambda p: p.x). The key tells sort which number to compare.

Practice

Each exercise has its answer under it. Try the work on paper first.

Exercise 1. Find the distance between (0, 0) and (3, 4). Use the Pythagoras steps from this lesson. Then check your answer with math.dist((0, 0), (3, 4)).

Show the answer

The differences are 3 and 4. Their squares are 9 and 16. The sum is 25. The square root of 25 is 5. The distance is 5. math.dist((0, 0), (3, 4)) returns 5.0.

Exercise 2. Four points are (1, 1), (4, 5), (5, 1), and (9, 9). Find the closest pair by checking pairs by eye, then give the distance.

Show the answer

The closest pair is (1, 1) and (5, 1). They have the same y, so their distance is just the x difference, which is 4. The next closest is (4, 5) and (5, 1), at about 4.123. The distance is 4.

Exercise 3. A method splits 8 points into two halves of 4 points. Each half returns a closest distance of 6.0. The strip holds 3 points. In the worst case, how many distance calculations can the strip step make for those 3 points?

Show the answer

The strip has 3 points, so there are at most 3 pairs to check, which is 3 times 2 divided by 2. Each pair is one distance calculation, so the worst case is 3. The early break can only make it smaller.

Exercise 4. The lesson says that one point in the strip is compared with at most 7 later points. Explain in your own words why 8 points is the most a strip box can hold, and why that gives 7.

Show the answer

The strip box is 2d wide and d tall. Cut it into 8 squares of side d / 2. Two points in one square would be less than d apart, but two points on the same side of the line are always at least d apart. So each square holds at most one point, and the box holds at most 8. A point uses one of those 8 slots, so at most 7 others can follow it in the strip.

Exercise 5. Points (0, 0), (4, 9), (5, 11), and (9, 0) were used in the first worked example. Suppose someone removes the strip step and returns only the smaller of the two half answers. What does the method return, and what is the true answer?

Show the answer

It returns 9.849, the closest distance inside the left half, which is the pair (0, 0) and (4, 9). The true answer is 2.236, the pair (4, 9) and (5, 11). The true pair crosses the dividing line, so only the strip step can find it.

Summary

• Brute force checks every pair, which takes O(n squared) time and grows too fast for large inputs.
• Divide and conquer sorts by x, solves the left and right halves, and checks a strip around the dividing line.
• The strip holds at most 8 points per box, so each point needs only a few comparisons, and the early break stops the loop sooner.
• Merging the halves by y keeps the strip sorted without sorting it again at each level, which gives O(n log n) total time.
• In Python, use // for positions, and pass a key function to sorted or sort when the objects have no <.