import matplotlib
if not hasattr(matplotlib.RcParams, "_get"):
    matplotlib.RcParams._get = dict.get

15.3 Lossless compression

In the previous section, our key question was how to remove perceptually redundant information. Here, the key question shifts: how can we leverage statistical redundancy to compress with absolutely no distortion?

Lossless codecs carry a robust guarantee of zero distortion:

\[\text{Decoder}(\text{Coder}(\mathbf{x})) = \mathbf{x}, \qquad \texttt{error}\big(\mathbf{x}, \text{Decoder}(\text{Coder}(\mathbf{x}))\big) = 0.\]

Our goal is to minimize the information rate, subject to this guarantee. We would also like the coder and decoder to be fast and simple to compute. For digital audio, lossless codecs like FLAC typically achieve compression factors of about \(1.4\times\) to \(2\times\), far less than lossy codecs, but with a perfect round trip.

The key idea of lossless compression is that the more predictable a signal is, the more we can compress it. We saw this idea in action with run-length encoding, which compressed strings with predictably long runs of \(0\)s. Let’s make this idea more precise.

Entropy

Imagine creating strings of bits by tossing a weighted coin, which comes up heads with probability \(p\). We write down \(1\) for heads and \(0\) for tails.

• If \(p = 0.01\), we get tails \(99\%\) of the time, so our strings have very long runs of \(0\)s. We already know of strategies to compress these strings, such as run-length encoding.

• If \(p = 0.5\) (a fair coin), the strings look random, and it is not obvious how we could do better than using \(1\) bit per toss.

What is the theoretical minimum number of bits per toss that we could use, as a function of \(p\)? Claude Shannon answered this question in 1948 [Sha48]. (You might remember him from the Nyquist-Shannon sampling theorem!) To answer it, he introduced a measure of the information in a probability distribution, which he called entropy (after a mathematically similar quantity in physics). His work founded the field of information theory:

Definition 52 (Entropy)

The entropy of a random variable \(X\) with alphabet \(\mathcal{X}\) is

\[H(X) \coloneqq -\sum_{x \in \mathcal{X}} P(x) \log P(x).\]

With a base-\(2\) logarithm, entropy is measured in bits.

Shannon’s source coding theorem states that no lossless codec can use fewer than \(H(X)\) bits per outcome, on average. This implies two goals: (1) estimate \(H(X)\) if possible, and (2) design a codec that can come arbitrarily close to this bound. Entropy measures how unpredictable, or surprising, a random variable is. The more predictable the outcomes, the lower the entropy, and the more we can compress them.

Let’s compute the entropy of our coin tosses, a random variable \(X\) with alphabet \(\{\text{heads}, \text{tails}\}\). For a fair coin, \(p = 0.5\):

\[H(X) = -\big(0.5 \log_2 0.5 + 0.5 \log_2 0.5\big) = -\log_2 2^{-1} = 1 \text{ bit per toss}.\]

As we suspected, we can’t do better than \(1\) bit per toss. For a coin with \(p = 0.01\):

\[H(X) = -\big(0.01 \log_2 0.01 + 0.99 \log_2 0.99\big) \approx 0.066 + 0.014 \approx 0.081 \text{ bits per toss}.\]

This implies that we could compress these strings with a compression factor of up to \(\frac{1}{0.081} \approx 12.4\times\)!

Fig. 126 The entropy of a coin toss, as a function of the probability of heads.

Huffman coding

How can we design a practical compression algorithm that approaches the theoretical bound imposed by entropy? One classic algorithm is Huffman coding, which David Huffman invented in 1952 [Huf52]. The idea is to assign each outcome a codeword (a string of bits), giving short codewords to likely outcomes and long codewords to unlikely ones.

With only two outcomes, heads and tails, the shortest possible codewords have \(1\) bit each, so we can’t do better than \(1\) bit per toss. Instead, we encode blocks of \(N\) tosses at a time. With \(N = 2\), there are four outcomes, HH, HT, TH, and TT. Let’s use a coin with \(p = 0.1\), whose entropy is about \(0.469\) bits per toss. Since the tosses are independent, the probability of each block is the product of the probabilities of its tosses, for example, \(P(\text{TT}) = 0.9 \cdot 0.9 = 0.81\).

Algorithm 1 (Huffman coding)

Given an alphabet and the probability of each outcome:

1. Create a node for each outcome, labeled with its probability.

2. While more than one node remains, merge the two least probable nodes into a new parent node, whose probability is the sum of theirs. Label the edges to its two children \(0\) and \(1\).

3. The codeword for each outcome is the sequence of edge labels on the path from the final (root) node to that outcome.

For our four outcomes, we first merge HH (\(0.01\)) with HT (\(0.09\)), then their parent (\(0.10\)) with TH (\(0.09\)), then their parent (\(0.19\)) with TT (\(0.81\)):

Fig. 127 The Huffman tree for blocks of \(N = 2\) tosses of a coin with \(P(\text{heads}) = 0.1\). Each codeword is the sequence of edge labels from the root to a leaf.

This gives the codewords TT → \(\texttt{0}\), TH → \(\texttt{11}\), HT → \(\texttt{100}\), and HH → \(\texttt{101}\). Notice that no codeword is the beginning (prefix) of another. This makes Huffman codes easy to decode, even though the codewords have different lengths. For example, \(\texttt{0110100}\) can only be read as \(\texttt{0}\,|\,\texttt{11}\,|\,\texttt{0}\,|\,\texttt{100}\), or TT TH TT HT.

The common case, a string of all tails, takes just \(1\) bit per pair, or \(\frac{1}{2}\) bit per toss. A rare string can take up to \(3\) bits per pair. In expectation, the code takes

\[0.81 \cdot 1 + 0.09 \cdot 2 + 0.09 \cdot 3 + 0.01 \cdot 3 = 1.29 \text{ bits per pair} = 0.645 \text{ bits per toss},\]

a compression factor of \(\frac{1}{0.645} \approx 1.55\times\), approaching the theoretical limit of \(\frac{1}{0.469} \approx 2.13\times\). As the block length \(N\) grows, the expected length of a Huffman code approaches the entropy. The widget below builds the Huffman tree for any block length \(N\) and probability \(p\):

# hide
# no-output
from IPython.utils.capture import capture_output
with capture_output():
    %pip install -q plotly anywidget

import heapq
import itertools
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import ipywidgets as widgets
import icm_plotly
from icm_plotly import RED, BLUE, GOLD, IRON, STEEL

# hide
# autorun
N_MAX = 6                                     # longest block of tosses
P_GRID = [round(0.01 * k, 2) for k in range(1, 51)]   # the P(heads) slider's stops
N_START, P_START = 2, 0.1
LABEL_MAX = 8                                 # label the tree only up to this many leaves


def block_probs(n: int, p: float) -> dict[str, float]:
    """Every block of n tosses ("HTT", ...) and its probability."""
    return {"".join(b): p ** b.count("H") * (1 - p) ** b.count("T")
            for b in itertools.product("HT", repeat=n)}


def huffman(probs: dict[str, float]) -> str | tuple:
    """The Huffman tree: a leaf is an outcome, a node is a pair (bit-0 child, bit-1 child)."""
    # the counter breaks ties between equal probabilities, so the tree is deterministic
    heap = [(prob, k, outcome) for k, (outcome, prob) in enumerate(probs.items())]
    heapq.heapify(heap)
    k = len(heap)
    while len(heap) > 1:
        p_a, _, a = heapq.heappop(heap)       # the least probable node
        p_b, _, b = heapq.heappop(heap)       # the next least probable
        heapq.heappush(heap, (p_a + p_b, k, (b, a)))   # the more probable child gets bit 0
        k += 1
    return heap[0][2]


def codewords(tree: str | tuple, prefix: str = "") -> dict[str, str]:
    """Read each outcome's codeword off the tree: the bits on the path from the root."""
    if isinstance(tree, str):
        return {tree: prefix}
    return {**codewords(tree[0], prefix + "0"), **codewords(tree[1], prefix + "1")}


def bits_per_toss(n: int, p: float) -> float:
    """Expected codeword length per block, divided by the n tosses in a block."""
    probs = block_probs(n, p)
    codes = codewords(huffman(probs))
    return sum(probs[b] * len(codes[b]) for b in probs) / n


def entropy(p: float) -> float:
    return -p * np.log2(p) - (1 - p) * np.log2(1 - p) if 0 < p < 1 else 0.0


def tree_layout(tree: str | tuple) -> dict[str, list]:
    """Plot coordinates: leaves left to right in order, each node above its children, y = depth."""
    lay = {"edge_x": [], "edge_y": [], "node_x": [], "node_y": [],
           "leaf_x": [], "leaf_y": [], "leaves": [], "bit_x": [], "bit_y": [], "bits": []}
    count = itertools.count()

    def place(node, depth):
        if isinstance(node, str):
            x = next(count)
            lay["leaf_x"].append(x)
            lay["leaf_y"].append(depth)
            lay["leaves"].append(node)
            return x
        xs = [place(child, depth + 1) for child in node]
        x = (xs[0] + xs[1]) / 2
        lay["node_x"].append(x)
        lay["node_y"].append(depth)
        for bit, x_child in enumerate(xs):
            lay["edge_x"] += [x, x_child, None]
            lay["edge_y"] += [depth, depth + 1, None]
            lay["bit_x"].append((x + x_child) / 2)
            lay["bit_y"].append(depth + 0.5)
            lay["bits"].append(str(bit))
        return x

    place(tree, 0)
    return lay


# the right panel's curves and each N's deepest tree (so the axes only move when N does)
BITS = {p: [bits_per_toss(n, p) for n in range(1, N_MAX + 1)] for p in P_GRID}
DEPTH = {n: max(len(c) for p in P_GRID for c in codewords(huffman(block_probs(n, p))).values())
         for n in range(1, N_MAX + 1)}


def tree_traces(n: int, p: float) -> dict:
    """Everything the left panel draws for one (N, p)."""
    probs = block_probs(n, p)
    tree = huffman(probs)
    codes = codewords(tree)
    lay = tree_layout(tree)
    labeled = len(probs) <= LABEL_MAX
    lay["labels"] = [f"{b}<br>{codes[b]}" if labeled else "" for b in lay["leaves"]]
    lay["hover"] = [f"{b}: P = {probs[b]:.3g}, codeword {codes[b]}" for b in lay["leaves"]]
    if not labeled:
        lay["bit_x"], lay["bit_y"], lay["bits"] = [], [], []
    lay["bit_pos"] = ["top left" if b == "0" else "top right" for b in lay["bits"]]
    lay["leaf_size"] = 8 if len(probs) <= 16 else 5
    lay["x_range"] = [-0.7, len(probs) - 0.3]
    lay["y_range"] = [DEPTH[n] + (1.1 if labeled else 0.5), -0.5]
    lay["y_dtick"] = 1 if DEPTH[n] <= 6 else None   # None: plotly picks the ticks
    return lay


def readout(n: int, p: float) -> str:
    per_toss = BITS[p][n - 1]
    h = entropy(p)
    return (f"N = {n}: {per_toss * n:.2f} bits per block → {per_toss:.3f} bits per toss "
            f"({1 / per_toss:.2f}× compression). "
            f"Entropy (theoretical limit): {h:.3f} bits per toss ({1 / h:.2f}×).")


def figure():
    fig = make_subplots(rows=1, cols=2, column_widths=[0.62, 0.38], horizontal_spacing=0.1)
    lay = tree_traces(N_START, P_START)
    # left: the tree (0 edges, 1 internal nodes, 2 leaves, 3 edge bits)
    fig.add_scatter(x=lay["edge_x"], y=lay["edge_y"], mode="lines",
                    line=dict(color=IRON, width=1.2), hoverinfo="skip", row=1, col=1)
    fig.add_scatter(x=lay["node_x"], y=lay["node_y"], mode="markers",
                    marker=dict(color=IRON, size=5), hoverinfo="skip", row=1, col=1)
    fig.add_scatter(x=lay["leaf_x"], y=lay["leaf_y"], mode="markers+text",
                    marker=dict(color=RED, size=lay["leaf_size"]), text=lay["labels"],
                    textposition="bottom center", textfont=dict(size=11),
                    hovertext=lay["hover"], hoverinfo="text", row=1, col=1)
    fig.add_scatter(x=lay["bit_x"], y=lay["bit_y"], mode="text", text=lay["bits"],
                    textposition=lay["bit_pos"], textfont=dict(size=11, color=IRON),
                    hoverinfo="skip", row=1, col=1)
    fig.update_xaxes(range=lay["x_range"], showline=False, ticks="", showticklabels=False,
                     fixedrange=True, row=1, col=1)
    fig.update_yaxes(range=lay["y_range"], dtick=lay["y_dtick"],
                     title_text="Codeword length (bits)", fixedrange=True, row=1, col=1)

    # right: bits per toss against block length (4 entropy, 5 Huffman, 6 the current N)
    ns = list(range(1, N_MAX + 1))
    h = entropy(P_START)
    fig.add_scatter(x=[0.5, 0.65, N_MAX + 0.5], y=[h, h, h], mode="lines+text",
                    line=dict(color=RED, width=1.6, dash="dash"), text=["", "entropy (theoretical limit)", ""],
                    textposition="bottom right", textfont=dict(color=RED),
                    hovertemplate="entropy (theoretical limit): %{y:.3f} bits per toss<extra></extra>",
                    row=1, col=2)
    fig.add_scatter(x=ns, y=BITS[P_START], mode="lines+markers",
                    line=dict(color=BLUE, width=2), marker=dict(color=BLUE, size=7),
                    hovertemplate="N = %{x}: %{y:.3f} bits per toss<extra>Huffman</extra>",
                    row=1, col=2)
    fig.add_scatter(x=[N_START], y=[BITS[P_START][N_START - 1]], mode="markers",
                    marker=dict(color=GOLD, size=13, line=dict(color=IRON, width=1)),
                    hoverinfo="skip", row=1, col=2)
    fig.update_xaxes(range=[0.5, N_MAX + 0.5], dtick=1, title_text="Block length N",
                     fixedrange=True, row=1, col=2)
    fig.update_yaxes(range=[0, 1.05], title_text="Bits per toss", fixedrange=True,
                     row=1, col=2)
    fig.update_layout(height=420)
    return fig


def controls(fig):
    n_slider = widgets.IntSlider(description="Block length $N$", min=1, max=N_MAX,
                                 value=N_START)
    p_slider = widgets.FloatSlider(description="P(heads) $p$", min=0.01, max=0.5, step=0.01,
                                   value=P_START, readout_format=".2f")
    text = widgets.HTML(f"<span style='font-size:0.9em'>{readout(N_START, P_START)}</span>")

    # the defaults snapshot the helpers; the page's notebooks share one kernel
    def update(n, p, tree_traces=tree_traces, readout=readout, entropy=entropy, BITS=BITS):
        p = round(p, 2)
        lay = tree_traces(n, p)
        h = entropy(p)
        with fig.batch_update():
            fig.data[0].x, fig.data[0].y = lay["edge_x"], lay["edge_y"]
            fig.data[1].x, fig.data[1].y = lay["node_x"], lay["node_y"]
            fig.data[2].x, fig.data[2].y = lay["leaf_x"], lay["leaf_y"]
            fig.data[2].text, fig.data[2].hovertext = lay["labels"], lay["hover"]
            fig.data[2].marker.size = lay["leaf_size"]
            fig.data[3].x, fig.data[3].y = lay["bit_x"], lay["bit_y"]
            fig.data[3].text, fig.data[3].textposition = lay["bits"], lay["bit_pos"]
            fig.layout.xaxis.range = lay["x_range"]
            fig.layout.yaxis.range = lay["y_range"]
            fig.layout.yaxis.dtick = lay["y_dtick"]
            fig.data[4].y = [h, h, h]
            fig.data[5].y = BITS[p]
            fig.data[6].x, fig.data[6].y = [n], [BITS[p][n - 1]]
        text.value = f"<span style='font-size:0.9em'>{readout(n, p)}</span>"

    widgets.interactive_output(update, {"n": n_slider, "p": p_slider})
    return widgets.VBox([n_slider, p_slider, text])


icm_plotly.show(figure, controls)

Block length \(N\)2

P(heads) \(p\)0.10

N = 2: 1.29 bits per block → 0.645 bits per toss (1.55× compression). Entropy (theoretical limit): 0.469 bits per toss (2.13×).

Huffman coding is how MP3 losslessly compresses its quantized spectrogram. The mostly-zero coefficients are highly predictable, so Huffman codes compress them efficiently. Ogg Vorbis uses closely related codebooks.

FLAC

The techniques we’ve studied so far work well for coin tosses, which are independent of one another. But for audio, we know that each sample is closely related to its neighbors. How can we exploit this?

The Free Lossless Audio Codec (FLAC), developed by Josh Coalson in 2001 and now maintained by the Xiph.Org Foundation, takes a two-step approach. It (1) uses a statistical model to predict each sample from the samples before it, removing the redundancy between neighboring samples, and (2) losslessly compresses whatever the model could not predict. FLAC splits audio into blocks of a few thousand samples, and for each block, the coder:

1. Fits a linear predictive coding (LPC) model of order \(K\), from Chapter 14, which predicts each sample as a linear combination of the \(K\) samples before it, with coefficients \(h_1, \ldots, h_K\).

2. Quantizes the coefficients to integers.

3. Computes the residual, the difference between the signal and the prediction.

4. Losslessly compresses the residual, using Rice coding.

The decoder reverses the process. Each step deserves a closer look.

Quantizing the coefficients. Our samples are integers (for example, \(16\)-bit PCM), but the optimal LPC coefficients are real numbers. To guarantee a perfect round trip, the coder and decoder must compute exactly the same predictions, so FLAC sticks to integer arithmetic. It scales each coefficient by \(2^s\) (for some shift \(s\), such as \(12\)) and rounds it to an integer, \(q_k = \text{round}(2^s h_k)\). Then the prediction, and the residual, are integers too:

\[\hat{x}[n] = \left\lfloor \frac{1}{2^s} \sum_{k=1}^{K} q_k \, x[n-k] \right\rfloor, \qquad e[n] = x[n] - \hat{x}[n].\]

The coder stores the \(K\) quantized coefficients, the first \(K\) samples of the block (verbatim), and the residual. Just as in Chapter 14, the decoder inverts this with a recursive filter, \(x[n] = e[n] + \hat{x}[n]\). Since \(\hat{x}[n]\) depends only on samples that the decoder has already reconstructed exactly, the round trip is perfect.

Why does this help? A good predictor leaves a residual that is much smaller than the signal, and concentrated near zero:

Fig. 128 A block of a drum loop as 16-bit integers, its prediction by a quantized LPC model with \(K = 8\), and the residual (left). The residual’s values are concentrated near zero (right).

Rice coding. Small residual values are much more common than large ones, often exponentially so. Rice coding is a simple and fast lossless code that is nearly optimal for such data [Ric79]. Unlike Huffman coding, which needs a table of probabilities for a fixed alphabet, Rice coding operates directly on integers, with a single parameter \(k\).

First, since residuals can be negative, Rice coding maps them to non-negative integers with a zigzag encoding, \(0, -1, 1, -2, 2, \ldots \mapsto 0, 1, 2, 3, 4, \ldots\):

\[\begin{split}u = \begin{cases} 2v & v \geq 0 \\ -2v - 1 & v < 0. \end{cases}\end{split}\]

Then, it splits \(u\) into a quotient and a remainder, \(q = \lfloor u / 2^k \rfloor\) and \(r = u \bmod 2^k\). It writes the quotient in unary, as \(q\) ones followed by a zero (\(0 \to \texttt{0}\), \(1 \to \texttt{10}\), \(2 \to \texttt{110}\), …), and the remainder in ordinary binary, with exactly \(k\) bits. Small values have small quotients, and therefore short codes. For example, with \(k = 4\):

• \(v = 5\): \(u = 10\), so \(q = 0\) and \(r = 10\). The code is \(\texttt{0}\,|\,\texttt{1010}\), \(5\) bits.

• \(v = -60\): \(u = 119\), so \(q = 7\) and \(r = 7\). The code is \(\texttt{11111110}\,|\,\texttt{0111}\), \(12\) bits.

The code below implements Rice coding for \(k = 4\):

def zigzag(v: int) -> int:
    """Interleave signed integers as 0, -1, 1, -2, 2, ... -> 0, 1, 2, 3, 4, ..."""
    return 2 * v if v >= 0 else -2 * v - 1


def rice_encode(v: int, k: int) -> str:
    """The Rice code of a signed integer v with parameter k, as a string of bits."""
    u = zigzag(v)
    q, r = u >> k, u & ((1 << k) - 1)        # quotient and remainder of u / 2^k
    unary = "1" * q + "0"                      # q ones, then a terminating zero
    binary = format(r, f"0{k}b") if k else ""  # the remainder, in exactly k bits
    return unary + binary


for v in [0, -1, 1, 5, -60]:
    code = rice_encode(v, k=4)
    print(f"{v:4d} -> {code} ({len(code)} bits)")

   0 -> 00000 (5 bits)
  -1 -> 00001 (5 bits)
   1 -> 00010 (5 bits)
   5 -> 01010 (5 bits)
 -60 -> 111111100111 (12 bits)

The FLAC coder chooses \(k\) for each block to minimize the total length of its codes. Larger residuals call for a larger \(k\). Putting all of these steps together, the code below implements a toy version of FLAC. It compresses our drum loop to fewer than \(10\) bits per sample (a compression factor of about \(1.8\times\)), and verifies that every block decodes perfectly:

# hide
import numpy as np
import pyquist as pq

Show setup

Hide setup

# collapse
def rice_bits(u: np.ndarray, k: int) -> int:
    """Total length of the Rice codes of non-negative integers u: a unary quotient
    (u >> k ones plus a zero), then a k-bit remainder."""
    return int(np.sum((u >> k) + 1 + k))


def zigzag(e: np.ndarray) -> np.ndarray:
    return np.where(e >= 0, 2 * e, -2 * e - 1)

def encode_block(x: np.ndarray, K: int = 8, shift: int = 12) -> tuple[np.ndarray, np.ndarray]:
    """Fit LPC coefficients to a block of integer samples, quantize them to integers
    (in units of 2^-shift), and return them along with the integer residual."""
    N = len(x)
    X = np.stack([x[K - k:N - k] for k in range(1, K + 1)], axis=1)  # row n: x[n-1], ..., x[n-K]
    h, *_ = np.linalg.lstsq(X.astype(float), x[K:].astype(float), rcond=None)
    q = np.round(h * 2 ** shift).astype(np.int64)                    # quantized coefficients
    prediction = (X @ q) >> shift                                    # integer arithmetic only
    return q, x[K:] - prediction


def decode_block(warmup: np.ndarray, q: np.ndarray, residual: np.ndarray,
                 shift: int = 12) -> np.ndarray:
    """Invert encode_block: x[n] = e[n] + prediction, one sample at a time (recursive)."""
    x = [int(v) for v in warmup]
    q = [int(c) for c in q]
    for e in residual:
        prediction = sum(c * x[-k] for k, c in enumerate(q, start=1)) >> shift
        x.append(int(e) + prediction)
    return np.array(x)


drums = pq.Audio.from_file("./assets/audio-drums.wav")
x = np.round(drums.samples[:, 0] * 32768).astype(np.int64)  # back to 16-bit integers
B, K = 4096, 8
bits = 0
for start in range(0, len(x) - K, B):
    block = x[start:start + B]
    q, residual = encode_block(block, K)
    k = min(range(16), key=lambda k: rice_bits(zigzag(residual), k))  # best k per block
    bits += 16 * K + 16 * K + 4 + rice_bits(zigzag(residual), k)      # warm-up, coefficients, k, residual
    assert np.array_equal(decode_block(block[:K], q, residual), block)  # lossless!

print(f"Uncompressed: 16 bits per sample")
print(f"Toy FLAC: {bits / len(x):.2f} bits per sample (compression factor {16 * len(x) / bits:.2f}×)")

Uncompressed: 16 bits per sample
Toy FLAC: 9.07 bits per sample (compression factor 1.76×)