Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images

A camera captures a scene, but noise corrupts every pixel. The observed image $y$ is a degraded version of the true image $x$: perhaps $y = x + \text{noise}$, or something more complex involving blur and nonlinear distortion. The restoration problem asks: given $y$, what was $x$?

This is an inverse problem, and it is ill-posed — many plausible $x$'s could produce the same $y$. To pick the best one, we need two ingredients from Bayesian inference:

• A prior $P(x)$ — what do clean images generally look like? (Neighboring pixels tend to be similar; edges are sparse.)

• A likelihood $P(y|x)$ — given a clean image $x$, how probable is the observed noisy image $y$?

Bayes' theorem combines them into the posterior $P(x|y) \propto P(y|x) \, P(x)$. The best restoration is the $x$ that maximizes this posterior — the MAP estimate.

The challenge: with a 256×256 image and 256 gray levels, the space of possible images has $256^{65536}$ states. No algorithm can examine them all. We need a way to explore this space efficiently, guided by the posterior. That's where Gibbs distributions and the Gibbs sampler enter.

The paper's first key insight is a deep connection between two seemingly different worlds.

On one side: Markov random fields (MRFs) — a way to model images where each pixel's value depends only on its immediate neighbors, not on distant pixels. This is a natural model for images: whether a pixel is bright or dark depends mostly on the pixels right around it.

On the other side: Gibbs distributions from statistical physics — probability distributions over configurations of a system, defined through an energy function. In physics, atoms in a crystal arrange themselves to minimize energy. Low-energy configurations (orderly crystal lattices) are exponentially more probable than high-energy ones (disordered states).

The Hammersley-Clifford theorem establishes that these are the same thing: every MRF is a Gibbs distribution, and vice versa. This means we can design image models by writing down an energy function — a much more intuitive task than specifying conditional probabilities for every pixel.

Think of it as a landscape: the energy function defines hills and valleys over the space of all possible images. Each valley is a plausible image. The Gibbs distribution tells you how probable each point in this landscape is — the deeper the valley, the more probable the configuration. Temperature controls how sharply probability concentrates in the valleys: at high temperature the system explores everywhere; at low temperature it settles into the deepest valleys.

The energy $U(x)$ is built by summing small, local contributions called clique potentials. A clique is a group of pixels that are all neighbors of each other — in the simplest case, a pair of adjacent pixels.

For image restoration, the energy has two parts:

• Prior energy — penalizes configurations that don't look like natural images. A simple choice: $\sum_{\text{neighbors } i,j} (x_i - x_j)^2$. This penalizes large differences between neighbors, encouraging smooth regions. More sophisticated priors add edge-preserving terms that allow sharp boundaries.

• Data energy — penalizes configurations that disagree with the observed noisy image. Typically: $\sum_i (y_i - x_i)^2$ for additive Gaussian noise.

The total energy $U(x) = \text{Prior energy} + \lambda \cdot \text{Data energy}$ balances smoothness against fidelity. The parameter $\lambda$ controls the tradeoff: large $\lambda$ trusts the data more; small $\lambda$ trusts the prior more.

The Gibbs sampler is elegantly simple. To sample from a complicated joint distribution $P(x_1, x_2, \ldots, x_n)$, you don't need to know the whole distribution at once. You just need to know how to sample each variable conditioned on all the others.

The algorithm cycles through all pixels. For each pixel $x_i$:

• Fix every other pixel at its current value.

• Compute the conditional distribution $P(x_i \mid x_{\text{neighbors}})$ — this depends only on the neighbors because of the Markov property.

• Draw a new value for $x_i$ from this conditional distribution.

After visiting every pixel once (a sweep), you've completed one iteration. The key theorem: as the number of sweeps grows, the distribution of the entire configuration converges to the target joint distribution $P(x)$. Each sweep brings the image closer to a sample from the posterior.

Gibbs sampling at a fixed temperature generates samples from the posterior — useful for computing averages and uncertainties. But for image restoration, Geman & Geman wanted the single best image: the MAP estimate. This is an optimization problem.

Here's the trick from physics: simulated annealing. At high temperature, the Gibbs distribution is nearly uniform — the sampler explores freely. As temperature decreases, probability concentrates on lower-energy states. At $T \to 0$, all probability sits on the global minimum.

The annealing schedule must cool slowly enough. Geman & Geman proved that if the temperature decreases as $T(t) = c / \log(1 + t)$ where $c$ is large enough, the algorithm converges to the global optimum with probability 1. This was the first rigorous convergence guarantee for simulated annealing on a general problem.

In practice, the logarithmic schedule is too slow. Practitioners use faster cooling (exponential or linear) and accept approximate solutions — which are usually excellent. The theoretical guarantee matters because it showed that stochastic methods can solve global optimization in combinatorial spaces, a result that influenced optimization across all of computer science.

The Gibbs sampler is a special case of the Metropolis-Hastings algorithm. Both are MCMC methods that build a Markov chain whose stationary distribution is the target. The difference lies in how they propose updates:

• Metropolis-Hastings proposes a random change and accepts or rejects it based on an acceptance ratio. Rejections waste computation — the chain stays put.

• Gibbs sampling proposes from the exact conditional distribution, so every proposal is accepted. No wasted steps, no tuning of proposal distributions.

The tradeoff: Gibbs sampling requires being able to compute and sample from each conditional $P(x_i \mid x_{-i})$ — easy for MRFs (thanks to the Markov property) but hard for models without conditional conjugacy. Metropolis-Hastings works for any model as long as you can evaluate the density ratio.

Two convergence results underpin the paper:

Ergodicity of the Gibbs sampler. At any fixed temperature $T > 0$, the Markov chain generated by Gibbs sampling is ergodic — it can reach any configuration from any other. This means it has a unique stationary distribution, and that distribution is exactly the Gibbs distribution. After enough sweeps, the samples are faithful representatives of the posterior.

Annealing convergence. If the temperature decreases slowly enough — specifically, $T(t) \geq c / \log(1+t)$ for a constant $c$ related to the energy barrier height — then the probability of being in a global minimum converges to 1. The sampler doesn't just explore: with annealing, it optimizes.

Together, these results mean: Gibbs sampling gives you correct samples (for Bayesian inference), and Gibbs sampling with annealing gives you the global optimum (for MAP estimation).

The paper's legacy extends far beyond image restoration. The Gibbs sampler became the computational engine behind Bayesian inference in thousands of applications. Every time a researcher fits a hierarchical model, runs a topic model, or does Bayesian structure learning, they are building on the foundations that Geman & Geman laid in 1984.