Numerical accuracy and speed in orm’s Fortran likelihood code

This document catalogs the changes made to src/orm_links.f90, src/ormll.f90 and src/ormeta.f90 for computational accuracy and speed, explains the underlying problem — computing differences of cumulative probabilities accurately — in detail, and describes how the design connects to the future addition of y-dependent effects (YDE / partial proportional odds).

None of this changes the statistical model. orm.fit.r, orm.rfit.r, init.c’s function bodies, and every .Fortran() call signature except ormeta’s (see below) are unchanged.

1. The problem: differencing cumulative probabilities

orm’s log-likelihood is built from terms of the form

d = F(x1) - F(x2)            (an interior or interval-censored category)
d = F(x1)  or  1 - F(x1)     (a boundary or singly-censored category)

where x1 = alpha(ia) + lp and x2 = alpha(ia2) + lp, F is the CDF for one of five link functions (logistic, probit, log-log, complementary log-log, Cauchy), and alpha holds the fitted intercepts.

With continuous or near-continuous Y, adjacent intercepts are close together — for N = 250,000, typical adjacent gaps are on the order of 1e-5 — so d is the difference of two numbers that agree to 4–5 significant digits. Two independent things can destroy that difference’s accuracy:

1. Cancellation in the subtraction itself. F(x1) and F(x2) each carry an absolute rounding error of about 1e-16 (or, in the upper tail, an error of order 1e-16·F, which stops shrinking once F saturates at 1 − 1e-16 in double precision). That fixed-size error does not shrink just because d is small, so the relative error of d grows in proportion to how narrow the cell is.
2. Loss of the cell width before the subtraction ever happens. x1 and x2 are each rounded to a double (alpha(ia) + lp, alpha(ia2) + lp) before F is ever evaluated. If lp is large (e.g. a strong linear predictor, or, in the future, a y-dependent shift), the rounding of x1 and x2 individually can already be larger than the true gap between them — at which point no amount of care in evaluating F can recover the correct d, because the input arguments have already lost the information.

Both problems compound for the tail-saturating links (log-log, complementary log-log): with the naive 1 - exp(-exp(x)) and exp(-exp(-x)) forms, F and 1 - F reach their float limits (0 or 1) at much smaller |x| than the logistic, so d collapses to exactly 0 — not just inaccurately small — well within the range of x actually encountered in fits, triggering spurious salloc = 999 failures.

Getting d right the first time is what matters. Once d has full relative accuracy, log(d) inherits it to within one rounding (about 1e-16), and a simple summation of wt * log(d) — subject to the ordinary accumulation error described in §2.5 — gives an accurate log-likelihood. No other transformation downstream of d can repair an inaccurate d.

1.1 Measured impact of the un-fixed code

Before this work, ormll.f90 computed d as p1 - p2 from two ordinary expit/pnorm/etc. evaluations. Measured against an arbitrary-precision (mpmath) reference at cell width 1e-5:

“Wrong” means d was off by order 1 relative error, or exactly 0 (triggering a spurious convergence failure), well inside the range of linear-predictor values a real fit can reach.

2. Changes made, in order

2.1 Accurate one-sided probabilities: never compute a small 1 - F (or F) by subtraction

linkpieces(x, link, need, f, g, p, t) returns f = F(x), g = 1 - F(x), and the pdf p, each computed so that the small one is never obtained by subtracting from 1:

• Logistic, probit, Cauchy: 1 - F(x) = F(-x), so the small tail is evaluated directly at the reflected argument (erfc(-x/√2) for probit, atan(1/x) for Cauchy) rather than as 1 - F(x).
• Log-log ↔︎ complementary log-log: the two link functions are reflections of each other (F_loglog(x) = exp(-exp(-x)), F_cloglog(x) = 1 - exp(-exp(x))), so each one’s complement is computed as the other link’s small-side value, again via expm1 rather than subtraction from 1.

t = exp(∓x) is also returned by linkpieces, since linkdpdf and celldiff need it and it should not be recomputed (§2.6).

2.2 Accurate two-point differences: celldiff uses the exact cell width, not two rounded arguments

celldiff(link, dl, x1, x2, f1, g1, f2, g2, p1, p2, t1, t2) takes the exact cell width dl = (alpha(ia) - alpha(ia2)) + (e1 - e2) (formed from small terms — see §2.4 and §3) as a first-class input, separate from x1 and x2 themselves, and uses a closed-form difference wherever one exists:

Every term in these formulas is non-negative, so there is no cancellation — this is the key structural idea, and it is what makes d accurate to full double precision (about 1e-15–1e-16 relative) however narrow the cell, for four of the five links.

For dl > 30 (only reachable via a wild Newton iterate, since a converged cell this wide has a probability below 1e-13), sinh/expm1 risk overflow; celldiff instead falls back to differencing the two accurate small-side tails (§2.1): G(x2) - G(x1) when x2 > 0, else F(x1) - F(x2). Both terms are then already tiny and well separated, so ordinary subtraction is safe there.

2.3 Probit: no closed form, so tail-switching is the best available

Probit has no algebraic identity for Φ(x1) - Φ(x2). celldiff differences whichever pair of one-sided values (§2.1) is smaller in magnitude (G(x2) - G(x1) when x2 > 0, else F(x1) - F(x2)), which removes the tail-saturation failure (probit no longer returns 999 or a wildly wrong d for moderately large |x|) but leaves a residual cancellation error of order eps / width — measured at about 4e-10 relative for a 1e-5-wide cell, versus complete failure in the original code at the same width once |x| exceeded about 8. This is the one link where narrow-cell accuracy is still limited by cancellation; there is no way to remove that within double precision without a different algorithm (e.g. an extended-precision or series-based Φ(x1) - Φ(x2) evaluation), which was judged not worth the complexity for this pass.

2.4 ormatoms: one link-agnostic per-observation routine

All link-specific arithmetic in the package now lives in orm_links.f90. ormatoms(link, n, k, alpha, ia, ia2, sgn, lp, wt, what, d, s1, s2, hee, h11, h22, h12, h1e, h2e, salloc, e1, e2) is the single per-observation loop that both ormll (the main Newton fit) and ormeta (per-cluster/per-node evaluation for random-intercept and future YDE fits) call; neither contains any link arithmetic of its own — they only accumulate what ormatoms returns into u/ha/hb/hab (ormll) or logd/g/h/s1/s2 (ormeta).

what controls how much is computed per observation:

• what = 1: d only (log-likelihood evaluation). Uses light routines (oned, celld) that call linkpieces only for the pieces (F, G, or the pdf) actually needed by that link’s formula in §2.2 — e.g. Cauchy needs none of them (the atan2 form uses x1, x2 directly), log-log needs only F(x1), cloglog only G(x2).
• what = 2: also s1, s2 (the score contributions to alpha(ia) and alpha(ia2) separately).
• what = 3: also hee, and, when the caller supplies storage for them, h11, h22, h12, h1e, h2e (the full second-derivative set).

For the logistic link, the second derivatives have closed forms with no division-heavy cancellation, verified against 450-digit arithmetic:

one intercept:   hee = h11 = h1e = -wt * pdf1
two intercepts:  hee = -wt*(pdf1+pdf2)
                 h1e = -wt*pdf1,  h2e = -wt*pdf2
                 h12 =  wt*pdf1*pdf2 / d^2
                 h11 = h1e - h12,  h22 = h2e - h12

using the identity d²(D/Dη) = D/Da1 + D/Da2, so h1e = h11+h12, h2e = h12+h22, hee = h11+2·h12+h22 hold exactly. Links 2–5 use the direct pdf/pdf-derivative formulas (still exact to the pdf’s own accuracy, just without the special-case cancellation-free logistic identities).

2.5 Summation: block the sum(wt * log(d)) in ormll

Fortran’s sum() (and a naive running total) accumulates N ≈ 250,000 terms of size ~10 sequentially, which alone contributes a rounding error of about 1.8e-8 to -2·logL at N = 250,000 (1.9e-7 at N = 1,000,000) — larger than the residual error left in d itself after §2.1–2.3. ormll’s new wlogsum accumulates in blocks of 512 (sum each block, then sum the block totals), which reduces this to about 1.2e-9 (250k) / 1.9e-9 (1M) — the resolution of a double at that magnitude — at negligible extra cost.

2.6 Removing redundant transcendental evaluations

• Shared exp. linkpieces computes t = exp(∓x) once and returns it for reuse by linkdpdf (the pdf-derivative) and celldiff (the exact-difference formulas for log-log/cloglog), instead of each recomputing its own exp.
• One erfc/atan instead of two. For probit and Cauchy, linkpieces computes only the smaller of F(x)/G(x)=1-F(x) directly and gets the larger one as 1 - (the smaller), which is exact to a single rounding because the larger value is always ≥ 0.5 (subtracting a number ≤ 0.5 from 1 loses no relative precision).
• expm1 with a series branch. Fortran has no expm1 intrinsic. expm1(x) = exp(x) - 1 is implemented as Kahan’s cancellation-free form (u-1)·x/log(u), u = exp(x), for general x, but for |x| < 0.02 (the common case: x here is typically the cell width or a small tail argument) a 9-term Taylor series is used instead, avoiding exp and log entirely. Verified against mpmath: relative error ≤ 2.1e-16 on both sides of the switch point, no discontinuity in accuracy at the switch.
• sinh with a series branch, same rationale. For the logistic half-width sinh(dl/2), |dl/2| < 0.02 uses a 4-term Taylor series instead of the sinh intrinsic. This is the single largest contributor to the logistic-link speedup, since a narrow cell (the common case with continuous Y) is exactly where the series applies.
• ormeta gained a what bit mask (1 = logd, 2 = g/s1/s2, 4 = h; any outputs not requested are left as passed in). The two R call sites that evaluate a trial log-likelihood only (agqStep’s inner AGQ-node loop, and orm.rfit.r’s Newton-step-size surrogate function) now pass what = 1L, skipping all derivative arithmetic; clusterModeFind passes what = 6L (needs only g, h); the sparse missing-information correction passes what = 3L (logd, g, s1, s2, not h). This required extending ormeta’s Fortran signature by one argument and updating src/init.c’s registration (19 → 20 args) and the four .Fortran(F_ormeta, ...) call sites in orm.rfit.r.

2.7 Net effect on speed (N = 250,000, single core, illustrative)

Timings vary run to run by roughly ±10-15% in the environment they were measured in and should be re-measured on the target machine; the pattern (large what=1 gains for the tail-heavy links; modest but real gains elsewhere) is the load-bearing conclusion, not the exact percentages.

Two structural options were evaluated and deliberately not adopted, for portability (plain, portable Fortran 2008 — no compiler-specific pragmas or non-default build flags):

• Manually inlining linkpieces/celldiff/linkdpdf’s bodies at their two call sites inside ormatoms (or forcing it via -finline-limit=/-finline-functions) recovers a further 5-20% on ormll what=3 and ormeta, because these routines are not inlined by gfortran at default -O2. Declined: doubles the size of the per-link branch in ormatoms’s main loop, or requires a non-default Makevars flag.
• Restructuring ormatoms around a fully-modular, argument-free internal procedure (Fortran’s analogue of an R closure — a CONTAINS-nested subroutine with no arguments, seeing the host’s local variables via host association) reproduces monolithic speed with a real per-observation subroutine boundary. This is a viable future refactor if ormatoms’s single loop needs to be split further (e.g. for YDE), but was not applied here since it wasn’t necessary yet.

2.8 Real-world confirmation: full orm fit, N = 250,000, 250,000 distinct Y levels, 10 covariates

A full orm.fit Newton fit on this design — N = 250,000, one intercept per distinct Y value (continuous Y), 10 covariates — was timed before and after this change on identical data:

Agreement between old and new fits on the same data:

• Coefficients (betas): agree to within 1e-12.
• Variance-covariance matrix elements (first intercept and all betas): agree to within 1e-11.

This is consistent with §2.7’s per-call timings (roughly 1.9x here, versus the isolated ormll/ormeta benchmarks) and with the expectation that a converged continuous-Y logistic fit should shift by an amount comparable to the accuracy improvement itself (§1.1’s ~1e-10-to-1e-2 errors in the old d, now ~1e-15), not by more. Coefficient- and covariance-level agreement at 1e-11–1e-12 confirms the new code is reproducing, not just running faster than, the old fit’s converged solution — the improved accuracy in d manifests as the fit converging to a very slightly different (more correct) optimum, at the precision level predicted, rather than as a different answer outright.

3. Differencing issues and their relevance to future YDE work

Y-dependent effects (YDE / partial proportional odds) let a subset of covariate effects vary with Y. Structurally, this means the linear predictor at each of the two intercepts in a cell becomes different:

x1 = alpha(ia)  + lp + e1,   e1 = z * gamma(ia)
x2 = alpha(ia2) + lp + e2,   e2 = z * gamma(ia2)

lp (offset + Xβ + any random effect) is common to both terms; e1, e2 are the y-dependent contributions specific to each intercept. ormatoms already has the interface this needs — this was implemented in this pass in anticipation of YDE, ahead of the modeling code itself:

• ormatoms’s signature takes lp and optional e1, e2 (default 0, i.e. today’s behavior exactly). Currently ormll passes only lp, and ormeta passes e2 = lp2 - lp1 only in the (currently never-exercised) case lp1 ≠ lp2.
• The cell width dl used in celldiff is formed inside ormatoms from (alpha(ia) - alpha(ia2)) + (e1 - e2) — i.e. from the small, y-dependent terms directly, never from x1 - x2 or from lp1 - lp2 after lp1, lp2 have already been separately rounded to doubles.

3.1 Why this matters: the caller cannot recover the width after the fact

This is the same cancellation problem as §1, one level up, and it is not fixable inside ormatoms if the caller has already destroyed the information before calling it:

• If a caller computes lp1 = lp + e1 and lp2 = lp + e2 as ordinary doubles and passes only those, the true difference e1 - e2 is gone: lp1 and lp2 are each rounded to about ulp(lp) (≈2e-15 for |lp| near 10), so lp1 - lp2 recovers e1 - e2 only to about that same absolute error — which, relative to a cell width of order 1e-5, is a ~2e-10 relative error in the width, and hence in d.
• Concretely: with lp = 10 and two candidate values of e1 differing by 8e-16 (i.e. two different, valid y-dependent effects), lp + e1 rounds to the identical double for both — the subroutine literally cannot distinguish the two cases from lp1, lp2 alone, regardless of what arithmetic it performs on them internally. No formula computed from already-rounded lp1, lp2 — including (lp1-lp2) - (alpha1-alpha2) — can recover the lost width; the information has to be preserved by the caller keeping e1, e2 separate all the way to the point where ormatoms forms alpha(ia) - alpha(ia2) + (e1 - e2) from the small numbers directly.
• Verified numerically: passing e1, e2 separately gives relative error in d of about 2-5e-15 (|lp| up to 20, cell width ~1e-5); pre-summing into lp1, lp2 and passing only those gives about 1-3e-10 at the same parameters — a ~1e5 amplification, exactly the 1/width factor from §1.

3.2 Consequence for the eventual YDE model code

Whatever R/Fortran code eventually builds lp1/lp2 (or calls ormatoms directly with e1/e2) for a YDE fit must keep the y-dependent contribution separate from the common linear predictor until it reaches ormatoms, i.e.:

• Compute lp = offset + X*beta (+ random effect) once, common to every observation and intercept.
• Compute e1 = z * gamma(ia), e2 = z * gamma(ia2) from the (typically small) y-dependent design and parameters, and pass them to ormatoms as separate arrays, not pre-added into lp.
• ormeta’s current lp1/lp2 R-facing signature is pre-summed (lp1 = lp+e1, lp2 = lp+e2 already added on the R side) — this is fine today only because e1 = e2 = 0 always holds (no YDE yet, so lp1 = lp2 exactly, and the cancellation issue is moot). Extending ormeta to real YDE use will require either (a) changing its R-facing signature to accept lp, e1, e2 separately (preferred, avoids the problem entirely), or (b) accepting the residual ~ulp(lp)-scale error that comes from ormeta internally re-deriving e2 = lp2 - lp1 from pre-summed inputs, which is what it does today as a stopgap.

3.3 What YDE additionally needs from the derivative structure

With e1 ≠ e2, x1 and x2 no longer move together under a YDE parameter’s gradient direction:

• s1, s2 (returned by ormatoms at what ≥ 2) are already exactly D(ell)/D(x1) and D(ell)/D(x2) individually — not just their sum — so they are usable as-is: a YDE score contribution is s1*z1 + s2*z2 for the observation’s y-dependent design vectors z1, z2 at the two intercepts.
• Likewise h11, h12, h22 (the second derivatives in x1, x2 separately) are what a YDE Hessian needs, rather than the combined hee, h1e, h2e, which describe a single, common shift of both x1 and x2 (valid for β and random-intercept parameters, where the same shift really does apply to both terms) and are not the right quantities once e1 ≠ e2. This distinction is called out explicitly in ormatoms’s header comment and in ormll’s Hessian-derivation comment, specifically so that whoever writes the YDE accumulation code does not reach for hee/h1e/h2e by analogy with the existing β/random-effect code and get a silently wrong Hessian.
• None of ormatoms’s per-link formulas (§2.1–2.3) needed to change to support this — linkpieces and celldiff already operate on x1, x2 as independent arguments; the only assumption removed was that they share one lp.

3.4 What remains unresolved for YDE, and is out of scope here

• Underflow, not just cancellation. For log-log and cloglog, a cell whose true probability underflows double precision (below ~1e-308) still returns salloc = 999 — accurate-but-tiny is not the same as representable. A YDE fit that pushes some cells into extreme tails (more likely with more free parameters) could hit this more often. A log-space version of celldiff’s formulas exists in principle but was not built, since it would also require the score/Hessian formulas to be re-derived in log space to be useful (the current s1 = pdf/d would still overflow/underflow even with an accurate log(d)).
• Probit’s residual cancellation (§2.3) is unchanged and would still apply, at the same magnitude, to YDE probit fits with narrow cells.
• No YDE modeling code (R or Fortran) exists yet. This document describes the interface and numerical groundwork laid in ormatoms in anticipation of that work, not an implementation of YDE itself.