Common Numerical Stability Mistakes in GLMs

Why Numerical Details Change a GLM Fit

A generalized linear model is defined by its probability model and its likelihood, not by the particular floating-point calculations used to fit it. A computation that looks harmless—such as replacing probabilities near zero by a small positive constant—can change the objective being optimized. Other shortcuts can instead produce overflow or an invalid value even when the exact likelihood is finite. These failures matter because an optimizer can only respond to the numbers it is given.

This tutorial focuses on four diagnostic questions: Has a probability been clipped in a way that changes the likelihood? Has an exponential been calculated before a logarithm that could have been evaluated stably? Does the likelihood have a meaningful limit at the boundary? And does a small coefficient update actually indicate that the fit is reliable?

Clipping Is Not an Exact Likelihood Calculation

For a Bernoulli observation \(y\in\{0,1\}\) with success probability \(p\), the log-likelihood contribution is \(y\log p+(1-y)\log(1-p)\). In a logistic model, \(p=1/(1+\exp(-\eta))\), where \(\eta\) is the linear predictor. Clipping replaces \(p\) in a calculation with a constrained value. For the symmetric clipping rule below, the restriction on \(\varepsilon\) is essential.

Definition: For \(0<\varepsilon<1/2\), probability clipping by \(\varepsilon\) replaces \(p\) by \(\widetilde p=\min(1-\varepsilon,\max(\varepsilon,p))\). Thus \(\widetilde p\in[\varepsilon,1-\varepsilon]\), a nonempty interval contained in \((0,1)\). Clipping may be useful for a display or a deliberately modified procedure, but it does not in general preserve the original likelihood.

The effect is not merely that the reported log-likelihood is slightly different. In a clipped region, the modified contribution can be constant as the model parameter changes, so its derivative no longer gives the score of the stated model.

Worked Example: Clipping Erases a Bernoulli Score

Consider one Bernoulli observation \(y=1\), a logistic predictor \(\eta=-40\), and clipping with \(\varepsilon=10^{-6}\). The model probability is

p=\frac{1}{1+\exp(40)}\approx 4.248354255\times 10^{-18}.

Without clipping, the log-likelihood is \(\log p=-\log(1+\exp(40))\), which is approximately \(-40\). Clipping replaces \(p\) with \(10^{-6}\), giving \(\log(10^{-6})\approx-13.81551056\). The clipped contribution is therefore about \(26.18448944\) larger than the exact contribution.

The score with respect to \(\eta\) for this observation is \(1-p\), approximately \(1\). But while \(p<\varepsilon\), the clipped probability is fixed at \(\varepsilon\), so the derivative of the clipped log contribution with respect to \(\eta\) is zero. An optimizer using that clipped objective receives no local signal to increase \(\eta\), even though the original model's score is nearly \(1\). The calculation has changed both the objective and its derivative.

Clipping can still be an explicit part of a chosen approximation or reporting convention; the important point is to distinguish that modified quantity from the likelihood of the original GLM. In particular, clipping should not silently replace a stable likelihood calculation.

Take Logarithms Without First Creating an Overflow

A standard source of overflow is the logistic log-likelihood written using \(\log(1+\exp(\eta))\). The exact expression is well defined for every finite \(\eta\), but evaluating \(\exp(\eta)\) directly can overflow when \(\eta\) is large. The function \(s(\eta)=\log(1+\exp(\eta))\), often called the softplus function, has an algebraically equivalent form that avoids this problem.

Result: For every real \(\eta\),

s(\eta)= \begin{cases} \eta+\log(1+\exp(-\eta)),&\eta\geq 0,\\ \log(1+\exp(\eta)),&\eta<0. \end{cases}

Derivation. If \(\eta\geq0\), factor \(\exp(\eta)\) from \(1+\exp(\eta)\):

\log(1+\exp(\eta)) =\log\!\left(\exp(\eta)(\exp(-\eta)+1)\right) =\eta+\log(1+\exp(-\eta)).

If \(\eta<0\), the original expression already has an exponential with argument less than zero. In the first branch, the exponential argument is \(-\eta\leq0\). Thus neither branch requires exponentiating a positive argument larger than one, which avoids the large-exponential overflow that the direct formula can cause. The branches agree at \(\eta=0\), where both give \(\log 2\). \(\square\)

For a Bernoulli logistic observation, the log-likelihood can be written \(y\eta-s(\eta)\). A further identity avoids subtracting two large, nearly equal quantities when \(y=1\): the contribution is \(-s(-\eta)\). When \(y=0\), it is \(-s(\eta)\). These forms also make the contribution’s sign and limiting behavior easier to inspect.

Worked Example: A Finite Logistic Log-Likelihood at an Extreme Predictor

Suppose \(y=0\) and \(\eta=1000\). The exact contribution is \(-\log(1+\exp(1000))\). Factoring as above gives

-\log(1+\exp(1000)) =-1000-\log(1+\exp(-1000)) \approx -1000.

The approximation is valid because \(\exp(-1000)\) is negligible at ordinary floating-point precision. By contrast, a direct calculation of \(\exp(1000)\) overflows in common floating-point arithmetic. The subsequent expression \(-\log(1+\exp(1000))\) can then become negative infinity, although the exact log-likelihood is finite and close to \(-1000\).

For \(y=1\) at the same predictor, use the equivalent contribution \(-s(-1000)=-\log(1+\exp(-1000))\), which is close to zero. Directly forming \(1000-\log(1+\exp(1000))\) can instead involve an infinite intermediate value and return negative infinity. The stable form avoids that failure. A floating-point result of zero for this last contribution is an approximation at machine precision, not a statement that the exact contribution is exactly zero.

Likelihood Limits at the Boundary

A probability parameter can approach zero or one even when the model’s finite parameter values yield probabilities strictly between them. The likelihood’s limit depends on the observed response. It is not correct to treat every logarithm of a boundary probability as an undifferentiated computational error.

Result: For a binomial observation \(Y=y\) from \(m\) trials, with integer \(0\leq y\leq m\), the log-likelihood as a function of success probability \(p\in(0,1)\) is

\ell(p)=\log {m\choose y}+y\log p+(m-y)\log(1-p).

At \(p\downarrow0\), its limit is finite if \(y=0\), and equals \(\log {m\choose 0}=0\); if \(y>0\), its limit is \(-\infty\). At \(p\uparrow1\), its limit is finite if \(y=m\), and equals \(\log {m\choose m}=0\); if \(y<m\), its limit is \(-\infty\).

Proof. As \(p\downarrow0\), \(\log(1-p)\to0\). If \(y=0\), the term \(y\log p\) is identically zero for \(p>0\), and \((m-y)\log(1-p)\to0\). If \(y>0\), then \(y\log p\to-\infty\), while the other terms have finite limits, so \(\ell(p)\to-\infty\). At the other endpoint, let \(q=1-p\). As \(q\downarrow0\), \(y\log p\to0\); if \(y=m\), the remaining logarithmic term has coefficient zero, while if \(y<m\), \((m-y)\log q\to-\infty\). This proves both limits. \(\square\)

The convention \(0\log 0=0\) is a concise way to represent the finite endpoint limits in expressions such as this one; it does not make \(\log 0\) finite. A positive coefficient multiplying \(\log 0\) still corresponds to a log-likelihood of negative infinity.

Worked Example: Distinguishing Two Poisson Boundary Limits

For a Poisson response \(Y=y\) with mean \(\mu>0\), the log-likelihood contribution is \(y\log\mu-\mu-\log(y!)\). If \(y=0\), this becomes \(-\mu\), so as \(\mu\downarrow0\) it tends to \(0\). At \(\mu=10^{-300}\), the contribution is \(-10^{-300}\), a finite value extremely close to zero.

If instead \(y=1\), the contribution is \(\log\mu-\mu\). At \(\mu=10^{-300}\), it is approximately \(-690.7755\), since \(\log(10^{-300})=-300\log 10\approx-690.7755\). As \(\mu\downarrow0\), this contribution tends to \(-\infty\), not to zero. Thus a boundary limit must be checked against the observed response before deciding whether a numerical value represents an error or a genuine likelihood limit.

For a GLM whose parameterization only permits interior means at finite coefficients, a boundary value may be approached without being attained by any finite coefficient vector. That distinction matters when interpreting an optimizer: a large coefficient or a boundary-like fitted mean does not, by itself, establish that a finite maximum has been found.

A Small Coefficient Step Is Not a Convergence Certificate

An iterative method can produce a small update for several reasons, including high curvature, rescaling, or a restricted step. The size of the update alone does not establish that the score is negligible, that the objective has been evaluated correctly, or that the fit satisfies the intended stopping criteria. The following binomial calculation isolates the high-curvature case.

Worked Example: A Tiny Newton Update with a Nonzero Score

Consider one binomial observation with \(m=10^{12}\) trials and \(y=500{,}000{,}000{,}001\) successes, modeled with logistic predictor \(\beta\). At \(\beta=0\), \(p=1/2\). The score and information for this observation are

U(\beta)=y-mp=1,\qquad I(\beta)=mp(1-p)=\frac{10^{12}}{4}=2.5\times10^{11}.

A Newton step for maximizing the log-likelihood is \(U/I\), so its magnitude is

\left|\frac{U}{I}\right| =\frac{1}{2.5\times10^{11}} =4\times10^{-12}.

That is an extremely small coefficient update, but the score at the current value is still \(1\), not zero. The example does not imply that a score of one is always practically important: its interpretation depends on scale and on the fitting procedure’s tolerances. It does show why a step-size test by itself cannot certify stationarity. A defensible diagnostic also checks an appropriately scaled score, objective behavior, parameter and predictor finiteness, and any constraints relevant to the fit.

Practical Checks and Common Mistakes

When a GLM fit behaves unexpectedly, trace the numerical calculation back to the mathematical quantity it is meant to represent. A stable evaluation should preserve the model’s likelihood and score, handle legitimate boundary limits explicitly, and distinguish a computational stopping signal from evidence that the fit is satisfactory.

Mistake: clipping and calling the result the original likelihood. Clipping can change objective values and make derivatives zero in regions where the original score is not zero. Use stable formulas for the original likelihood; if clipping defines a deliberate approximation, label and assess that approximation as such.
Mistake: calculating \(\exp(\eta)\) before taking a logarithm. An intermediate exponential may overflow even when the exact log-likelihood is finite. Use the branchwise softplus identity, or a trusted implementation of a stable log-one-plus-exponential function.
Mistake: treating every boundary value as the same kind of failure. For some responses, the likelihood has a finite boundary limit; for others, it tends to zero on the likelihood scale, or negative infinity on the log-likelihood scale. Check the response-specific limit.
Mistake: accepting a fit because coefficients barely changed. A small update can coexist with a nonzero score, and a numerically small update does not validate the objective calculation. Check complementary diagnostics rather than relying on step size alone.

Key takeaway: Preserve the intended likelihood, evaluate it without unstable intermediate arithmetic, and interpret endpoint limits and stopping diagnostics in context. Numerical stability is not a cosmetic refinement: it determines whether the computed objective and its derivatives correspond to the GLM being fitted.

Check Your Understanding

For each question, explain what the calculation represents, not only which formula you would use.

Why must the symmetric clipping rule require \(0<\varepsilon<1/2\), and how can clipping alter the score as well as the log-likelihood?
For a logistic Bernoulli observation with \(y=0\) and \(\eta=900\), explain why direct exponentiation can fail and state a stable expression for the contribution.
What is the limiting binomial log-likelihood as \(p\downarrow0\) when \(y=0\)? What changes when \(y>0\)?
For a Poisson response \(y=0\), what is the log-likelihood limit as \(\mu\downarrow0\), and why is that different from the case \(y=1\)?
Give one reason a Newton coefficient step can be tiny even when the current score is nonzero. Name at least two other diagnostics that should be checked before accepting a fit.