Two Backends, One Posterior

Four hours and twenty-one minutes after we launched it, the verification run on mac-247 finished. Two thousand five hundred draws under each of two backends from the same regime trading model, at the same seed, on a 2013 GeForce GT 650M paired with a FreeBSD 15 kernel. The output was a report of eight parameters and one line at the bottom. The line said:

{:fail, [
  {"logit_w1", %{rel_mean: 0.8046, sd_ratio: 1.611, ...}},
  {"theta_mr", %{rel_mean: 0.1472, sd_ratio: 8.296, ...}}
]}

A trading fleet ships when its posterior distributions agree across the machines that will compute them. The CPU BinaryBackend path and the GPU Vulkan chain shader had just failed to agree on two of eight random variables. The logit_w1 mean had shifted by four fifths of a posterior standard deviation. The theta_mr posterior sd was 8.3 times wider on the GPU than on the CPU. The chain shader, painstakingly synthesised from the same intermediate representation, dispatched through the same Vulkano-based NIF, running on hardware whose bring-up we had been through more than once, was producing measurably different samples.

This was not the first failure. The seed had been changed, the sample count had been raised from a hundred to five thousand, the compiler had been swapped between Nx.Vulkan and Nx.Defn.Evaluator, and the disagreement had survived. What made the failure interesting was that the shader was correct. What made it hard to prove was that the test was correct too. What made it publishable is that neither of these assertions was actually true.

What we thought we were testing

The Exmc.Trading.RegimeModel is a three-regime Bayesian model of asset returns. Eight free parameters: two regime-weight logits, a mean-reversion mean and speed, three scale parameters governing trend / mean-reversion / volatility regimes, and a trend drift. Sampled with NUTS through the eXMC sampler under one of two compilers. The compiler is a runtime choice.

Under compiler: :none, the leapfrog integrator runs through Nx.Defn.Evaluator, dispatching each Nx primitive to Nx.BinaryBackend— f64 arithmetic on the CPU. Under compiler: :vulkan, the sampler compiles a fused f64 “chain shader” whose GLSL was synthesised from the same IR that the CPU path executes. At K=32 that shader runs 32 leapfrog steps per vkQueueSubmit, which is the point.

For the trial fleet to work, both paths have to sample the same posterior. This is a stronger statement than it looks. Two implementations of the same Markov chain, seeded identically, produce different samples if they compute their leapfrog updates in any different order — and floating-point arithmetic is not associative. What we needed was a test that could detect “the chain shader is wrong” while tolerating “different arithmetic order.”

The wrong tolerance, first draft

The existing chain-shader coverage test compared 100-sample posteriors from CPU and GPU using two tolerances borrowed from the posteriordb validator:

rel_mean_tol: 0.5       # |mean_a - mean_b| / sd_a ≤ 0.5
sd_ratio ∈ [0.5, 2.0]   # sd_b/sd_a in the fifty percent band

On the regime model at seed=42, GPU posterior sd for theta_mr was 2.7 percent of CPU posterior sd. The ratio was 0.027, an order of magnitude outside the accepted band. The test failed loudly. The obvious hypothesis, and the one that shaped the first day of investigation, was that the GPU chain shader was underestimating something — possibly the tails of the HalfCauchy(1) prior used for the mean-reversion speed theta_mr. GLSL's log-standard instruction operates on floats, not doubles. The chain shader's “f64” log helper was double(log(float(x))) — a downcast to f32, apply the intrinsic, upcast back. Twenty-four bits of precision lost on every logarithm. This is what any careful reviewer sees first, in lib/exmc/nuts/custom_synth/multi_rv_custom_spec.ex, around line 404.

It was the wrong culprit.

Five seeds, three hosts, one refutation

The refutation came from what became “Phase 1” of the investigation: run the same 100-sample regime posterior five times at seeds 1 through 5, on each of three hosts: mac-247 (GT 650M), mac-248 (GT 750M), and super-io (Ampere RTX 3060 Ti). Record sd_ratio for the three heavy-tailed sigma parameters. The prediction: if the log-downcast was causing the collapse, Ampere — whose f32 pipeline is more accurate than the two Keplers' — should show tighter agreement.

Two things fell out of this run. The first was that mac-247 and mac-248 produced bit-identical ratios across all five seeds. Same GT 650M posterior samples as the GT 750M, byte for byte, seed after seed. The chain shader was deterministic across two generations of Kepler hardware.

The second was that super-io's Ampere ratios were worse than either Kepler's. On sigma_mr the Ampere sd_ratio ranged from 0.013 to 186 across the same five seeds — four orders of magnitude in a small sample. If the f32 log-downcast had been causing a Kepler-specific collapse, Ampere would have shown tighter agreement, not wider. It didn't. The chain shader's precision was not the culprit. The failing tolerance had been fitted to the wrong noise distribution.

HalfCauchy(1) is a heavy-tailed distribution with no finite variance. Its sample standard deviation, computed from any finite number of draws, is dominated by whether the chain happened to visit the far tail. At 100 samples the estimate of posterior sd for a HalfCauchy-prior parameter is essentially noise. The rel_mean_tol: 0.5, sd_ratio: [0.5, 2.0] tolerance was a good default for posteriors with bounded moments. For HalfCauchy tails at N=100 it was a coin flip.

This settled the Kepler-precision hypothesis. It also settled nothing about F1's failure at seed=42 and N=5000, which had been reproducible, and which still needed an explanation.

The shader, alone, is correct

The next question we asked the machinery was whether the chain shader math actually agreed with the CPU integrator when we compared them in isolation. One leapfrog step, matched initial state, no sampler around it. If the shader was drifting from CPU, the answer was there.

|Δq|_∞    = 2.24e-8    # position, essentially identical
|Δp|_∞    = 5.16e-6    # momentum, small drift
|Δgrad|_∞ = 2.24e-4    # gradient, ~1.5e-6 relative
|Δlogp|      = 323          # convention difference, not error

On Ampere and on the GT 750M, byte for byte the same. The 2.24e-4 max-abs difference on the gradient, whose largest element had magnitude 150, corresponds to a relative error of 1.5e-6 — consistent with the f32-cast on log and exp, completely inconsistent with the 0.8-standard-deviation gap in logit_w1 that F1 had reported after 5000 samples. Cumulative shader precision loss at 1.5e-6 relative per step cannot compound to 0.8 sigmas over any trajectory a real NUTS sampler will visit. It can't compound to that even if you tried. The 323 gap on logp was a convention difference: the chain shader outputs logp[k] as the log-density at the start of step k, not the end.

This left one class of explanation for F1: the sampler around the chain shader, not the chain shader itself. Something in the tree-building, mass-matrix adaptation, or step-size dual averaging was accumulating different arithmetic between the two backends.

The warmup that agreed with itself but not with the other

The instrumentation we added next was a trajectory-diff harness. Same seed, same fixture, 200 warmup steps, 50 sample steps under each compiler. Compare per-iteration draws.

End of warmup on CPU: eps = 0.283, one component of inv_mass at 3.95. End of warmup on GPU under Vulkan: eps = 0.239, same component at 1.33. Every one of the 50 sample iterations produced different draws. Not slightly different, in some cases. |Δq|_∞ reached 484 by iteration 39.

Reasonable next hypothesis: warmup routes differently under the two compilers. The D91 Option C fix, shipped a day earlier to address a K=32 warmup step-size collapse, forces warmup onto the per-leaf host-fallback path when the compiler is Nx.Vulkan. Under compiler: :none, the gate doesn't fire and warmup uses the K-fused speculative pre-compute path. Two different tree geometries during warmup, producing two different mass-matrix estimates, producing two different step sizes.

We turned off speculative pre-computation globally and re-ran, forcing both backends onto per-leaf warmup:

The GPU eps hadn't moved. Only the CPU changed. Whatever was happening in the CPU warmup path had responded to the speculative-pre-compute flag; the GPU had been unaffected because its warmup path was already the per-leaf one, thanks to D91.

And yet the divergence persisted. Warmup was not the issue. Sampling-phase K-fusion was not the issue. The two backends were producing different arithmetic at some layer we hadn't looked at yet.

The design decision that produces bit-drift

Exmc.Compiler.build_step_fn/3 is where the leapfrog integrator becomes an executable function. It calls Exmc.JIT.jit(fn, jit_opts), whose return value under compiler: :none is a function that executes through Nx.Defn.Evaluator against Nx.BinaryBackend, and under compiler: :vulkan is a function that executes through Nx.Vulkan's JIT (which host-falls-back to BinaryBackend for unsupported ops). Both paths ultimately do their arithmetic in the same f64 BinaryBackend. The paths differ in the intervening structure: the Evaluator's op-by-op interpretation versus Nx.Vulkan's JIT tracing and potential fusion.

Floating-point addition is commutative but not associative. A JIT that reorders a sum-of-products into a different tree produces different bit patterns than an interpreter that processes them left to right. Both are correct. Both are producing the requested computation to within the same guaranteed precision. Neither is the reference implementation of the other. This is not a bug. It is the reason the compiler exists.

The F1 divergence at N=5000 seed=42 was compiler-associativity drift, amplified by the sampler's sensitivity to initial state. Two backends, one seed. From iteration zero forward they diverge, not because either is wrong, but because they are arithmetically distinct implementations of the same mathematical operation. The chain shader we had suspected was correct. The tolerance we had used to check it was measuring floating-point non-associativity and calling it precision loss.

A verification standard for probabilistic backends

This is where the investigation stops being about a chain shader and starts being about how anyone should test that a probabilistic-programming backend samples the same posterior as a reference implementation.

Bit-identical comparison across backends is not achievable in this configuration and probably not achievable in most similar configurations. Once the compiler is allowed to reorder operations, the sample sequences will differ from iteration zero. This is by design. Any test that treats it as a defect will fail without illuminating anything.

Statistical parity is the correct contract. The test should not be “does implementation A produce the same samples as implementation B at the same seed.” It should be “do implementations A and B produce statistically- equivalent posterior distributions,” measured under tolerances sized to Monte Carlo error and augmented by whatever per-backend numerical drift exists in the arithmetic reordering.

Robust statistics beat moments for heavy-tailed posteriors. The first statistical-parity attempt used posterior mean and standard deviation with tolerance |Δmean|/ref_sd < 0.2 and sd_ratio ∈ [0.5, 2.0]. Three of the eight regime parameters failed. All three were parameters with HalfCauchy(1) priors — distributions whose population variance is undefined and whose sample sd is dominated by tail luck at any finite sample size. On theta_mr the CPU sd came out to 577; the GPU sd came out to 29; the two chains had happened to visit different tail excursions. Neither number was wrong. Neither number was converging to anything. The tolerance was measuring sample-variance-of-a-Cauchy, which is defined but not particularly informative.

The second attempt replaced mean and sd with median and interquartile range: |Δmedian|/ref_iqr < 0.5 and iqr_ratio ∈ [0.5, 2.0]. The sample median of a HalfCauchy is well-defined, converges at rate 1/√N, and is insensitive to tail excursions. IQR is the same story for scale. On mac-248 the report line came out to:

=== VERDICT ===
PASS — all 8 RVs match statistically on median+IQR.
Stage 2 prereq closes with the robust-stats contract.

theta_mr, the parameter whose mean/sd comparison had been failing at a factor of twenty, produced CPU median 1.09 with IQR 2.39 and GPU median 0.86 with IQR 2.25. Ratio 0.94. On mac-247, on entirely different Kepler silicon, the numbers were byte-identical to mac-248. Cross-Kepler bit-determinism, which Phase 1 had shown at 100 samples, held at 2000. Cross-backend statistical parity, which four days of investigation had implied should exist and which no test we'd written had been able to demonstrate, appeared as soon as we asked the correct question.

What the standard looks like

The verification-as-a-standard idea that emerged is a testing convention we would like to see in Nx.Backend-implementing libraries generally. The components:

N chains, not N samples. Draw N independent Markov chains under each backend, from different seeds. Compute per-backend posterior statistics on the concatenation. This turns per-chain variance into signal instead of noise.

Robust statistics for heavy-tailed priors. Median and IQR — or, more generally, quantiles at 25 / 50 / 75 — are the right summary for any parameter whose prior is HalfCauchy, Cauchy, Student-T with small degrees of freedom, LogNormal with heavy scale, or otherwise infinite-moment. For well-behaved parameters mean and sd are fine.

Tolerance sized to what matters. For a trading system the relevant tolerance is signal quality. Posterior means of trend and volatility parameters must agree to a small fraction of the trading-signal threshold. Heavy-tailed nuisance parameters can drift more, because they enter the signal only through their posterior median. For a scientific application the tolerance is whatever your downstream inference needs, and no smaller.

Cross-hardware determinism as a separate axis. The Kepler-Kepler byte-identity that fell out of the verification is a valuable property in its own right. It says “a reproducible run on our hardware will reproduce on your hardware.” That is not the same as cross-backend statistical parity, and it should be measured independently. In practice, an implementation that produces byte-identical output on two GPUs of the same architectural family, and statistically-equivalent output against a CPU reference on a third, has established both properties.

Chain-shader math verified in isolation. The one-step-leapfrog probe we built as part of the investigation is the diagnostic to keep. When a chain shader diverges from a CPU reference at the sampler level, running one step of both against a fixed initial state answers whether the shader math is at fault. In our case it wasn't. In another case, it might be, and the probe will find it in a few seconds instead of the several hours a full sampling run takes.

What was actually happening on the GT 650M

The regime model produces the same posterior distribution on FreeBSD 15 running vulkano against a 2013 GeForce GT 650M as it does on the CPU BinaryBackend reference. It also produces the same posterior on the 2013 GeForce GT 750M in a different 2013-vintage Mac Pro. The two Keplers agree byte for byte across seeds. The Ampere agrees statistically. The one thing they don't do is match each other bit-for-bit through the compiler, and they cannot, because Nx.Vulkan's JIT is not required to reorder operations the way Nx.Defn.Evaluator does.

Stage 2 of the D88 FreeBSD-fleet plan, which requires regime-model posterior parity between mac-248 and the super-io reference, closes for both Kepler hosts under the robust-stats contract. The chain shader was correct throughout. The tolerance we applied to it was wrong, then wrong differently, then finally right. The passage from a mean/sd failure with a factor-of-twenty gap to a median/IQR pass with a ratio of 0.94 happened without changing a line of shader code, or a line of sampler code, or a line of NIF code. What changed was our answer to the question: what do we owe the reference implementation.

Not bit-identity. That's a fantasy the compilers refuse to maintain. Statistical equivalence, with the right statistics, with the right tolerance, with the right null hypothesis.

Coda

Four hours and twenty-one minutes was the wall time of the F1 failure that started the investigation. Fifty minutes was the wall time of each of the median/IQR runs that closed it. The efficient path would have been to skip the intermediate steps and land directly on statistical parity with robust statistics. The efficient path was not available in advance. What was available was a shader that seemed to disagree with itself, a tolerance that seemed reasonable, and two Kepler GPUs from 2013 that had, over the course of a very long day, disagreed identically at every level of the investigation.

The verification standard survives because the shader passed. The shader passes because it computes the same posterior. It computes the same posterior because the leapfrog it implements matches the leapfrog the CPU implements, one f64 op at a time, byte for byte on the same silicon architecture and statistically across architectures. The last thing this investigation needed was a wrong answer at the shader level. That would have been another kind of story. This one is a story about tests: the ones we wrote first, the ones we should have written, the ones we propose that everyone write from now on.

The investigation notes live in exmc/research/, under filenames beginning kepler_precision_p1, kepler_shader_1step, kepler_trajectory_diff, and stage2_prereq_median_iqr. The USDT probes live in lib/exmc/nuts/leapfrog.ex and lib/exmc/nuts/vulkan/dispatch.ex; the DTrace consumer that pairs them is at trial/dtrace/leapfrog_compare.d. The verification proposal above is a discussion document, not a merged convention. It applies wherever an Nx.Backend-implementing library needs to check that its sampler agrees with a reference. If you write one, tell us.