Exponential Family

Knowledge Map Atlas
The Exponential Family MLE for the Exponential Family

The Exponential Family

Our treatment of Bayesian statistics showed that certain likelihood-prior pairs such as Binomial-Beta and Normal-Normal produce posteriors in the same family as the prior. This is not a coincidence. These distributions all belong to the exponential family, a unified framework that encompasses the vast majority of distributions used in statistics and machine learning, among them the Bernoulli, Poisson, Normal, Gamma, and Dirichlet distributions.

The exponential family is fundamental for three reasons:

Definition: Exponential Family

The exponential family is a family of probability distributions parameterized by natural parameters (or canonical parameters) \(\eta \in \mathbb{R}^K\) satisfying \(Z(\eta) \lt \infty\), with support over \(\mathcal{X} \subseteq \mathbb{R}^D\) such that \[ \begin{align*} p(x \mid \eta ) &= \frac{1}{Z(\eta)} h(x) \exp\{\eta^\top \mathcal{T}(x)\}\\\\ &= h(x) \exp\{\eta^\top \mathcal{T}(x) - A(\eta)\} \end{align*} \] where

  • \(h(x)\) is a base measure (or carrier measure): a function depending only on \(x\), not on \(\eta\), and often equal to \(1\).
  • \(\mathcal{T}(x) \in \mathbb{R}^K\) is the sufficient statistic vector. It is called sufficient because the likelihood depends on \(x\) only through \(\mathcal{T}(x)\). All information in \(x\) relevant to estimating \(\eta\) is captured by \(\mathcal{T}(x)\).
  • \(Z(\eta)\) is the normalization function, defined so that \(p(x \mid \eta)\) integrates to one over the support.

Each exponential family is defined by different \(h(x)\) and \(\mathcal{T}(x)\).

The normalization function \(Z(\eta)\) is often referred to as the partition function in statistical physics and machine learning. In particular, the log-partition function: \[ A(\eta) = \log Z(\eta) \] is convex over the convex set \(\Omega = \{\eta \in \mathbb{R}^K : A(\eta) \lt \infty\}\), and is strictly convex when the family is minimal. That \(\Omega\) is convex, and that \(A\) is convex on all of \(\Omega\), are consequences of Hölder's inequality that we take for granted here. On the interior of \(\Omega\), which is convex because \(\Omega\) is, the convexity of \(A\), and its strict convexity when the family is minimal, also follow from the Hessian formula \(\nabla^2 A(\eta) = \operatorname{Cov}[\mathcal{T}(x)]\) proved below.

In the context of optimization and information geometry, \(A(\eta)\) acts as a convex potential function that links natural parameters to moment parameters.

Definition: Minimal Representation

An exponential family is said to be minimal if there is no nonzero \(\eta \in \mathbb{R}^K\) such that \[ \eta^\top \mathcal{T}(x) = \text{const} \] holds for \(h(x)\,dx\)-almost every \(x\), that is, for all \(x\) outside a set on which \(h\) integrates to zero (for a discrete support, with the integral replaced by a sum).

Equivalently, the components of the sufficient statistic vector \(\mathcal{T}(x)\) (along with the constant function \(1\)) are linearly independent as functions of \(x\), where two functions that agree \(h(x)\,dx\)-almost everywhere are identified. This property ensures the identifiability of the model. There is a unique, one-to-one mapping between the natural parameters \(\eta\) and the resulting probability distribution.

Math vs. ML Perspectives: Handling the Constant

In many machine learning contexts, the condition is simplified to \[ \eta^\top \mathcal{T}(x) = 0. \] This assumes that the bias (constant) term has already been absorbed into the sufficient statistic vector \(\mathcal{T}(x)\) by augmenting it with a constant \(1\).

However, in our canonical definition where the log-partition function \(A(\eta)\) and base measure \(h(x)\) are explicitly separated, we must use the more rigorous \(\text{const}\) condition. If \(\eta^\top \mathcal{T}(x)\) were a constant, that value would be absorbed by \(A(\eta)\), meaning multiple \(\eta\) values could represent the same distribution. A minimal representation prevents this redundancy, ensuring that the Fisher information matrix is strictly positive definite and the optimization surface is well-behaved.

A common example of a non-minimal representation is the \(K\)-class multinomial distribution. Because the probabilities must sum to one (\(\sum p_i = 1\)), the natural parameters are over-specified. While we can work with this redundant form (as often done in Softmax regression), we can always reparameterize it into a minimal representation using \(K-1\) independent components to restore unique identifiability.

Let \(\eta = f(\phi)\), where \(\phi \in \mathbb{R}^M\) is some other possibly smaller set of parameters (\(M \leq K\)), and then \[ p(x \mid \phi ) = h(x) \exp\{ f(\phi)^\top \mathcal{T}(x) - A(f(\phi))\}. \] If \(M \lt K\) and the mapping \(\phi \to \eta\) is nonlinear, it is said to be a curved exponential family.

If \(\eta = f(\phi) = \phi\), the model is in canonical form and in addition, if \(\mathcal{T} =x\), we call it a natural exponential family (NEF): \[ p(x \mid \eta ) = h(x) \exp\{\eta^\top x - A(\eta)\}. \]

Finally, we define the moment parameters as follows: \[ m = \mathbb{E}[\mathcal{T}(x)] \in \mathbb{R}^K. \]

Example 1: Bernoulli Distribution

\[ \begin{align*} \operatorname{Ber}(x \mid \mu) &= \mu^x (1-\mu)^{1-x} \\\\ &= \exp\{x \log (\mu) + (1-x) \log (1-\mu)\}\\\\ &= \exp\{\mathcal{T}(x)^\top \eta\} \end{align*} \] where

  • \(\mathcal{T}(x) = [\mathbb{1}\{x=1\}, \, \mathbb{1}\{x=0\}]\).
  • \(\eta = [\log(\mu), \, \log(1-\mu)]\).
  • \(\mu\) is the mean parameter.

In this representation, the components \(\mathbb{1}\{x=1\}\) and \(\mathbb{1}\{x=0\}\) of \(\mathcal{T}(x)\) satisfy \(\mathbb{1}\{x=1\} + \mathbb{1}\{x=0\} = 1\) for all \(x\) in the support, violating the minimality condition with the nonzero vector \(v = (1, 1)\). Consequently, shifting \(\eta\) by any multiple of \(v\) yields the same distribution, so \(\eta\) is not uniquely determined. To restore uniqueness, we use a minimal representation: \[ \operatorname{Ber}(x \mid \mu) = \exp\left\{x \log \left(\frac{\mu}{1-\mu}\right) + \log (1-\mu)\right\} \] where

  • \(\mathcal{T}(x) = x\).
  • \(\eta = \log \left(\frac{\mu}{1-\mu}\right)\).
  • \(A(\eta) = -\log (1-\mu) = \log(1+ e^{\eta})\).
  • \(h(x) = 1\).

The Bernoulli example reveals a fundamental connection to classification models. The mean parameter \(\mu\) can be recovered from the canonical parameter \(\eta\) via: \[ \mu = \sigma(\eta) = \frac{1}{1+e^{-\eta}}. \] Here \(\sigma(\cdot)\) is the logistic function.

A direct computation gives \[ \begin{align*} \frac{dA}{d \eta} &= \frac{d}{d\eta} \log(1+ e^{\eta}) \\\\ &= \frac{e^{\eta}}{1 + e^{\eta}} \\\\ &= \frac{1}{1+e^{-\eta}} \\\\ &= \mu. \end{align*} \]

In general, the log-partition function \(A(\eta)\) acts as the cumulant generating function for the sufficient statistics \(\mathcal{T}(x)\). In other words, the derivatives of \(A(\eta)\) generate all the cumulants of \(\mathcal{T}(x)\). The first two are particularly important.

Theorem: Cumulants from the Log-Partition Function

Let \(\eta\) lie in the interior of \(\Omega = \{\eta \in \mathbb{R}^K : A(\eta) \lt \infty\}\). Then the gradient and Hessian of the log-partition function recover the first two cumulants of the sufficient statistic vector: \[ \nabla A(\eta) = \mathbb{E}[\mathcal{T}(x)], \quad \nabla^2 A(\eta) = \operatorname{Cov}[\mathcal{T}(x)]. \]

Proof:

Write \(Z(\eta) = \int h(x) \exp\{\eta^\top \mathcal{T}(x)\}\,dx\), so that \(A(\eta) = \log Z(\eta)\). For \(\eta\) in the interior of \(\Omega\) the integral converges on a neighbourhood of \(\eta\), and we take for granted that differentiation may be carried out under the integral sign there. Differentiating once gives \[ \begin{align*} \nabla A(\eta) &= \frac{\nabla Z(\eta)}{Z(\eta)} \\\\ &= \frac{1}{Z(\eta)} \int \mathcal{T}(x) h(x) \exp\{\eta^\top \mathcal{T}(x)\}\,dx \\\\ &= \int \mathcal{T}(x)\, p(x \mid \eta)\,dx = \mathbb{E}[\mathcal{T}(x)]. \end{align*} \] Differentiating a second time and applying the quotient rule to \(\nabla Z / Z\), \[ \begin{align*} \nabla^2 A(\eta) &= \frac{\nabla^2 Z(\eta)}{Z(\eta)} - \frac{\nabla Z(\eta)\,\nabla Z(\eta)^\top}{Z(\eta)^2} \\\\ &= \mathbb{E}[\mathcal{T}(x)\mathcal{T}(x)^\top] - \mathbb{E}[\mathcal{T}(x)]\,\mathbb{E}[\mathcal{T}(x)]^\top \\\\ &= \operatorname{Cov}[\mathcal{T}(x)]. \end{align*} \] The same computation applies verbatim to a discrete support, with the integral replaced by a sum over the support.

The first identity says the gradient of \(A\) returns the moment parameters \(m := \mathbb{E}[\mathcal{T}(x)]\). The second identity says the Hessian of \(A\) returns the covariance of the sufficient statistics.

As a consequence, for a minimal exponential family, the Hessian is strictly positive definite, and thus \(A(\eta)\) is strictly convex in \(\eta\). This convexity guarantees that the MLE (derived in the next section) is unique whenever it exists.

We now write two of the most widely used continuous distributions in machine learning in exponential family form.

Example 2: Normal Distribution

\[ \begin{align*} \mathcal{N}(x \mid \mu, \, \sigma^2) &= \frac{1}{\sigma\sqrt{2\pi}}\exp \left\{-\frac{1}{2\sigma^2}(x - \mu)^2 \right\} \\\\ &= \frac{1}{\sqrt{2\pi}} \exp \left\{ \frac{\mu}{\sigma^2}x - \frac{1}{2\sigma^2}x^2 -\frac{1}{2\sigma^2}\mu^2 -\log \sigma \right\} \end{align*} \] where

  • \(\mathcal{T}(x) = \begin{bmatrix}x \\ x^2 \end{bmatrix}\)
  • \(\eta = \begin{bmatrix} \frac{\mu}{\sigma^2} \\ -\frac{1}{2\sigma^2} \end{bmatrix} \)
  • \(A(\eta) = \frac{\mu^2}{2\sigma^2}+\log \sigma = -\frac{\eta_1^2}{4\eta_2}-\frac{1}{2}\log(-2\eta_2)\)
  • \(h(x) = \frac{1}{\sqrt{2\pi}}\).

Also, the moment parameters are given by: \[ m = \begin{bmatrix} \mu \\ \mu^2 + \sigma^2 \end{bmatrix}. \] Note. If \(\sigma = 1\), the distribution becomes a natural exponential family such that

  • \(\mathcal{T}(x) = x\)
  • \(\eta = \mu\)
  • \(A(\eta) = \frac{\mu^2}{2} = \frac{\eta^2}{2}\)
  • \(h(x) = \frac{1}{\sqrt{2\pi}}\exp\{-\frac{x^2}{2}\} = \mathcal{N}(x \mid 0, 1)\), which is not constant.

The univariate case extends naturally to higher dimensions. The multivariate normal distribution can be rewritten in the information form, parameterized by the precision matrix \(\Lambda = \Sigma^{-1}\) and \(\xi = \Lambda\mu\), which converts directly into canonical exponential-family form (Example 3 below). This form plays a central role in graphical models and message passing algorithms.

Example 3: Multivariate Normal Distribution (MVN)

\[ \begin{align*} \mathcal{N}(x \mid \mu, \Sigma) &= \frac{1}{(2\pi)^{\frac{D}{2}}\sqrt{\det(\Sigma)}} \exp \left\{ -\frac{1}{2}x^\top \Sigma^{-1}x + x^\top \Sigma^{-1}\mu -\frac{1}{2}\mu^\top \Sigma^{-1}\mu \right\}\\\\ &= c \exp\left \{x^\top \Sigma^{-1}\mu -\frac{1}{2}x^\top \Sigma^{-1}x \right\} \end{align*} \] where \[ c = \frac{\exp \left\{-\frac{1}{2}\mu^\top \Sigma^{-1} \mu\right\}}{(2\pi)^{\frac{D}{2}}\sqrt{\det(\Sigma)}} \] and \(\Sigma\) is a covariance matrix.

Now, we represent this model using canonical parameters. \[ \mathcal{N}_c (x \mid \xi, \Lambda) = c' \exp \left\{x^\top \xi - \frac{1}{2}x^\top \Lambda x \right\} \] where

  • \(\Lambda = \Sigma^{-1}\) is a precision matrix
  • \(\xi = \Sigma^{-1}\mu\) is a precision-weighted mean vector
  • \(c' = \frac{\exp \left\{-\frac{1}{2}\xi^\top \Lambda^{-1} \xi \right\}}{(2\pi)^{\frac{D}{2}}\sqrt{\det(\Lambda^{-1})}}\).

This representation is called information form and can be converted to exponential family notation as follows: \[ \begin{align*} \mathcal{N}_c (x \mid \xi, \Lambda) &= (2\pi)^{-\frac{D}{2}} \exp \left\{\frac{1}{2}\log | \Lambda | -\frac{1}{2}\xi^\top \Lambda^{-1}\xi \right\} \exp \left\{-\frac{1}{2}x^\top \Lambda x + x^\top \xi \right\} \\\\ &= h(x)g(\eta)\exp \left\{-\frac{1}{2}x^\top \Lambda x + x^\top \xi \right\} \\\\ &= h(x)g(\eta)\exp \left \{-\frac{1}{2}(\sum_{i, j}x_i x_j \Lambda_{ij}) + x^\top \xi \right\} \\\\ &= h(x)g(\eta)\exp \left\{-\frac{1}{2}\operatorname{vec}(\Lambda)^\top \operatorname{vec}(xx^\top) + x^\top \xi \right\} \\\\ &= h(x)\exp\{\eta^\top \mathcal{T}(x) - A(\eta)\} \end{align*} \] where

  • \(\mathcal{T}(x) = [x ; \operatorname{vec}(xx^\top)]\)
  • \(\eta = [\xi ; -\frac{1}{2}\operatorname{vec}(\Lambda)] = [\Sigma^{-1}\mu ; -\frac{1}{2}\operatorname{vec}(\Sigma^{-1})]\)
  • \(A(\eta) = -\log g(\eta) = -\frac{1}{2} \log | \Lambda | + \frac{1}{2}\xi^\top \Lambda^{-1} \xi \)
  • \(h(x) = (2\pi)^{-\frac{D}{2}}\).

The moment parameters are given by: \[ m = [\mu ; \operatorname{vec}(\mu\mu^\top + \Sigma)]. \]

Note. For \(D \geq 2\) this form is non-minimal: since \(x_i x_j = x_j x_i\), the entries of \(\operatorname{vec}(xx^\top)\) repeat, so the natural parameter \(\eta\) is over-specified by \(D(D-1)/2\) redundant components. A minimal representation would use only the upper-triangular (or lower-triangular) part of \(\Lambda\). However, in practice, the non-minimal representation is easier to plug into algorithms and stable for certain operations, while the minimal form is preferred for mathematical derivations.

MLE for the Exponential Family

One of the most useful consequences of the exponential family structure is that maximum likelihood estimation reduces to a simple condition: matching empirical moments to theoretical moments. This is because the gradient of the log-likelihood involves only the sufficient statistics and the log-partition function.

Assume \(\mathcal{D} = \{x_1, \ldots, x_N\}\) is an i.i.d. sample from an exponential family distribution with natural parameter \(\eta\). The likelihood is then given by: \[ \begin{align*} p(\mathcal{D} \mid \eta) &= \left\{\prod_{n=1}^N h(x_n)\right\} \exp \left\{\eta^\top \left[\sum_{n=1}^N \mathcal{T}(x_n)\right] - N A(\eta)\right\} \\\\ &\propto \exp\{\eta^\top \mathcal{T}(\mathcal{D}) - N A(\eta)\} \end{align*} \] where, dropping the \(\eta\)-independent factor \(\prod_n h(x_n)\), we define the aggregated sufficient statistic: \[ \mathcal{T}(\mathcal{D}) := \sum_{n=1}^N \mathcal{T}(x_n) \in \mathbb{R}^K. \]

The derivative of the log-partition function yields the expected value of the sufficient statistic vector: \[ \begin{align*} \nabla_{\eta} \log p(\mathcal{D} \mid \eta) &= \nabla_{\eta} \eta^\top \mathcal{T}(\mathcal{D}) - N \nabla_{\eta} A(\eta) \\\\ &= \mathcal{T}(\mathcal{D}) - N \mathbb{E}[\mathcal{T}(x)], \end{align*} \] where in the second line we used the first cumulant identity \(\nabla A(\eta) = \mathbb{E}[\mathcal{T}(x)]\) established in the previous section.

Setting this gradient to zero, we obtain the MLE \(\hat{\eta}\), characterized by \[ \mathbb{E}_{\hat{\eta}}[\mathcal{T}(x)] = \frac{1}{N}\sum_{n=1}^N \mathcal{T}(x_n). \] The left-hand side is the theoretical expectation of \(\mathcal{T}(x)\) under the model with natural parameter \(\hat{\eta}\). The right-hand side is the empirical average of the sufficient statistics. This principle is called moment matching.

Uniqueness does not give existence. A solution exists only when the empirical average on the right is an interior point of the set of attainable moments. For a Bernoulli sample in which every observation is \(0\), for instance, the condition forces \(\hat{\mu} = 0\), and no finite \(\hat{\eta}\) satisfies it.

For example, consider the univariate Gaussian \(\mathcal{N}(\mu, \sigma^2)\) with \(\mathcal{T}(x) = [x, x^2]^\top\) (Example 2 above). The moment matching condition becomes \[ \mathbb{E}_{\hat{\eta}}[x] = \bar{x}, \quad \mathbb{E}_{\hat{\eta}}[x^2] = \overline{x^2} := \frac{1}{N}\sum_{n=1}^N x_n^2. \]

Using \(\mathbb{E}[x] = \mu\) and \(\mathbb{E}[x^2] = \mu^2 + \sigma^2\), and provided the sample contains at least two distinct values (so that \(\overline{x^2} \gt \bar{x}^2\)), we obtain the closed-form MLE \[ \hat{\mu}_{\text{MLE}} = \bar{x}, \quad \widehat{\sigma^2}_{\text{MLE}} = \overline{x^2} - \bar{x}^2 = \frac{1}{N}\sum_{n=1}^N (x_n - \bar{x})^2. \]

In Variational Inference (VI), approximating an intractable posterior with an exponential family distribution is standard practice. Minimizing \(D_{\mathbb{KL}}(p \| q)\) over an exponential family \(q\) reproduces the moment matching condition of this section, with the expected sufficient statistics of \(p\) in place of the empirical average. Standard variational inference minimizes the reverse divergence \(D_{\mathbb{KL}}(q \| p)\), whose stationary condition is a different one. The log-partition function \(A(\eta)\) reappears as the Fenchel conjugate of the negative entropy, connecting the exponential family to convex duality and Information Geometry.

The exponential family structure reveals that the log-partition function \(A(\eta)\) encodes, through its derivatives on the interior of \(\Omega\), all the cumulants of the sufficient statistic \(\mathcal{T}(x)\). The support \(\mathcal{X}\) does not vary with \(\eta\), and for \(\eta\) in the interior of \(\Omega\), as in the proof of the cumulant identities above, we take for granted that differentiation may be carried out under the integral sign, so the expected negative Hessian form of the Fisher information matrix is available.

Differentiating \(\log p(x \mid \eta) = \log h(x) + \eta^\top \mathcal{T}(x) - A(\eta)\) twice in \(\eta\) leaves \(\nabla_\eta^2 \log p(x \mid \eta) = -\nabla^2 A(\eta)\), a quantity carrying no dependence on \(x\), so taking the expectation changes nothing and \(F(\eta) = \nabla^2 A(\eta)\). The second cumulant identity \(\nabla^2 A(\eta) = \operatorname{Cov}[\mathcal{T}(x)]\) established above is therefore a statement about the Fisher information itself, which measures how sensitive the distribution is to changes in the natural parameter and governs the geometry of the statistical model. That geometry is the subject of the next page.