Markov chain Monte Carlo

The posterior of a real model, a distribution over molecular configurations or an energy-based model of images has a density that can be evaluated only up to a constant, $p(x) = \tilde p(x) / Z$ with $Z$ an intractable integral, in many dimensions. The direct methods of the Monte Carlo note need a proposal close to $p$, and there is none. Markov chain Monte Carlo turns the ergodic theorem around: build a Markov chain whose stationary distribution is $p$, run it, and use its states as samples. The construction needs only ratios $\tilde p(y)/\tilde p(x)$, in which $Z$ cancels. This last note of the course derives the Metropolis–Hastings algorithm, shows how it fails, and ends with Langevin dynamics, which follows the gradient of the log-density and supplies the score that diffusion models reuse.

Notation. $p(x) = \tilde p(x)/Z$ is the target; $q(y \mid x)$ a proposal distribution; $X_0, X_1, \dots$ the chain. $\nabla \log p(x)$ is the score, which does not depend on $Z$.

Detailed balance#

A chain with transition probabilities $P(x \to y)$ satisfies detailed balance with respect to $p$ if

$$p(x)\, P(x \to y) = p(y)\, P(y \to x) \quad \text{for all } x, y .$$

In equilibrium, the flow of probability from $x$ to $y$ equals the flow back. Summing over $x$ gives $\sum_x p(x) P(x \to y) = p(y)$: $p$ is stationary. Detailed balance is a local condition that is easy to arrange, and if the chain is also irreducible and aperiodic, the Markov chains note guarantees convergence to $p$ from any start for finite state spaces; continuous spaces need analogous conditions.

Metropolis–Hastings#

From the current state $x$, propose $y \sim q(\cdot \mid x)$ and accept it with probability

$$\alpha(x, y) = \min\Big(1, \frac{\tilde p(y)\, q(x \mid y)}{\tilde p(x)\, q(y \mid x)}\Big) ;$$

otherwise stay at $x$. For $y \ne x$ the transition probability is $q(y \mid x)\,\alpha(x, y)$, and $p(x) q(y \mid x) \alpha(x, y) = \min\big(p(x) q(y \mid x),\ p(y) q(x \mid y)\big)$ is symmetric in $x$ and $y$: detailed balance holds. The constant $Z$ appears in both numerator and denominator and cancels. With a symmetric proposal, such as $y = x + \sigma \varepsilon$ with $\varepsilon \sim N(0, 1)$, the ratio is simply $\tilde p(y)/\tilde p(x)$: moves uphill are always accepted, moves downhill sometimes. This is random-walk Metropolis, the algorithm of Metropolis and colleagues from 1953.

The chain is valid for any step size $\sigma$, but $\sigma$ decides how well it works:

Figure

Metropolis with too small, good and too large steps

Over 20,000 steps on a standard normal target, a step of $0.1$ accepts 97% of proposals and is worth only 57 independent samples; a step of $25$ accepts 5% and is worth 708; a step of $2.4$ accepts 44% and is worth 4,815. For random-walk Metropolis in many dimensions, theory suggests tuning the step to an acceptance rate of about 0.23; in one dimension, about 0.44.

Mixing and diagnostics#

Successive states of a chain are correlated, so $n$ states usually carry less information than $n$ independent samples. The effective sample size belongs to a quantity of interest $g(X_t)$, not to the chain as a whole: with $\rho_k$ the autocorrelation of $g(X_t)$ at lag $k$ in the stationary chain,

$$\text{ESS} = \frac{n}{1 + 2\sum_{k \ge 1} \rho_k} ,$$

and it measures how many independent samples the $n$ states are worth for estimating $E_p[g(X)]$. The error of the average of $g$ over the chain is then about $\sigma / \sqrt{\text{ESS}}$, not $\sigma/\sqrt n$, with $\sigma^2 = \operatorname{Var}_p g(X)$. This needs a chain that has reached stationarity, autocorrelations whose sum converges, and a central limit theorem for the chain, which holds for well-behaved chains but not for every one. Positive correlations make the ESS smaller than $n$; negative ones can make it larger: a chain whose autocorrelations are $\rho_k = (-0.5)^k$ has $\text{ESS} = 3n$. The ESS values quoted in this note come from a simple estimator that stops summing where the autocorrelations drop below $0.05$; they are illustrative, and production tools use more careful estimators. Practical use adds three habits: discard an initial burn-in before the chain has reached its stationary regime, run several chains from dispersed starting points, and compare the variance between chains with the variance within them (the $\hat R$ statistic).

No diagnostic can detect what a chain has never visited. For a target with two separated modes, a chain with a small step stays in the mode it starts in:

Figure

A chain stuck in one mode

In 20,000 steps the small-step chain spends 0% of its time in the right mode, the large-step chain 50%. The first chain’s samples look perfectly healthy; only a second chain started in the other mode, or knowledge of the target, reveals the problem. Multimodality is one of the main practical failures of MCMC, and tempering methods, which run chains on flattened versions of the target, exist to cross the barriers.

Gibbs sampling#

When the conditional distribution of each coordinate given the others is easy to sample, Gibbs sampling updates one coordinate at a time from its conditional. It is Metropolis–Hastings with the conditional as the proposal, and the acceptance probability is always 1. For the bivariate normal of the multivariate Gaussian note with unit variances and correlation $\rho$, the conditionals are $X_1 \mid X_2 = x_2 \sim N(\rho x_2, 1 - \rho^2)$ and symmetrically.

Figure

Gibbs sampling moves in a zigzag

Strong correlation makes Gibbs slow: each move is along an axis, and the conditional variance $1 - \rho^2$ is small. Successive draws of $X_1$ have correlation $\rho^2$, which is $0.81$ for $\rho = 0.9$ and $0.98$ for $\rho = 0.99$; 20,000 sweeps are then worth about 2,300 and 190 independent samples. Reparametrizing to remove the correlation, or updating correlated coordinates in blocks, fixes it. Gibbs sampling is the workhorse for models with conjugate conditionals, such as topic models and many hierarchical Bayesian models.

Langevin dynamics#

Random-walk proposals ignore where the probability is. Langevin dynamics uses the score to push each step uphill and adds noise to keep exploring:

$$x_{t+1} = x_t + \frac{\varepsilon}{2} \nabla \log p(x_t) + \sqrt{\varepsilon}\, z_t, \qquad z_t \sim N(0, I) .$$

This is the unadjusted Langevin algorithm (ULA). For a smooth target and small enough steps its stationary distribution tends to $p$ as $\varepsilon \to 0$; for a finite step it is biased, and for too large a step, or a score that grows quickly in the tails, the chain can be unstable and have no stationary distribution at all. On a standard normal target, whose score is $-x$, the update is $x_{t+1} = (1 - \varepsilon/2)\, x_t + \sqrt{\varepsilon}\, z_t$, which is stable only for $0 < \varepsilon < 4$; then the stationary variance is $1/(1 - \varepsilon/4)$: $1.026$ for $\varepsilon = 0.1$, $1.143$ for $\varepsilon = 0.5$. For $\varepsilon \ge 4$ the chain diverges. Adding a Metropolis–Hastings acceptance step removes the discretization bias (MALA): $p$ is then exactly stationary for every step size, and the chain converges to it under the usual conditions, although a large step gets few acceptances; simulating Hamiltonian dynamics with a momentum variable gives Hamiltonian Monte Carlo, which in its adaptive NUTS form is the default sampler of modern probabilistic programming systems such as Stan. With minibatch gradients, the same update is stochastic gradient Langevin dynamics, which turns noisy gradient descent into approximate posterior sampling for neural networks.

The score needs neither $Z$ nor the density itself, only the gradient of its logarithm, and this is the bridge to generative models. The bridge carries the score, not the whole sampler, and three procedures should be kept apart:

• Equilibrium Langevin and MALA, this note: one fixed target $p$, a chain run long enough to forget its start, samples from (approximately, or for MALA exactly) its stationary distribution $p$.
• Annealed Langevin dynamics (Song and Ermon, 2019): a network learns the scores $\nabla \log p_\sigma$ of the data blurred by Gaussian noise of several sizes $\sigma$, and a few Langevin steps are run at each level, from the largest to the smallest. Every level has its own target, the chain is not run to equilibrium at any of them, and the result is not an MCMC sample of the data distribution (score matching note).
• Reverse diffusion (DDPM, score-based SDEs, Song et al., 2021): generation simulates the time reversal of a noising process, from pure noise at the final time back to data, in a fixed number of steps. Its drift contains the score $\nabla \log p_t$ of the noised data at time $t$, so the transition rule changes from step to step; it is a time-inhomogeneous Markov chain with no stationary target (DDPM sampling). Langevin steps can be added at each time as a corrector, but the reverse process itself is not MCMC.

What changed#

	Direct Monte Carlo	Markov chain Monte Carlo
Needs	samples from $p$ or a close proposal	$\tilde p$ up to a constant (or its score)
Samples	independent	correlated; count the effective sample size
High dimensions	proposals degrade exponentially	local moves still work, at the cost of mixing time
Typical algorithms	inversion, rejection, importance sampling	Metropolis–Hastings, Gibbs, Langevin, HMC
Main failure	infinite-variance weights	slow mixing, chains stuck in one mode

Limitations#

• Convergence is asymptotic and hard to verify. A chain can look stationary while it has explored only part of the target; diagnostics detect some failures, never all.
• Tuning matters. Step sizes, proposal shapes, parametrizations and the number of chains change efficiency by orders of magnitude.
• It is expensive for big data and big models. Every step evaluates the target, often on the whole dataset; for neural networks variational inference, ensembles or stochastic-gradient MCMC are cheaper approximations.

What you should remember#

• Detailed balance $p(x) P(x \to y) = p(y) P(y \to x)$ makes $p$ stationary.
• Metropolis–Hastings: propose $y \sim q(\cdot \mid x)$, accept with $\min\big(1, \tilde p(y) q(x \mid y) / (\tilde p(x) q(y \mid x))\big)$; only ratios of $\tilde p$ are needed.
• Too small a step crawls, too large a step is rejected; aim for acceptance of roughly 0.23 in many dimensions.
• For a quantity $g(X_t)$ of a stationary chain, correlated samples are worth an effective sample size $n/(1 + 2\sum \rho_k)$, usually less than $n$; run several chains; beware of modes a chain never visits.
• Gibbs sampling updates coordinates from their conditionals; strong correlation makes it slow.
• Langevin dynamics follows the score plus noise; without a correction it is biased by the step size and can be unstable ($0 < \varepsilon < 4$ on a standard normal), MALA removes the bias, and HMC uses the same score. Annealed Langevin and reverse diffusion reuse learned scores of noised data, but with targets that change across noise levels or time, so they are not equilibrium MCMC.

Code#

Random-walk Metropolis with three step sizes, a two-mode target, Gibbs sampling and unadjusted Langevin dynamics:

import numpy as np, math

rng = np.random.default_rng(0)

def metropolis(logp, x0, step, n):
    """Random-walk Metropolis: propose x + step * N(0, 1), accept with probability min(1, p(y) / p(x))."""
    x, lp, chain, accepted = x0, logp(x0), np.empty(n), 0
    for t in range(n):
        y = x + step * rng.normal()
        ly = logp(y)
        if math.log(rng.random()) < ly - lp:
            x, lp, accepted = y, ly, accepted + 1
        chain[t] = x
    return chain, accepted / n

def ess(chain):
    """A rough effective sample size for illustration: autocorrelations summed until they drop below 0.05 (production tools use more careful estimators)."""
    c = chain - chain.mean()
    f = np.fft.rfft(c, 2 * len(c))
    ac = np.fft.irfft(f * np.conj(f))[: len(c)]
    ac /= ac[0]
    cut = np.argmax(ac < 0.05)
    return len(c) / (1 + 2 * ac[1:cut].sum())

# a standard normal target with three step sizes, 20,000 steps each: acceptance rate and effective sample size
for step in (0.1, 2.4, 25):
    chain, rate = metropolis(lambda x: -x * x / 2, 0.0, step, 20_000)
    print(step, round(rate, 3), round(ess(chain)))
# 0.1 0.971 57 | 2.4 0.443 4815 | 25 0.053 708

# two modes at -3 and 3: the share of time in the right mode with a small and a large step
mix = lambda x: math.log(0.5 * math.exp(-(x + 3) ** 2 / 0.72) + 0.5 * math.exp(-(x - 3) ** 2 / 0.72))
for step in (0.5, 4.0):
    chain, rate = metropolis(mix, -3.0, step, 20_000)
    print(step, round(float((chain > 0).mean()), 3))
# 0.5 0.0 | 4.0 0.502

# Gibbs for a bivariate normal: successive draws of x1 have correlation rho^2
for rho in (0.5, 0.9, 0.99):
    x1 = x2 = 0.0
    s, xs = math.sqrt(1 - rho ** 2), np.empty(20_000)
    for t in range(20_000):
        x1 = rho * x2 + s * rng.normal()
        x2 = rho * x1 + s * rng.normal()
        xs[t] = x1
    print(rho, round(float(np.corrcoef(xs[:-1], xs[1:])[0, 1]), 3), round(ess(xs)))
# 0.5 0.264 12075 | 0.9 0.807 2292 | 0.99 0.978 188

# unadjusted Langevin on N(0, 1), whose score is -x: the step size biases the stationary variance to 1 / (1 - eps / 4)
for eps in (0.1, 0.5):
    x, xs = 0.0, np.empty(200_000)
    for t in range(200_000):
        x = x + eps / 2 * (-x) + math.sqrt(eps) * rng.normal()
        xs[t] = x
    print(eps, round(float(xs[1000:].var()), 3), round(1 / (1 - eps / 4), 3))
# 0.1 1.026 1.026 | 0.5 1.13 1.143

Next#

This is the last note of the course. It started from coins, dice and the question of what a probability is, and ends with the samplers behind Bayesian statistics and the score that diffusion models reuse. The generative models course builds on almost every note here: maximum likelihood and KL divergence, the multivariate Gaussian, the change-of-variables formula, Monte Carlo estimates of expectations, and Langevin dynamics with the score it follows.

References#

• N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, E. Teller. Equation of State Calculations by Fast Computing Machines. The Journal of Chemical Physics, 21(6), 1953.
• W. K. Hastings. Monte Carlo Sampling Methods Using Markov Chains and Their Applications. Biometrika, 57(1), 1970.
• S. Geman, D. Geman. Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 6(6), 1984.
• G. O. Roberts, A. Gelman, W. R. Gilks. Weak Convergence and Optimal Scaling of Random Walk Metropolis Algorithms. The Annals of Applied Probability, 7(1), 1997.
• R. M. Neal. MCMC Using Hamiltonian Dynamics. In Handbook of Markov Chain Monte Carlo, CRC Press, 2011.
• M. Welling, Y. W. Teh. Bayesian Learning via Stochastic Gradient Langevin Dynamics. ICML, 2011.
• Y. Song, S. Ermon. Generative Modeling by Estimating Gradients of the Data Distribution. NeurIPS, 2019. Annealed Langevin dynamics.
• Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, B. Poole. Score-Based Generative Modeling through Stochastic Differential Equations. ICLR, 2021. Reverse-time SDEs and predictor–corrector sampling.