Randomized & Approximation Algorithms

When we say an algorithm is randomized, we mean it makes random choices during execution — typically by sampling uniformly from some set or flipping coins. But there are two fundamentally different ways to use randomness, and confusing them is a common source of errors.

Las Vegas algorithms always produce the correct answer. The randomness only affects the running time. You cannot get a wrong answer from a Las Vegas algorithm; you can only be unlucky about how long it takes. Randomized quicksort is the canonical Las Vegas algorithm: the pivot choices are random, so the running time is a random variable, but the sorted output is always correct. When we analyze Las Vegas algorithms, we ask: what is the expected running time?

Monte Carlo algorithms run in fixed (or bounded) time but may produce an incorrect answer with some small probability. There is no randomness in the time budget — only in the correctness of the output. Miller-Rabin primality testing is the classic example: it runs in O(k log²n) time for k rounds, but declares a composite number prime with probability at most 4^{-k}. The guarantee is probabilistic, not absolute. When we analyze Monte Carlo algorithms, we ask: what is the probability of error?

The trade-off is fundamental. Las Vegas algorithms give you certainty of correctness at the cost of uncertain runtime. Monte Carlo algorithms give you a time budget at the cost of uncertain correctness. In practice, you can often convert one into the other: run a Monte Carlo algorithm until it agrees with itself k times (reducing error exponentially), or cut a Las Vegas algorithm short with a timeout (accepting a small failure probability).

A subtler taxonomy: Monte Carlo algorithms are one-sided if they can err in only one direction. A primality test that never calls a prime composite, but might call a composite prime, has one-sided error. Algorithms with two-sided error can be wrong in both directions. One-sided error is strictly weaker and often easier to handle: you can repeat and take the majority answer.

import random

# Las Vegas: randomized selection (always correct, random time)
def las_vegas_find(arr, k):
    """Find k-th smallest element. Always correct, expected O(n) time."""
    arr = arr[:]
    while True:
        pivot = random.choice(arr)
        lo = [x for x in arr if x < pivot]
        eq = [x for x in arr if x == pivot]
        hi = [x for x in arr if x > pivot]
        if k < len(lo):
            arr = lo
        elif k < len(lo) + len(eq):
            return pivot      # correct answer, guaranteed
        else:
            k -= len(lo) + len(eq)
            arr = hi

# Monte Carlo: estimate pi (fast, but approximate / probabilistically correct)
def monte_carlo_pi(n_samples):
    """Returns estimate of pi. More samples → smaller error."""
    inside = sum(
        1 for _ in range(n_samples)
        if random.random()**2 + random.random()**2 <= 1.0
    )
    return 4 * inside / n_samples

print(las_vegas_find([3,1,4,1,5,9,2,6], 4))  # always the 4th-smallest
print(monte_carlo_pi(1_000_000))              # ≈ 3.14159, not exact

Deterministic quicksort has a notorious Achilles' heel: if you always pick the first or last element as the pivot, a sorted or reverse-sorted array triggers O(n²) behavior. Adversarial inputs exist for any fixed pivot rule. The fix is elegant: choose the pivot uniformly at random at each recursive step. This single change guarantees expected O(n log n) time on every input, with no adversarial input possible — because the adversary doesn't know which random choices you will make.

The key analysis tool is Linearity of Expectation. Rather than tracking the entire recursion tree, we count comparisons directly.

Let the elements of the input be x₁ < x₂ < ... < xₙ (relabeled in sorted order). Define the indicator random variable:

X_{ij} = 1 if x_i and x_j are ever compared, 0 otherwise.

The total number of comparisons is C = Σᵢ<ⱼ X_{ij}. By Linearity of Expectation:

E[C] = Σᵢ<ⱼ Pr[x_i and x_j are compared].

The crucial observation: x_i and x_j are compared if and only if one of them is chosen as a pivot before any element in the set {x_i, x_{i+1}, ..., x_j} is chosen as a pivot. (Once any intermediate element becomes a pivot, x_i and x_j land in different subtrees and never meet.)

Since each of the (j-i+1) elements in {x_i, ..., x_j} is equally likely to be chosen first (uniform random pivots), the probability that x_i or x_j is chosen first is exactly 2/(j-i+1).

Therefore: E[C] = Σᵢ<ⱼ 2/(j-i+1) = 2 Σₖ₌₁ⁿ⁻¹ (n-k) · 1/(k+1) ≈ 2n ln n = O(n log n).

This is a beautiful example of the indicator random variable technique combined with Linearity of Expectation — two tools you will use again and again in randomized algorithm analysis.

For the failure probability: by Markov's inequality (or a more careful analysis using Chernoff bounds), the probability that randomized quicksort takes more than 10n ln n comparisons is less than 1/n. So the algorithm is not just fast in expectation — it is concentrated around its expectation.

import random

def randomized_quicksort(arr):
    if len(arr) <= 1:
        return arr
    pivot = random.choice(arr)   # uniform random pivot — the key move
    lo = [x for x in arr if x < pivot]
    eq = [x for x in arr if x == pivot]
    hi = [x for x in arr if x > pivot]
    return randomized_quicksort(lo) + eq + randomized_quicksort(hi)

# Demonstration: sorted input — the adversarial case for deterministic quicksort
n = 1000
sorted_input = list(range(n))
# Deterministic (first-element pivot) would be O(n^2) here.
# Randomized is O(n log n) in expectation — safe!
result = randomized_quicksort(sorted_input)
assert result == sorted_input   # always correct (Las Vegas)

A hash table promises expected O(1) insert/lookup. But the word expected hides a critical question: expected over what? If our hash function is fixed, an adversary can craft inputs that all map to the same bucket, forcing O(n) operations. The solution is universal hashing: randomize the hash function itself.

Definition (Universal Hash Family). A family ℋ of functions h: U → {0,1,...,m-1} is universal if for any two distinct keys x ≠ y: Pr_{h ∈ ℋ} [h(x) = h(y)] ≤ 1/m.

The simplest universal family for integer keys is the multiply-mod-prime family. Choose a prime p ≥ |U| and let: h_{a,b}(x) = ((ax + b) mod p) mod m, where a ∈ {1,...,p-1}, b ∈ {0,...,p-1}.

Why does universality give O(1) expected time? Fix any key x. For any other key y, the probability of collision is ≤ 1/m. By Linearity of Expectation, the expected number of keys that collide with x is at most n/m. If we keep the load factor n/m = O(1) (by resizing when n exceeds c·m for some constant c), each lookup examines O(1) keys in expectation.

Perfect Hashing (FKS scheme): Can we achieve worst-case O(1) lookups for a static set S of n keys? Yes! The Fredman-Komlós-Szemerédi (FKS) scheme works in two levels:

1. Map n keys to n buckets using a universal hash function. Some buckets get multiple keys.
2. For each bucket i with kᵢ keys, build a second-level hash table of size kᵢ². A random hash function in a universal family has no collisions at this second level with probability ≥ 1/2 (by a birthday-paradox argument: the expected number of collisions is at most kᵢ(kᵢ-1)/(2kᵢ²) = 1/2, so by Markov's inequality the probability of any collision is ≤ 1/2).

The total space is O(n) expected (since Σ kᵢ² = O(n) for a good first-level function). Lookups are always O(1) worst-case.

Cuckoo hashing is the modern favorite: two tables, two hash functions. Every key sits in one of two possible slots. Insert by displacing the current occupant ("cuckoo"-style) until everyone finds a home. Lookup is always O(1) worst-case (check two slots). Expected O(1) amortized insert.

import random

class UniversalHashTable:
    """Chaining hash table with multiply-mod-prime universal family."""
    def __init__(self, m=16):
        self.m = m
        self.p = 1_000_000_007  # prime larger than key universe
        self.a = random.randint(1, self.p - 1)
        self.b = random.randint(0, self.p - 1)
        self.buckets = [[] for _ in range(m)]
        self.n = 0

    def _hash(self, key):
        return ((self.a * hash(key) + self.b) % self.p) % self.m

    def insert(self, key, val):
        i = self._hash(key)
        for j, (k, _) in enumerate(self.buckets[i]):
            if k == key:
                self.buckets[i][j] = (key, val)
                return
        self.buckets[i].append((key, val))
        self.n += 1

    def lookup(self, key):
        for k, v in self.buckets[self._hash(key)]:
            if k == key:
                return v
        return None

t = UniversalHashTable()
t.insert("alice", 42)
t.insert("bob", 99)
print(t.lookup("alice"))  # 42

For NP-hard optimization problems, we cannot (under P ≠ NP) hope for a polynomial-time exact algorithm. But we can still say something rigorous. An approximation algorithm is a polynomial-time algorithm that produces a provably good-enough solution.

Definition (Approximation Ratio). For a minimization problem, an algorithm A is an α-approximation (where α ≥ 1) if for every instance I: A(I) ≤ α · OPT(I), where OPT(I) is the optimal solution value and A(I) is the cost of A's output. For maximization, the inequality flips: A(I) ≥ (1/α) · OPT(I), or equivalently we require A(I)/OPT(I) ≥ 1 - ε for approximation schemes.

The key challenge: how do you bound OPT(I) when you can't compute it? This is the central intellectual puzzle of approximation algorithm design. The standard technique is to find a lower bound on OPT (for minimization) that you can compute efficiently — typically a relaxation of the problem. If you can show that your solution costs at most α times this lower bound, and the lower bound is at most OPT, you get a provable α-approximation.

PTAS and FPTAS. An algorithm is a Polynomial-Time Approximation Scheme (PTAS) if for every fixed ε > 0 it achieves a (1+ε)-approximation in polynomial time (though the exponent may depend on 1/ε). A Fully Polynomial-Time Approximation Scheme (FPTAS) runs in time polynomial in both n and 1/ε — strictly stronger and rarer.

Inapproximability. Not all problems can be approximated well. The Inapproximability results — proved using the PCP Theorem and reductions — show, for example, that MAX-3SAT cannot be approximated beyond a factor of 7/8 in polynomial time (unless P=NP), and that Clique cannot be approximated within a factor of n^{1-ε} for any ε > 0. These results explain why certain approximation ratios are tight.

# Pseudocode: generic 2-approximation for weighted set cover
# (greedy: always pick the set with the best cost/coverage ratio)
def greedy_set_cover(universe, sets_with_costs):
    covered = set()
    chosen = []
    total_cost = 0
    remaining = set(universe)
    while remaining:
        # pick set maximizing |new elements covered| / cost
        best = max(
            sets_with_costs,
            key=lambda sc: len(sc[0] & remaining) / sc[1]
                          if len(sc[0] & remaining) > 0 else -1
        )
        chosen.append(best)
        remaining -= best[0]
        total_cost += best[1]
    return chosen, total_cost
# This achieves an H_n-approximation (H_n = 1 + 1/2 + ... + 1/n, the harmonic sum).

Vertex Cover is: given a graph G = (V, E), find a minimum-size subset S ⊆ V such that every edge has at least one endpoint in S. This problem is NP-hard in general, but it admits an elegant, provably correct 2-approximation.

The algorithm (Maximal Matching):

1. Find any maximal matching M (a set of edges no two of which share a vertex, and to which no more edges can be added).
2. Return all endpoints of edges in M as the vertex cover.

Why is the output always a valid cover? Suppose some edge (u,v) is not covered — neither u nor v is in the output. Then (u,v) ∉ M, and since M is maximal, some edge in M shares an endpoint with (u,v). But then that endpoint is in the cover, contradicting our assumption. Contradiction. So every edge is covered. ✓

Why is the approximation ratio exactly 2? Let M be the matching found. Our algorithm's output has size |output| = 2|M|. Now, in any vertex cover (including OPT), at least one endpoint of every matching edge must be included. Since matching edges are vertex-disjoint, OPT must contain at least |M| vertices just to cover the matching edges. Therefore: OPT ≥ |M|, and |output| = 2|M| ≤ 2·OPT.

This bound is tight: there exist graphs where the optimal vertex cover is exactly |M|, and the algorithm returns 2|M|. This is a beautiful, crisp proof: the lower bound on OPT (the size of the matching) is itself a subproblem easily solvable in polynomial time.

Can we do better than 2? This is a major open question! The Unique Games Conjecture (Khot 2002) implies that 2 is optimal — that no polynomial-time algorithm can achieve a (2-ε)-approximation for any ε > 0. This is one of the central open problems in approximation complexity.

Note: a maximal matching is easy to find greedily in O(E) time. (A maximum matching requires more sophisticated algorithms; we only need a maximal one here — any greedy traversal of edges works.)

def two_approx_vertex_cover(n, edges):
    """
    n: number of vertices (0-indexed)
    edges: list of (u, v) pairs
    Returns a vertex cover of size <= 2 * OPT.
    """
    in_matching = [False] * n
    cover = set()

    for u, v in edges:
        if not in_matching[u] and not in_matching[v]:
            # Add this edge to the maximal matching
            in_matching[u] = True
            in_matching[v] = True
            cover.add(u)
            cover.add(v)

    return cover

# Example: path graph 0-1-2-3
edges = [(0,1),(1,2),(2,3)]
print(two_approx_vertex_cover(4, edges))
# Could output {0,1,2,3} (size 4) while OPT={1,2} (size 2).
# But also might output {1,2} exactly. Worst case is 2*OPT.

The Travelling Salesman Problem (TSP) asks: given n cities with pairwise distances, find the shortest tour that visits every city exactly once and returns to the start. In the general case, TSP cannot be approximated within any constant factor in polynomial time (unless P=NP) — there exists a polynomial-time reduction showing that any constant-factor approximation would solve 3SAT.

But for Metric TSP — where distances satisfy the triangle inequality: d(u,w) ≤ d(u,v) + d(v,w) — we can do much better.

2-Approximation for Metric TSP (Double-Tree algorithm):

1. Compute a Minimum Spanning Tree T of the complete graph on the n cities. This costs OPT_MST ≤ OPT_TSP, because removing any edge from an optimal tour gives a spanning tree.
2. Perform a DFS traversal of T, listing cities in the order first visited. This gives a walk of length 2·OPT_MST ≤ 2·OPT_TSP (each edge traversed twice).
3. Shortcut the walk: whenever we revisit a city, skip it. The triangle inequality guarantees shortcutting never increases distance.

Why 2·OPT? OPT_MST ≤ OPT_TSP (step 1), the DFS walk has length ≤ 2·OPT_MST (step 2), and shortcutting only helps (step 3). So the final tour costs ≤ 2·OPT_TSP. ✓

Christofides Algorithm (3/2-approximation, 1976): Instead of the DFS walk, observe that the DFS traversal has a problem: the MST has odd-degree vertices, and traversing it requires some edges twice. Christofides' key insight:

1. Find MST T.
2. Let O = set of odd-degree vertices in T. By handshaking lemma, |O| is even.
3. Find a minimum-weight perfect matching M on the vertices in O. Since there are even many of them, this is possible in polynomial time.
4. Combine T ∪ M to get an Eulerian multigraph (every vertex now has even degree).
5. Find an Eulerian circuit, then shortcut.

Why 3/2? OPT_MST ≤ OPT_TSP. The matching M costs at most OPT_TSP/2 (consider the optimal TSP tour restricted to the odd-degree vertices — it gives two matchings, the cheaper of which costs ≤ OPT_TSP/2). Total cost ≤ OPT_MST + OPT_TSP/2 ≤ 3/2 · OPT_TSP. ✓

Recent breakthrough: In 2020, Karlin, Klein, and Oveis Gharan improved this to (3/2 - ε) for a small fixed ε > 0, breaking the 47-year-old Christofides barrier. As of 2025, the best known ratio is (3/2 - 1/600+), and it remains an open problem whether a (1+ε)-PTAS exists for Metric TSP.

# 2-approximation for Metric TSP: MST + DFS shortcut
import heapq

def prim_mst(dist):
    """Prim's algorithm. dist[i][j] = distance. Returns list of MST edges."""
    n = len(dist)
    in_mst = [False] * n
    key = [float('inf')] * n
    parent = [-1] * n
    key[0] = 0
    pq = [(0, 0)]
    mst_adj = [[] for _ in range(n)]
    while pq:
        d, u = heapq.heappop(pq)
        if in_mst[u]: continue
        in_mst[u] = True
        if parent[u] != -1:
            mst_adj[u].append(parent[u])
            mst_adj[parent[u]].append(u)
        for v in range(n):
            if not in_mst[v] and dist[u][v] < key[v]:
                key[v] = dist[u][v]
                parent[v] = u
                heapq.heappush(pq, (key[v], v))
    return mst_adj

def tsp_2approx(dist):
    n = len(dist)
    adj = prim_mst(dist)
    # DFS preorder traversal
    visited = [False] * n
    tour = []
    def dfs(u):
        visited[u] = True
        tour.append(u)
        for v in adj[u]:
            if not visited[v]:
                dfs(v)
    dfs(0)
    tour.append(tour[0])  # return to start
    cost = sum(dist[tour[i]][tour[i+1]] for i in range(n))
    return tour, cost

Two of the most practically impactful randomized data structures are built on hashing: Bloom filters and MinHash sketches.

Bloom Filters answer the question: "Have I seen this element before?" with bounded false-positive rate and zero false negatives, using only a fraction of the memory of a hash set. A Bloom filter for a universe U is a bit array B of m bits, initialized to 0, and k independent hash functions h₁,...,hₖ: U → {0,...,m-1}.

• Insert(x): Set B[hᵢ(x)] = 1 for all i = 1,...,k.
• Query(x): Return YES if B[hᵢ(x)] = 1 for all i = 1,...,k; else NO.

False negative rate: 0. If x was inserted, all its bits are set; every query returns YES. False positive rate (probability of wrongly saying YES for an uninserted x): After inserting n elements, the probability that any specific bit remains 0 is approximately (1 - 1/m)^{kn} ≈ e^{-kn/m}. The false positive probability is approximately (1 - e^{-kn/m})^k. Choosing k = (m/n)ln2 minimizes this, giving false positive probability ≈ (1/2)^k ≈ (0.6185)^{m/n}. With m = 10n bits, the false positive rate is below 1%.

Bloom filters are used in databases (avoid disk reads for missing keys), browsers (Chrome's Safe Browsing list), and distributed caches.

MinHash (Jaccard Similarity Sketching) solves: given two large sets A and B, estimate |A ∩ B|/|A ∪ B| (Jaccard similarity) without materializing the union. Define a random min-hash function h: pick a random hash h: U → [0,1] and let MinHash(S) = min_{x ∈ S} h(x). Key property:

Pr[MinHash(A) = MinHash(B)] = |A ∩ B| / |A ∪ B| = J(A,B).

By repeating with k independent hash functions and averaging the indicators, we estimate J(A,B) to within ε with O(1/ε²) hashes. MinHash powers document similarity detection (plagiarism detection, web crawl deduplication) and recommendation systems.

import hashlib, math

class BloomFilter:
    def __init__(self, capacity, fp_rate=0.01):
        self.m = math.ceil(-capacity * math.log(fp_rate) / (math.log(2)**2))
        self.k = max(1, round((self.m / capacity) * math.log(2)))
        self.bits = bytearray(math.ceil(self.m / 8))

    def _hashes(self, item):
        h = int(hashlib.md5(str(item).encode()).hexdigest(), 16)
        h2 = int(hashlib.sha1(str(item).encode()).hexdigest(), 16)
        return [(h + i * h2) % self.m for i in range(self.k)]

    def add(self, item):
        for bit in self._hashes(item):
            self.bits[bit // 8] |= (1 << (bit % 8))

    def __contains__(self, item):
        return all(
            self.bits[bit // 8] & (1 << (bit % 8))
            for bit in self._hashes(item)
        )

bf = BloomFilter(capacity=1000, fp_rate=0.01)
bf.add("hello")
print("hello" in bf)   # True  (never wrong)
print("world" in bf)   # False (unless false positive, < 1%)

Expected-value analysis tells you on average how a randomized algorithm behaves. But for real applications, you need to know: how often does it deviate badly from the average? Chernoff bounds give exponentially tight concentration inequalities for sums of independent random variables — the main tool for proving that randomized algorithms have good performance with high probability (w.h.p.), meaning probability at least 1 - 1/n^c for an arbitrarily large constant c.

Setup. Let X₁, ..., Xₙ be independent Bernoulli random variables (each is 0 or 1). Let X = Σ Xᵢ and μ = E[X].

Chernoff Upper Tail: For any δ > 0: Pr[X ≥ (1+δ)μ] ≤ (eᵟ / (1+δ)^{(1+δ)})^μ ≤ e^{-μδ²/3} for δ ≤ 1.

Chernoff Lower Tail: For 0 < δ < 1: Pr[X ≤ (1-δ)μ] ≤ e^{-μδ²/2}.

How to read this. The probability that X deviates from its expectation by more than a δ-fraction is exponentially small in μδ². This is spectacularly strong. For randomized quicksort, it tells us the depth of the recursion tree is O(log n) with high probability. For hashing, it tells us no bucket gets more than O(log n / log log n) elements. For random sampling, it tells us a sample of size O(1/ε² · log(1/δ)) gives an ε-accurate estimate with probability ≥ 1-δ.

Method of bounded differences (McDiarmid's inequality) extends Chernoff to non-sum random variables: if changing any single input affects the output by at most cᵢ, then Pr[|f - E[f]| ≥ t] ≤ 2 exp(-2t² / Σcᵢ²). This is the workhorse for analyzing more complex randomized algorithms.

Chernoff bounds are proved via the moment generating function (MGF) trick: for any t > 0, Pr[X ≥ a] = Pr[e^{tX} ≥ e^{ta}] ≤ E[e^{tX}] / e^{ta} by Markov's inequality, then expand the MGF using independence.

import random, math

def chernoff_demo(n=1000, p=0.5, delta=0.1, trials=100_000):
    """Empirically verify Chernoff bound for coin flips."""
    mu = n * p
    threshold = (1 + delta) * mu
    # Chernoff bound: Pr[X >= (1+delta)*mu] <= exp(-mu*delta^2/3)
    chernoff_bound = math.exp(-mu * delta**2 / 3)

    exceed = sum(
        1 for _ in range(trials)
        if sum(random.random() < p for _ in range(n)) >= threshold
    ) / trials

    print(f"Empirical Pr[X >= {threshold:.0f}]: {exceed:.6f}")
    print(f"Chernoff upper bound:           {chernoff_bound:.6f}")
    # The empirical value should be <= the bound

chernoff_demo()