Máquinas de Boltzmann de campo medio

Qué cubre este capítulo

Aquí tienes una imagen con el 18% de sus pixeles volteados al color equivocado, y aquí está la creencia que la arregla: un pixel probablemente se parece a sus vecinos. Ese único prior, escrito como una energía y llevado a su conclusión lógica, basta para regresar una imagen binaria ruidosa hacia la versión limpia de la que salió — sin conjunto de entrenamiento, sin etiquetas, sin descenso de gradiente sobre una pérdida. Este capítulo construye el modelo que lo logra: una máquina de Boltzmann sobre los pixeles de una imagen, el caballo de batalla de los modelos basados en energía, y la resuelve con inferencia variacional de campo medio, el truco barato de inferencia aproximada que convierte un posterior intratable en un puñado de llamadas a tanh.

El planteamiento es un modelo de Ising tomado directamente de la física. Cada pixel es un pequeño imán que puede apuntar hacia arriba o hacia abajo; los vecinos quieren coincidir, y cada pixel además siente un jalón hacia lo que la imagen ruidosa dice que debería ser. Esas dos fuerzas se negocian entre sí, y la configuración que las equilibra es una imagen limpia. Calcular la respuesta exacta es imposible — hay $2^{576}$ configuraciones de una cuadrícula de 24×24 — así que no lo hacemos. Aproximamos todo el posterior con un número por pixel, su media posterior, y actualizamos esos números uno a la vez hasta que dejan de moverse. Ese ciclo es todo el método, y vas a poder verlo correr: la animación a la mitad de esta página es un denoising real, barrido por barrido, con el ruido perdiendo la votación frente al consenso mientras un panel de energía libre se desliza cuesta abajo a un lado.

A lo que hay que ponerle el ojo es a esa energía libre. Igual que la log-verosimilitud en el capítulo de EM, es el número que certifica que el ciclo está funcionando — solo puede bajar, en cada barrido, y cuando se aplana la imagen está lista. EM y esto son primos cercanos; ambos son descenso por coordenadas sobre una cota variacional, y si aprendiste a confiar en la subida monótona de EM, ya sabes leer este descenso. Al final ponemos el resultado a competir contra el denoiser al que todo mundo recurre en la práctica, un filtro de mediana, y te voy a mostrar exactamente dónde el modelo de energía se gana su lugar y dónde no.

Un poco de historia

El modelo es más viejo que el machine learning. En 1920 Wilhelm Lenz le entregó a su estudiante Ernst Ising un juguete: una retícula de espines, cada uno ±1, cada uno prefiriendo alinearse con sus vecinos, y la pregunta de si esa preferencia local podía producir orden global — magnetismo. Ising resolvió el caso unidimensional en su tesis de 1925 y concluyó, erróneamente para dimensiones superiores, que nunca se magnetizaba. El modelo se quedó con su nombre de todos modos, y se convirtió en la mosca de la fruta de la física estadística: el sistema más simple donde las reglas locales suman comportamiento colectivo.

El salto de la física a la computación llegó a principios de los ochenta. El artículo de John Hopfield de 1982 mostró que una red de estas unidades ±1, conectada simétricamente, tiene una energía que solo puede disminuir conforme las unidades se voltean — así que la red se asienta en mínimos que puedes usar como memorias almacenadas. Dos años después, Geoffrey Hinton y Terry Sejnowski convirtieron la red determinista de Hopfield en una probabilística: su máquina de Boltzmann, bautizada en su artículo de 1985 con David Ackley, reemplazó el volteo duro por uno estocástico gobernado por la distribución de Boltzmann, lo que hizo de las unidades muestras de un modelo de probabilidad real y le dio a la red una regla de aprendizaje. Ese es el linaje en el que se sienta este capítulo — una energía sobre unidades binarias, leída como una probabilidad.

La otra mitad de la historia es el truco de inferencia. Muestrear una máquina de Boltzmann hasta la convergencia es lento, y en 1987 Carsten Peterson y James Anderson propusieron saltárselo: reemplaza el estado fluctuante de cada unidad por su promedio, y resuelve para los promedios directamente. Esa es la aproximación de campo medio, tomada tal cual de la teoría de Weiss del magnetismo, y la ecuación de punto fijo que escribieron es la actualización con tanh que construimos abajo. En paralelo, el artículo de 1984 de Stuart y Donald Geman puso exactamente este tipo de modelo a trabajar sobre imágenes — un prior de campo aleatorio de Markov para restauración de imágenes — y para finales de los noventa el programa de inferencia variacional de Michael Jordan, Zoubin Ghahramani y sus colaboradores había replanteado el campo medio como una instancia de una idea general: aproximar un posterior difícil con una familia simple y minimizar la distancia hacia él. Cada pieza de esa historia está en el código de esta página.

La intuición

Piensa en la imagen limpia como la verdad y en la imagen ruidosa como un testigo corrupto. El ruido de sal y pimienta voltea pixeles dispersos al color equivocado, así que el testigo dice la verdad la mayor parte del tiempo pero miente al azar sobre más o menos uno de cada cinco pixeles. Quieres recuperar la verdad, y tienes exactamente una pieza de conocimiento externo: las imágenes reales son suaves. La tinta se sienta junto a la tinta, el fondo junto al fondo. Un pixel negro solitario varado en un campo blanco es casi seguro una mentira.

Así que arma un estira y afloja. Dale a cada pixel dos razones para ser de un color o del otro. La primera es el acuerdo: un pixel que coincide con sus cuatro vecinos recibe recompensa, así que las regiones suaves salen baratas y los volteos aislados salen caros. La segunda es la lealtad al testigo: un pixel que coincide con lo que reportó la imagen ruidosa también recibe recompensa, para que la reconstrucción no pueda irse de pinta hacia una página en blanco. Ajusta las dos recompensas una contra la otra y la resolución es automática. Un pixel volteado rodeado de cuatro vecinos que no están de acuerdo con él siente cuatro unidades de presión para cambiar y solo una unidad de lealtad que lo mantiene en su lugar, así que cambia — el ruido pierde la votación. Un pixel en un borde real, donde los vecinos genuinamente discrepan, siente un voto dividido y le cede la decisión al testigo. El prior limpia los pixeles fáciles y no se mete con los difíciles.

Esa es toda la idea, y vale la pena ver la entrada antes que las matemáticas. A la izquierda está la imagen limpia de 24×24 — un dígito, binarizado a tinta negra sobre fondo blanco. A la derecha está la misma imagen después de que volteé el 18% de sus pixeles al azar. Ese revoltijo moteado de la derecha es todo lo que el algoritmo llega a ver.

Todavía puedes leer el dígito a través del moteado, y ese es justo el punto: tu sistema visual está haciendo inferencia de campo medio gratis, dejando que la mayoría intacta le gane la votación a las mentiras dispersas. Estamos a punto de escribir eso como aritmética.

Las matemáticas

Codifica la imagen como espines. Cada pixel $i$ recibe un valor oculto $x_i \in \{-1, +1\}$ — piensa $+1$ para tinta, $-1$ para fondo — y escribimos $y_i \in \{-1, +1\}$ para el pixel ruidoso que sí observamos. El modelo califica una configuración completa $x$ con una energía, baja para las configuraciones en las que creemos:

La primera suma recorre los pares de pixeles vecinos $\mathcal{E}$ (cada pixel unido a los cuatro que lo rodean); el producto $x_i x_j$ es $+1$ cuando dos vecinos coinciden y $-1$ cuando chocan, así que con $J > 0$ el acuerdo baja la energía — el prior de suavidad. La segunda suma es el término de datos: $x_i y_i$ es $+1$ cuando un pixel coincide con su observación ruidosa, así que con $h > 0$ coincidir con el testigo también baja la energía. $J$ fija cuánto confiamos en la suavidad, $h$ cuánto confiamos en la imagen ruidosa. Convierte la energía en una probabilidad con la distribución de Boltzmann, donde energía baja significa probabilidad alta:

Esa $Z$ es la trampa. Suma sobre las $2^{576}$ configuraciones de la cuadrícula, así que no podemos ni calcular $p(x)$ ni muestrear de ella directamente, y la marginal que de verdad queremos — la probabilidad de que un pixel dado sea tinta — está enterrada adentro. Aquí es donde el campo medio se gana su lugar. En vez del posterior verdadero usamos un sustituto $q$ que finge que los pixeles son independientes, un producto de distribuciones de un pixel:

Cada factor $q_i$ es una moneda sobre $\{-1, +1\}$, y queda completamente descrito por su media $\mu_i = \mathbb{E}_q[x_i] \in [-1, +1]$ — un pixel suave, de valor real, $+1$ para "seguro es tinta", $-1$ para "seguro es fondo", $0$ para "ni idea". Ajustar $q$ a $p$ significa elegir las $\mu_i$ que hacen a las dos distribuciones lo más cercanas posible en divergencia KL, y hacer ese cálculo para un pixel con los demás fijos da una ecuación de punto fijo limpia. La media de cada pixel es una tanh del campo total que lo presiona:

Léela directo de la energía: el campo sobre el pixel $i$ es $J$ por la suma de las medias actuales de sus vecinos más $h$ por su propia observación, y tanh aplasta ese campo a una media en $[-1, +1]$. Barre esto sobre cada pixel, una y otra vez, y las medias se asientan. ¿Al fondo de qué se están asentando? De la energía libre de campo medio, el objetivo que todo el ciclo minimiza:

Bajo la $q$ independiente, la energía esperada es simplemente la energía de Ising evaluada en las medias (porque $\mathbb{E}_q[x_i x_j] = \mu_i \mu_j$), y $H(q)$ es la entropía sumada de las monedas por pixel. La razón por la que $F$ es lo correcto a minimizar es que difiere del logaritmo de esa $Z$ intratable por exactamente el error de aproximación:

Como $\log Z$ es una constante que no podemos calcular pero tampoco necesitamos, empujar $F$ hacia abajo empuja $\mathrm{KL}(q \,\|\, p)$ hacia abajo — acerca nuestra aproximación factorizada al posterior verdadero tanto como lo permite el supuesto de independencia. Y cada actualización tanh de un solo pixel es el minimizador exacto de $F$ en esa coordenada, así que un barrido de ellas solo puede bajar $F$. Esa es la misma garantía de descenso-por-coordenadas-sobre-una-cota que hace confiable a EM, y la animación de abajo te hace verla cumplirse.

En qué es bueno, en qué no

La fortaleza es que el campo medio vuelve barato un cálculo imposible y lo convierte en algo que puedes leer. El posterior exacto sobre una imagen de 24×24 vive en un espacio de $2^{576}$ configuraciones con un normalizador que nadie puede sumar; el campo medio lo reemplaza con 576 números y una iteración de punto fijo que no es más que promedios de vecinos pasados por tanh. No hay fase de entrenamiento, no hay datos más allá de la única imagen, no hay learning rate, no hay paso que pueda divergir — solo un deslizamiento monótono cuesta abajo por la energía libre hasta un punto fijo. Cuando tu prior de verdad es local y por pares, que es exactamente el caso de la suavidad en una imagen, esta es más o menos la inferencia más barata que todavía respeta la estructura de toda la cuadrícula, y se paraleliza y escala como lo hace una convolución.

La debilidad es el supuesto que compró la velocidad. El campo medio declara a los pixeles independientes, y enfáticamente no lo son — ese es todo el contenido de un prior de suavidad. Así que la aproximación es sistemáticamente sobreconfiada: empuja cada $\mu_i$ hacia un $\pm 1$ duro más rápido de lo que lo haría el posterior verdadero, y no puede representar las correlaciones que cargan la estructura fina de una imagen. En un trazo delgado o una esquina afilada, donde la respuesta correcta depende de que los vecinos coincidan de manera coordinada — algo que la $q$ factorizada no puede expresar — el campo medio suaviza donde no debería y redondea la esquina. Además encuentra un mínimo local de $F$, no el global, así que el punto fijo depende de dónde empiezas y del orden en que barres. Y las dos perillas $J$ y $h$ te toca fijarlas a ti: sube $J$ demasiado y el modelo derrite toda la imagen en una mancha plana; bájala demasiado y copia el ruido tal cual. Es una aproximación rápida y honesta, no la respuesta exacta, y los lugares donde falla son los lugares donde los pixeles se estaban hablando a sus espaldas.

Los datos

Una imagen binaria de 24×24, y su gemela ruidosa. Construí la imagen limpia a partir de load_digits de scikit-learn — un "3" escrito a mano de 8×8 — umbralizándola a blanco y negro y luego escalando cada pixel a un bloque de 3×3 para obtener una cuadrícula de 24×24, 576 pixeles, 189 de ellos tinta. Eso da una forma con estructura genuina que proteger: trazos, curvas, un hueco en medio, el tipo de cosa con la que un prior de suavidad puede ayudar y también el tipo de cosa que puede sobre-suavizar. Luego la corrupción: sembré un generador aleatorio y volteé el 18% de los pixeles — 96 de los 576 — de su color verdadero al opuesto, ruido de sal y pimienta. La cuadrícula limpia y la ruidosa están ambas versionadas como CSV bajo data/, y la semilla está fija, así que cada número de esta página se reproduce exactamente. Al algoritmo se le entrega solo la cuadrícula ruidosa; la limpia existe únicamente para poder calificar el resultado al final.

Constrúyelo, una función a la vez

Seis funciones cortas, y cuatro de ellas son de tres líneas. El modelo es la energía, la inferencia es la actualización tanh barrida hasta la convergencia, y la lectura final es un umbral — nada en este archivo es más largo que un párrafo. Empieza con la geometría, el único lugar donde vive la estructura de la cuadrícula. Cada pixel necesita la suma de las medias de sus cuatro vecinos, y la manera honesta de obtenerla para toda la cuadrícula de una vez es rellenar el borde con ceros y sumar las cuatro copias desplazadas:

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

Esas son dos funciones. neighbor_sum es el acoplamiento por pares del campo aleatorio de Markov hecho concreto — "mira arriba, abajo, izquierda, derecha" — y energy califica una configuración dura con la energía de Ising de las matemáticas, el término de suavidad más el término de datos. Nunca minimizamos energy directamente; está aquí porque la energía libre se construye sobre ella y porque es aquello de lo que, en el fondo, trata todo el modelo. El motor de inferencia es una sola línea, la actualización de punto fijo para un 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 del campo, exactamente como se derivó — $J$ por la suma de vecinos más $h$ por la observación. Dale escalares y actualiza un pixel; dale cuadrículas y los actualiza todos. Ahora el barrido, y aquí está la única decisión de diseño que importa. Actualizo los pixeles uno a la vez en orden de rasterizado, cada uno usando las medias de sus vecinos tal como están en ese momento, de modo que un pixel corregido temprano en el barrido ayuda de inmediato a los que vienen después:

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

Esto es orden Gauss-Seidel, y no es un detalle — es lo que hace que la energía libre caiga monótonamente. Actualizar todos los pixeles a la vez con las medias del barrido anterior es más limpio de escribir y se vectoriza, pero puede pasarse de largo y oscilar, dos respuestas medio correctas parpadeando para siempre sin asentarse. El barrido de uno-a-la-vez es descenso por coordenadas de verdad, así que siempre va cuesta abajo. Para comprobar que así es, necesitamos el objetivo mismo:

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

Directo de las matemáticas: la energía esperada bajo la $q$ factorizada (la energía de Ising evaluada en las medias) menos la entropía de las monedas por pixel. Este es el número que nunca debe subir, y calcularlo en cada barrido es cómo sabemos que el ciclo está sano — si alguna vez sube un poco, el barrido tiene un bug. Leer una imagen limpia de vuelta es un umbral, la estimación de máxima marginal:

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)

Un pixel cuya media se inclina a positivo se declara tinta, el resto fondo. Finalmente el ciclo que lo amarra todo: arranca las medias en la imagen ruidosa, barre hasta que dejen de moverse, y opcionalmente graba toda la historia para poder reproducirla:

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))

Inicializa $\mu = y$ — la imagen ruidosa es nuestra primera conjetura de la limpia — luego barre, revisando después de cada pasada qué tanto se movió la media que más cambió y deteniéndote cuando cae por debajo de la tolerancia. Con record=True también guarda la cuadrícula completa de medias y la energía libre en cada barrido, que es exactamente lo que la animación reproduce. Los dos auxiliares de calificación, accuracy y PSNR, son para evaluar contra la imagen limpia que mantuvimos sellada.

Míralo trabajar

Este es el premio. Abajo hay una corrida real de la función denoise de arriba sobre el "3" ruidoso, un cuadro por barrido, empezando desde la imagen ruidosa misma. El panel superior es la cuadrícula actual de medias $\mu$ como un heatmap — grises suaves donde el modelo está indeciso, endureciéndose a blanco y negro conforme se compromete. El panel inferior traza la energía libre, el objetivo de las matemáticas, un punto por barrido. Dale play y mira cómo se come el moteado.

El primer barrido hace casi todo. La imagen llega con 96 pixeles equivocados — 83.3% correcta — y una sola pasada de descenso por coordenadas baja la energía libre de -1,010 a -1,312 y la imagen brinca de moteada a mayormente limpia, porque cada pixel volteado aislado está rodeado de vecinos que le ganan la votación de un solo golpe. De ahí en adelante es limpieza a lo largo de los bordes: el segundo barrido llega a -1,349, y luego una larga cola de mejoras menores a un entero mientras los últimos pixeles ambiguos cerca de los trazos se deciden. Para el barrido 11 la energía libre se ha asentado en -1,351.7 y la imagen ha dejado de cambiar; el ciclo corre un par de barridos más solo para confirmar que las medias dejaron de moverse, y converge después de 13. La estimación final tiene 19 pixeles equivocados de 576 — 96.7% correcta, desde el 83.3% inicial.

Ahora haz lo que el capítulo de EM te pidió hacer y escanea el panel inferior en busca de un solo barrido donde la energía suba. No hay ninguno. No es suerte de esta semilla; es la garantía de descenso por coordenadas de las matemáticas, y es la razón por la que puedes confiar en que el ciclo está mejorando incluso cuando la imagen ya dejó de cambiar visiblemente. Aquí está ese mismo descenso por su cuenta, cada barrido:

Un acantilado en el primer barrido, una rodilla, y luego plano. La misma forma monótona que la subida de la log-verosimilitud de EM, solo que apuntando hacia abajo en vez de hacia arriba porque estamos minimizando una energía libre en lugar de maximizar una verosimilitud — y por la misma razón: cada paso es el óptimo exacto de una cota variacional en una coordenada, así que el conjunto nunca se mueve en la dirección equivocada.

La implementación completa

El archivo entero, sin llamadas a librerías en el algoritmo mismo, de arriba a abajo. Esto es exactamente lo que corrió la animación:

"""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

La versión de librería

No hay un gemelo directo en scikit-learn para esto. Un denoiser MRF de campo medio no es un estimador con nombre que puedas importar y hacerle .fit(); es un modelo que ensamblas. El pariente con nombre más cercano en scikit-learn es sklearn.neural_network.BernoulliRBM, una máquina de Boltzmann restringida, y vale la pena saber por qué es un primo y no un sustituto. Un RBM es la misma familia Boltzmann/energía — unidades binarias, una energía, un p(x) ∝ exp(-E) — pero sus acoplamientos se aprenden de un corpus de imágenes por divergencia contrastiva, mientras que los nuestros están fijos a mano en $J = 1$. Es un aprendiz de features que descubre estructura a través de muchas imágenes; no es una herramienta para limpiar una sola imagen contra un prior de suavidad. Apúntalo a nuestro único "3" ruidoso y no tiene nada que aprender. Así que no es contra lo que hay que hacer el benchmark.

La comparación justa es contra el denoiser al que la gente de verdad recurre en una imagen binaria: un filtro de mediana de 3×3. Para una imagen de dos niveles, la mediana de nueve pixeles es simplemente su voto de mayoría, así que es la regla de "haz que cada pixel copie lo que diga la mayoría de su cajita" — local, sin parámetros, y completamente ciega a cualquier estructura más allá del vecindario inmediato de un pixel. scipy lo trae incluido, y ese es todo el 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)

Es un baseline bueno y honesto precisamente porque codifica el mismo instinto que nuestro modelo — confía en el vecindario — pero sin nada de la contabilidad global. Tanto el filtro de mediana como el MRF de campo medio creen que un pixel debería parecerse a sus alrededores. La diferencia es que el filtro de mediana decide cada pixel una sola vez, en aislamiento, desde una ventana fija de 3×3, mientras que el MRF deja que las decisiones se hablen entre sí e iteren hasta un punto fijo conjunto. La siguiente sección trata de cuánto vale esa diferencia.

Desde cero contra librería

Tres imágenes, calificadas contra el original limpio: la entrada ruidosa, el intento del filtro de mediana, y la reconstrucción de campo medio. Primero las imágenes, para que veas qué hizo cada una, luego los números.

Y las calificaciones, accuracy y PSNR lado a lado para las tres:

El modelo de campo medio gana, y no es por poquito. La imagen ruidosa arranca en 83.3% de accuracy por pixel y 7.78 dB de PSNR. El filtro de mediana la limpia a 92.0% y 10.98 dB — una mejora real, mata la mayor parte del moteado aislado exactamente como promete. Pero el MRF de campo medio llega a 96.7% y 14.82 dB, recortando el conteo de pixeles equivocados de los 96 de la imagen ruidosa a 19, contra los 46 del filtro de mediana. Casi cuatro decibeles de PSNR los separan. La razón se ve en las imágenes: donde dos o tres pixeles de ruido caen uno junto al otro, la ventana fija de 3×3 del filtro de mediana ve una mayoría local de ruido y preserva la mancha, mientras que la iteración del MRF deja que la región limpia circundante siga presionando hacia adentro durante varios barridos hasta que la mancha colapsa. La estructura global, resuelta por iteración, le gana a un voto local de un solo tiro.

Eso sí, voy a ser honesto sobre el techo. Ambos métodos pierden la misma batalla en el mismo lugar — las partes delgadas del trazo y las esquinas interiores afiladas del "3", donde la imagen verdadera tiene detalle de un pixel de ancho. El prior de suavidad no puede distinguir un trazo de un pixel de una pieza de ruido de un pixel; ambos se ven como un pixel en desacuerdo con sus vecinos, y el modelo redondea a los dos. Eso no es una falla de ajuste, es el supuesto de independencia asomándose: el campo medio es estructuralmente ciego al detalle coordinado y correlacionado del que está hecha la estructura fina. Gana de manera contundente con el moteado en las regiones planas y cede un poco en los bordes, y 96.7% es donde ese intercambio se asienta para esta imagen con $J = 1$, $h = 1.1$.

Conclusiones

Recurre a un modelo Boltzmann/MRF cuando puedas expresar lo que sabes como una preferencia local por pares — los vecinos deberían coincidir, las etiquetas adyacentes deberían ser consistentes, esto debería parecerse a aquello — y quieras una manera con principios de combinar ese prior con evidencia ruidosa. Escribirlo como una energía te da un modelo de probabilidad; el campo medio te da una manera de consultar realmente ese modelo sin pagar por inferencia exacta. Todo el aparato de esta página es una función de energía y una tanh barrida hasta la convergencia, y le ganó a un baseline local fuerte por cuatro decibeles sin ningún dato de entrenamiento. Cuando tu estructura es genuinamente local, eso es mucho kilometraje con muy poca maquinaria.

Guarda los dos avisos en el bolsillo, porque siempre son los mismos dos. Primero, el campo medio es una aproximación que asume que no existen justo las correlaciones de las que trata tu prior, así que es sobreconfiado y suaviza el detalle fino — cuando necesitas las correlaciones, echas mano de las herramientas más pesadas de las que el campo medio es la alternativa barata (propagación de creencias con ciclos, muestreo de Gibbs, o de plano un modelo más grande). Segundo, es descenso por coordenadas hacia un óptimo local de la energía libre, primo de EM exactamente en este aspecto: el descenso monótono está garantizado, el óptimo global no, y dónde aterrizas depende de dónde empiezas y de cómo barres.

Estas ideas no se quedaron en los ochenta. El marco basado en energía — califica configuraciones, convierte las calificaciones en probabilidades, haz inferencia aproximada — es la columna vertebral de los campos aleatorios condicionales, que fueron el estándar para segmentación de imágenes a nivel de pixel y etiquetado de secuencias durante una década, y la inferencia de campo medio para un CRF denso fue por años el paso de refinamiento atornillado a la salida de las redes de segmentación. La cota variacional que minimizamos aquí, KL desde una $q$ simple hacia una $p$ intratable, es la misma evidence lower bound que entrena a los autoencoders variacionales. Y la distribución de Boltzmann sobre una imagen corrupta, caminada de regreso hacia una muestra limpia, es precisamente la imagen detrás de los modelos de difusión modernos — un paisaje de energía aprendido, denoiseado por pasos. Construye la versión de 576 pixeles a mano una vez, mira caer la energía libre, y habrás construido la intuición para toda la familia basada en energía que vino después.