Sorting I: elementary sorts and merge sort

Sorting is the most studied problem in algorithms, for good reasons. It is useful on its own, it is a first step in many other algorithms, such as the duplicate check in lesson 1 and Kruskal's algorithm in lesson 10, and it shows almost every design technique in a small space. This lesson studies two simple quadratic sorts, selection sort and insertion sort, and then merge sort, the first divide-and-conquer algorithm in the course, which sorts in \(\Theta(n \log n)\) time. Along the way it introduces the vocabulary used to compare sorting algorithms: comparisons, inversions, stability and extra space.

The sorting problem and how we measure sorts

The input is a sequence of \(n\) items, each with a key drawn from a totally ordered set, and the output is a permutation of the input with keys in non-decreasing order. The algorithms in this lesson and the next only learn about keys by comparing two of them, so they are called comparison sorts, and we usually count key comparisons as the basic operation. Two other properties matter in practice. A sort is in place if it uses only \(O(1)\) extra memory beyond the input array, or \(O(\log n)\) for a recursion stack. A sort is stable if items with equal keys keep their original relative order, which matters when you sort records by one field after another.

Selection sort

Selection sort repeatedly finds the smallest remaining item and swaps it into place. After pass \(i\), positions 0 to \(i\) hold the \(i+1\) smallest items in order, which is the loop invariant. Finding the minimum of the \(n - i\) remaining items takes \(n - i - 1\) comparisons, so the total is \((n-1) + (n-2) + \cdots + 1 = n(n-1)/2\), whatever the input. Selection sort is therefore \(\Theta(n^2)\) in every case. Its one virtue is that it makes at most \(n - 1\) swaps, which matters only when moving items is far more expensive than comparing them.

Insertion sort

Insertion sort builds a sorted prefix one item at a time. It takes the next item and shifts larger items in the prefix one position right until it finds where the new item belongs, the way many people sort a hand of cards. Its invariant is that before processing position \(i\), the prefix \(a[0..i-1]\) contains the original first \(i\) items in sorted order. Insertion sort is stable, in place, simple, and very fast on small or nearly sorted inputs, which is why production sorting routines use it for short runs.

def insertion_sort(a):
    a = list(a)
    comps = shifts = 0
    for i in range(1, len(a)):
        key = a[i]
        j = i - 1
        while j >= 0:
            comps += 1
            if a[j] <= key:          # <= keeps equal keys in order: stable
                break
            a[j + 1] = a[j]          # shift right
            shifts += 1
            j -= 1
        a[j + 1] = key
    return a, comps, shifts

def inversions_brute(a):
    return sum(1 for i in range(len(a)) for j in range(i + 1, len(a)) if a[i] > a[j])

import random
random.seed(204)
n = 200
cases = {'sorted': list(range(n)), 'reversed': list(range(n, 0, -1)),
         'random': random.sample(range(n), n)}
for name, data in cases.items():
    out, c, s = insertion_sort(data)
    assert out == sorted(data)
    print(f'{name:>8}: comparisons {c:>6}, shifts {s:>6}, inversions {inversions_brute(data):>6}')

  sorted: comparisons    199, shifts      0, inversions      0
reversed: comparisons  19900, shifts  19900, inversions  19900
  random: comparisons   9753, shifts   9561, inversions   9561

The output reveals the exact cost of insertion sort. The number of shifts always equals the number of inversions in the input, where an inversion is a pair of positions \(i < j\) with \(a[i] > a[j]\). Each shift swaps one adjacent out-of-order pair and so removes exactly one inversion. A sorted list has no inversions and needs \(n - 1\) comparisons; a reversed list has the maximum \(n(n-1)/2\) inversions. A random permutation has, on average, half of all pairs inverted, \(n(n-1)/4\), so the average case is also \(\Theta(n^2)\). That argument proves something stronger: any sort that only swaps adjacent items, including bubble sort, needs \(\Omega(n^2)\) swaps on average.

Merge sort: divide and conquer

Merge sort follows the divide-and-conquer pattern from lesson 3. Divide the list into two halves. Conquer by sorting each half recursively. Combine by merging the two sorted halves into one sorted list. A list of length 0 or 1 is already sorted, which is the base case. All the real work happens in the merge step.

To merge two sorted lists \(L\) and \(R\), keep a pointer at the front of each. Repeatedly compare the two front items, move the smaller one to the output and advance its pointer. When one list runs out, append the rest of the other. Each comparison outputs one item, and the last item is output without a comparison, so merging lists of total length \(n\) takes at most \(n - 1\) comparisons. Taking from the left list when the keys are equal makes merge sort stable.

def merge(left, right, counter):
    out = []
    i = j = 0
    while i < len(left) and j < len(right):
        counter[0] += 1
        if left[i] <= right[j]:          # take from the left on ties: stable
            out.append(left[i]); i += 1
        else:
            out.append(right[j]); j += 1
    out.extend(left[i:])
    out.extend(right[j:])
    return out

def merge_sort(a, counter):
    if len(a) <= 1:
        return list(a)
    mid = len(a) // 2
    return merge(merge_sort(a[:mid], counter), merge_sort(a[mid:], counter), counter)

import math, random
random.seed(4)
print(merge_sort([38, 27, 43, 3, 9, 82, 10], [0]))
for n in (2**10, 2**14, 2**17):
    data = [random.random() for _ in range(n)]
    c = [0]
    assert merge_sort(data, c) == sorted(data)
    print(n, c[0], round(n * math.log2(n)), round(c[0] / (n * math.log2(n)), 3))

[3, 9, 10, 27, 38, 43, 82]
1024 8929 10240 0.872
16384 208469 229376 0.909
131072 2062401 2228224 0.926

On random inputs merge sort uses about \(n \log_2 n - 1.25n\) comparisons, which is 87 to 93 percent of \(n \log_2 n\) for these sizes and creeps toward 100 percent as \(n\) grows. The saving comes from merges that end with one list exhausted early. The worst case is at most \(n \lceil \log_2 n \rceil - n + 1\) comparisons, and even the best case, an already sorted list, needs about half of \(n \log_2 n\), because each merge still compares until one side runs out. Unlike insertion sort, merge sort does not adapt much to sorted input.

Why merge sort is correct

Correctness follows by strong induction on the length \(n\). Lists of length 0 or 1 are sorted, so the base case returns a correct result. For \(n \ge 2\), both halves are shorter than \(n\), so by the inductive hypothesis the recursive calls return them sorted. The merge then produces a sorted list containing exactly the items of both halves, by its own loop invariant: the output holds the smallest items seen so far in sorted order, and the front items of \(L\) and \(R\) are the smallest remaining ones in their lists. So merge sort returns a sorted permutation of its input.

Analysis and space

The running time satisfies \(T(n) = 2T(n/2) + \Theta(n)\), which lesson 3 solved with a recursion tree and the master theorem: \(T(n) \in \Theta(n \log n)\). This holds in every case, because the splitting does not depend on the data. The price is memory. The straightforward version creates new lists at every level, but even a careful array implementation needs a temporary buffer of \(n\) items for merging, so merge sort uses \(\Theta(n)\) extra space. In-place merging is possible but complicated and slow, which is one reason quicksort and heapsort, studied next, remain important.

Bottom-up merge sort and other variations

Merge sort need not be recursive. Bottom-up merge sort treats the input as \(n\) runs of length 1, merges adjacent pairs into runs of length 2, then 4, and so on, doubling the run length on each pass. There are \(\lceil \log_2 n \rceil\) passes, each costing \(\Theta(n)\), giving the same \(\Theta(n \log n)\) bound without any recursion. Natural merge sort goes further and starts from the runs already present in the data, so nearly sorted input needs very few passes. Python's built-in sort, Timsort, uses this idea together with insertion sort for short runs, which makes it \(\Theta(n)\) on already sorted data and \(\Theta(n \log n)\) in the worst case.

Merge sort is also the algorithm of choice for data that does not fit in memory. External merge sort reads chunks that fit in RAM, sorts each one, writes the sorted runs to disk and then merges many runs at once with a priority queue (lesson 5). Each pass reads and writes the data sequentially, which suits disks and networked storage far better than the random access other sorts need.

Counting inversions with merge sort

The number of inversions measures how different two rankings are, for example two people's rankings of the same films. Counting them by checking all pairs takes \(\Theta(n^2)\) time. A small change to merge sort counts them in \(\Theta(n \log n)\). During a merge, whenever an item from the right half is output before items remaining in the left half, it forms an inversion with every one of those remaining left items, because they are larger and came earlier. Adding that count at each step totals the inversions that cross between the halves; the recursive calls count those within each half.

def sort_count(a):
    """Return (sorted list, number of inversions) in Theta(n log n) time."""
    if len(a) <= 1:
        return list(a), 0
    mid = len(a) // 2
    left, x = sort_count(a[:mid])
    right, y = sort_count(a[mid:])
    out, i, j, cross = [], 0, 0, 0
    while i < len(left) and j < len(right):
        if left[i] <= right[j]:
            out.append(left[i]); i += 1
        else:
            out.append(right[j]); j += 1
            cross += len(left) - i       # right[j] precedes all remaining left items
    out += left[i:] + right[j:]
    return out, x + y + cross

import random
random.seed(7)
print(sort_count([3, 1, 4, 2])[1])
data = random.sample(range(2000), 2000)
fast = sort_count(data)[1]
slow = sum(1 for i in range(len(data)) for j in range(i + 1, len(data)) if data[i] > data[j])
print(fast, slow, fast == slow)

3
998169 998169 True

The count for 3, 1, 4, 2 agrees with the hand calculation above, and on 2,000 random items the fast method matches the brute-force count. The random count is close to the expected \(n(n-1)/4 \approx 999{,}500\). Recommendation systems use exactly this measure, often normalized as Kendall's tau, to compare rankings.