02 - Mutual Information and Kullback-Leibler Divergence

If you know the verb of a sentence, how much do you know about the subject? And if you describe a coin with the wrong distribution, how many bits do you pay for each toss? These are two different questions with the same mathematical answer. The mutual information I(X;Y) = H(X) − H(X|Y) is the average reduction in uncertainty about X that comes from knowing Y.

The Kullback-Leibler divergence D(p||q) = Σ p(x) log₂(p(x)/q(x)) measures how many bits per symbol are wasted by using a code built for q when the data really come from p.

The link is I(X;Y) = D(p(x,y) || p(x)p(y)): mutual information is the distance from independence.

This article computes both quantities, proves that the divergence is never negative, and tests on numbers three claims that are usually accepted on trust: the asymmetry, the fact that it is not a distance, and the difference between minimising D(p||q) and D(q||p). Every listing was run twice with identical output (Python 3.14.4, numpy 2.3.5, sympy 1.14.0, fixed seeds, no network). It presupposes article 01 of the series, on entropy and the chain rule.

Mutual information: three ways to compute it

MacKay defines mutual information with I(X;Y) = H(X) − H(X|Y) (equation 8.8 of §8.1) and states two properties: it is symmetric and never negative.

With the chain rule you get the equivalent forms H(X) + H(Y) − H(X,Y) and H(Y) − H(Y|X): the symmetry can be read in the second.

Shannon reaches the same quantity from the problem of noise: in section 12 of his paper, "Equivocation and channel capacity", he observes that for a noisy channel the transmission rate is obtained by subtracting from the entropy of the source the conditional entropy of the message given the received signal, R = H(x) − H_y(x). This is mutual information, which in this passage Shannon does not call by that name.

The listing computes it on the joint distribution of article 01 (subject cat or dog, verb runs or sleeps, with probabilities 0.4, 0.1, 0.3 and 0.2, invented for the example) in three ways: as a difference of entropies, as the divergence between the joint distribution and the product of the marginals (the formula in the opening) and as H(X) − H(X|Y).

Then it checks the two extremes: if Y is a copy of X, the mutual information equals H(X), that is, 1 bit for a fair coin; if the variables are independent, the joint distribution coincides with the product and the value is zero.

from math import log2

joint = {
    ("cat", "runs"): 0.4, ("cat", "sleeps"): 0.1,
    ("dog", "runs"): 0.3, ("dog", "sleeps"): 0.2,
}
px, py = {}, {}
for (x, y), p in joint.items():
    px[x] = px.get(x, 0) + p
    py[y] = py.get(y, 0) + p


def H(ps):
    return sum(p * log2(1 / p) for p in ps if p > 0)


# 1) as a difference of entropies: I = H(X) + H(Y) - H(X,Y)
I_entropie = H(px.values()) + H(py.values()) - H(joint.values())

# 2) as the KL divergence between the joint distribution and the product of the marginals
I_kl = sum(p * log2(p / (px[x] * py[y])) for (x, y), p in joint.items())

# 3) as a reduction of uncertainty: H(X) - H(X|Y)
H_x_dato_y = sum(
    py[y] * H([joint[(x, y)] / py[y] for x in px]) for y in py
)
I_riduzione = H(px.values()) - H_x_dato_y

print(f"I(X;Y) as a difference of entropies : {I_entropie:.6f} bit")
print(f"I(X;Y) as D_KL(p(x,y) || p(x)p(y))  : {I_kl:.6f} bit")
print(f"I(X;Y) = H(X) - H(X|Y)              : {I_riduzione:.6f} bit")
print(f"in nat: {I_kl * 0.6931471805599453:.6f}")

# extreme case: Y is the coin itself (Y determines X), I = H(X) = 1 bit
joint2 = {("H", "H"): 0.5, ("T", "T"): 0.5}
print("Y = X:", sum(p * log2(p / (0.5 * 0.5)) for p in joint2.values()), "bit")
# opposite case: independent, the joint distribution is the product
ind = {(x, y): px[x] * py[y] for x in px for y in py}
print("independent:", round(sum(p * log2(p / (px[x] * py[y])) for (x, y), p in ind.items()), 12), "bit")

# Output:
# I(X;Y) as a difference of entropies : 0.034852 bit
# I(X;Y) as D_KL(p(x,y) || p(x)p(y))  : 0.034852 bit
# I(X;Y) = H(X) - H(X|Y)              : 0.034852 bit
# in nat: 0.024157
# Y = X: 1.0 bit
# independent: 0.0 bit

The three calculations agree: 0.034852 bits, that is, 0.024157 nats. It is a small number, and consistent with article 01: the verb removes on average less than 0.035 bits of uncertainty about the subject.

A second independent check, with exact fractions and sympy symbolic logarithms, gives 0.0348515545597 bits.

To see the opposite case you need a strong dependence: a coin X and a die Y that shows 6 if the coin is heads and a uniform value between 1 and 5 if it is tails. Here Y determines X (6 means heads, anything else means tails), so H(X|Y) = 0 and I = H(X) = 1 bit; but I remains smaller than H(Y), because the die has an uncertainty of its own, 2.1610 bits, that the coin does not explain.

from math import log2
import sympy as sp

# Exact check with fractions: the joint distribution against the product of the marginals
R = sp.Rational
joint = {("cat", "runs"): R(2, 5), ("cat", "sleeps"): R(1, 10),
         ("dog", "runs"): R(3, 10), ("dog", "sleeps"): R(1, 5)}
px = {"cat": R(1, 2), "dog": R(1, 2)}
py = {"runs": R(7, 10), "sleeps": R(3, 10)}
I_esatta = sum(p * sp.log(p / (px[x] * py[y]), 2) for (x, y), p in joint.items())
print("exact I (sympy):", sp.N(I_esatta, 12), "bit")

# Coin X and die Y: if heads the die is 6, if tails it is uniform on 1..5
joint_md = {}
for y in range(1, 7):
    joint_md[("H", y)] = 0.5 if y == 6 else 0.0
    joint_md[("T", y)] = 0.0 if y == 6 else 0.1
px = {"H": 0.5, "T": 0.5}
py = {y: joint_md[("H", y)] + joint_md[("T", y)] for y in range(1, 7)}
I = sum(p * log2(p / (px[x] * py[y])) for (x, y), p in joint_md.items() if p > 0)
H = lambda ps: sum(p * log2(1 / p) for p in ps if p > 0)
print(f"H(X)={H(px.values()):.4f}  H(Y)={H(py.values()):.4f}  I(X;Y)={I:.4f} bit")
print(f"I = H(X): {abs(I - H(px.values())) < 1e-12}   I < H(Y): {I < H(py.values())}")

# Output:
# exact I (sympy): 0.0348515545597 bit
# H(X)=1.0000  H(Y)=2.1610  I(X;Y)=1.0000 bit
# I = H(X): True   I < H(Y): True

In general 0 ≤ I(X;Y) ≤ min(H(X), H(Y)): the upper bound follows from the fact that conditional entropies are not negative, and it is reached when one variable is a function of the other.

The Kullback-Leibler divergence and why it is never negative

Kullback and Leibler introduced the quantity in 1951, in the paper On Information and Sufficiency (Annals of Mathematical Statistics 22(1), pp. 79-86: metadata verified, text not read).

MacKay, §2.6, calls it relative entropy and gives the definition D(P||Q) = Σ P(x) log(P(x)/Q(x)) for two distributions on the same alphabet.

He adds two facts: Gibbs' inequality holds, D(P||Q) ≥ 0 with equality only if P = Q; and in general D(P||Q) ≠ D(Q||P), so, although it is sometimes called "KL distance", it is not properly a distance.

The proof uses Jensen's inequality: for a concave function such as the logarithm, the mean of the logarithms does not exceed the logarithm of the mean.

You sum only over the outcomes with p(x) > 0, and write D = −Σ p log(q/p). By Jensen, Σ p log(q/p) ≤ log Σ p · (q/p) = log Σ q, and the sum of q over the support of p is at most 1, so the logarithm is at most zero and D ≥ 0.

Equality requires that q/p be constant on the support and that the sum equal 1, that is, q = p. MacKay gives the same proof in the solution of exercise 2.26.

The case where q(x) = 0 where p(x) > 0 must be treated separately: the logarithm is infinite and the divergence is +∞.

Asymmetry, meaning in bits and cross-entropy

The listing computes the divergence between a fair coin p = (0.5, 0.5) and an unbalanced coin q = (0.9, 0.1) in both directions, in bits and in nats.

Then it gives the operational meaning: if for q you use an ideal code that assigns to each outcome log₂(1/q(x)) bits, but the data come from p, the average length is the cross-entropy H(p,q) = Σ p log₂(1/q) = H(p) + D(p||q).

MacKay states it in chapter 5: the divergence measures how many bits per symbol are wasted by using a code whose implicit probabilities are q while the true distribution is p. "Ideal" means non-integer lengths, which no real code has; article 03 shows what happens with integer lengths.

The simulation draws 200,000 outcomes from p with seed 7 and averages their surprise with respect to q: it is a second independent calculation of the formula, with the standard error stated. The last part shows the case that matters to anyone who trains classifiers: with p concentrated on a single class (one-hot) the entropy H(p) is zero, and the cross-entropy reduces to −ln q(correct class), identical to the divergence; the logits are three arbitrary numbers transformed by a softmax.

import numpy as np
from math import log, log2

rng = np.random.default_rng(7)


def kl(p, q):
    """D_KL(p||q) in bits; infinite if q is zero where p is not."""
    s = 0.0
    for a, b in zip(p, q):
        if a > 0:
            if b == 0:
                return float("inf")
            s += a * log2(a / b)
    return s


p, q = [0.5, 0.5], [0.9, 0.1]
print(f"D(p||q) = {kl(p, q):.4f} bit = {kl(p, q) * log(2):.4f} nat")
print(f"D(q||p) = {kl(q, p):.4f} bit = {kl(q, p) * log(2):.4f} nat")
print("different supports: D(p||r) =", kl(p, [1.0, 0.0]), "  D(r||p) =", kl([1.0, 0.0], p), "bit")

# Operational meaning: mean length of an ideal code built for q
H_p = sum(a * log2(1 / a) for a in p)
cross = sum(a * log2(1 / b) for a, b in zip(p, q))
print(f"H(p)={H_p:.4f}  H(p,q)={cross:.4f}  H(p)+D(p||q)={H_p + kl(p, q):.4f} bit")
x = rng.choice(2, size=200_000, p=p)
media = float(np.mean(-np.log2(np.array(q)[x])))
err = float(np.std(-np.log2(np.array(q)[x]), ddof=1) / np.sqrt(len(x)))
print(f"simulation, 200000 draws from p: {media:.4f} +/- {err:.4f} bits per symbol")

# With one-hot p the cross-entropy is -log q(correct class), equal to the KL
logits = np.array([2.0, 0.5, -1.0])
qs = np.exp(logits) / np.exp(logits).sum()
one_hot = np.array([0.0, 1.0, 0.0])
ce = -float(np.sum(one_hot[one_hot > 0] * np.log(qs[one_hot > 0])))
print(f"cross-entropy {ce:.6f} nat   -ln q(class) {-log(qs[1]):.6f}   KL {kl(one_hot, qs) * log(2):.6f}")

# Output:
# D(p||q) = 0.7370 bit = 0.5108 nat
# D(q||p) = 0.5310 bit = 0.3681 nat
# different supports: D(p||r) = inf   D(r||p) = 1.0 bit
# H(p)=1.0000  H(p,q)=1.7370  H(p)+D(p||q)=1.7370 bit
# simulation, 200000 draws from p: 1.7412 +/- 0.0035 bits per symbol
# cross-entropy 1.741311 nat   -ln q(class) 1.741311   KL 1.741311

The divergence depends on the direction: 0.7370 bits against 0.5310, that is, 0.5108 nats against 0.3681.

Mind the unit: 0.737 is bits, 0.511 is the same data in nats.

If q assigns probability zero to an outcome that p can produce, D(p||q) is infinite, while the opposite direction is 1 bit: the code of q cannot represent that outcome. The simulation gives 1.7412 ± 0.0035 bits against the 1.7370 of the formula, compatible. For a one-hot class, finally, cross-entropy and divergence coincide (1.741311 nats).

Since H(p) does not depend on q, minimising the cross-entropy with respect to q is the same thing as minimising D(p||q), and this is why that loss function has that name.

It is not a distance: the triangle breaks

A distance must be symmetric and satisfy the triangle inequality, d(a,c) ≤ d(a,b) + d(b,c). The asymmetry is already shown above; the triangle is tested with an exhaustive search.

The listing runs through all 729 triples of coins (probability of heads from 0.1 to 0.9 in steps of 0.1) and counts those in which going from P to R costs more than going via Q.

Then it looks at the Jensen-Shannon divergence, JS(p,q) = ½[D(p||m) + D(q||m)] with m = (p+q)/2, which is symmetric and bounded, and checks on 200,000 random triples with five outcomes (Dirichlet with parameter 0.3, seed 11) whether its square root respects the triangle.

import numpy as np
from math import log2, sqrt


def kl(p, q):
    return sum(a * log2(a / b) for a, b in zip(p, q) if a > 0)


def bern(a):
    return [a, 1 - a]


# Search a coarse grid for a triple of Bernoulli distributions that violates the triangle inequality
griglia = [i / 10 for i in range(1, 10)]
violazioni = []
for a in griglia:
    for b in griglia:
        for c in griglia:
            diretta = kl(bern(a), bern(c))
            via_b = kl(bern(a), bern(b)) + kl(bern(b), bern(c))
            if diretta > via_b + 1e-12:
                violazioni.append((diretta - via_b, a, b, c))
violazioni.sort(reverse=True)
print(f"violating triples out of {len(griglia) ** 3}: {len(violazioni)}")
_, a, b, c = violazioni[0]
print(f"example: P(1)={a}, Q(1)={b}, R(1)={c}")
print(f"D(P||R) = {kl(bern(a), bern(c)):.4f} bit  >  D(P||Q) + D(Q||R) = "
      f"{kl(bern(a), bern(b)):.4f} + {kl(bern(b), bern(c)):.4f} = "
      f"{kl(bern(a), bern(b)) + kl(bern(b), bern(c)):.4f} bit")

# Jensen-Shannon: symmetric, bounded by 1 bit; its square root, on 200000 triples with 5 outcomes, does not violate the triangle
def js(p, q):
    m = [(x + y) / 2 for x, y in zip(p, q)]
    return (kl(p, m) + kl(q, m)) / 2


rng = np.random.default_rng(11)
peggiore, massimo_js = np.inf, 0.0
for _ in range(200_000):
    p, q, r = (list(rng.dirichlet(np.ones(5) * 0.3)) for _ in range(3))
    peggiore = min(peggiore, sqrt(js(p, q)) + sqrt(js(q, r)) - sqrt(js(p, r)))
    massimo_js = max(massimo_js, js(p, q))
print(f"largest JS observed: {massimo_js:.4f} bit;  smallest triangle margin on sqrt(JS): {peggiore:.2e}")

# Output:
# violating triples out of 729: 168
# example: P(1)=0.1, Q(1)=0.6, R(1)=0.9
# D(P||R) = 2.5359 bit  >  D(P||Q) + D(Q||R) = 0.7944 + 0.4490 = 1.2435 bit
# largest JS observed: 0.9993 bit;  smallest triangle margin on sqrt(JS): 1.32e-02

In 168 triples out of 729 the triangle is violated. The printed example is clear-cut: going from 0.1 to 0.9 directly costs 2.5359 bits, going via 0.6 only 1.2435.

The Jensen-Shannon divergence stays below 1 bit (largest observed 0.9993) and is symmetric by construction. On the 200,000 triples the square root never violates the triangle, with a smallest margin of 0.013.

This is a numerical check, not a proof: that the root of JS is a metric is a result known in the literature, which I cite here but have not read, and for rigorous use it must be looked up at the source.

Which direction to minimise

The asymmetry changes the behaviour of whoever approximates. The listing takes a bimodal distribution p (two Gaussians at −3 and +3, discretised on a grid with step 0.1) and looks, with an exhaustive search over a grid of means and standard deviations, for the Gaussian q that minimises D(p||q) and the one that minimises D(q||p).

As a second independent calculation it compares with the moments of p: for the Gaussian family the minimiser of D(p||q) has the same mean and the same variance as p, and here the variance is 1 + 3² = 10.

import numpy as np

x = np.arange(-10, 10.0001, 0.1)
gauss = lambda m, s: np.exp(-0.5 * ((x - m) / s) ** 2)
p = 0.5 * gauss(-3, 1) + 0.5 * gauss(3, 1)
p /= p.sum()                                   # discretised bimodal distribution


def kl(a, b):
    ok = a > 0
    return float(np.sum(a[ok] * np.log2(a[ok] / np.maximum(b[ok], 1e-300))))


meglio = {"p||q": (np.inf,), "q||p": (np.inf,)}
for m in np.arange(-6, 6.01, 0.1):
    for s in np.arange(0.3, 6.01, 0.1):
        q = gauss(m, s)
        q /= q.sum()
        for chiave, valore in (("p||q", kl(p, q)), ("q||p", kl(q, p))):
            if valore < meglio[chiave][0]:
                meglio[chiave] = (valore, round(float(m), 1) + 0.0, round(float(s), 1))
for chiave, (v, m, s) in meglio.items():
    print(f"minimising D({chiave}): mu={m:+.1f}  sigma={s:.1f}  D={v:.4f} bit")
media = float(np.sum(x * p))
sigma = float(np.sqrt(np.sum((x - media) ** 2 * p)))
print(f"moments of p: mean {media:+.2f}, sigma {sigma:.3f}  (sqrt(1 + 3^2) = {np.sqrt(10):.3f})")

# Output:
# minimising D(p||q): mu=+0.0  sigma=3.2  D=0.6643 bit
# minimising D(q||p): mu=+3.0  sigma=1.0  D=0.9944 bit
# moments of p: mean +0.00, sigma 3.162  (sqrt(1 + 3^2) = 3.162)

The direction D(p||q) produces a wide Gaussian, with mean 0 and σ = 3.2, that covers both peaks and has its maximum where p is almost zero; the moments of p (σ = 3.162) confirm the result.

The direction D(q||p) picks a single peak, mean +3 and σ = 1 (by symmetry the minimum at −3 is worth the same as the one at +3; the listing reports one of them).

The explanation is in the formula: in D(p||q) the weight is p, so q is punished where p has mass and q has none; in D(q||p) the weight is q, and q is punished where it has mass and p does not, so it pays to concentrate on a single peak.

The result holds for this toy case (two symmetric peaks, a single family): it is not a rule about which direction is "better", and the choice depends on what an error costs in the application.

Estimating mutual information from counts

As for entropy, replacing probabilities with observed frequencies introduces a bias, and here it takes an even more insidious form: independent variables, for which I is zero, always come out as dependent. The listing draws 1000 independent pairs of two variables with 10 values each, computes the mutual information from the counts, and repeats 2000 times (seed 99).

import numpy as np

rng = np.random.default_rng(99)


def mi_bit(tabella):
    p = tabella / tabella.sum()
    px, py = p.sum(axis=1, keepdims=True), p.sum(axis=0, keepdims=True)
    ok = p > 0
    return float(np.sum(p[ok] * np.log2(p[ok] / (px @ py)[ok])))


k, N, prove = 10, 1000, 2000
stime = []
for _ in range(prove):
    t = np.zeros((k, k))
    np.add.at(t, (rng.integers(0, k, N), rng.integers(0, k, N)), 1)   # X and Y independent
    stime.append(mi_bit(t))
s = np.array(stime)
print(f"INDEPENDENT variables, {k}x{k} values, N={N}: true I = 0 bit")
print(f"I estimated from counts: {s.mean():.4f} +/- {s.std(ddof=1) / np.sqrt(prove):.4f} bit")
print(f"order of magnitude (k-1)^2 / (2 N ln 2) = {(k - 1) ** 2 / (2 * N * np.log(2)):.4f} bit")

# Output:
# INDEPENDENT variables, 10x10 values, N=1000: true I = 0 bit
# I estimated from counts: 0.0596 +/- 0.0002 bit
# order of magnitude (k-1)^2 / (2 N ln 2) = 0.0584 bit

The mean estimate is 0.0596 ± 0.0002 bits, not zero. Of the same order is (k−1)²/(2N ln 2) = 0.0584 bits, a first-order approximation that I use as an order of magnitude without having proved it. The practical meaning is that a "small" value of mutual information estimated on a sample is not evidence of dependence: it must be compared with what you get on shuffled data, and, according to that approximation, it grows with the number of values and decreases with the sample size.

Typical mistakes

• Swapping the two arguments of the divergence. D(p||q) is 0.7370 bits and D(q||p) 0.5310 bits for the same two coins.
• Calling it a distance. It is not symmetric and does not respect the triangle (168 triples out of 729).
• Changing unit without saying so. With the natural logarithm the values are in nats, with base 2 in bits, and dividing by ln 2 goes from the former to the latter.
• Forgetting the supports. If q has an outcome with probability zero where p has one, D(p||q) is infinite and an implementation that sums without checking gives inf or an error.
• Reading a small estimate as "independent". With 1000 samples and 10×10 values, independent variables give 0.06 bits.

The common thread is that the number you print depends on a choice that the formulas do not show: the direction, the unit, the support, the sample size. For every divergence or mutual information reported in a result it is worth writing these four things next to it.

What you can do now

You can compute mutual information in three equivalent ways and you know that it is a divergence from independence, you can prove with Jensen that the KL divergence is not negative and what is needed for it to be zero, and you can give it a meaning in wasted bits.

You can also show with numbers that it is not symmetric, that the triangle breaks and that the direction chosen changes the approximation, and you know that estimating these quantities from counts with few data gives values that are too high for mutual information and too low for entropy.

The next step is to turn entropy into code lengths: article 03 proves that no per-symbol code goes below H(X) and builds the optimal one; mutual information will return in article 04, where it measures how much gets through a noisy channel. The list of the articles is on the series page.

Later in the series you will also find article 05 on error-correcting codes and article 06 on lossy compression.