Median of Medians (BFPRT): Finding the k-th Smallest in Worst-Case Linear Time

Finding the k-th smallest of n values is called the selection problem, and the median is its best-known case. Sorting everything and reading off position k works, but it fixes the order of every element just to get one value.

In their 1973 paper “Time Bounds for Selection,” Blum, Floyd, Pratt, Rivest, and Tarjan gave an algorithm, PICK, whose number of comparisons is at most a linear function of n for any k. Because of how it picks its pivot it is usually called “median of medians,” and after its authors it is also known as BFPRT. This post follows that paper (Journal of Computer and System Sciences, vol. 7, pp. 448–461, in the copy published by co-author Rivest) to cover how PICK is built, why it runs in linear time, and the question “why groups of five?” Later improvements and specific library implementations are out of scope.

One pivot decides the cost

The obvious way to find the k-th value is to reuse quicksort’s partition step. Pick a pivot m, compare it with every other value, throw away the side that cannot contain the k-th value, and repeat on what is left. The paper describes PICK as quite similar to Hoare’s FIND algorithm, except that the element m used for partitioning is chosen more carefully.

The weak point of that approach is the pivot. What follows is this post’s own calculation. If the pivot is always the minimum or maximum of what remains, each partition removes only one element. The comparison count then becomes roughly n + (n − 1) + …, which grows to about n²/2. The post on quicksort’s worst case shows how pivot choice produces the worst case in quicksort. PICK spends extra comparisons on choosing the pivot so that every partition is guaranteed to remove at least a fixed fraction of the elements.

The steps of PICK

The paper describes PICK in general form with three parameters, c, d, and b. To select the i-th smallest element of a set S, one pass goes like this:

1. Split S into n/c columns of length c and sort each column. Let T be the set of the d-th smallest elements of the columns, and call PICK recursively to find the b-th element of T. That element is m.
2. Compare m with every element whose relation to m is not yet known, which gives m’s rank in S.
3. If the rank is i, m is the answer and the algorithm stops. If the rank is greater than i, discard the elements greater than or equal to m; if it is smaller, discard the elements less than or equal to m (and reduce i by the number discarded). Go back to step 1.

With d set to the middle of a column and b to the middle of T, m is “the median of the column medians,” which is where the name comes from. The paper states that each subset discarded this way contains at least one quarter of the remaining elements.

Why it is linear

The paper writes P(n) for the maximum number of comparisons, picks c = 21, d = 11, and b = n/42, and obtains this recurrence (equation 6):

P(n) ≤ 66n/21 + P(n/21) + n + P(31n/42)

Here is what each term means. 66 is h(21), the number of comparisons the Ford–Johnson method needs to sort 21 values, so the first term is the cost of sorting all n/21 columns. P(n/21) is the recursive call that finds the median of the n/21 column medians, n bounds the comparisons between m and the remaining elements, and P(31n/42) is the next pass on whatever is left after discarding. From this, the paper proves by induction that P(n) ≤ 58n/3.

Where 31n/42 comes from can be read off Figure 1 of the paper. The following is this post’s calculation based on that figure, ignoring columns that do not divide evenly and off-by-one differences (the paper also says it omits the details for n not divisible by c). Half of the column medians, meaning n/42 columns, are at most m. Each column is sorted, so in those columns the bottom 11 elements are all at most m, which gives at least 11n/42 elements known to be at most m. By the same reasoning, at least 11n/42 elements are known to be at least m. Whichever side is discarded, at least 11n/42 elements go, so at most 31n/42 remain.

What makes the total linear is the combined size of the two recursive calls: n/21 + 31n/42 = 33n/42, which is less than n. The total amount of work handed to the next level shrinks by a fixed ratio each time, so the costs over all levels add up to a convergent geometric series, and the whole thing is proportional to n. That explanation is this post’s; the paper sums it up as PICK discarding a significant fraction of S on each pass, at a cost proportional to the number of elements discarded.

Where the number five comes from

If you learned median of medians as “split into groups of five,” the numbers in the original paper may come as a surprise. The basic proof uses columns of 21, and the tuned version, PICK1, uses columns of 15. Five appears only in one sentence at the end of the proof: “Note that we must have c ≥ 5 for PICK to run in linear time.” The paper gives no separate proof for that sentence.

The reason shows up when you generalize the previous section to any c. This is also this post’s calculation. Take c to be odd, c = 2k + 1, so that each column has a single middle element. Half of the columns each contribute k + 1 elements known to be at most m, so about (k + 1)n/(2c) elements are known to be at most m, and at most n − (k + 1)n/(2c) remain. Writing the sizes of the two recursive calls as fractions of n:

Column size c	Median recursion	Next pass	Sum
3	1/3	2/3	1
5	1/5	7/10	9/10
7	1/7	5/7	6/7
21	1/21	31/42	11/14

With c = 3 the sum is exactly 1, so the total number of elements handled does not shrink from one level to the next, and the geometric-series argument above does not give a linear bound. Five is the smallest odd size for which the sum drops below 1, and explanations that use groups of five are using that minimum. This table is only a simple bound computed with PICK’s structure left as is; it says nothing about variants that change the structure.

You cannot make c as large as you like either. Sorting the columns in step 1 costs n·h(c)/c, and the paper notes that c must be bounded above by a constant for PICK to run in linear time. The paper picked c = 21 for the basic proof specifically to minimize the bound from its recurrence (5). Plugging c = 5 into the same recurrence (7 comparisons to sort five values, 7n/10 for the next pass) gives a bound of 24n, larger than the 58n/3 for c = 21. That comparison is this post’s calculation. Five is the minimum that works, not the best choice for keeping the constant small.

Constants and limitations stated in the paper

The paper’s conclusion is that selection can be done in linear time in the worst case. With the tuned PICK1, it proves that no more than about 5.43n comparisons are needed for any i (exactly 5.4305n, with the final 5 repeating), and PICK2 gives better bounds when i is near either end. On the lower-bound side, it shows that computing a median requires at least 3n/2 − 2 comparisons.

The paper also adds its own caveats about practice. It calls PICK1 “full of red tape” and says it could not in practice be implemented efficiently, and the authors expected the constants in both the upper and lower bounds to be improvable considerably. The paper measures cost in comparisons and assumes the n input numbers are distinct.

A Python translation

The code below translates PICK with groups of five into Python. The paper assumes distinct values, but this code counts the values equal to the pivot separately so that duplicates work too. k is zero-based.

def select(a, k):
    """Return the k-th smallest value in a (k is zero-based)."""
    if len(a) <= 5:
        return sorted(a)[k]
    # group into fives and collect each group's median (the last group may be shorter)
    groups = [sorted(a[i:i + 5]) for i in range(0, len(a), 5)]
    medians = [g[len(g) // 2] for g in groups]
    # recursively find the median of the medians and use it as the pivot
    pivot = select(medians, len(medians) // 2)
    lo = [x for x in a if x < pivot]
    hi = [x for x in a if x > pivot]
    n_eq = len(a) - len(lo) - len(hi)
    if k < len(lo):
        return select(lo, k)
    if k < len(lo) + n_eq:
        return pivot
    return select(hi, k - len(lo) - n_eq)


print(select([9, 1, 8, 2, 7, 3, 6, 4, 5, 0, 5], 5))  # 5

This code is meant to show the structure. It builds new lists at every step, and unlike the paper it does not skip comparisons whose outcome is already known. Its constant factor in comparisons differs from the paper’s analysis, but the structure in which the two recursive calls together cover fewer than n elements is the same.

Things to check

• Kind of guarantee: PICK’s linear bound is a deterministic worst-case guarantee, not a probabilistic or average-case one.
• Group size: the original paper does not fix five. Five is the smallest odd size for which the linear argument works with this structure; the paper’s basic proof uses 21 and the tuned version uses 15.
• Constants: the bound proved in the paper is about 5.43n comparisons, and the authors themselves thought the constant could be improved considerably.
• Duplicates: the paper assumes distinct values. An implementation has to handle values equal to the pivot separately.