R Python differences

R / Python numerical equivalence — known divergences

This page documents the cross-language numerical equivalence between the R smooth package and its Python port. It is the result of a detailed investigation that traced every cross-language gap down to its arithmetic origin. Where the two implementations diverge, this page explains why.

The short version:

Scenario	Equivalence
`OM` (any model / occurrence)	machine precision (~1e-9 or 0)
`OMG` (any combination)	machine precision (~1e-9 or 0)
`ADAM` with `initial="backcasting"` or `"complete"`	machine precision (0 or ~1e-11)
`ADAM` with `initial="optimal"`, multi-parameter scenarios	~1e-5 in FI, ~1e-3 in vcov — see below

The single remaining cross-language gap is in multi-parameter ADAM with initial="optimal" and is fully diagnosed below.

What is at machine precision

The following are bit-identical or agree to within a handful of ULPs between R and Python on the same input series:

Final coefficient vector coef(m) / m.coef
Log-likelihood logLik(m) / m.loglik
Loss value m$lossValue / m.loss_value
Fitted values m$fitted / m.fitted
Residuals m$residuals / m.residuals
Forecasts forecast(m, h=...) / m.predict(h=...)
Information criteria (AIC, AICc, BIC, BICc)
vcov(m), confint(m), and the contents of summary(m) in all scenarios except the one called out below.

What this required

Reaching machine-precision parity took several specific alignments, specifically:

Shared C++ numerical Hessian. Both languages now call into the same src/headers/hessianCore.h (via Rcpp on R, pybind11 on Python). R no longer depends on pracma::hessian. Same loop, same step size, same operation order — the Hessian wrapper itself contributes zero divergence.
bounds="none" during FI computation. Both implementations now disable the cost-function bounds penalty during the finite-difference Hessian. Without this, perturbations at boundary parameters (e.g. alpha=0) trip the 1e+300 infeasibility return and the inverse FI collapses to zero.
abs(diag(...)) of vcov. Mirrors R/adam.R:5226 ("just in case, take absolute values for the diagonal") — both sides agree on the sign convention for non-positive-semi-definite FIs.

The one remaining gap: ADAM with `initial="optimal"`

Symptom

For adam(y, model="AAN", initial="optimal") on a 120-observation series, R and Python produce:

Quantity                R                          Python                     Diff
coef (alpha)            0.96308367588546961        0.96308367588546961        0
coef (beta)             0.12202608739523153        0.12202608739523153        0
coef (level)            100.5313366407765          100.53133664077598        -5.26e-13
coef (trend)           -0.6389674596160431        -0.63896745961606904       -2.60e-14
logLik                 -168.84126131398833        -168.84126131398853         2.00e-13
vcov[0,0] (alpha,alpha) 0.0071360323…             0.0071358636…              ~1.7e-7
…                       …                          …                          …
max |vcov_py − vcov_r|                                                        ~5.06e-03

Same series, same model spec, same initial="optimal", but the level coef ends 1 ULP apart and the vcov on the worst entry differs by 5e-3.

The diagnosis

I traced this exhaustively. Every step is reproducible from the investigation log in the Python repo's git history; the summary follows.

Step 1 — Same optimiser (NLopt) but different versions

Python: bundled libnlopt in the wheel — version 2.10.0
R:      system /usr/lib/x86_64-linux-gnu/libnlopt.so.0 — version 2.7.1

The two NLopt versions are 3+ years of development apart. But: when given the same starting point x0, both versions converge to bit-identical final B. Version difference is not the cause.

Step 2 — The optimiser settings are identical

xtol_rel = 1e-6
xtol_abs = 1e-8
ftol_rel = 1e-8
ftol_abs = 0
maxeval  = length(B) * 40 = 160
algorithm = NLOPT_LN_NELDERMEAD

R's utils-adam.R:237-240 and Python's adam_general/core/estimator/optimization.py:180-182 agree to the literal value. Tolerances are not the cause.

Step 3 — The initial point `x0` already differs

Captured the exact x0 each implementation passes into nlopt::optimize:

alpha:  py=0.10000000000000001       r=0.10000000000000001       diff= 0
beta:   py=0.050000000000000003      r=0.050000000000000003      diff= 0
level:  py=96.513677541526278        r=96.513677541526292        diff=-1.42e-14   ← seed
trend:  py=0.11740463161952641       r=0.11740463161952630       diff=+1.11e-16

Before the optimiser runs, R and Python disagree on level by ~1 ULP at scale 100. That ULP is the seed of all downstream divergence.

Step 4 — Forcing Python to start from R's exact x0 produces bit-identical convergence

Monkey-patched nlopt.opt.optimize to replace Python's x0 with R's verbatim, then re-ran. The result:

                   Final coef diff vs R
Python default x0:       5.26e-13
Python with R's x0:      0.00e+00   ← BIT-IDENTICAL

This is the decisive evidence: NLopt is deterministic given the same x0. The cross-language divergence comes from the x0 generator, not the optimiser.

Step 5 — The x0 comes from `msdecompose`

R's R/utils-adam.R:666-678 and Python's adam_general/core/utils/utils.py both call msdecompose(y_in_sample, lags=1, type="additive", smoother="global") to seed level and trend. Compared the raw output:

Python msdecompose level: 96.513677541526278
R      msdecompose level: 96.513677541526292   diff = 1.42e-14

The exact same 1 ULP of divergence shows up at the msdecompose output, so the cause is upstream of msdecompose's post-processing. It is inside the smoother.

Step 6 — The smoother uses different least-squares solvers

The default smoother (smoother="global") fits a linear regression on the input series. The implementations:

R (R/msdecompose.R:96-113):

X <- cbind(1L, seq_len(n))
return(y - .lm.fit(X, y)$residuals)

.lm.fit is R's fast lean LM, backed by LINPACK's QR via dqrls.f.

Python (adam_general/core/utils/utils.py:283-300):

X = np.column_stack([np.ones(n), np.arange(1, n + 1)])
coef = np.linalg.lstsq(X, y, rcond=None)[0]
return X @ coef

np.linalg.lstsq is backed by LAPACK's dgelsd (SVD-based).

Both solve the same least-squares problem. Different numerical algorithms — LINPACK Householder QR vs LAPACK SVD — round at different bits during the elimination/back-substitution. Empirically:

Method                                    Max diff of fitted values vs R
np.linalg.lstsq (LAPACK SVD, current)             7.11e-14
np.linalg.lstsq → y − residuals (same path)       7.11e-14
np.linalg.solve(X.T @ X, X.T @ y) (normal eqs)    9.95e-14
scipy.linalg.lstsq(driver='gelsy') (QR+pivot)     4.26e-14  ← closest

None of NumPy / SciPy's LSQ paths match R's LINPACK QR bit-exactly. Even two different QR implementations (LINPACK Householder vs LAPACK QR with column pivoting) round at different bits. This is a property of the BLAS/LAPACK implementation, not an algorithmic bug.

The propagation chain

That ~7e-14 ULP per fitted value in the smoother propagates as follows:

LSQ solver (LINPACK QR  vs  LAPACK SVD)
        │  ~7e-14 ULP per fitted value
        ▼
smoothing_function_global  trend  array
        │  same ULP scale
        ▼
msdecompose  initial$nonseasonal (level / trend)
        │  partial cancellation in level = trend[0] − mean(diff(trend)) * lagsMax
        ▼  ~1.42e-14 in level,  ~1.11e-16 in trend
nlopt::optimize  starts Nelder-Mead from this x0
        │  ~160 iterations, ULP-different simplex trajectories
        ▼  ~5.26e-13 in final B  (level)
Finite-difference Hessian at non-identical B
        │  multiplied by 1/h² ≈ 6.7e7
        ▼  ~3e-5 in FI matrix entries
vcov  =  FI⁻¹    (with the 4×4 FI ill-conditioned for AAN)
        ▼  ~5e-3 max in the worst vcov entry

Each arrow is a real, identifiable amplification factor.

Why this only shows up in `initial="optimal"`

initial="backcasting" / "complete": level and trend are not in B. They are derived inside the C++ kernel by backcasting from the data each call. msdecompose is not the source of any value in B. Therefore the LSQ-solver-ULP never enters the optimiser's x0.
OM and OMG: same — B holds only smoothing parameters (α, β, γ, …), not initial states. msdecompose is not in the path.
initial="optimal": level and trend are estimated as elements of B. Their starting values come straight from msdecompose$initial$nonseasonal. The 1 ULP from the LSQ solver enters B[level] and propagates through the optimiser into the FI.

This explains why OM, OMG, and backcasting-ADAM are at machine precision while only initial="optimal" ADAM shows the residual.

A separate finding: R's reported optimum is sometimes a saddle

While verifying the explanation above, an additional finding emerged.

When Python's optimiser is started from R's converged model$B (not R's x0), Python's Nelder-Mead does not stay there. It walks away and finds a better local minimum:

                          loss
R   model$B (claimed optimum):     168.84126131398833
Python re-optimises from there:    164.75966898758691   ← better by ~4 units

The coefficient vectors at these two points differ by 0.72 in the worst entry — a real, multi-component move, not noise.

R itself emits the warning at vcov(m) time:

Observed Fisher Information is not positive semi-definite, which might
mean that the likelihood was not maximised properly.

That warning is consistent with what happened: R's nloptr settled at a saddle of the loss surface — a point where some Hessian eigenvalues are positive and others negative — rather than a strict minimum. Python's NLopt, given a "kicked" starting point at that saddle, escapes downhill along the negative-eigenvalue direction.

Both implementations encounter this saddle when running their own default workflow (starting from msdecompose x0). Neither finds the better minimum in normal use; both stop on the same suboptimal ridge. The "better minimum" is only reachable by deliberately displacing the starting point.

This is not directly related to the cross-language gap, but it explains why the FI computed at the converged point is non-PSD on this scenario. Once you compute the FI at a saddle, the diagonal becomes a mix of positive and negative curvature, the inverse is ill-conditioned, and small ULPs in the FI inputs amplify into large ULPs in vcov.

What could close the remaining gap

The root cause is the LSQ solver inside msdecompose's default global smoother. Options:

Option A — Closed-form linear regression (recommended)

For the default X = [1, t] case (intercept + time), the linear regression has a closed-form solution that involves only sums and products — no QR, no SVD, no BLAS dependency. The formula is the same on both sides and computable bit-identically:

t = np.arange(1, n + 1, dtype=float)
t_bar = (n + 1) / 2
y_bar = y.mean()
# Equivalent in R: same arithmetic, same order.
slope = ((t - t_bar) * (y - y_bar)).sum() / ((t - t_bar) ** 2).sum()
intercept = y_bar - slope * t_bar
fitted = intercept + slope * t

Implementing this in both R's msdecompose.R and Python's adam_general/core/utils/utils.py would make the smoother's output bit-identical across languages, and by Test B above, the entire ADAM optimiser convergence would follow. initial="optimal" would then reach machine precision parity.

The block-dummies fallback (n_groups > 1) covers a less-common code path and could either stay on the implementation-specific LSQ for now or be migrated to a shared C++ QR via the same pattern used for the Hessian.

Option B — Shared C++ QR via the existing Rcpp / pybind11 bridge

Same architecture as the shared hessianCore.h: write a small shared QR-based LSQ in src/headers/, expose it to both languages, and route smoothing_function_global through it on both sides. Guarantees bit-identical output across all column counts. More work than Option A for the same observable benefit on the default path; only worth doing if the block-dummies fallback is also a parity concern.

Option C — Accept the floor

Document the gap (this page) and the parity tests' tolerances around the LSQ-noise floor (atol=1e-3 on vcov for ADAM-with-initial="optimal"). Both implementations are independently correct under their respective numerical-algebra libraries; the gap is fundamental to LINPACK vs LAPACK/OpenBLAS rounding.

Tolerances in the parity test suite

The Python repo's tests/test_om_summary_r_comparison.py, test_om_r_comparison.py, test_omg_r_comparison.py, test_adam_summary_r_comparison.py, and test_fi_r_comparison.py carry the r_parity marker and are deselected in CI by default. The tolerances are calibrated to catch regressions, not to enforce bit-equivalence where the LSQ-solver floor is unavoidable:

Quantity	rtol	atol
Coefficients	1e-2	1e-3
Log-likelihood, loss, fitted	1e-12 (typically; sometimes tighter)	1e-13
OM / OMG / backcasting-ADAM vcov	1e-4	1e-6
ADAM-`optimal` vcov	0.30	1e-2
ADAM-`optimal` SE	0.20	2e-2
ADAM-`optimal` CI bounds	0.10	5e-2

If Option A or B is implemented in the future, the ADAM-optimal tolerances can be tightened to match the machine-precision band of the other rows. The tests will then guard the tighter floor.

Reproducibility / where to run the experiments

The investigation that produced this page used:

Local R checkout (devtools::load_all(".")) — necessary so that the smooth R code matches the current source rather than CRAN.
Local Python checkout (editable install via pip install -e .).
An ad-hoc nlopt.opt.optimize monkey-patch to inject specific starting points into Python's optimiser without changing its public API.
The _r_bridge.py test helper in the Python repo for round-trip JSON marshalling between Python and R.

If revisiting this analysis, the three decisive measurements are:

Capture x0 at the boundary between msdecompose output and nlopt::optimize input on both sides. Confirm any sub-ULP difference in level.
Inject R's x0 into Python's optimiser. Confirm bit-identical final B.
Compare msdecompose$initial$nonseasonal (R) to msdecompose(...)['initial']['nonseasonal'] (Python) on the same series. The 1 ULP in level reproduces every time.

If any of those three break in the future, the cause has moved and the chain in this document needs re-tracing.

Uh oh!

R Python differences

R / Python numerical equivalence — known divergences

What is at machine precision

What this required

The one remaining gap: ADAM with initial="optimal"

Symptom

The diagnosis

Step 1 — Same optimiser (NLopt) but different versions

Step 2 — The optimiser settings are identical

Step 3 — The initial point x0 already differs

Step 4 — Forcing Python to start from R's exact x0 produces bit-identical convergence

Step 5 — The x0 comes from msdecompose

Step 6 — The smoother uses different least-squares solvers

The propagation chain

Why this only shows up in initial="optimal"

A separate finding: R's reported optimum is sometimes a saddle

What could close the remaining gap

Option A — Closed-form linear regression (recommended)

Option B — Shared C++ QR via the existing Rcpp / pybind11 bridge

Option C — Accept the floor

Tolerances in the parity test suite

Reproducibility / where to run the experiments

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Clone this wiki locally

The one remaining gap: ADAM with `initial="optimal"`

Step 3 — The initial point `x0` already differs

Step 5 — The x0 comes from `msdecompose`

Why this only shows up in `initial="optimal"`