06 - Compressione con perdita e rate-distortion: il terzo teorema di Shannon e la quantizzazione

Hai un segnale, una fotografia o i pesi di una rete neurale e sei disposto a perdere un po' di qualità per risparmiare bit. La domanda che si pone subito è quanti bit servono davvero per restare sotto un errore medio D. La risposta di Shannon si chiama funzione rate-distortion R(D): il minimo di informazione mutua fra sorgente e ricostruzione, fra tutti i modi di ricostruire con distorsione media al più D. Per una sorgente gaussiana di varianza σ² e errore quadratico vale R(D) = ½ log2(σ²/D): ogni bit in più per simbolo divide la distorsione per 4, cioè 6,02 dB. I quantizzatori reali non lo raggiungono: un quantizzatore a un bit su una gaussiana ha distorsione 0,363 contro i 0,25 del limite, e uno uniforme con codifica entropica perde 0,2546 bit al limite di alta risoluzione. Misuro tutto questo con calcoli esatti e simulazioni a seme fisso, e chiudo la serie.

Che cosa misura R(D)

Sia X la sorgente, X̂ la sua ricostruzione e d(x, x̂) una funzione che dice quanto costa ricostruire x̂ quando la sorgente ha emesso x. La funzione rate-distortion è R(D) = min I(X; X̂) sulle ricostruzioni con E[d(X, X̂)] ≤ D. Shannon la definisce nel paragrafo 27 dell'articolo del 1948 (come minimo della velocità di trasmissione sotto un vincolo di qualità) e ne dà il significato operativo nel teorema 21: se la sorgente ha velocità R1 per una qualità v1, si può trasmetterla su un canale di capacità C con fedeltà vicina a v1 purché R1 ≤ C, e non si può se R1 > C. Il trattamento completo per sorgenti discrete è nell'articolo del 1959, che cito senza averlo letto.

La definizione è un problema di minimo: scegli la regola probabilistica p(x̂|x) che comprime di più, cioè che lascia meno informazione mutua, senza superare la distorsione concessa. Il valore operativo è che nessun codice con velocità sotto R(D) può restare entro D, mentre con velocità sopra R(D) esistono codici su blocchi abbastanza lunghi. Per D pari alla varianza (o più) basta non mandare nulla e ricostruire la media: R = 0.

La distorsione non deve essere una distanza

Negli appunti si legge spesso che la misura d deve essere simmetrica e rispettare la disuguaglianza triangolare. Shannon dice il contrario: nel paragrafo 26 definisce la valutazione della fedeltà come una media di una funzione ρ(x, y) che ha la natura generale di una «distanza», e precisa in nota che non è una metrica in senso stretto, perché in generale non soddisfa né ρ(x, y) = ρ(y, x) né la disuguaglianza triangolare. Basta che misuri quanto è sgradito ricevere y quando è stato emesso x. Gli esempi che usi in pratica sono due:

• l'errore quadratico, d(x, x̂) = (x - x̂)², per sorgenti continue come campioni audio, pixel e pesi;
• la distorsione di Hamming, che vale 0 se x = x̂ e 1 altrimenti, per sorgenti discrete: la sua media è la frazione di simboli sbagliati.

La scelta della distorsione è una scelta di che cosa ti interessa conservare: con la distorsione di Hamming tutti gli errori pesano uguale, con l'errore quadratico un errore doppio pesa quattro volte. La funzione R(D) cambia con la scelta, e la teoria non dice quale scegliere: lo decide l'uso del segnale ricostruito.

Il caso gaussiano: R(D) = ½ log₂(σ²/D)

Per una sorgente N(0, σ²) con errore quadratico, Polyanskiy e Wu (teorema 26.2) dimostrano R(D) = ½ log+(σ²/D), dove log+ vale zero quando l'argomento è minore di 1. Il listato la verifica due volte: con SymPy sui valori D = 1/4, 1/10, 1/100, e con l'algoritmo di Blahut-Arimoto per la rate-distortion, che minimizza l'informazione mutua su una gaussiana discretizzata con passo 0,05 senza conoscere la formula. Se il passo fosse troppo grosso il confronto sarebbe inquinato dalla discretizzazione, e infatti la tabella riporta lo scarto percentuale.

import math
import numpy as np
import sympy as sp

# 1) Forma chiusa R(D) = (1/2) log2(sigma^2 / D) per D < sigma^2 (Polyanskiy-Wu, teorema 26.2), con SymPy
s2, D = sp.symbols("sigma2 D", positive=True)
R = sp.log(s2 / D, 2) / 2
for d in (sp.Rational(1, 4), sp.Rational(1, 10), sp.Rational(1, 100)):
    val = R.subs({s2: 1, D: d})
    print(f"sigma^2 = 1, D = {d}:  R = {sp.N(val, 6)} bit/simbolo")
# Un bit in piu' per simbolo divide D per 4 (6,02 dB): D(R) = sigma^2 * 2^(-2R)
r = sp.Symbol("r")
print("D(R+1)/D(R) =", sp.simplify(2 ** (-2 * (r + 1)) / 2 ** (-2 * r)), "; in dB:", round(10 * math.log10(4), 4))

# 2) Controllo indipendente: algoritmo di Blahut-Arimoto per R(D) su una gaussiana discretizzata (passo 0.05)
x = np.linspace(-6, 6, 241)
p = np.exp(-x**2 / 2)
p /= p.sum()
d = (x[:, None] - x[None, :]) ** 2  # distorsione quadratica fra ingresso e ricostruzione


def blahut_arimoto_rd(s, iterazioni=3000):
    q = np.full(len(x), 1.0 / len(x))
    for _ in range(iterazioni):
        Q = q[None, :] * np.exp(-s * d)
        Q /= Q.sum(axis=1, keepdims=True)
        q = p @ Q
    Dm = float((p[:, None] * Q * d).sum())
    with np.errstate(divide="ignore", invalid="ignore"):
        Rm = float(np.nansum(p[:, None] * Q * np.log2(Q / q[None, :])))
    return Dm, Rm


print()
print("   s      D (BA)    R (BA)   R formula   scarto")
for s in (1.5, 3.0, 6.0, 25.0, 50.0):
    Dm, Rm = blahut_arimoto_rd(s)
    Rf = 0.5 * math.log2(1.0 / Dm)
    print(f"{s:5.1f}   {Dm:.5f}   {Rm:.5f}   {Rf:.5f}   {abs(Rm - Rf) / Rf * 100:5.2f}%")

# Uscita:
# sigma^2 = 1, D = 1/4:  R = 1.00000 bit/simbolo
# sigma^2 = 1, D = 1/10:  R = 1.66096 bit/simbolo
# sigma^2 = 1, D = 1/100:  R = 3.32193 bit/simbolo
# D(R+1)/D(R) = 1/4 ; in dB: 6.0206
#
#    s      D (BA)    R (BA)   R formula   scarto
#   1.5   0.33333   0.79248   0.79248    0.00%
#   3.0   0.16667   1.29248   1.29248    0.00%
#   6.0   0.08333   1.79248   1.79248    0.00%
#  25.0   0.02000   2.82193   2.82193    0.00%
#  50.0   0.01000   3.32193   3.32193    0.00%

Con σ² = 1, D = 1/4 richiede esattamente 1 bit e D = 1/100 richiede 3,32 bit: ogni decade di D costa ½ log2 10 = 1,66 bit, e ogni bit divide D per 4, cioè 6,02 dB. Nella tabella Blahut-Arimoto e formula coincidono alla quinta cifra, con scarto percentuale 0,00 in tutti i cinque punti: la discretizzazione a passo 0,05 non pesa. Nei punti calcolati la distorsione ottenuta vale 1/(2s), relazione che segue dalla struttura gaussiana del canale di prova e che qui è soltanto un'osservazione numerica.

Un codice perfetto che comprime

La rate-distortion non vale solo per segnali continui. Per una sorgente di bit equiprobabili con distorsione di Hamming, Polyanskiy e Wu (teorema 26.1) danno R(D) = h(p) - h(D) con p = 1/2, cioè 1 - h(D). Il loro «gioco dei 100 bit»: devi riassumere 100 bit casuali in 50. La strategia ingenua, salvare metà e tirare a indovinare l'altra, sbaglia il 25% dei bit; la rate-distortion promette circa l'11%. Il codice di Hamming (7,4), che conosci dall'articolo 5, è un quantizzatore concreto: sostituisci ogni blocco di 7 bit con la parola più vicina e memorizzi solo i 4 bit del messaggio, velocità 4/7.

import math
from fractions import Fraction
from itertools import product
import numpy as np


def h2(p):
    return -p * math.log2(p) - (1 - p) * math.log2(1 - p)


# Un codice perfetto come quantizzatore: ogni blocco di 7 bit e' sostituito dalla parola di Hamming (7,4) piu' vicina.
G = [(1, 0, 0, 0, 1, 1, 0), (0, 1, 0, 0, 1, 0, 1), (0, 0, 1, 0, 0, 1, 1), (0, 0, 0, 1, 1, 1, 1)]
parole = []
for m in product((0, 1), repeat=4):
    parole.append(tuple(sum(m[i] * G[i][j] for i in range(4)) % 2 for j in range(7)))
totale = sum(min(sum(a != b for a, b in zip(v, c)) for c in parole) for v in product((0, 1), repeat=7))
D_h = Fraction(totale, 128 * 7)  # distorsione media per bit su tutti i 128 vettori equiprobabili
print(f"Hamming (7,4) come quantizzatore: D = {totale}/{128 * 7} = {D_h} = {float(D_h):.4f} per bit; velocita' 4/7 = {4 / 7:.4f} bit per bit")
print(f"R(D) di una sorgente Ber(1/2) con questa D: 1 - h(D) = {1 - h2(float(D_h)):.4f}")

# Ancora la stessa velocita' 4/7: salvare 4 bit su 7 e tirare a indovinare gli altri 3 (seme fisso)
rng = np.random.default_rng(20261113)
n = 1_000_000
bit = rng.integers(0, 2, (n, 7))
ipotesi = rng.integers(0, 2, (n, 3))
D_naive = float((bit[:, 4:] != ipotesi).sum()) / (7 * n)
print(f"salvare 4 bit e tirare a indovinare 3: D = {D_naive:.4f} (teorico 3/14 = {3 / 14:.4f})")

# Il gioco dei 100 bit su 50: R(D) = 1 - h(D) con D = 0.11, contro il 25 per cento della strategia ingenua
print(f"R(0.11) = {1 - h2(0.11):.4f} bit/bit; strategia ingenua: meta' salvata, meta' indovinata: D = 0.25")

# Uscita:
# Hamming (7,4) come quantizzatore: D = 112/896 = 1/8 = 0.1250 per bit; velocita' 4/7 = 0.5714 bit per bit
# R(D) di una sorgente Ber(1/2) con questa D: 1 - h(D) = 0.4564
# salvare 4 bit e tirare a indovinare 3: D = 0.2142 (teorico 3/14 = 0.2143)
# R(0.11) = 0.5001 bit/bit; strategia ingenua: meta' salvata, meta' indovinata: D = 0.25

La distorsione è esattamente 1/8 per bit: 112 modifiche su 896 posizioni, perché ogni vettore dista 0 o 1 da una sola parola (il codice è perfetto) e solo 16 vettori su 128 sono già parole. A velocità 4/7 il limite teorico 1 - h(1/8) è 0,4564 bit per bit, quindi il codice di Hamming sta a 0,115 bit dal limite, mentre il quantizzatore ingenuo, che salva 4 bit e ne indovina 3, ha distorsione 0,2143 (3/14, simulato 0,2142). Con 50 bit su 100 il limite di D = 0,11 dà R = 0,5001, che conferma il numero del gioco, e la strategia ingenua dà 0,25. Questa prova mostra che un codice correttore è anche un compressore con perdita, ma non che un (7,4) raggiunga R(D): il limite si avvicina solo con blocchi molto più lunghi.

Reverse water-filling: dove spendere i bit

Se la sorgente ha componenti gaussiane indipendenti con varianze diverse σi², e la distorsione totale D si divide fra loro, il minimo di rate si ottiene con il reverse water-filling: si sceglie un livello λ tale che Σ min(σi², λ) = D, ogni componente con varianza sopra λ riceve distorsione λ, e quelle sotto λ non ricevono bit e restano alla loro varianza; il rate è ½ Σ log+(σi²/λ). È l'esercizio V.21(b) di Polyanskiy e Wu, che lo presentano come controparte del water-filling per la capacità. Attenzione al verso: la distorsione è min(σi², λ), non max, perché una componente più piccola del livello è già «sotto il livello» e non vale la pena codificarla.

import math
from fractions import Fraction
import numpy as np
from scipy.optimize import minimize

# Reverse water-filling: varianze (4, 1, 1/4), distorsione totale 3/2 (Polyanskiy-Wu, esercizio V.21b)
var = [Fraction(4), Fraction(1), Fraction(1, 4)]
D_tot = Fraction(3, 2)
lam = (D_tot - var[2]) / 2  # ipotesi: 1/4 < lam < 1, quindi min(4,lam) + min(1,lam) + min(1/4,lam) = 2 lam + 1/4
assert var[2] < lam < var[1]
D_i = [min(v, lam) for v in var]
assert sum(D_i) == D_tot
R_chiusa = sum(0.5 * math.log2(float(v / lam)) for v in var if v > lam)
print(f"lambda = {lam}; distorsioni per componente = {[str(d) for d in D_i]}; R = {R_chiusa:.6f} bit")


def rate(Dv):
    return sum(0.5 * math.log2(float(v) / d) for v, d in zip(var, Dv))


vincoli = [{"type": "eq", "fun": lambda Dv: sum(Dv) - float(D_tot)}]
limiti = [(1e-6, float(v)) for v in var]
risultato = minimize(rate, x0=[0.5, 0.5, 0.2], bounds=limiti, constraints=vincoli, method="SLSQP", options={"ftol": 1e-12, "maxiter": 500})
print(f"ottimizzatore SLSQP:   distorsioni {np.round(risultato.x, 4).tolist()}  R = {risultato.fun:.6f} bit")

# Uscita:
# lambda = 5/8; distorsioni per componente = ['5/8', '5/8', '1/4']; R = 1.678072 bit
# ottimizzatore SLSQP:   distorsioni [0.625, 0.625, 0.25]  R = 1.678072 bit

Con varianze 4, 1 e 1/4 e distorsione totale 3/2 il livello è λ = 5/8: le prime due componenti ricevono distorsione 5/8 e il rate vale ½ log2(4/(5/8)) + ½ log2(1/(5/8)) = 1,678 bit; la terza, con varianza 1/4 minore di λ, non riceve bit. L'ottimizzatore numerico SLSQP, che non conosce la regola, trova le stesse distorsioni e lo stesso rate. È la logica con cui un codec a trasformata decide quante cifre dare a ogni coefficiente, senza che io abbia verificato i dettagli di un codec in particolare.

Quantizzazione scalare e il prezzo del bit intero

Un quantizzatore scalare a M livelli suddivide l'asse in celle e sostituisce ogni valore con il livello della sua cella. Per un quantizzatore uniforme di passo Δ su una sorgente liscia l'errore in ogni cella è circa uniforme in [-Δ/2, Δ/2], con varianza Δ²/12. Su una sorgente uniforme in un intervallo di ampiezza L con M = 2b livelli, Δ = L/M e D = L²/(12 M²), come conferma il primo calcolo del listato con L = 1 e M = 8 (1/768): ogni bit in più divide D per 4, come per R(D). Per trovare i livelli ottimi su una sorgente non uniforme si usa l'algoritmo di Lloyd: assegna ogni campione al livello più vicino, sposta ogni livello sul centroide delle sue celle e ripeti finché gli assegnamenti non cambiano.

Il listato verifica la formula dell'errore uniforme con integrali esatti e applica Lloyd in tre situazioni: a sei campioni a mano, alla gaussiana con integrali esatti (senza campioni, per 2, 4, 8 e 16 livelli) e a un campione di 200.000 punti a seme fisso. Inoltre calcola in forma chiusa il quantizzatore a un bit con SymPy e misura, per tre passi, quanto un quantizzatore uniforme seguito da codifica entropica si discosta da R(D): la grandezza H(Q) è l'entropia dei livelli, cioè i bit medi che la codifica entropica spende.

import math
import numpy as np
import sympy as sp
from scipy import integrate
from scipy.stats import norm

# 1) Un bit su una gaussiana N(0,1): il centroide di x > 0 e la distorsione, in forma chiusa con SymPy
x = sp.symbols("x", real=True)
fx = sp.exp(-x**2 / 2) / sp.sqrt(2 * sp.pi)
centroide = sp.simplify(2 * sp.integrate(x * fx, (x, 0, sp.oo)))  # E[X | X > 0], perche' P(X > 0) = 1/2
dist = sp.simplify(1 - centroide**2)  # D = E[X^2] - c^2 con livelli +c e -c
print("livelli: +-", centroide, "=", sp.N(centroide, 6), "; D =", dist, "=", sp.N(dist, 6))
print("R(D) alla stessa distorsione:", round(0.5 * math.log2(1 / float(dist)), 4), "bit contro 1 bit del quantizzatore")

# 1-bis) Sorgente uniforme su [0, 1] con 8 livelli: D = L^2 / (12 M^2), con integrali esatti cella per cella
L, M = sp.Integer(1), 8
celle = [(L * i / M, L * (i + 1) / M) for i in range(M)]
D_unif = sum(sp.integrate((x - (a + b) / 2) ** 2, (x, a, b)) for a, b in celle)
print("sorgente uniforme su [0,1], 8 livelli: D =", D_unif, "; L^2/(12 M^2) =", L**2 / (12 * M**2))


# 2) Lloyd-Max "esatto" (integrali della gaussiana, nessun campione) per M = 2, 4, 8, 16 livelli
def lloyd_max_esatto(M, iterazioni=2000):
    livelli = np.linspace(-2, 2, M)
    for _ in range(iterazioni):
        soglie = np.concatenate(([-np.inf], (livelli[:-1] + livelli[1:]) / 2, [np.inf]))
        a, b = soglie[:-1], soglie[1:]
        massa = norm.cdf(b) - norm.cdf(a)
        livelli = (norm.pdf(a) - norm.pdf(b)) / massa  # centroide di una gaussiana troncata
    soglie = np.concatenate(([-np.inf], (livelli[:-1] + livelli[1:]) / 2, [np.inf]))
    Dm = 0.0
    for lv, a, b in zip(livelli, soglie[:-1], soglie[1:]):
        Dm += integrate.quad(lambda t: (t - lv) ** 2 * norm.pdf(t), a, b)[0]
    return Dm


print()
print("bit  livelli  D (Lloyd-Max)  D = 2^-2R (Shannon)  scarto in dB")
for bit in (1, 2, 3, 4):
    Dlm = lloyd_max_esatto(2**bit)
    Dsh = 2.0 ** (-2 * bit)
    print(f"{bit:3d}  {2**bit:7d}  {Dlm:13.5f}  {Dsh:19.5f}  {10 * math.log10(Dlm / Dsh):11.2f}")

# 3) Quantizzatore uniforme + codifica entropica: il guadagno asintotico di R(D) e' (1/2) log2(pi e / 6) = 0,2546 bit
print()
print("passo   D (esatta)   H(Q) bit   R(D) bit   H(Q) - R(D)")
for passo in (0.5, 0.25, 0.1):
    estremi = np.arange(-8, 8 + passo / 2, passo)
    centri = (estremi[:-1] + estremi[1:]) / 2
    prob = norm.cdf(estremi[1:]) - norm.cdf(estremi[:-1])
    H = -sum(q * math.log2(q) for q in prob if q > 0)
    Dq = sum(integrate.quad(lambda t: (t - c) ** 2 * norm.pdf(t), a, b)[0] for c, a, b in zip(centri, estremi[:-1], estremi[1:]))
    Rd = 0.5 * math.log2(1 / Dq)
    print(f"{passo:5.2f}   {Dq:.6f}   {H:8.4f}   {Rd:8.4f}   {H - Rd:10.4f}")
print("limite 0.5 * log2(pi e / 6) =", round(0.5 * math.log2(math.pi * math.e / 6), 4))

# 4) L'algoritmo di Lloyd su 6 campioni (passi a mano), con livelli iniziali -0.5 e +0.5
campioni = np.array([-1.2, -0.3, 0.1, 0.5, 0.8, 1.1])
mu = np.array([-0.5, 0.5])
print()
precedente = None
for it in range(1, 6):
    scelta = np.argmin((campioni[:, None] - mu[None, :]) ** 2, axis=1)
    if precedente is not None and np.array_equal(scelta, precedente):
        print(f"iterazione {it}: gli assegnamenti non cambiano, l'algoritmo si ferma")
        break
    precedente = scelta
    mu = np.array([campioni[scelta == j].mean() for j in range(2)])
    Dtot = float(((campioni - mu[scelta]) ** 2).mean())
    print(f"iterazione {it}: gruppi {[campioni[scelta == j].tolist() for j in range(2)]}  livelli {np.round(mu, 4).tolist()}  D medio {Dtot:.4f}")

# 5) Lloyd su un campione gaussiano a seme fisso, 4 livelli, contro il valore esatto
rng = np.random.default_rng(20261113)
campione = rng.standard_normal(200_000)
mu = np.array([-1.5, -0.5, 0.5, 1.5])
for _ in range(100):
    scelta = np.argmin((campione[:, None] - mu[None, :]) ** 2, axis=1)
    mu = np.array([campione[scelta == j].mean() for j in range(4)])
scelta = np.argmin((campione[:, None] - mu[None, :]) ** 2, axis=1)
Dc = float(((campione - mu[scelta]) ** 2).mean())
print(f"Lloyd su 200000 campioni, 4 livelli: {np.round(mu, 4).tolist()}  D = {Dc:.5f}  (esatto {lloyd_max_esatto(4):.5f})")

# Uscita:
# livelli: +- sqrt(2)/sqrt(pi) = 0.797885 ; D = (-2 + pi)/pi = 0.363380
# R(D) alla stessa distorsione: 0.7302 bit contro 1 bit del quantizzatore
# sorgente uniforme su [0,1], 8 livelli: D = 1/768 ; L^2/(12 M^2) = 1/768
#
# bit  livelli  D (Lloyd-Max)  D = 2^-2R (Shannon)  scarto in dB
#   1        2        0.36338              0.25000         1.62
#   2        4        0.11748              0.06250         2.74
#   3        8        0.03455              0.01562         3.45
#   4       16        0.00950              0.00391         3.86
#
# passo   D (esatta)   H(Q) bit   R(D) bit   H(Q) - R(D)
#  0.50   0.020833     3.0620     2.7925       0.2695
#  0.25   0.005208     4.0508     3.7925       0.2584
#  0.10   0.000833     5.3696     5.1144       0.2552
# limite 0.5 * log2(pi e / 6) = 0.2546
#
# iterazione 1: gruppi [[-1.2, -0.3], [0.1, 0.5, 0.8, 1.1]]  livelli [-0.75, 0.625]  D medio 0.1588
# iterazione 2: gli assegnamenti non cambiano, l'algoritmo si ferma
# Lloyd su 200000 campioni, 4 livelli: [-1.5152, -0.454, 0.4543, 1.5121]  D = 0.11736  (esatto 0.11748)

Il quantizzatore a un bit ha livelli ±√(2/π) = ±0,7979 e distorsione 1 - 2/π = 0,3634, valore che Polyanskiy e Wu riportano come (1 - 2/π)σ² ≈ 0,36σ². A quella distorsione R(D) vale 0,7302 bit, non 1: la differenza non si recupera con un quantizzatore a livelli fissi. La tabella per 2, 4, 8 e 16 livelli mostra un divario con il limite 2-2R che cresce da 1,62 a 3,86 dB. Con codifica entropica il divario in bit tende a ½ log2(πe/6) = 0,2546 bit, circa 1,53 dB: 0,2695, 0,2584 e 0,2552 per passi 0,5, 0,25 e 0,1. Il limite si ricava confrontando H(Q) ≈ h(X) - log2 Δ con R(D) a D = Δ²/12: è un mio calcolo, controllato qui con gli integrali esatti.

Sui sei campioni l'algoritmo converge dopo la prima iterazione (la seconda non cambia gli assegnamenti e l'algoritmo si ferma): 0,1 viene assegnato al livello +0,5 (distanza 0,4, contro 0,6 dal livello -0,5), non a quello negativo. I livelli finali sono −0,75 e 0,625, con distorsione media 0,1588. Con 200.000 campioni e quattro livelli, Lloyd arriva a D = 0,11736 contro il valore esatto 0,11748: lo scarto di 0,0001 è compatibile con l'errore di campionamento di un campione di quella dimensione, che non ho misurato a parte; Lloyd converge comunque a un minimo locale, non necessariamente globale.

I pesi di una rete: che cosa dice la teoria e che cosa no

Quantizzare i pesi di un modello è lo stesso problema: b bit per peso, distorsione quadratica media. Il listato usa pesi sintetici estratti da N(0, 0,02²) (non i pesi di una rete vera) e un quantizzatore uniforme simmetrico con una scala unica, scelta in due modi: dal massimo assoluto dei pesi, oppure ottimizzando la scala per minimizzare l'errore, che equivale a tagliare le code.

import math
import numpy as np

rng = np.random.default_rng(20261113)
# Quantizzazione di pesi SINTETICI N(0, 0.02^2) a 8 bit, scala per tensore: errore e confronto con il limite di Shannon
sigma = 0.02
w = rng.normal(0.0, sigma, 1_000_000)
print("bit  scala           MSE/sigma^2     SNR dB   limite Shannon dB")
for b in (4, 6, 8):
    qmax = 2 ** (b - 1) - 1
    scala_max = np.abs(w).max() / qmax
    candidati = np.linspace(0.5, 1.0, 51) * scala_max
    errori = [np.mean((w - np.clip(np.round(w / s), -qmax - 1, qmax) * s) ** 2) for s in candidati]
    s_best = candidati[int(np.argmin(errori))]
    for nome, s in (("max|w|", scala_max), ("ottimizzata", s_best)):
        mse = float(np.mean((w - np.clip(np.round(w / s), -qmax - 1, qmax) * s) ** 2))
        print(f"{b:3d}  {nome:11s}  {mse / sigma**2:12.3e}  {10 * math.log10(sigma**2 / mse):9.2f}  {10 * math.log10(4**b):13.2f}")
print("memoria: 32 bit -> 8 bit =", 32 // 8, "volte meno")

# Uscita:
# bit  scala           MSE/sigma^2     SNR dB   limite Shannon dB
#   4  max|w|          4.028e-02      13.95          24.08
#   4  ottimizzata     1.200e-02      19.21          24.08
#   6  max|w|          2.054e-03      26.87          36.12
#   6  ottimizzata     1.066e-03      29.72          36.12
#   8  max|w|          1.223e-04      39.13          48.16
#   8  ottimizzata     8.972e-05      40.47          48.16
# memoria: 32 bit -> 8 bit = 4 volte meno

A 8 bit, con scala dal massimo, il rapporto segnale-rumore vale 39,13 dB; con scala ottimizzata 40,47 dB; il limite di Shannon per pesi gaussiani con 8 bit fissi è 48,16 dB, cioè 6,02 · 8. A 4 bit i tre numeri sono 13,95, 19,21 e 24,08 dB: ottimizzare la scala vale più di 5 dB, che è quasi un bit. Il passaggio da 32 a 8 bit riduce la memoria di quattro volte per definizione. Non so, e la teoria non dice, quanta accuratezza perda un modello vero: dipende dalla rete e dai dati, e si misura sul modello. Il teorema 26.3 di Polyanskiy e Wu dà un altro dato utile: con media nulla, stessa varianza e errore quadratico, R(D) ≤ ½ log+(σ²/D) per ogni sorgente, quindi pesi non gaussiani sono almeno altrettanto comprimibili, e il limite gaussiano è prudente.

Dove si incontrano rate-distortion e capacità

La rate-distortion non è un argomento separato dalla capacità: il teorema 21 di Shannon dice che puoi trasmettere con fedeltà v1 su un canale di capacità C se e solo se R1 ≤ C, e Polyanskiy e Wu (teorema 26.6) lo riscrivono come velocità massima di codifica congiunta C / R(D). Il caso in cui la sorgente è un bit equiprobabile con distorsione di Hamming è la parte 2 del teorema dell'articolo 4: la velocità raggiungibile con un errore di bit pb è C / (1 - h(pb)), cioè C / R(pb) con R(pb) = 1 - h(pb) (osservazione 26.4 di Polyanskiy e Wu). Quindi il limite sopra la capacità, che l'articolo 4 chiama «accettare errori di bit», è un caso della compressione con perdita.

Nel machine learning la stessa forma compare in un altro posto. Alemi e altri (Fixing a Broken ELBO, arXiv 1711.00464) derivano limiti variazionali sull'informazione mutua fra ingresso e variabile latente di un autoencoder variazionale e ne ricavano una curva rate-distortion che descrive il compromesso fra compressione e accuratezza della ricostruzione. Hanno anche mostrato che esistono modelli con lo stesso ELBO e proprietà diverse: ho letto l'abstract, non l'articolo. Con la definizione standard dell'ELBO come ricostruzione meno KL fra posteriore e prior (la divergenza dell'articolo 2), il secondo termine fa da misura del «rate», e il primo da «distorsione».

Errori frequenti

• Pretendere che la distorsione sia una metrica. Shannon richiede soltanto una funzione che misuri quanto costa un errore; la simmetria non serve.
• Scambiare il segno nel water-filling inverso. La distorsione di una componente è min(σi², λ), non il massimo.
• Dimenticare che il limite è asintotico. R(D) vale per blocchi lunghi; con un bit per peso, il quantizzatore ottimo a un bit resta a 1,62 dB dal limite.
• Dedurre l'accuratezza di un modello dal solo errore quadratico dei pesi. La teoria quantifica la distorsione dei pesi; l'effetto sull'uscita della rete si misura sul modello.
• Dimenticare che la quantizzazione non è invertibile. Una volta quantizzato non recuperi l'originale, a differenza della compressione senza perdita dell'articolo 3.

Il tratto comune dei primi due e del quarto è dare a una formula un campo di validità più largo di quello che ha: il primo errore la restringe (metrica), il secondo ne sbaglia l'ordine (min), il quarto la estende dall'errore sui pesi all'accuratezza. Prima di applicarla, scrivi l'ipotesi che la rende vera.

Che cosa sai fare adesso, e che cosa hai in mano della serie

Sai calcolare R(D) per una sorgente gaussiana con errore quadratico e per una sorgente binaria con distorsione di Hamming, riconoscere che un quantizzatore a livelli fissi spreca da 1,6 a 3,9 dB e uno uniforme con codifica entropica 0,2546 bit, e allocare i bit fra componenti con il reverse water-filling. Sai anche collegare la compressione con perdita alla capacità. Con questo articolo la serie è completa. Si parte dall'entropia di Shannon, misura di incertezza di una sorgente; si passa a informazione mutua e divergenza di Kullback-Leibler, che servono a confrontare distribuzioni e a misurare ciò che un'osservazione rivela; si arriva alla codifica di sorgente con Huffman e il primo teorema di Shannon, cioè la compressione senza perdita; si affronta il canale rumoroso e la capacità, che fissa quanto si può comunicare in modo affidabile; si costruiscono i codici correttori con distanza di Hamming e i limiti di Hamming e Singleton; e si chiude qui con la compressione con perdita. L'indice della serie raccoglie tutti e sei gli articoli.

Domande frequenti