Skip to content

R Python differences

Jump to bottom
Ivan Svetunkov edited this page Jun 16, 2026 · 10 revisions

R vs Python differences

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) / SE — ADAM with initial="optimal" / "two-stage", evaluated at the same B ≤1e-3 (FD-Hessian discretisation floor)
coefbootstrap() — structural parity (same coef cardinality, PSD vcov) ✓
multicov(type="analytical") ≤1e-10 (was ≤5% — shared olsCore.h closed the optimiser-drift cascade)
multicov(type="empirical") — calls the same C++ ferrors backend on both sides ✓

The optimal / two-stage ADAM vcov used to carry a ~2% relative gap caused by the msdecompose global-smoother LSQ ULP. That gap has been closed by routing both languages through a shared src/headers/olsCore.h (pivoted QR + scale-invariant rank cutoff via Armadillo). See "What used to be the one remaining numerical gap" below.

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.

  8. Shared C++ OLS for msdecompose — src/headers/olsCore.h exposes olsCore(X, y, tol) (pivoted QR with a LINPACK-dqrls- style scale-invariant rank cutoff: a column is dropped when |R(i,i)| / ‖X.col(P(i))‖ < tol, default 1e-7). It is the single source of truth for the two OLS calls inside msdecompose (global-smoother trend fit and NA imputation). The R build wraps it as olsCpp (src/olsWrap.cpp, Rcpp) and the Python build wraps it as smooth.adam_general._ols.ols (src/python/olsWrap.cpp, pybind11 + carma, registered in python/CMakeLists.txt). This closes the historical LSQ-solver gap described in the next section.

What used to be the one remaining numerical gap

The ULP-level gap on ADAM with initial="optimal" (or "two-stage") on multi-parameter ETS models used to come from the linear-least- squares solver used inside msdecompose's default global smoother:

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

Both solve the same least-squares problem; they round at different bits during elimination / back-substitution. The resulting ULP entered the x0 passed to NLopt for initial="optimal" / "two-stage" (the backcasting / complete paths derive the level / trend inside the C++ kernel each call, so they were never affected). Downstream amplification — LSQ ULP → x0 ULP → final-B ULP → FI numerator → inverse-FI — used to compound this into a ~2% relative gap on the worst vcov / SE entry, even after the per-parameter relative Hessian step in hessianCore.h did most of the work.

That gap is now closed. The shared src/headers/olsCore.h runs the same pivoted QR + rank cutoff on both sides, so the OLS step itself produces the same level / trend seed to within a handful of ULPs (modulo BLAS-dispatch nuances inside Armadillo's geqp3). On the same five-scenario R-parity bench used by test_adam_summary_r_comparison.py (ANN, AAN, AAdN with bounds="admissible", ARIMA(1,0,1), and the two-stage AAN), inverting the FI at R's coefficients and comparing to R's vcov(m) gives:

Scenario max |ΔSE / SE| (was ≤2%) max |Δvcov / vcov| (was ≤2%)
ANN 6.4e-9 2.1e-7
AAN 1.6e-8 8.4e-8
AAdN (admissible) 5.0e-7 1.5e-6
ARIMA(1,0,1) 1.6e-4 3.3e-4
AAN, two-stage 3.7e-8 4.1e-7

That is a 4–7 order-of-magnitude tightening. The residual is now bounded by the FD-Hessian discretisation floor (step h ≈ ε^(1/4) ≈ 1.2e-4, central-difference error O(h²) ≈ 1e-8), not by the LSQ solver. The R-parity test (test_adam_summary_r_comparison.py) has been tightened from rtol=2e-2 → rtol=1e-3 on the at-the-same-B vcov / SE comparison to lock the new bound in.

msdecompose summation-order parity (ETS(M,M,M) on high-frequency seasonals)

A second msdecompose gap — unrelated to the OLS one above — caused Python's ETS(M,M,M) fit on taylor (n = 4032) to land in different NLopt basins from R on chaotic configurations (single seasonal lag in {48, 49, 50, 96}, and the combined lags=c(48, 336)). The root cause was that three reductions inside msdecompose did their summations in a different IEEE-754 order from R: the seasonal-centring mean and the trend-from-diff mean used NumPy's pairwise sum while R's mean() uses an 80-bit LDOUBLE accumulator; and — most importantly — the inner seasonal-mean was Python's np.mean, whereas R goes through stats::filter(weights = rep(1,N)/N) whose C-level cfilter walks the input from the last element down to the first, accumulating value · (1/N) in plain IEEE doubles. Those sub-ULP rounding-direction differences in the seasonal initials were amplified by the undamped multiplicative-trend recursion (4032 chained multiplicative steps) into NaN states on the Python side, tripping cost_functions.py's NaN → 1e+300 cap at the very first NLopt eval and blocking the optimiser from descending; with R's exact msdecompose bytes forced in, the same Python pipeline reached R's logLik byte-for-byte, proving msdecompose was the sole differentiator. The fix: two new helpers in python/src/smooth/adam_general/core/utils/utils.py — _fsum_nanmean (Shewchuk exact, matches R's LDOUBLE mean()) at the two centring sites, and _r_filter_mean (the backward value · inv_n accumulator) at the inner seasonal-mean site that mirrors R's cfilter byte-for-byte. After the fix every msdecompose output (level, trend, all seasonal patterns) is byte-identical to R on taylor and every ETS configuration previously affected now converges to R's logLik to 4 decimal places.

Configuration R native logLik Py before (np.mean) Py after (_r_filter_mean)
MMM, lags=[48] −29 096.16 −28 394.25 (stuck at x0) −29 096.16 ✓
MMM, lags=[49] −30 774.02 −33 478.48 −30 774.02 ✓
MMM, lags=[50] −30 748.68 −33 494.79 −30 748.68 ✓
MMM, lags=[96] −28 390.19 −28 387.40 −28 390.19 ✓
MMM, lags=[48, 336] −27 634.01 −73 254.19 (stuck at x0) −27 634.01 ✓
MMdM, lags=[48, 336] −27 633.876 −27 633.868 −27 633.876 ✓
MNM, lags=[48, 336] −27 634.44 −27 634.44 −27 634.44 ✓
Per-site fix R operation Python before Python after
Seasonal-centring (:402) mean(patterns[[i]][1:N], na.rm=T) (LDOUBLE) np.nanmean _fsum_nanmean
Initial-trend (:418) mean(diff(ySmooth), na.rm=T) (LDOUBLE) np.nanmean _fsum_nanmean
Inner seasonal-mean (:387) stats::filter(weights = rep(1,N)/N) (IEEE, backward scalar · 1/N) np.mean (pairwise) _r_filter_mean

Damped multiplicative trend (MMdM) and additive-trend-free (MNM) configurations were never sensitive enough to the ULP differences to diverge in the first place — only undamped MMM (φ = 1) amplifies the chaos at this scale. The _r_filter_mean helper runs in a Python loop but msdecompose only executes once at init time and each seasonal slice has ≤ ⌈n/lag⌉ elements, so the per-fit cost is negligible.

OM Bernoulli likelihood + FI Hessian (intermittent demand)

Two OM-side discrepancies caused the cross-language tests to disagree on intermittent demand fits where the Bernoulli log-likelihood is flat near the optimum. (1) om_cost.om_cf aggregated log(p) and log(1-p) via np.sum while R's om() uses an LDOUBLE-accumulator sum() — the 1-ULP CF difference at iter 0 steered NLopt into a different basin (e.g. on the seasonal MNM odds-ratio case R landed at alpha = 5.87e-5 while Py converged to the corner alpha = 0, producing 3 % relative differences in fitted probabilities). Fixed by routing both terms through math.fsum. (2) OM._fisher_information_matrix was running its central-difference Hessian off self._adam_created["mat_vt"] and self._profile["profiles_recent_table"] — both of which the C++ kernel's head-fill / backcasting passes mutate in place on every CF evaluation, so by the time vcov() runs both buffers reflect the post-fit state. R re-creates fresh matrices for FI via adam_creator + om_initial_transform on every Hessian probe (R/om.R:830-870). Python now snapshots the buffers before NLopt starts and restores those bytes inside the FI _cost closure on every probe; after the fix every OM r-parity test in test_om_r_comparison.py and test_om_summary_r_comparison.py passes at the existing rtol=1e-4 / atol=1e-6 tolerances.

CES "partial" / "full" seasonality (stats::lowess ULP gap)

CES with seasonality="none" and "simple" was bit-perfect against R from the start. "partial" and "full" used to diverge — on taylor-style data this manifested as a small ARMA-style coefficient gap; on the ces_quarterly benchmark full_quarterly was 2 logLik units off and the BOBYQA trajectory split visibly around the 25th evaluation. Traced step-by-step:

  1. The first ~24 BOBYQA evaluations were bit-identical between R and Python.
  2. At a degenerate B = [1.3, 0, 1.3, 1] the CF values diverged by 4 ULPs (1003.8731415564403 vs 1003.8731415564408). Every Python summation method (np.sum, math.fsum, np.longdouble) gave the same Py value, so the aggregator was innocent.
  3. Byte-comparing the kernel inputs at that B: matWt, matF, vecG, vectorYt, vectorOt were bit-identical. matVt and profilesRecent differed in the seasonal-component initial state rows by ~3.5e-14 (~25 ULPs).
  4. The seasonal-component seed in CES came from msdecompose(y, lags=…, type="additive")$seasonal with smoother="lowess". R uses stats::lowess (the original Fortran routine wrapped through C). Python uses greybox.lowess (a C++ port via pybind11). The two implementations agree at the algorithm level but diverge by ~3e-13 in their outputs on real-world data, with 46 of 80 cells differing on the ces_quarterly series.

Fix: switched CES (both R/adam-ces.R and Python core/ces/creator.py) to call msdecompose with smoother="global" instead of the historical default "lowess". The shared src/headers/olsCore.h backend behind smoother="global" produces byte-identical decompositions across languages (already verified for the ADAM / msdecompose r-parity work in the section above). After the swap the CES r-parity failure count drops from 84 to 18 — every airpassengers case becomes bit-perfect, and full_quarterly becomes bit-perfect. The remaining 18 failures are all quarterly cases where a residual 1-ULP cell-difference in msdecompose's global path cascades through the CES two-stage NLopt path into a ~2.7e-5 logLik gap — true optimiser-floor effect, well inside the realistic parity envelope.

Pinning down a fully bit-perfect lowess port (i.e., a faithful port of clowess.f rather than greybox's current C++ rewrite) remains tracked as a follow-up — useful if a future caller wants the lowess smoother explicitly — but is no longer on the CES critical path.

What still moves at the optimiser-convergence floor

confint() from Python's ADAM.fit() runs its own NLopt fit, so the reported bounds = py_coef ± py_SE × t_crit against R's r_coef ± r_SE × t_crit. Python's NLopt and R's nloptr converge to slightly different points on flat parts of the cost surface (worst on the ARIMA scenario, where the relative coefficient gap is up to ~1.5% on the AR / MA terms even at the same loss_value). This is an optimiser-floor effect, not an OLS effect — the per-scenario test_confint_matches_r therefore stays at rtol=2e-2 on the confidence bounds, while the at-R's-B vcov / SE comparison runs at the new rtol=1e-3.

Simulation functions and RNG-induced differences

This section covers the sim_* family (sim_es, sim_ssarima, sim_sma, sim_oes, sim_ces, sim_gum) and the simulate() methods (ADAM.simulate, OM.simulate, OMG.simulate).

Default seeds do NOT cross languages

R and Python use different random number generators:

  • R: the Mersenne Twister (set.seed(n)) for rnorm, rt, rgnorm, etc.
  • Python: NumPy's PCG64 default (np.random.default_rng(n)).

Even with identical seeds, the two languages produce completely different draws. The same sim_es(model="ANN", seed=42) call on either side gives an entirely different series. The same applies to ADAM.simulate(seed=42), OM.simulate(seed=42), and so on.

What is preserved across the R/Python boundary is the state-space machinery that consumes those draws — the matrix prep, the C++ adamCore::simulate kernel call, the post-processing (link functions for OM/OMG, scaling for ETS, etc.). These are bit- equivalent on both sides. The only divergence is the random draws the kernel ingests.

The "plug-in numbers" trick

To prove the state-space machinery itself matches across languages, all of the sim_* functions and the simulate() methods accept a callable in place of the distribution name:

# Python
errors = np.linspace(-1.5, 1.5, 60)
def feed(n):
    return errors[:n]
sim = sim_es(model="ANN", obs=60, randomizer=feed, ...)
# R — same vector
errors <- seq(-1.5, 1.5, length.out=60)
ourFunction <- function(n, ...) errors[1:n]
sim <- sim.es(model="ANN", obs=60, randomizer="ourFunction", dummy=1)

Both calls feed the same errors into the same C++ kernel, so the resulting series, states, and persistence vectors agree to floating-point. The r_parity-marked tests in the Python repo (test_sim_es_r_parity.py, test_sim_ssarima_r_parity.py, test_sim_sma_r_parity.py, test_sim_oes_r_parity.py, test_sim_ces_r_parity.py, test_sim_gum_r_parity.py) all use this trick and pass at atol=1e-9 or tighter:

Function Parity test Tolerance
sim_es ANN, AAA (with seasonal initial) 1e-10
sim_ssarima AR(1), IMA(1,1) 1e-9
sim_sma order 3 1e-10
sim_oes odds-ratio, inverse-odds-ratio, direct, general 1e-10
sim_ces none, full 1e-9
sim_gum local-level, 2-component 1e-10
OMG.simulate combined-probability 1e-9

Known divergences in OM.simulate / simulate.om

Three R-parity tests on the simulate-from-fit path are currently marked xfail with explanations:

  1. test_om_simulate_latent_matches_r — OM fits on the same y converge to slightly different smoothing parameters between Python (NLopt) and R for small intermittent series. With matched errors via the plug-in trick the latent series differ by the ratio alpha_python / alpha_R. This is the underlying OM-fit-parity issue, not a simulate-path divergence.

  2. test_om_simulate_seasonal_matches_r — multiplicative- seasonal OM (MNM, lags=[1, 12]) produces all-NaN latent on the Python side because the multiplicative seasonal state can collapse to zero under matched errors. R handles the same matrix prep without NaN. Tracked as a Python-side ADAM.simulate multiplicative-seasonal numerics follow-up.

  3. test_om_simulate_arima_matches_r — R's simulate.om on an OM with ARMA orders crashes inside the C++ kernel (element-wise multiplication: incompatible matrix dimensions: 2x1 and 1x1). An upstream R-side bug in simulateADAMCore's matrix prep for OM+ARIMA. Python's OM.simulate handles the same case cleanly (the smoke test passes); byte-equivalence is impossible until the R bug is fixed.

The OMG combined-probability case (test_omg_simulate_combined_probability_matches_r) does pass at atol=1e-9 because the OMG combination formula (omg_link_function) operates on the latent series and the optimizer divergence on a single OM fit gets averaged out when combining two independently-fit sub-models.

Practical guidance

When you want to reproduce a specific simulation across both languages:

  • Pre-generate the error vector in either language (it doesn't matter which — both have access to the other's outputs via CSV / JSON / R_data round-trip).
  • Pass it via the randomizer callable form on both sides.
  • The C++ kernel guarantees byte-equivalence; the msdecompose OLS path is now also shared (src/headers/olsCore.h), so the only residual divergence is BLAS-dispatch nuances inside Armadillo / LAPACK on the FD-Hessian pass.

When you want reproducible-within-language simulations:

  • Use the language's native seed argument. The output is deterministic given the same seed within that language but doesn't translate across.

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, test_sim_es_r_parity.py, test_sim_ssarima_r_parity.py, test_sim_sma_r_parity.py, test_sim_oes_r_parity.py, test_sim_ces_r_parity.py, test_sim_gum_r_parity.py, test_om_simulate_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 (at the same B) 1e-3 1e-6
ADAM optimal / two-stage confint bounds (each side fits) 2e-2 1e-2
multicov(type="analytical") 1e-10 1e-12
multicov(type="empirical") 1e-10 1e-12
coefbootstrap (stochastic — distributional only) structural
sim_* (plug-in-numbers trick) — 1e-9 to 1e-10
OMG.simulate (plug-in-numbers trick) — 1e-9
OM.simulate (plug-in-numbers trick) xfail (3 cases — see Simulation section above)

Clone this wiki locally