RNG parity for stochastic methods — feasibility analysis

Question (from maintainer): make the stochastic methods reproduce their external references exactly (Jaccard 1.0 for selectors). Do it properly: make the RNG an optional model parameter, implement additional RNG algorithms beyond splitmix64 to cover the tests, and initialise the tests with the right RNG + seed.

This document analyses impact, feasibility, and effort, and proposes a tiered plan. It is the basis for a Codex review.

Status update (2026-05-30)

This analysis is retained for rationale, but it is no longer the current task list. The additive RNG infrastructure has shipped (n4m_rng_kind_t with splitmix64, PCG64, R-MT, and numpy RandomState-MT), the UVE R-exact pilot is available behind rng_kind=MT_R, and the dashboard now renders documented selector RNG/noise/model mismatches as cross_check/BD J instead of red numeric divergences. See RNG_TIER0_INVENTORY.md for the authoritative post-investigation verdict table.

1. What we have today (verified in source)

• Two RNGs already exist in the core:

  • splitmix64 — drives binary-mask / index generators in the selectors (hardcoded inside each method).

  • n4m_rng_pcg64_* (public ABI, cpp/include/n4m/n4m.h:367+) — documented to reproduce numpy.random.default_rng(seed) bit-for-bit (uint64_fill ≡ default_rng(s).integers, standard_normal_fill ≡ ...).

• The config has NO rng-kind selector. cfg exposes solver, deflation, scaling, etc., but no n4m_config_set_rng_kind. RNG choice is currently baked into each method’s C++.

• ~14 seed-taking parameters across the public PLS/selector ABI (uint64_t seed | noise_seed | randomization_seed in pls.h): the stochastic surface = the variable selectors (uve, emcuve, cars, random_frog, spa-ish, ga, pso, vissa, irf, randomization, stability, scars, …) plus the ensemble PLS (bagging, boosting, random_subspace) and some augmenters.

2. The real root cause (uve_select, worked example)

The dashboard “Jaccard 0.75” for uve_select is NOT a binding bug:

Tier	What it actually runs	Selected indices
registry / cpp / python_tier1	call R plsVarSel::mcuve_pls (the canonical path is legacy=False → R)	[2,5,8,13]
ref_r_plsvarsel	R mcuve_pls (oracle)	[2,5,8,13]
r_tier1 / matlab_tier1 / python_tier2	the n4m C kernel n4m_uve_select	[2,5,8]

So the comparison is n4m’s own C kernel vs R mcuve_pls, and they differ by one boundary variable (index 13). Per the registry docstring (_uve_select_pls4all), the C kernel deliberately differs from R in three independent ways:

1. Noise model — n4m adds a fixed noise_features count of signed-uniform noise standardised per column; R adds ncol(X) uniform-noise columns.

2. RNG — n4m splitmix64 vs R Mersenne-Twister (set.seed(11)).

3. Fold / CV semantics — contiguous-fold coefficient sampling vs R’s.

Consequence: matching R bit-for-bit needs ALL THREE aligned, not just the RNG. With a different noise model or draw order, even an identical RNG stream diverges. “Make the RNG a parameter” is necessary but far from sufficient.

This generalises: each stochastic method’s reference is one of

• R (Mersenne-Twister, R’s set.seed scrambler, R’s rnorm inversion) — most selectors (plsVarSel, mdatools),

• numpy/sklearn (PCG64 — which n4m ALREADY reproduces) — the python-composite methods,

• auswahl (numpy under the hood) — random_frog / vissa / irf / vip_spa,

• deterministic (no RNG) — spa, vip, coefficient, interval.

3. Why exact repro is hard (the honest part)

Bit-exact match to R requires, per method:

1. R-compatible Mersenne-Twister. R’s RNG is MT19937 but with R’s own set.seed scrambling (Init_R_seed) and R’s unif_rand/norm_rand (inversion via qnorm). This is NOT the same byte stream as numpy’s RandomState (also MT but seeded differently) nor generic MT. It is well-documented (R src/main/RNG.c) but fiddly — ~a few hundred LOC + an exact-stream test against R.

2. Exact noise model — replicate how R builds its noise columns (count, distribution, placement).

3. Exact draw order — the sequence of RNG consumption must match R’s loop structure draw-for-draw; one extra/missing draw desynchronises everything.

4. Exact CV / fold semantics and the exact thresholding rule.

(2)–(4) amount to porting each R function’s algorithm into C++. That is the bulk of the effort and the risky part — some references (auswahl, mdatools) have undocumented internals.

There is also a design tension: n4m’s C kernels were deliberately built with their OWN deterministic, documented noise model (so they are reproducible and not hostage to R’s RNG). Making them mimic R means abandoning n4m’s own design to copy R bit-for-bit. That is a legitimate choice for “dashboard shows green”, but it is “n4m becomes R”, not “n4m is independently correct”.

3b. Concrete classification (verified by grep over registry source)

• 13 seed-taking methods in the public ABI (uint64_t seed|noise_seed| randomization_seed in pls.h).

• 15 registry helper functions route a selector to R (_*_indices_via_r): bve, cars(enpls), emcuve, ga, ipls/interval, ipw, randomization, rep, shaving, spa, sr, stability/mcuve, st, uve, wvc_threshold.

• 3 route to auswahl (numpy): random_frog, vip_spa, irf.

• 12 methods have their DEFAULT path (legacy=False) routing to R — i.e. for these 12 the “canonical n4m” value shown on the dashboard is R’s output, and the n4m C kernel is the legacy=True opt-in that the binding tiers exercise. This is the set where C-kernel-vs-R parity is measured.

• The numpy/sklearn composites (bagging/boosting/random_subspace/n_pls/ pls_qda/pls_glm/pls_cox/pls_logistic/group_sparse) BYPASS the C kernel entirely — NOT an RNG problem; correctly not_available for C-kernel bindings.

Net (corrected): the RNG-exact-repro effort is ~12–15 R-referenced selectors, not the 6 I first estimated — i.e. Tier C is larger than my first pass. The “cheap PCG64 win” (Tier A) is genuinely small, because the numpy-referenced stochastic surface is mostly composites the C kernel does not run. Tier B (RNG abstraction) is the enabler for the 12–15 R cases + 3 auswahl cases.

3c. Codex review corrections (verified vs source — authoritative over §3b)

A read-only Codex review (high reasoning) checked the analysis against the code and corrected three factual points + added two risks:

1. R-referenced set ≥ 12 and iriv_select is NOT R-backed. R-default selectors include stability, uve, spa, cars, ga, shaving, bve, emcuve, randomization, rep, ipw, st, wvc_threshold (registry.py:3588/3845/3957/4867/ 5001/5119/5244/5381). iriv_select routes through a NumPy port using default_rng (registry.py:4485, 7659), not R.

2. auswahl is NOT numpy-PCG64. It passes random_state=int(seed) (registry.py:1562/4014/7543); auswahl 0.9.0 uses sklearn check_random_state → legacy numpy.random.RandomState (MT19937), not default_rng/PCG64. The auswahl cases need a THIRD engine (numpy RandomState-MT); the “cheap PCG64 Tier A” is smaller/different than hoped.

3. n4m.h:344 overclaims PCG64 adoption — UVE/EMCUVE public APIs take only a seed (pls.h:1404/1605) and UVE uses file-local splitmix64 (uve_selection.cpp:144/215). The header comment needs correcting.

4. RISK (Tier B): symbol-additive ≠ behaviour-preserving. splitmix is duplicated file-local across ~10 selectors with method-specific seed perturbations (EMCUVE noise_seed + ensemble*golden emcuve_selection.cpp:87; randomization per-permutation seed :281; randomized-SVD context-seed offset model.cpp:766). Freeze current splitmix streams with golden tests BEFORE the vtable refactor; default rng_kind must = exact-current-splitmix.

5. RISK (Tier C): R-MT is bounded to set.seed+unif_rand(+rnorm). Methods that use R’s sample(), CV-segment generation or package internals also need sample.kind semantics + package draw order — NOT “a few hundred LOC”. The R-MT engine just landed covers the bounded case only.

Also: some references are “closest installable analog”, not the method proper (CARS vs enpls.fs, registry.py:6061) → keep those as documented cross_check/Jaccard, never a forced bit-match.

Revised verdict (Codex-aligned), supersedes §5:

• Tier 0 first — fix the method inventory; per method decide canonical / compat-mode / documented-cross_check (§3b was unreliable).

• Tier B only after golden stream tests freeze current splitmix; rng_kind additive, default = current splitmix.

• Implement only the RNGs actually needed: PCG64 exists; add numpy RandomState-MT for auswahl; R-MT (DONE) for true R-compat methods.

• Tier C only for high-value methods where the reference IS the target algorithm (UVE/MCUVE); cross_check/Jaccard for the analog-only ones.

4. Proposed tiered plan

Tier B (infrastructure — do this regardless; it is the right abstraction)

Make RNG a first-class, pluggable, configurable concept:

• typedef enum n4m_rng_kind_t { N4M_RNG_SPLITMIX64=0, N4M_RNG_PCG64=1, N4M_RNG_MT19937_R=2, ... } in a public header.

• A core RNG interface (vtable: next_u64, next_double, next_normal, reseed) so methods draw from an abstract RNG, not a hardcoded one.

• n4m_config_set_rng_kind(cfg, kind) + reuse the existing seed.

• Wire each of the ~14 stochastic methods to draw from the configured RNG.

• ABI impact: new public symbols → regenerate expected_symbols_* on all 3 platforms + bump N4M_ABI_VERSION_* + docs/abi/changes_log.md. The enum and setter are additive (backward-compatible), default = current behaviour (splitmix64) so nothing silently changes.

• Effort: ~1–2 weeks C++ + ABI gates + per-binding plumbing (R/MATLAB/Py/JS must expose the new setter). Self-contained, testable, valuable on its own.

Tier A (cheap wins — do first)

For methods whose reference is numpy/sklearn, n4m ALREADY has PCG64. If a method currently uses splitmix64 but its reference uses numpy, switching it to the existing PCG64 RNG (via Tier B’s config, or directly) can yield exact repro cheaply — no new RNG needed. Audit which of the ~14 fall in this bucket.

Tier C (expensive — per-method R alignment; decide method-by-method)

Implement N4M_RNG_MT19937_R (R-compatible MT + rnorm inversion), then for each R-referenced selector replicate R’s noise model + draw order + folds + threshold. Effort: days per method × ~10–12 methods = weeks→months, with no guarantee for the undocumented ones (auswahl/mdatools internals).

5. Recommendation (to be reviewed by Codex)

1. Tier B is worth doing on its own merits — “RNG as a configurable, pluggable model parameter” is the correct architecture and unblocks everything else. Ship it with default = splitmix64 (no behaviour change) and PCG64 already wired.

2. Tier A likely converts the numpy/sklearn-referenced stochastic methods to exact repro quickly once Tier B exists.

3. Tier C is a per-method decision, not a blanket commitment. For each R-referenced selector, weigh “port R’s algorithm into the C kernel” vs “accept it as a documented cross_check (n4m has its own valid algorithm)”. Bit-exact-to-R is achievable but is real reverse-engineering and partly makes n4m a re-implementation of R.

4. Independent of all the above, the dashboard must stop showing a Jaccard overlap as if it were an RMSE divergence (separate task #3, in progress) — a Jaccard 0.75 is “75% feature overlap”, not “0.75 numeric drift”.

6. Open questions for the maintainer

• For the R-referenced selectors, is the goal “n4m reproduces R exactly” (Tier C, n4m mimics R) or “n4m has a documented, correct selection that we display honestly as a cross-check”? This is the cost driver.

• Is the new RNG abstraction (Tier B) acceptable as an additive ABI change with an ABI-version bump and snapshot regen on all platforms?