Divide and conquer beyond sorting

Merge sort and quicksort are divide-and-conquer algorithms, but the technique reaches far beyond sorting. This lesson applies it to three classic problems. Selection finds the \(k\)-th smallest item without sorting, in linear time. The closest pair of points among \(n\) points in the plane is found in \(\Theta(n \log n)\) instead of \(\Theta(n^2)\). And Karatsuba's algorithm multiplies two \(n\)-digit numbers with fewer than the \(n^2\) digit multiplications taught in school. Each example shows a different way the combine or divide step can be made clever, and each analysis uses the recurrence tools from lesson 3.

The selection problem

Given \(n\) items and a rank \(k\), the selection problem asks for the item that would be at position \(k\) if the items were sorted. Finding the minimum or maximum takes \(n - 1\) comparisons, and finding the median, with \(k = \lceil n/2 \rceil\), is a common statistics task. Sorting first solves any selection in \(\Theta(n \log n)\), but that does more work than needed: we want one item, not a whole order. The question is whether selection can be done in linear time.

Quickselect

Quickselect borrows quicksort's partition. Choose a pivot, partition, and look at where the pivot lands. If its position is \(k\), it is the answer. If \(k\) is smaller, recurse only on the left side; if larger, recurse only on the right side, adjusting \(k\). Unlike quicksort, quickselect makes one recursive call, not two. With a random pivot, the expected size of the side we keep is at most \(3n/4\) about half the time, and the expected running time satisfies a recurrence whose solution is linear. A careful analysis shows the expected number of comparisons is at most \(4n\), and about \(3.4n\) for the median.

import random

def quickselect(a, k, counter):
    """Return the k-th smallest (0-based) of list a; count comparisons."""
    a = list(a)
    lo, hi = 0, len(a) - 1
    while True:
        if lo == hi:
            return a[lo]
        r = random.randint(lo, hi)
        a[r], a[hi] = a[hi], a[r]
        p, i = a[hi], lo
        for j in range(lo, hi):
            counter[0] += 1
            if a[j] < p:
                a[i], a[j] = a[j], a[i]
                i += 1
        a[i], a[hi] = a[hi], a[i]
        if k == i:
            return a[i]
        elif k < i:
            hi = i - 1
        else:
            lo = i + 1

random.seed(6)
for n in (1000, 10000, 100000):
    data = random.sample(range(10 * n), n)
    total = 0
    for trial in range(20):
        c = [0]
        med = quickselect(data, n // 2, c)
        total += c[0]
    assert med == sorted(data)[n // 2]
    print(f'n={n:>6}: median found, average comparisons {total / 20 / n:.2f} n')

n=  1000: median found, average comparisons 3.23 n
n= 10000: median found, average comparisons 3.20 n
n=100000: median found, average comparisons 3.44 n

The average comparisons per item stay roughly constant as \(n\) grows by a factor of 100, the signature of linear time. Quickselect is what most libraries use for selection; C++ calls it nth_element. Its worst case, like quicksort's, is \(\Theta(n^2)\), but only with astronomically bad luck.

Median of medians: linear in the worst case

Can we guarantee linear time with no randomness? In 1973 Blum, Floyd, Pratt, Rivest and Tarjan showed how, by choosing a pivot that is provably near the middle. Split the items into groups of 5 and find the median of each group with a constant number of comparisons. Recursively find the median of those \(\lceil n/5 \rceil\) medians and use it as the pivot. At least half of the group medians are at most the pivot, and each of those has two more items in its group that are at most it, so at least about \(3n/10\) items are at most the pivot. The same holds on the other side, so each recursive call keeps at most about \(7n/10\) items.

The master theorem does not apply to this recurrence, because the two subproblems have different sizes, so we use substitution. Guess \(T(n) \le d\,n\). Then \(T(n) \le d\,n/5 + 7d\,n/10 + c\,n = (9d/10 + c)\,n\), which is at most \(d\,n\) when \(d \ge 10c\). The key is that \(1/5 + 7/10 = 9/10 < 1\): the subproblems shrink in total, so the work forms a decreasing geometric series. With groups of 3, the fractions would be \(1/3 + 2/3 = 1\) and the bound would fail, which is why the algorithm uses groups of 5. In practice its constant is large, so it is mainly used as a fallback, as in introselect, and as a theoretical guarantee.

The closest pair of points

Given \(n\) points in the plane, find two at minimum distance. Checking every pair costs \(\Theta(n^2)\). The divide-and-conquer solution, due to Shamos and Hoey, runs in \(\Theta(n \log n)\). Sort the points by \(x\)-coordinate once. Split them by a vertical line into a left half and a right half of equal size and find the closest pair in each half recursively; let \(\delta\) be the smaller of the two distances. The closest pair overall is either one of those, or a pair with one point on each side. Any such pair closer than \(\delta\) must lie within the vertical strip of width \(2\delta\) around the dividing line.

The clever part is the combine step. Sort the points in the strip by \(y\)-coordinate. For each point, only the next few points in that order can be within \(\delta\) of it. Picture a \(2\delta \times \delta\) rectangle above the point, divided into eight squares of side \(\delta/2\). Two points in the same square would be at most \(\delta/\sqrt{2} < \delta\) apart and on the same side of the line, contradicting the definition of \(\delta\), so each square holds at most one point, and at most 7 other points need checking. The combine step is therefore linear, provided the points are kept sorted by \(y\), which the recursion can do by merging as in merge sort. The recurrence is \(T(n) = 2T(n/2) + \Theta(n)\), so \(T(n) \in \Theta(n \log n)\).

import math, random

def brute(pts):
    best = math.inf
    for i in range(len(pts)):
        for j in range(i + 1, len(pts)):
            best = min(best, math.dist(pts[i], pts[j]))
    return best

def closest(px, counter):
    """px sorted by x. Returns (distance, points sorted by y)."""
    n = len(px)
    if n <= 3:
        d = math.inf
        for i in range(n):
            for j in range(i + 1, n):
                counter[0] += 1
                d = min(d, math.dist(px[i], px[j]))
        return d, sorted(px, key=lambda p: p[1])
    mid = n // 2
    xline = px[mid][0]
    dl, ly = closest(px[:mid], counter)
    dr, ry = closest(px[mid:], counter)
    d = min(dl, dr)
    py, i, j = [], 0, 0                       # merge by y, as in merge sort
    while i < len(ly) and j < len(ry):
        if ly[i][1] <= ry[j][1]: py.append(ly[i]); i += 1
        else: py.append(ry[j]); j += 1
    py += ly[i:] + ry[j:]
    strip = [p for p in py if abs(p[0] - xline) < d]
    for a in range(len(strip)):
        for b in range(a + 1, min(a + 8, len(strip))):   # at most 7 following points
            counter[0] += 1
            d = min(d, math.dist(strip[a], strip[b]))
    return d, py

random.seed(42)
for n in (100, 1000, 4000):
    pts = [(random.random(), random.random()) for _ in range(n)]
    c = [0]
    d, _ = closest(sorted(pts), c)
    line = f'n={n:>5}: distance computations {c[0]:>6}'
    if n <= 1000:
        line += f', brute force {n * (n - 1) // 2:>6}, agree: {math.isclose(d, brute(pts))}'
    print(line)

n=  100: distance computations    614, brute force   4950, agree: True
n= 1000: distance computations  10898, brute force 499500, agree: True
n= 4000: distance computations  61753

For 1,000 points the divide-and-conquer version computes about 11,000 distances where brute force computes nearly half a million, and both find the same answer. Closest-pair computations appear in collision detection, clustering and air traffic control, and the strip argument is a model for many geometric algorithms.

Multiplying large integers: Karatsuba

The school method multiplies two \(n\)-digit numbers with \(n^2\) single-digit multiplications. For a long time this was believed optimal. In 1960, Anatoly Karatsuba, then a student, found a divide-and-conquer method that is faster. Write each number as a high half and a low half: \(x = x_1 B^{m} + x_0\) and \(y = y_1 B^{m} + y_0\), where \(B\) is the base and \(m = n/2\). Then \(xy = x_1 y_1 B^{2m} + (x_1 y_0 + x_0 y_1) B^m + x_0 y_0\), which seems to need four half-size products.

Karatsuba's insight is that the middle term can be obtained from one extra product: \(x_1 y_0 + x_0 y_1 = (x_1 + x_0)(y_1 + y_0) - x_1 y_1 - x_0 y_0\). So three half-size multiplications suffice, plus some additions and shifts, which cost \(\Theta(n)\). The recurrence \(T(n) = 3T(n/2) + \Theta(n)\) has \(a = 3\), \(b = 2\), \(d = 1\), and since \(\log_2 3 \approx 1.585 > 1\), the master theorem gives \(\Theta(n^{1.585})\).

def school(x, y):
    """Digit-by-digit multiplication; returns (product, single-digit multiplications)."""
    xs, ys = str(x), str(y)
    return x * y, len(xs) * len(ys)

def karatsuba(x, y, counter):
    if x < 10 or y < 10:
        counter[0] += 1                        # one single-digit multiplication
        return x * y
    m = max(len(str(x)), len(str(y))) // 2
    x1, x0 = divmod(x, 10 ** m)
    y1, y0 = divmod(y, 10 ** m)
    z2 = karatsuba(x1, y1, counter)
    z0 = karatsuba(x0, y0, counter)
    z1 = karatsuba(x1 + x0, y1 + y0, counter) - z2 - z0
    return z2 * 10 ** (2 * m) + z1 * 10 ** m + z0

import random
random.seed(3)
print(karatsuba(1234, 5678, [0]))
for digits in (64, 256, 1024):
    x = random.randrange(10 ** (digits - 1), 10 ** digits)
    y = random.randrange(10 ** (digits - 1), 10 ** digits)
    c = [0]
    assert karatsuba(x, y, c) == x * y
    print(f'{digits:>5} digits: school {school(x, y)[1]:>8} digit products, Karatsuba {c[0]:>7}, n^1.585 = {round(digits ** 1.585):>7}')

7006652
   64 digits: school     4096 digit products, Karatsuba    1003, n^1.585 =     729
  256 digits: school    65536 digit products, Karatsuba    9519, n^1.585 =    6562
 1024 digits: school  1048576 digit products, Karatsuba   85993, n^1.585 =   59064

At 1,024 digits, Karatsuba uses about a twelfth of the school method's digit products, and the gap widens with \(n\). Python's own integers switch to Karatsuba above a size threshold, and libraries such as GMP go further with Toom–Cook splitting into more parts and with fast Fourier transform methods, reaching \(O(n \log n)\) in a 2019 result by Harvey and van der Hoeven. The same idea of trading multiplications for additions gives Strassen's matrix algorithm from lesson 3.

Designing divide-and-conquer algorithms

The table suggests a design recipe. Find a way to split the problem into smaller instances of itself. Make the divide and combine steps as cheap as possible, ideally linear. Reduce the number of subproblems if you can, as Karatsuba did, or recurse on only one of them, as quickselect does. Then write the recurrence and solve it before you write the code, because the recurrence tells you whether the design is worth implementing. When subproblems overlap heavily, so the same subproblem is solved many times, divide and conquer becomes inefficient; dynamic programming, in CMPUT-304, handles that case by storing answers.