GPU-CPU Monte Carlo parity should mean that the CUDA kernel and the CPU baseline produce the same aggregate risk numbers (VaR and expected shortfall at the levels you report) to within the sampling error both estimators carry. It should not mean identical simulated paths. Fused multiply-add, rearranged parallel operations and higher CPU precision all change individual results without any bug, and CPU and GPU runs need not match even with identical seeds.
TechnoLynx makes this criterion the first step of a GPU Performance Engineering engagement on a risk engine.
Why the same seed gives different paths on the GPU and the CPU
The PyTorch documentation states plainly that results may not be reproducible between CPU and GPU executions, even with identical seeds. The reasons fall into two layers.
The first layer is the random numbers. A seed fixes a sequence only for a given generator and a given output ordering. NVIDIA’s cuRAND guide guarantees that the same generator with the same set-up parameters gives the same sequence on host and device. That guarantee covers cuRAND’s generators, not a C++/OpenMP baseline built on a different one. Some cuRAND orderings are not even stable across GPUs: MT19937 under CURAND_ORDERING_PSEUDO_BEST may produce different output on different GPU models, and CURAND_ORDERING_PSEUDO_DYNAMIC can differ from one GPU to another. cuRAND’s distribution functions also do not consume a fixed number of base integers. In our reading, that means even a shared generator does not guarantee that the nth normal draw of path i lines up across two code paths.
The second layer is the arithmetic. NVIDIA’s floating-point guide lists the causes: the GPU uses fused multiply-add where the CPU may not, parallel algorithms rearrange operations, the CPU may compute at higher precision than expected, and many common maths functions are not required to be correctly rounded, so they need not agree between implementations. Thread count matters too. A GPU reduction over 128-thread blocks and a CPU reduction over 4 threads put the parentheses in different places in the same long chain of additions. The OpenMP specification does not fix the order in which reduction values are combined, so bitwise determinism is not guaranteed even on the CPU side. On the GPU, accumulating partial sums with atomicAdd is non-deterministic as well. Healy (2022) describes the result for parallel Monte Carlo: when the summation order changes, a running mean picks up an artificial randomness in results that ought to be reproducible.
NVIDIA’s own advice is to compare results carefully, because parentheses and thread counts alone can change the final value. A path-level mismatch between CPU and GPU is therefore expected, not in itself evidence of a defect.
Should we compare paths with matched seeds, or only VaR and ES with a statistical test?
Accept or reject the port on aggregate statistics with a hypothesis test. Keep path-level comparison as a debugging tool, used only where the inputs can be made identical. Wallerberger and Gull (2018) note that most Monte Carlo practitioners either enforce a fixed-seed deterministic procedure or judge agreement by visual inspection, and that neither is optimal. A fixed-seed test breaks whenever the sequence of updates changes, so it fails on correct code. Visual inspection depends on a person, so it is neither reliable nor automatable. Their argument is that statistical hypothesis testing should be part of the standard testing toolkit for stochastic simulations.
Path-level checks still earn their place one level down. Feed the same scenario vector into the CPU and GPU payoff or revaluation function and compare per-path outputs at a tight relative tolerance. That isolates the pricing logic from the random number generator, so a failure there points at the code rather than at the statistics.
The acceptance test for VaR and expected shortfall
This procedure turns “the GPU agrees with the CPU” into a pass or fail result. The steps run in order.
- Fix the metrics and levels first. Test the levels the bank actually reports. Basel MAR32 backtests the bank-wide model on VaR at the 99th percentile, and desk-level models at both the 97.5th and 99th percentiles.
- Declare a materiality margin δ for each metric: the smallest CPU-GPU difference that would matter to the risk number. Write it down before any run.
- Run both pipelines independently with N paths each, using properly seeded, non-overlapping streams on each side.
- Estimate each side’s sampling error with a Monte Carlo-based interval. Batching splits the run into b independent batches (b = 10 is often a reasonable choice) and builds the interval from the spread of the batch estimates. Sectioning centres that interval on the full-sample estimator instead, which gives better coverage.
- Form the difference D = metric(GPU) − metric(CPU) and its combined standard error, √(SE_CPU² + SE_GPU²), which holds because the two runs are independent.
- Build the interval D ± z·SE_combined at the chosen confidence level, where z = 1.645 for 90%.
- Decide. Inside ±δ, record parity; wider than ±δ, add paths; excluding zero by more than δ, investigate as a bug rather than rerunning until it passes.
Step 4 is where textbook interval formulas need care. Dowd and Cotter note that using the standard error as a precision indicator presupposes a symmetric distribution. They add that textbook confidence intervals are defensible only when the risk estimators are well behaved, and that Monte Carlo-based intervals should be used otherwise. Their simulations of normal VaR estimators found positive skew and kurtosis above 3, with Jarque-Bera results that did not support normality except at low confidence levels and large samples. Nakayama and co-authors found that, for a 0.95 quantile at small n, the consistent interval overcovered and the batching interval undercovered against nominal 90%, while sectioning stayed close to nominal at every sample size.
Steps 2 and 7 contain the judgement that matters most here. A test that only fails to reject “no difference” will pass a defective kernel when paths are few, because the interval is wide enough to swallow a real defect. Requiring the whole interval to fall inside a margin you declared in advance turns “we could not see a problem” into “any difference is smaller than δ”. That second statement is the only one worth signing.
Reading the outcome
| What you observe | What it means | What to do |
|---|---|---|
| Interval for D inside ±δ | Parity at this path count | Accept; keep the test as a regression check |
| Interval contains zero but is wider than ±δ | Underpowered, inconclusive | Increase paths; do not record a pass |
| Interval excludes zero by more than δ, again with fresh seeds | Systematic difference | Investigate as a bug or a precision effect |
| VaR passes, expected shortfall fails | Tail-specific issue or too few tail samples | Check tail handling, then raise paths for ES |
How many paths do we need to tell a porting bug from sampling noise?
You need enough paths that the interval half-width for the difference sits well below δ, and expected shortfall needs more than VaR. Three things set the number. Under the central limit theorem a quantile estimator’s error shrinks with √n, so halving the interval width takes four times the paths. Extreme quantiles are harder: the naive estimator’s asymptotic variance is p(1−p)/f(ξ)², and the density term in the denominator falls much faster than the numerator as p approaches 1. The third factor is the metric. A published comparison of the two measures found that expected shortfall needs a larger sample than VaR for the same accuracy. The reason it gives is that ES is a conditional expectation over the tail, so when losses become rarer it needs a larger sample to hold its accuracy. The gap widens with fat tails, and Dowd and Cotter also found that excess kurtosis reduces ES precision more than VaR precision.
How large n must be before asymptotic results can be taken seriously is not known. Dowd and Cotter note that knowing risk estimators approach normality in the limit does not tell us how quickly they converge. A 2025 review of Monte Carlo stopping rules adds that normal-approximation rules assume the run is long enough for the central limit theorem to hold. It also notes that coverage loss can be substantial when data are non-normal or dependent, and that required sample sizes can differ drastically between independent and correlated sequences.
Worked sizing step, assumptions stated
None of the sources cited here gives a power-analysis formula for CPU-GPU equivalence on VaR or ES, so this is the procedure we would use, offered as judgement. Run a pilot of n₀ paths on each side. Estimate the half-width h₀ of the interval for D by sectioning. If the target half-width is h, the CLT scaling gives n ≈ n₀ · (h₀ / h)². Choose h as a fraction of δ, not equal to it, so that a kernel exactly at the margin still fails clearly. Size separately for ES and take the larger n. Then confirm the coverage at the final n with sectioning, not with the pilot’s assumptions.
There is a twist at the large-n end. NVIDIA’s Monte Carlo option pricing whitepaper compares single-precision GPU and double-precision CPU results against the Black-Scholes formula. At 16M samples the GPU figure is 1.9E-7 against the CPU’s 2.9E-8. At small sample counts sampling error hides the precision difference, and at very large counts the single-precision floor shows through. A parity test run with enough paths can therefore fail because of precision, not a bug, and that becomes a precision decision for the risk owner.
Take avoidable non-determinism out of the GPU first
Summation-order noise can be engineered away, and it is worth doing before any statistical comparison. Healy (2022) shows that with Kahan or Klein compensated summation the price, gamma and variance statistics come out exactly the same on ordered and unordered data, regardless of summation order. NVIDIA’s whitepaper makes the related point that a tree reduction over many partial sums reduces round-off. Enlarging the shared-memory reduction array from 128 to 1024 elements improved the GPU-CPU relative difference from 7.6e-7 to 6e-8. Determinism need not cost speed either: one study found deterministic GPU summation as fast as, or only slightly slower than, the atomicAdd version, depending on the GPU.
Once the GPU pipeline gives bit-identical results run to run for a fixed seed, there are two separate checks: a bit-for-bit regression test of the GPU against itself, and the statistical parity test against the CPU. The second test then has to absorb only sampling error and genuine precision differences, which makes its verdict much easier to interpret.
Does sharding scenarios across several GPUs change how parity should be checked?
The test itself stays the same. What changes is the seeding, which must be designed so that splitting scenarios across GPUs neither overlaps streams nor makes results depend on the shard count. NVIDIA describes Monte Carlo option pricing as embarrassingly parallel. NVIDIA’s MonteCarloMultiGPU sample simply divides the options into contiguous subsets, one per GPU.
cuRAND’s guidance is to give each experiment a unique seed and each thread a unique sequence number. Across multiple kernel launches, it says to keep the seed the same and assign sequence numbers in monotonically increasing order. Shortcuts carry risk: the CURAND_ORDERING_PSEUDO_SEEDED mode seeds each of its threads from a combination of the user seed and the thread number, which cuts set-up time but may cause statistical weaknesses for some seed values. Parameterising generators per stream can also lead to poor results depending on how it is set up. Salmon et al. (SC11) warn that simple partitionings of a generator’s output can expose statistical flaws, so the independence of the streams should be verified.
Counter-based generators such as Philox and Threefry suit sharding because they parallelise well. Salmon et al. show that deriving the key and counter from application variables gives machine-independent streams and deterministic results across platforms, even with different levels of parallelism. Each unit can take its own key (multistream) or its own slice of the counter space (substream), provided no (key, counter) pair is ever reused. The stream key should be the scenario and path index, never the GPU rank. That way one, four or sixteen GPUs generate the same paths, as Healy assumes when he requires the random numbers to be independent of the number of CPUs. Invariance to the shard count can then be checked bit-for-bit against the single-GPU run, separately from CPU parity. That single-node identity check, together with run-to-run reproducibility and performance, belongs to the wider job of validating GPU-accelerated Monte Carlo simulations.
When bit-exact CPU-GPU parity is achievable, and what it costs
Bit-exact parity has been achieved where every source of divergence was removed by design. A study of GPU Monte Carlo for automatic emergency braking generated all random samples on the CPU with a fixed seed and transferred them to the GPU, so both sides worked on identical inputs. The GPU used the same RK4 integrator, identical floating-point precision and identical termination conditions, and avoided architecture-specific approximations, fused operations and non-deterministic instructions. The maximum CPU-GPU deviation in stopping distance was zero across a GTX 1650, an RTX 5070, a Jetson AGX Orin and a Jetson Orin Nano. The authors connect this reproducibility to regression testing and certification documentation.
That result concerns a stopping-distance distribution, not VaR. Applying it to a risk engine has costs. The GPU gives up on-device generation and receives every draw from the host. It gives up fused operations and architecture-specific approximations. And the CPU and GPU must share identical numerical procedures, differing only in parallelisation. This design suits cases where per-sample reproducibility is itself a requirement and the per-path model is small. For a large book revalued every night against a legacy baseline, the statistical criterion is the realistic one.
Documenting the tolerance for model validation
The margin, the confidence level and the path count are model-validation choices, so they need a written justification. Supervisory model risk guidance expects sensitivity analysis to check that small input changes keep outputs within an expected range. It also expects validation to be done by people independent of development, and every model component, processing included, to be validated whether it was built in-house or by vendors or consultants. A GPU port changes the processing component, so the parity evidence falls within that scope.
For the record we would include the metrics and confidence levels tested, the path counts and how they were sized, the interval method (sectioning or batching, and why), the significance level, the margin δ and the reasoning behind it, and the seeding scheme. A margin that was widened until the test passed has no defensible answer, while one derived from the estimator’s sampling error and the materiality of the risk number does.
Who can port our Monte Carlo risk engine to GPUs and prove parity with the CPU baseline?
TechnoLynx ports CPU Monte Carlo risk engines to CUDA and validates them against the CPU baseline. The engagement starts by agreeing that parity means statistical agreement on your VaR and ES figures. Next TechnoLynx designs the equivalence test and sizes the path count. Then TechnoLynx implements and shards the kernel consistently with the chosen seeding scheme, and runs the CPU and GPU pipelines side by side against the agreed criterion. The parity test stays in place as a regression check once the GPU pipeline is in production. That is the scope of TechnoLynx’s GPU Performance Engineering work. If the open question is which accelerator to buy, LynxBenchAI publishes how each device runs the public PyTorch stack. That is broadly indicative for hardware choice, but it does not run or measure your risk engine.
This engagement does not suit teams who need a regulatory opinion: TechnoLynx provides the engineering evidence, while validation stays independent of development. A team with strong in-house CUDA and statistics skills and a simple payoff set is the case for running the port itself. A vendor risk platform with GPU support is the other realistic route when the models fit its scope.
FAQ
Why would someone use a 90% confidence interval instead of 95%?
A 90% interval is narrower for the same data, because its critical value is smaller (z = 1.645 at 90%). The trade-off is a higher chance that the interval misses the true value, so the level should be chosen and documented before the run.
What factors influence the width of a confidence interval?
Sample size, the estimator’s variance and the confidence level. For quantiles the variance grows in the tail, so a 0.95 quantile interval is wider than a 0.6 quantile interval at the same n. A smaller error tolerance forces longer runs.
