Chapter 3Probabilistic learning

In this chapter, we review how to perform probabilistic learning. We also introduce exponential family distributions, as they play a key role in this book.

3.1 Probability distributions

3.1.1 Discrete probability distributions

A discrete probability distribution over a set $\cY$ is specified by its probability mass function (PMF) $p \colon \cY \to [0,1]$. The probability of $\y \in \cY$ is then defined by $$\PP(Y = \y) \coloneqq p(\y),$$ where $Y$ denotes a random variable. When $Y$ follows a distribution $p$, we write $Y \sim p$ (with some abuse of notation, we use the same letter $p$ to denote the distribution and the PMF). The expectation of $\phi(Y)$, where $Y \sim p$ and $\phi \colon \cY \to \RR^M$, is then $$\EE[\phi(Y)] = \sum_{\y \in \cY} p(\y) \phi(\y),$$ its variance (for one-dimensional variables) is $$\VV[\phi(Y)] = \EE[(\phi(Y) - \EE[\phi(Y)])^2] = \sum_{\y \in \cY} p(\y) (\phi(\y) - \EE[\phi(Y)])^2$$ and its mode is $$\argmax_{\y \in \cY} p(\y).$$ The Kullback-Leibler (KL) divergence (also known as relative entropy) between two discrete distributions over $\cY$, with associated PMFs $p$ and $q$, is the statistical “distance” defined by $$\mathrm{KL}(p, q) \coloneqq \sum_{\y \in \cY} p(\y) \log\left(\frac{p(\y)}{q(\y)}\right) = \EE_{Y \sim p} \log\left(\frac{p(Y)}{q(Y)}\right).$$

3.1.2 Continuous probability distributions

A continuous probability distribution over $\cY$ is specified by its probability density function (PDF) $p \colon \cY \to \RR_+$. The probability of $\cA \subseteq \cY$ is then $$\PP(Y \in \cA) = \int_\cA p(y) dy.$$ The expectation, variance and KL divergence are defined analogously to the discrete setting, simply replacing $\sum_{\y \in \cY}$ with $\int_{\cY}$. Specifically, the expectation of $\phi(Y)$ is $$\EE[\phi(Y)] = \int_{\cY} p(y) \phi(y) dy,$$ the variance is $$\VV[\phi(Y)] = \EE[(\phi(Y) - \EE[\phi(Y)])^2] = \int_{\cY} p(y) (\phi(y) - \EE[\phi(Y)])^2 dy$$ and the KL divergence is $$\mathrm{KL}(p, q) \coloneqq \int_\cY p(y) \log\left(\frac{p(y)}{q(y)}\right) dy = \EE_{Y \sim p} \log\left(\frac{p(Y)}{q(Y)}\right).$$ The mode is defined as the argmax of the PDF.

When $\cY = \RR$, we can also define the cumulative distribution function (CDF) $$\PP(Y \le b) = \int_{-\infty}^b p(y) dy.$$ The probability of $Y$ lying in the semi-closed interval $(a, b]$ is then $$\PP(a < Y \le b) = \PP(Y \le b) - \PP(Y \le a).$$

3.2 Maximum likelihood estimation

3.2.1 Negative log-likelihood

We saw that a probability distribution over $\cY$ is specified by $p(\y)$, which is called the probability mass function (PMF) for discrete variables or the probability density function (PDF) for continuous variables. In practice, the true distribution $p$ generating the data is unknown and we wish to approximate it with a distribution $p_\lambdav$, with parameters $\lambdav \in \Lambda$. Given a finite set of i.i.d. observations $\y_1, \dots, \y_N$, how do we fit $\lambdav \in \Lambda$ to the data? This can be done by maximizing the likelihood of the data, i.e., we seek to solve $$\widehat \lambdav_N \coloneqq \argmax_{\lambdav \in \Lambda} \prod_{i=1}^N p_{\lambdav}(\y_i).$$ This is known as maximum likelihood estimation (MLE). Because the $\log$ function is monotonically increasing, this is equivalent to minimizing the negative log-likelihood, i.e., we have $$\widehat \lambdav_N = \argmin_{\lambdav \in \Lambda} -\sum_{i=1}^N \log p_{\lambdav}(\y_i).$$

3.2.2 Consistency w.r.t. the Kullback-Leibler divergence

It is well-known that the MLE is consistent, in the sense of the Kullback-Leibler divergence. That is, denoting the true distribution $p$ and $$\lambdav_\infty \coloneqq \argmin_{\lambdav \in \Lambda} \mathrm{KL}(p, p_\lambdav)$$ where $$\mathrm{KL}(p, p_\lambdav) = \EE_{Y \sim p} \log\left(\frac{p(Y)}{p_\lambdav(Y)}\right),$$ then $\widehat \lambdav_N \to \lambdav_\infty$ in expectation over the observations, as $N \to \infty$. This can be seen by using $$\mathrm{KL}(p, p_\lambdav) \approx \frac{1}{N} \sum_{i=1}^N \log\left(\frac{p(\y_i)}{p_\lambdav(\y_i)}\right) = \frac{1}{N} \sum_{i=1}^N \left(\log p(\y_i) - \log p_\lambdav(\y_i)\right)$$ and the law of large numbers.

3.3 Probabilistic supervised learning

3.3.1 Conditional probability distributions

Many times in machine learning, instead of a probability $\PP(Y = \y)$ for some $\y \in \cY$, we wish to define a conditional probability $\PP(Y = \y | X = \x)$, for some input $\x \in \cX$. This can be achieved by reduction to an unconditional probability distribution, $$\PP(Y = \y | X = \x) \coloneqq p_\lambdav(\y)$$ where $$\lambdav \coloneqq f(\x, \w)$$ and $f$ is a model function with model parameters $\w \in \cW$. That is, rather than being a deterministic function from $\cX$ to $\cY$, $f$ is a function from $\cX$ to $\Lambda$, the set of permissible distribution parameters of the output distribution associated with the input.

We emphasize that $\lambdav$ could be a single parameter or a collection of parameters. For instance, in the Bernoulli distribution, $\lambda = \pi$, while in the univariate normal distribution, $\lambdav = (\mu, \sigma)$.

In Section 3.4 and throughout this book, we will also use the notation $p_\thetav$ instead of $p_\lambdav$ when $\thetav$ are the canonical parameters of an exponential family distribution.

3.3.2 Inference

The main advantage of this probabilistic approach is that our prediction model is much richer than if we just learned a function from $\cX$ to $\cY$. We now have access to the whole distribution over possible outcomes in $\cY$ and can compute various statistics:

• Probability: $\PP(Y = \y | X = \x)$ or $\PP(Y \in \cA | X=\x)$,

• Expectation: $\EE[\phi(Y)|X=\x]$ for some function $\phi$,

• Variance: $\VV[\phi(Y)|X=\x]$,

• Mode: $\argmax_{\y \in \cY} p_\lambdav(\y)$.

We now review probability distributions useful for binary classification, multiclass classification, regression, multivariate regression, and integer regression. In the following, to make the notation more lightweight, we omit the dependence on $\x$.

3.3.3 Binary classification

For binary outcomes, where $\cY = \{0,1\}$, we can use a Bernoulli distribution with parameter $$\lambda \coloneqq \pi \in [0,1].$$ When a random variable $Y$ is distributed according to a Bernoulli distribution with parameter $\pi$, we write $$Y \sim \mathrm{Bernoulli}(\pi).$$ The PMF of this distribution is $$p_\pi(y) \coloneqq \begin{cases} \pi &\mbox{if } y = 1 \\ 1 - \pi &\mbox{if } y = 0 \end{cases}.$$ The Bernoulli distribution is a binomial distribution with a single trial. Since $y \in \{0, 1\}$, the PMF can be rewritten as $$p_\pi(y) = \pi^y (1-\pi)^{1 -y}.$$ The mean is $$\EE[Y] = \pi = \PP(Y = 1)$$ and the variance is $$\VV[Y] = \pi (1-\pi) = \PP(Y=1)\PP(Y=0).$$

Parameterization using a sigmoid

Since the parameter $\pi$ of a Bernoulli distribution needs to belong to $[0,1]$, we typically use a sigmoid function (Section 4.4.3), such as a logistic function as the output layer: $$\pi \coloneqq f(\x, \w) \coloneqq \logistic(g(\x, \w)),$$ where $g \colon \cX \times \cW \to \RR$ is for example a neural network and $$\mathrm{logistic}(a) \coloneqq \frac{1}{1 + \exp(-a)} \in (0,1).$$ When $g$ is linear in $\w$, this is known as binary logistic regression.

3.3.4 Multiclass classification

For categorical outcomes with $M$ possible choices, where $\cY = [M]$, we can use a categorical distribution with parameters $$\lambdav \coloneqq \piv \in \triangle^M,$$ where we define the probability simplex $$\triangle^M \coloneqq \{\piv \in \RR_+^M \colon \langle \piv, \ones \rangle = 1 \},$$ i.e., the set of valid discrete probability distributions. When $Y$ follows a categorical distribution with parameter $\piv$, we write $$Y \sim \mathrm{Categorical}(\piv).$$ The PMF of the categorical distribution is $$p_\piv(y) \coloneqq \langle \piv, \phi(y) \rangle = \pi_y,$$ where $$\phi(y) \coloneqq \e_y$$ is the standard basis vector for the coordinate $y \in [M]$.

Since $Y$ is a categorical variable, it does not make sense to compute the expectation of $Y$ but we can compute that of $\phi(Y) = \e_Y$, $$\EE_{Y \sim p_\piv}[\phi(Y)] = \piv.$$ Therefore, as was also the case for the Bernoulli distribution, the mean and the probability distribution (represented by the vector $\piv$) are the same in this case.

Parameterization using a softargmax

Since the parameter vector $\piv$ of a categorical distribution needs to belong to $\triangle^M$, we typically use a softargmax as the output layer: $$\piv \coloneqq f(\x, \w) \coloneqq \softargmax(g(\x, \w)),$$ where $g \colon \cX \times \cW \to \RR^M$ is for example a neural network and $$\softargmax(\u) \coloneqq \frac{\exp(\u)}{\sum_j \exp(u_j)} \in \mathrm{relint}(\triangle^M).$$ The output of the softargmax is in the relative interior of $\triangle^M$, $\mathrm{relint}(\triangle^M) = \triangle^M \cap \RR^M_{>0}$. That is, the produced probabilities are always strictly positive. The categorical distribution is a multinomial distribution with a single trial. When $g$ is linear in $\w$, this is therefore known as multiclass or multinomial logistic regression, though strictly speaking a multinomial distribution could use more than one trial.

3.3.5 Regression

For real outcomes, where $\cY = \RR$, we can use, among other choices, a normal distribution with parameters $$\lambdav \coloneqq (\mu, \sigma),$$ where $\mu \in \RR$ is the mean parameter and $\sigma \in \RR_+$ is the standard deviation parameter. When $Y$ follows a normal distribution with parameters $(\mu, \sigma)$, we write $$Y \sim \mathrm{Normal}(\mu, \sigma).$$ The PDF is $$p_{\mu,\sigma}(y) \coloneqq \frac{1}{\sigma \sqrt{2\pi} } \exp\left(-\frac{1}{2}\frac{(y-\mu)^2}{\sigma^2}\right).$$ The expectation is $$\EE_{Y \sim p_{\mu,\sigma}}[Y] = \mu.$$ One advantage of the probabilistic perspective is that we are not limited to predicting the mean. We can also compute the CDF $$\PP(Y \le y) = \frac{1}{2}\left[ 1 + \mathrm{erf}\left(\frac{y - \mu}{\sigma \sqrt{2}}\right) \right],$$ where we used the error function $$\operatorname{erf}(z) \coloneqq \frac{2}{\sqrt\pi}\int_0^z e^{-t^2}\,\mathrm{d}t.$$ This function is available in most scientific computing libraries, such as SciPy (Virtanen et al., 2020). We can also write $$\PP(Y \le y) = \Phi\left(\frac{y-\mu}{\sigma}\right)$$ where $$\Phi(z) \coloneqq \frac{1}{2}\left[ 1 + \mathrm{erf}\left(\frac{z}{\sqrt{2}}\right) \right] \tag{3.1}$$ is the CDF of the standard Gaussian distribution (with zero mean and unit variance). From the CDF, we also easily obtain $$\PP(a < Y \le b) = \frac{1}{2} \left[ \mathrm{erf}\left(\frac{b - \mu}{\sigma \sqrt{2}}\right) - \mathrm{erf}\left(\frac{a - \mu}{\sigma \sqrt{2}}\right) \right].$$

Parameterization

Typically, in regression, the mean is output by a model, while the standard deviation $\sigma$ is kept fixed (typically set to $1$). Since $\mu$ is unconstrained, we can simply set $$\mu \coloneqq f(\x, \w) \in \RR,$$ where $f \colon \cX \times \cW \to \RR$ is for example a neural network. That is, the output of $f$ is the mean of the distribution, $$\EE_{Y \sim p_{\mu,1}}[Y] = \mu = f(\x, \w).$$ We can also use $\mu$ to predict $\PP(Y \le y)$ or $\PP(a < Y \le b)$, as shown above.

3.3.6 Multivariate regression

More generally, for multivariate outcomes, where $\cY = \RR^M$, we can use a multivariate normal distribution with parameters $$\lambdav \coloneqq (\muv, \Sigmav),$$ where $\muv \in \RR^M$ is the mean and $\Sigmav \in \RR^{M \times M}$ is the covariance matrix. When $Y$ follows a multivariate normal distribution with parameters $(\muv, \Sigmav)$, we write $$Y \sim \mathrm{Normal}(\muv, \Sigmav).$$ The PDF is $$p_{\muv,\Sigmav}(\y) \coloneqq \frac{1}{\sqrt{(2 \pi)^M |\Sigmav|}} \exp\left(-\frac{1}{2} \langle \y - \muv, \Sigmav^{-1} (\y - \muv) \rangle\right).$$ Using a diagonal covariance matrix is equivalent to using $M$ independent normal distributions for each $Y_j$, for $j \in [M]$. The expectation is $$\EE_{Y \sim p_{\muv,\Sigmav}}[Y] = \muv.$$

Parameterization

Typically, in multivariate regression, the mean is output by a model, while the covariance matrix is kept fixed (typically set to the identity matrix). Since $\muv$ is again unconstrained, we can simply set $$\muv \coloneqq f(\x, \w) \in \RR^M.$$ More generally, we can parametrize the function $f$ so as to output both the mean $\muv$ and the covariance matrix $\Sigmav$, i.e., $$(\muv, \Sigmav) \coloneqq f(\x, \w) \in \RR^M \times \RR^{M \times M}.$$ The function $f$ must be designed such that $\Sigmav$ is symmetric and positive semi-definite. This is easy to achieve for instance by parametrizing $\Sigmav = \S \S^\top$ for some matrix $\S$.

3.3.7 Integer regression

For integer outcomes, where $\cY = \NN$, we can use, among other choices, a Poisson distribution with mean parameter $\lambda > 0$. When $Y$ follows a Poisson distribution with parameter $\lambda$, we write $$Y \sim \mathrm{Poisson}(\lambda).$$ The PMF is $$\PP(Y = y) = p_\lambda(y) \coloneqq \frac{\lambda^y \exp(-\lambda)}{y!}.$$ It is the probability of $y$ events occurring in an interval of time. The Poisson distribution is frequently used when there is a large number of possible events, each of which is rare.

The CDF is $$\PP(Y \le y) = \sum_{i=0}^y \PP(Y = i).$$ The Poisson distribution implies that the index of dispersion (the ratio between variance and mean) is $1$, since $$\EE[Y] = \VV[Y] = \lambda.$$ When this assumption is inappropriate, one can use generalized Poisson distributions (Satterthwaite, 1942).

Parameterization using an exponential

Since the parameter $\lambda$ of a Poisson distribution needs to be strictly positive, we typically use an exponential function as output layer: $$\lambda \coloneqq f(\x, \w) \coloneqq \exp(g(\x, \w)) > 0,$$ where $g \colon \cX \times \cW \to \RR$.

3.3.8 Loss functions

We now discuss how to learn the model parameters $\w \in \cW$ from input-output pairs $(\x_1, \y_1), \dots, (\x_N, \y_N)$.

Deterministic vs. probabilistic approaches

In a deterministic approach, if we used a mapping $f \colon \cX \times \cW \to \cY$, we could formulate an objective function of the form $$L(\w) \coloneqq \frac{1}{N} \sum_{i=1}^N \ell(f(\x_i, \w), \y_i),$$ where $\ell \colon \cY \times \cY \to \RR$ is a loss function. Unfortunately, $L(\w)$ would be typically discontinuous if $\cY$ is a discrete output space (as is the case in classification), making optimization difficult.

In contrast, in the probabilistic approach, we use a mapping $f \colon \cX \times \cW \to \Lambda$ to distribution parameters, and we can formulate an objective of the form $$L(\w) \coloneqq \frac{1}{N} \sum_{i=1}^N \ell(f(\x_i, \w), \y_i),$$ where $\ell \colon \Lambda \times \cY \to \RR$. This is typically a continuous objective, since $\Lambda$ is typically a continuous set and $p_\lambdav$ varies continuously w.r.t. $\lambdav$ even if $p_\lambdav$ is a distribution over a discrete set $\cY$. In other words, the probabilistic approach is not only powerful for the inference it allows us to do (probability, expectation, variance, mode), but also because it allows us to formulate a continuous and typically differentiable objective function!

This is summarized in Fig. 3.5.

Negative log-likelihood

In the conditional setting briefly reviewed in Section 3.3.1, we can use maximum likelihood estimation (MLE) to estimate the model parameters $\w \in \cW$ of $f$. Given a set of input-output pairs $(\x_1, \y_1), \dots, (\x_N, \y_N)$, we choose the model parameters that maximize the likelihood of the data, $$\widehat \w_N \coloneqq \argmax_{\w \in \cW} \prod_{i=1}^N p_{\lambdav_i}(\y_i),$$ where $$\lambdav_i \coloneqq f(\x_i, \w).$$ Again, this is equivalent to minimizing the negative log-likelihood, $$\widehat \w_N = \argmin_{\w \in \cW} -\sum_{i=1}^N \log p_{\lambdav_i}(\y_i).$$ In the notation above, this corresponds to defining the loss function $$\ell(\lambdav_i, \y_i) \coloneqq -\log p_{\lambdav_i}(\y_i).$$

Recovering well-known loss functions

Interestingly, MLE allows us to recover several popular loss functions.

• For the Bernoulli distribution with parameter $\lambda_i = \pi_i = \mathrm{logistic}(\theta_i)$, where $\theta_i \coloneqq g(\x_i, \w)$, we have $$\begin{align*} -\log p_{\lambda_i}(y_i) &= -\left[y_i \log \pi_i + (1-y_i) \log (1-\pi_i)\right] \\ &= \log(1 + \exp(\theta_i)) - y_i \theta_i \\ &= \mathrm{softplus}(\theta_i) - y_i \theta_i, \end{align*}$$ which is the binary logistic loss function.

• For the categorical distribution with parameters $\lambdav_i = \piv_i = \mathrm{softargmax}(\thetav_i)$, where $\thetav_i \coloneqq g(\x_i, \w)$, we have $$\begin{align*} -\log p_{\lambdav_i}(y_i) &= -\log \pi_{i,y_i} \\ &= \log \sum_{j=1}^M \exp(\theta_{i,j}) - \theta_{i,y_i} \\ &= \mathrm{logsumexp}(\thetav_i) - \langle \thetav_i, \e_{y_i} \rangle, \end{align*}$$ which is the multiclass logistic loss function, also known as cross-entropy loss.

• For the normal distribution with mean $\lambda_i = \mu_i = f(\x_i, \w)$ and fixed variance $\sigma_i^2$, we have $$-\log p_{\lambda_i}(y_i) = \frac{1}{2\sigma_i^2} (y_i - \mu_i)^2 +\frac{1}{2} \log \sigma_i^2 + \frac{1}{2} \log(2\pi),$$ which is, up to a constant and with unit variance, the squared loss function.

• For the Poisson distribution with mean $\lambda_i = \exp(\theta_i)$, where $\theta_i \coloneqq g(\x_i, \w)$, we have $$\begin{align*} -\log p_{\lambda_i}(y_i) &= -y_i \log(\lambda_i) + \lambda_i + \log(y_i!) \\ &= -y_i \theta_i + \exp(\theta_i) + \log(y_i!) \end{align*}$$ which is the Poisson loss function. The loss function is convex w.r.t. $\lambda_i$ and $\theta_i$ for $y_i \ge 0$.

3.4 Exponential family distributions

3.4.1 Definition

The exponential family is a class of probability distributions, whose PMF or PDF can be written in the form $$\begin{align*} p_\thetav(\y) &= \frac{h(\y) \exp\left[\langle \thetav, \phi(\y) \rangle\right]}{\exp(A(\thetav))} \\ &= h(\y) \exp\left[\langle \thetav, \phi(\y) \rangle - A(\thetav)\right], \end{align*}$$ where $\thetav$ are the natural or canonical parameters of the distribution. The function $h$ is known as the base measure. The function $\phi$ is the sufficient statistic: it holds all the information about $\y$ and is used to embed $\y$ in a vector space. The function $A$ is the log-partition or log-normalizer (see below for details). All the distributions we reviewed in Section 3.3 belong to the exponential family. With some abuse of notation, we use $p_\lambdav$ for the distribution in original form and $p_\thetav$ for the distribution in exponential family form. As we will see, we can go from $\thetav$ to $\lambdav$ and vice versa. We illustrate how to rewrite a distribution in exponential family form below.

We rewrite the previously-described distributions in exponential family form in Table 3.1. This list is non-exhaustive: there are many more distributions in the exponential family! (Barndorff-Nielsen, 2014)

3.4.2 The log-partition function

We focus our discussion in this section on discrete variables for simplicity. However, the results stated in this section also hold for continuous variables (Wainwright and Jordan, 2008). The log-partition function $A$ is the logarithm of the distribution’s normalization factor, $$A(\thetav) \coloneqq \log \sum_{\y \in \cY} h(\y) \exp\left[\langle \thetav, \phi(\y) \rangle\right].$$ We denote the set of valid parameters $$\Theta \coloneqq \{ \thetav \in \RR^M \colon A(\thetav) < +\infty\} \subseteq \RR^M.$$ We can conveniently rewrite $A(\thetav)$ as $$A(\thetav) = \mathrm{logsumexp}(B(\thetav)) \coloneqq \log \sum_{\y \in \cY} \exp([B(\thetav)]_\y),$$ where we defined the affine map $$B(\thetav) \coloneqq (\langle \thetav, \phi(\y) \rangle + \log h(\y))_{\y \in \cY}.$$ Since $A(\thetav)$ is the composition of $\mathrm{logsumexp}$, a convex function, and of $B$, an affine map, we immediately obtain the following proposition.

A major property of the log-partition function is that its gradient coincides with the expectation of $\phi(Y)$ according to $p_\thetav$.

The gradient $\nabla A(\thetav)$ is therefore often called the mean function. The set of achievable means $\muv(\thetav)$ is defined by $$\cM \coloneqq \mathrm{conv}(\phi(\cY)) \coloneqq \{ \EE_{p}[\phi(Y)] \colon p \in \cP(\cY)\},$$ where $\mathrm{conv}(\cS)$ is the convex hull of $\cS$ and $\cP(\cY)$ is the set of valid probability distributions over $\cY$.

Similarly, the Hessian $\nabla^2 A(\thetav)$ coincides with the covariance matrix of $\phi(Y)$ according to $p_\thetav$ (Wainwright and Jordan, 2008, 3).

When the exponential family is minimal, which means that the parameters $\thetav$ uniquely identify the distribution, it is known that $\nabla A$ is a one-to-one mapping from $\Theta$ to $\cM$. That is, $\muv(\thetav) = \nabla A(\thetav)$ and $\thetav = (\nabla A)^{-1}(\muv(\thetav))$.

3.4.3 Maximum entropy principle

Suppose we observe the empirical mean $\widehat \muv = \frac{1}{N} \sum_{i=1}^N \phi(\y_i) \in \cM$ of some observations $\y_1,\dots,\y_N$. How do we find a probability distribution achieving this mean? Clearly, such a distribution may not be unique. One way to choose among all possible distributions is by using the maximum entropy principle. Let us define the Shannon entropy by $$H(p) \coloneqq -\sum_{\y \in \cY} p(\y) \log p(\y)$$ for discrete variables and by $$H(p) \coloneqq -\int_{\cY} p(y) \log p(y) dy$$ for continuous variables. This captures the level of “uncertainty” in $p$, i.e., it is maximized when the distribution is uniform. Then, the maximum entropy distribution satisfying the first-order moment condition (i.e., whose expectation matches the empirical mean) is $$p^\star \coloneqq \argmax_{p \in \cP(\cY)} H(p) \quad \text{s.t.} \quad \EE_{Y \sim p}[\phi(Y)] = \widehat \muv.$$ It can be shown that the maximum entropy distribution is necessarily in the exponential family with sufficient statistic defined by $\phi$ and its canonical parameters $\thetav$ coincide with the Lagrange multipliers of the above constraint (Wainwright and Jordan, 2008, Section 3.1).

3.4.4 Maximum likelihood estimation

As in Section 3.2, to fit the parameters $\thetav \in \Theta$ of an exponential family distribution to some i.i.d. observations $\y_1, \dots, \y_N$, we can use the MLE principle, i.e., $$\widehat \thetav_N = \argmax_{\thetav \in \Theta} \prod_{i=1}^N p_\thetav(\y_i) = \argmin_{\thetav \in \Theta} -\sum_{i=1}^N \log p_\thetav(\y_i).$$ Fortunately, for exponential family distributions, the log probability/density enjoys a particularly simple form.

It follows from Proposition 3.1 that $\thetav \mapsto -\log p_\thetav(\y)$ is convex. Interestingly, we see that the gradient is the residual between the expectation of $\phi(Y)$ according to the model and the observed $\phi(\y)$. Therefore, the negative log-likelihood of an exponential family distribution can be seen as performing first moment matching.

3.4.5 Probabilistic learning with exponential families

In the supervised probabilistic learning setting, we wish to estimate a conditional distribution of the form $p_\w(\y \mid \x)$. Given a model function $f$, such as a neural network, a common approach for defining such a conditional distribution is by reduction to the unconditional setting, $$p_\w(\y \mid \x) \coloneqq p_\thetav(\y) \quad \text{where} \quad \thetav \coloneqq f(\x, \w).$$ In other words, the role of $f$ is to produce the parameters of $p_\thetav$ given some input $\x$. It is a function from $\cX \times \cW$ to $\Theta$. Note that $f$ must be designed such that it produces an output in $$\Theta \coloneqq \{ \thetav \in \RR^M \colon A(\thetav) < +\infty\}.$$ Many times, $\Theta$ will be the entire $\RR^M$ but this is not always the case. For instance, as we previously discussed, for a multivariate normal distribution, where $\thetav = (\muv, \Sigmav) = f(\x, \w)$, we need to ensure that $\Sigmav$ is a positive semidefinite matrix.

Training

Given input-output pairs $\{(\x_i, \y_i)\}_{i=1}^N$, we then seek to find the parameters $\w$ of $f(\x, \w)$ by minimizing the negative log-likelihood $$\argmin_{\w \in \cW} -\sum_{i=1}^N \log p_{\thetav_i}(\y_i) = \argmin_{\w \in \cW} \sum_{i=1}^N A(\thetav_i) - \langle \thetav_i, \phi(\y_i) \rangle$$ where $\thetav_i \coloneqq f(\x_i, \w)$. While $-\log p_\thetav(\y)$ is a convex function of $\thetav$ for exponential family distributions, we emphasize that $-\log p_{f(\x,\w)}(\y)$ is typically a nonconvex function of $\w$, when $f$ is a nonlinear function, such as a neural network.

Inference

Once we have found $\w$ by minimizing the objective function above, there are several possible strategies to perform inference for a new input $\x$.

• Expectation. When the goal is to compute the expectation of $\phi(Y)$, we can use $\nabla A(f(\x, \w))$. That is, we compute the distribution parameters associated with $\x$ by $\thetav = f(\x, \w)$ and then we compute the mean by $\muv = \nabla A(\thetav)$. When $f$ is linear in $\w$, the composition $\nabla A \circ f$ is called a generalized linear model. This is illustrated in Fig. 3.6.

• Probability. When the goal is to compute the probability of a certain $\y$, we can compute the distribution parameters associated with $\x$ by $\thetav = f(\x, \w)$ and then we can compute $\PP(Y = \y | X = \x) = p_\thetav(\y)$. In the particular case of the categorical distribution (of which the Bernoulli distribution is a special case), we point out again that the mean and the probability vector coincide: $$\muv = \piv = \nabla A(\thetav) = \softargmax(\thetav) \in \triangle^M.$$

• Other statistics. When the goal is to compute other quantities, such as the variance or the CDF, we can convert the natural parameters $\thetav$ to the original distribution parameters $\lambdav$ (see Table 3.1 for examples). Then, we can use established formulas for the distribution in original form, to compute the desired quantities.