Entropy Coding and Compression¶
Construct Huffman codes from symbol probabilities and check that encoding followed by decoding recovers the data exactly. Group symbols into blocks to see how average code lengths approach the entropy bound.
Run this tour¶
Run the cells in order with a Python 3 kernel. The first cell locates the companion data and toolbox and installs missing dependencies when needed. All worked examples include their implementation directly in this notebook. Random seeds make comparisons reproducible; you can change them to explore other samples.
# Locate the companion toolbox locally, or fetch it for a standalone/Colab copy.
from pathlib import Path
import importlib.util
import os
import subprocess
import sys
working = Path.cwd()
candidates = [working, working / "python", working.parent / "python"]
python_dir = next((p for p in candidates if (p / "nt_toolbox").is_dir()), None)
if python_dir is None:
checkout = working / "numerical-tours-support"
if not checkout.exists():
subprocess.run(
[
"git",
"clone",
"--depth",
"1",
"--branch",
"master",
"https://github.com/gpeyre/numerical-tours.git",
str(checkout),
],
check=True,
)
python_dir = checkout / "python"
os.chdir(python_dir)
if str(python_dir) not in sys.path:
sys.path.insert(0, str(python_dir))
requirements = python_dir / "requirements.txt"
if any(
importlib.util.find_spec(name) is None
for name in [
"numpy",
"scipy",
"matplotlib",
"skimage",
"sklearn",
"pywt",
"ipywidgets",
"cvxpy",
"skfmm",
"autograd",
"progressbar",
"celer",
]
):
subprocess.run(
[sys.executable, "-m", "pip", "install", "-r", str(requirements)], check=True
)
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(0)
plt.rcParams.update(
{
"figure.figsize": (8, 4),
"figure.dpi": 100,
"axes.spines.top": False,
"axes.spines.right": False,
"font.size": 11,
"image.cmap": "gray",
}
)
%matplotlib inline
$\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$
This numerical tour studies source coding using entropic coders (Huffman and arithmetic).
import numpy as np
import scipy as scp
import pylab as pyl
import matplotlib.pyplot as plt
from nt_toolbox.general import np, plt, pylab
from nt_toolbox.signal import np, plt, pylab
%matplotlib inline
Source Coding and Entropy¶
Entropic coding converts a vector $x$ of integers into a binary stream $y$. Entropic coding exploits the redundancies in the statistical distribution of the entries of $x$ to reduce as much as possible the size of $y$. The lower bound for the number of bits $p$ of $y$ is the Shannon bound :
$$p=-\sum_ih(i)\log_2(h(i))$$
where $h(i)$ is the probability of apparition of symbol $i$ in $x$.
Fist we generate a simple binary signal $x$ so that $0$ has a probability $p$ to appear in $x$.
Probability of 0.
p = 0.1
Size.
n = 512
Signal, should be with token 1,2.
from numpy import random
x = (random.rand(n) > p) + 1
One can check the probabilities by computing the empirical histogram.
h = [np.sum(x == 1), np.sum(x == 2)]
h = h / np.sum(h)
print("Empirical p = %.2f" % h[0])
We can compute the entropy of the distribution represented as a vector $h$ of proability that should sum to 1. We take a max to avoid problems with null probabilties.
e = -np.sum(h * np.log2([max(e, 1e-20) for e in h]))
print("Entropy = %.2f" % e)
Huffman Coding¶
A Hufman code $C$ associates to each symbol $i$ in $\{1,...,m\}$ a binary code $C_i$ whose length is as close as possible to the optimal bound $-\log_2\left(h(i)\right)$, where $h(i)$ is the probability of apparition of the symbol $i$.
We select a set of proabilities.
h = [0.1, 0.15, 0.4, 0.15, 0.2]
The tree $T$ contains the codes and is generated by an iterative algorithm. The initial "tree" is a collection of empty trees, pointing to the symbols numbers.
m = len(h)
T = [0] * m # create an empty tree
We build iteratively the Huffman tree by grouping together the two erees that have the smallest probabilities. The merged tree has a probability which is the sum of the two selected probabilities.
Initial probability.
# we use the symbols i = 0,1,2,3,4 (as strings) with the associated probabilities h(i)
for i in range(m):
T[i] = (h[i], str(i))
Iterative merging of the leading probabilities.
while len(T) > 1:
T.sort(
key=lambda node: node[0]
) # sort according to the first values of the tuples (the probabilities)
t = tuple(T[:2])
q = T[0][0] + T[1][0]
T = T[2:] + [(q, t)]
We trim the computed tree by removing the probabilities.
def trim(T):
T0 = T[1]
if type(T0) == str:
return T0
else:
return (trim(T0[0]), trim(T0[1]))
T = trim(T[0])
We display T using the ete3 package (install it in the terminal with "pip install ete3").
fig, ax = plt.subplots(figsize=(9, 4))
def draw_tree(tree, x=0.5, y=0.0, width=0.45):
if isinstance(tree, str):
ax.text(
x,
y,
f"symbol {tree}",
ha="center",
va="center",
bbox={"boxstyle": "round", "facecolor": "#dbeafe"},
)
return
for bit, child in enumerate(tree):
xc = x + (2 * bit - 1) * width / 2
ax.plot([x, xc], [y, y - 1], color="#334155")
ax.text((x + xc) / 2, y - 0.45, str(bit), ha="center")
draw_tree(child, xc, y - 1, width / 2)
draw_tree(T)
ax.axis("off")
plt.show()
Once the tree $T$ is computed, one can compute the code $C_{i}$ associated to each symbol $i$. This requires to perform a deep first search in the tree and stop at each node.
codes = {}
def huffman_gencode(T, codes, c):
if type(T) == str: # test if T is a leaf
codes[T] = c
else:
huffman_gencode(T[0], codes, c + "0")
huffman_gencode(T[1], codes, c + "1")
huffman_gencode(T, codes, "")
Display the code.
for e in codes:
print("Code of token " + e + ": " + codes[e])
We draw a vector $x$ according to the distribution $h$.
Size of the signal.
n = 1024
Randomization.
from numpy import random
def rand_discr(p, m=1):
"""
rand_discr - discrete random generator
y = rand_discr(p, n);
y is a random vector of length n drawn from
a variable X such that
p(i) = Prob( X=i )
Copyright (c) 2004 Gabriel Peyré
"""
# makes sure it sums to 1
p = p / np.sum(p)
n = len(p)
coin = random.rand(m)
cumprob = np.append(0, +np.cumsum(p))
sample = np.zeros(m)
for j in range(n):
ind = (coin > cumprob[j]) & (coin <= cumprob[j + 1])
sample[ind] = j
return sample
x = rand_discr(h, n)
Worked example 1
Implement the coding of the vector $x$ to obtain a binary vector $y$, which corresponds to replacing each sybmol $x(i)$ by the code $C_{x(i)}$.
y = ""
for e in x:
y = y + codes[str(int(e))]
Compare the length of the code with the entropy bound.
e = -np.sum(h * np.log2([max(e, 1e-20) for e in h]))
print("Entropy bound = %.2f" % (n * e))
print("Huffman code = %.2f" % len(y))
Decoding is more complicated, since it requires to iteratively parse the tree $T$.
Initial empty decoded stream.
x1 = []
Perform decoding.
T0 = T
for e in y:
if e == "0":
T0 = T0[0]
else:
T0 = T0[1]
if type(T0) == str:
i = i + 1
x1 += T0
T0 = T
We test if the decoding is correct.
from numpy import linalg
err = linalg.norm(np.subtract(x, [float(e) for e in x1]))
print("Error (should be zero) : %f " % err)
Huffman Block Coding¶
A Huffman coder is inefficient because it can distribute only an integer number of bit per symbol. In particular, distribution where one of the symbol has a large probability are not well coded using a Huffman code. This can be aleviated by replacing the set of $m$ symbols by $m^q$ symbols obtained by packing the symbols by blocks of $q$ (here we use $m=2$ for a binary alphabet). This breaks symbols with large probability into many symbols with smaller proablity, thus approaching the Shannon entropy bound.
Generate a binary vector with a high probability of having 1, so that the Huffman code is not very efficient (far from Shanon bound).
Proability of having 0.
t = 0.12
Probability distriution.
h = [t, 1 - t]
Generate signal.
from numpy import random
n = 4096 * 2
x = (random.rand(n) > t) + 1
For block of length $q=3$, create a new vector by coding each block with an integer in $\{1,...,m^q=2^3\}$. The new length of the vector is $n_1/q$ where $n_1=\lceil n/q\rceil q$.
Block size.
q = 3
Maximum token value.
m = 2
New size.
n1 = (n // q + 1) * q
New vector.
x1 = np.zeros(n1)
x1[: len(x)] = x
x1[len(x) :] = 1
x1 = x1 - 1
x2 = []
for i in range(0, n1, q):
mult = [m**j for j in range(q)]
x2.append(sum(x1[i : i + q] * mult))
We generate the probability table $H$ of $x_1$ that represents the probability of each new block symbols in $\{1,...,m^q\}$.
H = h
for i in range(q - 1):
Hold = H
H = []
for j in range(len(h)):
H = H + [e * h[j] for e in Hold]
A simpler way to compute this block-histogram is to use the Kronecker product.
H = h
for i in range(1, q):
H = np.kron(H, h)
Worked example 2
For various values of block size $k$, Perform the Huffman coding and compute the length of the code. Compare with the entropy lower bound.
e_bound = -np.sum(h * np.log2([max(e, 1e-20) for e in h]))
print("Entropy bound = %f" % e_bound)
print("---")
Err = []
for k in range(1, 11):
# define constants
m_token = 2
# new size
n1 = (n // k + 1) * k
# new vector
x1 = np.zeros(n1)
x1[: len(x)] = x
x1[len(x) :] = 1
x1 = x1 - 1
x2 = []
for i in range(0, n1, k):
mult = [m_token**i for i in range(k)]
x2.append(sum(x1[i : i + k] * mult))
# new probability distribution
H = h
for i in range(1, k):
H = np.kron(H, h)
# build Huffman tree
m = len(H)
T = [0] * m
for i in range(m):
T[i] = (H[i], str(i))
while len(T) > 1:
T.sort(key=lambda node: node[0])
t = tuple(T[:2])
q = T[0][0] + T[1][0]
T = T[2:] + [(q, t)]
T = trim(T[0])
# find the codes
codes = {}
huffman_gencode(T, codes, "")
# encode
y = ""
for e in x2:
y = y + codes[str(int(e))]
# append error
err = len(y) / len(x)
print("Huffman(block size = %i) = %f" % (k, err))
Err.append(err - e_bound)
plt.figure(figsize=(7, 5))
plt.plot(Err, linewidth=2)
plt.title("Huffman block coding performance")
plt.xlabel("Block size q")
plt.ylabel("Huffman error $-$ Entropy bound")
plt.show()
References and further reading¶
Thomas M. Cover and Joy A. Thomas. Elements of Information Theory. 2006, 2nd ed., Wiley. Entropy, source coding, and fundamental compression limits.
Claude E. Shannon. A Mathematical Theory of Communication. 1948, Bell System Technical Journal 27, 379–423 and 623–656. Entropy and the source coding theorem.
David A. Huffman. A Method for the Construction of Minimum-Redundancy Codes. 1952, Proceedings of the IRE 40(9), 1098–1101. The greedy construction of optimal binary prefix codes.
Stéphane Mallat. A Wavelet Tour of Signal Processing: The Sparse Way. 2009, 3rd ed., Academic Press. Multiresolution analysis, sparse approximation, and wavelet algorithms.
Gabriel Peyré. Advanced Signal, Image and Surface Processing. 2010, course notes. A mathematical companion to the Numerical Tours.