Information Theory for Computer Vision

Information theory measures uncertainty and dependence in bits. In computer vision it supplies the cross-entropy loss that trains classifiers, the entropy bound behind image compression, the mutual information that aligns a CT scan with an MR image, the divergences that compare histograms and fit generative models, and the information gain that chooses a robot’s next view. This article treats the quantities as one system: entropy, joint and conditional entropy, mutual information (MI), differential entropy and the Kullback–Leibler (KL) divergence, also called relative entropy or information gain. It covers how they relate, why minimizing the KL divergence in one direction or the other gives different fits, and what the data processing inequality and Fano’s inequality say about what processing an image can and cannot do. It then maps where each quantity is used in vision.

Definition

Information theory assigns numbers to probability distributions. The entropy $H(X)$ of a discrete random variable is the average surprise $-\log p(x)$ of its outcomes. The conditional entropy $H(X \mid Y)$ is the uncertainty left about $X$ once $Y$ is known. The mutual information $I(X; Y)$ is how much knowing one variable reduces uncertainty about the other. The KL divergence $\mathrm{KL}(p \,\Vert\, q)$ measures how far a distribution $q$ is from a distribution $p$ . For continuous variables, the differential entropy $h(X)$ replaces the sum with an integral over a density. It behaves differently from $H$ : it can be negative, and it changes when the variable is rescaled.

The base of the logarithm sets the unit. Base 2 gives bits and base $e$ gives nats (Murphy, 2022, §6.1.1). Every quantity in this article is in bits unless it is marked otherwise.

Intuition

Pick a pixel at random from an image of a dark background (three quarters of the area) and a bright object (one quarter). Before you look at the pixel, its region label $X$ has entropy $H(X) = -0.75 \log_2 0.75 - 0.25 \log_2 0.25 \approx 0.81$ bits, a little less than the one bit of a fair coin, because the background is more likely. Now read the pixel’s grey level $Y$ . If the two regions’ grey levels never overlap, $Y$ tells you the label exactly. Then $H(X \mid Y) = 0$ , and $Y$ carries all 0.81 bits about $X$ : $I(X; Y) = H(X)$ . If noise makes the grey levels overlap, some uncertainty remains, and $I(X; Y)$ falls below $H(X)$ . Mutual information is the part of the label’s uncertainty that the measurement removes.

Processing cannot recover lost information. Quantizing the grey level to three bits merges values. Whatever the quantized value tells about the label, the original value told at least as much. This is the data processing inequality. A blur behaves differently, because it reads more than the pixel. It averages the pixel’s neighbours, which usually share its label and carry independent noise. So the blurred value can tell more about the label than the pixel’s own value does. That does not break the inequality: the inequality concerns functions of what you observed, and the blur observed nine pixels. The example below computes all of these numbers exactly.

Formal Definition

Entropy, joint and conditional entropy

For a discrete $X$ with distribution $p$ over $K$ states, $H(X) = -\sum_k p(k) \log p(k)$ . It is largest for the uniform distribution, where $H = \log K$ , and zero for a distribution concentrated on one state (Murphy, 2022, §6.1.1). Cross-entropy loss gives this definition with the cross-entropy and KL divergence, and the identity $H(p, q) = H(p) + \mathrm{KL}(p \,\Vert\, q)$ that links them. For two variables (MacKay, 2003, eqs. 8.1 and 8.4),

H(X, Y) = -\sum_{x, y} p(x, y) \log p(x, y), \qquad H(X \mid Y) = -\sum_{x, y} p(x, y) \log p(x \mid y),

and the chain rule relates them (MacKay, 2003, eq. 8.7):

H(X, Y) = H(X) + H(Y \mid X) = H(Y) + H(X \mid Y).

Applied repeatedly, it gives $H(X_1, \dots, X_n) = \sum_i H(X_i \mid X_1, \dots, X_{i-1})$ (Murphy, 2022, eq. 6.18). The joint entropy is at most the sum of the marginal entropies, with equality exactly when the variables are independent, and at least the larger of the two (Murphy, 2022, §6.1.3). Conditioning reduces entropy on average, $H(X \mid Y) \le H(X)$ . One particular observation $y$ can still raise the uncertainty, $H(X \mid y) > H(X)$ ; only the average over $y$ cannot (Murphy, 2022, §6.1.4; MacKay, 2003, Exercise 8.2).

Mutual information

MacKay (2003, eq. 8.8) defines the mutual information as the reduction in uncertainty,

I(X; Y) = H(X) - H(X \mid Y),

and Murphy (2022, eqs. 6.47, 6.49 and 6.51) gives the equivalent forms

I(X; Y) = \mathrm{KL}\big(p(x, y) \,\Vert\, p(x)\, p(y)\big) = H(Y) - H(Y \mid X) = H(X) + H(Y) - H(X, Y).

The KL form shows that MI is non-negative and zero exactly when $X$ and $Y$ are independent. It is symmetric, and at most $\min(H(X), H(Y))$ for discrete variables (Murphy, 2022, eq. 6.73). The conditional mutual information $I(X; Y \mid Z) = H(X \mid Z) - H(X \mid Y, Z)$ is the information $Y$ adds about $X$ when $Z$ is already known (MacKay, 2003, eq. 8.10). It obeys its own chain rule, $I(Z, Y; X) = I(Z; X) + I(Y; X \mid Z)$ (Murphy, 2022, eq. 6.62).

The form $H(A) + H(B) - H(A, B)$ , with entropies estimated from a joint intensity histogram, is the one image registration maximizes. That article and the paper by Viola and Wells cover MI as an alignment cost.

Differential entropy

For a continuous variable with density $p$ , the differential entropy is $h(X) = -\int p(x) \log p(x)\, dx$ (Murphy, 2022, eq. 6.24). A uniform density on $(0, a)$ has $h = \log a$ . Densities can exceed 1, so $h$ can be negative: $U(0, 1/8)$ has $h = -3$ bits. A Gaussian with variance $\sigma^2$ has $h = \tfrac12 \log(2\pi e \sigma^2)$ (Murphy, 2022, eqs. 6.25–6.27).

Murphy (2022, §6.1.6), citing Cover and Thomas, connects $h$ with discrete entropy: an $n$ -bit quantization of $X$ has entropy of about $h(X) + n$ . Written for a bin width $\Delta$ , this is $H \approx h(X) - \log_2 \Delta$ (our restatement, with $\Delta = 2^{-n}$ ). So the entropy of a histogram depends on its bin width. It is not a property of the image alone until the bin width is fixed.

Two properties follow from the change-of-variables formula $p_{g(X)}(g(x)) = p_X(x) / \lvert g'(x) \rvert$ for an invertible, differentiable $g$ (our derivation). First, $h(g(X)) = h(X) + \mathbb{E}[\log \lvert g'(X) \rvert]$ , so $h(aX) = h(X) + \log \lvert a \rvert$ : contrast stretching by $a$ changes the differential entropy. Second, the Jacobian cancels in the ratio inside $I$ , so $I(g(X); Y) = I(X; Y)$ : mutual information does not change under any invertible map of either variable. MI between continuous variables is also non-negative, unlike $h$ (Murphy, 2022, eq. 6.48).

KL divergence

For distributions over the same states, $\mathrm{KL}(p \,\Vert\, q) = \sum_k p_k \log (p_k / q_k)$ , with the integral form for densities (Murphy, 2022, eqs. 6.29–6.30). It is non-negative and zero only when $p = q$ , which MacKay (2003, §2.6, eqs. 2.45–2.46) calls Gibbs’ inequality. In general $\mathrm{KL}(p \,\Vert\, q) \ne \mathrm{KL}(q \,\Vert\, p)$ , so although it is sometimes called the KL distance, it is not a distance (MacKay, 2003, §2.6). It is infinite wherever $p_k > 0$ and $q_k = 0$ . For two scalar Gaussians (Murphy, 2022, eq. 6.33),

\mathrm{KL}\big(\mathcal{N}(\mu_1, \sigma_1^2) \,\Vert\, \mathcal{N}(\mu_2, \sigma_2^2)\big) = \log \frac{\sigma_2}{\sigma_1} + \frac{\sigma_1^2 + (\mu_1 - \mu_2)^2}{2 \sigma_2^2} - \frac12 .

Minimizing $\mathrm{KL}(p_{\mathcal{D}} \,\Vert\, q)$ from the empirical distribution is maximum likelihood estimation (Murphy, 2022, §6.2.5), as that article and Cross-Entropy Loss explain.

Properties

Forward and reverse KL

To approximate a distribution $p$ by a simpler $q$ , one can minimize either $\mathrm{KL}(p \,\Vert\, q)$ or $\mathrm{KL}(q \,\Vert\, p)$ , and the two give different fits (Murphy, 2022, §6.2.6).

Forward KL, $\mathrm{KL}(p \,\Vert\, q)$ , also called the inclusive KL. Its minimization is an M-projection, or moment projection. The divergence is infinite wherever $p > 0$ and $q = 0$ , so the minimizer puts mass wherever $p$ has any. Murphy calls this zero-avoiding or mode-covering behaviour: a unimodal $q$ fitted to a bimodal $p$ spreads over both modes and typically over-estimates the support of $p$ .
Reverse KL, $\mathrm{KL}(q \,\Vert\, p)$ , also called the exclusive KL. Its minimization is an I-projection, or information projection. The divergence is infinite wherever $q > 0$ and $p = 0$ , so the minimizer avoids regions $p$ does not cover. Murphy calls this zero-forcing or mode-seeking: a unimodal $q$ locks onto one mode and typically under-estimates the support.

The name “moment projection” has a concrete reading for exponential families (our derivation). If $\log q = \eta^\top T(x) - A(\eta)$ , then $\mathrm{KL}(p \,\Vert\, q) = \text{const} - \eta^\top \mathbb{E}_p[T] + A(\eta)$ . Because $\nabla A(\eta) = \mathbb{E}_q[T]$ , the gradient vanishes where $\mathbb{E}_q[T] = \mathbb{E}_p[T]$ . For a Gaussian $q$ , $T = (x, x^2)$ , so the forward fit has $p$ ’s mean and variance. A second derivation of ours gives the value of the reverse fit at a mode. If $p = \sum_j w_j p_j$ has well-separated components and $q = p_j$ , then $p \approx w_j q$ where $q$ has its mass, and $\mathrm{KL}(q \,\Vert\, p) \approx -\log w_j$ . The larger mode therefore gives the smaller reverse divergence.

Which direction can be computed matters in practice. Torralba, Isola and Freeman (2024, §32.4.1) note that the forward direction $\mathrm{KL}(p_{\text{data}} \,\Vert\, p_\theta)$ is maximum likelihood, needing only samples from the data, while the reverse direction usually cannot be computed because the data density is not available. Variational methods work the other way round. MacKay (2003, §33.1) approximates a known but intractable $P$ by minimizing $\mathrm{KL}(Q \,\Vert\, P)$ over a simple $Q$ , the setting of the evidence lower bound, whose gap is a reverse KL. In GANs, the reverse-KL term in the generator’s gradient is the one Mode Collapse connects to dropped modes.

The data processing inequality

If $X \to Y \to Z$ forms a Markov chain, that is, $Z$ is independent of $X$ given $Y$ , then

I(X; Y) \ge I(X; Z)

(Murphy, 2022, Theorem 6.3.1). The proof expands $I(X; Y, Z)$ by the chain rule in two orders, $I(X; Z) + I(X; Y \mid Z) = I(X; Y) + I(X; Z \mid Y)$ . The Markov condition makes $I(X; Z \mid Y) = 0$ , and $I(X; Y \mid Z) \ge 0$ (Murphy, 2022, eqs. 6.77–6.79). MacKay (2003, Exercise 8.9) states the theorem for the state of the world $w$ , the data $d$ and the processed data $r$ , and remarks that it says as much about the definition of information as about processing. Equality holds when $Z$ is a sufficient statistic of $Y$ for $X$ (Murphy, 2022, §6.3.9). An invertible map is one case.

For images, the condition is the whole content of the theorem. Thresholding, quantizing or tone-mapping a pixel’s value is a function of $Y$ , so the chain holds and the pixel’s information about its label can only fall. A filter that reads neighbouring pixels is a function of the neighbourhood, not of the pixel’s value. The theorem then bounds the output by the information in the whole neighbourhood, which can exceed the pixel’s own. No vision pipeline adds information about the scene that is not in its input. In our reading, what processing can do is make information easier to use, which is a different property.

Fano’s inequality

Mutual information bounds classification error. For any estimator $\hat{Y} = f(X)$ with error probability $P_e = P(\hat{Y} \ne Y)$ , and with $H(E)$ the binary entropy of the error event,

H(Y \mid X) \le H(E) + P_e \log \lvert \mathcal{Y} \rvert ,

so $P_e \ge (H(Y \mid X) - 1)/\log \lvert \mathcal{Y} \rvert$ (Murphy, 2022, Theorem 6.3.2, eqs. 6.82–6.84). Murphy uses it to justify selecting features by their mutual information with the label: maximizing $I(X; Y)$ minimizes $H(Y \mid X)$ , which lowers the bound on the error of any classifier. For two classes, $\log_2 2 = 1$ and the weakened form is vacuous, but the first form still constrains the error.

Estimating information from images

Exact values need known distributions, which real images do not supply. For scalar variables a simple estimate bins the values into a histogram, and Murphy (2022, §6.3.5) warns that the number and placement of bins can affect the result considerably. The alternative he describes is the k-nearest-neighbour estimator of Kraskov, Stögbauer and Grassberger, which, as he notes, scikit-learn’s mutual_info_regression implements. In higher dimensions binning suffers from the curse of dimensionality (Murphy, 2022, §6.1.6.3). Normalized forms divide MI by an entropy so that it lies between 0 and 1. Murphy (2022, eq. 6.74) uses $\min(H(X), H(Y))$ , scikit-learn’s normalized_mutual_info_score defaults to the arithmetic mean (average_method="arithmetic"), and registration uses Studholme et al.’s $(H(A) + H(B))/H(A, B)$ (see Image Registration). Clustering covers the partition versions and their chance correction.

Conventions that change the numbers

Log base. scipy.stats.entropy uses the natural log unless base is given, and it normalizes both inputs to sum to one before computing (SciPy 1.16.3 source). scikit-learn’s mutual_info_score and mutual_info_classif return nats. mutual_info_classif uses n_neighbors=3, rescales continuous features and adds noise of relative size $10^{-10}$ , so it takes a random_state. scikit-learn’s decision trees use base 2, as Decision Tree records.
Zeros. SciPy computes the KL terms with rel_entr, which gives infinity where $p_k > 0$ and $q_k = 0$ . A histogram with an empty bin therefore has infinite KL divergence from any model that gives that bin mass. Similarity and Distance Measures describes the symmetric, finite alternatives used for histograms.
Bin width. Discrete entropy of a quantized continuous variable grows by one bit per halving of the bin width, as above. State the bin width whenever a density is histogrammed.

Examples

Forward and reverse KL on a grey-level histogram, and the data processing inequality on a two-region image. The script runs with NumPy 2.3.5, SciPy 1.16.3 and scikit-learn 1.7.2 on one thread in about one to two seconds. Two runs printed identical output.

Part (a) uses a 128 × 128 image with region A in columns 0–95 and region B in columns 96–127, an area ratio of 3 : 1. Grey levels are binned into 256 bins of width 1/256, with centres $x_k = (k + 0.5)/256$ . Each region’s grey levels follow a discretized Gaussian: the Gaussian density evaluated at the bin centres and normalized over the 256 bins, with mean 0.3 for A, 0.7 for B and standard deviation 0.05. The exact histogram of the image is then the mixture $p = 0.75\, c(0.3, 0.05) + 0.25\, c(0.7, 0.05)$ . A single discretized Gaussian $q$ from the same family is fitted by BFGS on $(\mu, \log \sigma)$ with analytic gradients:

once by forward KL, from $\mu = 0.5$ ;
nine times by reverse KL, from $\mu_0 = 0.1, 0.2, \dots, 0.9$ , all with $\sigma_0 = 0.1$ .

Further checks:

the gradients against central differences;
the forward fit against an independent root-finding solution of the moment equations;
scipy.stats.entropy against the hand sums;
the image itself, sampled pixel by pixel from its region’s distribution with seed 0, whose empirical histogram has empty bins.

Part (b) asks how much a pixel’s value tells about its region label $X$ . The label image is 130 × 130, with $X = 1$ (region B) in columns 97–129. The analysed pixels are the 128 × 128 whose 3 × 3 window lies inside the image, so $P(X = 1) = 0.25$ exactly. A pixel’s continuous value is $Y_c = 0.3 + 0.4X + 0.1\varepsilon$ with $\varepsilon$ standard normal. Three processed versions are derived from it:

$v = \mathrm{clip}(\mathrm{round}(255\, Y_c), 0, 255)$ , the stored 8-bit value;
$Z_3 = \lfloor v / 32 \rfloor$ , its top three bits;
$Z_b = \lfloor S / 288 \rfloor$ , where $S$ is the sum of $v$ over the pixel’s 3 × 3 window: a box mean followed by the top three bits, in integer arithmetic.

All three mutual informations are exact:

$P(v \mid x)$ from normal CDF differences;
the window sum’s distribution by convolving the nine per-pixel distributions for each count of B labels in the window;
$I(X; Y_c)$ by quadrature, computed two ways.

A Monte Carlo run passes 40 noisy images (seed 1) through the actual NumPy and SciPy pipeline and scores them with scikit-learn’s mutual_info_score. The script also checks Fano’s inequality with the Bayes error $P_e = \sum_z \min_x P(x, z)$ . It applies two invertible maps to $Y_c$ : a scaling by 1/4 and a tone curve $g(y) = 1/(1 + e^{-8(y - 0.5)})$ . Finally it compares the entropy of $\mathrm{round}(255\, Y_c)$ without clipping (bin width 1/255) with $h(Y_c) + \log_2 255$ .

import os
for var in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
    os.environ[var] = "1"                     # one BLAS thread, before NumPy loads
import numpy as np
from scipy import ndimage, optimize, special, stats
from sklearn.metrics import mutual_info_score

LN2 = np.log(2.0)

# ---------------------------------------------------------------- (a) forward and reverse KL
K = 256
x = (np.arange(K) + 0.5) / K                  # grey-level bin centres, bin width 1/256

def log_c(mu, sigma):
    """log of the discretized Gaussian: exp(-(x-mu)^2 / 2 sigma^2) normalized over the bins."""
    a = -0.5 * ((x - mu) / sigma) ** 2
    return a - special.logsumexp(a)

w = np.array([0.75, 0.25])                    # area fractions of regions A and B
logp = np.logaddexp(np.log(w[0]) + log_c(0.3, 0.05), np.log(w[1]) + log_c(0.7, 0.05))
p = np.exp(logp)                              # exact histogram of the image's grey levels

def kl_bits(a, b):
    """sum a log2(a/b), with 0 log 0 = 0 and inf where a > 0 and b = 0."""
    m = a > 0
    with np.errstate(divide="ignore"):
        return np.sum(a[m] * (np.log(a[m]) - np.log(b[m]))) / LN2

def moments(r):
    m = np.sum(r * x)
    return m, np.sqrt(np.sum(r * (x - m) ** 2))

def fwd(t):                                   # KL(p || q) in bits, t = (mu, log sigma)
    mu, s = t[0], np.exp(t[1]); lq = log_c(mu, s); q = np.exp(lq)
    val = np.sum(p * (logp - lq)) / LN2
    gm = (np.sum(q * x) - np.sum(p * x)) / s**2
    gs = (np.sum(q * (x - mu) ** 2) - np.sum(p * (x - mu) ** 2)) / s**2
    return val, np.array([gm, gs]) / LN2

def rev(t):                                   # KL(q || p) in bits
    mu, s = t[0], np.exp(t[1]); lq = log_c(mu, s); q = np.exp(lq)
    f = lq - logp
    sm, ss = (x - mu) / s**2, (x - mu) ** 2 / s**2    # d log q / d mu, d log q / d log sigma
    g = [np.sum(q * (u - np.sum(q * u)) * f) for u in (sm, ss)]
    return np.sum(q * f) / LN2, np.array(g) / LN2

def fit(fun, mu0, s0):
    r = optimize.minimize(fun, [mu0, np.log(s0)], jac=True, method="BFGS",
                          options={"gtol": 1e-8, "maxiter": 1000})
    return r.x[0], np.exp(r.x[1]), r.success, np.linalg.norm(fun(r.x)[1])

# C9: analytic gradients against central differences
t0, e = np.array([0.45, np.log(0.12)]), 1e-6
for name, fun in (("KL(p||q)", fwd), ("KL(q||p)", rev)):
    num = np.array([(fun(t0 + d)[0] - fun(t0 - d)[0]) / (2 * e) for d in e * np.eye(2)])
    print(f"C9 {name} gradient check: max rel error {np.max(np.abs(num - fun(t0)[1]) / np.abs(num)):.1e}")

mp, sp = moments(p)
print(f"p: mean {mp:.5f}, s.d. {sp:.5f}")
fits = [("forward", 0.5) + fit(fwd, 0.5, 0.1)]
fits += [("reverse", m0) + fit(rev, m0, 0.1) for m0 in np.arange(1, 10) / 10]
print("fit      start      mu   sigma  KL(p||q)  KL(q||p)  success  |grad|")
for kind, m0, mu, s, ok, g in fits:
    q = np.exp(log_c(mu, s))
    print(f"{kind:8s} {m0:5.1f} {mu:7.5f} {s:7.5f} {kl_bits(p, q):9.5f} {kl_bits(q, p):9.5f}"
          f"  {str(ok):7s}  {g:.1e}")

# C1: moment matching, and an independent root-finding solution of the moment equations
_, _, mu_f, s_f, *_ = fits[0]
qf = np.exp(log_c(mu_f, s_f))
mq, sq = moments(qf)
root = optimize.root(lambda t: np.array(moments(np.exp(log_c(t[0], np.exp(t[1]))))) - [mp, sp],
                     [mp, np.log(sp)], method="hybr", tol=1e-14)
print(f"C1 forward fit: |mean_q - mean_p| {abs(mq - mp):.1e}, |sd_q - sd_p| {abs(sq - sp):.1e}; "
      f"root-finding mu {root.x[0]:.6f}, sigma {np.exp(root.x[1]):.6f}, "
      f"max param diff {max(abs(root.x[0] - mu_f), abs(np.exp(root.x[1]) - s_f)):.1e}")
qn = np.exp(log_c(mp, sp))                    # naive closed form: mu, sigma = p's moments
print(f"C2 naive mu = mean_p, sigma = sd_p: KL(p||q) {kl_bits(p, qn):.6f}, "
      f"minus forward optimum {kl_bits(p, qn) - kl_bits(p, qf):.1e}")

# C3: scipy.stats.entropy against the hand sums
diffs = []
for _, _, m_, s_, *_ in fits:
    q = np.exp(log_c(m_, s_))
    diffs += [abs(stats.entropy(p, q, base=2) - kl_bits(p, q)),
              abs(stats.entropy(q, p, base=2) - kl_bits(q, p))]
print(f"C3 scipy.stats.entropy vs hand sums: max abs diff {max(diffs):.1e}")

# The image itself: region A in columns 0-95, region B in columns 96-127
rng = np.random.default_rng(0)
img = np.empty((128, 128), np.uint8)
img[:, :96] = rng.choice(K, size=(128, 96), p=np.exp(log_c(0.3, 0.05)))
img[:, 96:] = rng.choice(K, size=(128, 32), p=np.exp(log_c(0.7, 0.05)))
gA, gB = x[img[:, :96]], x[img[:, 96:]]
print(f"D1 region B pixels: {img[:, 96:].size} of {img.size}")
print(f"D2 region A mean {gA.mean():.4f} s.d. {gA.std():.4f}; region B mean {gB.mean():.4f} s.d. {gB.std():.4f}")
p_emp = np.bincount(img.ravel(), minlength=K) / img.size
print(f"C4 empirical histogram: {np.sum(p_emp == 0)} empty bins; KL(q_fwd||p_emp) hand {kl_bits(qf, p_emp)}, "
      f"scipy {stats.entropy(qf, p_emp, base=2)}; KL(p_emp||q_fwd) {kl_bits(p_emp, qf):.5f}")

# ---------------------------------------------------------------- (b) data processing inequality
lab = np.zeros((130, 130), np.int64)
lab[:, 97:] = 1                               # X = 1 in region B
valid = (slice(1, 129), slice(1, 129))        # pixels whose 3x3 window lies inside the image
X = lab[valid].ravel()
mu, sn = np.array([0.3, 0.7]), 0.1            # grey level of each region, noise s.d.
print(f"D3 analysed pixels with X = 1: {X.sum()} of {X.size}")
PX = np.bincount(X) / X.size

def H_bits(r):
    r = r[r > 0]
    return -np.sum(r * np.log2(r))

def mi_bits(joint):                           # joint table, rows x, columns z
    return H_bits(joint.sum(1)) + H_bits(joint.sum(0)) - H_bits(joint.ravel())

def pmf_v(m):                                 # P(v = k) for v = clip(round(255 y), 0, 255)
    cdf = stats.norm.cdf(((np.arange(256) + 0.5) / 255 - m) / sn)
    cdf[-1] = 1.0
    return np.diff(np.concatenate([[0.0], cdf]))

PV = np.array([pmf_v(m) for m in mu])         # P(v | x)
PZ3 = PV.reshape(2, 8, 32).sum(2)             # top 3 bits: v // 32
cnt = ndimage.correlate(lab, np.ones((3, 3), np.int64), mode="constant")[valid].ravel()
nxc = np.zeros((2, 10))                       # pixels with label x and c B-labels in the window
np.add.at(nxc, (X, cnt), 1.0)
PZb_c = np.zeros((10, 8))                     # P(Zb | c), Zb = (window sum of v) // 288
for c in range(10):
    s = np.array([1.0])
    for k in range(9):
        s = np.convolve(s, PV[1] if k < c else PV[0])
    PZb_c[c] = np.bincount(np.arange(s.size) // 288, weights=s, minlength=8)
J = {"v": PX[:, None] * PV, "Z3": PX[:, None] * PZ3, "Zb": nxc @ PZb_c / X.size}
I = {k: mi_bits(j) for k, j in J.items()}

# I(X; Y_c) for the continuous value, two ways on one grid
y, dy = np.linspace(-0.6, 1.6, 2_200_001, retstep=True)
fx = np.array([stats.norm.pdf(y, m, sn) for m in mu])
f = PX @ fx
def h_grid(dens, step):                       # differential entropy in bits by quadrature
    d = dens[dens > 0]
    return -np.sum(d * np.log2(d)) * step
hY, hYX = h_grid(f, dy), 0.5 * np.log2(2 * np.pi * np.e * sn**2)
post = PX[1] * fx[1] / f
Hb = lambda r: -(special.xlogy(r, r) + special.xlogy(1 - r, 1 - r)) / LN2
I_c = hY - hYX
print(f"C5 I(X;Yc): h(Y)-h(Y|X) {I_c:.9f}, H(X)-H(X|Y) {H_bits(PX) - np.sum(f * Hb(post)) * dy:.9f}")
print(f"H(X) = {H_bits(PX):.7f} bits")
print(f"I(X;Yc) {I_c:.7f} >= I(X;v) {I['v']:.7f} >= I(X;Z3) {I['Z3']:.7f};  I(X;Zb) {I['Zb']:.7f}")
print(f"  differences: Yc-v {I_c - I['v']:.2e}, v-Z3 {I['v'] - I['Z3']:.2e}, Zb-v {I['Zb'] - I['v']:.2e}")

# Monte Carlo through the actual pipeline, scored with scikit-learn (nats -> bits)
rng = np.random.default_rng(1)
mc = {"v": [], "Z3": [], "Zb": []}
for _ in range(40):
    Yc = mu[lab] + sn * rng.standard_normal(lab.shape)
    v = np.clip(np.round(255 * Yc), 0, 255).astype(np.int64)
    zb = ndimage.correlate(v, np.ones((3, 3), np.int64), mode="constant") // 288
    for k, z in (("v", v), ("Z3", v // 32), ("Zb", zb)):
        mc[k].append(z[valid].ravel())
for k in mc:
    est = mutual_info_score(np.tile(X, 40), np.concatenate(mc[k])) / LN2
    print(f"C6 I(X;{k}): exact {I[k]:.5f}, Monte Carlo {est:.5f}, diff {est - I[k]:+.1e}")

# Fano's inequality: H(X|Z) <= H_b(Pe) + Pe log2 2, Pe the Bayes error
for k, j in J.items():
    pe = j.min(0).sum()
    print(f"H8 {k:2s}: H(X|Z) {H_bits(PX) - I[k]:.5f} <= H_b(Pe) + Pe = {Hb(pe) + pe:.5f}  (Pe {pe:.5f})")

# Invertible maps of Yc: MI unchanged, differential entropy not
def mapped(ginv, dginv, lo, hi):              # densities of g(Yc) on a grid in the new variable
    u, du = np.linspace(lo, hi, 1_500_001, retstep=True)
    yy = ginv(u)
    with np.errstate(over="ignore"):
        gx = np.array([stats.norm.pdf(yy, m, sn) * np.abs(dginv(u)) for m in mu])
    gx = np.nan_to_num(gx)
    hg = h_grid(PX @ gx, du)
    return hg, hg - sum(PX[i] * h_grid(gx[i], du) for i in range(2))
h1, I1 = mapped(lambda u: 4 * u, lambda u: 4 + 0 * u, -0.15, 0.4)
e2 = 1e-15
h2, I2 = mapped(lambda u: 0.5 + np.log(u / (1 - u)) / 8, lambda u: 1 / (8 * u * (1 - u)), e2, 1 - e2)
sg = special.expit(8 * (y - 0.5))
Elog = np.sum(f * np.log2(8 * sg * (1 - sg))) * dy
print(f"H6 Yc/4: h {h1:.6f} (h(Yc) {hY:.6f}, change {h1 - hY:+.7f}), I diff {I1 - I_c:+.1e}")
print(f"H6 tone curve: h {h2:.6f} (change {h2 - hY:+.7f}, E[log2 g'] {Elog:+.7f}), I diff {I2 - I_c:+.1e}")

# Quantization: entropy of round(255 Yc) without clipping is about h(Yc) + log2 255
k = np.arange(-300, 556)
P = PX @ np.array([np.diff(stats.norm.cdf((np.r_[k - 0.5, k[-1] + 0.5] / 255 - m) / sn)) for m in mu])
print(f"H7 H(v unclipped) {H_bits(P):.6f}, h(Yc) + log2 255 {hY + np.log2(255):.6f}, "
      f"diff {H_bits(P) - hY - np.log2(255):+.1e}")

Output:

C9 KL(p||q) gradient check: max rel error 4.3e-11
C9 KL(q||p) gradient check: max rel error 3.8e-10
p: mean 0.40000, s.d. 0.18028
fit      start      mu   sigma  KL(p||q)  KL(q||p)  success  |grad|
forward    0.5 0.39074 0.19157   1.01461   2.65298  True     3.3e-10
reverse    0.1 0.30001 0.05002  10.71896   0.41496  True     1.2e-10
reverse    0.2 0.30001 0.05002  10.71896   0.41496  True     2.2e-10
reverse    0.3 0.30001 0.05002  10.71896   0.41496  True     5.4e-10
reverse    0.4 0.30001 0.05002  10.71896   0.41496  True     1.1e-10
reverse    0.5 0.30001 0.05002  10.71896   0.41496  True     1.1e-11
reverse    0.6 0.69997 0.05007  33.71518   1.99976  False    1.2e-07
reverse    0.7 0.69997 0.05007  33.71518   1.99976  True     7.8e-09
reverse    0.8 0.30001 0.05002  10.71896   0.41496  True     4.7e-10
reverse    0.9 0.30001 0.05002  10.71896   0.41496  True     2.7e-13
C1 forward fit: |mean_q - mean_p| 8.3e-12, |sd_q - sd_p| 2.3e-12; root-finding mu 0.390741, sigma 0.191575, max param diff 1.1e-11
C2 naive mu = mean_p, sigma = sd_p: KL(p||q) 1.019177, minus forward optimum 4.6e-03
C3 scipy.stats.entropy vs hand sums: max abs diff 1.8e-15
D1 region B pixels: 4096 of 16384
D2 region A mean 0.3001 s.d. 0.0502; region B mean 0.7008 s.d. 0.0501
C4 empirical histogram: 79 empty bins; KL(q_fwd||p_emp) hand inf, scipy inf; KL(p_emp||q_fwd) 1.02260
D3 analysed pixels with X = 1: 4096 of 16384
C5 I(X;Yc): h(Y)-h(Y|X) 0.737303341, H(X)-H(X|Y) 0.737303341
H(X) = 0.8112781 bits
I(X;Yc) 0.7373033 >= I(X;v) 0.7372819 >= I(X;Z3) 0.7160762;  I(X;Zb) 0.8085480
  differences: Yc-v 2.14e-05, v-Z3 2.12e-02, Zb-v 7.13e-02
C6 I(X;v): exact 0.73728, Monte Carlo 0.73719, diff -9.6e-05
C6 I(X;Z3): exact 0.71608, Monte Carlo 0.71543, diff -6.5e-04
C6 I(X;Zb): exact 0.80855, Monte Carlo 0.80865, diff +9.8e-05
H8 v : H(X|Z) 0.07400 <= H_b(Pe) + Pe = 0.15583  (Pe 0.01916)
H8 Z3: H(X|Z) 0.09520 <= H_b(Pe) + Pe = 0.17937  (Pe 0.02275)
H8 Zb: H(X|Z) 0.00273 <= H_b(Pe) + Pe = 0.00495  (Pe 0.00036)
H6 Yc/4: h -2.537529 (h(Yc) -0.537529, change -2.0000000), I diff -2.2e-16
H6 tone curve: h -0.507369 (change +0.0301603, E[log2 g'] +0.0301603), I diff +0.0e+00
H7 H(v unclipped) 7.456895, h(Yc) + log2 255 7.456824, diff +7.1e-05

The checks held. The data have the planned properties: region B holds 4,096 of the 16,384 pixels in both parts, and the sampled regions have the planned means and standard deviations to within 0.003. The analytic gradients agree with finite differences. The forward fit matches $p$ ’s moments, and its parameters agree with the root-finding solution. scipy.stats.entropy agrees with the hand sums, and the two formulas for $I(X; Y_c)$ agree. The Monte Carlo estimates differ from the exact values by at most $6.5 \times 10^{-4}$ bits. The naive parameter choice does not beat the optimizer, and the reverse divergence from the sampled image’s histogram is infinite both by hand and in SciPy.

Forward KL covers, reverse KL seeks. The forward fit sits between the modes with a standard deviation near four times a component’s, as hypothesized (held). Its own mean and standard deviation on the grid equal $p$ ’s, 0.40000 and 0.18028. Its parameters, $\mu = 0.39074$ and $\sigma = 0.19157$ , do not. In our reading, the reason is that the grid stops at 0 and 1, so a Gaussian with these parameters loses some tail, and the fitted parameters compensate for it. Setting the parameters to $p$ ‘s moments directly costs $4.6 \times 10^{-3}$ bits more. The reverse fits reached exactly two optima, one on each mode, with $\sigma$ within 0.2% of the components’ 0.05. At each optimum the divergence is close to $-\log_2 w$ : 0.41496 bits against $-\log_2 0.75 = 0.41504$ , and 1.99976 against 2. The larger mode gives the lower value; every part of this hypothesis held. Starts at 0.1–0.5, 0.8 and 0.9 reached the optimum at 0.3, and starts at 0.6 and 0.7 reached the one at 0.7.

Each fit is poor by the other measure, as hypothesized (held). The forward fit’s reverse divergence is 2.65 bits, more than either reverse optimum. The reverse fits’ forward divergences are 10.7 and 33.7 bits, against 1.01 for the forward fit, because each leaves a whole mode in its far tail. On the sampled image’s histogram, 79 of the 256 bins are empty. The reverse divergence from the forward fit is then infinite, both by hand and in SciPy, while the forward divergence stays near its exact value. Zero-forcing is not an abstraction for real histograms: any empty bin makes the reverse KL of a Gaussian infinite.

Processing a pixel loses information; pooling pixels does not contradict that. Along the chain $X \to Y_c \to v \to Z_3$ the information falls at every step (held). Storing 8 bits loses $2.1 \times 10^{-5}$ bits, and keeping three bits loses 0.021 bits, about 3% of what the pixel carries. $Z_3$ merges each run of 32 grey levels into one value. The blurred, quantized value $Z_b$ carries 0.809 bits, more than the stored pixel’s 0.737 and close to the label’s entropy of 0.811 (held). This is not a violation of the inequality (our derivation from the construction). $Z_b$ is computed from the nine values in the window, not from $v$ alone, so $X$ and $Z_b$ are not independent given $v$ , and Theorem 6.3.1’s condition fails for the chain $X \to v \to Z_b$ . The neighbours mostly share the pixel’s label, and their noise is independent, so their mean is a less noisy reading of the label. The inequality still holds for the chain through the whole window, whose information about $X$ is at least 0.809 bits.

Fano’s inequality held in all three cases. The bound is loose: for the stored pixel, $H(X \mid v) = 0.074$ bits against a bound of 0.156.

Differential entropy moves, mutual information does not. These two checks illustrate identities rather than test them: their residuals were known from a planning script before the run. Both held. The pixel’s differential entropy is already negative, −0.54 bits, which is possible because a density can exceed 1. Dividing by 4 lowers it by exactly 2 bits, and the tone curve raises it by $\mathbb{E}[\log_2 g'(Y_c)] = 0.030$ bits. The mutual information with the label is unchanged under both maps, to rounding error. Quantizing without clipping gives an entropy of 7.4569 bits, within $7.1 \times 10^{-5}$ of $h(Y_c) + \log_2 255$ , our restatement of Murphy’s §6.1.6 for this bin width.

Process note. The batch plan’s entry for part (b) reads “Region label X, noisy pixel Y and processed Z (Gaussian blur, then 3-bit quantization): exact I(X; Y) ≥ I(X; Z) (PML1 Thm 6.3.1)”. Before any script ran, we made five departures from it, logged them and had them approved at the pilot:

The blur is a 3 × 3 box mean on integers, whose distribution is an exact convolution.
Because the blur breaks the Markov chain when $Y$ is the pixel’s own value, the inequality is tested on $X \to Y_c \to v \to Z_3$ , and the blur is a separate case with the hypothesis $I(X; Z_b) > I(X; v)$ .
The noise standard deviation (0.1), the 130 × 130 image with its 128 × 128 analysed pixels, and the 8-bit rounding, which the batch plan leaves open, were fixed in our plan.
Checks were added: the root-finding reference, the naive moment-matched $q$ , the empirical histogram, the second MI formula, the Monte Carlo cross-check, the convergence and gradient checks, the tone curve, the quantization identity and Fano’s inequality.
The discretized Gaussian is the density at bin centres, normalized over the bins.

One bar was changed after a planning script printed a value that failed it. The convergence check’s first form was “every BFGS fit (one forward, nine reverse) reports success and ends with gradient norm < 1e-7”. The planning script printed “success False, |grad| 1.2e-07” for the reverse fit started at 0.6, and that check failed. The run above prints the same. The check was changed before the pilot to a gradient norm below $10^{-6}$ for every fit, with BFGS’s success flag (its line-search status) printed but not scored. The changed check held, with a largest norm of $1.2 \times 10^{-7}$ . Because it was set after the value was seen, it is an illustration, not a test. Also after that planning script, the quadrature grid for the mapped variables was changed from 2,200,001 to 1,500,001 points. With equal point counts, the sums for $Y_c/4$ were the sums for $Y_c$ scaled point by point, so they agreed by construction. The pilot was the full script, and the run shown is the same script with the same seeds; its output is identical to the pilot’s.

Common Misconceptions

“The KL divergence is a distance.” It is not symmetric (MacKay, 2003, §2.6). Murphy (2022, §6.2) introduces it as a divergence, which, unlike a metric, need not be symmetric or satisfy the triangle inequality, and Similarity and Distance Measures reports from Rubner et al. that it does not satisfy it. In the example, the forward fit is 1.01 bits from $p$ in one direction and 2.65 in the other. Which direction is minimized decides what the fit looks like.
“The data processing inequality means no filter can help.” It says that no function of the data adds information about the quantity of interest. A filter that pools neighbours uses more data than one pixel, and in the example it raises the pixel-level information from 0.737 to 0.809 bits. Processing can also make information usable, for example for a threshold, without adding any.
“Differential entropy is the information content of an image.” It can be negative, and it shifts with the units: the same pixel values have a differential entropy 2 bits lower when divided by 4. Murphy (2022, §6.1.6) shows that it can be negative. Mutual information between continuous variables has no such dependence (our derivation above, from the change of variables).
“Histogram entropy is a fixed property of the image.” It depends on the bin width, by about one bit per halving (Murphy, 2022, §6.1.6), and estimates of MI depend on the number and placement of bins (Murphy, 2022, §6.3.5).
“Three-variable information can be drawn as a Venn diagram.” MacKay (2003, Exercise 8.8) calls this representation misleading and points to the case $z = x + y \bmod 2$ with independent fair bits $x$ and $y$ . In our working of that case, $I(X; Y) = 0$ but $I(X; Y \mid Z) = 1$ bit: conditioning can increase mutual information, which overlapping circles cannot show.
“Conditioning on an observation always reduces uncertainty.” Only on average. A particular observation can leave the posterior more uncertain than the prior (MacKay, 2003, Exercise 8.2; Murphy, 2022, §6.1.4).

Where It Is Used

Training losses. The cross-entropy loss and label smoothing are cross-entropies. Knowledge distillation matches a student to a teacher by cross-entropy or KL. Maximum likelihood is forward KL from the data.
Compression. Image compression bounds the cost of a prefix code by the first-order entropy of its symbols.
Registration and stereo. Image registration, medical image registration and the Viola–Wells paper maximize MI between intensities. Semi-global matching can use MI as a matching cost.
Segmentation and features. Kapur’s method in thresholding chooses the threshold that maximizes the sum of the two classes’ entropies. Decision trees split by information gain. Fano’s inequality motivates choosing features by their MI with the label.
Comparing distributions. Similarity and distance measures uses KL and its symmetric relatives for histograms. Clustering scores partitions by (normalized, adjusted) MI. t-SNE minimizes a KL divergence between neighbourhood distributions and sets its scale by perplexity, an exponentiated entropy.
Generative models. Generative models are fitted by forward KL when their likelihood is tractable. The evidence lower bound of the variational autoencoder has a reverse-KL gap. GAN analyses connect the reverse KL to mode collapse. The Inception Score is the exponential of an expected KL divergence. Latent Space describes MI-based disentanglement scores.
Representation learning and uncertainty. The InfoNCE loss is a lower bound on MI between views, capped at $\log K$ , and underlies contrastive learning. Uncertainty estimation separates disagreement between models (an MI) from each model’s own entropy. Active perception chooses the next view by expected information gain.