RND Randomness Technique
Las Vegas and Monte Carlo
Also known as: Las Vegas and Monte Carlo algorithms · Las Vegas algorithm · Monte Carlo algorithm · Freivalds' algorithm · Schwartz–Zippel lemma · restarts · ZPP, RP, BPP
A randomised algorithm can pay for speed in two ways. Las Vegas always answers right, but its running time is random; Monte Carlo runs fast and predictably but sometimes errs, with a probability we choose. One turns into the other by a timer stop or by checking the answer, and the error can be made as small as you like by repetition.
A server has sent back a blurred photograph: it replaced each point with the sum of its nine neighbours along the row, multiplying the picture’s matrix by a blur matrix. Did it compute honestly? Recomputing is 16.8 million multiplications. There is a way to check in two hundred thousand, but it can be wrong: say “all right” when there is an error. It is wrong at most half the time, and twenty checks in a row will miss an error with a probability of less than one in a million. Spoil a point in the picture on the right and check for yourself.
One spoilt point is caught in about half the checks, and that is no coincidence but the exact probability: the check sees an error in column j only if the random vector is one in that column. Press “a sly damage”: two errors in one row with opposite signs that cancel in the sum of the row. An ordinary checksum will not notice them, while the check still catches them half the time.
The idea
Randomness in an algorithm can be spent in two ways. A Las Vegas algorithm always answers right, and only its running time is random. A Monte Carlo algorithm runs in predictable time but may err, with a probability we control. The two kinds turn into each other: Las Vegas stopped by a timer becomes Monte Carlo, and Monte Carlo whose answer can be checked quickly becomes Las Vegas. Repetition makes the Monte Carlo error as small as you like, and restarts cut off the long tails of Las Vegas.
The names were invented as a pair: “Monte Carlo” appeared in the 1940s among physicists modelling neutrons (see the article on the Monte Carlo method), and “Las Vegas” was proposed in 1979 by László Babai, one casino to match the other. The basics, randomised quicksort as Las Vegas and edge contraction as Monte Carlo, are taken apart in the course chapter “Let’s flip a coin”; here it is about what can be done with these algorithms.
A Monte Carlo error can have a direction. Freivalds’ check errs one way: if it said “there is an error”, there certainly is, while “all right” can be false. The Miller–Rabin primality test (its “composite” is always true) and the Bloom filter (its “no” is always true) are built the same way. There are also algorithms with two-sided error, which can be wrong either way, and they have to be amplified differently.
Freivalds’ check
Given n × n matrices A, B and C, is it true that AB = C? Multiplying A and B by the definition costs n³ multiplications. Rūsiņš Freivalds proposed in 1977 taking a random vector r of zeros and ones and comparing A(Br) with Cr. Multiplying a matrix by a vector costs n² operations, and there are three such multiplications.
If AB = C, the check always answers “right”. If AB ≠ C and each coordinate of r is independently 0 or 1 with probability ½, the check wrongly answers “right” with probability at most ½.
The first part is obvious: A(Br) = (AB)r = Cr. For the second write D = AB − C; the check errs exactly when Dr = 0. Since D ≠ 0, it has a non-zero element dij. Look at the i-th coordinate:
$$(Dr)_i = d_{ij}\,r_j + \sum_{k \ne j} d_{ik}\,r_k.$$Imagine that all rk with k ≠ j have already been chosen, but rj not yet; this is legitimate, since the coordinates are independent. The sum on the right has become some number s, and (Dr)i = dijrj + s. Of the two values rj ∈ {0, 1} at most one can make this expression zero: if both did, then s = 0 and dij + s = 0, whence dij = 0. So (Dr)i = 0 with probability at most ½ whatever the other coordinates, and Dr = 0 all the more so.
The proof also shows what the scene shows: one error in cell (i, j) is caught exactly when rj = 1, that is, exactly half the time, and it shows up in row i. In the listing below the matrix A is the link graph of three of the site’s textbooks (202 pages, 1624 links), and C = A·A is the number of two-click paths between pages. The correct product passes all 20 checks; in a product with one cell increased by one, 64 separate checks catch the error 34 times.
If an algorithm with one-sided error errs with probability at most q, then k independent runs, where the answer “there is an error” counts if any one of them gave it, err with probability at most qk.
A false answer “all right” is possible only if all k runs missed the error. The runs are independent, so the probability of this is a product of k probabilities, each at most q. And when AB = C, no run says “there is an error”, so it cannot err that way at all.
Twenty Freivalds checks cost 60n² operations and err with probability less than 2−20 ≈ 10−6. With n = 256, as for the photograph, that is four million operations against 16.8 million for a full recomputation, and with n = 10,000, 6·109 against 1012.
Las Vegas: right, but who knows when
Now the opposite case. Let us solve a sudoku from a puzzle bank with simple local search: fill each 3×3 box with the missing digits at random, and swap two free cells of one box if this does not increase the number of repeats in rows and columns. When no repeats are left, the sudoku is solved and the answer checked: it is always right. But how many steps will it take?
The histogram looks unusual. Many runs solve the puzzle in a few thousand steps, but for the first sudoku about half did not solve it even in a hundred thousand: such a run is stuck, with two or three repeats left that no single swap can remove without creating new ones. To get out, the descent needs rare luck, and its running time follows a “heavy-tailed” distribution: the mean time is enormous, although the median is small. Carla Gomes, Bart Selman, Nuno Crato and Henry Kautz showed in 2000 that many searches for satisfiability and scheduling problems behave this way.
What can be done with such an algorithm? It can be turned into Monte Carlo: stop it by a timer and say “don’t know”. Or it can be restarted: cut off a run that has dragged on and start again from a random filling. In the listing below, of eight runs without restarts four did not solve the puzzle in 20,000 steps, and the rest managed in 1303–4031 steps. With a restart every 2000 steps the puzzle is solved in 5001 steps on average.
Let the running time of a Las Vegas algorithm be a random variable T ≥ 0 with mean E[T]. If it is stopped after t steps with the answer “don’t know”, this happens with probability P(T > t) ≤ E[T] / t. In particular, with t = 2E[T] we get a Monte Carlo algorithm with failure probability at most ½.
T is no less than the variable equal to t when T > t and zero otherwise. Means preserve the inequality: E[T] ≥ t·P(T > t). It remains to divide by t.
On the sudoku Markov’s inequality is useless: half the runs did not fit even into a hundred thousand steps, and the mean time there is over fifty thousand, so the bound for t = 2000 comes out above 25, while a probability is never above one. But the real distribution can be measured, and here is what follows from it.
Let the runs be independent and identically distributed, and let each be cut off after t steps if it has not finished, after which a new one begins. If P(T ≤ t) > 0, the mean number of steps until an answer equals
$$\frac{\mathbb E\,[\min(T,\ t)]}{P(T \le t)}.$$Write p = P(T ≤ t) and X for the mean sought. The first run costs min(T, t) steps. With probability p it ended with an answer; with probability 1 − p it was cut off, and everything starts again, with the past having no influence. So X = E[min(T, t)] + (1 − p)·X, whence p·X = E[min(T, t)]. To subtract (1 − p)X from both sides legitimately, X must be finite; it is, because the number of runs is geometrically distributed with mean 1/p, and each costs at most t steps, so X ≤ t/p.
This formula is the curve under the histogram in the scene, and its minimum is the best cutoff. The benchmark below: 400 runs on the first sudoku, and from them the cost of a solution with a restart every t steps. With t = 2000 a solution costs 5635 steps on average; with t = 500, 33,257, because a run almost never makes it in 500 steps; with t = 100,000, 103,101, because the stuck runs drag on to the end.
The best cutoff depends on the problem, and the distribution is not known in advance. Michael Luby, Alistair Sinclair and David Zuckerman found a universal rule in 1993: restart after 1, 1, 2, 1, 1, 2, 4, 1, 1, 2, 1, 1, 2, 4, 8, … units of time. This strategy on average runs no longer than the best fixed cutoff times a logarithmic factor, and they proved that this factor is unavoidable: a strategy that does not know the distribution loses to the best cutoff by a logarithmic factor on some distribution. The MiniSat satisfiability solver, for example, restarts its search by this sequence.
The reverse move, from Monte Carlo to Las Vegas, is possible when the answer can be checked. If one attempt costs a steps and gives a right answer with probability p, and a check costs v, then “try and check until it works” is Las Vegas with mean time (a + v)/p: the number of attempts is geometric with mean 1/p, as in the proof above. Our sudoku solver is built just so: a cut-off run is a failed attempt, and “zero repeats” is the check.
Two-sided error and voting
If an algorithm can err both ways, “at least one said no” no longer works: there are false answers either way. Then it is run k times and the majority answer taken. If each run is right with probability above ½, the majority is rarely wrong, and exponentially more rarely as k grows.
Let each of k independent runs err with probability p < ½. Then the majority errs, that is, at least half the answers are wrong, with probability at most (4p(1 − p))k/2.
The probability that exactly j runs erred is C(k, j)·pj(1 − p)k−j. Rewrite the factor: pj(1 − p)k−j = (p(1 − p))k/2·(p/(1 − p))j − k/2. For j ≥ k/2 the exponent of the last factor is non-negative and its base p/(1 − p) is less than one, so it is at most 1. So the sum over all j ≥ k/2 is at most (p(1 − p))k/2 times the sum of the binomial coefficients, which is at most 2k. We get 2k(p(1 − p))k/2 = (4p(1 − p))k/2.
With p = ⅓ the base is 4·⅓·⅔ = 8/9, and a hundred runs err with probability less than (8/9)50 ≈ 0.003, a thousand less than 10−25. The bound is crude (the Chernoff inequality gives a sharper one), but the main point is visible from it: each new run multiplies the error probability by a constant less than one. So for Monte Carlo algorithms it does not matter whether they err a third of the time or one time in a hundred: a few dozen repetitions make the error smaller than the probability of a hardware fault.
Fingerprints: a polynomial at a random point
The most powerful source of Monte Carlo algorithms is a simple property of polynomials: a non-zero polynomial cannot vanish too often.
Let f(x1, …, xm) be a non-zero polynomial of degree d with coefficients from a field, and S a finite set of elements of that field. If each xi is chosen independently and uniformly from S, then f(x1, …, xm) = 0 with probability at most d / |S|.
Induction on the number of variables. For m = 1 a non-zero polynomial of degree d has at most d roots: if r is a root, f(x) is divisible by x − r without remainder (the remainder on division by x − r equals f(r) = 0), the quotient has degree d − 1, and the other roots of f are its roots, because in a field a product is zero only if one of the factors is. So at most d of the |S| equally likely values make f zero.
The step. Single out the highest power k of the variable x1: f = x1k·g(x2, …, xm) + (terms where the power of x1 is less than k), where g ≠ 0 and the degree of g is at most d − k. First choose x2, …, xm. By the induction hypothesis g vanishes with probability at most (d − k)/|S|. If g ≠ 0, then with the other variables chosen f is a non-zero polynomial in x1 alone of degree k, and it vanishes with probability at most k/|S|. The probability that f vanishes is at most the sum of these two probabilities: d/|S|.
An application that does not even need the many-variable case is comparing large collections. Do two web server logs contain the same requests, possibly in a different order? Write each log as the polynomial (z − x1)(z − x2)…(z − xn), where xi are the request numbers. It does not depend on the order, while different collections have different polynomials: their difference is non-zero and of degree at most n. Evaluate both polynomials at one random point z modulo the prime P = 261 − 1: this is the log’s fingerprint, one number instead of a hundred thousand, computed in one pass without sorting. By the lemma the fingerprints of different logs coincide with probability at most n / |S|.
In the listing is the log of the Kennedy Space Center server for 13 July 1995, 128,932 requests. A sorted copy gives the same fingerprint. A copy in which one request was lost and another written twice gives a different one. And so does a copy in which two numbers are shifted towards each other by one, although its simple sum of numbers is the same. The point z in the listing is taken from 253 values, so the probability of a false match is at most 128,932 / 253 ≈ 1.4·10−11. The same lemma underlies many other checks: the equality of polynomials given by formulas, the existence of a perfect matching (through the determinant of the Tutte matrix) and polynomial string hashing.
What comes from what
Complexity theory is interested in these transformations in general form. The class ZPP is problems solved by a Las Vegas algorithm in expected polynomial time. RP is those solved by a Monte Carlo algorithm in polynomial time with one-sided error: on the answer “no” it is always right, and on “yes” it errs with probability at most ½. coRP is the same with the error the other way. The classes were introduced by John Gill in the 1970s.
ZPP = RP ∩ coRP: a problem is solved by Las Vegas in expected polynomial time if and only if it has Monte Carlo algorithms with one-sided error in both directions.
From ZPP to RP. Let Las Vegas run t(n) steps on average. Stop it after 2t(n) steps; if it has not finished, answer “no”. By Markov’s inequality it fails to finish with probability at most ½, so on a true answer “yes” we wrongly say “no” with probability at most ½, and we say “yes” only when Las Vegas said it, so rightly. This is RP. Answering “yes” on failure gives coRP.
From RP ∩ coRP to ZPP. There is an algorithm R whose “yes” is always right, and an algorithm R′ whose “no” is always right. Run them in turn until R says “yes” or R′ says “no”, and output that answer, which is certainly right. If the right answer is “yes”, R says “yes” with probability at least ½ in each round; if “no”, the same holds for R′. So the number of rounds is no worse than geometric with probability ½, and on average there are at most two of them, each polynomial.
With two-sided error one gets the class BPP, and the majority vote makes its error as small as you like. Whether BPP equals P, that is, whether any randomness can be replaced by deterministic computation without a big loss of time, is an open question; most specialists believe it can.
Code
Freivalds’ check on the site’s link graph. A is stored as link lists, so A(Ar) costs two passes over the links, while the claimed product C is an ordinary dense n × n table, as if it had been sent from outside.
Log fingerprints. The product of two 61-bit numbers is an ordinary integer in Python, while in Go it does not fit in 64 bits and is computed with the function bits.Mul64, which returns the high and low halves, and the remainder is taken with bits.Rem64.
And Las Vegas with restarts on a sudoku from the puzzle bank. The function descent is one attempt with a ceiling, with_restarts restarts until an attempt succeeds; the answer checks itself: zero repeats.
What it costs
| What | Time | Error |
|---|---|---|
| check AB = C by recomputing | O(n³) (O(n2.37…) with fast algorithms) | none |
| Freivalds, k times | O(kn²) | at most 2−k, only “right” when wrong |
| fingerprints of two collections of n elements | O(n) modular multiplications | at most n / |S| |
| comparison by sorting | O(n log n) and memory for copies | none |
| Las Vegas stopped at c·E[T] | c·E[T] | failure at most 1/c of the time |
| Las Vegas with a restart every t | E[min(T, t)] / P(T ≤ t) on average | none |
| Monte Carlo with a check, success probability p | (attempt + check) / p on average | none |
| a vote of k runs | k runs | at most (4p(1 − p))k/2 |
The general rule reads off the table: repetition costs linearly, while the error falls exponentially. That is why a Monte Carlo algorithm with error probability ½ is practically as good as one with 10−6: the difference is twenty repetitions.
Where it works
- Checking others’ computations. Manuel Blum and Sampath Kannan proposed “program checkers” in 1989: instead of proving a program correct, check each of its answers with a fast randomised test. Freivalds’ check is the classic example.
- Comparing and synchronising data. Fingerprints modulo a prime let you compare files and collections without sending them; polynomial string hashes and the Rabin–Karp algorithm rest on the same idea.
- Restarts in solvers. Satisfiability solvers and search programs for timetables restart their search to cut off heavy tails; MiniSat does this by the Luby sequence by default.
- Classics of the catalogue: quicksort and quickselect with a random pivot are Las Vegas; the Miller–Rabin test, Karger’s minimum cut and the Bloom filter are Monte Carlo with one-sided error.
History
Randomised algorithms for “pure” problems rather than simulation appeared in the second half of the 1970s. In 1976–1977 Michael Rabin, Robert Solovay and Volker Strassen proposed probabilistic primality tests, in 1977 Rūsiņš Freivalds the check of matrix multiplication, and in 1978–1980 Richard DeMillo and Richard Lipton, Richard Zippel and Jacob Schwartz independently proved the lemma on the roots of a polynomial that now bears the names of the last two. John Gill at the same time described the classes of problems solved by randomised algorithms, and in 1979 László Babai, working on graph isomorphism testing, called algorithms that are never wrong but may run long Las Vegas algorithms.
Problems
In the first problem you have to check a matrix product too large to recompute. In the second, from a sample of the running times of a Las Vegas algorithm, choose the best cutoff for restarts; one of the tests is real runs of the sudoku solver from the article.
Checking the server When to start overData: Photographs by the site's author The link graph of the legost.in textbooks: the basic course, the mathematics course, the algorithms catalogue Kennedy Space Center web server log, 13 July 1995: NASA-HTTP, Internet Traffic Archive, freely redistributable Puzzles: the Sudoku Exchange puzzle bank, public domain; difficulty ratings by Sukaku Explainer