Logistic Regression as Conditional Bernoulli Likelihood

A linear model that predicts a real number extends to binary classification only if something maps the real line into the unit interval. The sigmoid is the standard choice, and it is usually introduced as an activation function borrowed from elsewhere. It is not borrowed; it is the inverse of the log-odds, and the log-odds is what makes the model a linear function in the parameters of a Bernoulli likelihood. Starting from the likelihood rather than from the sigmoid explains the loss, the gradient, the convexity, and the failure in the separable case as four views of one construction.

From log-odds to a probability

For a Bernoulli trial with success probability $\pi$, the odds are $\pi/(1-\pi)$, a positive number on the whole positive half-line. Its logarithm maps $(0,1)$ onto $\mathbb{R}$ bijectively, and it is the natural parameter of the Bernoulli family in exponential-family form: writing the distribution as $\exp(\eta y - A(\eta))$ with $y \in {0,1}$ gives $\eta = \log(\pi/(1-\pi))$ and log-partition $A(\eta) = \log(1 + e^{\eta})$.

Assuming the natural parameter is linear in the covariates, $\eta = x^\top \beta$, is the modelling choice that defines logistic regression. Inverting gives

$$\pi(x) = \frac{1}{1 + e^{-x^\top \beta}},$$

and the derivative of the log-partition, $A’(\eta) = \sigma(\eta)$, is the mean. The model is therefore a two-line construction: pick the exponential family, then declare its natural parameter linear. The probabilistic reading of $\beta_j$ is a coefficient on the change in log-odds per unit of $x_j$, holding the other covariates fixed, an additive statement about odds on a multiplicative scale, not about probabilities. Doubling an odds ratio of 1.5 to 3.0 corresponds to the same additive step in coefficients, which is why reporting a coefficient as a log-odds increment is the only interpretation that composes across features.

The intercept has a separate role. When the class balance is skewed, the intercept absorbs the log prior odds, and its magnitude can be large even when no covariate is informative: for a 99-to-1 split, the maximum-likelihood intercept alone is about −4.6, the natural log of the odds against. Regularising the intercept along with the coefficients shifts predictions away from the base rate, which is the usual explanation for coefficient-level regularisation degrading calibration on imbalanced data.

The negative log-likelihood is the cross-entropy

For independent observations, the likelihood is $\prod_i \pi_i^{y_i}(1-\pi_i)^{1-y_i}$, and its negative logarithm is

$$-\ell(\beta) = \sum_i \left[ \log(1 + e^{x_i^\top \beta}) - y_i x_i^\top \beta \right].$$

The same expression is the binary cross-entropy between the label and the predicted probability, which is not a coincidence: cross-entropy is the expected negative log-likelihood of the data distribution under the model, so identical objectives appear from two different derivations. Each term is the log-partition of the Bernoulli family minus the sufficient statistic times its coefficient, the canonical form for a linear exponential family with a linear predictor.

Numerical evaluation needs care because the two terms can overflow independently. For large negative $x^\top\beta$, $e^{x^\top\beta}$ underflows to zero and $\log(1+e^{x^\top\beta})$ is evaluated as $\log 0$ if the sum is computed naively; for large positive values the exponential itself overflows. The stable form is

$$\log(1 + e^{z}) = \max(z, 0) + \log\left(1 + e^{-|z|}\right),$$

which never exponentiates a positive number. It costs one abs, one max, and one log1p, and it is what every production implementation uses, because a loss that returns infinity for a confident correct prediction destroys the gradient for the entire batch.

import numpy as np

def bce_stable(Xb, y):
    z = Xb                       # linear predictor, shape (n,)
    # log(1 + exp(z)) - y*z, rearranged so exp() never sees a positive argument
    loss = np.maximum(z, 0.0) + np.log1p(np.exp(-np.abs(z))) - y * z
    return loss.mean()

Note that the instability is asymmetric in sign: a confidently correct prediction drives $e^{-|z|}$ to zero and the term to zero smoothly, while a confidently wrong prediction grows linearly in $|z|$ with no branch. That linear tail, rather than the flat tail of the squared loss, is what keeps a misclassified point contributing a gradient.

The gradient is a weighted residual sum

Differentiating each term with respect to $\beta$ uses $d A / d\eta = \sigma(\eta)$ and gives

$$\nabla_\beta \left[-\ell\right] = \sum_i \left(\sigma(x_i^\top\beta) - y_i\right) x_i = X^\top (\hat\pi - y).$$

The gradient vanishes when the predicted probabilities are orthogonal to the columns of the design, which is the score-equation condition for the exponential family: the sufficient statistic $X^\top y$ is matched by its expectation under the model. Least squares has the same shape, $X^\top(y - X\hat\beta) = 0$, with different residual meanings; the geometric derivation of that case is in least squares as orthogonal projection, and the structural similarity is why the two models share a fitting routine for the Hessian.

The Hessian is $X^\top W X$ with $W = \operatorname{diag}(\hat\pi_i(1-\hat\pi_i))$. Since each weight is positive and bounded above by $\tfrac14$, $W$ is positive definite, so the Hessian is positive semidefinite for every $\beta$ and the problem is convex. That is a strong claim: any local minimum is global, there are no saddle points to escape, and the only non-uniqueness comes from the null space of $X$, which is the same rank-deficiency story as in the linear case. Convexity is a consequence of $A”(\eta) = \sigma(\eta)(1-\sigma(\eta)) > 0$, the variance of the Bernoulli family; it is not an accident of how the loss was written.

The bound $w_i \le 1/4$ also bounds the curvature. Curvature is largest where the model is uncertain and the weights peak, and it collapses toward zero wherever predictions are saturated, which is exactly where the loss is nearly flat and the gradient is nearly zero. A first-order optimiser slows down precisely in the saturated region, and the constrained step sizes that make it converge are set by the $1/4$ bound on the maximum eigenvalue of $X^\top W X$, which is at most $\tfrac14 \lambda_{\max}(X^\top X)$.

Separable data drives the coefficients to infinity

If a hyperplane separates the classes exactly, the likelihood has no maximum at a finite $\beta$. Scaling any separating coefficient vector by $t$ multiplies every log-odds by $t$, drives each fitted probability toward 0 or 1 in the correct direction, and reduces the loss monotonically toward zero. The infimum is 0 and it is never attained. The optimiser therefore keeps moving, and the coefficients diverge.

This is not a numerical artifact to be worked around; it is the correct likelihood-answer to a question that has no interior solution. Three responses are common, and they differ in what they assume. Adding a penalty, an L2 term on the coefficients, or a Gaussian prior, since the two are the same computation, gives a finite optimum and is treated in ridge shrinkage in ill-conditioned designs. Capping the number of iterations truncates the divergence at an arbitrary point determined by the starting values and the step rule. Early stopping does the same thing but with the truncation point chosen by held-out loss, which at least makes the stopping rule a function of the data.

A partial-separation case is subtler. When a single feature separates a subset of the points, only the coefficient on that feature diverges, and the remaining coefficients converge to finite limits. Iteration limits hide this: the converging coefficients look stable to many digits while one of them marches outward, so an iteration limit that is adequate for a well-posed problem can silently truncate a coefficient that was still moving.

from scipy.optimize import minimize
import numpy as np

def fit_logistic(X, y, l2=1e-3, max_iter=100):
    n, p = X.shape
    def obj(beta):
        z = X @ beta
        loss = np.mean(np.maximum(z, 0) + np.log1p(np.exp(-np.abs(z))) - y * z)
        return loss + 0.5 * l2 * beta[1:] @ beta[1:]   # intercept unpenalised
    def grad(beta):
        pi = 1.0 / (1.0 + np.exp(-(X @ beta)))
        r = pi - y
        g = X.T @ r / n
        g[1:] += l2 * beta[1:]
        return g
    res = minimize(obj, np.zeros(p), jac=grad, method='L-BFGS-B',
                   options={'maxiter': max_iter})
    return res.x

The residual vector pi - y is also the right quantity to inspect after fitting: it is the per-observation contribution to the score, and its largest entries identify the observations the model cannot explain, which is the starting point for leverage and influence diagnostics.

What the coefficient does not tell you

The estimate is consistent under the assumed link, but the standard errors depend on the variance function $V(\pi) = \pi(1-\pi)$, which depends on the fitted probabilities. When the model is misspecified by omitting a covariate that interacts with an included one, the coefficient on the included feature converges not to zero but to a weighted average of the true conditional effects, because logistic regression is not collapsible, omitting an interacting variable biases the retained coefficient even when the omitted variable is independent of the retained one. The linear-model analogue has no such bias, which makes the difference easy to forget.

Perfect prediction of a subgroup also makes the corresponding coefficients unstable without any separation across the full sample, and a small number of influential observations can move the intercept substantially because the weights $\pi(1-\pi)$ are smallest exactly for those points, so their influence is damped, the opposite of the least-squares case where high leverage amplifies influence. Extending the model to more than two classes replaces the Bernoulli with a categorical family, and the choice between the two standard parameterisations is the subject of softmax versus one-vs-rest probabilities.