Monte Carlo Algorithms, in depth: one- and two-sided error, amplification, Freivalds' check, Schwartz-Zippel and estimation with confidence intervals

Monte Carlo versus Las Vegas, and the two kinds of error

Randomized algorithms split along one axis: what is random. A Las Vegas algorithm is always correct and its running time is random; randomized quicksort is the textbook case. A Monte Carlo algorithm has a bounded running time and its correctness is random. The probability is over the algorithm's own coin flips and holds for every input, which is the crucial difference from average-case analysis: no adversary can choose an input that makes the algorithm fail more often, as long as the adversary cannot see or predict the coins.

Monte Carlo decision algorithms come in two flavours. A one-sided error algorithm is never wrong on one of the two answers. Freivalds' check never says "AB ≠ C" when the product is right, and Miller-Rabin never calls a prime composite; when they say NO they hold a witness that proves it. Only the YES answer can be wrong. A two-sided error algorithm can be wrong either way, with probability at most some ε below 1/2 on every input. Complexity theory names the classes: RP and co-RP for one-sided polynomial-time tests, BPP for two-sided ones with ε ≤ 1/3.

Amplification: buying down the error

The error bound of a single run is rarely small enough to use as it stands. Amplification fixes that by repeating with fresh, independent randomness.

One-sided error. Run K times and answer NO if any run finds a witness. A wrong final YES needs every run to miss, so the error is at most ε^K. With ε = 1/2, 20 runs give about one in a million and 64 runs give 2^-64. Cost grows linearly in K while error falls exponentially, which is why one-sided tests are so pleasant.

Two-sided error. Any single run can lie, so you take a majority vote over K runs. If each run is right with probability p = 1/2 + δ, the Hoeffding bound says the majority is wrong with probability at most exp(-2Kδ²). For p = 2/3 that is exp(-K/18): 101 runs bound the error by about 0.004. The bound is loose. In our simulation of a 2/3-correct test the measured majority error was 0.333, 0.120, 0.027 and 0.0003 at K = 1, 11, 31 and 101, so the real error falls faster than the bound promises, but the bound is what you can prove.

The difference matters when δ is small. A test that is right 51% of the time needs about 46,000 votes before the bound reaches 10^-4, which is why practical two-sided algorithms work hard to make each run's advantage large.

Freivalds' matrix-product check

Freivalds' algorithm (1977) is the cleanest example in the class. Given n×n matrices A, B and C, decide whether AB = C. Recomputing the product costs O(n³) naively. Instead, pick a random vector r with entries in {0, 1} and compare A(Br) with Cr. That is three matrix-vector products, O(n²) work in total.

Why it works: let D = AB - C. If D = 0 the two sides agree for every r, so a correct product is never rejected. If D ≠ 0, some entry D[i][j] is non-zero. Fix all coordinates of r except r_j. Then (Dr)_i = D[i][j]·r_j + (everything else), and at most one of the two choices of r_j can make that zero. So each trial catches the error with probability at least 1/2, and a NO is backed by the vector r that exposed it.

import random

def matvec(M, v):
    return [sum(row[k] * v[k] for k in range(len(v))) for row in M]

def freivalds(A, B, C, trials=20, rng=random):
    """False is always right; True is wrong with probability <= 2**-trials."""
    n = len(C[0])
    for _ in range(trials):
        r = [rng.randrange(2) for _ in range(n)]
        if matvec(A, matvec(B, r)) != matvec(C, r):
            return False        # r is a witness: AB != C for certain
    return True                 # probably AB == C

Worked example and measured error

Worked example. Take A = [[1, 2], [3, 4]] and B = [[5, 6], [7, 8]]. The true product is [[19, 22], [43, 50]]. Suppose a buggy kernel returned C with the bottom right entry 51. Then D = AB - C is zero everywhere except D[1][1] = -1.

Draw r = (1, 0). Br = (5, 7), A(Br) = (19, 43), and Cr = (19, 43). They match: this trial missed, because the error lives in column 1 and r_1 = 0. Draw r = (0, 1). Br = (6, 8), A(Br) = (22, 50), Cr = (22, 51). Mismatch in row 1, so the product is wrong, and r = (0, 1) proves it. Of the four possible vectors, exactly the two with r_1 = 1 catch the error: the 1/2 bound is tight for a single wrong entry.

We checked this at scale. On a 60×60 integer product with one corrupted entry, the measured false-accept rate of a single trial was 0.503 over 20,000 trials, matching the analysis. Twenty trials cost 60n² multiply-adds against n³ for a naive product, so the check breaks even at n = 60 and costs about 1/67 of the product at n = 4,000.

The family: Schwartz-Zippel and friends

The same pattern of a cheap random probe plus an algebraic reason it rarely misses shows up across the field.

The Schwartz-Zippel lemma drives most of these. If P is a non-zero polynomial of total degree d and you evaluate it at a point drawn uniformly from S^n, for a finite set S, then P evaluates to zero with probability at most d/|S|. Working modulo a 61-bit prime with d = 1,000 gives error below 10^-15 per evaluation. That one lemma lets you test whether two arithmetic circuits compute the same polynomial, whether a bipartite graph has a perfect matching (via the Edmonds matrix determinant), and whether two strings are equal by comparing their polynomial fingerprints, which is the idea behind rolling hashes.

Monte Carlo estimation

The second meaning of "Monte Carlo" is numerical: estimate a quantity by averaging random samples. The guarantee is statistical rather than a yes/no error bound. If X has mean μ and variance σ², the average of N independent samples has standard error σ/√N, and for large N the central limit theorem gives an approximate 95% interval of mean ± 1.96·s/√N.

The classic demonstration estimates π: sample points in the unit square and count the fraction inside the quarter circle, times four. With a fixed seed we measured 3.152 ± 0.101 at N = 1,000 and 3.1422 ± 0.0102 at N = 100,000. A hundred times more samples bought ten times the precision. That square-root law is the defining fact of Monte Carlo estimation, and it holds in any dimension, which is why Monte Carlo beats grid quadrature for integrals over many variables, such as option pricing, light transport in rendering, and the expectations inside a policy-gradient update.

import math

def mc_mean(f, sampler, n, rng):
    """Estimate E[f(X)] with a 95% normal-approximation half-width."""
    s = s2 = 0.0
    for _ in range(n):
        y = f(sampler(rng))
        s += y
        s2 += y * y
    mean = s / n
    var = (s2 - n * mean * mean) / (n - 1)
    return mean, 1.96 * math.sqrt(var / n)

Because the error falls only as 1/√N, the useful lever is usually σ rather than N. Variance-reduction techniques such as control variates, antithetic sampling, stratification and importance sampling can cut σ by orders of magnitude. They are often worth more than a bigger cluster.

Converting between Monte Carlo and Las Vegas

Monte Carlo and Las Vegas convert into each other under the right conditions. If you can verify an answer cheaply, a Monte Carlo search becomes Las Vegas: repeat until the verifier accepts. The answer is then always right and only the running time is random. If the success probability per run is p, the expected number of runs is 1/p. In the other direction, any Las Vegas algorithm with expected time T becomes Monte Carlo by stopping it after 2T steps and answering arbitrarily. Markov's inequality bounds the failure probability by 1/2, and repetition shrinks it further.

This is a practical design choice. A factoring routine that checks its factors by multiplying has turned an uncertain step into a verified one. When a verifier is cheap, prefer it: a guarantee you can check beats one you have to trust.

Operational guidance

Choosing K is an engineering decision, so make it explicitly. Pick a target failure rate per call, multiply by the number of calls per year, and size K so the expected number of wrong answers across the fleet stays well below one. For one-sided ε = 1/2 tests, K = 64 is a common, easily defended choice.

• Seed from a real source. Use random.SystemRandom or secrets when an adversary chooses the inputs. Use a seeded random.Random(seed) when you need reproducible runs, and log the seed.
• Keep trials independent. Each trial needs fresh randomness. Reusing one r for K trials gives you one trial K times.
• Log the witness. When a one-sided test says NO, store the witness (r, the failing base, the cut). It turns a mysterious rejection into a reproducible bug report.
• Report intervals, not points. For estimators, ship the half-width with the estimate, and stop sampling when it falls below the tolerance instead of using a fixed N.
• Use exact arithmetic or a prime field for algebraic tests. Floating-point Freivalds on real matrices needs a tolerance, and a badly chosen one creates both kinds of error.

Failure modes

The error bound is rarely what breaks. These are the failures we see in practice.

• Predictable randomness. The bound assumes the input does not depend on the coins. Fix Miller-Rabin's bases and someone can construct composites that pass them; Arnault published such numbers in 1995. Hash-flooding attacks on fixed-seed hash tables are the same failure. If inputs can be adversarial, the coins must be secret.
• Correlated parallel streams. Workers seeded with the same value, or with consecutive seeds on a weak generator, produce correlated samples. The confidence interval then reports more precision than you have. Use a generator designed for independent streams, such as NumPy's SeedSequence.spawn.
• Heavy tails. If f(X) has infinite or huge variance, the sample mean converges slowly and the normal interval is fiction. Plot the samples and watch for a few values dominating the sum.
• Two-sided tests near 1/2. A small advantage δ needs on the order of 1/δ² votes. Budgets sized for δ = 1/6 collapse when the real advantage on some input class is 0.01.
• Overflow in the cheap check. Integer products in fixed-width types wrap. Do Freivalds modulo a prime, or in Python's arbitrary-precision integers.

Trade-offs

Monte Carlo trades certainty for speed and simplicity. Take the trade when a wrong answer costs about as much as a hardware fault, and the bound can be driven far below that. Avoid it when you need a proof for a third party, as in certified primality for a standards body or a legal audit. Avoid it too when the randomness cannot be kept secret from someone who benefits from a wrong answer. Prefer Las Vegas when a cheap verifier exists, and a deterministic algorithm when one is nearly as fast. The deterministic AKS primality test exists, yet libraries use Miller-Rabin because it is far faster.

What to do next

1. Implement Freivalds modulo a large prime, then check it against a GPU or BLAS matrix product you already trust. It is a cheap invariant for numerical kernels.
2. Write down the per-call failure budget for any probabilistic check you run, and derive K from it rather than copying a constant.
3. Read the Miller-Rabin article for a one-sided test with witnesses, and the Karger min-cut article for a Monte Carlo optimisation algorithm amplified by repetition.
4. Use polynomial hashing as a Schwartz-Zippel fingerprint, and compute its collision probability for your string lengths.
5. For estimation, add a stopping rule driven by the confidence half-width, then try one variance-reduction trick and measure the change in σ.
6. See Monte Carlo tree search for random rollouts used to guide decisions rather than estimate a single number.