Mean-field Boltzmann machines

What this chapter covers

Here is a picture with 18% of its pixels flipped to the wrong color, and here is the belief that fixes it: a pixel probably looks like its neighbors. That single prior, written down as an energy and pushed to its logical conclusion, is enough to walk a noisy binary image back toward the clean one it came from — no training set, no labels, no gradient descent on a loss. This chapter builds the model that does it: a Boltzmann machine over the pixels of an image, the workhorse energy-based model, and solves it with mean-field variational inference, the cheap approximate-inference trick that turns an intractable posterior into a handful of tanh calls.

The setup is an Ising model borrowed straight from physics. Every pixel is a little magnet that can point up or down; neighbors want to agree, and each pixel also feels a pull toward whatever the noisy image says it should be. Those two forces trade off, and the configuration that balances them is a clean image. Computing the exact answer is hopeless — there are $2^{576}$ configurations of a 24×24 grid — so we don't. We approximate the whole posterior with one number per pixel, its posterior mean, and update those numbers one at a time until they stop moving. That loop is the entire method, and you get to watch it run: the animation in the middle of this page is a real denoising, sweep by sweep, noise getting outvoted by agreement while a free-energy panel slides downhill beside it.

The thing to keep your eye on is that free energy. Like the log-likelihood in the EM chapter, it's the one number that certifies the loop is working — it only ever goes down, every sweep, and when it flattens the image is done. EM and this are close cousins; both are coordinate descent on a variational bound, and if you learned to trust EM's monotone climb you already know how to read this descent. At the end we race the result against the denoiser everyone actually reaches for, a median filter, and I'll show you exactly where the energy model earns its keep and where it doesn't.

A bit of history

The model is older than machine learning. In 1920 Wilhelm Lenz handed his student Ernst Ising a toy: a lattice of spins, each ±1, each preferring to align with its neighbors, and the question of whether that local preference could produce global order — magnetism. Ising solved the one-dimensional case in his 1925 thesis and concluded, wrongly for higher dimensions, that it never magnetized. The model kept his name anyway, and it became the fruit fly of statistical physics: the simplest system where local rules add up to collective behavior.

The leap from physics to computation came in the early 1980s. John Hopfield's 1982 paper showed that a network of these ±1 units, wired symmetrically, has an energy that can only decrease as the units flip — so the network settles into minima you can use as stored memories. Two years later Geoffrey Hinton and Terry Sejnowski turned Hopfield's deterministic net into a probabilistic one: their Boltzmann machine, named in their 1985 paper with David Ackley, replaced the hard flip with a stochastic one governed by the Boltzmann distribution, which made the units samples from a real probability model and gave the network a learning rule. That's the lineage this chapter sits in — an energy over binary units, read as a probability.

The other half of the story is the inference trick. Sampling a Boltzmann machine to convergence is slow, and in 1987 Carsten Peterson and James Anderson proposed skipping it: replace each unit's fluctuating state with its average, and solve for the averages directly. That's the mean-field approximation, lifted straight from the Weiss theory of magnetism, and the fixed-point equation they wrote down is the tanh update we build below. In parallel, Stuart and Donald Geman's 1984 paper put exactly this kind of model to work on images — a Markov random field prior for image restoration — and by the late 1990s the variational-inference program of Michael Jordan, Zoubin Ghahramani, and their collaborators had reframed mean field as one instance of a general idea: approximate a hard posterior with a simple family and minimize the distance to it. Every piece of that history is in the code on this page.

The intuition

Think of the clean image as the truth and the noisy image as a corrupted witness. Salt-and-pepper noise flips scattered pixels to the wrong color, so the witness is right most of the time but lies at random about roughly one pixel in five. You want to recover the truth, and you have exactly one piece of outside knowledge: real images are smooth. Ink sits next to ink, background next to background. A lone black pixel marooned in a white field is almost certainly a lie.

So set up a tug of war. Give every pixel two reasons to be one color or the other. The first is agreement: a pixel that matches its four neighbors is rewarded, so smooth regions are cheap and isolated flips are expensive. The second is loyalty to the witness: a pixel that matches what the noisy image reported is also rewarded, so the reconstruction can't just wander off into a blank page. Tune the two rewards against each other and the resolution is automatic. A flipped pixel surrounded by four neighbors that all disagree with it feels four units of pressure to switch and only one unit of loyalty holding it in place, so it switches — the noise gets outvoted. A pixel on a real edge, where neighbors genuinely disagree, feels a split vote and defers to the witness. The prior cleans the easy pixels and stays out of the way on the hard ones.

That is the whole idea, and it's worth seeing the input before the math. On the left is the clean 24×24 image — a digit, binarized to black ink on a white ground. On the right is the same image after I flipped 18% of its pixels at random. That speckled mess on the right is all the algorithm ever sees.

You can still read the digit through the speckle, which is the point: your visual system is doing mean-field inference for free, letting the intact majority outvote the scattered lies. We're about to write that down as arithmetic.

The math

Encode the image as spins. Each pixel $i$ gets a hidden value $x_i \in \{-1, +1\}$ — think $+1$ for ink, $-1$ for background — and we write $y_i \in \{-1, +1\}$ for the noisy pixel we actually observed. The model scores a whole configuration $x$ with an energy, low for the configurations we believe in:

The first sum runs over neighboring pixel pairs $\mathcal{E}$ (each pixel joined to the four around it); the product $x_i x_j$ is $+1$ when two neighbors agree and $-1$ when they clash, so with $J > 0$ agreement lowers the energy — the smoothness prior. The second sum is the data term: $x_i y_i$ is $+1$ when a pixel matches its noisy observation, so with $h > 0$ matching the witness lowers the energy too. $J$ sets how much we trust smoothness, $h$ how much we trust the noisy image. Turn the energy into a probability with the Boltzmann distribution, where low energy means high probability:

That $Z$ is the catch. It sums over all $2^{576}$ configurations of the grid, so we can neither compute $p(x)$ nor sample from it directly, and the marginal we actually want — the probability that a given pixel is ink — is buried inside it. This is where mean field earns its place. Instead of the true posterior we use a stand-in $q$ that pretends the pixels are independent, a product of one-pixel distributions:

Each factor $q_i$ is a coin over $\{-1, +1\}$, and it's fully described by its mean $\mu_i = \mathbb{E}_q[x_i] \in [-1, +1]$ — a soft, real-valued pixel, $+1$ for "surely ink," $-1$ for "surely background," $0$ for "no idea." Fitting $q$ to $p$ means choosing the $\mu_i$ that make the two distributions as close as possible in KL divergence, and doing that calculus for one pixel with the others held fixed gives a clean fixed-point equation. Each pixel's mean is a tanh of the total field pressing on it:

Read it straight off the energy: the field on pixel $i$ is $J$ times the sum of its neighbors' current means plus $h$ times its own observation, and tanh squashes that field into a mean in $[-1, +1]$. Sweep this over every pixel, again and again, and the means settle. What are they settling into the bottom of? The mean-field free energy, the objective the whole loop minimizes:

Under the independent $q$, the expected energy is just the Ising energy evaluated on the means (because $\mathbb{E}_q[x_i x_j] = \mu_i \mu_j$), and $H(q)$ is the summed entropy of the per-pixel coins. The reason $F$ is the right thing to minimize is that it differs from the log of that intractable $Z$ by exactly the approximation error:

Since $\log Z$ is a constant we can't compute but also don't need, driving $F$ down drives $\mathrm{KL}(q \,\|\, p)$ down — it pushes our factorized approximation as close to the true posterior as the independence assumption allows. And each single-pixel tanh update is the exact minimizer of $F$ in that one coordinate, so a sweep of them can only lower $F$. That is the same coordinate-descent-on-a-bound guarantee that makes EM trustworthy, and the animation below makes you watch it hold.

What it's good at, what it isn't

The strength is that mean field makes an impossible computation cheap and turns it into something you can read. The exact posterior over a 24×24 image lives in a space of $2^{576}$ configurations with a normalizer nobody can sum; mean field replaces it with 576 numbers and a fixed-point iteration that is nothing but neighbor averages passed through tanh. There's no training phase, no data beyond the single image, no learning rate, no step that can diverge — just a monotone slide down the free energy to a fixed point. When your prior really is local and pairwise, which is exactly the case for smoothness on an image, this is about the cheapest inference that still respects the whole grid's structure, and it parallelizes and scales the way a convolution does.

The weakness is the assumption that bought the speed. Mean field declares the pixels independent, and they are emphatically not — that's the entire content of a smoothness prior. So the approximation is systematically overconfident: it drives each $\mu_i$ toward a hard $\pm 1$ faster than the true posterior would, and it cannot represent the correlations that carry an image's fine structure. On a thin stroke or a sharp corner, where the right answer depends on neighbors agreeing in a coordinated way the factorized $q$ can't express, mean field smooths where it shouldn't and rounds the corner off. It also finds a local minimum of $F$, not the global one, so the fixed point depends on where you start and on the order you sweep. And the two knobs $J$ and $h$ are yours to set: crank $J$ too high and the model melts the whole image into one flat blob; too low and it copies the noise straight through. It's a fast, honest approximation, not the exact answer, and the places it fails are the places the pixels were talking to each other behind its back.

The data

One 24×24 binary image, and its noisy twin. I built the clean image from scikit-learn's load_digits — an 8×8 handwritten "3" — by thresholding it to black and white and then upsampling each pixel into a 3×3 block to get a 24×24 grid, 576 pixels, 189 of them ink. That gives a shape with genuine structure to protect: strokes, curves, a hole in the middle, the kind of thing a smoothness prior can help with and also the kind of thing it can over-smooth. Then the corruption: I seeded a random generator and flipped 18% of the pixels — 96 of the 576 — from their true color to the opposite, salt-and-pepper noise. The clean grid and the noisy grid are both committed as CSV under data/, and the seed is fixed, so every number on this page reproduces exactly. The algorithm is handed only the noisy grid; the clean one exists solely so we can grade the result at the end.

Build it, one function at a time

Six short functions, and four of them are three lines. The model is the energy, the inference is the tanh update swept to convergence, and the readout is a threshold — nothing in this file is longer than a paragraph. Start with the geometry, the one place the grid structure lives. Every pixel needs the sum of its four neighbors' means, and the honest way to get it for the whole grid at once is to pad the border with zeros and add the four shifted copies:

def neighbor_sum(mu):
    """Sum of each pixel's 4-connected neighbor means, as a full grid.

    Pad the grid with zeros so edge pixels simply have fewer neighbors (a
    free boundary), then add the four shifted copies: up, down, left, right.
    This is the only place the grid geometry lives — the pairwise coupling
    of the Markov random field is exactly "look at your four neighbors."
    """
    p = np.pad(mu, 1, mode="constant", constant_values=0.0)
    return p[:-2, 1:-1] + p[2:, 1:-1] + p[1:-1, :-2] + p[1:-1, 2:]


def energy(x, y, J, h):
    """The Ising / Boltzmann energy of a hard binary configuration x.

        E(x) = -J * sum over neighbor pairs (x_i x_j)  -  h * sum_i (x_i y_i)

    The first term is the smoothness prior: neighbors that agree (same sign)
    subtract J, so smooth images are low-energy. The second is the data term:
    a hidden pixel that matches its noisy observation subtracts h. Lower
    energy is more probable under p(x) ∝ exp(-E(x)). x and y are ±1 grids.
    """
    pair_term = 0.5 * float(np.sum(x * neighbor_sum(x)))  # each edge counted once
    data_term = float(np.sum(x * y))
    return -J * pair_term - h * data_term

That's two functions. neighbor_sum is the pairwise coupling of the Markov random field made concrete — "look up, down, left, right" — and energy scores a hard configuration with the Ising energy from the math, the smoothness term plus the data term. We never minimize energy directly; it's here because the free energy is built on it and because it's the thing the whole model is quietly about. The inference engine is a single line, the fixed-point update for one pixel:

def mean_field_pixel(nbr, y_i, J, h):
    """One pixel's mean-field fixed point: tanh of the incoming field.

    The field a pixel feels is J times the sum of its neighbors' means plus
    h times its own noisy observation. Passing that field through tanh gives
    the posterior mean the factorized q assigns to the pixel, in [-1, +1]:

        mu_i = tanh( J * nbr + h * y_i )

    Written with plain arithmetic, so `nbr` and `y_i` can be scalars (one
    pixel in a sweep) or whole grids. (In Bernoulli / sigmoid form this is the
    identical statement pi_i = sigmoid(2 * field) with pi_i = (1 + mu_i) / 2.)
    """
    return np.tanh(J * nbr + h * y_i)

tanh of the field, exactly as derived — $J$ times the neighbor sum plus $h$ times the observation. Feed it scalars and it updates one pixel; feed it grids and it updates all of them. Now the sweep, and here is the one design decision that matters. I update the pixels one at a time in raster order, each using its neighbors' means as they stand at that moment, so a pixel fixed early in the sweep immediately helps the ones after it:

def mean_field_sweep(mu, y, J, h):
    """One coordinate-descent sweep over the grid, one pixel at a time.

    Walk the pixels in raster order and update each from the means around it
    *as they stand right now* — a pixel updated early in the sweep already
    feeds the ones after it (Gauss-Seidel order). Each single-pixel tanh update
    is the exact coordinate-wise minimizer of the mean-field free energy, so
    sweeping this way drives F down monotonically. Updating every pixel at once
    from the previous sweep's means is cheaper but can oscillate instead of
    settle, so we take the honest coordinate descent. Returns the new grid.
    """
    mu = mu.astype(float).copy()
    rows, cols = mu.shape
    for i in range(rows):
        for j in range(cols):
            nbr = 0.0
            if i > 0:        nbr += mu[i - 1, j]
            if i + 1 < rows: nbr += mu[i + 1, j]
            if j > 0:        nbr += mu[i, j - 1]
            if j + 1 < cols: nbr += mu[i, j + 1]
            mu[i, j] = mean_field_pixel(nbr, y[i, j], J, h)
    return mu

This is Gauss-Seidel order, and it's not a detail — it's what makes the free energy fall monotonically. Updating every pixel at once from the previous sweep's means is tidier to write and it vectorizes, but it can overshoot and oscillate, two half-right answers flickering forever without settling. The one-at-a-time sweep is true coordinate descent, so it always goes downhill. To see that it does, we need the objective itself:

def free_energy(mu, y, J, h, eps=1e-12):
    """The mean-field free energy F(mu) that the sweeps minimize.

        F(mu) = E_q[E(x)]  -  H(q)

    Under the factorized q, E_q[x_i x_j] = mu_i mu_j, so the expected energy
    is the same Ising energy evaluated on the means. H(q) is the entropy of
    the independent Bernoullis. Minimizing F over mu is exactly minimizing
    KL(q || p) — the variational objective. Each tanh sweep is coordinate
    descent on this, so it decreases monotonically to a fixed point.
    """
    pair = 0.5 * float(np.sum(mu * neighbor_sum(mu)))
    expected_energy = -J * pair - h * float(np.sum(mu * y))
    p = np.clip((1.0 + mu) / 2.0, eps, 1.0 - eps)  # q(x_i = +1)
    entropy = float(-np.sum(p * np.log(p) + (1.0 - p) * np.log(1.0 - p)))
    return expected_energy - entropy

Straight from the math: the expected energy under the factorized $q$ (the Ising energy evaluated on the means) minus the entropy of the per-pixel coins. This is the number that must never rise, and computing it every sweep is how we know the loop is healthy — if it ever ticks up, the sweep has a bug. Reading a clean image back out is a threshold, the maximum-marginal estimate:

def to_binary(mu):
    """Threshold the means back to a clean ±1 image by their sign.

    A pixel whose posterior leans positive (mu_i >= 0) is read as +1, the
    rest as -1. This is the maximum-marginal-probability estimate of the
    clean image under the mean-field posterior.
    """
    return np.where(mu >= 0.0, 1, -1).astype(int)

A pixel whose mean leans positive is called ink, the rest background. Finally the loop that ties it together: start the means at the noisy image, sweep until they stop moving, and optionally record the whole history so we can replay it:

def denoise(y, J=1.0, h=1.1, max_sweeps=30, tol=1e-4, record=False):
    """Run mean-field inference to convergence, starting from the noisy image.

    Initialize the means at the observed pixels (mu = y), then sweep. Stop
    when the largest change in any mean drops below `tol`, or after
    `max_sweeps`. With record=True, also return the per-sweep history: the
    full grid of means, the free energy, and the pixel-agreement accuracy
    against the noisy input at each sweep, so the chapter can replay the
    denoising frame by frame.
    """
    mu = y.astype(float).copy()
    history = [mu.copy()] if record else None
    energies = [free_energy(mu, y, J, h)] if record else None
    for _ in range(max_sweeps):
        mu_new = mean_field_sweep(mu, y, J, h)
        delta = float(np.max(np.abs(mu_new - mu)))
        mu = mu_new
        if record:
            history.append(mu.copy())
            energies.append(free_energy(mu, y, J, h))
        if delta < tol:
            break
    if record:
        return mu, history, energies
    return mu


def accuracy(estimate, clean):
    """Fraction of pixels that match the clean image (both ±1 grids)."""
    return float((estimate == clean).mean())


def psnr(estimate, clean):
    """Peak signal-to-noise ratio between two ±1 images, in decibels.

    For a two-level image this is a monotone restatement of the pixel error
    rate — MSE on the ±1 scale is 4 times the fraction of wrong pixels — but
    it's the number image work usually quotes, so we report it too.
    """
    mse = float(np.mean((estimate.astype(float) - clean.astype(float)) ** 2))
    if mse == 0.0:
        return float("inf")
    peak = 2.0  # the ±1 range spans 2
    return float(10.0 * np.log10(peak ** 2 / mse))

Initialize $\mu = y$ — the noisy image is our first guess at the clean one — then sweep, checking after each pass how far the largest mean moved and stopping when it falls below tolerance. With record=True it also keeps the full grid of means and the free energy at every sweep, which is exactly what the animation plays back. The two scoring helpers, accuracy and PSNR, are for grading against the clean image we kept sealed.

Watch it work

This is the payoff. Below is a real run of the denoise function above on the noisy "3", one frame per sweep, starting from the noisy image itself. The top panel is the current grid of means $\mu$ as a heatmap — soft grays where the model is unsure, hardening to black and white as it commits. The bottom panel traces the free energy, the objective from the math, one point per sweep. Press play and watch the speckle get eaten.

The first sweep does almost everything. The image arrives with 96 wrong pixels — 83.3% correct — and a single coordinate-descent pass drops the free energy from -1,010 to -1,312 and the picture snaps from speckled to mostly clean, because every isolated flipped pixel is surrounded by neighbors that outvote it in one shot. From there it's cleanup along the edges: the second sweep reaches -1,349, and then a long tail of sub-integer improvements as the last few ambiguous pixels near the strokes make up their minds. By sweep 11 the free energy has settled at -1,351.7 and the image has stopped changing; the loop runs a couple more sweeps only to confirm the means have stopped moving, and converges after 13. The final estimate has 19 wrong pixels out of 576 — 96.7% correct, up from 83.3%.

Now do what the EM chapter asked you to do and scan the bottom panel for a single sweep where the energy goes up. There isn't one. That's not luck of this seed; it's the coordinate-descent guarantee from the math, and it's the reason you can trust the loop is improving even when the picture has stopped visibly changing. Here is that same descent on its own, every sweep:

A cliff on the first sweep, a knee, then flat. The same monotone shape as EM's log-likelihood climb, just pointing down instead of up because we're minimizing a free energy rather than maximizing a likelihood — and for the same reason: each step is the exact optimum of a variational bound in one coordinate, so the whole never moves the wrong way.

The full implementation

The entire file, no library calls in the algorithm itself, top to bottom. This is exactly what the animation ran:

"""A mean-field Boltzmann machine for binary image denoising, from scratch.

The model is an Ising-model / Markov random field over a grid of hidden
binary pixels x_i in {-1, +1}. Two things pull on each hidden pixel:

  * a smoothness term that couples it to its 4-connected neighbors (a pair
    of neighbors that agree lowers the energy), weight J > 0;
  * an observation term that couples it to the noisy pixel y_i we actually
    saw, strength h > 0.

Exact inference on this grid is intractable, so we approximate the posterior
over the clean image with a fully factorized distribution q(x) = prod_i q_i,
and summarize each q_i by its mean mu_i = E_q[x_i] in [-1, +1]. Mean-field
variational inference then reduces to a fixed-point update per pixel:

    mu_i <- tanh( J * (sum of neighbor means) + h * y_i )

Sweep that over every pixel until the means stop moving, threshold mu to a
sign, and you have a denoised image. That update is coordinate descent on the
mean-field free energy F(mu) = E_q[E(x)] - H(q); every sweep lowers F.

Pure NumPy. Every function below shows up in the chapter one step at a time;
the `# region:` markers are what the book's include directives pull in.
"""

import numpy as np


# region: energy
def neighbor_sum(mu):
    """Sum of each pixel's 4-connected neighbor means, as a full grid.

    Pad the grid with zeros so edge pixels simply have fewer neighbors (a
    free boundary), then add the four shifted copies: up, down, left, right.
    This is the only place the grid geometry lives — the pairwise coupling
    of the Markov random field is exactly "look at your four neighbors."
    """
    p = np.pad(mu, 1, mode="constant", constant_values=0.0)
    return p[:-2, 1:-1] + p[2:, 1:-1] + p[1:-1, :-2] + p[1:-1, 2:]


def energy(x, y, J, h):
    """The Ising / Boltzmann energy of a hard binary configuration x.

        E(x) = -J * sum over neighbor pairs (x_i x_j)  -  h * sum_i (x_i y_i)

    The first term is the smoothness prior: neighbors that agree (same sign)
    subtract J, so smooth images are low-energy. The second is the data term:
    a hidden pixel that matches its noisy observation subtracts h. Lower
    energy is more probable under p(x) ∝ exp(-E(x)). x and y are ±1 grids.
    """
    pair_term = 0.5 * float(np.sum(x * neighbor_sum(x)))  # each edge counted once
    data_term = float(np.sum(x * y))
    return -J * pair_term - h * data_term
# endregion


# region: update
def mean_field_pixel(nbr, y_i, J, h):
    """One pixel's mean-field fixed point: tanh of the incoming field.

    The field a pixel feels is J times the sum of its neighbors' means plus
    h times its own noisy observation. Passing that field through tanh gives
    the posterior mean the factorized q assigns to the pixel, in [-1, +1]:

        mu_i = tanh( J * nbr + h * y_i )

    Written with plain arithmetic, so `nbr` and `y_i` can be scalars (one
    pixel in a sweep) or whole grids. (In Bernoulli / sigmoid form this is the
    identical statement pi_i = sigmoid(2 * field) with pi_i = (1 + mu_i) / 2.)
    """
    return np.tanh(J * nbr + h * y_i)
# endregion


# region: sweep
def mean_field_sweep(mu, y, J, h):
    """One coordinate-descent sweep over the grid, one pixel at a time.

    Walk the pixels in raster order and update each from the means around it
    *as they stand right now* — a pixel updated early in the sweep already
    feeds the ones after it (Gauss-Seidel order). Each single-pixel tanh update
    is the exact coordinate-wise minimizer of the mean-field free energy, so
    sweeping this way drives F down monotonically. Updating every pixel at once
    from the previous sweep's means is cheaper but can oscillate instead of
    settle, so we take the honest coordinate descent. Returns the new grid.
    """
    mu = mu.astype(float).copy()
    rows, cols = mu.shape
    for i in range(rows):
        for j in range(cols):
            nbr = 0.0
            if i > 0:        nbr += mu[i - 1, j]
            if i + 1 < rows: nbr += mu[i + 1, j]
            if j > 0:        nbr += mu[i, j - 1]
            if j + 1 < cols: nbr += mu[i, j + 1]
            mu[i, j] = mean_field_pixel(nbr, y[i, j], J, h)
    return mu
# endregion


# region: free_energy
def free_energy(mu, y, J, h, eps=1e-12):
    """The mean-field free energy F(mu) that the sweeps minimize.

        F(mu) = E_q[E(x)]  -  H(q)

    Under the factorized q, E_q[x_i x_j] = mu_i mu_j, so the expected energy
    is the same Ising energy evaluated on the means. H(q) is the entropy of
    the independent Bernoullis. Minimizing F over mu is exactly minimizing
    KL(q || p) — the variational objective. Each tanh sweep is coordinate
    descent on this, so it decreases monotonically to a fixed point.
    """
    pair = 0.5 * float(np.sum(mu * neighbor_sum(mu)))
    expected_energy = -J * pair - h * float(np.sum(mu * y))
    p = np.clip((1.0 + mu) / 2.0, eps, 1.0 - eps)  # q(x_i = +1)
    entropy = float(-np.sum(p * np.log(p) + (1.0 - p) * np.log(1.0 - p)))
    return expected_energy - entropy
# endregion


# region: threshold
def to_binary(mu):
    """Threshold the means back to a clean ±1 image by their sign.

    A pixel whose posterior leans positive (mu_i >= 0) is read as +1, the
    rest as -1. This is the maximum-marginal-probability estimate of the
    clean image under the mean-field posterior.
    """
    return np.where(mu >= 0.0, 1, -1).astype(int)
# endregion


# region: denoise
def denoise(y, J=1.0, h=1.1, max_sweeps=30, tol=1e-4, record=False):
    """Run mean-field inference to convergence, starting from the noisy image.

    Initialize the means at the observed pixels (mu = y), then sweep. Stop
    when the largest change in any mean drops below `tol`, or after
    `max_sweeps`. With record=True, also return the per-sweep history: the
    full grid of means, the free energy, and the pixel-agreement accuracy
    against the noisy input at each sweep, so the chapter can replay the
    denoising frame by frame.
    """
    mu = y.astype(float).copy()
    history = [mu.copy()] if record else None
    energies = [free_energy(mu, y, J, h)] if record else None
    for _ in range(max_sweeps):
        mu_new = mean_field_sweep(mu, y, J, h)
        delta = float(np.max(np.abs(mu_new - mu)))
        mu = mu_new
        if record:
            history.append(mu.copy())
            energies.append(free_energy(mu, y, J, h))
        if delta < tol:
            break
    if record:
        return mu, history, energies
    return mu


def accuracy(estimate, clean):
    """Fraction of pixels that match the clean image (both ±1 grids)."""
    return float((estimate == clean).mean())


def psnr(estimate, clean):
    """Peak signal-to-noise ratio between two ±1 images, in decibels.

    For a two-level image this is a monotone restatement of the pixel error
    rate — MSE on the ±1 scale is 4 times the fraction of wrong pixels — but
    it's the number image work usually quotes, so we report it too.
    """
    mse = float(np.mean((estimate.astype(float) - clean.astype(float)) ** 2))
    if mse == 0.0:
        return float("inf")
    peak = 2.0  # the ±1 range spans 2
    return float(10.0 * np.log10(peak ** 2 / mse))
# endregion

The library version

There is no drop-in scikit-learn twin for this. A mean-field MRF denoiser isn't a named estimator you can import and .fit(); it's a model you assemble. The nearest named relative in scikit-learn is sklearn.neural_network.BernoulliRBM, a restricted Boltzmann machine, and it's worth knowing why it's a cousin and not a substitute. An RBM is the same Boltzmann/energy family — binary units, an energy, a p(x) ∝ exp(-E) — but its couplings are learned from a corpus of images by contrastive divergence, where ours are fixed by hand at $J = 1$. It's a feature learner that discovers structure across many images; it is not a tool for cleaning up one image against a smoothness prior. Point it at our single noisy "3" and there's nothing for it to learn. So it isn't the thing to benchmark against.

The fair comparison is against the denoiser people actually reach for on a binary image: a 3×3 median filter. For a two-level image the median of nine pixels is just their majority vote, so it's the "make each pixel copy whatever most of its little box says" rule — local, parameter-free, and completely blind to any structure past a single pixel's immediate neighborhood. scipy ships it, and that's the whole baseline:

def median_baseline(noisy):
    """Denoise a +/-1 image with a 3x3 median filter (majority vote).

    Slide a 3x3 window over the grid and replace each pixel with the median of
    the nine values under it. On a two-level image that median is the majority
    sign in the window, so an isolated flipped pixel surrounded by agreement
    gets overruled. It's the honest cheap baseline: local, parameter-free, and
    blind to any global structure beyond a single pixel's immediate box. Returns
    a +/-1 grid.
    """
    filtered = median_filter(noisy.astype(int), size=3, mode="nearest")
    return np.where(filtered >= 0, 1, -1).astype(int)

It's a good, honest baseline precisely because it encodes the same instinct as our model — trust the neighborhood — but with none of the global bookkeeping. Both the median filter and the mean-field MRF believe a pixel should look like its surroundings. The difference is that the median filter decides each pixel once, in isolation, from a fixed 3×3 window, while the MRF lets the decisions talk to each other and iterate to a joint fixed point. The next section is about how much that difference is worth.

Scratch versus library

Three images, graded against the clean original: the noisy input, the median filter's attempt, and the mean-field reconstruction. First the pictures, so you can see what each one did, then the numbers.

And the scores, accuracy and PSNR side by side for all three:

The mean-field model wins, and it's not close. The noisy image starts at 83.3% pixel accuracy and 7.78 dB PSNR. The median filter cleans it to 92.0% and 10.98 dB — a real improvement, it kills most of the isolated speckle exactly as advertised. But the mean-field MRF reaches 96.7% and 14.82 dB, cutting the wrong-pixel count from the noisy image's 96 down to 19, against the median filter's 46. Nearly four decibels of PSNR separate them. The reason is visible in the pictures: where two or three noise pixels happen to land next to each other, the median filter's fixed 3×3 window sees a local majority of noise and preserves the blob, while the MRF's iteration lets the surrounding clean region keep pressing inward over several sweeps until the blob collapses. Global structure, resolved by iteration, beats a one-shot local vote.

I'll be honest about the ceiling, though. Both methods lose the same battle in the same place — the thin parts of the stroke and the sharp inside corners of the "3," where the true image has detail one pixel wide. The smoothness prior can't tell a one-pixel stroke from a one-pixel piece of noise; both look like a pixel disagreeing with its neighbors, and the model rounds both off. That's not a tuning failure, it's the independence assumption showing through: mean field is structurally blind to the coordinated, correlated detail that fine structure is made of. It wins decisively on the speckle in the flat regions and gives back a little on the edges, and 96.7% is where that trade settles for this image at $J = 1$, $h = 1.1$.

Takeaways

Reach for a Boltzmann/MRF model when you can phrase what you know as a local, pairwise preference — neighbors should agree, adjacent labels should be consistent, this should look like that — and you want a principled way to fold that prior together with noisy evidence. Writing it as an energy gives you a probability model; mean field gives you a way to actually query that model without paying for exact inference. The whole apparatus on this page is one energy function and a tanh swept to convergence, and it beat a strong local baseline by four decibels with no training data at all. When your structure is genuinely local, that is a lot of mileage from very little machinery.

Keep the two caveats in your pocket, because they're the same two every time. First, mean field is an approximation that assumes away the very correlations your prior is about, so it's overconfident and it smooths away fine detail — when you need the correlations, you reach for the heavier tools mean field is the cheap alternative to (loopy belief propagation, Gibbs sampling, or just a bigger model). Second, it's coordinate descent to a local optimum of the free energy, cousin to EM in exactly this respect: the monotone descent is guaranteed, the global optimum is not, and where you land depends on where you start and how you sweep.

These ideas did not stay in the 1980s. The energy-based framing — score configurations, turn scores into probabilities, do approximate inference — is the backbone of conditional random fields, which were the standard for pixel-level image segmentation and sequence labeling for a decade, and mean-field inference for a dense CRF was for years the refinement step bolted onto the output of segmentation networks. The variational bound we minimized here, KL from a simple $q$ to an intractable $p$, is the same evidence lower bound that trains variational autoencoders. And the Boltzmann distribution over a corrupted image, walked back toward a clean sample, is precisely the picture behind modern diffusion models — a learned energy landscape, denoised in steps. Build the 576-pixel version by hand once, watch the free energy fall, and you've built the intuition for the whole energy-based family that came after it.