A walk that samples a distribution

Importance sampling assumed that a proposal distribution could be drawn from, and drawing from a distribution turns out to be the hard part. A density known by a formula is not a density anybody can sample: inverting its cumulative requires the cumulative, and computing that requires the very integral the sampling was meant to estimate.

The situation this rung is about is worse and more common. A weight $w(x)$ can be evaluated at any point, and the distribution wanted is $w$ divided by its total — a total nobody can compute, over a space too large to enumerate.

The rule #

Stand at state $i$. Propose a move to a neighbour $j$, chosen by any rule that is as likely to propose $j$ from $i$ as $i$ from $j$. Then:

• if $w(j) \geq w(i)$, move;
• otherwise move with probability $w(j)/w(i)$, and stay put otherwise.

That is the whole algorithm. It is due to Metropolis and four colleagues in 1953, and it is one of the ten or so algorithms that changed what could be computed.

The striking part is what is absent. The normalising constant never appears, because the rule only ever uses the ratio of two weights, and a constant cancels from a ratio. So a distribution can be sampled without anybody knowing what it is — the thing nobody can compute is the thing nobody needs.

Why it works #

The chain’s behaviour is settled by whether the target distribution $\pi$ is left alone by one step, and there is a condition sufficient for that which is far easier to check than the condition itself.

Detailed balance. For every pair of states,

$$\pi(i)P(i,j) = \pi(j)P(j,i),$$

where $P(i,j)$ is the chance of moving from $i$ to $j$. Read it as a statement about flows: in the long run as much probability crosses from $i$ to $j$ as crosses back.

If that holds, summing over $i$ gives $\sum_i \pi(i)P(i,j) = \pi(j)\sum_i P(j,i) = \pi(j)$, since a row of the transition matrix adds to 1. And $\sum_i \pi(i)P(i,j) = \pi(j)$ is exactly the statement that a step leaves $\pi$ where it was.

The Metropolis rule satisfies detailed balance by construction. Suppose $w(j) < w(i)$. Then $P(i,j) = q\cdot w(j)/w(i)$ and $P(j,i) = q$, where $q$ is the proposal probability, the same both ways. So

$$\pi(i)P(i,j) = \frac{w(i)}{Z}\cdot q\frac{w(j)}{w(i)} = \frac{q\,w(j)}{Z} = \pi(j)P(j,i),$$

and the $w(i)$ that would have made this messy has cancelled. The acceptance probability was chosen precisely so that it would.

The figure checks this rather than reciting it. It builds the whole transition matrix, verifies every row adds to one and holds no negative entry, and then asserts $\pi(i)P(i,j) = \pi(j)P(j,i)$ for every one of the pairs, to twelve decimal places. It then separately verifies that $\pi P = \pi$, which is the conclusion, computed independently of the argument that produces it.

The distribution as an eigenvector #

There is a second reading of $\pi P = \pi$ that connects this to a different part of the collection.

Written with $\pi$ as a row, that equation says $\pi$ is a left eigenvector of $P$ with eigenvalue 1. So the target distribution is an invariant direction of the transition matrix, and the whole apparatus of eigenvalues applies.

In particular the second eigenvalue decides how fast the chain forgets its start. A distribution that is not $\pi$ decomposes into $\pi$ plus a combination of the other eigenvectors, and each of those is multiplied by its own eigenvalue at every step. Since all the others are smaller than 1 in size, they die away — at a rate governed by the largest of them, which is called the spectral gap.

That is the power method, running by itself. Nobody chose to iterate a matrix; a random walk was set going, and the arithmetic of iterating a matrix is what governs it.

That reversal is worth dwelling on. The rung below in probability starts from a chain and asks where it goes; this one starts from where it should go and asks what chain would go there. The second question sounds harder and is easier, because detailed balance gives a sufficient condition that can be arranged by choosing one number — the acceptance probability — and nothing else.

What the samples are for #

A chain that visits states in proportion to their weight is an integrator. The average of any quantity over the states visited converges to the average under $\pi$, so anything that can be written as an expectation can be computed by running the chain and averaging.

That covers a great deal. In statistical physics the weight is $e^{-E/kT}$ — an exponential, which is why energy differences turn into probability ratios and the states are configurations of a system with $10^{23}$ parts; the quantity wanted is an average energy or magnetisation. In Bayesian statistics the weight is a likelihood times a prior and the states are parameter values; the quantity wanted is a posterior mean or a credible interval. In both cases the normalising constant is an integral over an enormous space, and in both cases it is never needed.

The proposal is free, and that matters #

Nothing in the argument constrained the proposal beyond symmetry: as likely to suggest $j$ from $i$ as $i$ from $j$. The figures use a step to a neighbour, and any symmetric rule would do — a jump to a random state, a step of a random size, a swap of two coordinates.

That freedom is where all the engineering is. A proposal making small steps is accepted often and moves slowly; one making large steps moves quickly when accepted and is rejected most of the time. There is an optimum between the two, and for a broad class of problems it is at an acceptance rate of about a quarter — a number that is genuinely surprising, since it says that throwing away three proposals in four is efficient.

Even the symmetry can be dropped, at the cost of one extra factor. The Metropolis–Hastings rule accepts with probability $\min\!\left(1, \frac{w(j)q(j,i)}{w(i)q(i,j)}\right)$, which reduces to the Metropolis rule when $q$ is symmetric and handles any proposal otherwise. Detailed balance still holds, by the same cancellation with one more term in it.

What is given up #

The samples are not independent, and that changes the arithmetic of how much a run is worth.

For $N$ independent draws the error is $\sigma/\sqrt N$. For a chain, consecutive states are correlated — the walk moves one step at a time — so $N$ steps of a chain are worth fewer than $N$ independent draws. The number they are worth is the effective sample size, $N$ divided by a factor built from the correlation between steps, and that factor can be enormous.

So the rate is still $1/\sqrt{N}$ in the effective count, and the effective count can be a tiny fraction of the actual one. All the practical difficulty of these methods is in that fraction.

There is also a start-up problem. The chain begins wherever it was put, which is generally not a typical state, and the early steps reflect that rather than the target. The usual response is to discard an initial segment — a burn-in — and the usual difficulty is that nobody can say how long it should be, because the honest answer depends on the spectral gap and the spectral gap is not available.

The failure that looks like success #

The characteristic way these methods go wrong is worth stating plainly, because it is the reason the field has a large literature on diagnostics.

A chain proposing small steps will explore its immediate neighbourhood thoroughly and may never cross a region of low weight to reach a second mode. Every diagnostic computed from such a run looks healthy: the samples are well mixed within the mode, the averages are stable, the estimated error is small. The result is confidently wrong, and the amount by which it is wrong is the weight of everything the chain never visited.

This is the same failure as the missed spike in importance sampling, one level worse. There, a badly chosen proposal might miss a region; here, the chain can be unable to reach it in any practical time, and the run gives no sign.

The mitigations are all about making the chain move further: larger proposals, several chains started in different places, tempering schemes that flatten the distribution to allow crossing and then sharpen it again. None of them is a guarantee, and none of them can be, since establishing that a chain has explored a space is at least as hard as exploring it. It is the shape of problem that an unprovable sentence has in a different setting: the check is of the same difficulty as the thing checked.

The comparison is the point of the whole construction. An unbiased walk spreads out forever and has no long-run distribution at all on an infinite line; confined to a finite set it settles into the uniform one. Weighting the acceptance is the only change, and it converts a process with one fixed distribution into a process with whatever fixed distribution is wanted.

Why a chain rather than rejection #

Rejection sampling produces genuinely independent draws and is exact, so it is worth being clear about why it is not simply used instead.

Rejection needs an easy density that everywhere exceeds the target after scaling, and keeps a draw with probability equal to their ratio. The fraction kept is the ratio of the two volumes — and in high dimensions a bounding density and a target density overlap in a fraction that shrinks exponentially with the dimension, however carefully the bound is chosen. At a hundred dimensions the acceptance rate is a number with a great many zeros after the point.

The chain has no such factor. Its cost is in correlation rather than in rejection, and correlation degrades a run by a constant factor while rejection degrades it exponentially. That is the entire trade, and it is why every method that works in high dimensions has given up independence.

Where the constant would have been #

It is worth being concrete about the size of the thing that cancels, because the cancellation is the reason the method exists and it is easy to read as a technicality.

For a system of $n$ particles each in one of two states, the normalising constant is a sum over $2^n$ configurations. At $n = 300$ that is more terms than there are atoms in the observable universe. No amount of computing power touches it, no clever rearrangement helps in general, and it is not a quantity anybody will ever have.

And it is required by every direct approach. Writing down the probability of a single configuration requires it. Drawing a sample by inverting a cumulative requires it. Computing an expectation as a weighted sum requires it.

The Metropolis rule requires the ratio of two configurations’ weights, which for a physical system is the difference of two energies, which is a local computation involving a handful of particles. A quantity nobody can compute has been replaced by one that costs almost nothing, and the replacement is exact rather than approximate.

That is the whole reason the algorithm mattered enough to change what physics could calculate, and it is why its authors were computing the equation of state of a hard-sphere gas on a machine with a few kilobytes of memory in 1953.

What the picture cannot show #

The figures use a dozen or so states so that the transition matrix can be written down and every pair checked. Nothing that is hard about these methods is visible at that size: with twelve states the chain mixes in a few hundred steps, every mode is reachable, and there is no question of a diagnostic being misleading. The genuine subject is a space of $10^{100}$ configurations, and its difficulties are entirely absent here.

That is a deliberate trade rather than an oversight. On twelve states the two claims the method rests on — detailed balance, and the target being fixed by a step — can be checked exhaustively, pair by pair, which is worth more than a picture of a case where they cannot. What is shown is that the argument is correct; what is not shown is that the algorithm is practical.

Nor does the trace of the walk carry any information about correlation, which is the quantity that decides how much a run is worth. A trace of three hundred steps looks like a trace of three hundred steps whether the effective sample size is three hundred or three.

The ladder from here #

Rungs above: the Gibbs sampler, which updates one coordinate at a time and needs no acceptance step. Hamiltonian Monte Carlo, which proposes long moves by simulating a physical trajectory and is why high-dimensional Bayesian computation became practical. The spectral gap and mixing times, and the coupling arguments that bound them. Simulated annealing, which is this algorithm with the temperature lowered and is an optimiser rather than a sampler. And the convergence diagnostics, and the standing difficulty that none of them can prove what it is asked to.

The idea in one sentence #

A distribution nobody can compute was sampled by a process that never asks for it.

That is worth holding on to as a pattern rather than an algorithm. The quantity that is out of reach — the normalising constant — was arranged to cancel, by building a rule that uses only ratios; and the property that was wanted — visiting states in the right proportion — was obtained by imposing a local condition, detailed balance, that could be checked one pair at a time.

A global property secured by a local condition is the shape of a great many arguments, and it is what makes them checkable. Nothing in the figures verifies that the chain visits states in the right proportion by watching it; what is verified is a condition on every pair of states, from which the global behaviour follows.