Metropolis-Hastings algorithm

The Metropolis-Hastings algorithm is the foundational method of Markov Chain Monte Carlo (MCMC) sampling — a technique for drawing samples from probability distributions that are too complex to sample from directly. Developed by Nicholas Metropolis and colleagues in 1953 and generalized by W. K. Hastings in 1970, the algorithm constructs a Markov chain whose stationary distribution is the target distribution of interest. By simulating this chain and collecting samples after a burn-in period, one can approximate expectations, quantiles, and other properties of distributions that resist analytic treatment.

The algorithm's core insight is deceptively simple: propose a new state, decide whether to accept it based on how probable it is relative to the current state, and repeat. This local exploration strategy, repeated thousands or millions of times, produces a random walk through the state space that eventually visits regions in proportion to their probability mass. The algorithm requires only that one can evaluate the target distribution up to a proportionality constant — a minimal requirement that makes it applicable to Bayesian inference, statistical physics, combinatorial optimization, and any domain where probabilities are known structurally but not normalizing constants.

The Algorithm

Given a target distribution π(x) and a proposal distribution q(x'|x), the Metropolis-Hastings algorithm proceeds as follows:

1. Initialize at some state x₀.
2. Propose a new state x' by sampling from q(x'|xₜ).
3. Compute the acceptance ratio α = min(1, [π(x')q(xₜ|x')] / [π(xₜ)q(x'|xₜ)]).
4. With probability α, accept the proposal and set xₜ₊₁ = x'. Otherwise, reject and set xₜ₊₁ = xₜ.
5. Repeat from step 2.

The acceptance ratio ensures detailed balance with respect to π: the chain is as likely to move from x to x' as from x' to x, when weighted by the stationary probabilities. This balance, combined with ergodicity conditions on the proposal, guarantees that the chain converges to π regardless of the initial state.

The critical design choice is the proposal distribution. A random walk proposal — adding Gaussian noise to the current state — is simple but can mix slowly in high dimensions or multimodal landscapes. An independent proposal — sampling from a fixed distribution — can jump between modes but requires a good approximation of the target. Adaptive proposals that tune their spread based on past samples have been developed, but they sacrifice the Markov property and require careful theoretical justification.

Why It Works: The Systems View

The Metropolis-Hastings algorithm is not merely a statistical tool. It is a dynamical system that solves an inference problem by local exploration. Viewed through the lens of dynamical systems theory, the algorithm constructs a Markov chain whose transition matrix has π as its dominant eigenvector — the unique fixed point toward which all initial distributions converge. The rate of convergence is governed by the spectral gap: the difference between the two largest eigenvalues of the transition matrix. A large spectral gap means fast mixing; a small gap means the chain lingers in local regions for long periods, producing highly correlated samples and unreliable estimates.

This dynamical perspective reveals why the algorithm struggles in high dimensions. As the dimensionality of the state space grows, the volume of the space grows exponentially, while the proposal distribution typically explores only a local neighborhood. The chain becomes trapped in regions of high probability separated by vast low-probability deserts — a phenomenon analogous to glassy dynamics in statistical physics. The effective sample size — the number of effectively independent samples produced by the chain — collapses, and the computational cost of reliable inference explodes.

The connection to self-organized criticality is subtle but real. A well-tuned Metropolis-Hastings chain operates near a critical point in its parameter space: the proposal variance must be large enough to escape local modes but small enough to maintain a reasonable acceptance rate. Too large, and proposals are almost always rejected; too small, and the chain diffuses too slowly. The optimal variance balances exploration and exploitation — a trade-off that appears in evolutionary algorithms, reinforcement learning, and any system that must search a complex landscape.

Applications and Extensions

In Bayesian inference, the Metropolis-Hastings algorithm enables posterior computation for models with thousands of parameters and complex hierarchical structures. It underlies probabilistic programming languages like BUGS and early versions of Stan, and it remains the default sampler when gradients are unavailable or expensive. The algorithm's flexibility — it requires only likelihood evaluations, not gradients or conjugacy — makes it indispensable for structural equation models, mixture models, and latent variable frameworks.

In statistical physics, the algorithm samples configurations from the Boltzmann distribution, enabling the study of phase transitions, critical phenomena, and disordered systems. The Ising model, spin glass models, and polymer simulations all rely on MCMC methods descended from Metropolis-Hastings.

Extensions of the algorithm address its limitations. Hamiltonian Monte Carlo uses gradient information to propose states that follow the geometry of the target, dramatically improving mixing in high dimensions. Parallel tempering runs multiple chains at different temperatures, allowing swaps that escape local modes. Reversible jump MCMC extends the framework to variable-dimensional spaces, enabling model comparison and selection.

Limitations and Open Problems

Despite its generality, the Metropolis-Hastings algorithm has fundamental limitations. The mixing time — the number of iterations required for the chain to approximate its stationary distribution — can be exponentially large in the problem size for complex targets. Diagnosing convergence remains an art: standard diagnostics like the Gelman-Rubin statistic detect obvious failures but can miss subtle forms of non-convergence.

The algorithm also embodies a form of representational debt. By constructing a Markov chain that explores the state space locally, it discards global structural information about the target distribution. When the target has complex geometry — narrow ridges, strong correlations, multimodal structure — the local exploration strategy is intrinsically inefficient. The information needed for efficient sampling is present in the target's geometry, but the algorithm has no mechanism to extract and use it. Hamiltonian Monte Carlo partially repays this debt by using gradients; variational methods repay it by using parametric approximations. But for targets that are neither gradient-friendly nor parametrically approximable, the debt remains.

The open frontier is the design of proposal mechanisms that adapt to the target's geometry without sacrificing theoretical guarantees. Recent work on Langevin dynamics, Riemannian manifold MCMC, and normalizing flows as proposals represents attempts to build smarter explorers. Whether these methods can be unified into a single framework — or whether the problem of efficient sampling is fundamentally as hard as the inference problem it solves — remains open.

The Metropolis-Hastings algorithm is the original sin of computational inference: it solves the problem of sampling from complex distributions by making a problem that is hard in principle easy in practice — at the cost of making it slow in practice. Its enduring power is not its efficiency but its universality: it works for any target, any proposal, any dimension. But this universality is also its weakness. The algorithm does not learn from the target's structure; it blindly proposes and cautiously accepts, like a bureaucrat who processes every form correctly but never asks whether the process could be streamlined. In an era where machine learning has made adaptivity the default, the Metropolis-Hastings algorithm's stubborn non-adaptivity is both a philosophical stance — all distributions are equal — and a computational dead end. The next generation of samplers will not be universal. They will be structural. They will know what they are sampling from, and they will use that knowledge to explore with purpose rather than persistence.