R Python differences

R / Python numerical equivalence

This page documents where the R smooth package and the Python port produce bit-identical results and where they don't. It is the running record of the cross-language status as the two implementations converge.

What is identical (machine precision)

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

Quantity	Across all scenarios
coef(m) / m.coef — fitted coefficient vector	✓
m$lossValue / m.loss_value — cost-function value	✓
m$logLik / m.loglik	✓
m$fitted / m.fitted — fitted values	✓
m$residuals / m.residuals	✓
forecast(m, h=…) / m.predict(h=…) — point forecasts	✓
sigma(m) / m.sigma — empirical residual std (R formula)	✓
AIC / AICc / BIC / BICc	✓
vcov(m), confint(m), summary(m) — OM, OMG, and ADAM with initial="backcasting" / "complete"	✓
vcov(m) / confint(m) — ADAM with initial="optimal" / "two-stage"	≤2% — see below
coefbootstrap() — structural parity (same coef cardinality, PSD vcov)	✓
multicov(type="analytical")	≤5%
multicov(type="empirical") — calls the same C++ ferrors backend on both sides	✓

What had to be aligned

The shared C++ core and the matching call sites:

1. Shared finite-difference Hessian — src/headers/hessianCore.h is the single source of truth for the FD Hessian used by both R (vcov.adam / vcov.om / vcov.omg) and Python (numerical_hessian). Per-parameter relative step hᵢ = ε^(1/4) · max(|xᵢ|, 1) so large-magnitude parameters (initial level / trend / seasonal in B when initial="optimal" or "two-stage") get a meaningful perturbation.

2. bounds="none" during FI computation. Both implementations disable the cost-function bounds penalty during the FD Hessian so perturbations at boundary parameters don't trip the 1e+300 infeasibility return.

3. abs(diag(...)) of vcov. Mirrors R/adam.R:5226 — both sides agree on the sign convention for non-PSD FIs.

4. sigma() formula. Python ADAM.sigma computes sqrt(SS / df) exactly per R's sigma.adam (R/adam.R:4625-4658), with a distribution-specific SS and df = nobs - nparam (Python's nparam already excludes the scale parameter that R counts and then subtracts). The internal optimisation scale R calls m$scale is exposed as m.scale (separate from m.sigma).

5. OM.sigma / sigma.om. OM has no scale parameter on the probability axis; both languages report sqrt(mean(residuals²)) on the link-transformed scale — formula from oes_old / oesg_old.

6. dgnorm shape lower bound. Python's optimiser lower bound on the generalised-normal shape is 1e-10, matching R's adam_checkOptimizer. A tighter bound (0.25) was previously used and changed the NLopt Nelder-Mead simplex behaviour, producing different local minima on the same cost surface.

7. two-stage ≡ optimal for FI. Two-stage produces the same B shape as optimal; both R (R/adam.R near the initialTypeFI block) and Python treat them identically when computing the Hessian.

The one remaining numerical gap

The ULP-level gap on ADAM with initial="optimal" (or "two-stage") on multi-parameter ETS models comes from the linear-least-squares solver used inside msdecompose's default global smoother — and nowhere else.

• R uses .lm.fit → LINPACK Householder QR (dqrls.f).
• Python uses numpy.linalg.lstsq → LAPACK dgelsd (SVD-based).

Both solve the same least-squares problem; they round at different bits during elimination/back-substitution. Result: msdecompose's level (and to a much smaller extent trend) differs by ~1 ULP on the data scale (level ≈ 100 for AirPassengers → ~1.4e-14 absolute difference).

That ULP enters the x0 passed to NLopt only when level / trend are in B — i.e. for initial="optimal" and "two-stage". The backcasting and complete paths derive these inside the C++ kernel from the data each call, so the LSQ-solver ULP never reaches the optimiser. OM and OMG likewise don't carry initial states in B. The optimiser is deterministic given the same x0, so this single ULP is the entire source of downstream divergence in vcov / confint / summary for optimal / two-stage.

The downstream amplification chain — LSQ ULP → x0 ULP → final-B ULP → FI numerator → inverse-FI — used to compound this into ~5e-3 diffs on the worst vcov entry. With the per-parameter relative Hessian step (point 1 above) the FI amplification is now small enough that the gap collapses to roughly 2% relative on vcov / SE / confint for the optimal-initials path.

This is a property of the BLAS/LAPACK implementation, not an algorithmic bug. Two different LSQ solvers will always round differently. Closing the residual gap would require either (a) using a closed-form formula for the X = [1, t] case inside msdecompose (intercept + time — trivially bit-identical across implementations) or (b) routing msdecompose's smoother through a shared C++ QR via the same Rcpp / pybind11 bridge used for the Hessian.

Current parity-test tolerances

The Python repo's r_parity tests (test_om_summary_r_comparison.py, test_omg_r_comparison.py, test_adam_summary_r_comparison.py, test_fi_r_comparison.py, test_multicov_r_parity.py, test_coefbootstrap_r_parity.py) carry the r_parity marker and are deselected in CI by default.

Quantity	rtol	atol
Coefficients (all scenarios)	machine precision	—
Log-likelihood, loss, fitted, residuals	machine precision	—
OM / OMG / backcasting-ADAM vcov / SE / CI	1e-4	1e-6
ADAM optimal / two-stage vcov, SE, CI bounds	2e-2	1e-3
multicov(type="analytical")	5e-2	1e-3
multicov(type="empirical")	5e-2	1e-3
coefbootstrap (stochastic — distributional only)	structural