Logistic regression

Problem

Many targets are yes-or-no: will this student pass, is this email spam, will this customer churn. Linear regression can't model them well. Its predictions run off past 0 and 1, squared error punishes a confident correct answer for being "too far" from the label, and one extreme point can tip the whole line.

What you want is a probability of the positive class, between 0 and 1, that still comes from an interpretable weighted sum of features. Logistic regression is exactly that, and it is far more than a baseline: it is the output layer of nearly every neural classifier, the model behind most click-through and credit-risk scores, and the template for the cross-entropy loss that trains language models.

Intuition

Keep the linear score z = w·x + b, but stop reading it as the answer. Read it as evidence: each feature pushes the score up or down, and z = 0 means "no evidence either way". Then squash the score into a probability with the sigmoid

σ(z) = 1 / (1 + e^−z)

which maps 0 to ½, large positive scores toward 1 and large negative ones toward 0. Run the sigmoid backwards and you find what z really is: the log-odds, log(p / (1 − p)). The model is linear in log-odds, so each unit of a feature multiplies the odds of the positive class by e^w, wherever you start.

Twelve students' hours studied and whether they passed: the fitted model is z = 1.424·hours − 4.630. Each extra hour multiplies the odds of passing by e^1.424 ≈ 4.16, and the probability crosses ½ where z = 0, at 3.25 hours.

Left: probabilities. Right: the same model as log-odds, where it is a plain line. The sigmoid only bends that line into the interval (0, 1).

Mechanics

The model

For an example x with d features, logistic regression predicts

p(y = 1 | x) = σ(w·x + b)    p(y = 0 | x) = 1 − σ(w·x + b) = σ(−(w·x + b))

Three sigmoid facts carry the rest of the page: σ(−z) = 1 − σ(z); its inverse is the logit, log(p / (1 − p)); and its derivative is σ′(z) = σ(z)(1 − σ(z)), largest (¼) at z = 0 and vanishing for confident predictions.

The cross-entropy loss and its gradient

Choose w and b to maximize the probability the model assigns to the labels it saw. With p_i = σ(z_i), the likelihood of one example is p_i^y_i(1 − p_i)^1−y_i; take minus the log and average to get the log loss, the binary cross-entropy between the labels and the predictions:

L(w, b) = −(1/n) Σ_i [ y_i log p_i + (1 − y_i) log(1 − p_i) ]

A confident wrong answer costs a lot (−log 0.01 ≈ 4.6) and a confident right one almost nothing; the loss-functions page computes it straight from logits, as log(1 + e^−z) for a positive label, so it never takes the log of a rounded-off zero.

Differentiate with the chain rule. For one example, ∂ℓ/∂p = −y/p + (1 − y)/(1 − p) and ∂p/∂z = p(1 − p); their product collapses to

∂ℓ/∂z = p − y    ∇_wL = (1/n) Xᵀ(p − y)    ∂L/∂b = (1/n) Σ (p_i − y_i)

The error signal is just "prediction minus label", the same form as linear regression's residual. That is not a coincidence: both are generalized linear models with their natural link, and the same p − y appears at the output of every softmax classifier.

Setting the bias gradient to zero gives a check you can run on any unpenalized fit: the predicted probabilities sum to the number of positives. And unlike squared error, there is no closed form; ∇L = 0 is nonlinear in w. But L is convex (its Hessian, below, is positive semidefinite), so any method that reaches a stationary point has found the global minimum.

Newton's method, a.k.a. IRLS

Because the loss is smooth and convex with a cheap Hessian, Newton's method is the classic solver. Differentiate the gradient once more:

H = (1/n) Xᵀ S X,   S = diag(p_i(1 − p_i))    w ← w − H⁻¹∇L

Each step solves a weighted least-squares problem, with weight p_i(1 − p_i) on example i, which is why statisticians call it iteratively reweighted least squares. Examples the model is already sure about get almost no weight; the uncertain ones near the boundary drive the update. Near the optimum Newton converges quadratically: the number of correct digits roughly doubles per step.

Gradient descent uses only ∇L. It is cheaper per step (no d × d solve) but its speed depends on the conditioning, just as on the linear-regression page, and uncentered features make it slow.

Four Newton steps beat two thousand gradient steps. Each Newton step costs a d × d solve, which is trivial here and prohibitive with a million features; that tradeoff is why scikit-learn defaults to L-BFGS, which approximates the Hessian from recent gradients.

The decision boundary

Predict the positive class when p ≥ ½, that is when w·x + b ≥ 0. The set w·x + b = 0 is a line in two dimensions and a hyperplane in d, perpendicular to w. So logistic regression is a linear classifier: its contours of equal probability are parallel hyperplanes, and ‖w‖ sets how quickly probability changes as you cross them. The p = 0.9 contour sits where w·x + b = log 9 ≈ 2.2, a distance 2.2 / ‖w‖ from the boundary.

The threshold ½ is a choice, not part of the model. If a missed positive costs ten times a false alarm, predict positive above p = 1/11; the probabilities are what you threshold, which only makes sense if they are calibrated. Unpenalized logistic regression tends to be well calibrated on the data it was fitted to; strong penalties pull probabilities toward ½.

Regularization and C

If the classes are linearly separable, the unpenalized loss has no minimum. Any separating w can be doubled to push every p_i closer to its label and the loss closer to zero, so ‖w‖ grows forever and the probabilities become 0s and 1s. With many features, separability is the norm, not the exception. A penalty fixes it.

scikit-learn writes the penalty through C, an inverse strength, in the SVM style:

min_w,b   C Σ_i ℓ(y_i, w·x_i + b) + ½‖w‖²

Small C means strong regularization, small weights and soft probabilities; large C trusts the data. The default is C = 1, so LogisticRegression() is penalized unless you say otherwise (pass C=np.inf for the plain maximum-likelihood fit). The intercept is not penalized. Because the data term is a sum, the same C is effectively weaker penalty on a bigger dataset. An L1 penalty (l1_ratio=1, with the liblinear or saga solver) zeroes coefficients exactly, as lasso does. As with ridge, scale the features first, or the penalty falls unevenly on them.

Left: only C = 100 separates every training point. The stronger penalties give up the class 1 point at (3, 3.5), predicting 0.48 and 0.42 for it, in exchange for a much smaller ‖w‖, so their probability ramps gently (the dashed contours for C = 1 sit far apart). Right: without a penalty the separable data drives the weights off to infinity.

Many classes: softmax regression

With K classes, give each class its own weight vector and score, z_k = w_k·x + b_k, and turn the scores into a distribution with the softmax:

p(y = k | x) = e^z_k / Σ_j e^z_j    ℓ = −log p(y | x)    ∂ℓ/∂z_k = p_k − [k = y]

The loss is categorical cross-entropy and the gradient is again "prediction minus one-hot label". With K = 2 it reduces to the sigmoid on z₁ − z₀: adding the same number to every score changes nothing, so only differences matter. This is also how a neural classifier ends: an nn.Linear layer producing K logits, then CrossEntropyLoss. In matrix form the scores are XWᵀ + b with W of shape (K, d), which is exactly coef_'s shape in scikit-learn.

The alternative is one-vs-rest: K separate binary models, each "class k versus everything else", with probabilities renormalized afterwards. Since version 1.5 scikit-learn's LogisticRegression always fits the true multinomial model when K ≥ 3 (except with liblinear, which is binary only); wrap it in OneVsRestClassifier to get one-vs-rest.

Solvers

All of scikit-learn's solvers minimize the same objective; they differ in how they get there and which penalties they support.

• lbfgs (default): quasi-Newton, builds a Hessian approximation from the last few gradients. A good general choice for dense data and L2 or no penalty.
• newton-cholesky: exact Newton steps, factoring the d × d Hessian. Best when n is much larger than d and d is modest (say, under a thousand), including one-hot features.
• newton-cg: Newton steps solved approximately by conjugate gradient, using only Hessian-vector products.
• liblinear: coordinate descent; binary only; supports L1. Unusually, it penalizes the intercept too, so set intercept_scaling high if that matters.
• sag and saga: stochastic average gradient, for very large n; saga also handles L1 and elastic net. Both need standardized features to converge quickly.

Worked example

Implementation

import torch
import torch.nn.functional as F

hours = torch.tensor([0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0]).unsqueeze(1)
passed = torch.tensor([0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1], dtype=torch.float32).unsqueeze(1)

model = torch.nn.Linear(1, 1)                 # z = w * hours + b, a logit
opt = torch.optim.LBFGS(model.parameters(), line_search_fn="strong_wolfe")

def closure():
    opt.zero_grad()
    loss = F.binary_cross_entropy_with_logits(model(hours), passed)
    loss.backward()
    return loss

for _ in range(10):
    opt.step(closure)
model.weight, model.bias                      # ≈ 1.4245 and -4.6296
torch.sigmoid(model(torch.tensor([[2.0], [4.0]])))   # ≈ 0.14 and 0.74

# K classes: K logits, then cross-entropy on the integer label
logits = torch.tensor([[2.0, 1.0, -1.0]])
F.cross_entropy(logits, torch.tensor([0]))    # 0.3490

import numpy as np
from sklearn.linear_model import LogisticRegression

X = np.array([0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6]).reshape(-1, 1)
y = np.array([0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1])

clf = LogisticRegression(C=np.inf).fit(X, y)  # C=inf: no penalty
clf.coef_, clf.intercept_                     # ≈ [[1.4245]], [-4.6296]
clf.predict_proba([[2.0], [4.0]])[:, 1]       # ≈ [0.14, 0.74]
clf.decision_function([[3.25]])               # ≈ [0.0], the boundary

LogisticRegression().fit(X, y).coef_          # C=1 by default: a smaller slope

Newton's method in plain Python:

import math

def sigmoid(z):
    return 1 / (1 + math.exp(-z)) if z >= 0 else math.exp(z) / (1 + math.exp(z))

def fit_logistic_1d(xs, ys, steps=8):
    w = b = 0.0
    for _ in range(steps):
        ps = [sigmoid(w * x + b) for x in xs]
        gw = sum((p - y) * x for p, x, y in zip(ps, xs, ys))
        gb = sum(p - y for p, y in zip(ps, ys))
        s = [p * (1 - p) for p in ps]
        hww = sum(si * x * x for si, x in zip(s, xs))
        hwb = sum(si * x for si, x in zip(s, xs))
        hbb = sum(s)
        det = hww * hbb - hwb * hwb             # solve the 2x2 Newton system
        w -= (hbb * gw - hwb * gb) / det
        b -= (hww * gb - hwb * gw) / det
    return w, b

fit_logistic_1d([0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6],
                [0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1])   # (1.4245, -4.6296)

Tradeoffs

• Linear boundary. The model can only separate classes with a hyperplane in its features. Interactions, polynomial terms or splines bend it, at the cost of more columns; past that point, trees or neural networks learn the features instead.
• Interpretable, with caveats. Coefficients are log odds ratios, holding the other features fixed, and carry every caution from linear regression: units, correlated features and confounders. Penalties also shrink them, so a regularized coefficient understates the effect it estimates.
• Probabilities. A fit by maximum likelihood gives fairly well calibrated probabilities, which most classifiers (SVMs, boosted trees) do not; check with a reliability plot before trusting them, especially under strong regularization or class imbalance.
• Imbalance. With rare positives, the intercept absorbs the base rate and accuracy is meaningless. Use class_weight, a moved threshold, and metrics like precision, recall or log loss rather than accuracy.
• Solver choice. Newton-type methods converge in a handful of steps but scale as d² memory and d³ time per step; first-order and stochastic methods scale to huge data but need scaled features and more iterations. A ConvergenceWarning usually means unscaled features, not too few iterations.

Check yourself

Connections

References

• Hastie, Tibshirani and Friedman. The Elements of Statistical Learning, 2nd ed., section 4.4. Springer, 2009.
• Bishop. Pattern Recognition and Machine Learning, sections 4.3.2–4.3.4 (logistic regression, IRLS, multiclass). Springer, 2006.
• Cox. The Regression Analysis of Binary Sequences. Journal of the Royal Statistical Society B, 1958.
• Nelder and Wedderburn. Generalized Linear Models. Journal of the Royal Statistical Society A, 1972.
• Albert and Anderson. On the Existence of Maximum Likelihood Estimates in Logistic Regression Models. Biometrika, 1984.
• Fan, Chang, Hsieh, Wang and Lin. LIBLINEAR: A Library for Large Linear Classification. JMLR, 2008.
• scikit-learn 1.9.1: _logistic.py, _linear_loss.py and _newton_solver.py.
• PyTorch v2.15.0: torch/nn/modules/linear.py and torch/nn/modules/loss.py.