Information Theory

Before we can define entropy — the entropy of a distribution — we must define the information content of a single outcome. This is called self-information or surprisal.

The idea is intuitive: learning that a rare event happened gives you more information than learning that a common event happened. If you learn that it rained in Seattle today (high probability), you are barely surprised. If you learn it snowed in Bangalore in June (extremely low probability), you are shocked — and that shock is information.

Formally, for an event with probability p, its self-information is: $$I(x) = \log_2 \frac{1}{p} = -\log_2 p$$

Why the log? There are three axiomatic reasons:

1. Certainty axiom: If p = 1, I = 0. Certain events carry no information.
2. Surprise monotonicity: If p₁ < p₂, then I(x₁) > I(x₂). Rarer events are more informative.
3. Additivity for independence: If two events A and B are independent, I(A ∩ B) = I(A) + I(B). This forces the logarithm. Intuitively: learning two independent facts gives you the sum of their information.

The unit depends on the base:

• Base 2 → bits (the standard in CS)
• Base e → nats (common in statistics and ML)
• Base 10 → hartleys (rare)

Example: A fair coin. p(heads) = 0.5. I(heads) = -log₂(0.5) = 1 bit. This is where the word 'bit' comes from — it is literally one binary digit of surprise, the information content of one fair coin flip.

import math

def self_information(p: float, base: float = 2) -> float:
    """Self-information (surprisal) of an event with probability p."""
    if p <= 0 or p > 1:
        raise ValueError("Probability must be in (0, 1]")
    if p == 1:
        return 0.0
    return -math.log(p, base)

# Examples
print(self_information(0.5))    # 1.0 bit  (fair coin)
print(self_information(0.25))   # 2.0 bits (fair 4-sided die outcome)
print(self_information(1/52))   # ~5.7 bits (drawing a specific card)
print(self_information(1))      # 0.0 bits  (certain event)

Entropy H(X) is the expected self-information of a random variable X. It is the average surprise you experience when you observe an outcome drawn from the distribution.

$$H(X) = -\sum_{x \in \mathcal{X}} p(x) \log_2 p(x)$$

(By convention, 0 · log(0) = 0 since lim_{p→0} p log p = 0.)

Intuition via extremes:

• A deterministic variable (one outcome has probability 1) has H = 0 bits. You always know what will happen, so observing it tells you nothing.
• A uniform distribution over n outcomes has H = log₂(n) bits — this is the maximum possible entropy for n outcomes. Uncertainty is maximized when all outcomes are equally likely.

Properties:

1. H(X) ≥ 0 always.
2. H(X) ≤ log₂|𝒳|, with equality iff X is uniform. (Proved by Jensen's inequality applied to the concave log function.)
3. H(X) is a concave function of the probability vector.
4. H is continuous in the probabilities.

The binary entropy function H_b(p) = -p log₂ p - (1-p) log₂(1-p) is special — it is the entropy of a Bernoulli(p) variable. It equals 0 at p=0 and p=1 (no uncertainty), and peaks at 1 bit when p=0.5 (maximum uncertainty).

Analogy: Entropy is like the average number of yes/no questions you need to ask to identify the outcome. If one outcome is overwhelmingly likely, you need very few questions. If all outcomes are equally likely, you need log₂(n) questions. This is not just an analogy — the source coding theorem makes it precise.

import math
from collections import Counter

def entropy(probs: list[float], base: float = 2) -> float:
    """Shannon entropy of a discrete distribution."""
    return -sum(p * math.log(p, base) for p in probs if p > 0)

# From a sequence of symbols
def empirical_entropy(data) -> float:
    n = len(data)
    counts = Counter(data)
    probs = [c / n for c in counts.values()]
    return entropy(probs)

# A loaded die heavily favoring 1
probs_loaded = [0.9, 0.02, 0.02, 0.02, 0.02, 0.02]
print(f"Loaded die entropy:  {entropy(probs_loaded):.4f} bits")  # ~0.57

# A fair die
probs_fair = [1/6] * 6
print(f"Fair die entropy:    {entropy(probs_fair):.4f} bits")  # ~2.585
print(f"log2(6) =            {math.log2(6):.4f} bits")          # same

# English text has ~1 bit/char of true entropy (Markov model),
# but 0th-order entropy is ~4.2 bits
print(empirical_entropy(list("hello world")))

For two random variables X and Y:

Joint entropy is entropy of the pair: $$H(X,Y) = -\sum_{x,y} p(x,y) \log_2 p(x,y)$$

Conditional entropy H(Y|X) is the average entropy of Y given knowledge of X: $$H(Y|X) = \sum_x p(x) H(Y|X=x) = -\sum_{x,y} p(x,y) \log_2 p(y|x)$$

Conditional entropy measures the remaining uncertainty in Y after you have observed X. It is not the same as H(Y|X=x) for a specific x — it is the expectation over x.

The chain rule for entropy: $$H(X,Y) = H(X) + H(Y|X)$$

This is intuitive: the information in (X,Y) together equals the information in X alone, plus the additional information in Y once you know X. It generalizes to any number of variables: $$H(X_1, X_2, \ldots, X_n) = \sum_{i=1}^n H(X_i | X_1, \ldots, X_{i-1})$$

Key inequality: Conditioning never increases entropy: $$H(Y|X) \leq H(Y)$$ with equality iff X and Y are independent. Knowing something about the world can only reduce (or maintain) uncertainty — it never adds uncertainty. This is one of the most important inequalities in all of information theory.

Subadditivity: H(X,Y) ≤ H(X) + H(Y), with equality iff X and Y are independent. Joint information is at most the sum of individual information.

import numpy as np
from itertools import product

def joint_entropy(joint_pmf: np.ndarray) -> float:
    """Entropy of the joint distribution given as a 2D array."""
    flat = joint_pmf.flatten()
    return -sum(p * np.log2(p) for p in flat if p > 0)

def conditional_entropy(joint_pmf: np.ndarray) -> float:
    """H(Y|X) where axis 0 is X, axis 1 is Y."""
    px = joint_pmf.sum(axis=1)
    result = 0.0
    for i, pxi in enumerate(px):
        if pxi > 0:
            py_given_x = joint_pmf[i] / pxi
            result += pxi * (-sum(p * np.log2(p) for p in py_given_x if p > 0))
    return result

# Example: X = weather (sunny/rainy), Y = umbrella (yes/no)
# P(X=sunny, Y=no)=0.5, P(X=sunny, Y=yes)=0.1,
# P(X=rainy, Y=no)=0.1, P(X=rainy, Y=yes)=0.3
joint = np.array([[0.5, 0.1],
                  [0.1, 0.3]])

hx = -sum(p * np.log2(p) for p in joint.sum(axis=1) if p > 0)
hy = -sum(p * np.log2(p) for p in joint.sum(axis=0) if p > 0)
hxy = joint_entropy(joint)
hy_given_x = conditional_entropy(joint)

print(f"H(X)      = {hx:.4f}")
print(f"H(Y)      = {hy:.4f}")
print(f"H(X,Y)    = {hxy:.4f}")
print(f"H(Y|X)    = {hy_given_x:.4f}")
print(f"Chain rule check H(X)+H(Y|X) = {hx + hy_given_x:.4f} == H(X,Y) = {hxy:.4f}")

Mutual information I(X;Y) measures how much knowing X reduces uncertainty about Y (and vice versa — it is symmetric):

$$I(X;Y) = H(X) + H(Y) - H(X,Y) = H(X) - H(X|Y) = H(Y) - H(Y|X)$$

Equivalently: $$I(X;Y) = \sum_{x,y} p(x,y) \log_2 \frac{p(x,y)}{p(x)p(y)}$$

This last form is the KL divergence between the joint distribution and the product of marginals — it measures how far the joint is from independence.

Key properties:

• I(X;Y) ≥ 0, with equality iff X and Y are independent.
• I(X;Y) = I(Y;X) — mutual information is symmetric.
• I(X;X) = H(X) — a variable has maximum mutual information with itself.
• I(X;Y) ≤ min(H(X), H(Y)) — you can't extract more information about Y from X than is contained in either variable alone.

The information diagram: Draw H(X) and H(Y) as overlapping circles. The overlap is I(X;Y). The non-overlapping part of H(X) is H(X|Y), and vice versa. The total area is H(X,Y). This Venn-diagram intuition is exact mathematically.

Applications:

• Feature selection: High I(feature; label) means the feature is informative. This is the foundation of information-gain-based decision trees (CART, ID3, C4.5).
• Channel capacity: The maximum mutual information over input distributions is the channel capacity.
• Neuroscience: Used to measure how much a neuron's firing pattern encodes about the stimulus.
• Compression: I(X;Y) = H(X) - H(X|Y) tells you how much a correlated variable Y compresses the transmission of X.

import numpy as np

def mutual_information(joint_pmf: np.ndarray) -> float:
    """I(X;Y) from joint distribution (axis 0 = X, axis 1 = Y)."""
    px = joint_pmf.sum(axis=1, keepdims=True)
    py = joint_pmf.sum(axis=0, keepdims=True)
    result = 0.0
    for i in range(joint_pmf.shape[0]):
        for j in range(joint_pmf.shape[1]):
            pxy = joint_pmf[i, j]
            if pxy > 0:
                result += pxy * np.log2(pxy / (px[i, 0] * py[0, j]))
    return result

# Perfect correlation: knowing X tells you exactly Y
joint_perfect = np.array([[0.5, 0.0],
                           [0.0, 0.5]])
print(f"Perfect corr I(X;Y) = {mutual_information(joint_perfect):.4f}")  # 1.0

# Independence: X and Y are uncorrelated
joint_indep = np.array([[0.25, 0.25],
                         [0.25, 0.25]])
print(f"Independence I(X;Y) = {mutual_information(joint_indep):.4f}")    # 0.0

# Partial correlation
joint_partial = np.array([[0.4, 0.1],
                           [0.1, 0.4]])
print(f"Partial corr I(X;Y) = {mutual_information(joint_partial):.4f}")  # ~0.278

KL divergence (Kullback-Leibler divergence), also called relative entropy, measures how different one probability distribution Q is from a reference distribution P:

$$D_{KL}(P | Q) = \sum_x p(x) \log_2 \frac{p(x)}{q(x)}$$

It answers: if nature generates data from P but you are using a code (model) optimized for Q, how many extra bits per symbol do you waste?

Critical properties:

• D_KL(P‖Q) ≥ 0 always (Gibbs' inequality / information inequality). Proved via Jensen's inequality: -log is convex, so E[-log(Q/P)] ≥ -log(E[Q/P]) = 0.
• D_KL(P‖Q) = 0 iff P = Q everywhere.
• D_KL is NOT symmetric: D_KL(P‖Q) ≠ D_KL(Q‖P) in general.
• D_KL is NOT a metric (no triangle inequality), despite being a measure of "distance."

Cross-entropy H(P, Q) is defined as: $$H(P, Q) = -\sum_x p(x) \log_2 q(x) = H(P) + D_{KL}(P | Q)$$

Cross-entropy is the expected number of bits used to encode events from P using a code optimized for Q. When Q = P, this reduces to H(P) — perfect compression. The excess H(P,Q) - H(P) = D_KL(P‖Q) is the wasted bits due to model mismatch.

Why cross-entropy is the canonical ML loss: In supervised learning, P is the true label distribution (a one-hot vector for classification), and Q is the model's predicted distribution (softmax output). Minimizing cross-entropy H(P,Q) is equivalent to minimizing KL divergence D_KL(P‖Q) (since H(P) is constant w.r.t. model parameters), which is equivalent to maximum likelihood estimation. This is not a heuristic choice — it is the information-theoretically correct loss.

Forward vs reverse KL:

• D_KL(P‖Q): 'forward KL'. When Q ≈ 0 but P > 0, the term blows up to infinity — Q is forced to cover all regions where P has mass (mode-covering behavior). Used in ML training (MLE).
• D_KL(Q‖P): 'reverse KL'. When P ≈ 0 but Q > 0, the term blows up — Q is forced to avoid regions where P has no mass (mode-seeking behavior). Used in variational inference.

import numpy as np

def kl_divergence(p: np.ndarray, q: np.ndarray) -> float:
    """D_KL(P || Q) in bits. Assumes same support."""
    # Clip to avoid log(0); p*log(p/q) -> 0 when p -> 0
    mask = p > 0
    return float(np.sum(p[mask] * np.log2(p[mask] / q[mask])))

def cross_entropy(p: np.ndarray, q: np.ndarray) -> float:
    """H(P, Q) in bits."""
    mask = p > 0
    return float(-np.sum(p[mask] * np.log2(q[mask])))

# True distribution: mostly class 0
p = np.array([0.9, 0.05, 0.05])

# Model 1: almost correct
q1 = np.array([0.85, 0.1, 0.05])
# Model 2: totally wrong
q2 = np.array([0.1, 0.1, 0.8])

print(f"H(P)         = {cross_entropy(p, p):.4f} bits")  # entropy of P
print(f"H(P, Q1)     = {cross_entropy(p, q1):.4f} bits")  # low waste
print(f"H(P, Q2)     = {cross_entropy(p, q2):.4f} bits")  # high waste
print(f"KL(P||Q1)    = {kl_divergence(p, q1):.4f} bits")  # small mismatch
print(f"KL(P||Q2)    = {kl_divergence(p, q2):.4f} bits")  # large mismatch
print(f"KL(Q2||P)    = {kl_divergence(q2, p):.4f} bits")  # asymmetric!

# In PyTorch, cross_entropy loss is exactly this:
# loss = -sum(p * log_softmax(logits)) = H(one_hot, softmax(logits))

Shannon's source coding theorem is one of the most profound results in all of science. It answers: what is the fundamental limit of lossless data compression?

Setup: You have a discrete memoryless source (DMS) that produces symbols i.i.d. from distribution P over alphabet 𝒳. You want to encode long sequences of n symbols into binary strings, then perfectly decode them.

Shannon's Source Coding Theorem: For any ε > 0 and any rate R > H(X) bits/symbol, there exists an encoder/decoder pair such that for large enough n, the code works with probability ≥ 1-ε using at most R bits/symbol. Conversely, if R < H(X), any code must fail on a non-negligible fraction of inputs.

In words: H(X) bits/symbol is the sharp threshold for lossless compression. You cannot do better on average; you can always approach it.

The proof idea — typical sequences: The key is the asymptotic equipartition property (AEP): By the law of large numbers, for large n: $$-\frac{1}{n} \log p(x_1, \ldots, x_n) \to H(X) \quad \text{in probability}$$

This means that with high probability, the sequence (x₁,...,xₙ) is in the typical set 𝒜_ε^(n), whose members all have probability approximately 2^(-nH(X)). The typical set has:

• At most 2^(n(H+ε)) members.
• Probability ≥ 1 - ε for large n.

The scheme: Assign short codes (nR bits) to the 2^(nR) ≈ 2^(nH) typical sequences. Everything else (atypical, probability < ε) gets a special 'error' flag. This achieves rate H(X) with vanishing error probability.

What it means for practical compression: ZIP, gzip, zstd, bzip2 are all trying to approach this limit. They use learned models of the source statistics (Lempel-Ziv algorithms approach the entropy of the stationary process). The theorem tells us when we have 'solved' compression — and the answer is not 'compress until unreadable,' but 'compress until we reach H(X).'

# Illustration of the AEP: empirical entropy rate converges to H(X)
import numpy as np
import math

def demo_aep(probs: list, n_trials: int = 100_000, n: int = 1000, seed=42):
    """
    Draw n-length sequences from the DMS with given probs.
    Check that -1/n * log2 P(sequence) ≈ H(X) for most sequences.
    """
    rng = np.random.default_rng(seed)
    k = len(probs)
    H = -sum(p * math.log2(p) for p in probs if p > 0)

    # Sample n_trials sequences of length n
    log_probs_per_seq = []
    log_p = np.log2(probs)
    for _ in range(n_trials):
        seq = rng.choice(k, size=n, p=probs)
        log_prob = np.sum(log_p[seq])   # log P(x_1,...,x_n)
        log_probs_per_seq.append(-log_prob / n) # should ≈ H

    empirical = np.mean(log_probs_per_seq)
    std_dev   = np.std(log_probs_per_seq)
    print(f"True H(X)          = {H:.4f} bits")
    print(f"Mean -1/n*logP(seq) = {empirical:.4f} bits  (should ≈ H)")
    print(f"Std dev            = {std_dev:.5f} bits  (→0 as n→∞)")
    # Fraction in ε-typical set (within ε=0.05 of H)
    eps = 0.05
    typical = sum(abs(v - H) < eps for v in log_probs_per_seq)
    print(f"Fraction typical (ε={eps}): {typical/n_trials:.3f}  (→1 as n→∞)")

demo_aep(probs=[0.5, 0.25, 0.125, 0.125])

The source coding theorem guarantees H(X) is the limit, but does not give us a practical code. Huffman coding is an algorithm that constructs an optimal prefix-free (instantaneous) binary code for a known distribution P.

Prefix-free codes: A code is prefix-free if no codeword is a prefix of another. This allows unambiguous decoding without separators. Prefix-free codes correspond bijectively to binary trees (leaves = symbols, path from root = codeword). The expected code length is: $$L = \sum_i p_i \ell_i$$ where ℓᵢ is the length of codeword for symbol i.

Kraft's inequality: A prefix-free code over a binary alphabet exists iff $\sum_i 2^{-\ell_i} \leq 1$. Huffman codes achieve equality.

Huffman's algorithm (1952, David Huffman as a student assignment):

1. Create a leaf node for each symbol with its probability.
2. Build a min-heap (priority queue) by probability.
3. Repeatedly: extract the two lowest-probability nodes, create a new internal node with their sum as probability, add it back to the heap.
4. When one node remains, it is the root. The tree is complete.

Optimality theorem: Huffman coding achieves the minimum expected length among all prefix-free codes. The expected length L satisfies: $$H(X) \leq L < H(X) + 1$$

So Huffman wastes at most 1 bit per symbol. For block coding (encoding k symbols at a time), the waste drops to 1/k bits/symbol — approaching H(X) from above.

Limitation: Huffman coding requires knowing the distribution P exactly in advance. Adaptive and universal codes (Lempel-Ziv, arithmetic coding) relax this requirement. Arithmetic coding can approach H(X) more closely than Huffman by using fractional bits per symbol.

import heapq
from dataclasses import dataclass, field
from typing import Optional

@dataclass(order=True)
class HNode:
    prob: float
    symbol: Optional[str] = field(default=None, compare=False)
    left:  Optional['HNode'] = field(default=None, compare=False)
    right: Optional['HNode'] = field(default=None, compare=False)

def huffman_tree(probs: dict[str, float]) -> HNode:
    heap = [HNode(p, s) for s, p in probs.items()]
    heapq.heapify(heap)
    while len(heap) > 1:
        a = heapq.heappop(heap)
        b = heapq.heappop(heap)
        parent = HNode(a.prob + b.prob, left=a, right=b)
        heapq.heappush(heap, parent)
    return heap[0]

def build_codebook(node: HNode, prefix='') -> dict[str, str]:
    if node.symbol is not None:   # leaf
        return {node.symbol: prefix or '0'}  # single symbol edge case
    codes = {}
    if node.left:  codes.update(build_codebook(node.left,  prefix + '0'))
    if node.right: codes.update(build_codebook(node.right, prefix + '1'))
    return codes

# Classic example from Shannon's paper
probs = {'a': 0.5, 'b': 0.25, 'c': 0.125, 'd': 0.125}
root = huffman_tree(probs)
codes = build_codebook(root)

print("Symbol | Prob  | Code | Length")
for sym, p in sorted(probs.items()):
    print(f"  {sym}    | {p:.3f} | {codes[sym]:<4} | {len(codes[sym])}")

L = sum(probs[s] * len(codes[s]) for s in probs)
import math
H = -sum(p * math.log2(p) for p in probs.values())
print(f"\nH(X) = {H:.4f} bits,  E[L] = {L:.4f} bits,  waste = {L-H:.4f} bits")
# Expected: H=1.75, L=1.75 (this distribution is 'dyadic' — Huffman is exact!)

Shannon's channel capacity theorem answers a different question: given a noisy communication channel, what is the maximum reliable information rate?

The model: A discrete memoryless channel (DMC) is defined by a conditional distribution P(Y|X) — given input symbol X, the output Y is randomly determined by the channel. For example, a binary symmetric channel (BSC) with error probability p flips each bit independently with probability p.

Channel capacity C is defined as: $$C = \max_{P(X)} I(X; Y)$$

This is the maximum mutual information between input and output, maximized over all input distributions P(X). It is measured in bits per channel use.

Shannon's Noisy Channel Coding Theorem (1948): For any rate R < C, there exists a sequence of codes such that the probability of error → 0 as block length n → ∞. For any R > C, the probability of error → 1.

C is the sharp threshold. Below it, perfect reliable communication is possible. Above it, it is not.

Examples:

• BSC(p): C = 1 - H_b(p). At p=0.5 (completely random), C=0. At p=0, C=1 (perfect channel).
• Binary erasure channel (BEC) with erasure prob ε: C = 1 - ε.
• AWGN channel with bandwidth B and SNR: C = B · log₂(1 + SNR) — the Shannon-Hartley theorem. This is why doubling power gives only log gain in capacity, while doubling bandwidth gives linear gain.

The proof sketch: As before, random codes over the typical set work. With high probability, a random codebook of size 2^(nR) for R < C will have the property that different codewords produce distinguishable output sequences. Jointly typical decoding achieves vanishing error.

Practical impact: This theorem told engineers in 1948 that reliable digital communication was possible at all! Before Shannon, engineers believed you had to trade off reliability against speed. Shannon showed there is a threshold C — below it you get both. Modern codes (turbo codes, LDPC, polar codes) essentially achieve C.

import math
import numpy as np
from scipy.optimize import minimize_scalar

def binary_entropy(p: float) -> float:
    if p == 0 or p == 1: return 0.0
    return -p * math.log2(p) - (1-p) * math.log2(1-p)

def bsc_capacity(p: float) -> float:
    """Capacity of Binary Symmetric Channel with crossover prob p."""
    return 1.0 - binary_entropy(p)

def bec_capacity(epsilon: float) -> float:
    """Capacity of Binary Erasure Channel with erasure prob epsilon."""
    return 1.0 - epsilon

def awgn_capacity_bits(snr_linear: float, bandwidth_hz: float = 1.0) -> float:
    """Shannon-Hartley: C = B * log2(1 + SNR), in bits/sec."""
    return bandwidth_hz * math.log2(1 + snr_linear)

# BSC: capacity vs error rate
print("BSC capacity:")
for p in [0.0, 0.05, 0.1, 0.2, 0.5]:
    print(f"  p={p:.2f}: C = {bsc_capacity(p):.4f} bits/use")

# AWGN: bandwidth vs power tradeoff
print("\nAWGN (SNR=10dB=10x linear):")
for B in [1, 2, 4, 8]:
    C = awgn_capacity_bits(10.0, B)
    print(f"  B={B} Hz: C = {C:.4f} bits/sec")

print("\nAWGN (B=1 Hz, varying SNR):")
for snr_db in [0, 3, 10, 20, 30]:
    snr_lin = 10 ** (snr_db / 10)
    C = awgn_capacity_bits(snr_lin)
    print(f"  SNR={snr_db:2d} dB: C = {C:.4f} bits/use")