Floating-Point Error in GLM Computation

Why Numerical Precision Matters

Likelihood calculations can be mathematically well defined and still be unreliable when carried out on a computer. A likelihood may be a product of many small probabilities; a log-likelihood may be a sum of many terms; and an optimizer may compare two log-likelihood values that are nearly equal. Finite-precision arithmetic can underflow, round intermediate results, or erase a small difference against a much larger number.

These effects matter in generalized linear model (GLM) computation because optimization decisions depend on numerical quantities: whether a proposed step increases the log-likelihood, whether an update is sufficiently small, or whether two parameter values fit almost equally well. A computed value is not automatically an accurate value. The first step is to understand what the computer can represent and how arithmetic changes the result.

Floating-Point Numbers and Rounding

A floating-point system represents real numbers using a finite number of binary digits. In the widely used IEEE 754 binary64 format, numbers near 1 have 53 significant binary digits, including the leading digit. Arithmetic results are rounded to a representable number; the usual default is rounding to nearest, with ties rounded to an even final significand digit.

For a basic operation whose exact result is in the normal range and does not overflow, a common rounding model is \(\operatorname{fl}(a\mathbin{\circ}b)=(a\mathbin{\circ}b)(1+\delta)\), where \(\circ\) denotes addition, subtraction, multiplication, or division and \(|\delta|\leq u\). This is a relative-error description of one operation, not a guarantee that an entire long calculation has relative error at most \(u\). Errors can accumulate, and the model does not apply unchanged to results that are subnormal or underflow to zero.

Binary64 also has a finite range. Its smallest positive normal number is \(2^{-1022}\), and its smallest positive subnormal number is \(2^{-1074}\). Subnormal numbers extend the range toward zero, but their spacing is fixed rather than proportional to their size. A calculation that continues below this range may round to zero. At the other end, a result larger than the finite range overflows.

Products, Likelihoods, and Accumulated Rounding

For independent observations, a likelihood is often expressed as a product of individual probability or density contributions. Each contribution may be representable, yet their product can be too small to represent. In exact arithmetic, taking logarithms turns a product into a sum: \(\log\{\prod_{i=1}^{n} f_i\}=\sum_{i=1}^{n}\log(f_i)\), when all \(f_i>0\). This identity explains why log-likelihoods are widely used, but it does not eliminate rounding: each logarithm and each addition is still computed with finite precision.

A simple bound describes how rounding can accumulate in a product. Assume each multiplication is rounded according to the relative-error model, and assume no intermediate result overflows or underflows. If \(p\) is the exact product of \(n\) factors and \(\widehat p\) is the computed product, then

Proof. Write each rounded multiplication as the exact multiplication followed by a factor \(1+\delta_j\), with \(|\delta_j|\leq u\). Repeated substitution gives \[ \widehat p=p\prod_{j=1}^{n}(1+\delta_j). \] Consequently, \[ \frac{|\widehat p-p|}{|p|} \leq \prod_{j=1}^{n}(1+|\delta_j|)-1 \leq (1+u)^n-1. \] By the binomial theorem, \((1+u)^n-1=\sum_{k=1}^{n}{n\choose k}u^k\). Since \({n\choose k}\leq n^k\), this is at most \(\sum_{k=1}^{n}(nu)^k\). For \(nu<1\), the geometric-series bound gives \[ \sum_{k=1}^{n}(nu)^k < \frac{nu}{1-nu}=\gamma_n. \] This proves the stated bound. \(\square\)

This result is a bound under explicit assumptions, not a prediction that every product incurs the maximum error. It also cannot rescue a calculation that underflows: the relative-error model used in the proof excludes that case.

Cancellation When Comparing Log-Likelihoods

Cancellation occurs when nearly equal numbers are subtracted. Even if the subtraction of the represented numbers is carried out accurately, rounding may already have discarded the small difference in the inputs. This is especially relevant when comparing log-likelihoods that share a large common component.

Proof. Subtract the exact difference from the perturbed difference: \[ \big[(x+\Delta x)-(y+\Delta y)\big]-(x-y)=\Delta x-\Delta y. \] The triangle inequality gives \[ |\Delta x-\Delta y|\leq |\Delta x|+|\Delta y| \leq \varepsilon(|x|+|y|). \] Dividing by \(|x-y|\), which is positive because \(x\neq y\), gives the relative-error bound. When \(|x-y|\) is small compared with \(|x|\) and \(|y|\), this bound can be large. \(\square\)

In some computations, an algebraic expression for a difference can avoid carrying large shared terms through two separate evaluations. Whether that is possible depends on the likelihood and the calculation being performed. The general diagnostic is to compare the size of the difference with the scale and precision of the values from which it was computed.

Reading Numerical Results Carefully

A computed log-likelihood should be interpreted with awareness of the operations that produced it. If a calculation sums many terms, rounding may accumulate; if it multiplies many contributions, underflow may erase the product; and if it subtracts nearby values, cancellation may make the difference inaccurate. These are distinct mechanisms, so a single statement such as “double precision is accurate to 16 digits” is not a sufficient diagnosis.

In optimization, numerical effects can interact with convergence criteria. For example, a reported likelihood change of zero could mean that the objective truly did not change, or that the change was smaller than the resolution of the computed values. Likewise, an apparent improvement smaller than the rounding uncertainty should not be treated as decisive evidence that one parameter value is better. Interpret the objective alongside the size of the step, the scale of the values, and the behavior of the score and information calculations.

Common pitfall: “The computer printed the same log-likelihood for both parameter values, so the two likelihoods must be equal.” This conclusion confuses equality of rounded outputs with equality of exact quantities. If the values are large and their difference is smaller than their local floating-point spacing, both may round to the same number. The correction is to treat the printed equality as a statement about numerical resolution, not a proof of mathematical equality. In particular, do not interpret more displayed decimal places as evidence of more underlying precision.