Logistic Regression — Theory

0. Notation

All symbols used below — defined once.

Symbol	Type	Meaning
\(n\)	scalar	number of training examples
\(p\)	scalar	number of features (including bias if used)
\(x_i\)	vector of length \(p\)	feature vector of the \(i\)-th example (a column of \(X^T\))
\(X\)	matrix of size \(n \times p\)	design matrix; row \(i\) is \(x_i^T\)
\(y_i\)	scalar in \(\lbrace0, 1\rbrace\)	binary label of the \(i\)-th example
\(y\)	vector of length \(n\)	stacked labels
\(\theta\)	vector of length \(p\)	model parameters
\(z_i\)	scalar	linear score for example \(i\): \(z_i := x_i^T \theta\)
\(\sigma(z)\)	scalar in \((0, 1)\)	sigmoid function: \(\sigma(z) := 1 / (1 + e^{-z})\)
\(p_i\)	scalar in \((0, 1)\)	model probability \(P(y_i = 1 \mid x_i, \theta) = \sigma(z_i)\)
\(p\)	vector of length \(n\)	stacked probabilities \((p_1, \dots, p_n)\); distinguished from the feature count \(p\) by context and type
\(L(\theta)\)	scalar	the cross-entropy / NLL loss; defined in §3
\(\lambda\)	scalar \(\ge 0\)	regularisation strength
\(I_p\)	matrix of size \(p \times p\)	identity matrix
\(W\)	diagonal matrix of size \(n \times n\)	weights \(W_{ii} := p_i (1 - p_i)\) used in the Hessian (§5)
\(K\)	scalar	number of classes in the multi-class extension (§11)

Bias convention. We prepend a constant column of ones to \(X\) so the intercept is absorbed into \(\theta_0\). Penalties (when added in §10) act on \(\theta_1, \dots, \theta_{p-1}\) only, not on the intercept.

Vector convention. All vectors are column vectors. Lowercase Latin / Greek = vector; uppercase = matrix.

Norms carry explicit subscripts: \(\Vert\theta\Vert_1\), \(\Vert\theta\Vert_2^2\).

1. WHY — From Linear Scores to Probabilities

A linear score \(z = \theta^T x\) takes values on \((-\infty, +\infty)\). When the target is a binary label \(y \in \lbrace0, 1\rbrace\), using \(z\) directly as a prediction fails:

Out-of-range predictions. A prediction of \(\hat{y} = 1.4\) has no meaning as a probability — there is no "140 percent chance."
Outlier sensitivity. Observations far from the boundary distort the fitted line's slope, dragging the decision threshold in the wrong direction.
Bad probability model. The truth is closer to: low \(x\) → near-certain class 0 (probability near 0), high \(x\) → near-certain class 1 (probability near 1), with a smooth transition in between. A straight line cannot capture this saturation.

The fix: take the linear score and squash it through a function that maps \((-\infty, +\infty) \to (0, 1)\). That function is the sigmoid.

1.1 Assumptions

Bernoulli response. Each \(y_i \in \lbrace0, 1\rbrace\) is Bernoulli with success probability \(p_i = \sigma(x_i^T \theta)\).
Conditional independence. Given the features and \(\theta\), the labels \(y_1, \dots, y_n\) are independent (equivalently: the pairs \((x_i, y_i)\) are i.i.d.).
Linear log-odds. The log-odds \(\log(p_i / (1 - p_i)) = x_i^T \theta\) is a linear function of features.
No perfect multicollinearity. \(X\) has full column rank (needed for the MLE — when it exists — to be unique, and for finite standard errors).
No separation. No hyperplane separates the classes completely or quasi-completely; otherwise the MLE does not exist (Albert & Anderson, 1984; see §5.3).

2. The Sigmoid Function and Its Identities

2.1 Definition

The sigmoid (also called the logistic function) is

\[\sigma(z) := \frac{1}{1 + e^{-z}}, \qquad z \in \mathbb{R}.\]

Range \((0, 1)\), strictly increasing, infinitely differentiable, \(\sigma(0) = 1/2\).

2.2 Identities

All three follow from direct calculation:

(I1) Reflection. \(\sigma(-z) = 1 - \sigma(z)\).

(I2) Derivative. \(\sigma'(z) = \sigma(z) (1 - \sigma(z))\).

(I3) Logit inverse. If \(p = \sigma(z)\) with \(p \in (0, 1)\), then \(z = \log(p / (1 - p))\).

Proof of (I2). Write \(\sigma(z) = (1 + e^{-z})^{-1}\) and apply the chain rule:

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

\(\blacksquare\)

Reading. Identity (I2) is the algebraic reason logistic regression has clean closed-form gradients and Hessian — the derivative of the sigmoid is itself expressible as a product of sigmoids. Identity (I3) is what lets us interpret \(\theta^T x\) as a log-odds (§8).

3. Likelihood and Cross-Entropy Loss

3.1 Model

Conditional on \(x_i\) and \(\theta\), model \(y_i\) as a Bernoulli random variable with success probability \(p_i := \sigma(x_i^T \theta)\):

\[y_i \mid x_i, \theta \sim \text{Bernoulli}(p_i), \qquad i = 1, \dots, n.\]

Assume training examples are i.i.d. given \(\theta\).

3.2 Likelihood

The probability mass function of a Bernoulli is \(P(y_i \mid p_i) = p_i^{y_i} (1 - p_i)^{1 - y_i}\). The full data likelihood is the product:

\[\mathcal{L}(\theta) = \prod_{i=1}^n p_i^{y_i} (1 - p_i)^{1 - y_i}.\]

The exponents act as a switch: example \(i\) contributes the factor \(p_i\) when \(y_i = 1\), and the factor \(1 - p_i\) when \(y_i = 0\).

3.3 Negative log-likelihood — the cross-entropy loss

Take the negative log of the likelihood and divide by \(n\):

\[L(\theta) := -\frac{1}{n} \log \mathcal{L}(\theta) = -\frac{1}{n} \sum_{i=1}^n \Big[y_i \log p_i + (1-y_i)\log(1-p_i)\Big]. \qquad (3.1)\]

This is the cross-entropy loss for binary classification (also called the binary log-loss or logistic loss).

Three key observations:

Non-negative. \(p_i \in (0, 1)\) implies \(\log p_i \le 0\) and \(\log(1 - p_i) \le 0\). So \(L(\theta) \ge 0\).
Asymmetric per example. For \(y_i = 1\) only \(\log p_i\) appears; for \(y_i = 0\) only \(\log(1 - p_i)\). The loss for one example is large exactly when the model puts low probability on the true class.
Blows up on confident wrong predictions. If \(y_i = 1\) but \(p_i \to 0\), then \(-\log p_i \to +\infty\): the per-example loss is unbounded (the sigmoid saturates; the loss does not).

Cross-entropy aggressively punishes confident-but-wrong predictions — precisely why softmax + cross-entropy is the default loss in deep learning.

3.4 Rewriting the loss with \(z_i\) — the margin form

Replace \(p_i = \sigma(z_i)\) in (3.1) and use \(1 - \sigma(z) = \sigma(-z)\) from (I1):

\[L(\theta) = -\frac{1}{n} \sum_{i=1}^n \Big[ y_i \log \sigma(z_i) + (1 - y_i) \log \sigma(-z_i) \Big].\]

Use the identity \(\log \sigma(z) = -\log(1 + e^{-z})\). The single-example loss simplifies into a clean unified form:

\[ \ell_i(\theta) = \log(1 + e^{-y_i' z_i}) \quad \text{where } y_i' := 2 y_i - 1 \in \{-1, +1\}. \qquad (3.2) \]

Proof. For \(y_i = 1\), so that \(y_i' = +1\):

\[ -\log \sigma(z_i) = \log(1 + e^{-z_i}) = \log(1 + e^{-y_i' z_i}) \]

For \(y_i = 0\), so that \(y_i' = -1\):

\[ -\log \sigma(-z_i) = \log(1 + e^{z_i}) = \log(1 + e^{-y_i' z_i}) \]

\(\blacksquare\)

Result: The cross-entropy loss is a smooth, monotone function of the margin \(y_i' z_i\). A large positive margin (correct sign, high magnitude) → small loss. A large negative margin (wrong sign) → large loss.

This is the same margin concept that drives SVMs, but with the smooth log-loss instead of the kinked hinge loss.

4. Gradient of the Cross-Entropy Loss

4.1 Theorem (gradient)

Theorem 4.1. The gradient of the cross-entropy loss (3.1) is

\[ \nabla L(\theta) = \frac{1}{n} \sum_{i=1}^n (p_i - y_i) \, x_i = \frac{1}{n} X^T (p - y), \qquad (4.1) \]

where \(p \in \mathbb{R}^n\) is the vector with entries \(p_i = \sigma(x_i^T \theta)\).

Proof. Differentiate the single-example loss with respect to \(\theta\). From (3.1):

\[\frac{\partial \ell_i}{\partial \theta} = -\Big[ y_i \cdot \frac{\sigma'(z_i)}{\sigma(z_i)} - (1 - y_i) \cdot \frac{\sigma'(z_i)}{1 - \sigma(z_i)} \Big] \cdot \frac{\partial z_i}{\partial \theta}.\]

Use (I2), \(\sigma'(z_i) = \sigma(z_i)(1 - \sigma(z_i)) = p_i(1 - p_i)\), and \(\partial z_i / \partial \theta = x_i\):

\[\frac{\partial \ell_i}{\partial \theta} = -\Big[ y_i (1 - p_i) - (1 - y_i) p_i \Big] x_i = -\Big[ y_i - p_i \Big] x_i = (p_i - y_i) x_i.\]

Average over \(i\) to get (4.1). \(\blacksquare\)

Result: \(\nabla L(\theta) = \frac{1}{n} X^T (p - y)\).

4.2 Comparison with OLS

Model	Gradient	Residual
Linear regression	\((1/n) \cdot X^T r\)	\(r = X\theta - y\)
Logistic regression	\((1/n) \cdot X^T r\)	\(r = p - y\), with \(p_i = \sigma(x_i^T \theta)\)

The structure is identical — sum over examples of (residual) \(\cdot\) (feature). The only change is that the residual is now \(p_i - y_i\) (predicted probability minus label). This is why every linear-regression algorithm (GD, SGD, coordinate descent) generalises straight to logistic regression with almost no change.

Convention note. The linear-regression row uses the loss \(\frac{1}{2n}\Vert X\theta - y\Vert_2^2\); with the plain MSE \(\frac{1}{n}\Vert X\theta - y\Vert_2^2\) the gradient is \((2/n) X^T r\) instead — a cosmetic constant that does not change the minimiser.

4.3 First-order condition

Setting \(\nabla L(\theta) = 0\) gives

\[X^T (p - y) = 0. \qquad (4.2)\]

This says the residual vector \(p - y\) must be orthogonal to every column of \(X\) — the same geometric condition as OLS, but for the probability residual. The crucial difference: \(p\) depends on \(\theta\) through the sigmoid, so (4.2) is a non-linear equation in \(\theta\).

5. Hessian and Convexity

5.1 Theorem (Hessian)

Theorem 5.1. The Hessian of the cross-entropy loss (3.1) is

\[\nabla^2 L(\theta) = \frac{1}{n} X^T W X, \qquad (5.1)\]

where \(W \in \mathbb{R}^{n \times n}\) is the diagonal matrix with entries \(W_{ii} = p_i (1 - p_i)\).

Proof. Differentiate (4.1) once more. Each gradient entry is \(\frac{1}{n} \sum_i (p_i - y_i) x_{ij}\).

The derivative of \(p_i\) with respect to \(\theta_k\) is \(\sigma'(z_i) \cdot x_{ik} = p_i(1 - p_i) x_{ik}\) by (I2). So

\[\frac{\partial^2 L}{\partial \theta_j \partial \theta_k} = \frac{1}{n} \sum_i p_i(1 - p_i) x_{ij} x_{ik}.\]

In matrix form this is exactly \((1/n) \cdot X^T W X\). \(\blacksquare\)

Result: \(\nabla^2 L(\theta) = \frac{1}{n} X^T W X\) with \(W_{ii} = p_i(1 - p_i)\).

5.2 Theorem (convexity)

Theorem 5.2. The cross-entropy loss \(L(\theta)\) is convex on \(\mathbb{R}^p\). If \(X\) has full column rank, then \(L\) is strictly convex. (No condition on \(W\) is needed: \(W_{ii} = p_i(1 - p_i) > 0\) automatically, since \(\sigma\) never reaches \(0\) or \(1\) at finite \(\theta\).)

Proof. A twice-differentiable function is convex iff its Hessian is positive semidefinite everywhere. For any vector \(v \in \mathbb{R}^p\):

\[v^T \nabla^2 L(\theta) v = \frac{1}{n} v^T X^T W X v = \frac{1}{n} (X v)^T W (X v) = \frac{1}{n} \sum_i W_{ii} (X v)_i^2.\]

Each \(W_{ii} = p_i (1 - p_i) > 0\) for \(p_i \in (0, 1)\), and \((X v)_i^2 \ge 0\). So the sum is \(\ge 0\) → positive semidefinite → convex.

If \(X\) has full column rank, then \(Xv = 0\) implies \(v = 0\), so the sum is strictly positive for \(v \ne 0\) → positive definite → strictly convex. \(\blacksquare\)

Result: The cross-entropy loss is convex — strictly convex when \(X\) has full column rank — so it has at most one global minimiser and no spurious local minima. This is what separates logistic regression from neural networks (non-convex, multiple local minima). Caution: strict convexity guarantees uniqueness, not existence — on separable data no minimiser exists (§5.3).

5.3 Failure case — perfect separability

If the two classes are perfectly linearly separable, then \(L(\theta)\) has no finite minimiser: shrinking the loss toward zero requires \(\Vert\theta\Vert \to \infty\) to push the sigmoid to 0/1 saturation on every example. The MLE does not exist.

The same failure occurs under quasi-complete separation (a separating hyperplane with some points lying exactly on it); Albert & Anderson (1984) show the MLE exists iff neither kind of separation is present.

Cure: Add a coercive penalty such as the ridge term \(\lambda \Vert\theta\Vert_2^2\) (§10). The penalised loss then tends to \(+\infty\) as \(\Vert\theta\Vert_2 \to \infty\), so a finite minimiser exists even on perfectly separable data — and it is unique, since the ridge term is strictly convex.

(Strict convexity alone would not be enough: it guarantees uniqueness, not existence; coercivity is what restores existence.)

6. No Closed-Form Solution

For OLS we had \(\hat{\theta}_{\text{OLS}} = (X^T X)^{-1} X^T y\) in closed form. Why does this not work for logistic regression?

The first-order condition (4.2) is \(X^T (p - y) = 0\), which expands to

\[\sum_{i=1}^n \Big( \sigma(x_i^T \theta) - y_i \Big) x_i = 0.\]

Each term contains \(\sigma(x_i^T \theta)\) — a non-linear function of \(\theta\). There is no way to factor \(\theta\) out and isolate it. The equation is transcendental.

We are forced to solve iteratively. The good news (Theorem 5.2): the loss is convex, so there are no bad local minima — provided a minimiser exists (no separation, §5.3) and step sizes are chosen appropriately (§9), standard descent methods converge to the global optimum from any starting point.

7. Decision Boundary Geometry

In \(d\) dimensions, the logistic model predicts

\[P(y = 1 \mid x) = \sigma(\theta_0 + \theta_1 x_1 + \cdots + \theta_d x_d).\]

The decision boundary is the set of points where the model is on the fence — \(P(y = 1 \mid x) = 0.5\) — which (since \(\sigma(0) = 1/2\) and \(\sigma\) is strictly increasing) is exactly

\[\theta_0 + \theta_1 x_1 + \cdots + \theta_d x_d = 0.\]

This is a hyperplane. So logistic regression is a linear classifier, even though its probability output is a non-linear S-curve.

What varies across space is how confident the model is — not the shape of the boundary:

Far from the hyperplane, on the "1" side, probability tends to 1.
Far from the hyperplane, on the "0" side, probability tends to 0.
Right on the hyperplane, probability is 0.5.
The transition band where the probability swings from 0.1 to 0.9 has a width proportional to \(1 / \Vert w\Vert_2\), where \(w := (\theta_1, \dots, \theta_d)\) is the weight part of \(\theta\).
The intercept \(\theta_0\) shifts the boundary but not its sharpness.
A larger \(\Vert w\Vert_2\) means a sharper transition.

Practical consequence. Logistic regression cannot separate XOR-shaped data with a single boundary. When a non-linear classifier is needed, either engineer non-linear features (polynomial trick) or use a non-linear model (decision trees, neural networks).

8. Log-Odds Interpretation

8.1 Odds and the logit

Definition (odds). The odds of an event with probability \(p\) is \(\text{odds}(p) := p / (1 - p)\).

\(p = 0.5\) → odds = 1 ("50:50").
\(p = 0.9\) → odds = 9 ("9 to 1 in favour").
\(p = 0.1\) → odds \(\approx\) 0.111 ("1 to 9").

The function \(p \mapsto \log(p / (1 - p))\) is called the logit — and by (I3) it is the inverse of the sigmoid. So if \(p = \sigma(z)\), then \(z = \log(p / (1 - p))\).

8.2 Coefficient interpretation

The model's linear score \(z = \theta^T x\) is the log-odds of class 1. Interpretation by coefficient:

Increasing feature \(x_j\) by 1 unit (other features fixed) multiplies the odds of class 1 by \(e^{\theta_j}\).
\(\theta_j > 0\) → feature is associated with class 1.
\(\theta_j < 0\) → feature is associated with class 0.
\(|\theta_j|\) large → feature has a strong effect.

This is exactly the kind of interpretability that makes logistic regression a workhorse in medicine, social science, and credit risk: the coefficients have a direct meaning on the odds-ratio scale.

9. Optimisation Methods

Since no closed-form solution exists (§6), we must solve the first-order condition iteratively. Four methods, in increasing sophistication.

9.1 Gradient descent

Update: \(\theta_{k+1} \leftarrow \theta_k - \eta \cdot \nabla L(\theta_k)\) with \(\nabla L(\theta_k) = (1/n) X^T(p_k - y)\).

Step size. The Lipschitz constant of \(\nabla L\) is bounded by

\[ L_{\text{smooth}} \le \frac{1}{4 n} \, \lambda_{\max}(X^T X), \]

because \(W_{ii} = p_i (1 - p_i) \le 1/4\) for every \(p_i \in (0, 1)\) (maximum at \(p_i = 1/2\)). With step \(\eta = 1 / L_{\text{smooth}}\) (and assuming a minimiser \(\theta^\ast\) exists — no separation, §5.3), convergence is

\[ L(\theta_k) - L(\theta^\ast) \le \frac{\|\theta_0 - \theta^\ast\|_2^2}{2 \eta \cdot k} = O(1 / k). \]

Result: Sub-linear convergence — simple but slow.

9.2 Newton's method — IRLS

Newton's method replaces the constant step direction \(-g_k\) by \(-(\nabla^2 L(\theta_k))^{-1} g_k\), using local curvature:

\[\theta_{k+1} = \theta_k - \Big(\frac{1}{n} X^T W_k X\Big)^{-1} \cdot \frac{1}{n} X^T(p_k - y).\]

Why IRLS. Define the "working response" \(z_k := X \theta_k + W_k^{-1} (y - p_k)\). The Newton update is equivalent to the weighted least squares solve:

\[\theta_{k+1} = (X^T W_k X)^{-1} X^T W_k z_k.\]

At each step we solve a weighted OLS problem with weights \(W_k\) and response \(z_k\). The weights depend on \(\theta_k\), so we re-weight and re-solve at every iteration — hence iteratively reweighted least squares (IRLS).

Convergence. Locally quadratic — once \(\theta_k\) is close enough to \(\theta^\ast\):

\[ \Vert\theta_{k+1} - \theta^\ast\Vert_2 \le C \cdot \Vert\theta_k - \theta^\ast\Vert_2^2 \]

So the number of correct digits roughly doubles per iteration. (Far from the optimum, a damped / line-search Newton step is used to guarantee progress.) 5–15 iterations reach machine precision on typical problems.

Cost per step. \(\Theta(n p^2 + p^3)\) — one \(p \times p\) matrix solve per iteration. Fine when \(p \lesssim 10^3\).

Result: Newton/IRLS converges quadratically — vastly faster than GD per iteration, at the cost of \(O(p^3)\) per step.

9.3 Stochastic gradient descent (SGD)

Compute the gradient on a random mini-batch of size \(b \ll n\), paying \(O(b p)\) per step instead of \(O(n p)\). With a decaying step size \(\eta_k \propto 1/\sqrt{k}\), SGD converges at rate \(O(1/\sqrt{k})\) in expectation (for the averaged iterate) in the convex case.

When to use: \(n\) very large (millions of examples), online learning (data arrives in a stream), or deep learning (logistic regression is the smallest such network).

9.4 L-BFGS — the practical default

L-BFGS (Limited-memory BFGS) is a quasi-Newton method that approximates the inverse Hessian using only the last \(m\) gradient differences (typically \(m = 10\text{–}20\)). It gets:

Fast convergence in practice (full-memory BFGS is provably super-linear; L-BFGS itself is guaranteed only linear, yet is typically far faster than GD).
GD-like memory cost (\(O(mp)\), no \(p \times p\) Hessian).
No second derivatives needed.

This is the practical default. scikit-learn's LogisticRegression(solver="lbfgs") uses L-BFGS.

Result: Use L-BFGS in practice; understand Newton/IRLS conceptually.

10. Regularised Logistic Regression

10.1 The penalised loss

Same recipe as Ridge/Lasso for regression: add a penalty on \(\theta\) (excluding the intercept):

\[ L_{\text{ridge}}(\theta) := L(\theta) + \lambda \|\theta\|_2^2, \qquad L_{\text{lasso}}(\theta) := L(\theta) + \lambda \|\theta\|_1. \qquad (10.1) \]

(Per §0, the norms in (10.1) run over \(\theta_1, \dots, \theta_{p-1}\) only — the intercept \(\theta_0\) is never penalised.)

Both are convex (sum of two convex functions). The L2-penalised objective is strictly convex and coercive, so it has a unique finite minimiser even on linearly separable data — the regulariser fixes the failure case of §5.3.

10.2 Bayesian interpretation

Theorem 10.2. Under the Bernoulli likelihood from §3 and the prior

\(\theta_j \sim \text{Normal}(0, \tau^2)\) independently → MAP equals L2-regularised logistic regression with \(\lambda = 1 / (2 n \tau^2)\).

\(\theta_j \sim \text{Laplace}(0, b)\) independently → MAP equals L1-regularised logistic regression with \(\lambda = 1 / (n b)\).

Proof sketch. Bayes' rule:

\[ \log p(\theta \mid y, X) = \log \mathcal{L}(\theta) + \log p(\theta) + \text{const} \]

Negate and divide by \(n\) — the data term is \(L(\theta)\), and the prior adds \(\frac{1}{2 n \tau^2} \Vert\theta\Vert_2^2\) (Normal) or \(\frac{1}{n b} \Vert\theta\Vert_1\) (Laplace). Identifying the multiplier with \(\lambda\) gives the claim. \(\blacksquare\)

10.3 Gradient of the regularised loss

L2. \(\nabla L_{\text{ridge}}(\theta) = \frac{1}{n} X^T (p - y) + 2 \lambda \theta\).
L1. \(L_{\text{lasso}}\) is not differentiable where any \(\theta_j = 0\); work with the subdifferential \(\partial L_{\text{lasso}}(\theta) = \big\lbrace \frac{1}{n} X^T (p - y) + \lambda s \thinspace:\thinspace s \in \partial \Vert\theta\Vert_1 \big\rbrace\).
Subgradient entries. Here \(s_j = \operatorname{sign}(\theta_j)\) if \(\theta_j \ne 0\) and \(s_j \in [-1, 1]\) if \(\theta_j = 0\).

10.4 Algorithmic changes

For L2: gradient and Hessian pick up extra terms:

\[\nabla^2 L_{\text{ridge}}(\theta) = \frac{1}{n} X^T W X + 2 \lambda I_p.\]

The added \(2\lambda I_p\) makes the system strictly positive definite even when \(X\) is rank-deficient — numerically stable Newton solve.

For L1: use the proximal-gradient (ISTA) recipe — alternate a gradient step on the smooth data term with a soft-thresholding step for the L1 penalty. FISTA acceleration applies unchanged.

Result: L1 logistic regression produces sparse coefficients (variable selection for classification). L2 logistic regression shrinks coefficients smoothly and is the default in most software.

In scikit-learn's LogisticRegression, C is the inverse regularisation strength (larger C = weaker penalty); it minimises \(\frac{1}{2}\Vert w\Vert_2^2 + C \sum_i \log(1 + e^{-y_i' z_i})\), so under the conventions of (10.1) the exact correspondence is \(C = 1/(2 n \lambda)\).

11. Multi-Class Extension: Softmax Regression

For \(K \ge 3\) classes, replace the sigmoid by the softmax function and the Bernoulli by the categorical distribution.

11.1 Definition (softmax)

Given a vector of scores \(z \in \mathbb{R}^K\):

\[\text{softmax}(z)_k = \frac{e^{z_k}}{\sum_{j=1}^K e^{z_j}}, \qquad k = 1, \dots, K.\]

By construction \(\text{softmax}(z)_k \in (0, 1)\) and \(\sum_k \text{softmax}(z)_k = 1\). When \(K = 2\), the softmax reduces to the sigmoid (up to re-parameterisation \(z = z_1 - z_2\)).

11.2 Model

One parameter vector \(\theta_k \in \mathbb{R}^p\) per class. Stack as \(\Theta \in \mathbb{R}^{K \times p}\):

\[P(y_i = k \mid x_i, \Theta) = \text{softmax}(\Theta x_i)_k = \frac{e^{\theta_k^T x_i}}{\sum_{j=1}^K e^{\theta_j^T x_i}}.\]

Each class carries its own parameter vector, and the softmax turns the \(K\) scores \(\theta_k^T x_i\) into probabilities that sum to one across classes.

11.3 Loss

Encode each label as a one-hot vector \(y_i \in \lbrace0, 1\rbrace^K\). The multinomial cross-entropy loss is:

\[L(\Theta) = -\frac{1}{n} \sum_{i=1}^n \sum_{k=1}^K y_{ik} \log P(y_i = k \mid x_i, \Theta). \qquad (11.1)\]

One-hot encoding kills every term but the one for the true class, so the loss is the average negative log-probability assigned to the correct label.

11.4 Gradient

Same shape as Theorem 4.1, per class:

\[ \frac{\partial L}{\partial \theta_k} = \frac{1}{n} \sum_{i=1}^n (P(y_i = k \mid x_i, \Theta) - y_{ik}) \, x_i. \qquad (11.2) \]

Result: The multinomial gradient has exactly the binary structure \((1/n) X^T(P - Y)\), repeated per class.

11.5 Identifiability

Softmax has a redundancy: shifting every \(\theta_k\) by the same vector \(c\) leaves the probabilities unchanged. So \(\Theta\) is identifiable only up to a constant shift — fixed by either (i) setting one class as reference (\(\theta_K = 0\)), or (ii) adding an L2 penalty.

12. Statistical Properties

12.1 Consistency

Everything in §12.1–12.3 assumes the model is correctly specified: there exists a true \(\theta^\ast\) with \(P(y = 1 \mid x) = \sigma(x^T \theta^\ast)\). Under this and mild regularity (\(X\) full rank with high probability, \(E[\Vert x_i\Vert^4] < \infty\)), the MLE is consistent: \(\hat{\theta}_n \to \theta^\ast\) in probability as \(n \to \infty\).

12.2 Asymptotic normality

\[\sqrt{n} (\hat{\theta}_n - \theta^\ast) \xrightarrow{d} \text{Normal}(0, \mathcal{I}(\theta^\ast)^{-1}),\]

where \(\mathcal{I}(\theta^\ast)\) is the Fisher information matrix:

\[ \mathcal{I}(\theta^\ast) := \mathbb{E}_x [\sigma(x^T \theta^\ast) (1 - \sigma(x^T \theta^\ast)) \, x x^T]. \]

The empirical analogue is \((1/n) X^T W X\) — the Hessian from Theorem 5.1.

12.3 Standard errors and confidence intervals

\[\widehat{\text{Var}}(\hat{\theta}_n) \approx (X^T \widehat{W} X)^{-1},\]

where \(\widehat{W}\) uses plug-in probabilities \(\hat{p}_i = \sigma(x_i^T \hat{\theta}_n)\). The square root of the \(j\)-th diagonal entry gives \(\widehat{SE}(\hat{\theta}_j)\); a 95% CI is \(\hat{\theta}_j \pm 1.96 \cdot \widehat{SE}(\hat{\theta}_j)\).

12.4 Threshold selection

The model outputs \(p(x) = P(y = 1 \mid x)\). To produce a class label, pick threshold \(\tau\):

0-1 loss. If \(p(x)\) is the true conditional probability, the Bayes-optimal threshold is \(\tau = 0.5\).
Asymmetric costs. If false positive costs \(c_{10}\) and false negative costs \(c_{01}\):

\[\tau^\ast = \frac{c_{10}}{c_{10} + c_{01}}.\]

A costlier false positive pushes \(\tau^\ast\) up, so the model must be more confident before it predicts class 1; equal costs give back the 0-1 threshold \(\tau = 0.5\).

12.5 Calibration

A classifier is calibrated if \(P(y = 1 \mid p(x) = q) = q\) for every \(q \in (0, 1)\). Logistic regression is usually well-calibrated when the model is correctly specified — because its loss is the cross-entropy (Bernoulli likelihood).

12.6 ROC and AUC

Sweep \(\tau\) from 1 to 0. At each \(\tau\): TPR = \(P(\hat{y}=1 \mid y=1)\), FPR = \(P(\hat{y}=1 \mid y=0)\). The ROC curve plots TPR vs FPR.

AUC (Area Under the Curve) is threshold-free and equals \(P(p(x_+) > p(x_-)) + \tfrac{1}{2} P(p(x_+) = p(x_-))\) — the probability that the model ranks a random positive above a random negative (ties count half).

13. Connections

Linear regression. Same gradient structure \((1/n) X^T r\); residual changes from \(X\theta - y\) to \(\sigma(X\theta) - y\).
Ridge / Lasso regression. Same penalty terms, same Bayesian interpretation; only the data term changes from squared error to cross-entropy.
Smallest neural network. A logistic regression model is exactly one linear layer → sigmoid activation → output. Stack two and you have a 2-layer MLP.
Training algorithm. The training algorithm (gradient descent on cross-entropy) is the same algorithm used to train every modern deep network.
Support vector machines. Same margin concept (\(y' z\)), but with the smooth log-loss \(\log(1 + e^{-y' z})\) instead of the kinked hinge loss \(\max(0, 1 - y' z)\).
Generalised linear models (GLM). IRLS generalises to any exponential family (Poisson regression, gamma regression) by varying \(W_k\) and \(z_k\).
Probability
Gradient Descent
Loss Functions

14. Failure Modes

Non-linear decision boundary. Logistic regression's boundary is a hyperplane (§7). XOR-shaped data, images, high-order interactions → needs feature engineering or non-linear models.
Perfect separability. Unregularised MLE doesn't exist (§5.3). Cure: L2 regularisation.
Many irrelevant features. Variance scales with \(p/n\). Cure: L1 regularisation for variable selection.
Class imbalance. Default threshold 0.5 predicts majority class. Cure: tune threshold from cost ratio (§12.4), or use class weighting.
Model mis-specification. When the true decision function is far from linear, probabilities can be over-confident. Cure: post-hoc calibration or a more flexible model.

15. References

Hastie, T., Tibshirani, R., & Friedman, J. (2009). The Elements of Statistical Learning: Data Mining, Inference, and Prediction (2nd ed.). Springer. Chapter 4.4: Logistic Regression.
Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer. Chapter 4.3: Probabilistic Discriminative Models.
McCullagh, P., & Nelder, J. A. (1989). Generalized Linear Models (2nd ed.). Chapman and Hall/CRC.
Albert, A., & Anderson, J. A. (1984). On the existence of maximum likelihood estimates in logistic regression models. Biometrika, 71(1), 1–10. (Shows the MLE exists iff the data exhibit neither complete nor quasi-complete separation — the theory behind §5.3.)