dnorm() operation ordering differs between RTMB and TMB

[divit...@gmail.com]

Hi all,

While porting the Gaussian likelihood in glmmTMB from TMB to RTMB, I noticed a small but reproducible numerical difference between otherwise equivalent models. After some debugging, I traced it to the order of operations used in dnorm().

RTMB currently computes the log-density as

ans <- -0.5 * r * r - log(sqrt(2 * pi)) - log(sd)

whereas TMB computes

-log(sqrt(2*M_PI)) - log(sd) - 0.5 * resid * resid

These expressions are mathematically equivalent, but because floating-point subtraction is not associative, they can differ by about 1e-13 in the objective. In my case, that was enough for nlminb() to follow a slightly different optimization path, leading to slightly different fitted coefficients.

Changing my RTMB implementation to match TMB's operation ordering eliminated the discrepancy and produced identical objectives, gradients, coefficients, and log-likelihoods.

My questions are:
1. Is the difference in operation ordering between RTMB::dnorm() and TMB::dnorm() intentional?
2. If so, is there a particular reason for preferring the current RTMB ordering?
3. Would it make sense for RTMB::dnorm() to use the same operation order as TMB::dnorm() to improve numerical reproducibility between equivalent RTMB and TMB models?

I've written up the details, including the relevant source code, numerical example, and test results, here:

https://github.com/glmmTMB/glmmRTMB/issues/4

Thanks!

Best regards,
Divit Purwar