Maximum Entropy

Motivation: the temperature of a gas

The temperature of a gas corresponds to the average kinetic energy of its molecules. What can we say about the distribution of velocities at a given temperature? Physics tells us it is the maximum entropy distribution under the temperature constraint — the Maxwell–Boltzmann distribution.

IntuitionMacrostates and microstates

The maximum entropy distribution corresponds to the macrostate (indexed by the empirical distribution) that has the most microstates (the individual gas velocities). Implicit in the use of maximum entropy methods in physics is a sort of AEP: all microstates are equally probable, so the macrostate with the most of them is the one you will observe.

This is the same counting we did in the method of types: polynomially many types, exponentially many sequences per type, and $\lvert T(P)\rvert \doteq 2^{nH(P)}$ — so the type with the largest entropy has the most sequences. Maximum entropy is a counting principle before it is a modelling principle.

ImportantWhat this chapter adds

The conditional limit theorem already showed that maximum entropy is a theorem about conditioning, not just a modelling choice. This chapter gives the general solution: given any set of moment constraints, what exactly is the maximizing density, and when does it fail to exist?

The maximum entropy distribution

The problem. Maximize $h(f)$ over all probability densities $f$ satisfying

$f(x) \ge 0$ , with equality outside the support set $S$
$\int_S f(x)\,dx = 1$
$\int_S f(x) r_i(x)\,dx = \alpha_i$ for $1 \le i \le m$

So $f$ is a density on $S$ meeting moment constraints $\alpha_1, \ldots, \alpha_m$ .

Approach 1: calculus (suggests the form)

DerivationLagrange multipliers on the entropy functional

$h(f)$ is a concave function over a convex set. Form the functional

$J(f) = -\int f\ln f + \lambda_0\int f + \sum_{i=1}^m \lambda_i\int f r_i$

and “differentiate” with respect to $f(x)$ , the $x$ th component of $f$ :

$\frac{\partial J}{\partial f(x)} = -\ln f(x) - 1 + \lambda_0 + \sum_{i=1}^m \lambda_i r_i(x)$

Setting this to zero and solving for $f(x)$ gives the exponential form below.

FormulaThe maximizing density

$f(x) = e^{\lambda_0 - 1 + \sum_{i=1}^m \lambda_i r_i(x)}, \qquad x \in S$

with $\lambda_0, \lambda_1, \ldots, \lambda_m$ chosen so that $f$ satisfies the constraints.

WarningCalculus only suggests the form

This argument produces a candidate. To prove it is the maximum we would have to take the second variation. It is far simpler to use the information inequality $D(g\Vert f)\ge0$ , which is Approach 2.

Approach 2: the information inequality (proves it)

TheoremMaximum entropy distribution

Let $f^*(x) = f_\lambda(x) = e^{\lambda_0 + \sum_{i=1}^m\lambda_i r_i(x)}$ , $x\in S$ , with the $\lambda_i$ chosen so $f^*$ satisfies the constraints. Then $f^*$ uniquely maximizes $h(f)$ over all densities satisfying the constraints.

ProofInsert

f^*

into the entropy, then use

D \ge 0

Let $g$ satisfy the constraints. Multiply and divide inside the logarithm by $f^*$ :

$h(g) = -\int_S g\ln g = -\int_S g\ln\frac{g}{f^*}f^* = -D(g\Vert f^*) - \int_S g\ln f^*$

$\overset{(a)}{\le} -\int_S g\ln f^* \overset{(b)}{=} -\int_S g\left(\lambda_0 + \sum_i\lambda_i r_i\right) \overset{(c)}{=} -\int_S f^*\left(\lambda_0 + \sum_i \lambda_i r_i\right) = -\int_S f^*\ln f^* = h(f^*)$

where (a) is non-negativity of relative entropy, (b) is the definition of $f^*$ , and (c) is the fact that $g$ and $f^*$ satisfy the same constraints — so the integral of any linear combination of $1$ and the $r_i$ is the same under both. Equality in (a) holds iff $g = f^*$ almost everywhere, proving uniqueness.

∎

IntuitionWhy the proof works so cleanly

The whole argument turns on step (c): $\ln f^*$ is, by construction, a linear combination of exactly the functions the constraints pin down. So the cross term $\int g\ln f^*$ cannot tell $g$ from $f^*$ , and the only surviving difference is $-D(g\Vert f^*) \le 0$ .

This is the identical trick used in the Gaussian max-entropy theorem, where $\log\phi_K$ was a quadratic form and the constraint was a covariance. That theorem is the special case $r_i(x) = x_ix_j$ . The same approach holds for discrete entropies and for multivariate distributions.

Examples

ExampleOne-dimensional gas with a temperature constraint

Constraints $\mathbb{E}X = 0$ , $\mathbb{E}X^2 = \sigma^2$ give the form $f(x) = e^{\lambda_0 + \lambda_1 x + \lambda_2 x^2}$ . The exponent is a quadratic, so by inspection this is a normal density; matching the constraints:

$f(x) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{x^2}{2\sigma^2}}$

ExampleDice, no constraints

$S = \{1,\ldots,6\}$ with no moment constraints: the maximizer is the uniform distribution $p(x) = \tfrac16$ . (Setting all $\lambda_i = 0$ leaves a constant.)

ExampleBounded support

$S = [a,b]$ , no other constraints: the uniform distribution on that range.

ExamplePositive half-line with a mean constraint — the atmosphere

$S = [0,\infty)$ and $\mathbb{E}X = \mu$ give the exponential density

$f(x) = \frac{1}{\mu}e^{-x/\mu}, \qquad x \ge 0$

Physical reading: consider the height $X$ of molecules in the atmosphere. The average potential energy is fixed, so the gas tends to the maximum entropy distribution subject to $\mathbb{E}(mgX)$ fixed — an exponential. The density of the atmosphere does indeed have this distribution.

WarningA mean constraint on the whole line gives no maximizer

$S = (-\infty,\infty)$ with only $\mathbb{E}X = \mu$ : the maximum entropy is infinite and there is no maximum entropy distribution. (Consider normal distributions with larger and larger variances — all satisfy the mean constraint, with entropy growing without bound.) A first moment alone does not confine the distribution to a finite effective volume.

ExampleWhole line, first and second moments

$\mathbb{E}X = \alpha_1$ , $\mathbb{E}X^2 = \alpha_2$ give $\mathcal{N}(\alpha_1, \alpha_2 - \alpha_1^2)$ .

ExampleMultivariate covariance constraint

$S = \mathbb{R}^n$ , $\mathbb{E}X_iX_j = K_{ij}$ : the density is $f(x) = e^{\lambda_0 + \sum_{i,j}\lambda_{ij}x_ix_j}$ . The exponent is a quadratic form, so by inspection it is a multivariate normal with zero mean; matching the second moments fixes the covariance:

$f(x) = \frac{1}{\sqrt{2\pi}^{\,n}\lvert K\rvert^{1/2}}e^{-\frac12 x^TK^{-1}x}, \qquad h = \frac12\log(2\pi e)^n\lvert K\rvert$

ExamplePartial covariance constraints force zeros in $K^{-1}$

Suppose we know $K_{ij}$ only for $(i,j)\in A$ — say only for $i = j\pm2$ . Comparing the general form $e^{\lambda_0 + \sum\lambda_{ij}x_ix_j}$ with the Gaussian density $e^{-\frac12 x^TK^{-1}x}$ , the $\lambda_{ij}$ are the entries of $-\tfrac12 K^{-1}$ — and unconstrained pairs have $\lambda_{ij} = 0$ . Hence

$\big(K^{-1}\big)_{ij} = 0 \quad \text{for } (i,j)\in A^c$

The entries of the inverse covariance vanish outside the constraint set. This is the information-theoretic origin of sparse precision matrices and of conditional-independence structure in Gaussian graphical models.

ExampleDice with a mean constraint — Boltzmann’s example

Throw $n$ dice and be told the total number of spots is $n\alpha$ . What proportion show face $i$ ?

The counting approach. There are $\binom{n}{n_1,\ldots,n_6}$ ways for $n_i$ dice to show face $i$ — a macrostate indexed by $(n_1,\ldots,n_6)$ corresponding to that many microstates, each of probability $6^{-n}$ . To find the most probable macrostate, maximize the multinomial subject to $\sum_{i=1}^6 in_i = n\alpha$ .

Using a crude Stirling approximation $n!\approx(n/e)^n$ :

$\binom{n}{n_1,\ldots,n_6} \approx \frac{(n/e)^n}{\prod_{i=1}^6 (n_i/e)^{n_i}} = \prod_{i=1}^6\left(\frac{n}{n_i}\right)^{n_i} = e^{nH\left(\frac{n_1}{n},\ldots,\frac{n_6}{n}\right)}$

So maximizing the multinomial under the constraint is almost equivalent to maximizing $H(p_1,\ldots,p_6)$ under $\sum ip_i = \alpha$ . By the theorem,

$p_i^* = \frac{e^{\lambda i}}{\sum_{j=1}^6 e^{\lambda j}}$

with $\lambda$ chosen so $\sum ip_i^* = \alpha$ . The most probable macrostate is $(np_1^*, \ldots, np_6^*)$ .

The approximations are justified. By the conditional limit theorem, not only is the maximum entropy macrostate the most likely, it contains almost all the probability: for rational $\alpha$ ,

$\Pr\left\{\left\lvert\frac{N_i}{n} - p_i^*\right\rvert < \epsilon,\ i = 1,\ldots,6 \ \middle|\ \frac1n\sum_i X_i = \alpha\right\} \to 1$

as $n\to\infty$ along the subsequence where $n\alpha$ is an integer.

The anomalous maximum entropy problem

The problem. We proved the maximizer is $f(x) = e^{\lambda_0 + \sum\lambda_i h_i(x)}$ if the $\lambda_i$ satisfying the constraints exist. What if they don’t?

Consider maximizing entropy subject to four constraints on the whole line:

$\int f = 1, \qquad \int xf = \alpha_1, \qquad \int x^2 f = \alpha_2, \qquad \int x^3 f = \alpha_3$

WarningThe method appears to fail

The maximizer, if it exists, must be $f(x) = e^{\lambda_0 + \lambda_1 x + \lambda_2 x^2 + \lambda_3 x^3}$ . But if $\lambda_3 \ne 0$ then $\int_{-\infty}^\infty f = \infty$ and the density cannot be normalized — a cubic exponent blows up on one side. So $\lambda_3$ must be 0. But then we have four equations and only three free variables, so in general the constants cannot be chosen.

ImportantThe resolution: the supremum is not attained

The entropy has a least upper bound under these constraints, but it is not possible to attain it.

With only first and second moment constraints, the maximizer is the normal with those moments. Adding a third moment constraint cannot raise the maximum. Can we achieve that same value? We cannot achieve it, but we can come arbitrarily close.

IntuitionThe wiggle construction

Take a normal distribution and add a small “wiggle” at a very high value of $x$ . Because $x^3$ grows so fast out there, a tiny amount of probability mass placed far enough out moves the third moment a lot while barely disturbing the first and second — and barely disturbing the entropy. Any residual change to the first two moments can be cancelled by adding further wiggles to balance them.

By choosing the positions of the wiggles we can hit any value of the third moment without reducing the entropy significantly below that of the associated normal. Hence

$\sup_f h(f) = h\big(\mathcal{N}(0, \alpha_2 - \alpha_1^2)\big) = \frac12\ln 2\pi e(\alpha_2 - \alpha_1^2)$

NoteThe lesson

The maximum entropy may be only $\epsilon$ -achievable. The theorem’s hypothesis — “if $\lambda_0,\ldots,\lambda_m$ satisfying the constraints exist” — is not a technicality to be waved through. When a constraint function grows faster than the exponential tail can be normalized against, the supremum is approached but never reached.

Spectrum estimation

The problem. Given a stationary zero-mean process $\{X_i\}$ , the autocorrelation function is $R(k) = \mathbb{E}X_iX_{i+k}$ , and its Fourier transform is the power spectral density

$S(\lambda) = \sum_{m=-\infty}^{\infty}R(m)e^{-im\lambda}, \qquad -\pi < \lambda \le \pi$

Since the spectrum indicates the structure of the process, we want to estimate it from a sample.

WarningWhy the naive (periodogram) method fails

The simplest approach estimates the autocorrelation by sample averages,

$\hat R(k) = \frac{1}{n-k}\sum_{i=1}^{n-k}X_iX_{i+k}$

and feeds all of them into the transform. The resulting estimate does not converge to the true spectrum for large $n$ , so this periodogram method is rarely used.

The reason: the estimates have very different accuracies. Low lags $k$ are based on many samples; high $k$ on very few. Truncating — setting high-lag autocorrelations to 0 — introduces artifacts from the sudden transition, and various windowing schemes that smooth the transition reduce spectral resolution and can produce negative power estimates.

ImportantBurg’s idea

In the late 1960s, working on geophysical spectral estimation, Burg suggested: instead of setting the high-lag autocorrelations to zero, set them to the values that make the fewest assumptions about the data — i.e. the values that maximize the entropy rate of the process.

Setting them to zero is an assumption (that the process decorrelates abruptly); maximizing entropy is the absence of one. This is the maximum entropy principle as articulated by Jaynes.

Burg assumed the process stationary and Gaussian, and found that the maximizer subject to the correlation constraints is an autoregressive Gaussian process of the appropriate order. Where an underlying autoregressive model is plausible, this determines the model parameters — as in linear predictive coding for speech. The method is known as the maximum entropy method or Burg’s method.

Entropy rates of a Gaussian process

DefinitionDifferential entropy rate

$h(\mathcal{X}) = \lim_{n\to\infty}\frac{h(X_1,\ldots,X_n)}{n}$

if the limit exists. As in the discrete case, for stationary processes the limit exists and equals both

$h(\mathcal{X}) = \lim_{n\to\infty}\frac{h(X_1,\ldots,X_n)}{n} = \lim_{n\to\infty}h(X_n\mid X_{n-1},\ldots,X_1)$

For a stationary Gaussian process,

$h(X_1,\ldots,X_n) = \frac12\log(2\pi e)^n\lvert K^{(n)}\rvert$

where $K^{(n)}$ is Toeplitz with $R(0), R(1), \ldots, R(n-1)$ along the top row, i.e. $K^{(n)}_{ij} = R(i-j)$ .

ImportantKolmogorov’s formula

As $n\to\infty$ the density of eigenvalues of the covariance matrix tends to the spectrum of the process, and the entropy rate is

$h(\mathcal{X}) = \frac12\log2\pi e + \frac{1}{4\pi}\int_{-\pi}^{\pi}\log S(\lambda)\,d\lambda$

IntuitionSame Toeplitz-to-spectrum move as in the Gaussian channel

The step “eigenvalue density of a Toeplitz matrix $\to$ power spectrum” is exactly the one that turned eigenvalue water-filling into water-filling in the spectral domain. The entropy rate is (up to constants) the average of $\log S$ over frequency — a geometric mean of the spectrum.

FormulaEntropy rate as prediction error

The entropy rate is also $\lim_n h(X_n\mid X^{n-1})$ . Since the process is Gaussian the conditional distribution is Gaussian, so this conditional entropy is $\frac12\log2\pi e\sigma_\infty^2$ , where $\sigma_\infty^2$ is the variance of the error in the best estimate of $X_n$ given the infinite past. Hence

$\sigma_\infty^2 = \frac{1}{2\pi e}2^{2h(\mathcal{X})}$

IntuitionWhat this identity says

The entropy rate corresponds to the minimum mean-squared error of the best predictor of the process given its infinite past. A process you can predict well has low entropy rate. Note this is precisely the estimation bound $\mathbb{E}(X-\hat X)^2 \ge \frac{1}{2\pi e}e^{2h}$ — holding with equality here, because the process is Gaussian and the conditional mean is the estimator.

Burg’s maximum entropy theorem

TheoremBurg

The maximum entropy rate stochastic process $\{X_i\}$ satisfying

$\mathbb{E}X_iX_{i+k} = \alpha_k, \qquad k = 0,1,\ldots,p, \quad \text{for all } i$

is the $p$ th-order Gauss–Markov process

$X_i = -\sum_{k=1}^p a_kX_{i-k} + Z_i$

where the $Z_i$ are i.i.d. $\sim\mathcal{N}(0,\sigma^2)$ and $a_1,\ldots,a_p,\sigma^2$ are chosen to satisfy the constraints.

ImportantWhat is not assumed

We do not assume $\{X_i\}$ is (a) zero mean, (b) Gaussian, or (c) wide-sense stationary. All three are conclusions: the maximizer turns out to be Gaussian and autoregressive, over the class of all processes meeting the correlation constraints.

ProofTwo bounds: Gaussian beats anything, and

p

th-order Markov beats Gaussian

Let $X_1,\ldots,X_n$ be any process satisfying the constraints, and let $Z_1,\ldots,Z_n$ be a Gaussian process with the same covariance matrix.

Step 1: Gaussian is the best for a given covariance. By the multivariate max-entropy theorem,

$h(X_1,\ldots,X_n) \le h(Z_1,\ldots,Z_n)$

Step 2: Chain rule, then truncate the conditioning.

$h(Z_1,\ldots,Z_n) = h(Z_1,\ldots,Z_p) + \sum_{i=p+1}^n h(Z_i\mid Z_{i-1},\ldots,Z_1) \le h(Z_1,\ldots,Z_p) + \sum_{i=p+1}^n h(Z_i\mid Z_{i-1},\ldots,Z_{i-p})$

by the chain rule and conditioning-reduces-entropy. Truncating the past to $p$ terms is where the $p$ th-order structure enters — and it is an upper bound, so it points the right way.

Step 3: Replace by a genuine $p$ th-order process. Define $Z_1',\ldots,Z_n'$ as a $p$ th-order Gauss–Markov process with the same distribution as $Z_1,\ldots,Z_n$ for all orders up to $p$ (existence verified by the Yule–Walker equations below). Since $h(Z_i\mid Z_{i-1},\ldots,Z_{i-p})$ depends only on the $p$ th-order distribution, these conditional entropies are equal:

$h(X_1,\ldots,X_n) \le h(Z_1',\ldots,Z_p') + \sum_{i=p+1}^n h(Z_i'\mid Z_{i-1}',\ldots,Z_{i-p}') = h(Z_1',\ldots,Z_n')$

where the last equality is the $p$ th-order Markovity of $\{Z_i'\}$ — for this process the truncated conditioning is not an approximation but exact.

Step 4: Divide and take the limit.

$\lim\frac1n h(X_1,\ldots,X_n) \le \lim\frac1n h(Z_1',\ldots,Z_n') = h^* = \frac12\log2\pi e\sigma^2$

∎

IntuitionThe proof in one line

The entropy of a finite segment is bounded above by that of a Gaussian process with the same covariance, which is in turn bounded above by that of the minimal order Gauss–Markov process satisfying the same covariance constraints. Two successive relaxations, each of which the candidate process attains with equality.

The Yule–Walker equations

The question left open. Given covariances $R(0),\ldots,R(p)$ , does a $p$ th-order Gauss–Markov process with those covariances exist, and can we find its coefficients?

DerivationMultiply the recursion by

X_{i-l}

and take expectations

From $X_i = -\sum_{k=1}^p a_kX_{i-k} + Z_i$ , multiplying by $X_{i-l}$ and using $R(k) = R(-k)$ :

$R(0) = -\sum_{k=1}^p a_kR(-k) + \sigma^2, \qquad R(l) = -\sum_{k=1}^p a_kR(l-k), \quad l = 1, 2, \ldots$

These are the Yule–Walker equations: $p+1$ equations in the $p+1$ unknowns $a_1,\ldots,a_p,\sigma^2$ . So we can solve for the process parameters from the covariances — answering the existence question.

ImportantYule–Walker extension

Beyond providing a linear system for the $a_k$ and $\sigma^2$ , these equations indicate how the autocorrelations behave for lags greater than $p$ : the high-lag autocorrelations are an extension of the values for lags below $p$ , called the Yule–Walker extension. This is precisely Burg’s answer to the truncation problem — rather than zeroing the unknown lags, extend them by the recursion.

Fast algorithms (Levinson, Durbin) exploit the special structure of these equations to compute the coefficients efficiently. (Set $a_0 = 1$ for consistent notation.)

FormulaThe maximum entropy spectrum

$S(\lambda) = \sum_{m=-\infty}^\infty R(m)e^{-im\lambda} = \frac{\sigma^2}{\left\lvert 1 + \sum_{k=1}^p a_ke^{-ik\lambda}\right\rvert^2}, \qquad -\pi\le\lambda\le\pi$

This is the maximum entropy spectral density subject to $R(0), R(1),\ldots,R(p)$ .

NoteComputing the entropy rate without the

a_i

Let $K_p$ be the autocorrelation matrix with $R_0,\ldots,R_p$ along the top row. Then

$h^* = h(X_p\mid X_{p-1},\ldots,X_0) = h(X_0,\ldots,X_p) - h(X_0,\ldots,X_{p-1})$

$= \frac12\log(2\pi e)^{p+1}\lvert K_p\rvert - \frac12\log(2\pi e)^p\lvert K_{p-1}\rvert = \frac12\log(2\pi e)\frac{\lvert K_p\rvert}{\lvert K_{p-1}\rvert}$

— a ratio of determinants, no coefficients required.

WarningHow many lags? The model-order question

In practice we are given a sample $X_1,\ldots,X_n$ and compute autocorrelations from it. How many lags should we use — what is the optimum $p$ ? Nothing in Burg’s theorem answers this: it tells you the best process for a given $p$ , and larger $p$ always allows higher entropy… on the sample.

A logically sound method is to choose the $p$ minimizing the total description length in a two-stage description of the data (Rissanen, Barron) — pay for the model, then pay for the data given the model. This is closely related to Kolmogorov complexity, and is the principled answer to a question that maximum entropy alone cannot settle.

Quick reference

Constraint set	Maximum entropy distribution
Support $S = [a,b]$ , nothing else	Uniform on $[a,b]$
Finite alphabet, nothing else	Uniform
$S = [0,\infty)$ , $\mathbb{E}X = \mu$	Exponential, $\frac1\mu e^{-x/\mu}$
$S = \mathbb{R}$ , $\mathbb{E}X = \mu$ only	None — entropy unbounded
$S = \mathbb{R}$ , $\mathbb{E}X$ and $\mathbb{E}X^2$	$\mathcal{N}(\alpha_1, \alpha_2-\alpha_1^2)$
$S = \mathbb{R}^n$ , $\mathbb{E}X_iX_j = K_{ij}$	$\mathcal{N}(0,K)$ , $h = \frac12\log(2\pi e)^n\lvert K\rvert$
Covariance known only on $A$	$\mathcal{N}$ with $(K^{-1})_{ij} = 0$ off $A$
Finite alphabet, $\mathbb{E}X = \alpha$	$p_i^*\propto e^{\lambda i}$ (exponential tilt)
$S = \mathbb{R}$ , first three moments	None attained; $\sup h = h(\mathcal{N}(0,\alpha_2-\alpha_1^2))$
Process, $\mathbb{E}X_iX_{i+k} = \alpha_k$ , $k\le p$	$p$ th-order Gauss–Markov (Burg)

Facts to keep to hand:

The general form is always $f(x) = e^{\lambda_0 + \sum_i\lambda_ir_i(x)}$ — an exponential family whose sufficient statistics are exactly the constrained moments.
The proof is one move: $\ln f^*$ is a linear combination of the constrained functions, so $\int g\ln f^* = \int f^*\ln f^*$ , leaving only $-D(g\Vert f^*)\le0$ .
Existence is a real hypothesis. If the required $\lambda_i$ do not exist (cubic moment on the whole line), the supremum is only $\epsilon$ -achievable via distant wiggles.
Gaussian process facts: $h(\mathcal{X}) = \frac12\log2\pi e + \frac{1}{4\pi}\int\log S(\lambda)d\lambda$ (Kolmogorov), and $\sigma_\infty^2 = \frac{1}{2\pi e}2^{2h(\mathcal{X})}$ — entropy rate is best-prediction error.
Burg: max entropy rate under $p$ correlation constraints $\Rightarrow$ $p$ th-order Gauss–Markov, with parameters from Yule–Walker and high-lag correlations given by the Yule–Walker extension rather than set to zero.
Maximum entropy does not choose the model order $p$ ; description length does.

References

[[lit/cover2006-elements|Cover & Thomas, Elements of Information Theory]] — Chapter 12

Related