Reliable Stopping Rules for Iterative GLM Fits

Why a GLM Fit Needs More Than One Stopping Check

An iterative fitting algorithm must decide when to stop updating the coefficient vector. A small update is not, by itself, evidence that the likelihood has been optimized: the algorithm may have taken a very cautious step. A small change in log-likelihood is not enough either, because the objective can be locally flat while the score remains appreciable. Conversely, a score that looks small in the raw coefficient units can be misleading when predictors are measured on very different scales.

A reliable stopping decision checks three different features of the fit: whether the coefficient update is small relative to meaningful coefficient scales, whether the score is small relative to local curvature, and whether the objective has stopped changing on a per-observation scale. These are complementary diagnostics, not three interchangeable ways to measure the same thing.

Define the Three Criteria

Consider independent observations with a differentiable GLM log-likelihood \(\ell(\beta)\), where \(\beta\) is a \(p\)-dimensional coefficient vector. At iteration \(k\), let \(\beta^{(k)}\) be the current accepted iterate and \(\beta^{(k+1)}\) the proposed accepted update. Write \(U(\beta)=\partial\ell(\beta)/\partial\beta\) for the score and \(J(\beta)=-\partial^2\ell(\beta)/\partial\beta\partial\beta^\mathsf{T}\) for the observed information when these derivatives exist.

For each coefficient, choose a positive reference scale \(s_j\), expressed in the units of \(\beta_j\). It should reflect a meaningful change in that coefficient, not merely a convenient computer representation. A mixed absolute-relative scaled step is

R_{\mathrm{step}}= \max_{1\leq j\leq p} \frac{\left|\beta_j^{(k+1)}-\beta_j^{(k)}\right|} {s_j+\max\!\left(\left|\beta_j^{(k)}\right|,\left|\beta_j^{(k+1)}\right|\right)}.

The scale \(s_j\) prevents division by a value near zero; the coefficient-dependent term also makes the check relative for large coefficients. The step check passes when \(R_{\mathrm{step}}\leq\varepsilon_{\mathrm{step}}\), for a chosen dimensionless tolerance \(\varepsilon_{\mathrm{step}}\).

To assess the score, use a positive-definite local curvature matrix \(M\) in coefficient coordinates, such as the information matrix used by the fitting algorithm when it is positive definite. Define

R_{\mathrm{score}}= \sqrt{U(\beta^{(k+1)})^\mathsf{T}M^{-1}U(\beta^{(k+1)})}.

This score measure accounts for the scale of local curvature. In computation, evaluate it by solving a linear system rather than forming a matrix inverse. If the available curvature matrix is singular or seriously ill-conditioned, this diagnostic is not trustworthy; that is a reason to investigate the fit, not a reason to declare convergence.

For \(n\) independent observations, the objective check uses the average log-likelihood change:

R_{\mathrm{obj}}= \frac{\left|\ell(\beta^{(k+1)})-\ell(\beta^{(k)})\right|}{n}.

The absolute value detects a small change in either direction; it does not assert that a decrease is acceptable. A maximization algorithm should still reject or investigate updates that decrease the objective contrary to its update rule. The division by \(n\) makes the tolerance refer to average change per observation. Unlike a relative change divided by the absolute value of the log-likelihood, this check is unchanged if a constant is added to the log-likelihood.

Stopping rule: Declare numerical convergence only when the scaled step, the information-scaled score, and the per-observation objective change all pass their specified tolerances, the objective and derivatives are finite, and the curvature calculation is usable. Passing these checks is a numerical stopping decision; it is not, by itself, proof that a global maximum has been found.

What the Score Check Measures

When \(M=J(\beta)\) is positive definite, the score criterion is the length of the score measured in the local information metric. It is invariant to a nonsingular linear change of coefficient coordinates. To see this, let \(\beta=A\gamma\) for an invertible matrix \(A\). By the chain rule, the score and information in \(\gamma\) coordinates are \(U_\gamma=A^\mathsf{T}U_\beta\) and \(J_\gamma=A^\mathsf{T}J_\beta A\). Therefore

U_\gamma^\mathsf{T}J_\gamma^{-1}U_\gamma = U_\beta^\mathsf{T}A \left(A^\mathsf{T}J_\beta A\right)^{-1} A^\mathsf{T}U_\beta = U_\beta^\mathsf{T}J_\beta^{-1}U_\beta.

The equality follows from \((A^\mathsf{T}J_\beta A)^{-1}=A^{-1}J_\beta^{-1}(A^\mathsf{T})^{-1}\). Thus, changing coefficient units does not change this score measure. This invariance does not remove the need for a usable positive-definite curvature matrix: if the matrix cannot be solved against reliably, the score check cannot be relied on.

A Local Result Connecting Score and Objective Change

The score criterion has a direct interpretation for a full Newton step. The following result makes that interpretation precise under an explicit local curvature condition.

Theorem: Suppose \(\ell\) is twice continuously differentiable on a neighborhood of \(\beta\), \(J(\beta)\) is positive definite, and the full Newton step for maximizing \(\ell\) is \(\delta=J(\beta)^{-1}U(\beta)\). Suppose, for every \(t\in[0,1]\), the information along the step satisfies \[ (1-\rho)J(\beta)\preceq J(\beta+t\delta)\preceq(1+\rho)J(\beta) \] for some \(0\leq\rho<1\). Then, with \(\lambda^2=U(\beta)^\mathsf{T}J(\beta)^{-1}U(\beta)\), \[ \frac{1-\rho}{2}\lambda^2 \leq \ell(\beta+\delta)-\ell(\beta) \leq \frac{1+\rho}{2}\lambda^2. \]

Proof. Define \(g(t)=\ell(\beta+t\delta)\) for \(0\leq t\leq1\). The chain rule gives \(g'(0)=U(\beta)^\mathsf{T}\delta=\lambda^2\) and \(g''(t)=-\delta^\mathsf{T}J(\beta+t\delta)\delta\). Twice integrating the derivative yields the exact identity

\ell(\beta+\delta)-\ell(\beta) = \lambda^2-\int_0^1(1-t)\, \delta^\mathsf{T}J(\beta+t\delta)\delta\,dt.

The assumed matrix inequalities imply that the integrand's quadratic form lies between \((1-\rho)\delta^\mathsf{T}J(\beta)\delta\) and \((1+\rho)\delta^\mathsf{T}J(\beta)\delta\). Since \(\delta=J(\beta)^{-1}U(\beta)\), we have \(\delta^\mathsf{T}J(\beta)\delta=\lambda^2\). Also, \(\int_0^1(1-t)\,dt=1/2\). Substitution gives the stated lower and upper bounds. \(\square\)

This result explains why a small information-scaled score is informative about a Newton step's local predicted improvement when the curvature is positive definite and does not vary too much along that step. It does not say that every algorithm takes a full Newton step: damping, step control, or other update choices can make the actual coefficient change much smaller.

Worked Examples

Worked Example: All Three Checks Pass Near a Logistic Fit

In an intercept-only logistic model with \(n=100\) Bernoulli observations, suppose 50 responses are successes. Omitting the constant binomial coefficient, the log-likelihood is \(\ell(\beta)=50\beta-100\log(1+\exp(\beta))\). Its score is \(U(\beta)=50-100p(\beta)\), and its observed information is \(J(\beta)=100p(\beta)(1-p(\beta))\), where \(p(\beta)=1/(1+\exp(-\beta))\).

Suppose an accepted update goes from \(\beta^{(k)}=0.001\) to \(\beta^{(k+1)}=0.0001\), with reference scale \(s_1=1\). Then

R_{\mathrm{step}}= \frac{0.0009}{1+0.001} \approx 0.0008991.

At the new value, \(p(0.0001)\approx0.500025\), so \(U(0.0001)\approx-0.0025\) and \(J(0.0001)\approx25\). The score measure is approximately \(|-0.0025|/\sqrt{25}=0.0005\).

Because the data are balanced, \(\ell(\beta)-\ell(0)=-100\log\cosh(\beta/2)\). Thus the objective increase from \(0.001\) to \(0.0001\) is \[ 100\{\log\cosh(0.0005)-\log\cosh(0.00005)\} \approx 0.000012375. \] The per-observation change is about \(1.2375\times10^{-7}\). With tolerances \(0.001\) for the step, \(0.01\) for the score, and \(10^{-6}\) for the per-observation objective change, all three checks pass.

Worked Example: Rescaling a Predictor Changes Raw Steps

Suppose a linear predictor contains \(\beta x\), and a predictor is re-expressed as \(z=100x\). The equivalent coefficient is \(\gamma=\beta/100\), since \(\beta x=\gamma z\). A particular update changes \(\beta\) from \(0.50\) to \(0.49\), and therefore changes \(\gamma\) from \(0.005\) to \(0.0049\). Both updates change the linear predictor by the same amount: at \(x=2\), the change is \(0.01(2)=0.02\); at \(z=200\), it is \(0.0001(200)=0.02\).

A raw absolute-step rule with tolerance \(0.001\) says the \(\gamma\) update, of size \(0.0001\), passes, but the \(\beta\) update, of size \(0.01\), fails. That disagreement is due only to the units. Instead, choose \(s_\beta=1\) and \(s_\gamma=0.01\), which transform in the same ratio as the coefficients. The scaled steps are \[ \frac{0.01}{1+0.50}\approx0.00667 \quad\text{and}\quad \frac{0.0001}{0.01+0.005}\approx0.00667. \] The scaled criterion gives the same assessment in both parameterizations.

Worked Example: A Tiny Objective Change Can Hide a Nontrivial Score

Consider one Gaussian observation \(Y=0.1\) from a normal model with known variance \(1\) and fitted mean \(\beta\). Ignoring an additive constant, \(\ell(\beta)=-\tfrac12(0.1-\beta)^2\), \(U(\beta)=0.1-\beta\), and \(J(\beta)=1\). Suppose an algorithm moves from \(\beta=0\) to \(\beta=0.001\), with \(s_1=1\).

The scaled step is \(0.001/(1+0.001)\approx0.000999\). The log-likelihood change is \[ -\tfrac12(0.099)^2+\tfrac12(0.1)^2 =0.0000995, \] so the per-observation change is also \(0.0000995\). Both pass tolerances of \(0.002\) for the step and \(0.001\) for objective change. But at the new iterate the score measure is \(\sqrt{0.099^2/1}=0.099\), which fails a score tolerance of \(0.01\). The update and objective change are small because the step is tiny; the fit is not yet stationary according to the score check.

Practical Use and Common Pitfalls

Evaluate all three checks at the accepted iterate, after verifying that the objective and the quantities needed for the score are finite. Choose tolerances before interpreting a fit, and choose coefficient reference scales to match the coefficient units. If the model has predictors on very different scales, raw coefficient changes can be especially hard to compare; predictor scaling and a matching transformation of reference scales make the step criterion more interpretable.

A common mistake is to stop as soon as the objective changes by less than a relative tolerance such as \(|\Delta\ell|/|\ell|\). This can be misleading because adding a constant to the log-likelihood changes the denominator without changing the fit or \(\Delta\ell\). It can also allow a heavily damped, tiny update to pass while the score remains large. Use a declared objective scale such as per-observation change, and require the score and step checks as well.

Another mistake is to treat three passing tolerances as proof of a maximum. They show that the algorithm has made little progress, the local score is small under the chosen curvature scale, and the objective is changing little at the stated numerical resolution. They do not establish that the likelihood has a global maximum, or that the numerical derivatives and objective evaluations are accurate enough. A poorly conditioned or unusable curvature calculation, non-finite values, or an objective decrease that conflicts with the algorithm's acceptance rule calls for investigation rather than an automatic convergence label.

Key takeaway: A dependable GLM stopping rule combines a coefficient-scale-aware step, a curvature-scaled score, and an objective change on a declared per-observation scale. A single small raw number is not a convergence diagnosis.

Check Your Understanding

Use the definitions and examples above to assess each situation.

Why can the same predictor effect produce different raw coefficient-step sizes after changing the predictor's units?
What does the information-scaled score measure, and what condition on its curvature matrix is required?
Why is a per-observation objective change unaffected by adding a constant to the log-likelihood?
An algorithm reports a very small objective change but a score measure above tolerance. Should it stop under the combined rule? Explain.
Do three passing stopping checks prove that the fitted coefficients maximize the likelihood globally? Why or why not?