Logistic Regression: Deriving the Sigmoid and the Cross Entropy from Maximum Likelihood

0. Key points

• In a classification problem the output is a label, $0$ or $1$. Fitting least squares directly to such data produces predictions outside the range of a probability, and lets far-away points that ought to be irrelevant to the decision drag the boundary around.
• To turn the linear score $\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}$ into a probability we use the sigmoid function $\sigma(z) = 1/(1+e^{-z})$. This is not an arbitrary choice: it is equivalent to assuming that the log odds (the logit) is linear.
• The loss function is not something we pick by hand. Writing down the maximum likelihood estimate for a Bernoulli model, the negative log likelihood is the cross entropy error.
• The gradient takes the surprisingly simple form $\nabla L(\boldsymbol{w}) = \sum_i (\mu_i - y_i)\boldsymbol{x}_i$. The derivative of the sigmoid cancels against the derivative of the cross entropy, and this cancellation is what keeps learning fast.
• The Hessian is $X^{\mathsf{T}}SX \succeq O$, so $L$ is convex. But the stationarity condition is a transcendental equation: unlike linear regression, there is no closed-form solution. Hence one has no option but to search numerically using derivatives, which leads directly to gradient descent in the next chapter.
• When the data are linearly separable the maximum likelihood estimate fails to exist and the weights diverge. Adding $L^2$ regularisation restores existence and uniqueness of the minimiser, and this coincides with MAP estimation under a Gaussian prior.

1. Motivation: what breaks when we fit a line to labels 0 and 1

1.1. The classification problem

In Linear regression and least squares the output was a real number, as when predicting weight from height. Most problems one actually wants to solve, however, are not of that kind.

• Is this email spam or not?
• Does a patient with this set of test values have the disease or not?
• Is the animal in this image a cat or not?

In each case the output is a yes/no choice. Such a problem is called a binary classification problem. Mathematically, we formulate it as predicting a label $y \in \{0, 1\}$ from a feature vector $\boldsymbol{x} \in \mathbb{R}^{d}$.

A naive question arises at once. Labels $y$ are, after all, numbers, so why not simply use linear regression? Fit $\hat{y} = \boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}$ by least squares and declare “yes” when $\hat{y} \ge 0.5$ and “no” otherwise. This looks plausible. Trying it out makes clear what goes wrong.

1.2. Least squares applied to $0/1$ labels

Consider predicting the pass/fail outcome $y$ of an examination ($1$ means pass) from the study time $t$ in hours, with the following four data points.

$t$	1	2	3	4
$y$	0	0	1	1

Fit $\hat{y} = at + b$ by least squares. Since $\bar{t} = 2.5$, $\bar{y} = 0.5$, $\sum_i (t_i - \bar{t})^2 = 2.25 + 0.25 + 0.25 + 2.25 = 5$ and $\sum_i (t_i - \bar{t})(y_i - \bar{y}) = 0.75 + 0.25 + 0.25 + 0.75 = 2$, we get

$$a = \frac{2}{5} = 0.4, \qquad b = 0.5 - 0.4 \times 2.5 = -0.5 .$$

The decision boundary is $\hat{y} = 0.5$, that is $t = 2.5$, and all four points are classified correctly. So far so good. But look at the predicted values themselves: $\hat{y}(1) = -0.1$ and $\hat{y}(4) = 1.1$. We have obtained a negative probability and a probability exceeding $1$. “Your probability of passing is $-10\%$” carries no meaning.

The trouble is not merely cosmetic. Add the point $(t, y) = (20, 1)$ — someone studied for 20 hours and passed, an entirely unremarkable observation. Now $\bar{t} = 6$, $\bar{y} = 0.6$, $\sum_i (t_i - \bar{t})^2 = 25 + 16 + 9 + 4 + 196 = 250$ and $\sum_i (t_i - \bar{t})(y_i - \bar{y}) = 3 + 2.4 - 1.2 - 0.8 + 5.6 = 9$, so

$$a = \frac{9}{250} = 0.036, \qquad b = 0.6 - 0.036 \times 6 = 0.384 .$$

The boundary moves to $0.036t + 0.384 = 0.5$, that is $t = 3.22\ldots$. Consequently the point $t = 3$ now has $\hat{y}(3) = 0.492 < 0.5$, so a point that had been classified correctly becomes misclassified.

Why does this happen? The squared error $(\hat{y}_i - y_i)^2$ penalises $\hat{y}_i$ for being far from $y_i$. From the point of view of classification, however, the point at $t = 20$ is equally correct whether $\hat{y} = 1.1$ or $\hat{y} = 5$: both are comfortably on the “pass” side. The squared error counts being too deep on the correct side as an error, and flattens the line in order to reduce it. In short, the squared error is the wrong objective for classification.

1.3. What we need

Two requirements now stand out.

1. A mechanism squeezing the output into $(0,1)$, so that predictions can be read as probabilities.
2. A loss function derived from a probability model, so that the difference between “answered $0.9$ and was right” and “answered $0.55$ and was right” is measured by a non-arbitrary criterion.

Logistic regression supplies both. The sigmoid function (Definition 3.1) answers the first, and the cross entropy error coming from maximum likelihood (Definition 4.1) answers the second. Only at the last stage, when we come to minimise that loss, does differentiation become genuinely indispensable: linear regression needed only the normal equations, a system of linear equations, whereas logistic regression admits no closed-form solution (Remark 5.4).

flowchart LR
A["feature vector x"] --> B["linear score z = w·x"]
B --> C["probability p = sigmoid(z)"]
C --> D["negative log likelihood = cross entropy error L"]
D --> E["gradient grad L = sum of (p - y) x"]
E --> F["update weights w"]
F -.-> B

The overall picture of logistic regression: a linear score becomes a probability, the probability becomes a loss, and the gradient of the loss corrects the weights.

2. Preliminaries: notation and assumptions

Throughout, the data consist of $n$ pairs $(\boldsymbol{x}_1, y_1), \ldots, (\boldsymbol{x}_n, y_n)$ with $\boldsymbol{x}_i \in \mathbb{R}^{d}$ and $y_i \in \{0,1\}$. The intercept (bias) is absorbed into the feature vector: the first component of every $\boldsymbol{x}_i$ is taken to be $1$, so that the corresponding weight $w_1$ plays the role of the intercept. With this convention no intercept appears explicitly and every formula below is uniformly of the form $\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}$.

Let the design matrix $X$ be the $n \times d$ matrix whose $i$-th row is $\boldsymbol{x}_i^{\mathsf{T}}$ (the same notation as the design matrix（Definition 2.1）［Linear Regression and Least Squares］ of Linear regression and least squares). We write $\boldsymbol{y} = (y_1, \ldots, y_n)^{\mathsf{T}} \in \mathbb{R}^{n}$ for the vector of labels.

Probabilistically, we assume that the $y_i$ are conditionally independent given the $\boldsymbol{x}_i$. We do not model the distribution of $\boldsymbol{x}_i$ itself in any way (we return to this point in Example 3.4). For the basics of random variables and expectation see Random variables and expectation.

Differentiation with respect to a vector means collecting the partial derivatives, $\nabla f(\boldsymbol{w}) = (\partial f/\partial w_1, \ldots, \partial f/\partial w_d)^{\mathsf{T}}$, and the Hessian is $(\nabla^2 f)_{jk} = \partial^2 f / \partial w_j \partial w_k$. For details see the definition of the Hessian（Definition 7.3）［多変数関数の微分と偏微分］ in Differentiation of functions of several variables. For a symmetric matrix $A$, we write $A \succeq O$ for positive semidefiniteness and $A \succ O$ for positive definiteness.

3. Turning outputs into probabilities: the sigmoid and the logit

3.1. The sigmoid function

Definition 3.1（Sigmoid function (standard logistic function)）

Define $\sigma : \mathbb{R} \to \mathbb{R}$ by

$$\sigma(z) = \frac{1}{1 + e^{-z}} .$$

This function is called the sigmoid function, or the standard logistic function.

The name comes from the shape of the graph, an $S$ (from the stem of the Greek letter sigma, plus eides, “shaped like”). It was originally introduced by Verhulst in 1838 as the solution of the differential equation $\frac{dp}{dt} = p(1-p)$ describing population growth. Note that this differential equation is precisely property (3) below.

The graph of the sigmoid function. It passes through 0.5 at z = 0 and approaches 0 and 1 at the two ends.

Proposition 3.2（Basic properties of the sigmoid）

For the function $\sigma$ of Definition 3.1, the following hold.

1. For every $z \in \mathbb{R}$ we have $0 < \sigma(z) < 1$; moreover $\sigma$ is of class $C^{\infty}$ on $\mathbb{R}$, strictly increasing, and $\lim_{z \to -\infty} \sigma(z) = 0$, $\lim_{z \to +\infty} \sigma(z) = 1$.
2. For every $z$, $\sigma(-z) = 1 - \sigma(z)$. In particular $\sigma(0) = 1/2$.
3. For every $z$, $\sigma'(z) = \sigma(z)\bigl(1 - \sigma(z)\bigr) = \sigma(z)\,\sigma(-z) > 0$.
4. $\sigma : \mathbb{R} \to (0,1)$ is a bijection, and its inverse is $\sigma^{-1}(p) = \log \dfrac{p}{1-p}$ for $0 < p < 1$.

Proof（Proposition 3.2）

(1) For every $z$ we have $e^{-z} > 0$, so $1 + e^{-z} > 1 > 0$ and therefore $0 < \sigma(z) = 1/(1+e^{-z}) < 1$. The map $z \mapsto e^{-z}$ is $C^{\infty}$ and the denominator $1+e^{-z}$ never vanishes, so the quotient $\sigma$ is $C^{\infty}$ as well. Strict monotonicity follows from $\sigma' > 0$, proved in (3). As for the limits: as $z \to -\infty$ we have $e^{-z} \to +\infty$, hence $\sigma(z) \to 0$; as $z \to +\infty$ we have $e^{-z} \to 0$, hence $\sigma(z) \to 1$.

(2) By definition $\sigma(-z) = 1/(1+e^{z})$. On the other hand

$$1 - \sigma(z) = 1 - \frac{1}{1+e^{-z}} = \frac{(1+e^{-z}) - 1}{1+e^{-z}} = \frac{e^{-z}}{1+e^{-z}} ,$$

and multiplying numerator and denominator by $e^{z}$ turns this into $\dfrac{1}{e^{z}+1}$, which is $\sigma(-z)$. Setting $z = 0$ gives $\sigma(0) = 1 - \sigma(0)$, that is $\sigma(0) = 1/2$.

(3) Apply the chain rule to $\sigma(z) = (1+e^{-z})^{-1}$. The outer derivative is $-(1+e^{-z})^{-2}$ and the derivative of the inner function $1+e^{-z}$ is $-e^{-z}$, so

$$\sigma'(z) = -(1+e^{-z})^{-2} \cdot (-e^{-z}) = \frac{e^{-z}}{(1+e^{-z})^{2}} .$$

On the other hand, the intermediate computation in (2) gives $1 - \sigma(z) = \dfrac{e^{-z}}{1+e^{-z}}$, whence

$$\sigma(z)\bigl(1-\sigma(z)\bigr) = \frac{1}{1+e^{-z}} \cdot \frac{e^{-z}}{1+e^{-z}} = \frac{e^{-z}}{(1+e^{-z})^{2}} ,$$

and the two agree. Finally, $1 - \sigma(z) = \sigma(-z)$ by (2), so $\sigma'(z) = \sigma(z)\sigma(-z)$. By (1) both $\sigma(z) > 0$ and $\sigma(-z) > 0$, so $\sigma'(z) > 0$.

(4) By (3) the function $\sigma$ is strictly increasing, hence injective. By (1) it is continuous, its range is contained in $(0,1)$, and its limits at the two ends are $0$ and $1$, so by the intermediate value theorem it attains every value in $(0,1)$. Hence $\sigma : \mathbb{R} \to (0,1)$ is a bijection. The inverse is found by solving $p = 1/(1+e^{-z})$ for $z$. Taking reciprocals gives $1 + e^{-z} = 1/p$, that is $e^{-z} = (1-p)/p$. Taking logarithms gives $-z = \log\dfrac{1-p}{p}$, hence $z = \log\dfrac{p}{1-p}$.

∎

Property (3) is the identity used most often in this article. The fact that the derivative of the sigmoid is a polynomial in the sigmoid itself is what will make the gradient computation come out clean (Theorem 5.1).

3.2. Why this function — the logit and Bayes’ theorem

There are plenty of smooth increasing functions with values in $(0,1)$; the cumulative distribution function $\Phi$ of the standard normal would do, and the resulting model is called probit regression. So why the sigmoid? The answer lies in part (4) of Proposition 3.2.

Definition 3.3（Odds and logit）

For $0 < p < 1$, the quantity $\dfrac{p}{1-p}$ is called the odds of the probability $p$, and its logarithm

$$\operatorname{logit}(p) = \log \frac{p}{1-p}$$

is called the logit, or the log odds. By Proposition 3.2 (4) we have $\operatorname{logit} = \sigma^{-1}$.

The odds is the ratio of the number of favourable cases to the number of unfavourable ones. It is the same usage as “3 to 1” in horse racing or sport: for $p = 0.75$ the odds are $3$, that is “3 to 1”. A probability is confined to the bounded interval $[0,1]$, whereas the odds ranges over $(0, \infty)$ and its logarithm, the logit, over all of $\mathbb{R}$. The role of the logit is to convert a probability into a quantity one may move linearly.

Consequently, setting $\mu = \sigma(\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x})$ is exactly equivalent to applying the logit to both sides:

$$\log \frac{\mu}{1-\mu} = \boldsymbol{w}^{\mathsf{T}}\boldsymbol{x} .$$

So the assumption behind logistic regression is neither that the probability is linear nor that it has a sigmoid shape; it is the single statement that the log odds is a linear combination of the features. The sigmoid is nothing but that assumption solved for the probability.

The assumption does hold in natural situations, as the following example shows.

Example 3.4（The sigmoid emerges from two normal distributions）

Let the prior probabilities of the classes $y \in \{0,1\}$ be $\pi_1 = P(y=1)$ and $\pi_0 = 1-\pi_1$, and let $p(\boldsymbol{x} \mid y=k)$ be the class-conditional densities of the features. By Bayes' theorem（Theorem 2.2）［確率論とベイズ統計］,

$$P(y=1 \mid \boldsymbol{x}) = \frac{p(\boldsymbol{x}\mid y=1)\pi_1}{p(\boldsymbol{x}\mid y=1)\pi_1 + p(\boldsymbol{x}\mid y=0)\pi_0} = \frac{1}{1 + \exp(-a)}, \qquad a = \log \frac{p(\boldsymbol{x}\mid y=1)\pi_1}{p(\boldsymbol{x}\mid y=0)\pi_0}.$$

The middle step is nothing more than dividing numerator and denominator by $p(\boldsymbol{x}\mid y=1)\pi_1$ and using $\dfrac{p(\boldsymbol{x}\mid y=0)\pi_0}{p(\boldsymbol{x}\mid y=1)\pi_1} = e^{-a}$. Thus the sigmoid appears with no assumption at all, because $a$ is precisely the log odds.

What remains is the question whether $a$ is an affine function of $\boldsymbol{x}$. If both classes are normal with a common (invertible) covariance matrix $\Sigma$, say $N(\boldsymbol{\mu}_1, \Sigma)$ and $N(\boldsymbol{\mu}_0, \Sigma)$, the normalising constants cancel and

$$\begin{aligned} a &= \log\frac{\pi_1}{\pi_0} - \tfrac12 (\boldsymbol{x}-\boldsymbol{\mu}_1)^{\mathsf{T}}\Sigma^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_1) + \tfrac12 (\boldsymbol{x}-\boldsymbol{\mu}_0)^{\mathsf{T}}\Sigma^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_0) \\[2pt] &= (\boldsymbol{\mu}_1 - \boldsymbol{\mu}_0)^{\mathsf{T}}\Sigma^{-1}\boldsymbol{x} \;-\; \tfrac12 \boldsymbol{\mu}_1^{\mathsf{T}}\Sigma^{-1}\boldsymbol{\mu}_1 + \tfrac12 \boldsymbol{\mu}_0^{\mathsf{T}}\Sigma^{-1}\boldsymbol{\mu}_0 + \log\frac{\pi_1}{\pi_0} \end{aligned}$$

In the second line we used that the quadratic term $-\tfrac12\boldsymbol{x}^{\mathsf{T}}\Sigma^{-1}\boldsymbol{x}$ comes out of both brackets and cancels (it would not cancel if the covariance matrices differed). What is left is affine in $\boldsymbol{x}$.

Let us put in numbers in one dimension. With $\mu_0 = 0$, $\mu_1 = 2$, variance $1$ and $\pi_1 = \pi_0 = 1/2$,

$$a = -\tfrac12 (x-2)^2 + \tfrac12 x^2 = 2x - 2, \qquad P(y=1\mid x) = \sigma(2x-2).$$

The boundary $P = 1/2$ is at $x = 1$, the midpoint of the two means. Logistic regression may be read as the model that estimates the coefficients $(-2, 2)$ of this $a$ directly, without passing through $\boldsymbol{\mu}_k$ or $\Sigma$.

Remark 3.5

Example 3.4 says that assuming normal distributions leads to logistic regression, but the converse fails. Many class-conditional distributions besides the normal give linear log odds (most of the exponential family does), and logistic regression covers all of them at once. The stance of modelling $P(y \mid \boldsymbol{x})$ directly without ever building $p(\boldsymbol{x}\mid y)$ is called a discriminative model. The comparison with generative models is taken up in The role of probability and Bayesian statistics.

3.3. The logistic regression model and how to read its coefficients

Definition 3.6（Logistic regression model）

For a parameter $\boldsymbol{w} \in \mathbb{R}^{d}$, the model defining the conditional distribution of the label $y \in \{0,1\}$ given a feature $\boldsymbol{x} \in \mathbb{R}^{d}$ by

$$P(y = 1 \mid \boldsymbol{x};\boldsymbol{w}) = \sigma(\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}), \qquad P(y = 0 \mid \boldsymbol{x};\boldsymbol{w}) = 1 - \sigma(\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x})$$

is called the logistic regression model. We call $z = \boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}$ the logit or the score, and $\mu = \sigma(z)$ the predicted probability. The two formulas combine into

$$P(y \mid \boldsymbol{x};\boldsymbol{w}) = \mu^{y}(1-\mu)^{1-y}, \qquad \mu = \sigma(\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x})$$

(substituting $y=1$ gives $\mu$, and $y=0$ gives $1-\mu$). That is, $y$ follows a Bernoulli distribution with success probability $\mu$.

Example 3.7（A coefficient is a multiplier of the odds）

Suppose a model predicting the probability of passing from study time $t$ has been estimated as $\log\dfrac{\mu}{1-\mu} = -3 + 1.2\,t$. How should the coefficient $1.2$ be read?

The correct reading is: increasing $t$ by $1$ increases the log odds by $1.2$, that is, multiplies the odds by $e^{1.2} = 3.32\ldots$. Let us check.

• $t = 2$: $z = -0.6$, $\mu = \sigma(-0.6) = 0.3543$, odds $= e^{-0.6} = 0.5488$.
• $t = 3$: $z = 0.6$, $\mu = \sigma(0.6) = 0.6457$, odds $= e^{0.6} = 1.8221$. The odds ratio is $1.8221/0.5488 = 3.320 = e^{1.2}$.
• $t = 5$: $z = 3$, $\mu = 0.9526$, odds $= e^{3} = 20.09$.
• $t = 6$: $z = 4.2$, $\mu = 0.9852$, odds $= e^{4.2} = 66.69$. The odds ratio is again $66.69/20.09 = 3.320 = e^{1.2}$.

The odds ratio is constant everywhere, but the increase in probability is not. From $t : 2 \to 3$ the probability moves a great deal, $0.354 \to 0.646$ (a gain of $+0.29$), whereas from $t : 5 \to 6$ it moves only $0.953 \to 0.985$ ($+0.03$). Where the probability is already close to $1$, tripling the odds barely raises it. Reading the coefficient $1.2$ as “raises the probability by $1.2$” is simply wrong.

4. Where the loss comes from: maximum likelihood and the cross entropy error

4.1. The likelihood

Definition 3.6 writes the data-generating rule as a probability. The standard principle for fixing the parameters of such a model is maximum likelihood: choose the parameter that makes the data at hand as probable as possible.

Write $\mu_i(\boldsymbol{w}) = \sigma(\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i)$. By the conditional independence assumption of §2, the joint probability of the observed labels $y_1, \ldots, y_n$ — the likelihood — is a product:

$$\mathcal{L}(\boldsymbol{w}) = \prod_{i=1}^{n} P(y_i \mid \boldsymbol{x}_i;\boldsymbol{w}) = \prod_{i=1}^{n} \mu_i^{\,y_i}(1-\mu_i)^{1-y_i}.$$

A product is awkward, so we take logarithms. Since $\log$ is strictly increasing, the maximisers of $\mathcal{L}$ and of $\log\mathcal{L}$ coincide exactly. Flipping the sign, as is customary in optimisation, produces the following quantity.

Definition 4.1（Cross entropy error (negative log likelihood)）

For data $(\boldsymbol{x}_i, y_i)_{i=1}^{n}$ and the model of Definition 3.6, writing $\mu_i = \sigma(\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i)$, the quantity

$$L(\boldsymbol{w}) = -\sum_{i=1}^{n} \Bigl[\, y_i \log \mu_i + (1-y_i)\log(1-\mu_i) \,\Bigr]$$

is called the cross entropy error, or the negative log likelihood. Since $0 < \mu_i < 1$ (Proposition 3.2 (1)) the logarithms are always defined, and $L(\boldsymbol{w}) > 0$.

4.2. The name “cross entropy”

The name comes from information theory. For two probability distributions $p, q$ on a finite set, $H(p,q) = -\sum_{k} p_k \log q_k$ is called their cross entropy. The $i$-th term of Definition 4.1 is exactly the cross entropy between the distribution determined by the label, $p = (1-y_i,\; y_i)$ (since $y_i$ is $0$ or $1$, this is a distribution concentrated at a point), and the model distribution $q = (1-\mu_i,\; \mu_i)$.

Moreover there is the decomposition $H(p,q) = H(p) + D_{\mathrm{KL}}(p \,\|\, q)$, and here $p$ is a point mass, so its entropy is $H(p) = 0$. Therefore

$$L(\boldsymbol{w}) = \sum_{i=1}^{n} D_{\mathrm{KL}}\bigl(\,\delta_{y_i} \,\big\|\, \mathrm{Ber}(\mu_i)\,\bigr)$$

The information-theoretic reading is thus that lowering the cross entropy error is the same as bringing the model’s predictive distribution closer to the distribution of the observed labels.

Proposition 4.2（Equivalence of maximum likelihood and cross entropy minimisation）

Let $\mathcal{L}$ be the likelihood above and $L$ the cross entropy error of Definition 4.1. For every $\boldsymbol{w} \in \mathbb{R}^{d}$ we have $L(\boldsymbol{w}) = -\log\mathcal{L}(\boldsymbol{w})$, and consequently, as sets,

$$\operatorname*{arg\,max}_{\boldsymbol{w}\in\mathbb{R}^{d}} \mathcal{L}(\boldsymbol{w}) = \operatorname*{arg\,min}_{\boldsymbol{w}\in\mathbb{R}^{d}} L(\boldsymbol{w})$$

(the equality holds including the case where both sides are empty).

Proof（Proposition 4.2）

Every factor of $\mathcal{L}(\boldsymbol{w}) = \prod_i \mu_i^{y_i}(1-\mu_i)^{1-y_i}$ is strictly positive by Proposition 3.2 (1), so $\mathcal{L}(\boldsymbol{w}) > 0$ and its logarithm exists. The logarithm of a product is the sum of the logarithms, so

$$\log \mathcal{L}(\boldsymbol{w}) = \sum_{i=1}^{n}\Bigl[\, y_i \log\mu_i + (1-y_i)\log(1-\mu_i) \,\Bigr] = -L(\boldsymbol{w}).$$

The map $t \mapsto -t$ is a strictly decreasing bijection of $\mathbb{R}$, so (together with the monotonicity of $\log$) the inequalities $\mathcal{L}(\boldsymbol{w}) \ge \mathcal{L}(\boldsymbol{w}')$ and $L(\boldsymbol{w}) \le L(\boldsymbol{w}')$ are equivalent. Hence the set of maximisers of $\mathcal{L}$ and the set of minimisers of $L$ coincide.

∎

In other words, a loss function is not designed but derived from a probability model. One may restate the failure of the squared error in §1.2 by saying that the squared error is the negative log likelihood of a different model, one in which the output is normally distributed (maximum likelihood under Gaussian noise（Proposition 3.3）［確率論とベイズ統計］). Labels taking the values $0$ and $1$ are not normally distributed.

Tip

In an implementation, computing $\log(1+e^{z})$ naively overflows when $z$ is large. The $i$-th term of Definition 4.1 can be collapsed into the single expression $\log(1+e^{z_i}) - y_i z_i$ (this follows from $\log\mu = -\log(1+e^{-z})$, $\log(1-\mu) = -\log(1+e^{z})$ and $\log(1+e^{-z}) = \log(1+e^{z}) - z$), and rewriting it further as $\log(1+e^{z}) = \max(z,0) + \log\bigl(1+e^{-|z|}\bigr)$ keeps the argument of the exponential at most $0$, so the computation is safe.

4.3. Solving the simplest case

Example 4.3（The intercept-only model can be solved explicitly）

Consider a model with no features, only an intercept: $d = 1$ with $\boldsymbol{x}_i = (1)$. Then $\mu_i = \sigma(b)$ is a constant independent of $i$. If $k$ of the $n$ observations have $y_i = 1$, then

$$L(b) = -\bigl[\, k \log\sigma(b) + (n-k)\log(1-\sigma(b)) \,\bigr].$$

Differentiate. By Proposition 3.2 (3) we have $\dfrac{d}{db}\log\sigma(b) = \dfrac{\sigma'(b)}{\sigma(b)} = 1-\sigma(b)$, and likewise $\dfrac{d}{db}\log(1-\sigma(b)) = \dfrac{-\sigma'(b)}{1-\sigma(b)} = -\sigma(b)$. Hence

$$L'(b) = -\bigl[\, k(1-\sigma(b)) - (n-k)\sigma(b) \,\bigr] = n\,\sigma(b) - k .$$

The equation $L'(b) = 0$ is equivalent to $\sigma(b) = k/n$. If $0 < k < n$ then $k/n \in (0,1)$, so by Proposition 3.2 (4) there is a unique solution,

$$\hat{b} = \operatorname{logit}\!\left(\frac{k}{n}\right) = \log\frac{k}{n-k}.$$

For instance with $n = 100$ and $k = 30$ we get $\hat{b} = \log(30/70) = \log(3/7) = -0.8473$, so the predicted probability is $\sigma(-0.8473) = 0.30$, which is the empirical rate of positives itself. Maximum likelihood returns the obvious answer, as it should.

If, on the other hand, $k = 0$ or $k = n$, then $k/n$ lies outside $(0,1)$ and $L'(b) = 0$ has no solution. For $k = n$ the loss $L(b) = n\log(1+e^{-b})$ tends to $0$ as $b \to +\infty$ but never attains it: the maximum likelihood estimate does not exist. This is the simplest instance of Theorem 6.3.

5. Differentiation becomes necessary: the gradient and its meaning

5.1. The gradient formula

In Example 4.3 there was a single parameter, so we could differentiate and solve. What happens for general $d$? First we compute the gradient.

Theorem 5.1（Gradient of the cross entropy error）

Fix $\boldsymbol{x}_1,\ldots,\boldsymbol{x}_n \in \mathbb{R}^{d}$ and $y_1,\ldots,y_n \in \{0,1\}$ arbitrarily, let $L$ be the cross entropy error of Definition 4.1, and put $\mu_i(\boldsymbol{w}) = \sigma(\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i)$. Then $L$ is of class $C^{\infty}$ on $\mathbb{R}^{d}$ and

$$\nabla L(\boldsymbol{w}) = \sum_{i=1}^{n} \bigl(\mu_i(\boldsymbol{w}) - y_i\bigr)\,\boldsymbol{x}_i = X^{\mathsf{T}}\bigl(\boldsymbol{\mu}(\boldsymbol{w}) - \boldsymbol{y}\bigr) ,$$

where $X$ is the $n\times d$ design matrix whose $i$-th row is $\boldsymbol{x}_i^{\mathsf{T}}$ and $\boldsymbol{\mu}(\boldsymbol{w}) = (\mu_1,\ldots,\mu_n)^{\mathsf{T}}$.

Proof（Theorem 5.1）

Put $z_i = \boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i = \sum_{j=1}^{d} w_j x_{ij}$ and $\mu_i = \sigma(z_i)$. Each $z_i$ is affine in $\boldsymbol{w}$ hence $C^{\infty}$, $\sigma$ is $C^{\infty}$ (Proposition 3.2 (1)), and $\log$ is $C^{\infty}$ on $(0,\infty)$ with $\mu_i, 1-\mu_i \in (0,1)$; so $L$, being a finite sum of compositions of these, is $C^{\infty}$.

Write the $i$-th term as $\ell_i = -\bigl[y_i\log\mu_i + (1-y_i)\log(1-\mu_i)\bigr]$ and apply the chain rule along $\mu_i \to z_i \to w_j$.

Step 1: differentiate with respect to $\mu_i$.

$$\frac{\partial \ell_i}{\partial \mu_i} = -\frac{y_i}{\mu_i} + \frac{1-y_i}{1-\mu_i} = \frac{-y_i(1-\mu_i) + \mu_i(1-y_i)}{\mu_i(1-\mu_i)} = \frac{-y_i + y_i\mu_i + \mu_i - \mu_i y_i}{\mu_i(1-\mu_i)} = \frac{\mu_i - y_i}{\mu_i(1-\mu_i)} .$$

Along the way we put the fractions over the common denominator $\mu_i(1-\mu_i)$ and used that $y_i\mu_i$ and $-\mu_i y_i$ cancel in the numerator.

Step 2: differentiate with respect to $z_i$. By Proposition 3.2 (3), $\dfrac{d\mu_i}{dz_i} = \sigma'(z_i) = \mu_i(1-\mu_i)$.

Step 3: differentiate with respect to $w_j$. From $z_i = \sum_{j} w_j x_{ij}$ we get $\dfrac{\partial z_i}{\partial w_j} = x_{ij}$.

Multiplying the three, the denominator $\mu_i(1-\mu_i)$ of Step 1 cancels against the factor $\mu_i(1-\mu_i)$ of Step 2:

$$\frac{\partial \ell_i}{\partial w_j} = \frac{\mu_i - y_i}{\mu_i(1-\mu_i)} \cdot \mu_i(1-\mu_i) \cdot x_{ij} = (\mu_i - y_i)\,x_{ij}.$$

The cancellation is legitimate because $\mu_i(1-\mu_i) \ne 0$, which follows from $0 < \mu_i < 1$ (Proposition 3.2 (1)). Summing over $i$ and collecting $j = 1,\ldots,d$ gives

$$\nabla L(\boldsymbol{w}) = \sum_{i=1}^{n}(\mu_i - y_i)\boldsymbol{x}_i .$$

Finally, since $\boldsymbol{x}_i$ is the $i$-th row of $X$, we have $\sum_i (\mu_i - y_i)\boldsymbol{x}_i = X^{\mathsf{T}}(\boldsymbol{\mu}-\boldsymbol{y})$ (the columns of $X^{\mathsf{T}}$ are the $\boldsymbol{x}_i$, so multiplying $X^{\mathsf{T}}$ by a vector forms a linear combination of those columns).

∎

This formula looks just like the one for linear regression, whose least squares gradient was $X^{\mathsf{T}}(X\boldsymbol{w} - \boldsymbol{y})$. The only difference is that the prediction has changed from $X\boldsymbol{w}$ to $\sigma(X\boldsymbol{w})$. The structure — weight the residual (prediction minus observation) by the features and add up — is shared.

Corollary 5.2（Mean calibration of a model with an intercept）

Suppose the model contains an intercept, that is, for some $j_0$ we have $x_{i j_0} = 1$ for all $i$. Then every $\boldsymbol{w}^{*}$ satisfying $\nabla L(\boldsymbol{w}^{*}) = \boldsymbol{0}$ obeys

$$\frac{1}{n}\sum_{i=1}^{n} \mu_i(\boldsymbol{w}^{*}) = \frac{1}{n}\sum_{i=1}^{n} y_i .$$

That is, the mean predicted probability equals the proportion of positives in the data.

Proof（Corollary 5.2）

By Theorem 5.1 the $j_0$-th component of $\nabla L$ is $\sum_{i}(\mu_i - y_i)x_{ij_0}$. By hypothesis $x_{ij_0} = 1$, so this equals $\sum_i (\mu_i - y_i)$. Since $\nabla L(\boldsymbol{w}^{*}) = \boldsymbol{0}$, this component vanishes as well, that is $\sum_i \mu_i = \sum_i y_i$. Dividing both sides by $n$ gives the claim.

∎

Corollary 5.2 guarantees that a maximum-likelihood logistic regression is “right on average”. If it assigns an average probability of $0.3$ to 100 people, then exactly 30 of them were positives. Example 4.3 is nothing but the case $d=1$ of this corollary.

5.2. What the clean cancellation means

The cancellation in the proof of Theorem 5.1 is no accident: the sigmoid and the cross entropy are a pair chosen to be combined that way. Seeing what happens with the squared error instead makes the point clear.

Example 5.3（With squared error the gradient vanishes, and convexity is lost too）

Using the squared error $E(\boldsymbol{w}) = \frac12\sum_i (\mu_i - y_i)^2$ with the same model $\mu_i = \sigma(\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i)$ changes only Step 1 of the proof of Theorem 5.1, to $\partial E/\partial\mu_i = \mu_i - y_i$; the factor $\mu_i(1-\mu_i)$ from Step 2 then survives uncancelled:

$$\frac{\partial E}{\partial w_j} = \sum_{i=1}^{n} (\mu_i - y_i)\,\mu_i(1-\mu_i)\,x_{ij}.$$

To see what the extra factor $\mu_i(1-\mu_i)$ does, take a single point that is “confidently wrong”: $d = 1$, $x = 1$, $y = 1$, $w = -10$. Then $\mu = \sigma(-10) = 4.5398\times 10^{-5}$, so

• cross entropy gradient: $\mu - y = -0.99995$;
• squared error gradient: $(\mu-y)\mu(1-\mu) = -4.5394\times 10^{-5}$.

The ratio is about $22028$. At the point where the model is most badly wrong, the squared error learns essentially nothing. The reason is that the sigmoid saturates and $\sigma' \approx 0$; this is the simplest form of the phenomenon known as vanishing gradients.

Worse still, this $E$ is not even convex. In the same one-point setting, put $s = \sigma(-w) = 1-\mu$, so that $E(w) = \frac12 s^2$. Using $\dfrac{ds}{dw} = -\sigma'(-w) = -s(1-s)$ (Proposition 3.2 (2),(3)),

$$E'(w) = s\cdot\frac{ds}{dw} = -s^{2}(1-s), \qquad E''(w) = \bigl(-2s + 3s^{2}\bigr)\cdot\frac{ds}{dw} = s^{2}(1-s)(2-3s).$$

Since $s \in (0,1)$, the sign of $E''$ is the sign of $2-3s$, that is, it depends on whether $s < 2/3$. At $s = 1/2$ (i.e. $w=0$) we get $E'' = 0.0625 > 0$, whereas at $s = 0.9$ (i.e. $w = -\log 9 = -2.197$) we get $E'' = -0.0567 < 0$. Convexity flips across the inflection point $w = -\log 2$.

At the same single point the cross entropy is $L(w) = -\log\sigma(w) = \log(1+e^{-w})$, with $L'(w) = -(1-\sigma(w))$ and $L''(w) = \sigma(w)(1-\sigma(w)) > 0$, so it is strictly convex. Moreover $L'(w) \to -1$ as $w \to -\infty$: the gradient does not vanish.

5.3. There is no closed-form solution

Now that we have the gradient, the maximum likelihood estimate must satisfy the stationarity condition

$$X^{\mathsf{T}}\bigl(\sigma(X\boldsymbol{w}) - \boldsymbol{y}\bigr) = \boldsymbol{0}$$

(where $\sigma$ acts componentwise). This is the decisive parting of the ways from linear regression.

Remark 5.4

The stationarity condition for linear regression is the normal equation $X^{\mathsf{T}}X\boldsymbol{w} = X^{\mathsf{T}}\boldsymbol{y}$ (the normal equation（Theorem 3.3）［Linear Regression and Least Squares］), a system of linear equations in $\boldsymbol{w}$. If $X^{\mathsf{T}}X$ is invertible, we may write $\boldsymbol{w} = (X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\boldsymbol{y}$, obtained by finitely many arithmetic operations.

The stationarity condition for logistic regression, by contrast, is a transcendental equation mixing exponentials with polynomials. Even in the case $d=1$ with $\boldsymbol{x}_i = (t_i)$ it reads

$$\sum_{i=1}^{n} \frac{t_i}{1+e^{-w t_i}} = \sum_{i=1}^{n} y_i t_i ,$$

whose left-hand side is an elementary function of $w$, yet no general formula solving it for $w$ in elementary functions is known (apart from special cases that reduce to $n=1$ and can be solved as in Example 4.3). A rigorous proof that no elementary closed form exists belongs to differential Galois theory and we do not enter into it here, but the practical consequence is plain: we must give up on solving by symbolic manipulation and search numerically instead.

And when searching numerically, the only local information telling us, from the current point $\boldsymbol{w}$, which way to move so that $L$ decreases is the gradient $\nabla L(\boldsymbol{w})$. Here differentiation ceases to be a computational device and becomes a compass for search. The method that actually runs this search is gradient descent (Definition 4.1［勾配降下法］), and the machinery for computing gradients efficiently in multilayer models is described in Neural networks and backpropagation.

5.4. A numerical example

Example 5.5（Running one gradient step by hand）

Fit the data of §1.2 ($t = 1,2,3,4$ with $y = 0,0,1,1$) with an intercept, so that $\boldsymbol{x}_i = (1, t_i)^{\mathsf{T}}$ and $\boldsymbol{w} = (b, a)^{\mathsf{T}}$.

Initial point $\boldsymbol{w} = (0,0)$. Since $z_i = 0$, we have $\mu_i = \sigma(0) = 0.5$ (Proposition 3.2 (2)). The loss is

$$L(\boldsymbol{0}) = -\sum_{i=1}^{4}\log 0.5 = 4\log 2 = 2.7726 .$$

The residuals are $\boldsymbol{\mu}-\boldsymbol{y} = (0.5,\, 0.5,\, -0.5,\, -0.5)$, so by Theorem 5.1

$$\nabla L(\boldsymbol{0}) = \begin{pmatrix} 0.5+0.5-0.5-0.5 \\ 0.5\cdot 1 + 0.5\cdot 2 - 0.5\cdot 3 - 0.5\cdot 4\end{pmatrix} = \begin{pmatrix} 0 \\ -2 \end{pmatrix}.$$

The intercept component vanishes exactly as Corollary 5.2 predicts, because $\bar{\mu} = 0.5 = \bar{y}$. The slope component is negative, so increasing $a$ decreases the loss.

One step. With learning rate $\eta = 0.1$ we set $\boldsymbol{w} \leftarrow \boldsymbol{w} - \eta\nabla L(\boldsymbol{w}) = (0,\, 0.2)$. Then $z_i = 0.2, 0.4, 0.6, 0.8$ and $\mu_i = 0.5498, 0.5987, 0.6457, 0.6900$, giving

$$L = -\bigl[\log 0.4502 + \log 0.4013 + \log 0.6457 + \log 0.6900\bigr] = 0.7981+0.9130+0.4375+0.3711 = 2.5197 .$$

Indeed the loss has decreased from $2.7726$. The new gradient, from $\boldsymbol{\mu}-\boldsymbol{y} = (0.5498,\,0.5987,\,-0.3543,\,-0.3100)$, is

$$\nabla L = \begin{pmatrix} 0.5498+0.5987-0.3543-0.3100 \\ 0.5498+1.1974-1.0630-1.2401 \end{pmatrix} = \begin{pmatrix} 0.4842 \\ -0.5559 \end{pmatrix} .$$

This time the intercept component is positive. Raising only the slope pushed all the predictions upward and broke the mean calibration. The next step will therefore lower the intercept while raising the slope.

Turning the computation above into code gives the following. Following the remark in §4.2, the loss is collapsed into the form $\log(1+e^{z}) - yz$ and then rewritten as $\log(1+e^{z}) = \max(z,0)+\log(1+e^{-|z|})$ to avoid overflow.

import numpy as np

def softplus(z):                      # a safe computation of log(1 + exp(z))
    return np.maximum(z, 0.0) + np.log1p(np.exp(-np.abs(z)))

def loss(w, X, y):
    z = X @ w
    return float(np.sum(softplus(z) - y * z))

def grad(w, X, y):                    # the gradient formula verbatim
    mu = 1.0 / (1.0 + np.exp(-(X @ w)))
    return X.T @ (mu - y)

X = np.array([[1.0, 1.0], [1.0, 2.0], [1.0, 3.0], [1.0, 4.0]])
y = np.array([0.0, 0.0, 1.0, 1.0])

w = np.zeros(2)
print(loss(w, X, y), grad(w, X, y))   # 2.772588722239781  [ 0. -2.]

for _ in range(3):
    w = w - 0.1 * grad(w, X, y)
    print(w, loss(w, X, y))

The losses printed decrease monotonically: $2.7726 \to 2.5197 \to 2.4706 \to 2.4305$.

6. Convexity: why the search works, and when it does not

Having decided to search numerically, the next thing to check is whether the search can find anything. For a general function, a point where the gradient vanishes may be a local minimum, a local maximum or a saddle point. The cross entropy error, however, has a good property.

6.1. The Hessian and convexity

Theorem 6.1（Convexity of the cross entropy error）

In the setting of Theorem 5.1, put $S(\boldsymbol{w}) = \operatorname{diag}\bigl(\mu_1(1-\mu_1), \ldots, \mu_n(1-\mu_n)\bigr)$. Then

$$\nabla^{2} L(\boldsymbol{w}) = \sum_{i=1}^{n} \mu_i(1-\mu_i)\,\boldsymbol{x}_i\boldsymbol{x}_i^{\mathsf{T}} = X^{\mathsf{T}}S(\boldsymbol{w})X .$$

This matrix is positive semidefinite for every $\boldsymbol{w}$, and consequently $L$ is a convex function on $\mathbb{R}^{d}$. If moreover $\operatorname{rank} X = d$ (the columns of $X$ are linearly independent), then $\nabla^{2}L(\boldsymbol{w}) \succ O$ for every $\boldsymbol{w}$ and $L$ is strictly convex.

Proof（Theorem 6.1）

Computation of the Hessian. By Theorem 5.1, $\dfrac{\partial L}{\partial w_j} = \sum_i (\mu_i - y_i)x_{ij}$. The $y_i$ are constants, so differentiating once more with respect to $w_k$ only the $\mu_i$ contribute:

$$\frac{\partial^{2} L}{\partial w_j \partial w_k} = \sum_{i=1}^{n} \frac{\partial \mu_i}{\partial w_k}\, x_{ij} = \sum_{i=1}^{n} \sigma'(z_i)\,\frac{\partial z_i}{\partial w_k}\, x_{ij} = \sum_{i=1}^{n} \mu_i(1-\mu_i)\, x_{ik} x_{ij} .$$

The second equality is the chain rule and the third uses Proposition 3.2 (3) together with $\partial z_i/\partial w_k = x_{ik}$. Since $x_{ij}x_{ik}$ is the $(j,k)$ entry of the matrix $\boldsymbol{x}_i\boldsymbol{x}_i^{\mathsf{T}}$, in matrix form $\nabla^2 L = \sum_i \mu_i(1-\mu_i)\boldsymbol{x}_i\boldsymbol{x}_i^{\mathsf{T}}$. As the $i$-th row of $X$ is $\boldsymbol{x}_i^{\mathsf{T}}$, this equals $X^{\mathsf{T}}S X$.

Positive semidefiniteness. For any $\boldsymbol{v}\in\mathbb{R}^{d}$,

$$\boldsymbol{v}^{\mathsf{T}}\nabla^{2}L(\boldsymbol{w})\boldsymbol{v} = \sum_{i=1}^{n}\mu_i(1-\mu_i)\,\boldsymbol{v}^{\mathsf{T}}\boldsymbol{x}_i\boldsymbol{x}_i^{\mathsf{T}}\boldsymbol{v} = \sum_{i=1}^{n}\mu_i(1-\mu_i)\,(\boldsymbol{x}_i^{\mathsf{T}}\boldsymbol{v})^{2} \;\ge\; 0 .$$

Each summand is nonnegative because $0 < \mu_i < 1$ by Proposition 3.2 (1), hence $\mu_i(1-\mu_i) > 0$, and $(\boldsymbol{x}_i^{\mathsf{T}}\boldsymbol{v})^2 \ge 0$.

Convexity. Take arbitrary $\boldsymbol{w}_0, \boldsymbol{w}_1 \in \mathbb{R}^{d}$ and put $\boldsymbol{h} = \boldsymbol{w}_1 - \boldsymbol{w}_0$ and $g(t) = L(\boldsymbol{w}_0 + t\boldsymbol{h})$. Since $L$ is $C^{\infty}$ (Theorem 5.1), $g$ is $C^{2}$ on $\mathbb{R}$, and by the chain rule $g''(t) = \boldsymbol{h}^{\mathsf{T}}\nabla^{2}L(\boldsymbol{w}_0+t\boldsymbol{h})\boldsymbol{h} \ge 0$. A function of one variable with nonnegative second derivative is convex, so $g$ is convex on $[0,1]$ and $g(t) \le (1-t)g(0) + t\,g(1)$, that is,

$$L\bigl((1-t)\boldsymbol{w}_0 + t\boldsymbol{w}_1\bigr) \le (1-t)L(\boldsymbol{w}_0) + t\,L(\boldsymbol{w}_1) \qquad (0\le t\le 1).$$

As $\boldsymbol{w}_0,\boldsymbol{w}_1$ were arbitrary, $L$ is convex.

Strict convexity. Assume $\operatorname{rank}X = d$ and let $\boldsymbol{v}\ne\boldsymbol{0}$. Suppose $\boldsymbol{v}^{\mathsf{T}}\nabla^{2}L\boldsymbol{v} = 0$ in the identity above. Since all summands are nonnegative, each must vanish, that is $\mu_i(1-\mu_i)(\boldsymbol{x}_i^{\mathsf{T}}\boldsymbol{v})^2 = 0$. As $\mu_i(1-\mu_i) > 0$, we get $\boldsymbol{x}_i^{\mathsf{T}}\boldsymbol{v} = 0$ for every $i$, which means $X\boldsymbol{v} = \boldsymbol{0}$. Since $\operatorname{rank}X = d$, the kernel of $X$ is $\{\boldsymbol{0}\}$, so $\boldsymbol{v} = \boldsymbol{0}$, a contradiction. Hence $\boldsymbol{v}\ne\boldsymbol{0}$ implies $\boldsymbol{v}^{\mathsf{T}}\nabla^{2}L\boldsymbol{v} > 0$, that is $\nabla^2 L \succ O$. In that case the function $g$ above has $g'' > 0$, so $g$ is strictly convex and therefore so is $L$.

∎

Corollary 6.2（Stationary points are global minima）

In the setting of Theorem 6.1, if $\boldsymbol{w}^{*}\in\mathbb{R}^{d}$ satisfies $\nabla L(\boldsymbol{w}^{*}) = \boldsymbol{0}$, then $\boldsymbol{w}^{*}$ is a global minimiser of $L$. Conversely, every global minimiser is a stationary point.

Proof（Corollary 6.2）

Take any $\boldsymbol{w}\in\mathbb{R}^{d}$ and put $\boldsymbol{h} = \boldsymbol{w}-\boldsymbol{w}^{*}$ and $g(t) = L(\boldsymbol{w}^{*}+t\boldsymbol{h})$. Since $g$ is $C^{2}$, Taylor’s theorem in one variable with Lagrange remainder provides $\theta\in(0,1)$ with

$$g(1) = g(0) + g'(0) + \tfrac12 g''(\theta).$$

Here $g(1) = L(\boldsymbol{w})$, $g(0) = L(\boldsymbol{w}^{*})$, $g'(0) = \nabla L(\boldsymbol{w}^{*})^{\mathsf{T}}\boldsymbol{h} = 0$ by hypothesis, and $g''(\theta) = \boldsymbol{h}^{\mathsf{T}}\nabla^{2}L(\boldsymbol{w}^{*}+\theta\boldsymbol{h})\boldsymbol{h} \ge 0$ by the positive semidefiniteness in Theorem 6.1. Therefore $L(\boldsymbol{w}) \ge L(\boldsymbol{w}^{*})$ for every $\boldsymbol{w}$. The converse follows because $L$ is differentiable, so its gradient vanishes at a global minimum (Fermat’s theorem). For Taylor’s theorem see Theorem 5.3［Mean Value Theorems and Taylor's Theorem］ in The mean value theorem and Taylor’s theorem.

∎

This is why it is legitimate to search using the gradient alone: there is no risk of being trapped in a local minimum, and wherever the gradient vanishes is the answer. Loss functions in deep learning are generally not convex, so this guarantee is a considerable advantage of logistic regression.

6.2. When the data are linearly separable the estimate does not exist

Convexity guarantees that whatever is found is globally optimal; it does not guarantee that anything is found. Indeed, in the following common situation there is no minimiser.

Theorem 6.3（No maximum likelihood estimate under linear separability）

Suppose the data $(\boldsymbol{x}_i, y_i)_{i=1}^{n}$ with $n \ge 1$ are strictly linearly separable, that is, there exists a vector $\boldsymbol{v}\in\mathbb{R}^{d}$ such that

$$y_i = 1 \implies \boldsymbol{x}_i^{\mathsf{T}}\boldsymbol{v} > 0, \qquad y_i = 0 \implies \boldsymbol{x}_i^{\mathsf{T}}\boldsymbol{v} < 0$$

for every $i$. Then, for the function $L$ of Definition 4.1,

$$\inf_{\boldsymbol{w}\in\mathbb{R}^{d}} L(\boldsymbol{w}) = 0 ,$$

but this infimum is not attained. Moreover every sequence $(\boldsymbol{w}_k)$ with $L(\boldsymbol{w}_k)\to 0$ satisfies $\|\boldsymbol{w}_k\| \to \infty$.

Proof（Theorem 6.3）

(a) $L > 0$. By Proposition 3.2 (1) we have $0 < \mu_i < 1$, so a term with $y_i = 1$, namely $-\log\mu_i$, is strictly positive because $\mu_i < 1$, and a term with $y_i = 0$, namely $-\log(1-\mu_i)$, is strictly positive because $1-\mu_i < 1$. Being a sum of $n \ge 1$ strictly positive numbers, $L(\boldsymbol{w}) > 0$ for every $\boldsymbol{w}$.

(b) $L(t\boldsymbol{v})\to 0$. Let $t > 0$ and put $c_i = \boldsymbol{x}_i^{\mathsf{T}}\boldsymbol{v}$, so that $\boldsymbol{w} = t\boldsymbol{v}$ gives $z_i = tc_i$. For a term with $y_i = 1$, Definition 3.1 gives $\log\sigma(z) = -\log(1+e^{-z})$, hence

$$-\log\sigma(tc_i) = \log\bigl(1+e^{-tc_i}\bigr) \xrightarrow{\;t\to\infty\;} \log 1 = 0$$

(by hypothesis $c_i > 0$, so $e^{-tc_i}\to 0$). For a term with $y_i = 0$, Proposition 3.2 (2) gives $1-\sigma(tc_i) = \sigma(-tc_i)$, hence

$$-\log\bigl(1-\sigma(tc_i)\bigr) = -\log\sigma(-tc_i) = \log\bigl(1+e^{tc_i}\bigr) \xrightarrow{\;t\to\infty\;} 0$$

(by hypothesis $c_i < 0$, so $e^{tc_i}\to 0$). Being a finite sum, $L(t\boldsymbol{v})\to 0$.

(c) The infimum and its non-attainment. By (a), $L > 0$; by (b), $L$ comes arbitrarily close to $0$; hence $\inf L = 0$. But by (a) no $\boldsymbol{w}$ gives $L(\boldsymbol{w}) = 0$, so the infimum is not attained.

(d) Divergence. Suppose $L(\boldsymbol{w}_k)\to 0$ and $(\boldsymbol{w}_k)$ were bounded. By the Bolzano–Weierstrass theorem some subsequence converges, $\boldsymbol{w}_{k_m}\to\boldsymbol{w}_{\infty}$. Since $L$ is continuous (indeed $C^{\infty}$ by Theorem 5.1), $L(\boldsymbol{w}_{\infty}) = \lim_m L(\boldsymbol{w}_{k_m}) = 0$, contradicting (a). Hence $(\boldsymbol{w}_k)$ is unbounded. Furthermore, if it had a bounded subsequence, that subsequence would also satisfy $L \to 0$ and the same argument would give a contradiction. Having no bounded subsequence is precisely $\|\boldsymbol{w}_k\|\to\infty$.

∎

Example 6.4（Watching the weights diverge）

The data of §1.2 ($t=1,2$ with $y=0$, and $t=3,4$ with $y=1$) are separated at $t = 2.5$. Taking $\boldsymbol{v} = (-2.5,\, 1)^{\mathsf{T}}$ in Theorem 6.3 gives $c_i = t_i - 2.5 = -1.5, -0.5, 0.5, 1.5$, so the sign conditions hold. Computing the loss along $\boldsymbol{w} = \alpha\boldsymbol{v}$ gives the following.

$\alpha$	1	2	5	10	20
$\|\boldsymbol{w}\|$	2.69	5.39	13.46	26.93	53.85
$L(\boldsymbol{w})$	1.3510	0.7237	0.1589	0.01343	0.0000908

Let us verify the value at $\alpha = 1$ by hand. Here $z_i = -1.5, -0.5, 0.5, 1.5$, the two points with $y=0$ contribute $\log(1+e^{z})$ and the two with $y=1$ contribute $\log(1+e^{-z})$, so

$$L = \log(1+e^{-1.5}) + \log(1+e^{-0.5}) + \log(1+e^{-0.5}) + \log(1+e^{-1.5}) = 2(0.2014 + 0.4741) = 1.3510 .$$

The loss decreases monotonically towards $0$ while $\|\boldsymbol{w}\|$ grows without bound. This is why running the gradient descent of Example 5.5 indefinitely makes the weights grow forever. The practical nuisance is that on separable data all predicted probabilities stick to $0$ or $1$, so the information about how confident the model is gets lost. When the number of features $d$ exceeds the number of data points $n$ the data are almost always separable, so this is no exotic scenario.

6.3. Fixing it by regularisation

Theorem 6.5（Existence and uniqueness for L2-regularised maximum likelihood）

Let $\lambda > 0$ and, for the function $L$ of Definition 4.1, put

$$L_{\lambda}(\boldsymbol{w}) = L(\boldsymbol{w}) + \frac{\lambda}{2}\|\boldsymbol{w}\|^{2} .$$

No condition whatsoever is imposed on the data $(\boldsymbol{x}_i,y_i)_{i=1}^n$ (they may be separable, and $X$ need not have full column rank). Then $L_{\lambda}$ has exactly one global minimiser $\hat{\boldsymbol{w}}_{\lambda}$ on $\mathbb{R}^{d}$, and it is the unique solution of the equation

$$X^{\mathsf{T}}\bigl(\sigma(X\hat{\boldsymbol{w}}_{\lambda}) - \boldsymbol{y}\bigr) + \lambda\,\hat{\boldsymbol{w}}_{\lambda} = \boldsymbol{0} .$$

Proof（Theorem 6.5）

Existence. By part (a) of Theorem 6.3 we have $L \ge 0$, so $L_{\lambda}(\boldsymbol{w}) \ge \frac{\lambda}{2}\|\boldsymbol{w}\|^{2}$. On the other hand, at $\boldsymbol{w} = \boldsymbol{0}$ we have $\mu_i = 1/2$ and hence $L_{\lambda}(\boldsymbol{0}) = L(\boldsymbol{0}) = n\log 2$. Taking $R = \sqrt{2n\log 2/\lambda} + 1$, for $\|\boldsymbol{w}\| > R$ we get

$$L_{\lambda}(\boldsymbol{w}) \ge \frac{\lambda}{2}\|\boldsymbol{w}\|^{2} > \frac{\lambda}{2}R^{2} > n\log 2 = L_{\lambda}(\boldsymbol{0}) .$$

Therefore the infimum of $L_{\lambda}$ over $\mathbb{R}^{d}$ coincides with its infimum over the closed ball $\bar{B}(\boldsymbol{0},R) = \{\boldsymbol{w} : \|\boldsymbol{w}\|\le R\}$. That ball is a bounded closed subset of $\mathbb{R}^d$, hence compact, and $L_{\lambda}$ is continuous, so by the Weierstrass extreme value theorem the minimum is attained at some point $\hat{\boldsymbol{w}}_{\lambda}$. This is a global minimiser over all of $\mathbb{R}^{d}$.

Uniqueness. The Hessian of $\|\boldsymbol{w}\|^{2}$ is $2I$, so $\nabla^{2}L_{\lambda}(\boldsymbol{w}) = X^{\mathsf{T}}S(\boldsymbol{w})X + \lambda I$. For any $\boldsymbol{v}\ne\boldsymbol{0}$, positive semidefiniteness from Theorem 6.1 gives

$$\boldsymbol{v}^{\mathsf{T}}\nabla^{2}L_{\lambda}\boldsymbol{v} = \boldsymbol{v}^{\mathsf{T}}X^{\mathsf{T}}SX\boldsymbol{v} + \lambda\|\boldsymbol{v}\|^{2} \ge \lambda\|\boldsymbol{v}\|^{2} > 0 .$$

Now suppose $\boldsymbol{w}_1 \ne \boldsymbol{w}_2$ were both global minimisers. Both are stationary, so $\nabla L_{\lambda}(\boldsymbol{w}_1) = \boldsymbol{0}$. Carrying out the same Taylor expansion as in the proof of Corollary 6.2 with $\boldsymbol{h} = \boldsymbol{w}_2-\boldsymbol{w}_1 \ne \boldsymbol{0}$, there is $\theta\in(0,1)$ with

$$L_{\lambda}(\boldsymbol{w}_2) = L_{\lambda}(\boldsymbol{w}_1) + 0 + \tfrac12\boldsymbol{h}^{\mathsf{T}}\nabla^{2}L_{\lambda}(\boldsymbol{w}_1+\theta\boldsymbol{h})\boldsymbol{h} \ge L_{\lambda}(\boldsymbol{w}_1) + \frac{\lambda}{2}\|\boldsymbol{h}\|^{2} > L_{\lambda}(\boldsymbol{w}_1) ,$$

contradicting the minimality of $\boldsymbol{w}_2$. Hence the minimiser is unique.

The equation. By Theorem 5.1 and $\nabla\bigl(\frac{\lambda}{2}\|\boldsymbol{w}\|^{2}\bigr) = \lambda\boldsymbol{w}$ we have $\nabla L_{\lambda}(\boldsymbol{w}) = X^{\mathsf{T}}(\sigma(X\boldsymbol{w})-\boldsymbol{y}) + \lambda\boldsymbol{w}$. As $L_{\lambda}$ is convex (a sum of the convex function $L$ and the convex function $\frac{\lambda}{2}\|\boldsymbol{w}\|^2$), the argument of Corollary 6.2 shows that being stationary and being a global minimiser are equivalent. Since the minimiser is unique, so is the solution of the stationarity equation.

∎

Remark 6.6（Regularisation is a Gaussian prior）

The term in $\lambda$ looks like an engineering trick that penalises large weights, but in the language of probability it has a natural interpretation. Regard $\boldsymbol{w}$ itself as a random variable with prior $\boldsymbol{w} \sim N(\boldsymbol{0}, \tau^{2}I)$. By Bayes’ theorem the posterior satisfies $p(\boldsymbol{w}\mid \text{data}) \propto \mathcal{L}(\boldsymbol{w})\,p(\boldsymbol{w})$, so its negative logarithm is

$$-\log p(\boldsymbol{w}\mid \text{data}) = L(\boldsymbol{w}) + \frac{1}{2\tau^{2}}\|\boldsymbol{w}\|^{2} + \text{const} .$$

This is exactly $L_{\lambda}$ with $\lambda = 1/\tau^{2}$. Minimising with $L^2$ regularisation is the same as MAP estimation (maximising the posterior) under a Gaussian prior. The correspondence also matches intuition: the smaller $\tau$ is — the more strongly we believe the weights lie near $\boldsymbol{0}$ — the larger $\lambda$ becomes. For details see L2 regularisation as MAP estimation under a Gaussian prior（Theorem 5.1）［確率論とベイズ統計］ in The role of probability and Bayesian statistics.

7. Exercises

Exercise 7.1Easy

For the function $\sigma$ of Definition 3.1, show that $\sigma''(z) = \sigma'(z)\bigl(1-2\sigma(z)\bigr)$ and verify that $\sigma$ has an inflection point at $z = 0$.

Solution

By Proposition 3.2 (3), $\sigma' = \sigma(1-\sigma) = \sigma - \sigma^{2}$. Differentiating with respect to $z$, by the product rule (or the chain rule),

$$\sigma'' = \sigma' - 2\sigma\sigma' = \sigma'(1 - 2\sigma) .$$

By Proposition 3.2 (3) we have $\sigma' > 0$, so the sign of $\sigma''$ is determined solely by the sign of $1-2\sigma(z)$. Since $\sigma$ is strictly increasing with $\sigma(0) = 1/2$ (same proposition, (2)):

• for $z < 0$ we have $\sigma(z) < 1/2$, so $1-2\sigma(z) > 0$, that is $\sigma'' > 0$ (convex);
• at $z = 0$ we have $\sigma'' = 0$;
• for $z > 0$ we have $\sigma(z) > 1/2$, so $\sigma'' < 0$ (concave).

The concavity changes across $z = 0$, so $z = 0$ is an inflection point. Since $\sigma'(0) = \frac12\cdot\frac12 = \frac14$, the tangent there has slope $1/4$, the maximal steepness of the sigmoid.

Exercise 7.2Standard

A model for the probability of contracting a certain disease has been estimated as $\log\dfrac{\mu}{1-\mu} = -4 + 0.8\,x_1 + 1.5\,x_2$, where $x_1$ is age measured in units of 10 years and $x_2$ equals $1$ for a smoker and $0$ for a non-smoker.

1. Find the probability of contracting the disease for a 50-year-old ($x_1 = 5$) non-smoker.
2. Holding everything else fixed, by what factor are the odds for a smoker larger than for a non-smoker?
3. Under this model, being a smoker raises the odds by as much as how many years of ageing?

Solution

1. We have $z = -4 + 0.8\times 5 + 1.5\times 0 = -4 + 4 = 0$, so by Proposition 3.2 (2) the probability is $\mu = \sigma(0) = 0.5$, that is $50\%$.

2. Changing $x_2$ from $0$ to $1$ increases the log odds by $1.5$, so the odds are multiplied by $e^{1.5} = 4.4817$. As in Example 3.7, this factor does not depend on the value of $x_1$. The factor by which the probability changes, however, does depend on $x_1$: at $x_1 = 5$ the probability only goes from $0.5$ to $\sigma(1.5) = 0.8176$, a factor of $1.635$.

3. We look for $\Delta$ making the log odds of a smoker $(x_1, 1)$ equal to those of an older non-smoker $(x_1 + \Delta, 0)$:

$$-4 + 0.8x_1 + 1.5 = -4 + 0.8(x_1 + \Delta) \iff 1.5 = 0.8\,\Delta \iff \Delta = 1.875 .$$

Since $x_1$ is measured in units of 10 years, this is 18.75 years. Checking numerically: a smoker with $x_1 = 5$ (aged 50) has $z = -4 + 4 + 1.5 = 1.5$, and a non-smoker with $x_1 = 6.875$ (aged 68.75) has $z = -4 + 0.8\times 6.875 = -4 + 5.5 = 1.5$; they agree. Both have probability $\sigma(1.5) = 0.8176$.

The reason this computation does not depend on the value of $x_1$ is that the log odds is a linear combination of $x_1$ and $x_2$. In a model with an interaction term $x_1x_2$, the conversion depends on age and can no longer be expressed as a fixed number of years.

Exercise 7.3Standard

Relabel by $\tilde{y}_i = 2y_i - 1 \in \{-1, +1\}$. With $z_i = \boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i$, show that the function $L$ of Definition 4.1 can be written

$$L(\boldsymbol{w}) = \sum_{i=1}^{n} \log\bigl(1 + e^{-\tilde{y}_i z_i}\bigr) ,$$

then compute the gradient from this expression and check that it agrees with Theorem 5.1.

Solution

The expression. Treat the $i$-th term case by case. If $y_i = 1$ (so $\tilde{y}_i = +1$), the term in Definition 4.1 is $-\log\mu_i = -\log\sigma(z_i)$. From $\sigma(z) = 1/(1+e^{-z})$ in Definition 3.1 we get $-\log\sigma(z_i) = \log(1+e^{-z_i}) = \log(1+e^{-\tilde{y}_i z_i})$.

If $y_i = 0$ (so $\tilde{y}_i = -1$), the term is $-\log(1-\mu_i)$. By Proposition 3.2 (2) we have $1-\sigma(z_i) = \sigma(-z_i)$, so

$$-\log(1-\mu_i) = -\log\sigma(-z_i) = \log\bigl(1+e^{z_i}\bigr) = \log\bigl(1+e^{-(-1)z_i}\bigr) = \log\bigl(1+e^{-\tilde{y}_i z_i}\bigr) ,$$

which is the same formula in both cases.

The gradient. Put $u_i = -\tilde{y}_i z_i$, so the $i$-th term is $\log(1+e^{u_i})$, and $\dfrac{d}{du}\log(1+e^{u}) = \dfrac{e^{u}}{1+e^{u}} = \sigma(u)$ (divide numerator and denominator by $e^{u}$ to recover the definition of $\sigma$). By the chain rule, $\dfrac{\partial u_i}{\partial \boldsymbol{w}} = -\tilde{y}_i\boldsymbol{x}_i$, so

$$\nabla L(\boldsymbol{w}) = -\sum_{i=1}^{n} \tilde{y}_i\,\sigma(-\tilde{y}_i z_i)\,\boldsymbol{x}_i .$$

Agreement with Theorem 5.1 is checked case by case. If $y_i = 1$: $-\tilde{y}_i\sigma(-\tilde{y}_iz_i) = -\sigma(-z_i) = -(1-\mu_i) = \mu_i - 1 = \mu_i - y_i$. If $y_i = 0$: $-\tilde{y}_i\sigma(-\tilde{y}_iz_i) = +\sigma(z_i) = \mu_i = \mu_i - y_i$. Both equal $\mu_i - y_i$, so the two expressions agree.

This form exhibits the structure “the larger the margin $\tilde{y}_i z_i$, the smaller the loss”, and puts the loss into a shape directly comparable with the hinge loss $\max(0, 1-\tilde{y}_iz_i)$ of support vector machines.

Exercise 7.4Hard

Let $\lambda > 0$ and consider $L_{\lambda}$ of Theorem 6.5 together with its unique minimiser $\hat{\boldsymbol{w}}_{\lambda}$. Show that for every $\boldsymbol{w}\in\mathbb{R}^{d}$

$$L_{\lambda}(\boldsymbol{w}) \;\ge\; L_{\lambda}(\hat{\boldsymbol{w}}_{\lambda}) + \frac{\lambda}{2}\bigl\|\boldsymbol{w}-\hat{\boldsymbol{w}}_{\lambda}\bigr\|^{2} ,$$

and use this to deduce $\|\hat{\boldsymbol{w}}_{\lambda}\| \le \sqrt{2n\log 2/\lambda}$.

Solution

The inequality. Put $\boldsymbol{h} = \boldsymbol{w} - \hat{\boldsymbol{w}}_{\lambda}$ and $g(t) = L_{\lambda}(\hat{\boldsymbol{w}}_{\lambda} + t\boldsymbol{h})$. Since $L$ is $C^{\infty}$ (Theorem 5.1) and $\frac{\lambda}{2}\|\boldsymbol{w}\|^{2}$ is a polynomial, $L_{\lambda}$ is $C^{\infty}$ and $g$ is $C^{2}$. By Taylor’s theorem there is $\theta\in(0,1)$ with

$$L_{\lambda}(\boldsymbol{w}) = g(1) = g(0) + g'(0) + \tfrac12 g''(\theta).$$

Since $\hat{\boldsymbol{w}}_{\lambda}$ is a minimiser, $\nabla L_{\lambda}(\hat{\boldsymbol{w}}_{\lambda}) = \boldsymbol{0}$ and therefore $g'(0) = \nabla L_{\lambda}(\hat{\boldsymbol{w}}_{\lambda})^{\mathsf{T}}\boldsymbol{h} = 0$. Also, as shown in the proof of Theorem 6.5, $\nabla^{2}L_{\lambda} = X^{\mathsf{T}}SX + \lambda I$, and $X^{\mathsf{T}}SX \succeq O$ by Theorem 6.1, so

$$g''(\theta) = \boldsymbol{h}^{\mathsf{T}}\bigl(X^{\mathsf{T}}S X + \lambda I\bigr)\boldsymbol{h} \ge \lambda\|\boldsymbol{h}\|^{2}.$$

Substituting gives $L_{\lambda}(\boldsymbol{w}) \ge L_{\lambda}(\hat{\boldsymbol{w}}_{\lambda}) + \frac{\lambda}{2}\|\boldsymbol{h}\|^{2}$. This property is called $\lambda$-strong convexity; it is stronger than strict convexity, asserting that the function is bounded below by a quadratic.

The upper bound. Take $\boldsymbol{w} = \boldsymbol{0}$ in the inequality. As seen in the proof of Theorem 6.5, $L_{\lambda}(\boldsymbol{0}) = L(\boldsymbol{0}) = n\log 2$ (there are $n$ terms with $\mu_i = 1/2$), so

$$n\log 2 \;\ge\; L_{\lambda}(\hat{\boldsymbol{w}}_{\lambda}) + \frac{\lambda}{2}\|\hat{\boldsymbol{w}}_{\lambda}\|^{2} \;\ge\; \frac{\lambda}{2}\|\hat{\boldsymbol{w}}_{\lambda}\|^{2} .$$

The last inequality uses $L_{\lambda}(\hat{\boldsymbol{w}}_{\lambda}) \ge 0$, which follows from part (a) of Theorem 6.3 together with $\frac{\lambda}{2}\|\hat{\boldsymbol{w}}_\lambda\|^2 \ge 0$. Rearranging, $\|\hat{\boldsymbol{w}}_{\lambda}\|^{2} \le 2n\log 2/\lambda$, that is $\|\hat{\boldsymbol{w}}_{\lambda}\| \le \sqrt{2n\log 2/\lambda}$.

That the weights stay within this range even for linearly separable data is quantitative evidence that regularisation really does stop the divergence of Theorem 6.3. The bound tending to $\infty$ as $\lambda \to 0$ is likewise consistent with the unregularised situation.

References

• C. M. Bishop, Pattern Recognition and Machine Learning, Springer, 2006 — Chapter 4, “Linear Models for Classification”. §4.2 contains the derivation from a generative model given in Example 3.4, and §4.3 treats logistic regression and IRLS.
• T. Hastie, R. Tibshirani, J. Friedman, The Elements of Statistical Learning, 2nd ed., Springer, 2009 — Chapter 4, “Linear Methods for Classification”, including the unboundedness in the linearly separable case and the treatment of regularisation. Version made available by the authors.
• S. Boyd, L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004 — Chapter 3 (convex functions and strong convexity) and Chapter 7 (maximum likelihood estimation as convex optimisation). Version made available by the authors.
• I. Goodfellow, Y. Bengio, A. Courville, Deep Learning, MIT Press, 2016 — Chapter 6, on the correspondence between the choice of output unit and the cross entropy loss, and on the gradient saturation discussed in Example 5.3. Public version.
• J. Berkson, “Application of the Logistic Function to Bio-Assay”, Journal of the American Statistical Association 39 (1944), 357–365 — the paper introducing the word “logit”.
• J. A. Nelder, R. W. M. Wedderburn, “Generalized Linear Models”, Journal of the Royal Statistical Society, Series A 135 (1972), 370–384 — the paper placing logistic regression within the framework of generalised linear models.
• Takuya Kubo, Data Kaiseki no tame no Tokei Modeling Nyumon (Introduction to Statistical Modelling for Data Analysis), Iwanami Shoten, 2012 (in Japanese) — a chapter treating logistic regression from the standpoint of generalised linear models. An accessible introduction from the practical side of statistical modelling.

Appendix: generalisation to several classes, and Newton’s method

Softmax regression. When there are $K$ classes, provide a weight vector $\boldsymbol{w}_k \in \mathbb{R}^{d}$ for each class and set

$$P(y = k \mid \boldsymbol{x}) = \frac{\exp(\boldsymbol{w}_k^{\mathsf{T}}\boldsymbol{x})}{\sum_{l=1}^{K}\exp(\boldsymbol{w}_l^{\mathsf{T}}\boldsymbol{x})} .$$

This is the softmax function. For $K = 2$, dividing numerator and denominator by $\exp(\boldsymbol{w}_0^{\mathsf{T}}\boldsymbol{x})$ gives $P(y=1\mid\boldsymbol{x}) = \sigma\bigl((\boldsymbol{w}_1-\boldsymbol{w}_0)^{\mathsf{T}}\boldsymbol{x}\bigr)$, recovering the sigmoid (only differences of weights matter, so the parameters carry only $K-1$ blocks’ worth of freedom). Writing the label as a one-hot vector $\boldsymbol{t}_i$ (with a $1$ in the $y_i$-th component only), the loss is $L = -\sum_i \sum_k t_{ik}\log \mu_{ik}$, again a cross entropy. The gradient has the same shape as in Theorem 5.1:

$$\frac{\partial L}{\partial \boldsymbol{w}_k} = \sum_{i=1}^{n} (\mu_{ik} - t_{ik})\,\boldsymbol{x}_i$$

(see the gradient of softmax with cross entropy（Proposition 7.1）［ニューラルネットワークと逆伝播］). The point is that the structure producing the cancellation is preserved intact.

Newton’s method and IRLS. Since Theorem 6.1 gave us the Hessian as well, we can use second-order information rather than the gradient alone. The Newton update

$$\boldsymbol{w} \leftarrow \boldsymbol{w} - \bigl(X^{\mathsf{T}}SX\bigr)^{-1}X^{\mathsf{T}}(\boldsymbol{\mu}-\boldsymbol{y})$$

can, after rearranging the right-hand side, be rewritten as $\boldsymbol{w} \leftarrow (X^{\mathsf{T}}SX)^{-1}X^{\mathsf{T}}S\boldsymbol{z}$ with $\boldsymbol{z} = X\boldsymbol{w} + S^{-1}(\boldsymbol{y}-\boldsymbol{\mu})$, which is the form of a weighted least squares problem. Because the weights $S$ are updated at every iteration, this is called IRLS (iteratively reweighted least squares). It converges quickly, but each iteration handles the inverse of a $d\times d$ matrix (in practice, a system of linear equations), which becomes heavy for large $d$. That difference in cost is one reason why first-order gradient descent is used in deep learning.