Logistic regression

What it is and when you reach for it

Logistic regression is a linear model for the probability of a binary outcome: it scores an input with w·x + b and squashes that score through the sigmoid to a number in (0, 1) that is trained to be the actual probability of the positive class. It occupies a special position in the model hierarchy. Below it sits linear regression, which it repairs for classification (least squares on 0/1 labels penalizes confident correct predictions and produces values outside [0, 1]); above it sits every neural classifier, whose final layer is exactly a logistic (or multinomial logistic) regression on learned features. You reach for it first on any new classification problem: it trains in seconds, its coefficients can be read as log-odds contributions, its probabilities are honest enough to threshold against real costs, and it establishes the baseline any heavier model must beat. In large-scale production settings with sparse features it is often not just the baseline but the deployed model.

The math

From Bernoulli likelihood to cross-entropy

Model the label as a coin flip whose bias depends on x: y ∼ Bernoulli(p) with p = σ(w·x + b). For one example the likelihood is p^y(1 − p)^1−y, a single expression that reads p when y = 1 and 1 − p when y = 0. The dataset likelihood is the product over examples, so its negative log, which we minimize, is the sum

L(w, b) = −(1/n) Σᵢ [ yᵢ log pᵢ + (1 − yᵢ) log(1 − pᵢ) ]

and this is precisely binary cross-entropy. Nothing was chosen here; cross-entropy is not a loss someone picked because it worked, it is the unique consequence of saying "the label is Bernoulli and I will do maximum likelihood." That pedigree is why the minimizer of this loss is calibrated in-sample: at the optimum, among examples scored around 0.7, about 70 percent are positive, because any systematic miscalibration would leave likelihood on the table.

Sigmoid and logit are inverses

The sigmoid σ(z) = 1 / (1 + e^−z) maps the real line to (0, 1); its inverse is the logit, log(p / (1 − p)), the log of the odds. So the model statement "p = σ(w·x + b)" and the statement "the log-odds of the positive class are linear in x" are the same sentence read in opposite directions, and the second one explains the coefficients: increasing feature j by one unit adds w_j to the log-odds, multiplying the odds by e^w_j, regardless of where you started. This is the interpretation epidemiologists and credit modelers rely on. Two identities used constantly below: σ(−z) = 1 − σ(z), and dσ/dz = σ(z)(1 − σ(z)), which you can verify by differentiating (1 + e^−z)⁻¹ directly.

The gradient, fully derived

Write z = w·x + b and p = σ(z), and take one example's loss ℓ = −y log p − (1 − y) log(1 − p). Chain rule through p first:

∂ℓ/∂p = −y/p + (1 − y)/(1 − p)

∂ℓ/∂z = ∂ℓ/∂p · dp/dz
       = [ −y/p + (1 − y)/(1 − p) ] · p(1 − p)
       = −y(1 − p) + (1 − y)p
       = −y + yp + p − yp
       = p − y

Every messy term cancels and the derivative of the loss with respect to the score is simply p − y, the prediction error. The last hop is linear: ∂z/∂w = x and ∂z/∂b = 1, so for the full dataset, stacking the examples into X,

∂L/∂w = (1/n) Xᵀ (p − y)        ∂L/∂b = mean(p − y)

This is the same "error times input" form as linear regression, which is not a coincidence: both are generalized linear models whose link function is matched to their likelihood, and for every such matched pair the gradient collapses to residual times features. It is also the form that survives into deep learning: softmax plus cross-entropy at the top of any network backpropagates exactly p − y into the logits, one reason that pairing is universal (the mechanics of the softmax half are on the softmax page).

A worked numeric step

Take w = (1, −1), b = 0, one positive example x = (2, 1), y = 1, learning rate 0.5. The score is z = 1·2 + (−1)·1 = 1, so p = σ(1) ≈ 0.731 and the loss is −log 0.731 ≈ 0.313. The error is p − y ≈ −0.269, so the gradients are ∂L/∂w = −0.269 · (2, 1) = (−0.538, −0.269) and ∂L/∂b = −0.269. One step: w ← (1.269, −0.866), b ← 0.135. The new score is 1.269·2 − 0.866·1 + 0.135 ≈ 1.807, p rises to σ(1.807) ≈ 0.859, and the loss falls to about 0.152. The mislabeled-side pull is proportional to how wrong the probability was, and a perfectly confident correct prediction (p = y) contributes exactly zero gradient, which is the sense in which the loss knows when it is done.

Why the loss is convex

Per example, write the loss in terms of the score using log(1 + e^z) = logaddexp(0, z): ℓ(z) = logaddexp(0, z) − yz. Its first derivative is p − y as derived, and its second derivative is dp/dz = p(1 − p), which is positive everywhere, so ℓ is convex in z. The score is an affine function of (w, b), and composing a convex function with an affine map preserves convexity, so each example's loss is convex in the parameters and the average keeps it that way. In matrix form the Hessian is (1/n) Xᵀ S X with S = diag(p_i(1 − p_i)), which is positive semidefinite by construction. Convexity means there are no local minima to fall into: every method from plain gradient descent to L-BFGS lands at the same global optimum, differing only in how fast, and that is why comparing your implementation against a library is a sharp test rather than a shrug about different basins.

Regularization

The semidefinite Hessian hides one real failure mode: if the data is linearly separable, the likelihood can always be improved by scaling w up, pushing every p toward 0 or 1, so the weights diverge and the optimum does not exist. Adding an L2 penalty (λ/2)‖w‖² makes the Hessian positive definite, the objective strictly convex, and the optimum finite and unique; statistically it is a Gaussian prior on the weights. L1 instead drives coordinates exactly to zero and doubles as feature selection, at the cost of a non-smooth objective that needs solvers built for it. As with the SVM, the bias should be exempt: shrinking b toward zero encodes a belief that classes are balanced that you probably do not hold. One convention trap: scikit-learn parameterizes regularization by C = 1/(nλ) in the mean-loss convention used here, and larger C means less regularization.

Implementation, twice

The first pair is the standard formulation: PyTorch with nn.Linear and BCEWithLogitsLoss, which fuses the sigmoid into the loss with the logaddexp identity so large scores never overflow, and JAX with the same stable loss written out explicitly, value_and_grad, and the whole training loop compiled as a single lax.scan. One optimizer detail worth copying: weight decay is applied to the weights through a parameter group and explicitly not to the bias.

import torch
from torch import nn

def train_logreg(X, y, lam=1e-3, lr=0.5, epochs=300, batch_size=None):
    """Logistic regression via nn.Linear + BCEWithLogitsLoss.

    X: (n, d) float tensor.  y: (n,) float tensor of 0/1 labels.
    batch_size=None trains full batch; an int trains minibatch SGD.
    BCEWithLogitsLoss takes raw scores: the sigmoid is fused into
    the loss for numerical stability, so the model has no sigmoid.
    """
    n, d = X.shape
    model = nn.Linear(d, 1)
    loss_fn = nn.BCEWithLogitsLoss()
    # decay the weights, never the bias
    opt = torch.optim.SGD(
        [{"params": [model.weight], "weight_decay": lam},
         {"params": [model.bias], "weight_decay": 0.0}],
        lr=lr,
    )

    def step(xb, yb):
        opt.zero_grad()
        loss = loss_fn(model(xb).squeeze(-1), yb)
        loss.backward()
        opt.step()
        return loss.item()

    losses = []
    for _ in range(epochs):
        if batch_size is None:
            losses.append(step(X, y))
        else:
            perm = torch.randperm(n)
            for i in range(0, n, batch_size):
                idx = perm[i : i + batch_size]
                losses.append(step(X[idx], y[idx]))
    return model, losses

def predict_proba(model, X):
    with torch.no_grad():
        return torch.sigmoid(model(X).squeeze(-1))

import jax
import jax.numpy as jnp

def loss_fn(params, X, y, lam):
    """Stable BCE on logits: logaddexp(0, z) - y*z == -[y log p + (1-y) log(1-p)]."""
    w, b = params
    z = X @ w + b
    nll = jnp.mean(jnp.logaddexp(0.0, z) - y * z)
    return nll + 0.5 * lam * jnp.dot(w, w)   # bias never regularized

def train_logreg(X, y, lam=1e-3, lr=0.5, epochs=300):
    """Full-batch gradient descent; the whole loop is one compiled scan."""
    grad_fn = jax.value_and_grad(loss_fn)

    def step(params, _):
        loss, g = grad_fn(params, X, y, lam)
        params = jax.tree.map(lambda p, gi: p - lr * gi, params, g)
        return params, loss

    params0 = (jnp.zeros(X.shape[1]), jnp.array(0.0))
    params, losses = jax.lax.scan(step, params0, None, length=epochs)
    return params, losses

def predict_proba(params, X):
    w, b = params
    return jax.nn.sigmoid(X @ w + b)

The second pair strips the abstraction away. The PyTorch tab implements the derived gradient Xᵀ(p − y)/n by hand, no autograd, which is the version to write when you want to prove to yourself the calculus above is the whole story. The JAX tab does minibatch training the JAX way: reshuffle each epoch with an explicit PRNG key, reshape the epoch into a (steps, batch, d) tensor, and scan over the batches inside a scan over epochs, so the entire multi-epoch minibatch schedule compiles to one XLA program with no Python in the hot path.

import torch

def train_logreg_manual(X, y, lam=1e-3, lr=0.5, epochs=300):
    """The derivation, verbatim: grad_w = X^T (p - y) / n + lam * w.

    No autograd anywhere; compare against the nn version to confirm
    they converge to the same (w, b) on the same data.
    """
    n, d = X.shape
    w = torch.zeros(d)
    b = torch.zeros(())
    losses = []
    for _ in range(epochs):
        p = torch.sigmoid(X @ w + b)       # (n,)
        err = p - y                        # exactly dL/dz from the math
        w -= lr * (X.T @ err / n + lam * w)
        b -= lr * err.mean()
        # loss via the stable identity, for monitoring only
        z = X @ w + b
        nll = (torch.clamp(z, min=0) - y * z + torch.log1p(torch.exp(-z.abs()))).mean()
        losses.append((nll + 0.5 * lam * w.dot(w)).item())
    return w, b, losses

import jax
import jax.numpy as jnp

def train_logreg_minibatch(X, y, key, lam=1e-3, lr=0.5,
                           epochs=50, batch=256):
    """Minibatch SGD, fully compiled: scan over epochs, scan over batches.

    Drops the ragged tail each epoch (n % batch examples), which the
    fresh shuffle makes harmless when n is much larger than batch.
    """
    n, d = X.shape
    steps = n // batch
    grad_fn = jax.value_and_grad(loss_fn)   # loss_fn from the previous block

    def batch_step(params, xy):
        Xb, yb = xy
        loss, g = grad_fn(params, Xb, yb, lam)
        params = jax.tree.map(lambda p, gi: p - lr * gi, params, g)
        return params, loss

    def epoch_step(params, key):
        perm = jax.random.permutation(key, n)[: steps * batch]
        Xb = X[perm].reshape(steps, batch, d)
        yb = y[perm].reshape(steps, batch)
        params, losses = jax.lax.scan(batch_step, params, (Xb, yb))
        return params, losses.mean()

    params0 = (jnp.zeros(d), jnp.array(0.0))
    keys = jax.random.split(key, epochs)
    params, epoch_losses = jax.lax.scan(epoch_step, params0, keys)
    return params, epoch_losses

The multinomial extension

Nothing structural changes with K classes. The weight vector becomes a matrix W with one column of scores per class, sigmoid becomes softmax, the Bernoulli likelihood becomes categorical, and the gradient of the loss with respect to the logits becomes p − onehot(y), the same error form with the same cancellation. Binary logistic regression is the K = 2 case with one redundant column removed.

import torch
from torch import nn

def train_multinomial(X, y, num_classes, lam=1e-3, lr=0.5, epochs=300):
    """y: (n,) int64 class indices.  CrossEntropyLoss = log_softmax + NLL,
    fused and stable, so the model outputs raw logits."""
    model = nn.Linear(X.shape[1], num_classes)
    loss_fn = nn.CrossEntropyLoss()
    opt = torch.optim.SGD(
        [{"params": [model.weight], "weight_decay": lam},
         {"params": [model.bias], "weight_decay": 0.0}],
        lr=lr,
    )
    for _ in range(epochs):
        opt.zero_grad()
        loss = loss_fn(model(X), y)
        loss.backward()
        opt.step()
    return model

import jax
import jax.numpy as jnp

def multinomial_loss(params, X, y, lam):
    """y: (n,) int class indices.  log_softmax keeps everything stable."""
    W, b = params                            # W: (d, K), b: (K,)
    logits = X @ W + b
    logp = jax.nn.log_softmax(logits, axis=-1)
    nll = -logp[jnp.arange(X.shape[0]), y].mean()
    return nll + 0.5 * lam * jnp.sum(W * W)

def train_multinomial(X, y, num_classes, lam=1e-3, lr=0.5, epochs=300):
    grad_fn = jax.value_and_grad(multinomial_loss)

    def step(params, _):
        loss, g = grad_fn(params, X, y, lam)
        params = jax.tree.map(lambda p, gi: p - lr * gi, params, g)
        return params, loss

    params0 = (jnp.zeros((X.shape[1], num_classes)),
               jnp.zeros(num_classes))
    params, losses = jax.lax.scan(step, params0, None, length=epochs)
    return params, losses

Using it on a real shape of problem

A representative shape: 5,000 examples, 50 features, generated from an actual logistic model so we know the truth the fit should recover. Expect the full-batch loss to decay smoothly toward the Bayes floor (it will not reach zero, because the labels are genuinely noisy near the boundary), the minibatch loss to follow the same trend with jitter, and accuracy to settle in the 80s for this noise level. Because the objective is strictly convex under L2, the manual-gradient version, the nn version, and the JAX versions must all agree on the final weights to several decimals; if they do not, one of them has a bug, and that agreement is the first test worth running. Exact figures depend on seed and hardware, the shapes do not.

import torch

torch.manual_seed(0)
n, d = 5000, 50
X = torch.randn(n, d)
w_true = torch.randn(d) / d ** 0.5
y = (torch.rand(n) < torch.sigmoid(X @ w_true * 3)).float()

model, losses = train_logreg(X, y, lam=1e-3, lr=0.5,
                             epochs=100, batch_size=256)
p = predict_proba(model, X)
acc = ((p > 0.5).float() == y).float().mean()
print(f"loss {losses[0]:.3f} -> {losses[-1]:.3f}, acc {acc:.3f}")
# smooth-ish decay with minibatch jitter; acc ~0.81 (seed-dependent)

import jax
import jax.numpy as jnp

key = jax.random.PRNGKey(0)
kx, kw, ky, kt = jax.random.split(key, 4)
n, d = 5000, 50
X = jax.random.normal(kx, (n, d))
w_true = jax.random.normal(kw, (d,)) / d ** 0.5
y = (jax.random.uniform(ky, (n,)) < jax.nn.sigmoid(X @ w_true * 3)
     ).astype(jnp.float32)

params, epoch_losses = train_logreg_minibatch(X, y, kt, lam=1e-3,
                                              lr=0.5, epochs=100)
p = predict_proba(params, X)
acc = jnp.mean((p > 0.5).astype(jnp.float32) == y)
print(f"loss {epoch_losses[0]:.3f} -> {epoch_losses[-1]:.3f}, acc {acc:.3f}")
# smooth-ish decay with minibatch jitter; acc ~0.86 (seed-dependent)

Applications

The flagship production application is click-through rate prediction. For years the ad systems at Google, Facebook, and Microsoft ran logistic regression over billions of sparse cross-features, because the model trains online one example at a time, its probabilities feed directly into expected-value bidding where calibration is money, and the p − y gradient costs almost nothing per impression. Google's "Ad Click Prediction: a View from the Trenches" paper describes the FTRL-Proximal variant they ran at that scale, and the ad click prediction system design page on this site walks through where the model sits in the full serving pipeline. Even after deep models took over ranking, logistic regression persists as the calibration layer and as the shadow baseline that new models must beat by enough to justify their serving cost.

Beyond ads, logistic regression is the standard model wherever the coefficients themselves are the product: credit scoring, where regulators require explainable adverse-action reasons; epidemiology and clinical research, where odds ratios per risk factor are the published result; and churn and conversion modeling, where a calibrated probability is fed into an expected value decision. And there is a sense in which every deep classifier ships a logistic regression: the final linear-plus-softmax layer of any neural network is exactly multinomial logistic regression on the features the rest of the network learned, which is why linear probes, the standard tool for asking what a representation contains, are just this model fit on frozen activations, and why everything on this page about calibration, convexity, and p − y transfers upward intact.

Against the real libraries

The production reference is scikit-learn's LogisticRegression, and its solver parameter is a tour of convex optimization trade-offs. The default, lbfgs, is a quasi-Newton method that builds a low-rank curvature estimate from gradient history; on dense small-to-medium problems with L2 it converges in tens of iterations to precision that first-order methods reach slowly or never, and it should be your default too. liblinear wraps the same LIBLINEAR coordinate-descent library the linear SVM uses; it supports L1, is very fast on small and sparse problems, but is inherently one-vs-rest for multiclass and, a real gotcha, regularizes the intercept. saga is the incremental-gradient option: it processes one example at a time with variance-reduced updates, is the only solver supporting elastic net, and wins when n is in the millions and sparse, where forming full gradients is the bottleneck. The practical rule: lbfgs unless you need L1 or elastic net or have huge sparse data, then saga; liblinear for small sparse L1 problems.

What the library adds over this page's implementations is precision and bookkeeping rather than different mathematics: line searches and convergence tolerances that squeeze the convex problem to its floor, class weighting, warm starts across a regularization path, and robust multiclass handling. The from-scratch version is enough when the model is a component of a larger differentiable system, when you need a custom loss or training schedule, or when you want streaming updates the way the CTR systems do. Verification is sharp because the optimum is unique: fit LogisticRegression(C=1/(n*lam), solver="lbfgs", tol=1e-8) on the same standardized data, run your own trainer to convergence, and the weight vectors should match to three or four decimals; any stubborn gap traces to the intercept convention, unconverged SGD, or a mismatched C.

The second library is statsmodels, and it answers a different question. scikit-learn tells you what the model predicts; statsmodels' Logit tells you what the model means: standard errors on each coefficient, z-statistics and p-values, confidence intervals, and likelihood-ratio comparisons between nested models, all from the same fit. If the deliverable is "feature j is associated with the outcome, odds ratio 1.4, 95 percent CI [1.2, 1.7]", that is statsmodels' job, and no amount of scikit-learn will produce it. The two fit the same objective (mind that statsmodels does not regularize by default while scikit-learn always does, so set penalty=None to compare), and running both on one dataset is a good habit: prediction metrics from one, inference from the other.

Traps and misconceptions

Sigmoid then BCELoss. Computing the probability and passing it to a log loss re-exposes the overflow that the fused version was built to avoid: for z = −100 the probability underflows to 0 and log(0) is −inf. Always give raw logits to BCEWithLogitsLoss or use the logaddexp form; the JAX implementations above never materialize p during training at all.

Forgetting that separable data breaks the model. On linearly separable data the unregularized weights diverge and the probabilities saturate to 0 and 1. If your training loss keeps creeping down forever and the weight norm keeps growing, you are watching this happen; any nonzero L2 fixes it, which is one reason scikit-learn regularizes by default.

Reading coefficients off unscaled features. Coefficient magnitude mixes the effect size with the feature's unit; a weight on "income in dollars" will look microscopic next to one on "age in decades" regardless of importance. Standardize before comparing coefficients, and remember that even then the coefficients describe association under the model, not causation.

Accuracy as the training signal. The model is fit to log loss and its output is a probability; thresholding at 0.5 to report accuracy throws away exactly the calibration that makes it valuable and misleads badly under class imbalance. Evaluate with log loss or calibration curves plus a threshold-free metric, and choose the operating threshold from the actual costs of the two error types.

Expecting p-values from scikit-learn. There is no significance testing anywhere in sklearn's LogisticRegression, and improvising it from the coefficients is easy to get wrong because the default L2 penalty biases the estimates. When inference is the point, fit unpenalized in statsmodels, which computes the standard errors from the Hessian for you.