02 - Informazione mutua e divergenza di Kullback-Leibler

Se conosci il verbo di una frase, quanto sai del soggetto? E se descrivi una moneta con la distribuzione sbagliata, quanto paghi in bit per ogni lancio? Sono due domande diverse con la stessa risposta matematica. L'informazione mutua I(X;Y) = H(X) − H(X|Y) è la riduzione media di incertezza su X che deriva dal conoscere Y.

La divergenza di Kullback-Leibler D(p||q) = Σ p(x) log₂(p(x)/q(x)) misura quanti bit per simbolo si sprecano usando un codice costruito per q quando i dati vengono davvero da p.

Il legame è I(X;Y) = D(p(x,y) || p(x)p(y)): l'informazione mutua è la distanza dall'indipendenza.

Questo articolo calcola entrambe le quantità, dimostra che la divergenza non è mai negativa, e mette alla prova su numeri tre affermazioni che di solito si accettano sulla fiducia: l'asimmetria, il fatto che non sia una distanza, e la differenza fra minimizzare D(p||q) e D(q||p). Ogni listato è stato eseguito due volte con uscita identica (Python 3.14.4, numpy 2.3.5, sympy 1.14.0, semi fissi, nessuna rete). Presuppone l'articolo 01 della serie, sull'entropia e la regola della catena.

Informazione mutua: tre modi di calcolarla

MacKay definisce l'informazione mutua con I(X;Y) = H(X) − H(X|Y) (equazione 8.8 del §8.1) e ne enuncia due proprietà: è simmetrica e non è mai negativa.

Con la regola della catena si ottengono le forme equivalenti H(X) + H(Y) − H(X,Y) e H(Y) − H(Y|X): la simmetria si legge nella seconda.

Shannon arriva alla stessa quantità dal problema del rumore: nella sezione 12 del suo articolo, «Equivocation and channel capacity», osserva che per un canale rumoroso la velocità di trasmissione si ottiene sottraendo all'entropia della sorgente l'entropia condizionata del messaggio dato il segnale ricevuto, R = H(x) − H_y(x). È l'informazione mutua, che in questo passo Shannon non chiama così.

Il listato la calcola sulla congiunta dell'articolo 01 (soggetto gatto o cane, verbo corre o dorme, con probabilità 0,4, 0,1, 0,3 e 0,2, inventate per l'esempio) in tre modi: come differenza di entropie, come divergenza fra la congiunta e il prodotto delle marginali (la formula nell'apertura) e come H(X) − H(X|Y).

Poi controlla i due estremi: se Y è una copia di X l'informazione mutua vale H(X), cioè 1 bit per una moneta equa; se le variabili sono indipendenti, la congiunta coincide con il prodotto e il valore è zero.

from math import log2

joint = {
    ("gatto", "corre"): 0.4, ("gatto", "dorme"): 0.1,
    ("cane", "corre"): 0.3, ("cane", "dorme"): 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) per differenza di entropie: I = H(X) + H(Y) - H(X,Y)
I_entropie = H(px.values()) + H(py.values()) - H(joint.values())

# 2) come divergenza KL fra la congiunta e il prodotto delle marginali
I_kl = sum(p * log2(p / (px[x] * py[y])) for (x, y), p in joint.items())

# 3) come riduzione d'incertezza: 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) per differenza di entropie : {I_entropie:.6f} bit")
print(f"I(X;Y) come 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}")

# caso estremo: Y è la moneta stessa (Y determina X), I = H(X) = 1 bit
joint2 = {("T", "T"): 0.5, ("C", "C"): 0.5}
print("Y = X:", sum(p * log2(p / (0.5 * 0.5)) for p in joint2.values()), "bit")
# caso opposto: indipendenti, la congiunta è il prodotto
ind = {(x, y): px[x] * py[y] for x in px for y in py}
print("indipendenti:", round(sum(p * log2(p / (px[x] * py[y])) for (x, y), p in ind.items()), 12), "bit")

# Output:
# I(X;Y) per differenza di entropie : 0.034852 bit
# I(X;Y) come 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
# indipendenti: 0.0 bit

I tre calcoli coincidono: 0,034852 bit, cioè 0,024157 nat. È un numero piccolo, e coerente con l'articolo 01: il verbo toglie in media meno di 0,035 bit di incertezza sul soggetto.

Un secondo controllo indipendente, con frazioni esatte e logaritmi simbolici di sympy, dà 0,0348515545597 bit.

Per vedere il caso opposto serve una dipendenza forte: una moneta X e un dado Y che esce 6 se la moneta è testa e un valore uniforme fra 1 e 5 se è croce. Qui Y determina X (6 vuol dire testa, altro vuol dire croce), quindi H(X|Y) = 0 e I = H(X) = 1 bit; ma I resta minore di H(Y), perché il dado ha un'incertezza propria, 2,1610 bit, che la moneta non spiega.

from math import log2
import sympy as sp

# Verifica esatta con frazioni: la congiunta e il prodotto delle marginali
R = sp.Rational
joint = {("gatto", "corre"): R(2, 5), ("gatto", "dorme"): R(1, 10),
         ("cane", "corre"): R(3, 10), ("cane", "dorme"): R(1, 5)}
px = {"gatto": R(1, 2), "cane": R(1, 2)}
py = {"corre": R(7, 10), "dorme": R(3, 10)}
I_esatta = sum(p * sp.log(p / (px[x] * py[y]), 2) for (x, y), p in joint.items())
print("I esatta (sympy):", sp.N(I_esatta, 12), "bit")

# Moneta X e dado Y: se esce testa il dado è 6, se esce croce è uniforme su 1..5
joint_md = {}
for y in range(1, 7):
    joint_md[("T", y)] = 0.5 if y == 6 else 0.0
    joint_md[("C", y)] = 0.0 if y == 6 else 0.1
px = {"T": 0.5, "C": 0.5}
py = {y: joint_md[("T", y)] + joint_md[("C", 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:
# I esatta (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 generale 0 ≤ I(X;Y) ≤ min(H(X), H(Y)): il limite superiore discende dal fatto che le entropie condizionate non sono negative, ed è raggiunto quando una variabile è funzione dell'altra.

La divergenza di Kullback-Leibler e perché non è mai negativa

Kullback e Leibler introducono la quantità nel 1951, nell'articolo On Information and Sufficiency (Annals of Mathematical Statistics 22(1), pp. 79-86: metadati verificati, testo non letto).

MacKay, §2.6, la chiama entropia relativa e ne dà la definizione D(P||Q) = Σ P(x) log(P(x)/Q(x)) per due distribuzioni sullo stesso alfabeto.

Aggiunge due fatti: vale la disuguaglianza di Gibbs, D(P||Q) ≥ 0 con uguaglianza solo se P = Q; e in generale D(P||Q) ≠ D(Q||P), quindi, pur chiamata talvolta «distanza KL», non è propriamente una distanza.

La dimostrazione usa la disuguaglianza di Jensen: per una funzione concava come il logaritmo, la media dei logaritmi non supera il logaritmo della media.

Si somma solo sugli esiti con p(x) > 0, e si scrive D = −Σ p log(q/p). Per Jensen, Σ p log(q/p) ≤ log Σ p · (q/p) = log Σ q, e la somma di q sul supporto di p è al più 1, quindi il logaritmo è al più zero e D ≥ 0.

L'uguaglianza richiede che q/p sia costante sul supporto e che la somma valga 1, cioè q = p. MacKay propone la stessa dimostrazione nella soluzione dell'esercizio 2.26.

Il caso in cui q(x) = 0 dove p(x) > 0 va trattato a parte: il logaritmo è infinito e la divergenza vale +∞.

Asimmetria, significato in bit e cross-entropy

Il listato calcola la divergenza fra una moneta equa p = (0,5; 0,5) e una moneta sbilanciata q = (0,9; 0,1) nei due versi, in bit e in nat.

Poi dà il significato operativo: se si usa per q un codice ideale che assegna a ogni esito log₂(1/q(x)) bit, ma i dati vengono da p, la lunghezza media è la cross-entropy H(p,q) = Σ p log₂(1/q) = H(p) + D(p||q).

MacKay lo enuncia nel capitolo 5: la divergenza misura quanti bit per simbolo si sprecano usando un codice le cui probabilità implicite sono q mentre la distribuzione vera è p. «Ideale» significa lunghezze non intere, che nessun codice reale ha; l'articolo 03 mostra che cosa succede con lunghezze intere.

La simulazione estrae 200.000 esiti da p con seme 7 e ne media la sorpresa rispetto a q: è un secondo calcolo indipendente della formula, con l'errore standard dichiarato. L'ultima parte mostra il caso che interessa a chi addestra classificatori: con p concentrata su una sola classe (one-hot) l'entropia H(p) vale zero, e la cross-entropy si riduce a −ln q(classe giusta), identica alla divergenza; i logit sono tre numeri qualsiasi trasformati da una 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 bit; infinito se q è zero dove p non lo è."""
    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("supporti diversi: D(p||r) =", kl(p, [1.0, 0.0]), "  D(r||p) =", kl([1.0, 0.0], p), "bit")

# Significato operativo: lunghezza media di un codice ideale costruito per 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"simulazione, 200000 estrazioni da p: {media:.4f} +/- {err:.4f} bit per simbolo")

# Con p one-hot la cross-entropy è -log q(classe giusta), uguale alla 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(classe) {-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
# supporti diversi: 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
# simulazione, 200000 estrazioni da p: 1.7412 +/- 0.0035 bit per simbolo
# cross-entropy 1.741311 nat   -ln q(classe) 1.741311   KL 1.741311

La divergenza dipende dal verso: 0,7370 bit contro 0,5310, cioè 0,5108 nat contro 0,3681.

Attenzione all'unità: 0,737 sono bit, 0,511 sono gli stessi dati in nat.

Se q assegna probabilità zero a un esito che p può produrre, D(p||q) è infinita, mentre il verso opposto vale 1 bit: il codice di q non può rappresentare quell'esito. La simulazione dà 1,7412 ± 0,0035 bit contro i 1,7370 della formula, compatibili. Per una classe one-hot, infine, cross-entropy e divergenza coincidono (1,741311 nat).

Siccome H(p) non dipende da q, minimizzare la cross-entropy rispetto a q è la stessa cosa che minimizzare D(p||q), ed è per questo che quella funzione di perdita ha quel nome.

Non è una distanza: il triangolo si rompe

Una distanza deve essere simmetrica e soddisfare la disuguaglianza triangolare, d(a,c) ≤ d(a,b) + d(b,c). L'asimmetria è già dimostrata sopra; il triangolo si prova con una ricerca esaustiva.

Il listato percorre tutte le 729 terne di monete (probabilità di testa da 0,1 a 0,9 a passi di 0,1) e conta quelle in cui il passaggio da P a R costa più del passaggio via Q.

Poi guarda la divergenza di Jensen-Shannon, JS(p,q) = ½[D(p||m) + D(q||m)] con m = (p+q)/2, che è simmetrica e limitata, e controlla su 200.000 terne casuali a cinque esiti (Dirichlet con parametro 0,3, seme 11) se la sua radice quadrata rispetta il triangolo.

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]


# Cerca su una griglia grossolana una terna di Bernoulli che viola la disuguaglianza triangolare
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"terne violanti su {len(griglia) ** 3}: {len(violazioni)}")
_, a, b, c = violazioni[0]
print(f"esempio: 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: simmetrica, limitata da 1 bit; la sua radice, su 200000 terne a 5 esiti, non viola il triangolo
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"JS massima osservata: {massimo_js:.4f} bit;  margine minimo del triangolo su sqrt(JS): {peggiore:.2e}")

# Output:
# terne violanti su 729: 168
# esempio: 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
# JS massima osservata: 0.9993 bit;  margine minimo del triangolo su sqrt(JS): 1.32e-02

In 168 terne su 729 il triangolo è violato. L'esempio stampato è netto: da 0,1 a 0,9 direttamente costa 2,5359 bit, passando per 0,6 solo 1,2435.

La divergenza di Jensen-Shannon resta sotto 1 bit (massimo osservato 0,9993) ed è simmetrica per costruzione. Sulle 200.000 terne la radice quadrata non viola mai il triangolo, con margine minimo 0,013.

Questo è un controllo numerico, non una dimostrazione: che la radice di JS sia una metrica è un risultato noto in letteratura, che qui cito ma non ho letto, e per un uso rigoroso va cercato nella fonte.

Quale verso minimizzare

L'asimmetria cambia il comportamento di chi approssima. Il listato prende una distribuzione bimodale p (due gaussiane in −3 e +3, discretizzate su una griglia con passo 0,1) e cerca, con una ricerca esaustiva su una griglia di medie e deviazioni standard, la gaussiana q che minimizza D(p||q) e quella che minimizza D(q||p).

Come secondo calcolo indipendente confronta con i momenti di p: per la famiglia gaussiana il minimizzatore di D(p||q) ha la stessa media e la stessa varianza di p, e qui la varianza è 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()                                   # distribuzione bimodale discretizzata


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"minimizzando 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"momenti di p: media {media:+.2f}, sigma {sigma:.3f}  (sqrt(1 + 3^2) = {np.sqrt(10):.3f})")

# Output:
# minimizzando D(p||q): mu=+0.0  sigma=3.2  D=0.6643 bit
# minimizzando D(q||p): mu=+3.0  sigma=1.0  D=0.9944 bit
# momenti di p: media +0.00, sigma 3.162  (sqrt(1 + 3^2) = 3.162)

Il verso D(p||q) produce una gaussiana larga, con media 0 e σ = 3,2, che copre entrambi i picchi e ha il suo massimo dove p è quasi zero; i momenti di p (σ = 3,162) confermano il risultato.

Il verso D(q||p) sceglie un solo picco, media +3 e σ = 1 (per simmetria il minimo a −3 vale quanto quello a +3; il listato ne riporta uno).

La spiegazione è nella formula: in D(p||q) il peso è p, quindi q viene punita dove p ha massa e q non ne ha; in D(q||p) il peso è q, e q viene punita dove ha massa e p no, per cui le conviene concentrarsi su un picco.

Il risultato vale per questo giocattolo (due picchi simmetrici, una sola famiglia): non è una regola su quale verso sia «migliore», e la scelta dipende da che cosa costa un errore nell'applicazione.

Stimare l'informazione mutua dai conteggi

Come per l'entropia, sostituire le probabilità con le frequenze osservate introduce una distorsione, e qui ha una forma ancora più insidiosa: variabili indipendenti, per le quali I vale zero, risultano sempre dipendenti. Il listato estrae 1000 coppie indipendenti di due variabili con 10 valori ciascuna, calcola l'informazione mutua dai conteggi, e ripete 2000 volte (seme 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 e Y indipendenti
    stime.append(mi_bit(t))
s = np.array(stime)
print(f"variabili INDIPENDENTI, {k}x{k} valori, N={N}: I vera = 0 bit")
print(f"I stimata dai conteggi: {s.mean():.4f} +/- {s.std(ddof=1) / np.sqrt(prove):.4f} bit")
print(f"ordine di grandezza (k-1)^2 / (2 N ln 2) = {(k - 1) ** 2 / (2 * N * np.log(2)):.4f} bit")

# Output:
# variabili INDIPENDENTI, 10x10 valori, N=1000: I vera = 0 bit
# I stimata dai conteggi: 0.0596 +/- 0.0002 bit
# ordine di grandezza (k-1)^2 / (2 N ln 2) = 0.0584 bit

La stima media è 0,0596 ± 0,0002 bit, non zero. Dello stesso ordine è (k−1)²/(2N ln 2) = 0,0584 bit, un'approssimazione di primo ordine che uso come ordine di grandezza senza averla dimostrata. Il significato pratico è che un valore «piccolo» di informazione mutua stimato su un campione non è un'evidenza di dipendenza: va confrontato con quello che si ottiene su dati mescolati, e, secondo quell'approssimazione, cresce con il numero di valori e decresce con la dimensione del campione.

Errori tipici

• Scambiare i due argomenti della divergenza. D(p||q) vale 0,7370 bit e D(q||p) 0,5310 bit per le stesse due monete.
• Chiamarla distanza. Non è simmetrica e non rispetta il triangolo (168 terne su 729).
• Cambiare unità senza dirlo. Con il logaritmo naturale i valori sono in nat, con base 2 in bit, e dividere per ln 2 passa dai primi ai secondi.
• Dimenticare i supporti. Se q ha un esito a probabilità zero dove p ne ha, D(p||q) è infinita e un'implementazione che somma senza controllare dà inf o un errore.
• Leggere una stima piccola come «indipendenti». Con 1000 campioni e 10×10 valori, variabili indipendenti danno 0,06 bit.

Il filo comune è che il numero che si stampa dipende da una scelta che le formule non mostrano: il verso, l'unità, il supporto, la dimensione del campione. Per ogni divergenza o informazione mutua riportata in un risultato vale la pena scrivere accanto queste quattro cose.

Che cosa sai fare adesso

Sai calcolare l'informazione mutua in tre modi equivalenti e sai che è una divergenza dall'indipendenza, sai dimostrare con Jensen che la divergenza KL non è negativa e che cosa serve perché sia zero, e sai attribuirle un significato in bit sprecati.

Sai anche mostrare con numeri che non è simmetrica, che il triangolo si rompe e che il verso scelto cambia l'approssimazione, e sai che stimare queste grandezze dai conteggi con pochi dati dà valori troppo alti per l'informazione mutua e troppo bassi per l'entropia.

Il passo seguente è tradurre l'entropia in lunghezze di codice: l'articolo 03 dimostra che nessun codice per simbolo scende sotto H(X) e costruisce quello ottimale; l'informazione mutua tornerà nell'articolo 04, dove misura quanto passa in un canale rumoroso. L'elenco degli articoli è nella pagina della serie.

Più avanti nella serie trovi anche l'articolo 05 sui codici correttori e l'articolo 06 sulla compressione con perdita.