Random Spanning Tree, in depth: why MST and DFS are biased, Aldous-Broder, the Wilson loop-erased walk, effective resistance and testing a sampler

What uniform means and how many trees there are

Uniform means probability 1/T for each tree, where T is the number of spanning trees. Kirchhoff's matrix-tree theorem gives T exactly: build the Laplacian L = D - A (degrees on the diagonal, -1 for each edge), delete any one row and the matching column, and take the determinant. The complete graph on 4 vertices has 16 spanning trees (Cayley's formula n to the power n-2), a 3 by 3 grid has 192, and the diamond graph used below, two triangles sharing an edge, has 8. The numbers in this article are computed by an exact rational determinant when the page is built, not copied from a table.

T grows exponentially with graph size, so you can never list the trees and pick one for a real graph. The algorithms below sample a uniform tree without ever counting, in time related to how quickly a random walk explores the graph. The counting tools are still useful as test oracles on small graphs, and the same Laplacian structure appears in PageRank and in the spanning-tree counts for cactus graphs.

Why the obvious shortcuts are biased

Two shortcuts are tempting. The first: give every edge an independent random weight and take the minimum spanning tree with Kruskal. Since random weights only matter through their order, this is Kruskal on a uniformly random edge order, so it can be enumerated exactly over all 120 orders of the diamond's five edges. The result is not uniform: trees take probability 7/60 or 2/15 instead of 1/8, and the shared edge appears with probability 8/15 instead of 1/2. The cycle property explains it: an edge is left out exactly when it comes last on some cycle through it. The shared edge lies on two triangles, so it is left out with probability 1/3 + 1/3 - 1/5 = 7/15. An outer edge lies on one triangle and the 4-cycle, so it is kept with probability 1 - (1/3 + 1/4 - 1/5) = 37/60 rather than 5/8. Nothing in that cycle arithmetic reproduces the uniform marginals.

The second shortcut, a randomised depth-first search (the classic recursive backtracker for mazes), is worse: on the same graph its trees have probabilities ranging from 1/12 to 7/48. On grids the bias is visible by eye: DFS mazes have long winding corridors and few dead ends, while uniform mazes have many short dead ends. Neither shortcut is wrong for a game, but if your application relies on the uniform distribution, such as an MCMC proposal, a statistical test, or a sparsifier with guarantees, the bias silently invalidates the result.

Aldous-Broder: first-entrance edges

The Aldous-Broder algorithm, found independently by Aldous and Broder around 1990, is the simplest exact sampler. Start a simple random walk at any vertex. Every time the walk enters a vertex it has never visited, add the edge it arrived on to the tree. Stop when every vertex has been visited. The resulting tree is uniform. The proof runs the walk backwards in time and uses the stationarity of the walk, but the procedure is three lines:

def aldous_broder(adj, rng, start=0):
    visited, tree, u = {start}, [], start
    while len(visited) < len(adj):
        v = rng.choice(adj[u])
        if v not in visited:
            visited.add(v)
            tree.append((u, v))      # first-entrance edge
        u = v
    return tree

Its running time is the cover time of the graph, the expected number of steps for a random walk to visit every vertex. That is fine on expanders and complete graphs but poor on graphs with bottlenecks or long paths: on a path or a cycle the cover time is quadratic in the number of vertices, and late in the run the walk spends most of its time revisiting vertices that are already in the tree.

Wilson algorithm: loop-erased random walks

David Wilson's 1996 algorithm is the standard choice. Pick a root and put it in the tree. For each vertex not yet in the tree, run a random walk from it until the walk hits the tree. Erase the loops from that walk in the order they formed, which leaves a simple path from the start vertex to the tree, and add the path to the tree. Repeat until every vertex is in.

Loop erasure sounds fiddly, but there is a neat trick: store, for each vertex, only the direction the walk most recently left it. Following those pointers from the start vertex gives exactly the loop-erased path, because any loop the walk made was later exited by a different edge and the pointer was overwritten.

import random

def wilson(adj, rng, root=0):
    """Uniform spanning tree by loop-erased random walks. adj: dict vertex -> list of neighbours."""
    in_tree = {root}
    nxt = {}
    for start in adj:
        # 1. random walk from start until it hits the tree, remembering only the LAST exit
        u = start
        while u not in in_tree:
            nxt[u] = rng.choice(adj[u])   # overwriting erases loops implicitly
            u = nxt[u]
        # 2. retrace the loop-erased path and add it to the tree
        u = start
        while u not in in_tree:
            in_tree.add(u)
            u = nxt[u]
    return [(u, nxt[u]) for u in nxt]     # every non-root vertex points to its parent

edges = [(0, 1), (0, 2), (1, 2), (1, 3), (2, 3)]
adj = {v: [] for v in range(4)}
for a, b in edges:
    adj[a].append(b)
    adj[b].append(a)
print(wilson(adj, random.Random(7)))

The output is uniform over all spanning trees, regardless of the root and of the order in which start vertices are processed. Its running time depends on the root, but Wilson showed the expected number of walk steps is governed by the mean hitting time of the graph, which is never larger than the cover time and is usually far smaller, because each walk stops as soon as it touches the growing tree instead of having to find the last unvisited vertex. On a 1000 by 1000 grid that is the difference between a run you wait for and one you do not.

Worked example on the diamond

Run Wilson by hand on the diamond with root 0. Vertex 0 is in the tree. Start at vertex 1: suppose the walk goes 1 to 3, 3 to 2, 2 to 1, 1 to 0. The pointers record nxt[3] = 2, nxt[2] = 1, and nxt[1] = 0, because vertex 1's pointer was overwritten by its second exit. Retracing from 1 gives the path 1 to 0 only: the loop 1, 3, 2, 1 has been erased. The tree is now {0, 1}. Start at vertex 2: its stale pointer is replaced when the new walk leaves it. Suppose the walk goes 2 to 3, 3 to 1. Vertex 1 is in the tree, so the walk stops, and the path 2, 3, 1 is added. Vertex 3 is now in the tree as well, so the loop over start vertices skips it. The final tree is {0-1, 1-3, 3-2}, one of the 8, and each of the 8 is reached with probability exactly 1/8. A quick simulation of 80,000 runs gives frequencies between 0.123 and 0.127 for every tree.

Weighted trees and effective resistance

For a weighted graph, the weighted UST picks each tree with probability proportional to the product of its edge weights. Wilson's algorithm handles this by making the walk step to a neighbour with probability proportional to the edge weight, and the matrix-tree theorem still counts the total weight if the Laplacian uses weighted degrees.

The deepest fact about USTs connects them to electrical networks. Treat each edge as a resistor of one ohm (or of resistance 1/weight). Then the probability that edge (u, v) belongs to a uniform spanning tree equals the effective resistance between u and v. On the diamond, the shared edge 1-2 is in parallel with two two-ohm paths, so its effective resistance is 1 / (1 + 1/2 + 1/2) = 1/2, and indeed 4 of the 8 trees contain it. This is Kirchhoff's theorem, and it is why USTs appear in spectral sparsification: sampling edges by effective resistance keeps the graph's cut and Laplacian structure, and a UST is a natural way to sample by that measure.

The relationship also works in reverse as a sanity check and as an estimator. If you need effective resistances on a graph too large for an exact Laplacian solve, sample a few hundred USTs with Wilson's algorithm and count how often each edge appears; the frequencies converge to the resistances. Conversely, if your sampler's edge frequencies disagree with resistances computed exactly on a medium-sized graph, the sampler is biased, even when it always returns valid trees. Every valid tree has exactly n - 1 edges, so the edge probabilities of a uniform tree always sum to n - 1, which is Foster's theorem: the effective resistances of all edges sum to n - 1. On the diamond the five resistances are 1/2 for the shared edge and 5/8 for each outer edge, and 1/2 + 4 times 5/8 = 3.

Engineering notes

• Iterate, never recurse. Walks and retracing are loops; a recursive implementation overflows the stack on large grids.
• Use arrays, not dictionaries, for large graphs. Store adjacency in CSR form and the pointer and in-tree flags in integer and boolean arrays.
• Choose the root sensibly. The distribution does not depend on the root, but the running time does. A high-degree or central root makes walks hit the tree sooner.
• Seed and record the RNG. Reproducible samples make debugging and experiments possible. Use a generator with good statistical quality, not a short linear congruential generator.
• Mazes. On a grid graph a UST is a perfect maze with a uniform distribution over all perfect mazes; draw walls for every grid edge not in the tree.
• Disconnected graphs. Run the algorithm per connected component, found with union-find or BFS; otherwise the walk from an unreachable vertex never terminates.

Testing a sampler

Samplers fail silently, so test the distribution, not just the output. On small graphs, count the trees with the matrix-tree theorem, draw many samples, and run a chi-square test against uniform. With 8 trees there are 7 degrees of freedom, and a statistic above about 24 is suspicious at the 0.1% level. The random-weight MST sampler fails this on the diamond with a large enough number of runs; Wilson passes.

from collections import Counter
import random

def chi_square_uniform(sampler, adj, n_trees, runs=80_000, seed=1):
    rng = random.Random(seed)
    counts = Counter(frozenset(tuple(sorted(e)) for e in sampler(adj, rng))
                     for _ in range(runs))
    assert len(counts) == n_trees, "some trees never appear"
    expected = runs / n_trees
    return sum((c - expected) ** 2 / expected for c in counts.values())   # df = n_trees - 1

Also check structural invariants on every output: exactly n - 1 edges, all vertices connected, no cycles. And check edge marginals on a larger graph against effective resistances computed by solving a Laplacian system, which catches subtler bias than per-tree counts can on graphs too big to enumerate.

Failure modes

• Using MST with random weights or a randomised DFS where uniformity matters.
• Forgetting loop erasure, adding every vertex the walk visits; the result can contain cycles or be biased.
• Walking on a directed or disconnected graph. Wilson's algorithm on a directed graph samples arborescences towards the root, a different object; on a disconnected graph it never terminates.
• Aldous-Broder on a long path or bottlenecked graph, paying the full cover time.
• Weights interpreted the wrong way round. In a weighted UST the step probability is proportional to conductance (weight), not to length.

What to do next

1. Implement Wilson's algorithm from the code above and run it on the diamond graph.
2. Compute the tree count with a Laplacian determinant and run the chi-square test.
3. Repeat with random-weight Kruskal and see it fail the same test.
4. Generate a 30 by 30 grid maze with Wilson and compare its dead-end count to a DFS maze.
5. Add weights and check one edge marginal against its effective resistance.
6. Read about simulated annealing and other MCMC methods, where uniform samplers serve as proposals: simulated annealing.