Randomised algorithms

Sometimes a random choice works as well as a clever one. Randomised algorithms make random choices as they run, and their guarantees are about every input: on any input, the time is small on average, or the chance of a wrong answer is tiny and can be made as small as you like. This is the last lesson of the course.

Two kinds of random algorithm

• Always right, time random (called Las Vegas): quicksort and quickselect with random pivots (SORT-3, SORT-8). The answer is always correct; the guarantee is on the expected time, the average over the random choices (MATH-7).
• Always fast, rarely wrong (called Monte Carlo): a Bloom filter (HASH-8), Freivalds’s check, the Miller–Rabin test, Karger’s cut. The time is fixed; the guarantee is a bound on the chance of a wrong answer.

The randomness is in the method, not in the input: no input is bad for it, because the random choices are made after the input is fixed. A Monte Carlo method run k times independently has a chance of error that shrinks like the k-th power of one run’s chance. This is called amplification.

Random pivots

Why does quicksort with a random pivot make about 2n ln n comparisons on average, for every input? Name the values by their sorted order, z₁ < z₂ < … < zₙ. Two values zᵢ and zⱼ (i < j) are compared only if one of them is the first pivot chosen among zᵢ, …, zⱼ: until then they stay in the same part, and once a value between them is chosen, they go to different parts. Each of these j − i + 1 values is equally likely to be first, so the chance is 2/(j − i + 1). Adding up over all pairs (the expected count is the sum of the chances, MATH-7) gives less than 2n(1 + 1/2 + … + 1/n), which is about 2n ln n.

Checking a product

Is A · B = C for three n × n matrices? Multiplying A by B costs about n³ steps (n · n · n). Freivalds (1977): pick a random vector r of 0s and 1s and compare A · (B · r) with C · r. Each is two matrix-vector products, O(n²).

def same_product(A, B, C, rounds, rand):
    # rand() returns 0 or 1; a wrong C escapes each round with chance at most 1/2
    n = len(A)
    for _ in range(rounds):
        r = [rand() for _ in range(n)]
        Br = [sum(B[i][k] * r[k] for k in range(n)) for i in range(n)]
        ABr = [sum(A[i][k] * Br[k] for k in range(n)) for i in range(n)]
        Cr = [sum(C[i][k] * r[k] for k in range(n)) for i in range(n)]
        if ABr != Cr:
            return False                 # certainly A·B != C
    return True                          # A·B == C, wrong with chance <= 2**-rounds

Why at most 1/2. If A · B ≠ C, the difference D = A · B − C has some entry dᵢⱼ ≠ 0. Entry i of D · r is dᵢⱼ · rⱼ plus a sum s that does not involve rⱼ. Fix every other entry of r; of the two choices of rⱼ, at most one makes dᵢⱼ · rⱼ + s = 0. So D · r = 0 has chance at most 1/2. A “no” answer is always right; a “yes” after 20 rounds is wrong with chance at most 2⁻²⁰, about one in a million.

Prime or not: Miller–Rabin

MATH-5’s Fermat test checks aⁿ⁻¹ mod n = 1: true for every a when n is prime. But Carmichael numbers such as 561 pass it for every a that shares no factor with them. Miller (1976) and Rabin (1980) added a second fact about primes:

(If x² − 1 = (x − 1)(x + 1) is a multiple of a prime, one of the two factors is.) So write n − 1 = d · 2ˢ with d odd, and look at aᵈ, a²ᵈ, a⁴ᵈ, …, aⁿ⁻¹, each the square of the one before, all modulo n. For a prime, the list either starts at 1 or reaches n − 1 before it reaches 1. If instead a value other than ±1 squares to 1, or the last value is not 1, then n is certainly not prime and a is a witness. The figure shows a prime, the Carmichael number 561, and 2047, which fools base 2.

How often it errs. For an odd composite n, at most a quarter of the bases from 1 to n − 1 are strong liars (bases that pass). So k random bases all pass with chance at most (1/4)ᵏ: 20 bases give less than one chance in 10¹² (a million million). Each base costs one pow(a, d, n) and s squarings, O(log n) multiplications (MATH-3).

For numbers of a bounded size, no luck is needed at all. Searches by computer found the smallest composite that fools every one of the first k prime bases. Below 3.3 × 10²⁴ (33 followed by 23 zeros), the 13 bases 2, 3, 5, …, 41 catch every composite, so the test becomes exact. This is derandomisation: a fixed set of choices proved to work for every input in range. The practice problem uses it.

Karger’s minimum cut

A cut of an undirected graph splits its vertices into two non-empty sides; its size is the number of edges between them. FLOW-4 finds the smallest cut between two given vertices with a flow. Karger (1993) finds the smallest cut overall with random choices:

1. pick an edge at random and contract it: merge its two ends into one vertex, keep all other edges (repeats too), drop loops;
2. repeat until two vertices are left; the edges between them are a cut.

Why it often works. Fix a smallest cut, of size c. Every vertex has at least c edges (else its own edges would be a smaller cut), so with v vertices left there are at least c · v / 2 edges, and the chance of contracting one of the cut’s c edges is at most 2/v. The cut survives all contractions with chance at least (1 − 2/n)(1 − 2/(n − 1)) … (1 − 2/3), and this product equals 2 / (n(n − 1)). Repeating about n² ln n times and keeping the smallest cut found misses with chance about 1/n². Karger and Stein made it faster by sharing the early contractions between runs.

Testing without luck

A test that passes or fails by chance is useless (MATH-7). The practice pages of this course test random code in three ways: pass the random function in, so a test can script its choices (MATH-P7); check a property that holds for every outcome; or, as here, remove the randomness by a proof. In your own code, random.seed(…) makes runs repeatable, random.Random(seed) gives a generator of your own, and secrets gives numbers nobody can predict, for passwords and keys.

Now try HARD-P5, Prime or not.