Newton and IRLS Fitting for Logistic Regression

The logistic negative log-likelihood has no closed-form minimiser, but it does have a closed-form Hessian, and the function it defines is convex with a curvature that is cheap to evaluate. Newton’s method exploits both. For a generalised linear model with a canonical link, the resulting algorithm turns out to be exactly a sequence of weighted linear regressions, which is why it is written as iteratively reweighted least squares rather than as a second-order optimiser. The identity is not a coincidence; it follows from the Hessian of a canonical-link GLM being the weighted cross-product matrix of the design.

The Newton step for a convex objective

Newton’s method replaces the objective at the current iterate with its second-order Taylor expansion and jumps to that quadratic’s minimiser. Writing the negative log-likelihood as $f(\beta)$ with gradient $g$ and Hessian $H$,

$$\beta_{k+1} = \beta_k - H^{-1} g .$$

For a strictly convex objective with a Lipschitz Hessian, the method converges quadratically near the optimum: the number of correct digits roughly doubles each step. Far from the optimum it can overshoot, because the quadratic model is only local, and full Newton steps are then usually replaced by a damped version that halves the step until the objective decreases.

For the logistic model the gradient is $g = X^\top(\hat\pi - y)$ and the Hessian is $H = X^\top W X$ with $W = \operatorname{diag}(\hat\pi_i(1-\hat\pi_i))$. Substituting,

$$\beta_{k+1} = \beta_k + (X^\top W X)^{-1}X^\top(y - \hat\pi).$$

The term $y - \hat\pi$ is the residual on the probability scale, so the update is a weighted least-squares regression of that residual on $X$, using the current variances as weights. Rewriting the update as a weighted regression of a working response $z_i = x_i^\top\beta_k + (y_i - \hat\pi_i)/w_i$ on $X$ with weights $w_i$ gives the classical IRLS formulation: at each iteration, solve

$$\hat\beta = (X^\top W X)^{-1} X^\top W z .$$

Both forms are the same computation. The working-response form is the one that connects to a single weighted-least-squares routine, and it makes the dependence on the current fit visible: $z$ moves with $\beta$, and $W$ changes with $\hat\pi$, so the problem being solved changes every iteration.

import numpy as np

def irls(X, y, iters=25, tol=1e-8):
    beta = np.zeros(X.shape[1])
    for _ in range(iters):
        eta = X @ beta
        pi = 1.0 / (1.0 + np.exp(-eta))
        w = np.clip(pi * (1.0 - pi), 1e-10, None)   # weights vanish at saturation
        z = eta + (y - pi) / w                      # working response
        sw = np.sqrt(w)                             # absorb W into the design
        beta_new = np.linalg.lstsq(X * sw[:, None], z * sw, rcond=None)[0]
        if np.max(np.abs(beta_new - beta)) < tol:
            return beta_new
        beta = beta_new
    return beta

Scaling rows by $\sqrt{w}$ is numerically preferable to forming $X^\top W X$ explicitly, for the reason set out in ridge shrinkage in ill-conditioned designs: the cross-product squares the condition number of the design, and the weights can make that worse by spanning several orders of magnitude within one iteration. The square-root formulation keeps the problem in the design matrix’s own conditioning and lets a QR-based solver see it directly.

The weights are the information content

Each weight $w_i = \pi_i(1-\pi_i)$ is the Fisher information for observation $i$ about the linear predictor. It peaks at $\pi = 0.5$ with value $0.25$ and goes to zero as the probability saturates. The Hessian is therefore a sum of rank-one contributions $w_i x_i x_i^\top$, and an observation whose prediction is near 0 or 1 contributes almost nothing to the curvature of the objective.

This has a direct reading. A well-predicted observation carries little information about where the decision boundary is, because moving the boundary slightly would not change its prediction. An uncertain observation carries a lot, because it sits near the boundary and any movement of the boundary changes its fitted value. The IRLS iteration weights observations accordingly, and this is the mechanism behind the claim that logistic regression ignores examples it has already explained, not through any filtering step but because their Hessian contribution vanishes.

The same fact explains the algorithm’s behaviour on separable data. As the fit diverges, every $\pi_i$ approaches 0 or 1, so every weight approaches zero and $H$ becomes progressively closer to singular. The Newton step length, which involves $H^{-1}$, grows without bound, while the actual decrease in the objective shrinks. The iteration takes enormous steps in parameter space for negligible gains in loss, and eventually the linear solve loses all precision. That degenerate case is not repairable by more iterations; the unregularised maximum likelihood estimate does not exist, as discussed in logistic regression as conditional likelihood.

The practical fixes attach at different points. Adding a ridge penalty $\lambda I$ to $X^\top W X$ bounds the inverse below and makes every step finite, at the cost of a biased estimate. Clipping the weights to a floor, as in the code above, prevents division by zero in the working response but does not prevent divergence; it only makes the arithmetic finite. Detecting divergence and stopping is honest, provided the stopping criterion is a function of the data rather than an iteration count chosen in advance.

Damping and the trust region

A full Newton step is based on a quadratic approximation, and for a logistic model with a large coefficient the approximation can be poor. The standard safeguard is a line search along the Newton direction: accept the step only if it decreases the objective and, for a stronger guarantee, if it satisfies the Armijo condition

$$f(\beta + \alpha d) \le f(\beta) + c,\alpha g^\top d,$$

with $c$ a small number such as $10^{-4}$ and $\alpha$ the step length found by halving from one. The Newton direction satisfies $g^\top d < 0$ whenever $H$ is positive definite, so a decrease along the direction always exists, and backtracking terminates because the derivative at $\alpha = 0$ is strictly negative when the gradient is nonzero.

The alternative to a line search is a trust region, which bounds the step length and adapts the bound based on how well the quadratic model predicted the observed decrease. Trust regions are more robust when the Hessian is nearly singular, because the bound keeps the step finite even when $H^{-1}g$ is huge. For logistic regression with a shrinking Hessian, this is exactly the situation the method has to survive, and a trust-region Newton method is the principled version of what weight clipping approximates.

Convergence is quadratic in the final phase and roughly linear in the early phase. Each iteration costs one Hessian assembly of $O(np^2)$ work and one solve of $O(p^3)$, so with moderate $p$ and large $n$ the cost per iteration is dominated by forming $X^\top W X$, and the total cost is a small multiple of a single least-squares solve. For $p$ in the thousands and $n$ in the millions, the Hessian assembly becomes the bottleneck and methods that never form it, stochastic gradient variants, or quasi-Newton updates that build an approximate inverse from gradient differences, take over despite their slower convergence per iteration.

What the standard errors cost

The Hessian at convergence is the observed information matrix, and its inverse is the asymptotic covariance of the maximum likelihood estimate. When the model is correctly specified, the observed and expected information coincide, so $\operatorname{Cov}(\hat\beta) \approx (X^\top W X)^{-1}$ with $W$ evaluated at the fitted probabilities. The square roots of the diagonal are the standard errors, and the whole machinery is already computed as a byproduct of the last Newton step.

When the model is misspecified by a missing interaction, a wrong link, or clustered observations, the inverse information is not a consistent estimator of the sampling variance. The sandwich estimator replaces it with

$$\hat V = (X^\top W X)^{-1} \left(\sum_i \hat r_i^2 x_i x_i^\top\right) (X^\top W X)^{-1},$$

where $\hat r_i$ is the score contribution of observation $i$. For clustered data the middle term sums scores within clusters before squaring, and the difference between the two estimates is often a factor of two or more on real clustered survey or panel data.

import numpy as np

def robust_cov(X, y, beta):
    pi = 1.0 / (1.0 + np.exp(-(X @ beta)))
    W = pi * (1.0 - pi)
    A = X.T @ (X * W[:, None])            # observed information
    A_inv = np.linalg.inv(A)
    r = (y - pi)[:, None] * X             # per-observation score
    B = r.T @ r                           # no clustering here
    return A_inv @ B @ A_inv

Treating the converged weights as fixed constants is the approximation behind every “solve once and read the diagonal” approach. It is accurate to the extent that the fitted probabilities are close to their true values, which is generally the case at convergence, but it degrades when many observations have extreme fitted values, because a handful of weights then dominate the information and the effective sample size is much smaller than $n$.

Iteration count is also a diagnostic. A fit that converges in five iterations on data with a single informative direction and twenty on data that is nearly separable has a different explanation for its coefficient values, and logging the iteration count alongside the estimate makes that difference visible.

The weights also define a leverage for each observation in the weighted regression, $h_i = w_i x_i^\top (X^\top W X)^{-1} x_i$, and the sum of the leverages equals the number of parameters. Since $w_i$ is at most $0.25$, an observation can have high leverage only by sitting far out in the design, meaning a feature vector that few other rows resemble, rather than by being uncertain. The two sources of influence are therefore separable here: uncertainty enters through the weights, and unusual position enters through the design, and an observation that combines both dominates the fit while an observation with only one of them has the effect a linear model would give it. That separation has no counterpart in least squares, where the hat matrix diagonal sums to the parameter count on its own and leverage alone determines influence.