Logistic Regression, in depth: log-odds, log loss, Newton and gradient solvers, separation, regularisation and calibration

The model: linear log-odds

Start from odds. If p is the probability of the positive class, the odds are p/(1 - p) and the log-odds, or logit, is log(p/(1 - p)), which ranges over the whole real line. Logistic regression assumes that the log-odds are linear in the features: log(p/(1 - p)) = w.x + b = z. Solve for p and you get p = 1/(1 + e^-z), the sigmoid. So the sigmoid is not an arbitrary squashing function: it is the inverse of the assumption that evidence adds up on the log-odds scale.

The decision boundary, where p = 0.5, is the hyperplane w.x + b = 0. The model is linear in its features, and you get curved boundaries only by adding features such as interactions, splines or embeddings. The magnitude of z is a measure of confidence: at z = 2, p = 0.881, and at z = 4, p = 0.982.

The loss and its gradient

Treat each label as a Bernoulli draw with probability p_i. The negative log-likelihood of the data, averaged over the n examples, is the log loss, or binary cross-entropy: L = -(1/n) sum of [y_i log p_i + (1 - y_i) log(1 - p_i)]. Cross-entropy, derived walks through the same loss in the multiclass setting used by language models.

Two properties make it pleasant to optimise. First, the derivative of the sigmoid is p(1 - p), which cancels neatly with the loss, so the gradient is just X^T(p - y)/n: the prediction error weighted by the features. It has the same form as least squares in linear regression. Second, the loss is convex in w, with Hessian X^TSX/n where S = diag(p_i(1 - p_i)). Any local minimum is therefore global, and Newton's method converges in a handful of steps.

Compute the loss in a numerically stable way. Evaluating log(sigmoid(z)) literally gives log(0) = -infinity for large negative z. Use the identity that the per-example loss equals log(1 + e^z) - yz, and evaluate log(1 + e^z) with np.logaddexp(0, z).

Two solvers in NumPy

import numpy as np
from scipy.special import expit          # overflow-safe sigmoid

def loss_and_grad(w, X, y, lam):
    """X has a leading column of ones; w[0] is the intercept (not penalised)."""
    z = X @ w
    p = expit(z)
    loss = np.mean(np.logaddexp(0, z) - y * z) + 0.5 * lam * w[1:] @ w[1:]
    g = X.T @ (p - y) / len(y)
    g[1:] += lam * w[1:]
    return loss, g, p

def fit_gd(X, y, lam=1e-2, lr=0.5, steps=2000):
    X = np.c_[np.ones(len(X)), X]
    w = np.zeros(X.shape[1])
    for _ in range(steps):
        _, g, _ = loss_and_grad(w, X, y, lam)
        w -= lr * g
    return w

def fit_newton(X, y, lam=1e-2, iters=50, tol=1e-10):   # also called IRLS
    X = np.c_[np.ones(len(X)), X]
    n, d = X.shape
    w = np.zeros(d)
    R = lam * np.eye(d)
    R[0, 0] = 0.0
    for _ in range(iters):
        _, g, p = loss_and_grad(w, X, y, lam)
        H = (X * (p * (1 - p))[:, None]).T @ X / n + R
        step = np.linalg.solve(H, g)
        w -= step
        if np.abs(step).max() < tol:
            break
    return w

Gradient descent costs O(nd) per step and scales to any n, especially with mini-batches. The trade-offs of step size and convergence are covered in gradient descent theory. Newton's method builds and solves a d by d system, O(nd² + d³) per step, but typically needs fewer than 20 steps. It is the right tool for up to a few thousand features. Each Newton step is a weighted least-squares solve, which is why statisticians call it iteratively reweighted least squares (IRLS). Quasi-Newton L-BFGS sits between the two, and it is what most libraries use by default.

Worked example, and what separation does

Take four one-feature examples: x = 1, 2, 3, 4 with labels y = 0, 0, 1, 1. Start at w = 0, b = 0. Every prediction is p = 0.5, so the loss is log 2 = 0.6931. The gradient with respect to w is the mean of (p - y)x = (0.5 + 1.0 - 1.5 - 2.0)/4 = -0.5. The gradient with respect to b is the mean of (p - y) = 0. One gradient step with learning rate 1 gives w = 0.5, b = 0, and the loss falls to 0.6539. The boundary is still at x = 0, so all four points are predicted positive. Later steps push b negative until the boundary reaches x = 2.5, the midpoint between the classes.

Now the trap. This data is separable: a threshold at 2.5 classifies every point correctly. Any boundary at 2.5 can be made sharper by scaling w and b up together, which pushes every p closer to its label and lowers the loss. So the unregularised loss has no minimum. Running the Newton code above with lam = 0 confirms it. The weight goes 1.6, 2.84, 4.53, 6.58, and then grows by about 2 per iteration. By iteration 38 it reaches w = 74.7, the loss is about 10^-17, and the Hessian becomes numerically singular, so the solve fails. Gradient descent fails more slowly: w = 5.8 after 1,000 steps and 14.9 after 100,000, still rising.

With a penalty the problem has a unique finite solution. The boundary stays at 2.5, and lambda now sets how confident the probabilities are. In real data with many features, perfect or quasi-perfect separation is common: one rare category that is always positive is enough. This is why regularisation should be on by default.

Regularisation and scikit-learn

L2 regularisation adds lambda/2 times |w|². It shrinks all coefficients smoothly and keeps correlated features together. L1 adds lambda times |w|₁, which drives some coefficients exactly to zero and so selects features. Elastic net mixes the two. Regularisation compares them across model families. Always standardise features before penalising, or the penalty falls unevenly on features measured in large and small units.

In scikit-learn, LogisticRegression is regularised by default. It minimises C times the summed log loss plus half the squared L2 norm, with C=1.0, so smaller C means stronger regularisation. For n examples, the page's lambda corresponds to 1/(C n). Since version 1.8 the penalty argument is deprecated (removal is scheduled for 1.10), and the penalty type is set with l1_ratio: 0.0 for L2 (the new default), 1.0 for L1, and values in between for elastic net, which needs the saga solver. Code that must run on older versions should keep using penalty. The multi_class argument is gone from the 1.8 signature, and multiclass problems are fitted as a multinomial model.

from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegressionCV

model = make_pipeline(
    StandardScaler(),
    LogisticRegressionCV(Cs=10, cv=5, scoring="neg_log_loss", max_iter=1000),
)
model.fit(X_train, y_train)
p_valid = model.predict_proba(X_valid)[:, 1]

Score the search on log loss, not accuracy, when you need good probabilities. Raise max_iter if you see convergence warnings, rather than ignoring them. On 1.8, LogisticRegressionCV also warns that l1_ratios=None is deprecated; pass l1_ratios=(0.0,) there to keep plain L2. One solver detail: liblinear penalises the intercept, while lbfgs does not.

Reading coefficients

Each coefficient is a change in log-odds per unit of its feature, holding the others fixed. Exponentiate it to get an odds ratio: with standardised inputs, a coefficient of 0.7 means that one standard deviation more of that feature multiplies the odds by e^0.7 = 2.01. Three caveats apply before anyone reads causation into the table.

• Coefficients are shrunk by the penalty, so their size reflects regularisation as well as the data.
• With correlated features, credit is split arbitrarily between them, and signs can flip from one refit to the next. Check stability across bootstrap refits before reporting individual coefficients.
• "Holding the others fixed" may describe a combination that never occurs in the data.

Probabilities, thresholds and imbalance

Logistic regression fitted by maximum likelihood is often reasonably well calibrated on data like its training set. Several common practices break that, and the threshold for a decision should come from costs, not from 0.5.

• Choosing a threshold. If a false negative costs c_FN and a false positive costs c_FP, act when p >= c_FP/(c_FP + c_FN), assuming the probabilities are calibrated.
• class_weight="balanced" reweights examples by n/(classes x class count). It can improve recall on a rare class, but it deliberately inflates its probabilities, so predicted rates no longer match reality.
• Downsampling negatives. If you keep a fraction r of negatives, the fitted log-odds are too high by -log r. Add log r to the intercept to recover population probabilities. This works because only the intercept absorbs the sampling change.
• Strong regularisation pulls probabilities toward the base rate. Check reliability curves and apply Platt or isotonic recalibration as described in model calibration.

More than two classes

For K classes, multinomial logistic regression keeps one weight vector per class and replaces the sigmoid with the softmax: p_k = e^z_k divided by the sum of e^z_j. The gradient keeps the same form, X^T(P - Y)/n with one-hot Y. One class's weights are redundant, since adding the same vector to every class leaves the softmax unchanged, and regularisation is what pins them down. One-vs-rest, meaning K separate binary models, is the alternative when classes are not mutually exclusive. In scikit-learn, wrap the estimator in OneVsRestClassifier.

Operational guidance

• Ship the pipeline, not the weights. Scaling, encoding, clipping and the intercept correction all belong to the model artefact. Version them together.
• Unseen categories. Configure the encoder to ignore unknown categories, which then contribute zero, rather than crash at serving time.
• Monitor the mean predicted probability against the observed positive rate, the distribution of z, and the log loss on labelled feedback. Calibration drift usually appears before accuracy drops.
• Serving is cheap. A prediction is one dot product, O(d), and it is fast enough for request paths, in-database scoring and edge devices.
• Retrain on a schedule, and alert when a refit moves the coefficients sharply. That signals drift or a broken feature.

Failure modes

• Separation without a penalty. Weights diverge, and probabilities of 0 and 1 appear, as in the worked example.
• Unscaled features, which slow L-BFGS and distort the penalty.
• Leakage features, such as a field filled in after the outcome, which produce suspiciously near-perfect separation.
• Treating balanced-weight probabilities as real rates in pricing or capacity planning.
• Expecting non-linear boundaries. If the signal is an interaction or a threshold effect, add the feature or use a tree model.

Trade-offs

What to do next

1. Run fit_newton on the four-point example with lam = 0, then 0.01, and reproduce the table.
2. On your own data, build the scaler plus LogisticRegressionCV pipeline scored on log loss, and record the chosen C.
3. Plot a reliability curve on validation data, and recalibrate if the points leave the diagonal.
4. Set the decision threshold from your false-positive and false-negative costs.
5. If you downsampled or used class weights, correct the intercept or recalibrate before reporting probabilities.
6. Check coefficient stability with bootstrap refits before explaining any single coefficient to stakeholders.