One of the very first models studied in any introductory machine learning course is the logistic regression model for binary classification:

\[\hat y_i = \sigma(X_i \beta)\]

Here, \(X_i = [ x_{i1} x_{i2} \dots x_{id} ]\) are the features (also called covariates) of an input we want to classify. For example, these could be the pixels of an image, or a vector representation of some text. We’ll also assume \(x_{i1} = 1\), so as to include a bias term in our model. \(\beta\) is a \(d\)-dimensional column vector of our learnable model parameters, and \(\hat y_i\) is the model output. The function \(\sigma\) is the sigmoid function, and it is given by:

\[\sigma(z) = \frac{1}{1 + \exp(-z)} \tag{1}\label{eq1}\]

In an introductory ML course, the sigmoid function is produced seemingly out of nowhere, and it happens to be exactly what we need to make our model work:

  • The range of \(\sigma\) is \((0, 1)\), which means our model output \(\hat y_i\) can be interpreted as a probability distribution (more on this later). Hence our model can make classification predictions simply by rounding \(\hat y_i\) to the nearest integer, either 0 or 1.
  • \(\sigma\) is differentiable everywhere, and its derivative can be elegantly expressed as \(\sigma'(z) = \sigma(z)(1 - \sigma(z))\). In particular the derivative is nowhere 0, because \(\sigma(z)\) is within \((0,1)\) for any \(z\). This is a useful property if we are training our model using gradient descent.

But certainly sigmoid isn’t the only differentiable function to possess these desirable properties. What makes sigmoid so well-suited for binary classification? Where does the expression for sigmoid \(\eqref{eq1}\) come from? In the following sections I will present an answer to these questions from a statistical perspective, and explain why sigmoid is in fact the most “natural” choice for our binary classification model.

The Exponential Family

We begin with an overview of the exponential family of distributions, perhaps a seemingly unrelated topic. A probability distribution with density \(p(x \vert \theta)\) is said to belong to the exponential family if \(p(x \vert \theta)\) can be written in the following form:

\[p(x \vert \theta) = h(x) \exp \Big( \eta(\theta)^\top T(x) - A(\eta(\theta))\Big) \tag{2}\label{eq2}\]
  • \(\theta\) could be any set of parameters, and \(\eta\) is a reparametrization function. The parameters \(\eta(\theta)\) are called the natural parameters (or canonical parameters) of the distribution.
  • \(T(x)\) is called the sufficient statistic for the distribution (we’ll see why).
  • \(A(\eta(\theta))\) is called the cumulant function. It can be interpreted as a normalization constant for the PDF.
  • \(h(x)\) could be any nonnegative function.

It turns out that many common distributions belong to the exponential family, including the Bernoulli, Gaussian, Poisson, Gamma, and Beta distributions.

Example: Bernoulli Distribution

As an example, let’s see how the PDF of the Bernoulli distribution can be expressed in the form of \(\eqref{eq2}\). A common expression for the PDF of the Bernoulli is:

\[p(x \vert \theta) = \theta^x (1 - \theta)^{1-x}\]

where \(\theta\) is a parameter in \([0, 1]\). This gives:

\[p(x \vert \theta) = \Big(\frac{\theta}{1-\theta}\Big)^x (1 - \theta)\] \[= \exp \Big(x \log\Big(\frac{\theta}{1-\theta}\Big) + \log(1 - \theta) \Big) \tag{3}\label{eq3}\]

Note that \(\eqref{eq3}\) is in the desired form \(\eqref{eq2}\):

  • \(\eta(\theta) = \log\big(\frac{\theta}{1-\theta}\big)\) is the natural parameter
  • \(T(x) = x\) is the sufficient statistic
  • \(A(\eta(\theta)) = - \log (1 - \theta)\) is a normalizing constant
  • $h(x) = 1$

In this example, \(\theta\) is called the mean parameter of the Bernoulli distribution. In general, distributions in the exponential family have mean parameters which are not necessarily equal to their natural parameters. Formally, the mean parameters of a distribution in the exponential family are defined as \(\theta := \mathbb{E}_{x \sim p}[T(x)]\) , the expected value of the sufficient statistic. Note that since \(T(x) = x\) for the Bernoulli distribution, this definition coincides with our intuition for the mean parameter of a Bernoulli distribution being the expected success rate: \(\theta = \mathbb{E}_{x \sim p}[x]\). As another example, the mean parameters of the Gaussian distribution are the familiar \(\theta = [ \mu \sigma^2 ]\), while its natural parameters are the (perhaps not-so-natural) \([\frac{\mu }{\sigma ^{2}}{\frac {1}{2\sigma ^{2}}}]\).

Sufficient Statistics

Why should we care about the exponential family? One reason is that their sufficient statistics are easy to determine, which makes them valuable for parametric inference. Recall that a statistic is any function of random samples \(X\) from a distribution. Informally, a statistic \(T(X)\) is sufficient for a distribution \(p(x \vert \theta)\) if we can estimate parameters \(\theta\) using \(T(X)\), and if there is no further information we can obtain from \(X\) about \(\theta\) that is not already included in \(T(X)\).

As an example, consider the Bernoulli distribution once again, parametrized by \(\theta\). If we sample \(x_1,...x_n \stackrel{iid}{\sim} Bernoulli(\theta)\), then knowing the sample mean \(\frac{1}{n} \sum_{i=1}^n x_i\) is enough to estimate \(\theta\). Indeed, the MLE estimate for \(\theta\) is \(\hat \theta = \frac{1}{n} \sum_{i=1}^n x_i\). Knowing anything more about \(x_1,...,x_n\) could not give us a better estimate for \(\theta\). We say that \(\sum_{i=1}^n x_i\) is a sufficient statistic for the Bernoulli distribution (assuming \(n\) is known beforehand).

Notice that the joint distribution of samples \(x_1,...,x_n\) drawn \(iid\) from an arbitrary distribution in the exponential family is:

\[p(x_1,...,x_n \vert \theta) = \Big(\prod_{i=1}^n h(x_i) \Big) \exp\Big( \eta(\theta)^\top \underbrace{\big(\sum_{i=1}^n T(x_i)\big)}_{\text{sufficient!}} - n A(\eta(\theta))\Big) \tag{4}\label{eq4}\]

In the previous section we showed \(T(x) = x\) for the Bernoulli distribution. It follows that \(\sum_{i=1}^n T(x_i) = \sum_{i=1}^n x_i\) is precisely the sufficient statistic for the Bernoulli distribution! It can be shown that in the general case as well, \(\sum_{i=1}^n T(x_i)\) is a sufficient statistic for any distribution in the exponential family.

Hence, all the information needed to estimate the parameters of an exponential-family distribution from \(iid\) samples \(x_1,...,x_n\) is contained in a single vector \(\sum_{i=1}^n T(x_i)\) of fixed dimension (as opposed to dimension that grows with sample size \(n\)), and this vector can simply be read off the PDF \(\eqref{eq4}\) of the joint distribution. As an aside, it turns out that much more is true. The Pitman-Koopman-Darmois theorem states that the exponential family of distributions is the only family of distributions (under a few mild assumptions) for which the sufficient statistics are finite-dimensional for arbitrarily large sample sizes.

Reparametrizing

Earlier, we computed a reparametrization \(\eta(\theta) = \log\big(\frac{\theta}{1-\theta}\big)\) of the Bernoulli PDF which mapped its mean parameter \(\theta\) to its natural parameter. Notice furthermore that the mapping is invertible, with inverse given by (does this look familiar?):

\[\theta = \frac{1}{1 + \exp(-\eta(\theta))}\]

It is now natural to ask whether it is always possible to reparametrize from the natural parameters \(\eta\) to the mean parameters \(\theta = \mathbb{E}_{x \sim p}[T(x)]\) for any distribution in the exponential family. If \(\theta\) are the mean parameters, can we always find an \(\eta(\theta)\) that is one-to-one? As it turns out, we can. Yet another beautiful property satisfied by any distribution in the exponential family is the following:

\[\frac{\partial A(\eta)}{\partial \eta} = \mathbb{E}_{x \sim p}[T(x)] \tag{5}\label{eq5}\] \[\frac{\partial^2 A(\eta)}{\partial \eta ^2} = \text{Var}_{x \sim p}[T(x)] \tag{6}\label{eq6}\]

Here, \(\eta := \eta(\theta)\) denotes the natural parameters and \(A(\eta)\) is the cumulant function. Note that \(\eqref{eq5}\) is precisely the expression for the mean parameters \(\theta\), so we can write \(\eqref{eq6}\) as:

\[\frac{\partial^2 A(\eta)}{\partial \eta ^2} = \frac{\partial \theta(\eta) }{\partial \eta} = \text{Var}_{x \sim p}[T(x)]\]

Since the variance \(\text{Var}_{x \sim p}[T(x)]\) is a strictly positive quantity, it follows that \(\dfrac{\partial \theta(\eta)}{\partial \eta} > 0\).

Therefore, we can always express the mean parameters \(\theta\) as a function of \(\eta\) (the function given by \(\eqref{eq5}\)), and this function is guaranteed to be strictly increasing, hence invertible.

The Generalized Linear Model

Classification, Revisited

In the introduction, we considered the logistic regression model from a discriminative perspective: our model was a function of input features, parametrized by \(\beta\), that would make predictions to distinguish (“discriminate”) between inputs from different classes. Suppose we now consider the same model, from a generative perspective. We assume that conditioned on the input features \(X_i\), the true class \(y_i\) of an input was drawn by some generative process from a Bernoulli distribution with unknown mean parameter \(\theta_i\). We have observed \(X_1,...,X_n\) and \(y_1,...,y_n\). Our goal is to learn \(\beta\) so we can predict a likely estimate for each \(\theta_i\). Here is this model:

\[y_i \sim Bernoulli(\theta_i)\] \[\theta_i = \sigma(X_i\beta) \tag{7}\label{eq7}\]

Generalized Linear Model

As we now know quite well, the Bernoulli distribution belongs to the more general exponential family. We have seen how numerous properties of the Bernoulli distribution could be generalized to arbitrary distributions in the exponential family, including the formula for a sufficient statistic, and the invertibility of \(\eta(\theta)\) when \(\theta\) are the mean parameters. Can we generalize our classification model as well?

In \(\eqref{eq7}\) we are modeling the mean parameter of the Bernoulli distribution. Instead, let’s choose an arbitrary distribution \(p(x \vert \theta)\) in the exponential family, parametrized by its mean parameters \(\theta\), and model these. For simplicity we assume the mean parameter \(\theta\) is one-dimensional, as in the Bernoulli distribution. We also replace \(\sigma\) in \(\eqref{eq7}\) with an arbitrary invertible function \(g^{-1}\). Here is the resulting model:

\[y \sim p(\cdot \vert \theta)\] \[g(\theta) = X\beta\]

We have just constructed the Generalized Linear Model! The Generalized Linear Model (GLM) is a model for the mean parameters of an arbitrary distribution in the exponential family. The model is specified by three components:

  • The distribution \(p(\cdot \vert \theta)\): a distribution from the exponential family parametrized by its mean parameters \(\theta\)
  • The link function \(g(\theta)\): an invertible function of the mean parameters \(\theta\)
  • Covariates \(X\) and a linear predictor \(X\beta\), where \(\beta\) are the model parameters

Notice that our model “works” precisely because of the elegant properties of the exponential family described earlier. We always know what a sufficient statistic looks like for a distribution in the exponential family, and in particular we know that it will always be of a fixed dimension. Therefore it is feasible to model the mean parameters \(\theta\), which correspond to the expectation of the sufficient statistic. We also know that there is an invertible parametrization \(\theta(\eta)\) of the mean parameters, so we will be able to recover an estimate for the natural parameters \(\eta\) if we can estimate \(\theta\).

The power of the GLM is in its generality. While the logistic regression model could only model binary data, the GLM provides a framework to model response data drawn from any distribution in the exponential family (and by the Pitman-Koopman-Darmois theorem, drawn from any distribution with sufficient statistics of fixed dimension). We can now model count data with the Poisson distribution, wait times with Gamma, multi-class data with the Multinomial. In addition, the GLM is consistent with other well-known models: with the Gaussian distribution, for instance, we can recover traditional linear regression.

On the flipside, the generality of the GLM may also seem daunting. We now have so much more freedom in specifying our model. For example, how do we choose a link function \(g\)? In general there are several possible options for the link function for a given distribution, and each has their own advantages. In the next section we will examine a special choice for the link function.

The canonical link is the link function that results in the natural parameters \(\eta\) themselves being modeled as a linear function \(X\beta\) of the covariates. To be precise, the canonical link function is simply \(g(\theta) = \eta(\theta)\).

This is a cute definition, but it might seem like another trick pulled out of a hat. What makes this choice of link function so special? We will prove the following result:

Proposition: Suppose we have a GLM where the canonical link function \(g\) is used, and the specified distribution satisfies \(T(y) = y\). Then the MLE estimate \(\hat \beta\) satisfies the orthogonality condition \(\sum_{i=1}^n (y_i - \hat \theta_i)X_i = 0\), where \(\hat \theta_i = g^{-1}(X_i \hat \beta)\).

Proof. Given \(iid\) samples \((X_1, y_1), ..., (X_n, y_n)\) where the \(y_i\)’s are distributed according to the specified exponential-family distribution, the log-likelihood for this data is:

\[\sum_{i=1}^n \log(p(y_i \vert \hat \theta_i)) = \sum_{i=1}^n \log (h(y_i)) + \sum_{i=1}^n \eta(\hat \theta_i) T(y_i) - A(\eta(\hat \theta_i))\] \[= \sum_{i=1}^n \log (h(y_i)) + \sum_{i=1}^n \eta(\hat \theta_i) y_i - A(\eta(\hat \theta_i)) \qquad \text{since } T(y) = y\]

The canonical link function is \(g(\theta) = \eta(\theta)\), and under our model we have \(g(\hat \theta_i) = X_i \hat \beta\). Plugging this in gives:

\[\sum_{i=1}^n \log(p(y_i \vert \hat \theta_i)) = \sum_{i=1}^n \log (h(y_i)) + \sum_{i=1}^n (X_i \hat \beta) y_i - A(X_i \hat \beta) \tag{8}\label{eq8}\]

Now we differentiate \(\eqref{eq8}\) with respect to \(\beta\) and set the derivative to 0 (using \(\eqref{eq5}\)), since this condition must be satisfied by the MLE \(\hat \beta\):

\[\sum_{i=1}^n X_iy_i - X_i \hat \theta_i = 0\] \[\sum_{i=1}^n (y_i - \hat \theta_i)X_i = 0 \tag{9}\label{eq9}\]

This is the orthogonality condition, as required. \(\blacksquare\)

Note that with the assumption \(T(y) = y\), the mean parameter \(\theta_i\) is precisely the expected value of \(y_i\). Let \(\epsilon_i := y_i - \hat \theta_i\). This corresponds to the leftover noise in our observations \(y_i\) that has not been explained by our model of their means \(\theta_i\). Now recall that the feature vectors \(X_i = [x_{i1} x_{i2} \dots x_{id}]\) have \(x_{i1} = 1\), so as to include a bias term in our linear predictor. Examining the vector equality \(\eqref{eq9}\) along the first component, it follows that \(\bar \epsilon := \sum_{i=1}^n \epsilon_i = 0\). Hence, with \(\bar x_{\cdot j} := \sum_{i=1}^n x_{ij}\), for each covariate \(x_{\cdot j}\) we have:

\[\text{Cov}(\epsilon_i, x_{ij}) = \frac{1}{n-1} \Big(\sum_{i=1}^n \epsilon_i x_{ij} - \bar \epsilon \bar x_{\cdot j}\Big) = 0\]

(Here \(\text{Cov}(\cdot)\) denotes the sample covariance.)

The unexplained noise \(\epsilon_i\) and each covariate \(x_{ij}\) are uncorrelated! In other words, all the linear dependence of \(\theta_i\) on \(X_i\) has been “explained away” by our model under the assumptions we made. This is intuitively a very natural and desirable property for our optimal model parameters to satisfy. As a final note, the assumption \(T(y) = y\) that we needed for this proof holds for many common distributions, including Bernoulli (as we have seen), Normal, Poisson, and Gamma distributions.

Why Sigmoid?

We now return to the world of machine learning and the question that motivated our study of GLMs: why sigmoid? As you may have noticed, the sigmoid function is the inverse of the canonical link function for the Bernoulli distribution. Earlier in this section we saw that we can cast the binary classification task as a parametric inference problem for a GLM. In the latter setting, we estimate our parameters \(\beta\) by maximizing the joint likelihood of our observations. In a machine learning setting, when training a binary logistic classifier, we typically minimize the cross-entropy loss (dependence on \(\beta\) is implicit in the formula):

\[\mathcal{L}_{CE}(\beta) = \frac{1}{n} \sum_{i=1}^n - y_i \log(\hat \theta_i) - (1 - y_i) \log(1 - \hat \theta_i)\]

It is easy to see that \(\mathcal{L}_{CE}\) is exactly the average negative log-likelihood of observations \((X_i, y_i)\) under a Bernoulli distribution parametrized by \(\hat \theta_i = \sigma(X_i\beta)\). Hence, gradient descent on this objective is equivalent to iteratively maximizing the log-likelihood of a Bernoulli GLM with canonical link function! Using the canonical link, namely the sigmoid function, ensures the optimal parameters \(\hat \beta\) will satisfy the desirable orthogonality condition. Moreover, the expression of the gradient of \(\mathcal{L}_{CE}\) is as given in \(\eqref{eq9}\). As an added bonus, we can thus interpret gradient descent on \(\mathcal{L}_{CE}\) as an iteratively re-weighted least squares algorithm. This is in fact identical to the procedure that statisticians use to fit GLMs; see Nelder and Wedderburn (1972) for more on this.

Recap

From a machine learning perspective, we can intuitively convince ourselves that sigmoid and cross-entropy loss possess some desirable properties for classification and learning using gradient descent. As it happens, all our intuition can be formalized through the theory of Generalized Linear Models, and we find that sigmoid and cross-entropy loss are truly the most natural choices. Here is a summary of what we have seen:

  1. A distribution (under some additional conditions) belongs to the exponential family if and only if there is a vector of sufficient statistics that has fixed dimension when sample size increases. We can always reparametrize a distribution in the exponential family with the mean parameters instead of the natural parameters. The mean parameters are defined as the expected value of the sufficient statistic, and are therefore finite-dimensional.
  2. Due to these nice properties of the exponential family, we can construct the Generalized Linear Model, which models the mean parameters of a distribution from the exponential family as a linear function of input features passed through a possibly nonlinear link function. In general there can be several possible choices for the link function for a given distribution.
  3. The canonical link is the link function which results in the natural parameters themselves being modeled as a linear function of the input features. When using the canonical link function, the MLE estimate for the model parameters satisfies some convenient and desirable properties.
  4. The logistic regression model is a special case of a GLM with Bernoulli distribution. Using the canonical link for this model is equivalent to using sigmoid activation, and cross-entropy loss is exactly equal to the negative log-likelihood for the GLM. Hence our standard binary classification model enjoys all the aforementioned desirable properties of the GLM with canonical link.

References

  • J. A. Nelder and R. W. M. Wedderburn (1972), Generalized Linear Models, Journal of the Royal Statistical Society. (pdf)
  • Kevin P. Murphy (2013) Machine Learning: A Probabilistic Perspective. (Ch. 9)
  • Probabilistic Graphical Model lecture slides, CMU (pdf)
  • David Blei, “Exponential Families” (pdf)
  • CSC412 Winter 2020 Course Notes, University of Toronto (link)
  • Probabilistic ML lecture notes, Princeton University (pdf)
  • Generalized Linear Models lecture slides, University of Michigan (pdf)