Logistic Regression

Logistic regression is a linear classifier that outputs a probability. It computes a weighted sum of an image’s features, adds a bias, and passes the result through the logistic sigmoid, so that the score becomes a number between 0 and 1: the predicted probability that the image, a patch or a pixel belongs to the positive class. Despite its name it is a classification method, not a regression in the machine-learning sense (Bishop, 2006, §4.3.2). It is the standard baseline for machine learning on fixed features, the model a linear probe fits on top of a frozen network, and, with a softmax in place of the sigmoid, the output layer of most image classification networks.

The model is fitted by maximum likelihood: its negative log-likelihood is the binary cross-entropy loss, which is convex in the weights, so a global minimum is found by Newton’s method or by gradient descent.

The logistic function itself is older than the classifier. According to Cramer’s (2002) history, Pierre-François Verhulst used it for population growth in three papers from 1838 to 1847 and named it the logistic in the second of them, published in 1845. In Cramer’s account, Joseph Berkson proposed it in 1944 as an alternative to the probit for bio-assay, the analysis of how the probability of a response depends on a dose, and coined the term logit by analogy with the probit; Wilson and Worcester (1943) were, Cramer judges, probably the first to publish an application in bio-assay, just before Berkson. Cramer names Cox, in articles of the 1960s and a 1969 book, as one of the first to recognize its wider uses, and dates the multinomial generalization to Cox (1966) and, independently, Theil (1969).

Problem

Given images described by feature vectors $\mathbf{x} \in \mathbb{R}^d$ and binary labels $y \in \{0, 1\}$, estimate the probability $P(y = 1 \mid \mathbf{x})$ for a new image, and with it a decision. The features can be designed (the gradient energy of a patch, a HOG descriptor), the raw pixel values, or the embedding a pretrained network computes. This article uses labels in $\{0, 1\}$, as Hastie, Tibshirani and Friedman (2009, §4.4.1), Bishop (§4.3.2) and scikit-learn do; Murphy (2022, §10.2) also gives the form for $\tilde{y} \in \{-1, +1\}$, the convention of support vector machines.

Inputs and Outputs

• Input: a training set of $n$ images, stored as a data matrix $X \in \mathbb{R}^{n \times d}$ with one feature vector per row and a label vector $\mathbf{y} \in \{0, 1\}^n$; optionally a regularization strength and class or sample weights.
• Output: a weight vector $\mathbf{w} \in \mathbb{R}^d$ and a bias $b$. For a new image, $p = \sigma(\mathbf{w}^\top \mathbf{x} + b)$ is the predicted probability of class 1, and thresholding it (at 0.5 by default) gives the predicted label. For $K$ classes, the output is a weight matrix and a bias vector, and the prediction a probability for each class.

Intuition

Each feature contributes evidence. A patch whose gradient energy is high is more likely to contain an edge; a patch whose gradients all point the same way is more likely still. Logistic regression adds the contributions up, each scaled by a learned weight, into a single score $a$, and reads the score as log-odds: $a = 0$ means even odds, each unit of $a$ multiplies the odds by $e \approx 2.7$, and the sigmoid converts the odds into a probability. Because the evidence is added, the boundary between the regions labelled 0 and 1 is flat, a line in a two-feature plot and a hyperplane in general.

Fitting asks which weights make the observed labels most probable. Images far on the correct side of the boundary are already explained and contribute little; images near the boundary or on the wrong side pull the weights. Larger weights make the sigmoid steeper and every prediction more confident. That lowers the loss only while the confident predictions are right, so when the classes overlap the fit settles at finite weights, and when they can be separated it does not settle at all.

Algorithm

Fitting by Newton’s method, written as iteratively reweighted least squares (IRLS; Hastie et al., §4.4.1; Murphy, §10.2.6, Algorithm 10.1). Prepend a constant 1 to each feature vector so that the bias becomes the first weight.

1. Start from $\mathbf{w} = \mathbf{0}$, where every probability is $0.5$.
2. Compute the scores $a_i = \mathbf{w}^\top \mathbf{x}_i$ and probabilities $p_i = \sigma(a_i)$.
3. Compute the weights $s_i = p_i(1 - p_i)$ and the working response $z_i = a_i + (y_i - p_i)/s_i$.
4. Solve the weighted least-squares problem $\min_{\mathbf{w}} \sum_i s_i (z_i - \mathbf{w}^\top \mathbf{x}_i)^2$ for the new $\mathbf{w}$.
5. Repeat from step 2 until the gradient or the change in the negative log-likelihood is negligible.

Large data sets are fitted instead by first-order or quasi-Newton methods (gradient descent, stochastic gradient descent, L-BFGS), which need only the gradient below.

Mathematical Formulation

The model and the log-odds

The binary model is (Murphy, 2022, Eq. 10.1–10.2)

$$P(y = 1 \mid \mathbf{x}) = \sigma(\mathbf{w}^\top \mathbf{x} + b), \qquad \sigma(a) = \frac{1}{1 + e^{-a}},$$

and inverting the sigmoid shows what the score means:

$$a = \mathbf{w}^\top \mathbf{x} + b = \log \frac{P(y = 1 \mid \mathbf{x})}{P(y = 0 \mid \mathbf{x})}.$$

The score is the logit, the log-odds of class 1 (Murphy calls it the logit or pre-activation). Increasing feature $j$ by one unit adds $w_j$ to the log-odds and multiplies the odds by $e^{w_j}$, whatever the other features are. Two identities are used throughout: $\sigma(-a) = 1 - \sigma(a)$, and $\sigma'(a) = \sigma(a)\,(1 - \sigma(a))$ (Bishop, Eq. 4.88; Murphy, Eq. 10.15).

The model is linear in the log-odds, not in the probability. One setting in which it is exactly right: if each class’s features are Gaussian with a shared covariance matrix $\Sigma$, the posterior of class 1 is $\sigma(\mathbf{w}^\top \mathbf{x} + w_0)$ with $\mathbf{w} = \Sigma^{-1}(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2)$, because the quadratic terms in $\mathbf{x}$ cancel (Bishop, §4.2.1). Logistic regression estimates the same form directly, without modelling how the features are distributed.

The decision boundary

With equal costs for the two kinds of error, the best decision is class 1 when $P(y = 1 \mid \mathbf{x}) > 0.5$, that is when $a > 0$ (Murphy, §10.2.1). The decision boundary $\mathbf{w}^\top \mathbf{x} + b = 0$ is a hyperplane with normal $\mathbf{w}$; its orientation is set by the direction of $\mathbf{w}$, and the norm $\lVert \mathbf{w} \rVert$ sets how steeply the probability changes across it and hence how confident the predictions are. The signed distance of $\mathbf{x}$ from the boundary is $a / \lVert \mathbf{w} \rVert$. Other thresholds give other operating points on the classifier’s ROC curve; nonlinear boundaries come from nonlinear features, such as polynomial terms (Murphy, §10.2.2) or a network’s learned representation.

Fitting by maximum likelihood

Absorb the bias into $\mathbf{w}$ by giving $X$ a leading column of ones, and write $p_i = \sigma(\mathbf{w}^\top \mathbf{x}_i)$ and $\mathbf{p}$ for the vector of the $p_i$. Each label is a Bernoulli outcome with probability $p_i$, so the negative log-likelihood is the binary cross-entropy (Hastie et al., Eq. 4.20; Murphy, Eq. 10.8, who divides it and its derivatives by $n$):

$$E(\mathbf{w}) = -\sum_{i=1}^{n} \big[ y_i \log p_i + (1 - y_i) \log (1 - p_i) \big] = \sum_{i=1}^{n} \big[ \log\big(1 + e^{\mathbf{w}^\top \mathbf{x}_i}\big) - y_i\, \mathbf{w}^\top \mathbf{x}_i \big].$$

The maximum likelihood and cross-entropy loss articles explain why the negative log-likelihood of a categorical label is the cross-entropy. Its gradient and Hessian are (Hastie et al., Eq. 4.24–4.25; Murphy, Eq. 10.21–10.24)

$$\nabla E = X^\top (\mathbf{p} - \mathbf{y}), \qquad \nabla^2 E = X^\top S X, \qquad S = \operatorname{diag}\big(p_1(1 - p_1), \dots, p_n(1 - p_n)\big).$$

The gradient weights each image by its error $p_i - y_i$; the derivative of the sigmoid has cancelled against the logarithm, as Bishop (§4.3.2) points out, so the gradient has the same form as for least squares. Its component for the bias gives the score equation $\sum_i p_i = \sum_i y_i$ at the optimum: the fitted model predicts as many class-1 images in the training set as there are (Hastie et al., §4.4.1).

Since every $s_i > 0$, $\mathbf{v}^\top X^\top S X \mathbf{v} = \lVert S^{1/2} X \mathbf{v} \rVert^2 \ge 0$, so the Hessian is positive semidefinite and $E$ is convex (Murphy, §10.2.3.4). It is positive definite, and $E$ strictly convex, when $X$ has full column rank; Murphy states the strict inequality without that condition, which fails whenever there are fewer images than weights. Every stationary point is therefore a global minimum, by the argument in optimization. Murphy notes that probabilities close to 0 or 1 make the Hessian nearly singular, which an $L^2$ penalty prevents.

Newton’s method as iteratively reweighted least squares

There is no closed-form solution (Bishop, §4.3.3). Newton’s step $\mathbf{w}^{\text{new}} = \mathbf{w} - (X^\top S X)^{-1} X^\top (\mathbf{p} - \mathbf{y})$ can be rearranged as (Hastie et al., Eq. 4.26–4.28; Murphy, Eq. 10.36–10.43)

$$\mathbf{w}^{\text{new}} = (X^\top S X)^{-1} X^\top S \mathbf{z}, \qquad \mathbf{z} = X\mathbf{w} + S^{-1}(\mathbf{y} - \mathbf{p}),$$

the solution of a weighted least-squares problem with weights $s_i$ and response $\mathbf{z}$. Because $S$ and $\mathbf{z}$ change with $\mathbf{w}$, the problem is solved repeatedly, hence the name. Hastie and colleagues suggest $\mathbf{w} = \mathbf{0}$ as a starting value and note that convergence is never guaranteed but typical, since the log-likelihood is concave; in the rare cases that a step overshoots and the log-likelihood decreases, halving the step guarantees convergence. Close to the optimum the convergence is quadratic, as for Newton’s method in general (see optimization).

Regularization and scikit-learn’s C

Adding a penalty to $E$ gives MAP estimation under a prior on the weights (Murphy, §10.2.7). With an $L^2$ penalty,

$$E_\lambda(\mathbf{w}) = E(\mathbf{w}) + \frac{\lambda}{2} \lVert \mathbf{w} \rVert^2, \qquad \nabla E_\lambda = X^\top(\mathbf{p} - \mathbf{y}) + \lambda \mathbf{w}, \qquad \nabla^2 E_\lambda = X^\top S X + \lambda I,$$

where the bias is usually left out of the penalty (Hastie et al., §4.4.4). The Hessian is now positive definite for any $\lambda > 0$, so the minimizer is unique. The larger $\lambda$, the less flexible the model; Murphy (§10.2.7, Figure 10.6) illustrates this on a logistic regression with degree-4 polynomial features. An $L^1$ penalty $\lambda \lVert \mathbf{w} \rVert_1$ drives some weights to exactly zero and selects features; Hastie and colleagues (§4.4.4) recommend standardizing the features first so that the penalty treats them equally. The regularization article compares the two penalties.

scikit-learn writes the same objective with the inverse strength $C$. Its user guide for version 1.7.2 gives the binary objective as the weighted average of the cross-entropy over the training images, with sample and class weights, plus the penalty $r(\mathbf{w})$ divided by $C$ times the total weight; $r(\mathbf{w}) = \tfrac{1}{2}\lVert \mathbf{w} \rVert^2$ for the default $L^2$ penalty, and the intercept is not in it. Multiplying through shows that, without sample or class weights, this is $C \cdot E(\mathbf{w}) + \tfrac{1}{2}\lVert \mathbf{w} \rVert^2$, so $\lambda = 1/C$: a smaller $C$ means stronger regularization, as the documentation says. Because $C$ multiplies the summed loss, a fixed $C$ penalizes relatively less as the training set grows; Murphy’s $\lambda$ (§10.2.7) instead multiplies $\lVert \mathbf{w} \rVert^2$ added to the mean loss.

Several classes: multinomial regression and one-vs-rest

Multinomial logistic regression (softmax regression) gives each class $k$ a weight vector $\mathbf{w}_k$ and bias $b_k$, and turns the $K$ scores into probabilities with the softmax (Murphy, Eq. 10.54–10.55):

$$P(y = k \mid \mathbf{x}) = \frac{\exp(\mathbf{w}_k^\top \mathbf{x} + b_k)}{\sum_{j=1}^{K} \exp(\mathbf{w}_j^\top \mathbf{x} + b_j)}.$$

The negative log-likelihood is the categorical cross-entropy, again convex, with gradient $X^\top (P - Y)$ for the matrices of predicted probabilities $P$ and one-hot labels $Y$, and a Hessian of size $Kd \times Kd$ (Murphy, §10.3.2). The boundary between classes $k$ and $j$ is the hyperplane where their scores are equal. Because adding the same vector to every $\mathbf{w}_k$ changes nothing, $K - 1$ weight vectors suffice, which is how Hastie and colleagues write the model (Eq. 4.17), with the last class as reference; for $K = 2$ it reduces to the binary model, the two-class identity of the softmax. scikit-learn keeps all $K$ vectors, a choice its user guide makes for ease of implementation and so that no class is singled out, and notes that the solution may then not be unique without a penalty. With an $L^2$ penalty it is unique: at the optimum, the weights of each feature sum to zero over the classes (Murphy, §10.3.5).

One-vs-rest (OvR) instead fits $K$ independent binary models, class $k$ against all others, and predicts the class with the highest score. Bishop (§4.1.2) shows that combining two-class decisions by their signs leaves regions of input space claimed by no class or by several, and that a single $K$-class discriminant that picks the largest of $K$ linear scores avoids them; one-vs-rest with an arg max over scores is of this form, but its scores were fitted separately, so their probabilities do not form a distribution (our reasoning; the fourth example measures the consequence).

Parameters

Each parameter below is named as in scikit-learn 1.7.2’s LogisticRegression, with its default as the 1.7.2 documentation states it.

• Regularization (penalty='l2', C=1.0). The documentation stresses that regularization is applied by default, which its user guide calls common in machine learning but not in statistics; penalty=None gives the plain maximum likelihood fit, and 'l1' and 'elasticnet' are available with some solvers. Choose $C$ on a validation set or by cross-validation (LogisticRegressionCV builds this in), usually over a logarithmic grid. Standardize the features first: the penalty treats all weights alike, which presumes features of similar scale (Murphy, §10.2.8).
• Solver (solver='lbfgs', max_iter=100, tol=1e-4). L-BFGS supports the $L^2$ penalty or none, and the multinomial loss. The alternatives are 'newton-cg', 'newton-cholesky' (which computes the full Hessian, with memory quadratic in the number of features times classes), 'sag' and 'saga' (stochastic average gradient methods for large data; only 'saga' supports elastic net), and 'liblinear', which supports $L^1$ and $L^2$ but handles only binary problems natively.
• Class weights (class_weight=None). A dictionary or 'balanced', which weights each class by $n / (K\, n_k)$ for $n_k$ images in class $k$, inversely to its frequency. The third example shows what this does to the probabilities.
• Intercept (fit_intercept=True). The intercept is outside the penalty in the objective above. With liblinear it is the weight of a synthetic constant feature and is penalized like any other weight; intercept_scaling lessens that.
• Multiclass handling. In 1.7.2, every solver except liblinear minimizes the multinomial loss when there are three or more classes. The multi_class parameter, whose description still lists 'auto', 'ovr' and 'multinomial' with default 'auto' while the signature’s default is 'deprecated', was deprecated in version 1.5 and is due for removal in 1.8; the documentation recommends wrapping the model in OneVsRestClassifier for one-vs-rest. Binary problems are fitted as a single binary model.

Complexity

One gradient costs $O(nd)$ for $n$ images and $d$ features, and one IRLS iteration $O(nd^2 + d^3)$: forming $X^\top S X$ and solving with it. Newton’s method needs few iterations (the first example needs nine), against hundreds for gradient descent on the same problem, so in our reading it is the better choice while $d$ is small enough for a $d \times d$ solve to be cheap. For the multinomial model the Hessian has $(Kd)^2$ entries, which Murphy (§10.3.3) gives as the reason quasi-Newton methods such as L-BFGS are common; for very many images, stochastic gradient methods process them in mini-batches. Prediction costs $O(d)$ per image, or $O(Kd)$ with $K$ classes, the cost of one matrix–vector product. Memory is that of the data plus $O(d^2)$ for a Newton method or $O(md)$ for L-BFGS with $m$ stored updates.

Implementation

The examples classify synthetic $16 \times 16$ grayscale patches with exactly known ground truth: a flat patch of random brightness, a patch crossed by a straight edge of random orientation, offset, polarity and contrast (0.05 to 0.5), or a corner of a brighter or darker region bounded by two edges, all with Gaussian noise of standard deviation 0.15. Two designed features describe each patch: the logarithm of its mean gradient energy, and the coherence of its structure tensor, $\sqrt{(J_{xx} - J_{yy})^2 + 4J_{xy}^2}/(J_{xx} + J_{yy})$, which is near 1 when the gradients share one orientation and near 0 when they have none. Here $J_{xx}$, $J_{yy}$ and $J_{xy}$ are the patch averages of $g_x^2$, $g_y^2$ and $g_x g_y$, with the gradients $g_x, g_y$ taken by derivatives of a Gaussian of standard deviation 1.5 pixels. The outputs came from NumPy 2.3.5, SciPy 1.16.3 and scikit-learn 1.7.2; each example runs in a few seconds.

Newton, gradient descent and scikit-learn

The first block implements IRLS and gradient descent in NumPy for edge versus flat patches, using 400 training patches and 20,000 test patches, and compares both with scikit-learn, first without a penalty and then with the $L^2$ penalty at $C = 1$. A tight tolerance is passed to scikit-learn so that the comparison measures the solution rather than the stopping rule.

import numpy as np
from scipy import ndimage
from sklearn.linear_model import LogisticRegression

S = 16                                            # patch size (pixels)
yy, xx = np.mgrid[0:S, 0:S] - (S - 1) / 2         # coordinates centred on the patch

def render(kind, rng):
    """kind 0: flat patch; 1: straight step edge; 2: corner of a bright or dark region."""
    img = np.full((S, S), rng.uniform(0.2, 0.8))              # random mean brightness
    if kind > 0:
        c = rng.uniform(0.05, 0.5) * rng.choice([-1, 1])      # random contrast and polarity
        t1 = rng.uniform(0, np.pi)
        side = xx * np.cos(t1) + yy * np.sin(t1) > rng.uniform(-3, 3)
        if kind == 2:                                          # second edge at 60-120 degrees
            t2 = t1 + rng.uniform(np.pi / 3, 2 * np.pi / 3)
            side &= xx * np.cos(t2) + yy * np.sin(t2) > rng.uniform(-3, 3)
        img += c * side
    return img + rng.normal(0, 0.15, (S, S))                   # Gaussian noise

def dataset(n, kinds, rng, p=None):
    y = rng.choice(len(kinds), size=n, p=p)                   # label = index into kinds
    P = np.array([render(kinds[k], rng) for k in y])
    return P, y

def features(P):
    """Log mean gradient energy and orientation coherence of the structure tensor."""
    gx = ndimage.gaussian_filter(P, (0, 1.5, 1.5), order=(0, 0, 1))
    gy = ndimage.gaussian_filter(P, (0, 1.5, 1.5), order=(0, 1, 0))
    jxx, jyy, jxy = (gx**2).mean((1, 2)), (gy**2).mean((1, 2)), (gx * gy).mean((1, 2))
    energy = jxx + jyy
    coherence = np.sqrt((jxx - jyy) ** 2 + 4 * jxy**2) / energy   # 0: no dominant orientation, 1: one
    return np.column_stack([np.log(energy), coherence])

def sigmoid(a):
    return 0.5 * (1 + np.tanh(a / 2))                  # equals 1/(1+exp(-a)), without overflow

def nll(w, A, y, lam=0.0):
    a = A @ w                                          # logits
    loss = np.sum(np.logaddexp(0, a) - y * a)          # -sum[y log p + (1-y) log(1-p)]
    return loss + 0.5 * lam * w[1:] @ w[1:]            # intercept w[0] is not penalized

def newton_irls(A, y, lam=0.0, iters=20):
    """Newton's method = iteratively reweighted least squares, from w = 0."""
    w = np.zeros(A.shape[1]); R = lam * np.eye(A.shape[1]); R[0, 0] = 0
    history = [nll(w, A, y, lam)]
    for _ in range(iters):
        p = sigmoid(A @ w)
        s = p * (1 - p)                                # IRLS weights
        g = A.T @ (p - y) + R @ w                      # gradient
        H = A.T @ (s[:, None] * A) + R                 # Hessian X^T S X (+ penalty)
        w = w - np.linalg.solve(H, g)
        history.append(nll(w, A, y, lam))
        if np.abs(g).max() < 1e-10:
            break
    return w, history

rng = np.random.default_rng(0)
Ptr, ytr = dataset(400, [0, 1], rng)                 # 400 training patches: flat (0) or edge (1)
Pte, yte = dataset(20000, [0, 1], rng)               # 20,000 test patches
Ftr, Fte = features(Ptr), features(Pte)
mu, sd = Ftr.mean(0), Ftr.std(0)                     # standardize with training statistics
Atr = np.column_stack([np.ones(len(ytr)), (Ftr - mu) / sd])   # first column: intercept
Ate = np.column_stack([np.ones(len(yte)), (Fte - mu) / sd])

w_nt, hist = newton_irls(Atr, ytr)
print("Newton/IRLS, mean NLL per iteration:", np.round(np.array(hist) / len(ytr), 6))

# Gradient descent on the mean NLL with step 1/L, L = largest eigenvalue of X^T X / (4n),
# an upper bound on the Hessian because p(1-p) <= 1/4.
L = np.linalg.eigvalsh(Atr.T @ Atr / (4 * len(ytr))).max()
w_gd = np.zeros(3); target = hist[-1] + 1e-6 * len(ytr)
for it in range(1, 100001):
    w_gd -= (Atr.T @ (sigmoid(Atr @ w_gd) - ytr) / len(ytr)) / L
    if nll(w_gd, Atr, ytr) < target:
        break
print(f"gradient descent: {it} iterations to come within 1e-6 of the minimum mean NLL")

sk = LogisticRegression(penalty=None, tol=1e-10, max_iter=1000).fit(Atr[:, 1:], ytr)
w_sk = np.r_[sk.intercept_, sk.coef_.ravel()]
print("weights [b, w_energy, w_coherence]:")
print("  Newton/IRLS     ", np.round(w_nt, 4))
print("  gradient descent", np.round(w_gd, 4))
print("  scikit-learn    ", np.round(w_sk, 4), "(penalty=None)")

# scikit-learn's default penalty: L2 with C = 1.0, i.e. minimize C * sum(NLL) + ||w||^2 / 2
sk1 = LogisticRegression(tol=1e-10).fit(Atr[:, 1:], ytr)
w_nt1, _ = newton_irls(Atr, ytr, lam=1 / 1.0)
print("  Newton, lam=1/C ", np.round(w_nt1, 4))
print("  scikit-learn    ", np.round(np.r_[sk1.intercept_, sk1.coef_.ravel()], 4), "(L2, C=1.0)")
err = np.mean((Ate @ w_nt > 0) != yte)
print(f"test error at p = 0.5: {err:.3f}")

Output:

Newton/IRLS, mean NLL per iteration: [0.693147 0.346037 0.270007 0.237095 0.226016 0.224217 0.224158 0.224158
 0.224158 0.224158]
gradient descent: 443 iterations to come within 1e-6 of the minimum mean NLL
weights [b, w_energy, w_coherence]:
  Newton/IRLS      [2.4018 4.3606 1.2164]
  gradient descent [2.3919 4.3432 1.2168]
  scikit-learn     [2.4018 4.3606 1.2164] (penalty=None)
  Newton, lam=1/C  [1.8602 3.3473 1.2248]
  scikit-learn     [1.8602 3.3473 1.2248] (L2, C=1.0)
test error at p = 0.5: 0.087

Newton’s method reaches the minimum in nine steps, and the mean negative log-likelihood stops changing in its sixth decimal after the sixth: the last steps gain several correct digits each, the fast local convergence of Newton’s method. Gradient descent with a safe fixed step needs several hundred iterations to get within $10^{-6}$ of the same value, and stops with weights that still differ in the second decimal: near the minimum the objective is flat along some directions, so a small loss gap allows a visible weight error. Newton’s method and scikit-learn agree to all four printed decimals, and Newton with $\lambda = 1/C$ reproduces scikit-learn’s penalized fit, which confirms the reading of $C$ above.

The weights are read as log-odds per standard deviation of each feature. One standard deviation more log gradient energy adds 4.36 to the log-odds, multiplying the odds of an edge by about $e^{4.36} \approx 78$; one standard deviation more coherence multiplies them by about $e^{1.22} \approx 3.4$. The penalty at $C = 1$ shrinks the energy weight by about a quarter and leaves the coherence weight almost unchanged. The test error of 0.087 is not a fitting failure: the faintest edges, with a contrast as low as a third of the noise level, have features that overlap those of flat patches, so no boundary in this feature plane can classify every patch correctly.

Separable data: the weights diverge

The second block, which continues from the first, describes each patch by its 256 raw pixel values and trains on only 40 patches. With more weights than images the two classes can be separated by a hyperplane whatever the labels are, and the maximum likelihood estimate does not exist: as Bishop (§4.3.2) explains, the likelihood keeps increasing as a separating hyperplane’s $\lVert \mathbf{w} \rVert$ grows, the sigmoid tends to a step and every training image gets probability 1 for its own class; Hastie and colleagues (§4.4.5) call the estimates undefined, that is, infinite. The block runs unregularized gradient descent and then scikit-learn along a path of increasing $C$, and compares the direction of $\mathbf{w}$ with that of a linear support vector machine with a very large $C$, which approximates the separating hyperplane with the largest margin.

from sklearn.svm import SVC

# 40 training patches described by their 256 raw pixel values: more weights than images,
# so the two classes are linearly separable whatever the labels.
rng = np.random.default_rng(1)
Ptr, ytr = dataset(40, [0, 1], rng)
Pte, yte = dataset(5000, [0, 1], rng)
Xtr = Ptr.reshape(len(ytr), -1); Xte = Pte.reshape(len(yte), -1)
mu = Xtr.mean(0); Xtr = Xtr - mu; Xte = Xte - mu
Atr = np.column_stack([np.ones(len(ytr)), Xtr])

# Unregularized gradient descent: the loss keeps falling and the weights keep growing.
L = np.linalg.eigvalsh(Atr.T @ Atr / (4 * len(ytr))).max()
w = np.zeros(Atr.shape[1]); checkpoints = [10, 100, 1000, 10000, 100000]
print("iteration   mean NLL   ||w||    smallest margin y*a")
for it in range(1, checkpoints[-1] + 1):
    w -= (Atr.T @ (sigmoid(Atr @ w) - ytr) / len(ytr)) / L
    if it in checkpoints:
        a = Atr @ w
        print(f"{it:9d}   {nll(w, Atr, ytr) / len(ytr):.2e}   {np.linalg.norm(w[1:]):6.2f}"
              f"   {np.min((2 * ytr - 1) * a):6.2f}")

# The regularization path: as C grows the norm diverges, but the direction converges
# to that of the maximum-margin separating hyperplane (hard-margin linear SVM).
svm = SVC(kernel="linear", C=1e10).fit(Xtr, ytr)
u = svm.coef_.ravel() / np.linalg.norm(svm.coef_)
print("\n      C      ||w||   cos(w, SVM)   train NLL   test error")
for C in [1e-2, 1e0, 1e2, 1e4, 1e6, 1e8]:
    m = LogisticRegression(C=C, tol=1e-12, max_iter=100000).fit(Xtr, ytr)
    v = m.coef_.ravel()
    p = m.predict_proba(Xtr)[:, 1]
    tr_nll = -np.mean(ytr * np.log(p) + (1 - ytr) * np.log1p(-p))
    print(f"{C:8.0e}  {np.linalg.norm(v):8.2f}   {v @ u / np.linalg.norm(v):9.4f}"
          f"   {tr_nll:9.2e}   {np.mean(m.predict(Xte) != yte):.3f}")

Output:

iteration   mean NLL   ||w||    smallest margin y*a
       10   5.56e-01     0.70    -0.03
      100   1.77e-01     4.09     1.14
     1000   1.86e-02    10.06     3.48
    10000   1.80e-03    16.07     5.83
   100000   1.79e-04    21.99     8.15

      C      ||w||   cos(w, SVM)   train NLL   test error
   1e-02      0.08      0.9185    6.71e-01   0.491
   1e+00      3.29      0.9855    2.35e-01   0.459
   1e+02     12.27      0.9976    7.89e-03   0.464
   1e+04     22.53      0.9990    1.44e-04   0.466
   1e+06     33.34      0.9995    2.14e-06   0.467
   1e+08     44.41      0.9997    2.85e-08   0.468

Gradient descent never stops making progress. Each tenfold increase in iterations divides the training loss by about ten, adds about six to $\lVert \mathbf{w} \rVert$ and widens the smallest margin $y \cdot a$ (with $y$ coded as $\pm 1$): the weights grow without bound, only slowly, and every training patch is pushed toward certainty. Newton’s method fails sooner: with 40 images and 257 weights, $X^\top S X$ has rank at most 40 and cannot be inverted at any step (our reasoning). Along the regularization path the norm likewise grows without limit as $C$ increases, while the direction settles: its cosine with the maximum-margin direction rises from 0.92 to 0.9997. This is the behaviour Hastie and colleagues note in a footnote to §12.3.2: for separable data the penalized estimate diverges as the penalty vanishes, but its direction converges to the optimal separating hyperplane’s.

None of these fits generalizes; the test error stays between 0.46 and 0.49, close to guessing. In our reading the cause is the representation rather than the fitting: an edge of either polarity at any orientation and position is not a linear function of pixel values, since reversing an edge’s polarity reverses the sign of the edge’s contribution to every pixel and with it to the score. Regularization keeps the weights finite and the problem well posed; it cannot supply features the model lacks. The same divergence happens in feature spaces of modest dimension whenever the training classes happen to be separable, which Bishop notes can occur even with many more images than parameters.

Calibration, class weights and overfitting

The third block, which continues from the first, makes edges rare, 10% of all patches, as edge-containing windows might be in a scan of mostly uniform regions. It compares three models: the two-feature model with scikit-learn’s defaults; the same model with class_weight='balanced'; and a model on raw pixels from only 300 patches with almost no penalty ($C = 10^4$). For each it prints a reliability table over 20,000 test patches, comparing the mean predicted probability of an edge in each of ten equal-width bins with the fraction of edges actually in that bin, and the binary expected calibration error computed from those bins. It repeats the experiment for five seeds and prints the range of the ECE.

def reliability(p, y, bins=10):
    """Per equal-width bin: count, mean predicted probability, observed frequency; and the ECE."""
    idx = np.minimum((p * bins).astype(int), bins - 1)
    rows, ece = [], 0.0
    for b in range(bins):
        sel = idx == b
        if sel.any():
            rows.append((b / bins, sel.sum(), p[sel].mean(), y[sel].mean()))
            ece += sel.mean() * abs(p[sel].mean() - y[sel].mean())
    return rows, ece

def experiment(seed, show=False):
    rng = np.random.default_rng(seed)
    prior = [0.9, 0.1]                                  # edges are rare: 10% of patches
    Ptr, ytr = dataset(2000, [0, 1], rng, prior)
    Pte, yte = dataset(20000, [0, 1], rng, prior)
    Ftr, Fte = features(Ptr), features(Pte)
    mu, sd = Ftr.mean(0), Ftr.std(0)
    Ztr, Zte = (Ftr - mu) / sd, (Fte - mu) / sd
    models = {
        "features": (LogisticRegression().fit(Ztr, ytr), Zte),
        "features, balanced": (LogisticRegression(class_weight="balanced").fit(Ztr, ytr), Zte),
    }
    if show:
        plain, bal = models["features"][0], models["features, balanced"][0]
        print(f"training set: edge fraction {ytr.mean():.4f}, "
              f"mean predicted p {plain.predict_proba(Ztr)[:, 1].mean():.4f}")
        print(f"intercept shift from balanced weights: {bal.intercept_[0] - plain.intercept_[0]:.3f}"
              f"  (log of the class ratio: {np.log((ytr == 0).sum() / (ytr == 1).sum()):.3f})")
        print(f"weights [w_energy, w_coherence]: plain {np.round(plain.coef_.ravel(), 3)},"
              f" balanced {np.round(bal.coef_.ravel(), 3)}")
    # Raw pixels, few images, almost no penalty: the model can nearly separate its training set.
    Xtr = Ptr[:300].reshape(300, -1); Xte = Pte.reshape(len(yte), -1)
    m = Xtr.mean(0)
    models["pixels, C=1e4"] = (LogisticRegression(C=1e4, max_iter=5000).fit(Xtr - m, ytr[:300]), Xte - m)
    out = {}
    for name, (model, Z) in models.items():
        p = model.predict_proba(Z)[:, 1]
        rows, ece = reliability(p, yte)
        pred = p > 0.5
        conf = np.maximum(p, 1 - p).mean()             # mean confidence in the predicted label
        out[name] = (ece, conf, np.mean(pred == yte), pred[yte == 1].mean(), yte[pred].mean())
        if show and name != "pixels, C=1e4":
            print(f"{name}: reliability on 20,000 test patches (10% edges)")
            print("   bin    count   mean p   edge freq")
            for lo, n, mp, fr in rows:
                print(f"   {lo:.1f}  {n:7d}   {mp:.3f}    {fr:.3f}")
    return out

res = experiment(0, show=True)
all_runs = [experiment(s) for s in range(5)]
print("\nmodel                ECE    ECE, seeds 0-4  confidence  accuracy  recall  precision")
for name in res:
    eces = [r[name][0] for r in all_runs]
    ece, conf, acc, rec, prec = res[name]
    print(f"{name:20s} {ece:.3f}  {min(eces):.3f}-{max(eces):.3f}     {conf:.3f}     {acc:.3f}   {rec:.3f}   {prec:.3f}")

Output:

training set: edge fraction 0.0995, mean predicted p 0.0995
intercept shift from balanced weights: 2.191  (log of the class ratio: 2.203)
weights [w_energy, w_coherence]: plain [3.048 1.049], balanced [2.859 0.826]
features: reliability on 20,000 test patches (10% edges)
   bin    count   mean p   edge freq
   0.0    17458   0.015    0.014
   0.1      630   0.138    0.108
   0.2      185   0.244    0.222
   0.3       82   0.344    0.512
   0.4       64   0.451    0.609
   0.5       36   0.545    0.778
   0.6       47   0.651    0.872
   0.7       45   0.758    0.978
   0.8       84   0.852    1.000
   0.9     1369   0.993    1.000
features, balanced: reliability on 20,000 test patches (10% edges)
   bin    count   mean p   edge freq
   0.0    10993   0.040    0.007
   0.1     3420   0.143    0.020
   0.2     1673   0.245    0.028
   0.3      969   0.346    0.036
   0.4      568   0.446    0.065
   0.5      352   0.547    0.088
   0.6      236   0.647    0.195
   0.7      131   0.746    0.389
   0.8      116   0.846    0.647
   0.9     1542   0.992    0.995

model                ECE    ECE, seeds 0-4  confidence  accuracy  recall  precision
features             0.006  0.005-0.008     0.974     0.977   0.782   0.991
features, balanced   0.100  0.100-0.118     0.889     0.955   0.868   0.731
pixels, C=1e4        0.143  0.088-0.171     0.982     0.848   0.144   0.180

The plain two-feature model is well calibrated where most patches fall: in the lowest and highest bins, which hold over 90% of the test set, the predicted probabilities match the observed edge frequencies to within about 0.01, and the ECE is under 0.01 for every seed. On the training set its mean predicted probability equals the fraction of edges exactly, the score equation for the intercept; the penalty does not disturb it because the intercept is not penalized. Calibration is not perfect everywhere, though: in the sparsely populated bins from 0.3 to 0.8, edges are more frequent than predicted, by up to 0.23. In our reading the cause is the model’s form. The true log-odds of an edge given these two features is not exactly linear in them, and the maximum likelihood fit matches the frequencies on average, as the score equations require, but not in every region of feature space.

Balanced class weights make the training set behave as if edges and flat patches were equally common. The intercept rises by 2.19, close to the logarithm of the training ratio of flat patches to edges, 2.20. By our reasoning this is what a correctly specified model would give in the limit of many images: weighting one class by a factor $r$ relative to the other multiplies the odds that minimize the weighted loss by exactly $r$ at every $\mathbf{x}$, which adds $\log r$ to the intercept and leaves the feature weights unchanged. With a finite sample and a model that is only approximately right, the shift is an approximation: here the feature weights also change, the coherence weight by about a fifth. The probabilities therefore overestimate edges in every bin below 0.9, the ECE rises more than tenfold, and the 0.5 threshold moves to where the unweighted model gives an edge a probability of about 0.1, raising recall from 0.78 to 0.87 and lowering precision from 0.99 to 0.73. By the same reasoning, nearly the same operating point could be had from the unweighted model by lowering its threshold to about 0.1, keeping its probabilities calibrated; the ROC curve and precision and recall articles describe how to choose it.

The pixel model overfits its 300 training patches. Its mean confidence, 0.98, is far above its accuracy, 0.85, lower than the roughly 0.90 of always answering “flat”; across seeds its ECE ranges from 0.09 to 0.17.

Multinomial versus one-vs-rest

The fourth block, which continues from the first, has three classes, flat, edge and corner, described by the same two features. It fits scikit-learn’s multinomial model, checks it against a softmax regression minimized with SciPy’s L-BFGS, and compares it with one-vs-rest on 20,000 test patches over five seeds.

from scipy.optimize import minimize
from scipy.special import logsumexp
from sklearn.multiclass import OneVsRestClassifier
from sklearn.metrics import log_loss

def softmax_objective(theta, A, Y, C):
    """C * sum of cross-entropies + ||W||^2 / 2 (intercepts unpenalized), and its gradient."""
    W = theta.reshape(A.shape[1], Y.shape[1])
    Z = A @ W                                          # logits, one column per class
    logP = Z - logsumexp(Z, axis=1, keepdims=True)
    f = -C * np.sum(Y * logP) + 0.5 * np.sum(W[1:] ** 2)
    G = C * A.T @ (np.exp(logP) - Y)                   # A^T (P - Y)
    G[1:] += W[1:]
    return f, G.ravel()

def run(seed):
    rng = np.random.default_rng(seed)
    Ptr, ytr = dataset(600, [0, 1, 2], rng)           # flat, edge, corner
    Pte, yte = dataset(20000, [0, 1, 2], rng)
    Ftr, Fte = features(Ptr), features(Pte)
    mu, sd = Ftr.mean(0), Ftr.std(0)
    Ztr, Zte = (Ftr - mu) / sd, (Fte - mu) / sd
    multi = LogisticRegression(tol=1e-10).fit(Ztr, ytr)                     # multinomial (softmax)
    ovr = OneVsRestClassifier(LogisticRegression(tol=1e-10)).fit(Ztr, ytr)  # three binary models
    A = np.column_stack([np.ones(len(ytr)), Ztr]); Y = np.eye(3)[ytr]
    sol = minimize(softmax_objective, np.zeros(9), args=(A, Y, 1.0), jac=True,
                   method="L-BFGS-B", options={"gtol": 1e-10, "maxiter": 10000})
    W = sol.x.reshape(3, 3)
    W_sk = np.vstack([multi.intercept_, multi.coef_.T])
    raw = np.column_stack([e.predict_proba(Zte)[:, 1] for e in ovr.estimators_])
    pm, po = multi.predict_proba(Zte), ovr.predict_proba(Zte)
    return dict(diff=np.abs(W - W_sk).max(), colsum=np.abs(multi.coef_.sum(0)).max(),
                err_m=np.mean(pm.argmax(1) != yte), err_o=np.mean(po.argmax(1) != yte),
                nll_m=log_loss(yte, pm), nll_o=log_loss(yte, po),
                disagree=np.mean(pm.argmax(1) != po.argmax(1)),
                rawsum=(raw.sum(1).min(), raw.sum(1).max()))

r = run(0)
print(f"own softmax regression vs scikit-learn, largest weight difference: {r['diff']:.1e}")
print(f"sum of each feature's weights over the classes (L2 makes it zero): {r['colsum']:.1e}")
print(f"one-vs-rest: the three binary probabilities sum to {r['rawsum'][0]:.2f}-{r['rawsum'][1]:.2f}"
      " before normalization")
print("\nseed   error multinomial  error OvR   log loss multinomial  log loss OvR   disagree")
for s in range(5):
    r = run(s)
    print(f"{s:4d}   {r['err_m']:.4f}             {r['err_o']:.4f}      {r['nll_m']:.4f}"
          f"                {r['nll_o']:.4f}         {r['disagree']:.4f}")

Output:

own softmax regression vs scikit-learn, largest weight difference: 1.8e-05
sum of each feature's weights over the classes (L2 makes it zero): 4.0e-15
one-vs-rest: the three binary probabilities sum to 0.56-1.75 before normalization

seed   error multinomial  error OvR   log loss multinomial  log loss OvR   disagree
   0   0.1641             0.1648      0.4457                0.4904         0.0204
   1   0.1674             0.1693      0.4450                0.4915         0.0217
   2   0.1639             0.1652      0.4539                0.5018         0.0192
   3   0.1623             0.1671      0.4425                0.4873         0.0232
   4   0.1648             0.1676      0.4558                0.5029         0.0219

The hand-written softmax regression and scikit-learn’s default multiclass model are the same estimator, agreeing to about $10^{-5}$, and the $L^2$ penalty has made each feature’s weights sum to zero over the three classes, as Murphy’s identifiability argument predicts. The three one-vs-rest models were never fitted to form a distribution: their probabilities for one patch sum to anything from about 0.56 to 1.75, and scikit-learn’s OneVsRestClassifier divides them by their sum. The two schemes label about 2% of test patches differently and have nearly the same error, the multinomial model’s being lower for all five seeds by up to 0.005. The difference is clearer in the probabilities: the multinomial model’s test log loss is 0.045 to 0.048 lower for every seed, because only it was fitted as one distribution over the three classes.

Properties and Behavior

• Convex fit, unique answer. The negative log-likelihood is convex, so any solver that converges finds the global minimum, and with an $L^2$ penalty or a full-rank, non-separable training set that minimum is unique. Different solvers differ in speed, not in the answer, as the first example shows.
• Separation breaks maximum likelihood. On linearly separable training data the maximum likelihood estimate does not exist, as the second example shows. Bishop (§4.3.2) adds that a continuum of separating hyperplanes then gives the same training probabilities, so the solution found depends on the optimizer and the initialization; any positive penalty restores a finite, unique solution.
• Calibration under conditions. The score equations make the fitted probabilities match the observed class frequencies on average over the training set, and, without a penalty, also when weighted by each feature, $\sum_i x_{ij} p_i = \sum_i x_{ij} y_i$ (Hastie et al., Eq. 4.21); with an $L^2$ penalty only the intercept’s equation still holds exactly. Calibration at every probability level additionally needs a correct linear log-odds model, enough data for the number of weights, unweighted classes and a test distribution like the training one; the third example breaks the first three. On eight binary classification problems, Niculescu-Mizil and Caruana (2005) found logistic regression, like neural networks and bagged trees, well calibrated before any correction, with post-hoc calibration of little help on most problems, whereas the maximum-margin methods (boosted trees, boosted stumps and SVMs) pushed their predictions away from 0 and 1.
• Class imbalance shifts the intercept. Calibrated probabilities of a rare class are small, so the 0.5 threshold misses many of its members. Reweighting the classes, or resampling the data to balance them (Murphy, §10.3.8.2), mainly moves the boundary’s offset, a threshold change paid for with calibration, as the third example shows; choosing the threshold on the calibrated model keeps the probabilities interpretable.
• Inference on the weights. If the model is correct, the maximum likelihood weights are consistent and approximately Gaussian with covariance $(X^\top S X)^{-1}$, which gives standard errors and tests for whether a feature matters (Hastie et al., §4.4.3). Statistics uses logistic regression mostly this way, as a tool for understanding the role of each input (Hastie et al., §4.4.1); vision uses it mostly to predict.
• A one-layer neural network. A network with no hidden layer, a sigmoid output and the binary cross-entropy loss is binary logistic regression, and a softmax output with the categorical cross-entropy gives multinomial logistic regression, the “softmax classifier with no hidden layer” of the cross-entropy loss example. In a deep classifier the last layer is a logistic regression on the features the earlier layers compute; Niculescu-Mizil and Caruana describe a network with a single sigmoid output as a linear classifier in the span of its hidden units, with a sigmoid that calibrates its predictions. The neural network and multilayer perceptron articles add the hidden layers. Murphy (§10.2.5) shows the perceptron as a limiting case that replaces the sigmoid by a step function.
• Linear probing. A logistic regression on the frozen features of a pretrained network is a common choice of linear probe in transfer learning, and a linear classifier of this kind, trained on the full labelled training set, is a standard protocol for evaluating representations learned by self-supervised learning. By our reasoning, convexity is what makes the probe a fair measure of the features: its result depends on the features, the labels and the penalty, not on an initialization. With embeddings of hundreds or thousands of dimensions and fewer labelled images than dimensions, the training classes are separable whatever the labels, as in the second example, so the penalty matters.
• Discriminative, not generative. Logistic regression models $P(y \mid \mathbf{x})$ and ignores how the features are distributed. Linear discriminant analysis gives posteriors of exactly the same linear-logit form by fitting a Gaussian to each class with a covariance matrix shared by all classes. Hastie and colleagues (§4.4.5) compare the two: if the classes really are Gaussian with shared covariance, ignoring that loses, in the worst case, about 30% efficiency asymptotically in the error rate (Efron, 1975, as they cite him); but LDA is not robust to gross outliers, logistic regression relies on fewer assumptions and is generally felt to be the safer bet, and in their experience the two give very similar results.

Limitations

• Linear boundaries only. The decision boundary is a hyperplane in the feature space, so whatever is not linearly separable there must be made so by the features. Raw pixels are a poor space for it, as the second example shows; polynomial features (Murphy, §10.3.1), kernels, or features learned by a convolutional network are the usual remedies.
• Needs a penalty, and the penalty needs tuning. Without one the fit diverges on separable data; with one, the result depends on $C$ and on the scale of the features.
• Probabilities are only as good as the model. Outside the conditions above, and especially under distribution shift, the predicted probabilities can be badly miscalibrated.
• Mutually exclusive labels for the softmax. Images that can carry several labels at once need one sigmoid per label, a multi-label model (Murphy, §10.3), as the cross-entropy loss article describes.
• Mislabelled images. Every training image contributes to the likelihood, and a mislabelled image deep in the other class’s region pulls on the weights; Murphy (§10.4) describes robust variants, such as a mixture likelihood that models label noise explicitly.

Variants

• $L^1$ and elastic-net logistic regression, which select features by driving weights to zero (Hastie et al., §4.4.4, for the $L^1$ penalty); scikit-learn provides $L^1$ through the liblinear and saga solvers and elastic net through saga.
• Platt scaling fits a one-dimensional logistic regression to the scores of another classifier to turn them into probabilities; Platt proposed it for support vector machines in a book chapter published in 2000. As Niculescu-Mizil and Caruana (2005) describe it, the sigmoid’s two parameters are fitted by maximum likelihood on an independent calibration set, with slightly softened targets to limit overfitting. Their experiments found it most effective for methods whose miscalibration has a sigmoid shape, such as SVMs and boosted trees, and of little benefit for already calibrated models such as logistic regression itself. Temperature scaling, described in the cross-entropy loss article, does a similar job for a network’s softmax with a single parameter that rescales the logits.
• Probit regression replaces the sigmoid by the standard normal distribution function; the two curves are very similar over a wide central range (Cramer, 2002), and the probit is the older of the two in bio-assay.
• Logistic regression trained by stochastic gradient descent, as in scikit-learn’s SGDClassifier with the log loss, for data sets too large for full-batch solvers; Murphy (§10.2.4) notes that with a suitably decaying step size it converges to the global optimum, since the objective is convex.
• Bayesian logistic regression keeps a posterior over the weights instead of a point estimate, to capture the uncertainty that matters when data are few; the posterior has no closed form and is approximated, for example by a Gaussian at the MAP estimate (the Laplace approximation), and predictions average over it (Murphy, §10.5).
• Support vector machines replace the logistic loss by the hinge loss. Hastie and colleagues (§12.3.2, Table 12.1) compare the two: both give no or little penalty to points well inside their margin and a linear one to points far on the wrong side, but the hinge loss estimates only the class while the logistic loss estimates the log-odds.