Median of Medians, in depth: the 30 percent guarantee, the recurrence with constants, why groups of five, and introselect in real libraries

Why selection needs a guarantee

Quickselect picks a pivot, partitions the array into elements below, equal to and above it, and recurses into only the side that contains rank k. With a random pivot the expected work is linear, about 3.4n comparisons for the median in the textbook analysis. The worst case is not. If every pivot lands near an end, each round removes only a few elements and the total work is n + (n - 1) + (n - 2) + ..., which is quadratic. A fixed rule such as first element or median of three can be defeated by a crafted input, and a random rule can still be unlucky.

Whether this matters depends on who controls the input. A batch job over your own data can live with expected bounds. A service computing percentiles over user-supplied values, or a library used by unknown callers, cannot: the worst case is exactly what an attacker will find. The fix is a pivot whose quality is proved rather than hoped for.

The algorithm

The algorithm to select rank k from an array A of n elements is:

1. If n is small (say at most five), sort A directly and return A[k].
2. Split A into ceil(n/5) groups of five, the last group possibly shorter.
3. Sort each group (insertion sort, at most 10 comparisons for five elements; 7 is optimal) and take its median.
4. Recursively call select on the list of group medians to find its median, M. This is the pivot.
5. Partition A three ways around M into L (less), E (equal) and G (greater).
6. If k < |L| recurse into L with k; if k < |L| + |E| return M; otherwise recurse into G with k - |L| - |E|.

There are two recursive calls, and that is the subtle part. The first, on the medians, is on an input of size n/5 and exists only to choose the pivot. The second is the usual quickselect step. The total stays linear only because both calls are on inputs that are a constant fraction smaller than n and the two fractions add up to less than one.

Why the pivot is good

Picture the groups as columns, each sorted top to bottom, and the columns ordered by their medians, as in the figure. The pivot M is the median of the middle row. Every column to the left of M's column, including M's own, has a median at most M, and in each of those columns the median and the two elements above it are at most M. So at least three elements per column are at most M, for about half the columns.

Made precise: there are g = ceil(n/5) groups. At least ceil(g/2) group medians are at most M. Discard M's own group and the possibly short last group; each remaining group contributes at least three elements at most M. That gives at least 3(ceil(g/2) - 2), which is at least 3n/10 - 6 elements at most M, and by symmetry at least 3n/10 - 6 elements at least M. Whichever side the recursion enters therefore has at most n - (3n/10 - 6) = 7n/10 + 6 elements.

This is a worst-case statement. It holds for every input and every k, which is what quickselect with a random pivot cannot offer. Note what it does not say: the pivot is not near the true median, only guaranteed to sit somewhere between roughly the 30th and 70th percentile.

The recurrence, solved

Let T(n) be the worst-case number of steps. Sorting the groups and partitioning are linear, say at most an for a constant a. The two recursive calls give:

T(n) <= T(ceil(n/5)) + T(7n/10 + 6) + a*n        for n >= 140
T(n) <= constant                                  for n <  140

Guess T(n) <= cn and substitute. Using ceil(n/5) <= n/5 + 1, the right side is at most c(n/5 + 1) + c(7n/10 + 6) + an = 9cn/10 + 7c + an. That is at most cn when cn/10 - 7c >= an, that is c >= 10an / (n - 70). For n >= 140 the fraction n / (n - 70) is at most 2, so c = 20a works, and choosing c large enough also covers the base cases. Hence T(n) = O(n).

The key quantity is 1/5 + 7/10 = 9/10. Each level of the recursion tree does at most 9/10 of the work of the level above, so the total is a geometric series bounded by ten times the top level. The large constant hidden there is the price of the guarantee, and it is why nobody uses this algorithm as their first choice.

Why groups of five

Nothing in the method fixes the group size. With groups of size g (odd), the guaranteed fraction on each side is about (g + 1)/(4g) of the elements, and the recurrence becomes T(n) <= T(n/g) + T(n(3g - 1)/(4g)) + an.

Five is the smallest odd size for which the straightforward analysis gives linear time, and larger groups trade a smaller recursion on the medians for more comparisons to sort each group. The straightforward analysis is not the last word: Chen and Dumitrescu (2015) showed that groups of three or four can also run in linear time if the algorithm is modified, for example by applying the grouping step twice before recursing. Five is the default because its proof is short, not because it is a law.

An implementation that counts its work

The implementation below is written for clarity and measurement rather than speed: it uses lists instead of in-place index ranges, and routes every comparison through a counter so the cost can be measured. The same select function runs with either pivot rule.

import random

COMPARISONS = 0

def less(a, b):
    global COMPARISONS
    COMPARISONS += 1
    return a < b

def small_sort(xs):
    xs = list(xs)
    for i in range(1, len(xs)):
        x, j = xs[i], i - 1
        while j >= 0 and less(x, xs[j]):
            xs[j + 1] = xs[j]
            j -= 1
        xs[j + 1] = x
    return xs

def partition3(xs, pivot):
    lo, eq, hi = [], [], []
    for x in xs:
        if less(x, pivot):
            lo.append(x)
        elif less(pivot, x):
            hi.append(x)
        else:
            eq.append(x)
    return lo, eq, hi

def mom_pivot(xs):
    medians = [small_sort(xs[i:i + 5])[(len(xs[i:i + 5]) - 1) // 2]
               for i in range(0, len(xs), 5)]
    return select(medians, (len(medians) - 1) // 2, mom_pivot)

def random_pivot(xs):
    return random.choice(xs)

def select(xs, k, choose_pivot):
    # Return the k-th smallest (0-based) element of xs.
    while True:
        if len(xs) <= 5:
            return small_sort(xs)[k]
        p = choose_pivot(xs)
        lo, eq, hi = partition3(xs, p)
        if k < len(lo):
            xs = lo
        elif k < len(lo) + len(eq):
            return p
        else:
            k -= len(lo) + len(eq)
            xs = hi

Two details matter. The partition is three-way: if it only split into less-than and greater-or-equal, an array of identical values would never shrink and the loop would not terminate. And the second recursion is written as a loop, so the stack depth comes only from the pivot recursion on the medians, which shrinks by a factor of five each time and is therefore logarithmic.

Worked example: 25 values

Take the 25 values 51, 29, 60, 93, 16, 19, 78, 22, 56, 84, 17, 74, 37, 14, 21, 65, 63, 18, 40, 85, 64, 89, 25, 38, 96 and ask for the median, rank k = 12 counting from zero.

The medians are 51, 56, 21, 63, 64. Five values is a base case, so the recursive call sorts them to 21, 51, 56, 63, 64 and returns 56. Partitioning all 25 values around 56 gives L with 13 elements (every value below 56), E = {56} and G with 11 elements. Since k = 12 is less than |L| = 13, the answer is in L, and the next round selects rank 12 of those 13 values, which is their maximum, 51.

The structure promised at least 8 elements on each side of the pivot; the actual split was 13 and 11. The guarantee is a floor, not a forecast.

Measured cost against random pivots

Running the code above on uniformly random floats, selecting the median, gave these comparisons per element. Note that partition3 spends one comparison on elements below the pivot and two on the rest, about 1.5 per element.

Median of medians settles at about ten comparisons per element whatever the input, which is what a linear worst-case bound looks like. Random pivots average about 5n here (the textbook 3.4n times 1.5 for this partition), but single runs scatter widely around that. Wall-clock favours random pivots further, because median of medians copies more data.

Introselect and what libraries ship

Production libraries do not choose between the two; they combine them. Musser's introselect (1997) runs quickselect with a cheap pivot rule, watches how fast the problem shrinks, and switches to a guaranteed method if it is not shrinking fast enough. Typical inputs pay the quickselect price and adversarial ones pay a bounded penalty.

• Rust's slice::select_nth_unstable documents an introselect based on ipnsort whose fallback is median of medians using Tukey's ninther for pivot selection, which, in the words of the documentation, guarantees linear runtime for all inputs.
• NumPy's numpy.partition uses kind='introselect' by default, and its documentation lists the worst case as O(n).
• C++ std::nth_element is only required by the standard to be linear on average. libstdc++ implements introselect with a heap-based fallback, which bounds the worst case at O(n log n) rather than O(n). If you need a hard linear bound in C++, you have to supply it.

The general pattern, cheap path first with a proven fallback behind a trip-wire, is the same one introsort uses for sorting, described in Quicksort, in depth. For the related problem of returning the top k items in order, see Partial Sort, in depth.

Failure modes

• Two-way partition with duplicates. Splitting into less-than and greater-or-equal never shrinks an array of equal keys. Always partition three ways or use a scheme that splits equal keys.
• Off-by-one in the median of a short group. The last group may have one to four elements. Take index (len - 1) // 2 consistently; a mismatch silently breaks the 3n/10 bound for small n.
• Recursion depth. Recursing on both the medians and the chosen side doubles stack use; turn the second recursion into a loop as above.
• Using it as the main path. At ten-plus comparisons per element plus extra data movement it loses to quickselect on almost every real input. Use it as the fallback in introselect.
• Assuming the library is worst-case linear. Check the documented guarantee; nth_element in common C++ implementations is not.
• NaN in floats. NaN is neither less nor greater than anything, so it lands in E and corrupts the result. Filter or order NaN first.

Trade-offs

What to do next

1. Decide whether your selection inputs can be adversarial or pathological; if so, an expected bound is not enough.
2. Look up the documented worst case of the selection routine your language ships, and note it next to the call site.
3. If you write your own, implement introselect: quickselect with a shrink check, falling back to median of medians when the check fails.
4. Use a three-way partition and test arrays of all-equal values, sorted, reverse-sorted and organ-pipe inputs.
5. Add a property test comparing your selector with sorting on random small arrays, including duplicates and every k.
6. Count comparisons on your real data, as above, before deciding the guaranteed path is too slow.