Validating custom CUDA kernels against cuBLAS or cuSOLVER starts from one fact. Two correct implementations of the same maths can disagree in the low-order bits, because floating-point summation order, fused multiply-add contraction, TF32 tensor-core paths and split-K reductions all change where rounding falls. So “does it match bit for bit?” is a test floating-point code cannot always pass.

Numerical validation belongs inside AI engineering for finance, not after the optimisation work.

Why do custom CUDA kernels give slightly different results from cuBLAS or cuSOLVER?

Floating-point addition and multiplication are not associative: every intermediate result is rounded, and the order of evaluation decides how those rounding errors accumulate. NVIDIA’s CUDA Best Practices Guide puts it directly: (A+B)+C is not guaranteed to equal A+(B+C), and parallelising a computation changes the order of operations, so parallel results need not match sequential ones. NVIDIA’s floating-point white paper shows that seemingly identical computations give different results even when every basic operation is IEEE 754 compliant.

Several distinct mechanisms produce the reordering:

• Reduction layout. A typical GPU reduction has the threads in a block do a parallel tree reduction, followed by a serial reduction across blocks, so changing the threads per block rearranges the parentheses in the sum.
• FMA contraction. A fused multiply-add rounds once where a separate multiply and add round twice. Whether the compiler fuses depends on the platform and the build, which matters whenever CPU and GPU results are compared, and nvcc allows contractions and re-associations by default.
• Tensor cores and TF32. Since version 11.0.0, cuBLAS may use Tensor Cores automatically unless pedantic modes disable them. Tensor cores accumulate in their own precision and order, so a matrix multiply on them does not agree bit for bit with the same multiply built from ordinary FMAs.
• Split-K. Some cuBLAS GEMM algorithms split the K dimension to raise occupancy and then sum the chunks. Because the library may pick different kernels or reduction schemes for different shapes, splitting a problem along M or N can give results that are not bitwise identical to the unsplit problem. Split-K is also why matmuls lose batch invariance. The bits are fixed by where the contracted axis is cut, how long each part is, and the types the partial sums are stored and merged in.
• Host precision. x86 processors can use 80-bit extended precision, which often differs from pure 64-bit device arithmetic.

Add to this that a vendor library does not fully disclose what arithmetic it performed, so a kernel written elsewhere has trouble matching it. NVIDIA’s own guidance follows from all of this: differences between CPU and GPU results are not automatically evidence that the GPU result is wrong.

Before any diff means anything, find out which internal path the library took for your exact call: the math mode and compute type on the cuBLAS handle, and on the cuSOLVER side the factorisation routine and whether a mixed-precision iterative-refinement solver is running with a lower internal precision. TechnoLynx’s numerical provenance audit gives you a documented explanation for each difference: which routine and precision mode your pipeline’s calls actually take, settled before anyone argues about tolerances.

Validating custom CUDA kernels operator by operator

A comparison of final prices alone does not show where divergence entered the pipeline. To locate it, emit intermediate results from both implementations at matching checkpoints (after pivoting, after each reduction, after each solver stage), diff each against a trusted reference, and find the first point where the difference becomes material. The audit should then classify that divergence as systematic bias tied to an algorithmic choice or as rounding noise, and check whether it grows with grid size or time-step count.

TechnoLynx builds this as an operator-level diff harness that runs on the target GPU and precision, sweeps problem sizes to classify bias against noise, and is packaged so it can be rerun whenever a kernel or library version changes.

The procedure, in the order it has to run

1. Fix the reference. Keep known-good outputs from representative inputs and compare every change against them using criteria suited to the algorithm, since bitwise identity is often not achievable with floating point.
2. Pin down what the library did. Runtime profiling shows how many kernels a cuBLAS call launched, over what grid and whether it took a workspace, which narrows a shape to an algorithm family.
3. Probe the black box. The 2026 study that reconstructed cuBLAS GEMM bit for bit settled the arithmetic with numerical experiments: extreme values of opposite sign plus one tiny value placed on an otherwise empty axis show where the rounding groups fall.
4. Use what the library exposes. cuSOLVER’s expert iterative-refinement interface returns the residual norm at each iteration and the iteration count, which gives you a checkpoint inside a routine you otherwise cannot open.
5. Diff element-wise, then by size. Report the maximum absolute difference for every quantity and add relative error for reduction-heavy paths, the paths most sensitive to accumulation order.
6. Classify. Noise from accumulation order grows gradually with accumulation depth; the Feltor developers report that even a tree summation’s error grows with the square root of the array length. A difference with a consistent sign that does not follow that pattern points to an algorithmic difference or a bug.

Two limits apply. A vendor library will not fully disclose the arithmetic it performed, so a kernel written elsewhere has difficulty matching it. And telling bias from noise needs enough varied runs to have statistical power; sparse sampling can misclassify either way.

Can cuBLAS be stopped from using TF32 so it runs full FP64 and matches our reference?

Yes, and the reference run should do exactly that. CUBLAS_PEDANTIC_MATH uses the prescribed precision and standardised arithmetic for every phase of the calculation and is intended for numerical robustness studies, testing and debugging. CUBLAS_COMPUTE_64F_PEDANTIC uses FP64 throughout and also disables algorithmic optimisations such as Gauss complexity reduction (3M). Setting the environment variable NVIDIA_TF32_OVERRIDE=0 overrides any default or programmatic configuration, so cuBLAS will not accelerate single-precision work with TF32. For mixed-precision routines, CUBLAS_MATH_DISALLOW_REDUCED_PRECISION_REDUCTION forces reductions to use the compute type rather than the lower-precision output type.

One more trap: on recent toolkits, “FP64” GEMM is not always native FP64. cuBLAS in CUDA 13.0 Update 2 added FP64 emulation on Blackwell INT8 tensor cores, with an automatic framework that picks emulation or native execution. For a reference run, request the pedantic FP64 compute type explicitly rather than trusting the default path. On the cuSOLVER side, math modes are set per handle, and an iterative-refinement solve whose lowest precision equals the input/output type falls back to a main-precision factorisation. On the compiler side, NVIDIA recommends the default -ftz=false -prec-div=true -prec-sqrt=true and no --use_fast_math for the most portable results.

NVIDIA’s cuSOLVER documentation expects at most a 2x speed-up for large problems when an FP64 solve may factorise in FP32, and 3x to 4x with FP16, so pinning everything to FP64 gives that up. The decision is straightforward: where outputs must match a legacy kernel within a tight bound, pay for FP64 in the reference and audit builds; where some divergence is acceptable, keep the faster default path in production and document the divergence the audit measured.

How do you set an error budget for kernel-versus-library differences that auditors will accept?

The error budget we would defend combines an analytical bound for each operation with the measured distribution of differences over repeated, varied runs. Divergence inside the budget is classified as expected rounding; divergence outside it, or with a consistent sign, is investigated as a bug. Two kernels running the same formulas can differ because the order of summation and the use of FMA change the result.

TechnoLynx delivers this through its GPU Performance Engineering work: it defines the bounding method (ULP analysis, repeated runs), builds the harness that compares outputs across precisions and inputs, confirms the budget holds across realistic portfolios and scenarios, and packages it as auditor-facing documentation rerun after kernel or library updates.

Rump’s rigorous bound for a floating-point dot product is at least (n+2)·eps times the unit in the first place of the computed sum, so it grows with the vector length n. Rump also notes that the worst case is very rare, so such bounds usually grossly overestimate the true error. Treat the analytical bound as the line that must never be crossed, not as the expected size of the difference.

cuSOLVER contains a tolerance of this kind. Its iterative-refinement solver stops when the residual norm falls below sqrt(N) × ‖x‖ × ‖A‖ × eps, with eps the machine epsilon of the input/output type. That is a sensible template: the tolerance grows with problem size and with the magnitudes involved, and is anchored to the precision of the output.

Worked example: what an accumulation-order budget looks like

A 2026 arXiv study of optimised depthwise-convolution CUDA kernels compared every output element-wise against a PyTorch reference. Forward outputs matched within float32 precision. The weight gradient, the path most sensitive to accumulation order, showed differences that grew with problem size, reaching a maximum absolute difference of 4.96×10⁻⁴, or about 1.05×10⁻⁶ relative, at the largest configuration. The authors judged this to be within single-precision tolerance and consistent with known properties of floating-point reductions. Three things carry over to a pricing budget: state it in relative terms, expect it to grow with accumulation depth, and set it separately for each output path, since the forward and reduction paths behaved differently.

Why the tolerance itself needs defending

The authors of RESOLVE, a 2026 GPU-kernel validation framework, note that a permissive tolerance can hide bugs while a restrictive one rejects valid optimisations. Their caution is concrete: two KernelBench candidates that had passed the benchmark’s tolerance-based tests were still incorrect. A budget alone is therefore not enough; it needs the determinism and operator-level checks alongside it.

For auditors, supervisory guidance already frames validation this way. The Federal Reserve’s SR 11-7 guidance (2011, US) says validation can set thresholds for acceptable error by analysing the distribution of outcomes around expected values, and that vendor and third-party products need validation even when components are proprietary.

We would document three caveats: re-derive the budget when the precision mode (FP32, FP64, TF32) or the hardware changes; investigate systematic divergence beyond rounding noise as a bug; and be ready to tie the budget to a stated risk tolerance. The parity protocol for a full port (choosing the tolerance, running the parity book, signing off) is covered in our guide to proving kernel equivalence after porting or precision changes.

Proving a faster kernel has not moved the numbers a decision depends on

Do not ship a kernel or library optimisation until two independent checks pass: an error-bound analysis for the specific arithmetic it changed, and a regression suite run against a disclosed, maintained baseline showing no change in the decision-relevant outputs beyond tolerance; if either fails, revert or narrow the optimisation. For a pricing desk, those outputs are per-trade prices, Greeks and anything that triggers a hedge or a limit, not just portfolio aggregates.

TechnoLynx runs this validation inside its GPU Performance Engineering work, starting by agreeing the tolerance for decision-driving outputs, then choosing the mixed-precision strategy, checking the faster kernel against the baseline and error bound, and leaving drift monitoring in place.

An analogy from molecular dynamics: Amber validates its GPU code against its CPU code. Amber keeps a full double-precision model, DPFP, specifically for careful regression testing against its CPU code, while its default mixed model, SPFP, is used for speed. Amber’s developers note that SPFP diverges from CPU trajectories rapidly, yet statistical properties converge to the same values. Agreement in aggregate is still not proof: a GPU free-energy validation in Amber observed that two codes can agree within statistical error while one uses a subtly wrong potential, so energies and gradients should be compared directly. In that comparison, energies agreed to at least six significant figures between the CPU double-precision and GPU mixed-precision models. The pricing equivalent is to compare per-trade values directly between optimised and baseline runs, even when a portfolio-level reconciliation already passes.

A bitwise determinism check is worth running first. RESOLVE requires that repeated executions of a kernel on a fixed input give identical bits, which flags schedule-dependent behaviour with no tolerance needed. Its determinism tester caught injected synchronisation faults and found four previously unreported issues in production mega-kernels, including two clear bugs.

Error bounds are worst-case or statistical estimates, not guarantees for every input, and the baseline has to be maintained or the suite drifts with it. For high-stakes outputs, an independent review on top of the suite is reasonable.

Who can validate our custom CUDA kernels against cuBLAS and cuSOLVER?

TechnoLynx builds the validation layer quant teams need before a faster kernel goes live: a traced account of what cuBLAS and cuSOLVER actually compute, an error budget you can show auditors, and regression gates that catch library changes. This work sits in TechnoLynx’s GPU Performance Engineering engagement, which starts with framing which kernels, precisions and outputs must reconcile. Where the open question is which GPU to buy, LynxBenchAI published results show how devices run the public PyTorch stack; they do not measure your own pricing code.

Building this in-house works when someone on the team already knows cuBLAS math modes, reduction order and statistical classification of differences, and has time free from production support. When choosing a partner, ask for the harness and the documentation as deliverables you keep, not just a verdict. If the need is a one-off check of a single routine, NVIDIA’s cuBLAS and cuSOLVER documentation plus a pedantic FP64 reference run may be enough for a numerically experienced team.

How do you stop a cuBLAS or cuSOLVER upgrade from silently changing rounding behaviour?

Pin the qualified library version and re-certify on every upgrade with a back-to-back regression suite and a numerical-tolerance acceptance bar. Vendor documentation on internal reduction order is incomplete, so empirical regression testing is the real safeguard. NVIDIA’s cuBLAS documentation is explicit that bitwise reproducibility is not guaranteed across toolkit versions, because the implementation may change.

TechnoLynx sets this up as a certification gate: tolerance bars per operation and end to end, a regression harness comparing library outputs with legacy kernel outputs on representative and edge-case inputs, the pinned version with its re-certification trigger written into the upgrade pipeline, and a record of the precision sensitivities found, for the risk sign-off.

Version changes alter more than the last bit. The 2026 study that reconstructed cuBLAS GEMM bitwise found that whether cuBLASLt splits K at all depends on both the architecture and the library version. It also located a defect in cuBLAS itself: at very deep contraction lengths the library summed a whole number of blocks and dropped the tail. A regression suite with deep-K shapes would catch that; a version-to-version spot check on typical shapes would not.

cuBLAS’s fixed-point FP64 emulation may allocate extra workspace dynamically when the per-handle workspace is insufficient, and allocation failures fall back to non-emulated routines; providing enough workspace through cublasSetWorkspace() avoids this. The automatic dynamic precision framework decides per operation whether emulation or native hardware runs. A certification record should therefore name the toolkit and driver versions, GPU architecture, math mode, compute type and workspace size, since NVIDIA’s CUDA guide lists hardware, compiler version and flags, runtime conditions and thread configuration among the conditions for a deterministic result. The same 2026 study offers a stricter model: a descriptor that records only the parameters fixing reduction order (such as where split-K cuts the K axis), so two machines given the same descriptor owe each other the same bits.

The decision rule: if the library documents its numerical behaviour and the version can be pinned, pin it and re-certify on each upgrade; if it does not, treat every update as untrusted until it passes the suite. Bit identity across versions is not the goal, because non-associativity rules it out in general; bounded, documented tolerance is.

Is a GPU-versus-CPU gap non-associativity in the reduction, or a real bug?

Compare the observed deviation with the expected bound under controlled conditions. If it stays within the bound once seeds are matched and deterministic settings are on, it is benign parallel reduction numerical variance; if it exceeds the bound, or persists under matched seeds and deterministic mode, it is a bug to chase in the code. NVIDIA’s CUDA programming guide states that differences between host and device do not necessarily mean the GPU result is incorrect.

TechnoLynx handles this diagnosis within GPU Performance Engineering: it pins down which kernels, precisions and reduction orders are involved, builds deterministic configurations and per-trade instrumentation on the target GPU, and reports a benign-versus-bug classification against the CPU engine.

A difference in the 7th decimal of an FP32 result sits close to the limit of single precision, so summation order is the first suspect. That makes it a hypothesis to test, not a conclusion; Amber’s maintainers, for comparison, expect minor test differences from undefined GPU operation order, limited precision and RNG variation, yet still require the diff to be inspected manually.

Separating Monte Carlo noise from implementation error

TechnoLynx’s approach is to feed both engines identical random inputs: generate the normals once, or use the same quasi-random sequence, and run both CPU and GPU paths over them. Any remaining difference is then implementation or arithmetic, not path variance. Amber’s maintainers list variation in random number generators across GPU hardware among the causes of test differences. Tests for a subtle precision or RNG bug in a cuRAND port should compare per-path distribution statistics as well as the final price, and include a deliberately small path count where a bug cannot hide in the noise.

Can deterministic cuBLAS settings give bit-for-bit results for audit, and what does that cost?

Yes, within limits. NVIDIA states that all cuBLAS routines from a given toolkit version give the same bits on every run on GPUs with the same architecture and number of SMs, though not across toolkit versions. The guarantee no longer holds with multiple active CUDA streams or fixed-point emulation, and routines with atomics-based alternates reproduce run to run when atomics mode is disabled. cuSOLVER’s dense routines are deterministic by default, and cusolverDnSetDeterministicMode() allows non-deterministic results for performance. Shape dependence is a separate issue: cuBLASLt’s batch-invariance mode makes results independent of how a problem is split along M or N, but it is experimental, restricts algorithm choice, disables split-K and may reduce performance, and it is currently supported only on compute capability 10.x GPUs.

On cost, the 2026 cuBLAS reconstruction study cites fixing the reduction order at up to 20% of matrix-multiply throughput. A 2024 study of floating-point non-associativity on GPUs also measured the performance impact of replacing atomics with deterministic alternatives. Where audit builds use deterministic settings that production fast paths do not, the audit has to record which configuration produced which numbers.

Determinism is not correctness, either: the Feltor developers made their parallel simulations bitwise reproducible and still concluded that reproducibility and accuracy alone do not show correct behaviour. Bitwise determinism checks belong alongside a reference comparison and an error budget, never in place of them.

Triage table: benign variance or a bug?

Observation	Likely reading	Next step
Difference within the stated budget, near-zero mean, grows gradually with accumulation depth	Reduction-order rounding	Record it in the budget documentation
Repeated GPU runs on the same input give different bits	Schedule-dependent behaviour: atomics, streams or a race	Disable atomics, single stream, rerun the determinism check
Difference vanishes under pedantic FP64 and TF32 off	Precision path, not logic	Decide whether production may keep the fast path
Difference persists under matched random inputs and deterministic mode, beyond budget	Implementation bug	Operator-level diff to find the first divergent checkpoint
Difference far above rounding on a few products	Ill-conditioned problem amplifying round-off	Trade-level root cause, below

When the largest deviations cluster on a few trades

When the outliers concentrate on a small subset of trades, the question changes from “is this noise?” to “why these trades?”. Rank per-trade deviation, group the top tail by product type and model path, reproduce those trades on the GPU in deterministic mode, and escalate anything that stays outside tolerance as a real bug. The reason to look for structure is that ill-conditioned problems amplify small rounding differences into large output differences, so clustering on one product family may reflect that family’s numerics rather than a defect, and only the deterministic reproduction settles which.

What about steep skew and low-liquidity strikes, where the numerics diverge more?

We found no source that measures kernel-versus-library divergence in these regimes specifically, so what follows is engineering judgement, not a finding. The mechanism is textbook: inverting a price to an implied volatility divides any price error by vega, so where vega is small (far wings, short-dated far strikes) a rounding-level price difference becomes a much larger volatility difference. NVIDIA applies the same logic in cuSOLVER, recommending one SVD routine only for well-conditioned data because its formulation can greatly amplify errors. We would set tolerances per region of the surface rather than one global number, include deliberately ill-conditioned strikes in the regression set, and accept that a match in calm regimes says nothing about stressed ones.

Did the speed-up come from the kernel change, or from quietly dropped precision?

Attribute a throughput gain by pairing a kernel-level profile with a per-change divergence breakdown, then re-measure with each risky optimisation reverted in turn. If the gain survives the revert and divergence stays inside budget, the speed-up is genuine; if removing one fused op or precision cast collapses both the speed-up and the divergence, that change was buying speed with accuracy.

TechnoLynx runs this as an ablation study: device-level profiling of per-kernel time, selective reverts of each precision cast, fusion or layout change on the target hardware, throughput and divergence measured per configuration, and a written attribution of where the speed came from and what residual risk remains.

Some optimisations trade precision by design. nvcc’s --use_fast_math replaces math calls with intrinsics, running faster at reduced precision and accuracy, and NVIDIA recommends replacing calls selectively, only where the gain justifies the changed accuracy. cuBLASLt’s FP8 fast-accumulation flag skips periodic promotion of intermediates to higher precision on Ada and Hopper. In NVIDIA’s own emulation case study, cutting FP64 emulation to 39 mantissa bits gave nearly 3x end-to-end speed-up while output stayed consistent to 12 significant digits. Each is legitimate when it is chosen knowingly and covered by the budget, and a problem when it arrives as an unexamined default.

Measure the gain where it counts, too. In the 2026 depthwise-convolution study, a 3.26x kernel speed-up shrank to 1.29x per training epoch because framework overhead, synchronisation and other operations took a growing share of runtime. For per-kernel attribution, see our guide to GPU profiler tools such as Nsight Systems and Nsight Compute, and profiling-based validation for GPU workloads on checking that a profile reflects the real workload. A weak reference can hide real precision loss, so an attribution is only as good as its baseline.

Whether the speed-up is still acceptable after reverting the riskiest change is a business judgement, not a numerical one. A faster kernel is only worth keeping if you can name which change bought the speed and show that its numerical cost sits inside a budget someone outside the team can check.