06 - Lossy compression and rate-distortion: Shannon's third theorem and model quantization

You have a signal, a photograph or the weights of a neural network and you are willing to lose a little quality to save bits. The question that comes up at once is how many bits you really need to stay under a mean error D. Shannon's answer is called the rate-distortion function R(D): the minimum of the mutual information between source and reconstruction, over all the ways of reconstructing with mean distortion at most D. For a Gaussian source of variance σ² and squared error it is R(D) = ½ log2(σ²/D): every extra bit per symbol divides the distortion by 4, that is 6.02 dB. Real quantizers do not reach it: a one-bit quantizer on a Gaussian has distortion 0.363 against the limit's 0.25, and a uniform one with entropy coding loses 0.2546 bits at the high-resolution limit. I measure all this with exact calculations and fixed-seed simulations, and I close the series.

What R(D) measures

Let X be the source, X̂ its reconstruction and d(x, x̂) a function that says how much it costs to reconstruct x̂ when the source emitted x. The rate-distortion function is R(D) = min I(X; X̂) over the reconstructions with E[d(X, X̂)] ≤ D. Shannon defines it in section 27 of the 1948 paper (as the minimum of the transmission rate under a quality constraint) and gives its operational meaning in theorem 21: if the source has rate R1 for a quality v1, it can be transmitted over a channel of capacity C with fidelity close to v1 provided R1 ≤ C, and it cannot if R1 > C. The full treatment for discrete sources is in the 1959 paper, which I cite without having read it.

The definition is a minimization problem: you choose the probabilistic rule p(x̂|x) that compresses the most, that is that leaves the least mutual information, without exceeding the allowed distortion. The operational value is that no code with rate below R(D) can stay within D, while with rate above R(D) there exist codes on blocks long enough. For D equal to the variance (or more) it is enough to send nothing and reconstruct the mean: R = 0.

Distortion need not be a distance

In lecture notes you often read that the measure d must be symmetric and respect the triangle inequality. Shannon says the opposite: in section 26 he defines the evaluation of fidelity as an average of a function ρ(x, y) that has the general nature of a "distance", and he states in a footnote that it is not a "metric" in the strict sense, because in general it satisfies neither ρ(x, y) = ρ(y, x) nor the triangle inequality. It is enough that it measures how bad it is to receive y when x was emitted. The two examples you use in practice are:

• the squared error, d(x, x̂) = (x - x̂)², for continuous sources such as audio samples, pixels and weights;
• the Hamming distortion, which is 0 if x = x̂ and 1 otherwise, for discrete sources: its mean is the fraction of wrong symbols.

The choice of distortion is a choice of what you care about preserving: with the Hamming distortion all errors weigh the same, with the squared error a double error weighs four times as much. The function R(D) changes with the choice, and the theory does not say which to choose: the use of the reconstructed signal decides.

The Gaussian case: R(D) = ½ log₂(σ²/D)

For a source N(0, σ²) with squared error, Polyanskiy and Wu (theorem 26.2) prove R(D) = ½ log+(σ²/D), where log+ is zero when the argument is less than 1. The listing checks it twice: with SymPy on the values D = 1/4, 1/10, 1/100, and with the Blahut-Arimoto algorithm for rate-distortion, which minimizes the mutual information on a Gaussian discretized with step 0.05 without knowing the formula. If the step were too coarse the comparison would be polluted by the discretization, and indeed the table reports the percentage gap.

import math
import numpy as np
import sympy as sp

# 1) Closed form R(D) = (1/2) log2(sigma^2 / D) for D < sigma^2 (Polyanskiy-Wu, theorem 26.2), with 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/symbol")
# One more bit per symbol divides D by 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) Independent check: Blahut-Arimoto algorithm for R(D) on a discretized Gaussian (step 0.05)
x = np.linspace(-6, 6, 241)
p = np.exp(-x**2 / 2)
p /= p.sum()
d = (x[:, None] - x[None, :]) ** 2  # squared distortion between input and reconstruction


def blahut_arimoto_rd(s, iterations=3000):
    q = np.full(len(x), 1.0 / len(x))
    for _ in range(iterations):
        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   gap")
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}%")

# Output:
# sigma^2 = 1, D = 1/4:  R = 1.00000 bit/symbol
# sigma^2 = 1, D = 1/10:  R = 1.66096 bit/symbol
# sigma^2 = 1, D = 1/100:  R = 3.32193 bit/symbol
# D(R+1)/D(R) = 1/4 ; in dB: 6.0206
#
#    s      D (BA)    R (BA)   R formula   gap
#   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%

With σ² = 1, D = 1/4 requires exactly 1 bit and D = 1/100 requires 3.32 bits: every decade of D costs ½ log2 10 = 1.66 bits, and every bit divides D by 4, that is 6.02 dB. In the table Blahut-Arimoto and the formula coincide to the fifth digit, with percentage gap 0.00 at all five points: the discretization at step 0.05 does not weigh. At the computed points the distortion obtained is 1/(2s), a relation that follows from the Gaussian structure of the test channel and that here is only a numerical observation.

A perfect code that compresses

Rate-distortion does not hold only for continuous signals. For a source of equiprobable bits with Hamming distortion, Polyanskiy and Wu (theorem 26.1) give R(D) = h(p) - h(D) with p = 1/2, that is 1 - h(D). Their "game of 100 bits": you must summarize 100 random bits in 50. The naive strategy, saving half and guessing the other half, gets 25% of the bits wrong; rate-distortion promises about 11%. The Hamming (7,4) code, which you know from article 5, is a concrete quantizer: you replace every block of 7 bits with the nearest codeword and store only the 4 message bits, rate 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)


# A perfect code as a quantizer: every block of 7 bits is replaced by the nearest Hamming (7,4) codeword.
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)]
words = []
for m in product((0, 1), repeat=4):
    words.append(tuple(sum(m[i] * G[i][j] for i in range(4)) % 2 for j in range(7)))
total = sum(min(sum(a != b for a, b in zip(v, c)) for c in words) for v in product((0, 1), repeat=7))
D_h = Fraction(total, 128 * 7)  # mean distortion per bit over all 128 equiprobable vectors
print(f"Hamming (7,4) as a quantizer: D = {total}/{128 * 7} = {D_h} = {float(D_h):.4f} per bit; rate 4/7 = {4 / 7:.4f} bits per bit")
print(f"R(D) of a Ber(1/2) source with this D: 1 - h(D) = {1 - h2(float(D_h)):.4f}")

# Same rate 4/7 again: save 4 bits out of 7 and guess the other 3 (fixed seed)
rng = np.random.default_rng(20261113)
n = 1_000_000
bit = rng.integers(0, 2, (n, 7))
guesses = rng.integers(0, 2, (n, 3))
D_naive = float((bit[:, 4:] != guesses).sum()) / (7 * n)
print(f"save 4 bits and guess 3: D = {D_naive:.4f} (theoretical 3/14 = {3 / 14:.4f})")

# The game of 100 bits in 50: R(D) = 1 - h(D) with D = 0.11, against the 25 percent of the naive strategy
print(f"R(0.11) = {1 - h2(0.11):.4f} bit/bit; naive strategy: half saved, half guessed: D = 0.25")

# Output:
# Hamming (7,4) as a quantizer: D = 112/896 = 1/8 = 0.1250 per bit; rate 4/7 = 0.5714 bits per bit
# R(D) of a Ber(1/2) source with this D: 1 - h(D) = 0.4564
# save 4 bits and guess 3: D = 0.2142 (theoretical 3/14 = 0.2143)
# R(0.11) = 0.5001 bit/bit; naive strategy: half saved, half guessed: D = 0.25

The distortion is exactly 1/8 per bit: 112 changes over 896 positions, because every vector is at distance 0 or 1 from exactly one codeword (the code is perfect) and only 16 vectors out of 128 are already codewords. At rate 4/7 the theoretical limit 1 - h(1/8) is 0.4564 bits per bit, so the Hamming code is 0.115 bits from the limit, while the naive quantizer, which saves 4 bits and guesses 3, has distortion 0.2143 (3/14, simulated 0.2142). With 50 bits out of 100 the limit at D = 0.11 gives R = 0.5001, which confirms the number of the game, and the naive strategy gives 0.25. This test shows that an error-correcting code is also a lossy compressor, but not that a (7,4) reaches R(D): the limit is approached only with much longer blocks.

Reverse water-filling: where to spend the bits

If the source has independent Gaussian components with different variances σi², and the total distortion D is divided among them, the minimum rate is obtained with reverse water-filling: you choose a level λ such that Σ min(σi², λ) = D, every component with variance above λ receives distortion λ, and those below λ receive no bits and stay at their variance; the rate is ½ Σ log+(σi²/λ). It is exercise V.21(b) of Polyanskiy and Wu, who present it as the counterpart of water-filling for capacity. Mind the direction: the distortion is min(σi², λ), not max, because a component smaller than the level is already "below the level" and it is not worth encoding.

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

# Reverse water-filling: variances (4, 1, 1/4), total distortion 3/2 (Polyanskiy-Wu, exercise V.21b)
var = [Fraction(4), Fraction(1), Fraction(1, 4)]
D_tot = Fraction(3, 2)
lam = (D_tot - var[2]) / 2  # assumption: 1/4 < lam < 1, hence 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_closed = sum(0.5 * math.log2(float(v / lam)) for v in var if v > lam)
print(f"lambda = {lam}; distortions per component = {[str(d) for d in D_i]}; R = {R_closed:.6f} bits")


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


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

# Output:
# lambda = 5/8; distortions per component = ['5/8', '5/8', '1/4']; R = 1.678072 bits
# SLSQP optimizer:       distortions [0.625, 0.625, 0.25]  R = 1.678072 bits

With variances 4, 1 and 1/4 and total distortion 3/2 the level is λ = 5/8: the first two components receive distortion 5/8 and the rate is ½ log2(4/(5/8)) + ½ log2(1/(5/8)) = 1.678 bits; the third, with variance 1/4 less than λ, receives no bits. The numerical SLSQP optimizer, which does not know the rule, finds the same distortions and the same rate. It is the logic by which a transform codec decides how many digits to give to each coefficient, without my having verified the details of any particular codec.

Scalar quantization and the price of the whole bit

A scalar quantizer with M levels divides the axis into cells and replaces every value with the level of its cell. For a uniform quantizer of step Δ on a smooth source the error in each cell is roughly uniform in [-Δ/2, Δ/2], with variance Δ²/12. On a source uniform over an interval of width L with M = 2b levels, Δ = L/M and D = L²/(12 M²), as the first calculation of the listing confirms with L = 1 and M = 8 (1/768): every extra bit divides D by 4, as for R(D). To find the optimal levels on a non-uniform source you use Lloyd's algorithm: assign every sample to the nearest level, move every level to the centroid of its cells and repeat until the assignments no longer change.

The listing checks the formula of the uniform error with exact integrals and applies Lloyd in three situations: to six samples by hand, to the Gaussian with exact integrals (no samples, for 2, 4, 8 and 16 levels) and to a sample of 200,000 points with a fixed seed. It also computes in closed form the one-bit quantizer with SymPy and measures, for three steps, how far a uniform quantizer followed by entropy coding is from R(D): the quantity H(Q) is the entropy of the levels, that is the mean bits that entropy coding spends.

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

# 1) One bit on a Gaussian N(0,1): the centroid of x > 0 and the distortion, in closed form with SymPy
x = sp.symbols("x", real=True)
fx = sp.exp(-x**2 / 2) / sp.sqrt(2 * sp.pi)
centroid = sp.simplify(2 * sp.integrate(x * fx, (x, 0, sp.oo)))  # E[X | X > 0], because P(X > 0) = 1/2
dist = sp.simplify(1 - centroid**2)  # D = E[X^2] - c^2 with levels +c and -c
print("levels: +-", centroid, "=", sp.N(centroid, 6), "; D =", dist, "=", sp.N(dist, 6))
print("R(D) at the same distortion:", round(0.5 * math.log2(1 / float(dist)), 4), "bits against 1 bit of the quantizer")

# 1-bis) Uniform source on [0, 1] with 8 levels: D = L^2 / (12 M^2), with exact integrals cell by cell
L, M = sp.Integer(1), 8
cells = [(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 cells)
print("uniform source on [0,1], 8 levels: D =", D_unif, "; L^2/(12 M^2) =", L**2 / (12 * M**2))


# 2) "Exact" Lloyd-Max (Gaussian integrals, no samples) for M = 2, 4, 8, 16 levels
def lloyd_max_exact(M, iterations=2000):
    levels = np.linspace(-2, 2, M)
    for _ in range(iterations):
        thresholds = np.concatenate(([-np.inf], (levels[:-1] + levels[1:]) / 2, [np.inf]))
        a, b = thresholds[:-1], thresholds[1:]
        mass = norm.cdf(b) - norm.cdf(a)
        levels = (norm.pdf(a) - norm.pdf(b)) / mass  # centroid of a truncated Gaussian
    thresholds = np.concatenate(([-np.inf], (levels[:-1] + levels[1:]) / 2, [np.inf]))
    Dm = 0.0
    for lv, a, b in zip(levels, thresholds[:-1], thresholds[1:]):
        Dm += integrate.quad(lambda t: (t - lv) ** 2 * norm.pdf(t), a, b)[0]
    return Dm


print()
print("bit  levels   D (Lloyd-Max)  D = 2^-2R (Shannon)  gap in dB")
for bit in (1, 2, 3, 4):
    Dlm = lloyd_max_exact(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) Uniform quantizer + entropy coding: the asymptotic gap from R(D) is (1/2) log2(pi e / 6) = 0.2546 bits
print()
print("step    D (exact)    H(Q) bit   R(D) bit   H(Q) - R(D)")
for step in (0.5, 0.25, 0.1):
    edges = np.arange(-8, 8 + step / 2, step)
    centers = (edges[:-1] + edges[1:]) / 2
    prob = norm.cdf(edges[1:]) - norm.cdf(edges[:-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(centers, edges[:-1], edges[1:]))
    Rd = 0.5 * math.log2(1 / Dq)
    print(f"{step:5.2f}   {Dq:.6f}   {H:8.4f}   {Rd:8.4f}   {H - Rd:10.4f}")
print("limit 0.5 * log2(pi e / 6) =", round(0.5 * math.log2(math.pi * math.e / 6), 4))

# 4) Lloyd's algorithm on 6 samples (steps by hand), with initial levels -0.5 and +0.5
samples = np.array([-1.2, -0.3, 0.1, 0.5, 0.8, 1.1])
mu = np.array([-0.5, 0.5])
print()
previous = None
for it in range(1, 6):
    choice = np.argmin((samples[:, None] - mu[None, :]) ** 2, axis=1)
    if previous is not None and np.array_equal(choice, previous):
        print(f"iteration {it}: the assignments do not change, the algorithm stops")
        break
    previous = choice
    mu = np.array([samples[choice == j].mean() for j in range(2)])
    Dtot = float(((samples - mu[choice]) ** 2).mean())
    print(f"iteration {it}: groups {[samples[choice == j].tolist() for j in range(2)]}  levels {np.round(mu, 4).tolist()}  mean D {Dtot:.4f}")

# 5) Lloyd on a Gaussian sample with a fixed seed, 4 levels, against the exact value
rng = np.random.default_rng(20261113)
sample = rng.standard_normal(200_000)
mu = np.array([-1.5, -0.5, 0.5, 1.5])
for _ in range(100):
    choice = np.argmin((sample[:, None] - mu[None, :]) ** 2, axis=1)
    mu = np.array([sample[choice == j].mean() for j in range(4)])
choice = np.argmin((sample[:, None] - mu[None, :]) ** 2, axis=1)
Dc = float(((sample - mu[choice]) ** 2).mean())
print(f"Lloyd on 200000 samples, 4 levels: {np.round(mu, 4).tolist()}  D = {Dc:.5f}  (exact {lloyd_max_exact(4):.5f})")

# Output:
# levels: +- sqrt(2)/sqrt(pi) = 0.797885 ; D = (-2 + pi)/pi = 0.363380
# R(D) at the same distortion: 0.7302 bits against 1 bit of the quantizer
# uniform source on [0,1], 8 levels: D = 1/768 ; L^2/(12 M^2) = 1/768
#
# bit  levels   D (Lloyd-Max)  D = 2^-2R (Shannon)  gap 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
#
# step    D (exact)    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
# limit 0.5 * log2(pi e / 6) = 0.2546
#
# iteration 1: groups [[-1.2, -0.3], [0.1, 0.5, 0.8, 1.1]]  levels [-0.75, 0.625]  mean D 0.1588
# iteration 2: the assignments do not change, the algorithm stops
# Lloyd on 200000 samples, 4 levels: [-1.5152, -0.454, 0.4543, 1.5121]  D = 0.11736  (exact 0.11748)

The one-bit quantizer has levels ±√(2/π) = ±0.7979 and distortion 1 - 2/π = 0.3634, a value that Polyanskiy and Wu report as (1 - 2/π)σ² ≈ 0.36σ². At that distortion R(D) is 0.7302 bits, not 1: the difference cannot be recovered with a quantizer with fixed levels. The table for 2, 4, 8 and 16 levels shows a gap from the limit 2-2R that grows from 1.62 to 3.86 dB. With entropy coding the gap in bits tends to ½ log2(πe/6) = 0.2546 bits, about 1.53 dB: 0.2695, 0.2584 and 0.2552 for steps 0.5, 0.25 and 0.1. The limit is obtained by comparing H(Q) ≈ h(X) - log2 Δ with R(D) at D = Δ²/12: it is my own calculation, checked here with the exact integrals.

On the six samples the algorithm converges after the first iteration (the second changes no assignment and the algorithm stops): 0.1 is assigned to the level +0.5 (distance 0.4, against 0.6 from the level -0.5), not to the negative one. The final levels are −0.75 and 0.625, with mean distortion 0.1588. With 200,000 samples and four levels, Lloyd reaches D = 0.11736 against the exact value 0.11748: the gap of 0.0001 is compatible with the sampling error of a sample of that size, which I did not measure separately; Lloyd in any case converges to a local minimum, not necessarily a global one.

The weights of a network: what the theory says and what it does not

Quantizing the weights of a model is the same problem: b bits per weight, mean squared distortion. The listing uses synthetic weights drawn from N(0, 0.02²) (not the weights of a real network) and a symmetric uniform quantizer with a single scale, chosen in two ways: from the absolute maximum of the weights, or by optimizing the scale to minimize the error, which amounts to clipping the tails.

import math
import numpy as np

rng = np.random.default_rng(20261113)
# Quantization of SYNTHETIC weights N(0, 0.02^2) at 4, 6 and 8 bits, one scale per tensor: error and comparison with the Shannon limit
sigma = 0.02
w = rng.normal(0.0, sigma, 1_000_000)
print("bit  scale           MSE/sigma^2     SNR dB   Shannon limit dB")
for b in (4, 6, 8):
    qmax = 2 ** (b - 1) - 1
    scale_max = np.abs(w).max() / qmax
    candidates = np.linspace(0.5, 1.0, 51) * scale_max
    errors = [np.mean((w - np.clip(np.round(w / s), -qmax - 1, qmax) * s) ** 2) for s in candidates]
    s_best = candidates[int(np.argmin(errors))]
    for name, s in (("max|w|", scale_max), ("optimized", s_best)):
        mse = float(np.mean((w - np.clip(np.round(w / s), -qmax - 1, qmax) * s) ** 2))
        print(f"{b:3d}  {name:11s}  {mse / sigma**2:12.3e}  {10 * math.log10(sigma**2 / mse):9.2f}  {10 * math.log10(4**b):13.2f}")
print("memory: 32 bit -> 8 bit =", 32 // 8, "times less")

# Output:
# bit  scale           MSE/sigma^2     SNR dB   Shannon limit dB
#   4  max|w|          4.028e-02      13.95          24.08
#   4  optimized       1.200e-02      19.21          24.08
#   6  max|w|          2.054e-03      26.87          36.12
#   6  optimized       1.066e-03      29.72          36.12
#   8  max|w|          1.223e-04      39.13          48.16
#   8  optimized       8.972e-05      40.47          48.16
# memory: 32 bit -> 8 bit = 4 times less

At 8 bits, with the scale from the maximum, the signal-to-noise ratio is 39.13 dB; with the optimized scale 40.47 dB; Shannon's limit for Gaussian weights with 8 fixed bits is 48.16 dB, that is 6.02 · 8. At 4 bits the three numbers are 13.95, 19.21 and 24.08 dB: optimizing the scale is worth more than 5 dB, which is almost one bit. The move from 32 to 8 bits reduces memory fourfold by definition. I do not know, and the theory does not say, how much accuracy a real model loses: it depends on the network and on the data, and it is measured on the model. Theorem 26.3 of Polyanskiy and Wu gives another useful fact: with zero mean, the same variance and squared error, R(D) ≤ ½ log+(σ²/D) for every source, so non-Gaussian weights are at least as compressible, and the Gaussian limit is conservative.

Where rate-distortion and capacity meet

Rate-distortion is not a topic separate from capacity: Shannon's theorem 21 says that you can transmit with fidelity v1 over a channel of capacity C if and only if R1 ≤ C, and Polyanskiy and Wu (theorem 26.6) rewrite it as the maximum joint coding rate C / R(D). The case in which the source is an equiprobable bit with Hamming distortion is part 2 of the theorem in article 4: the rate achievable with a bit error pb is C / (1 - h(pb)), that is C / R(pb) with R(pb) = 1 - h(pb) (remark 26.4 of Polyanskiy and Wu). So the limit above capacity, which article 4 calls "accepting bit errors", is a case of lossy compression.

In machine learning the same form appears elsewhere. Alemi et al. (Fixing a Broken ELBO, arXiv 1711.00464) derive variational bounds on the mutual information between the input and the latent variable of a variational autoencoder and obtain from them a rate-distortion curve that describes the trade-off between compression and reconstruction accuracy. They also showed that there are models with the same ELBO and different properties: I read the abstract, not the paper. With the standard definition of the ELBO as reconstruction minus the KL between posterior and prior (the divergence of article 2), the second term acts as the measure of the "rate", and the first as the "distortion".

Common mistakes

• Insisting that distortion be a metric. Shannon requires only a function that measures how much an error costs; symmetry is not needed.
• Swapping the sign in inverse water-filling. The distortion of a component is min(σi², λ), not the maximum.
• Forgetting that the limit is asymptotic. R(D) holds for long blocks; with one bit per weight, the optimal one-bit quantizer stays 1.62 dB from the limit.
• Inferring a model's accuracy from the squared error of the weights alone. The theory quantifies the distortion of the weights; the effect on the output of the network is measured on the model.
• Forgetting that quantization is not invertible. Once quantized you do not recover the original, unlike the lossless compression of article 3.

The common trait of the first two and the fourth is giving a formula a wider field of validity than it has: the first mistake narrows it (metric), the second gets its order wrong (min), the fourth extends it from the error on the weights to accuracy. Before applying it, write down the hypothesis that makes it true.

What you can do now, and what you hold of the series

You can compute R(D) for a Gaussian source with squared error and for a binary source with Hamming distortion, recognise that a fixed-level quantizer wastes from 1.6 to 3.9 dB and a uniform one with entropy coding 0.2546 bits, and allocate bits among components with reverse water-filling. You can also connect lossy compression to capacity. With this article the series is complete. It starts from Shannon entropy, the measure of a source's uncertainty; it goes on to mutual information and Kullback-Leibler divergence, which serve to compare distributions and to measure what an observation reveals; it arrives at source coding with Huffman and Shannon's first theorem, that is lossless compression; it tackles the noisy channel and capacity, which fixes how much can be communicated reliably; it builds the error-correcting codes with Hamming distance and the Hamming and Singleton bounds; and it closes here with lossy compression. The series index gathers all six articles.

Frequently asked questions