The GPU Precision Trap: How an fp32 Backtest on Apple Metal Silently Returns Garbage

Part of the "Backtests Without Illusions" series.

We ported our parameter-sweep backtest to the GPU and it got 2,796× faster. The numbers it returned were plausible. They were also, on the first working version, complete garbage — off by a factor of two hundred — and nothing crashed, nothing warned, nothing looked wrong.

This is the trap that makes GPU backtesting on Apple Silicon dangerous in a way CPU code is not: Metal has no float64. Every number your backtest touches on an Apple GPU is a 32-bit float, whether you asked for that or not. And the single most tempting way to vectorize a moving average — the O(n) prefix-sum trick every performance-minded quant reaches for — is precisely the formulation that fp32 cannot survive at price scale. It does not error out. It runs at full speed and hands you a plausible-looking equity curve built on a moving average that is wrong by a factor of 211.

The fix is the interesting part, because it is not "use more bits" (you can't) and it is not "be more careful" (the naive code is already careful). The fix is to compute the same mathematical quantity by a different sum — a direct windowed convolution — that keeps every intermediate value small enough for fp32 to represent exactly. Same WMA, same result to seven significant figures, and 55.9× faster than single-thread numba on the CPU. This article is the autopsy: why Metal forces fp32 on you, exactly where the obvious formulation overflows, why the correct one doesn't, and how we proved — by counting trades, not by eyeballing curves — that the fast version and the honest fp64 version agree.

All numbers here are measured on an Apple M2 Max, 150,000 bars × 80 parameter combinations, best-of-3, from the repo's scripts/bench_param_sweep.py (the M5 GPU method, commit 97eadaf), scripts/engine_multitf_gpu.py (04d71e8), and the design notes in scripts/GPU_NOTES.md.

Metal has no float64

On a CPU, a Python/numpy backtest runs in double precision by default. float64 gives you a 52-bit mantissa: integers are represented exactly up to 2⁵³ ≈ 9×10¹⁵, and relative precision is about 1.1×10⁻¹⁶. You almost never think about it, because at price scale — a BTC close around 30,000, a cumulative sum reaching 10¹⁴ — fp64 has margin to spare.

Apple's GPU does not offer that margin, because it does not offer float64 at all. The Metal Shading Language has no double type; MLX, PyTorch-MPS, and every other framework that targets Apple Silicon inherits this. There is no flag to turn it on and no slow-but-correct fallback. If your computation touches the GPU, it happens in fp32:

• Mantissa: 23 bits (24 with the implicit leading bit). Integers are exact only up to 2²⁴ = 16,777,216 ≈ 1.6×10⁷.
• Relative precision: ~1.2×10⁻⁷. Roughly seven significant decimal digits, and no more.

That 1.6×10⁷ ceiling is the whole story. It sounds generous — sixteen million — until you notice that a backtest routinely builds intermediate quantities far larger than that, and the moment an intermediate crosses 1.6×10⁷, fp32 can no longer even represent consecutive integers, let alone the fractional structure you need. The precision doesn't degrade gracefully; it falls off a cliff, and the cliff is at a number your data sails past without a second thought.

The trap is that fp32 is almost always fine. Most of a backtest — prices, returns, PnL, Sharpe — lives comfortably in the range where seven digits is plenty. So the naive port works, passes a smoke test, produces sane-looking output. The failure is localized to exactly one operation, the one place where an intermediate blows past 10⁷, and that operation is the one everybody vectorizes first.

The tempting formulation: WMA in one prefix-sum sweep

Our strategy leans on Hull moving averages. An HMA is three weighted moving averages composed together; the HMA3 variant is four. A weighted moving average with a linear kernel over a window of length p is

$\mathrm{WMA}_t = \frac{\sum_{k=1}^{p} k \cdot x_{t-p+k}}{p(p+1)/2}.$

Sweeping thousands of parameter combinations over 150k bars, the WMA convolutions are the cost. So the instinct — the correct instinct, on a CPU — is to make each WMA O(n) instead of O(n·p) using prefix sums. You precompute two cumulative sums once,

$S_1[t] = \sum_{j \le t} x_j, \qquad S_2[t] = \sum_{j \le t} j \cdot x_j,$

and then any window's linearly-weighted sum collapses to a handful of differences and index shifts of S1 and S2. No per-window loop, no per-window reduction — two cumsum passes and the whole matrix of WMAs falls out of array arithmetic. It vectorizes beautifully, it maps perfectly onto a GPU's parallel-scan primitive, and on fp64 it is exactly right.

It is also the single worst thing you can do in fp32, and the reason is hiding in S2.

S2 = mx.cumsum(j * price)     # j is the global bar index: 0, 1, 2, ... , n-1

The term j · price is the problem. With j running up to 150,000 and price around 30,000, the last term alone is 4.5×10⁹, and S2 is the running sum of 150,000 such terms. It does not stay at price scale. It climbs into the territory where fp32 has already stopped being able to count.

Where it overflows: the arithmetic of the trap

Let's put the orders of magnitude next to each other, because that is where the whole failure lives.

S2 = cumsum(j · price) reaches roughly price · n²/2 ≈ 30,000 · (150,000)²/2 ≈ 3×10¹⁴. Call it ~10¹⁴. Now recall the fp32 exact-integer ceiling: ~1.6×10⁷. The running sum overshoots the last integer fp32 can represent exactly by seven orders of magnitude.

What does that mean concretely? Near 10¹⁴, the gap between two representable fp32 numbers — one unit in the last place, the ULP — is about 2²³ ≈ 8×10⁶. So once S2 is up in the 10¹⁴ range, it is only known to within ±8 million. Every value it stores has been rounded to the nearest multiple of ~8×10⁶.

Now watch what the WMA recovery does with that. To extract a single window's weighted sum, you subtract two neighboring S2 values (plus S1 corrections). Those two S2 values are each ~10¹⁴, each carrying ±8×10⁶ of rounding noise. Their true difference — the windowed quantity you actually want — corresponds, after normalization, to a WMA on the order of the price itself, ~3×10⁴. So the arithmetic is:

$\underbrace{(10^{14} \pm 8\times10^6)}_{\text{fp32}} - \underbrace{(10^{14} \pm 8\times10^6)}_{\text{fp32}} = \underbrace{\text{a quantity you need to a few parts in } 10^4}_{\text{true answer}}$

This is catastrophic cancellation in its purest form: the rounding error of each operand (±8×10⁶) is larger than the answer you are trying to recover. The signal is smaller than the noise floor of the numbers it is extracted from. It is not that you lose a few digits — you lose all of them, and what comes back is dominated by the accumulated rounding of the cumsum.

The measured consequence, from GPU_NOTES.md: for a WMA computed this way on 150k bars at price ~30,000, the maximum relative error against fp64 reaches ~211. Not 211 percent — 211×. The computed moving average can be two orders of magnitude away from the true one. And here is the part that makes it a trap rather than a bug: it runs to completion and returns finite, plausible numbers. No overflow to infinity, no NaN, no exception. A moving average that is wrong by 211× still looks like a moving average — it is smooth, it is finite, it is in roughly the right ballpark on the bars where cancellation happens to be mild — so it sails through every sanity check that isn't a direct comparison against a trusted reference. You get a full backtest, a full equity curve, a full set of "optimal" parameters, all built on an indicator that is fiction.

The fix is not more precision — it is a different sum

The reflex, once you see the error, is to reach for more precision — accumulate in fp64, or use compensated (Kahan) summation. On Metal, the first is simply unavailable. But you don't need either, because the problem was never the number of bits. The problem was the formulation. The prefix-sum trick manufactures 10¹⁴-scale intermediates and then subtracts them back down; the magnitudes it creates are an artifact of the algorithm, not of the answer. Choose a formulation that never creates them and fp32 is fine.

That formulation is the definition itself: a direct windowed convolution. Instead of two global cumulative sums, slide the length-p linear kernel across the series and sum it in place. Each output is a sum of at most p ≈ 200 terms, and every term is weight × price where the weights are normalized to sum to 1 — so each term is on the order of price / p, every partial sum stays around price scale (~3×10⁴), and no intermediate ever comes within six orders of magnitude of the fp32 ceiling. There is nothing to cancel because nothing ever inflated.

In MLX this is one primitive — mx.conv1d — which is exactly what GPUs are built to do fast:

def _mx_wma_valid(x, period):
    w = mx.arange(1, period + 1, dtype=mx.float32) / (period * (period + 1) / 2.0)
    return mx.conv1d(x.reshape(1, -1, 1), w.reshape(1, period, 1), padding=0).reshape(-1)

Same WMA, mathematically identical to the prefix-sum version and to the CPU's fp64 vec_wma/nb_wma. But now the measured maximum relative error against fp64 is 8.2×10⁻⁷ — right at the fp32 noise floor of ~1.2×10⁻⁷, seven significant figures of agreement. The formulation that looks slower on paper (O(n·p) instead of O(n)) is the only one that is correct, and — because it is a dense convolution the GPU parallelizes across both bars and windows at once — it is also blisteringly fast. We went from a relative error of 211 to 8×10⁻⁷ by changing how we sum, not how many bits we sum in.

Two practical notes that fall out of doing it this way. First, MLX won't propagate NaNs through a conv1d the way numpy does, so the warm-up region (the first p−1 bars, where a windowed average isn't defined) can't be marked with NaN on the GPU. We don't need it to be: the valid start of every series is known analytically, invalid prefixes are filled with zeros that are never read, and the NaN padding is restored on the CPU afterward — bit-for-bit the same validity mask as the vectorized and numba versions. Second, the whole sweep shares one cand_close series and reuses windows heavily across combos, so a single batched conv1d with many output channels computes every unique WMA the sweep needs in one GPU call, materialized with one mx.eval().

Proving you didn't fall in: parity by trade count

Here is the uncomfortable question the last section should raise: if a WMA wrong by 211× still looks like a WMA, how do you know the 8×10⁻⁷ version is actually right and not just wrong more subtly? You cannot eyeball it. You need an invariant that a downstream, discrete part of the pipeline exposes — and a backtest hands you a perfect one: the trades.

The other GPU methods in the ladder (M0–M4) run entirely in fp64, so we hold them to a strict equivalence assert — identical trade counts, PnL matching to atol=1e-6. The fp32 GPU method (M5) cannot pass that by construction, and quietly loosening the assert for everyone to accommodate it would be exactly the kind of dishonesty this series exists to fight. So M5 gets its own quantitative parity report, report_equiv_fp32, comparing its extracted trades against the fp64 reference.

The mechanism of any residual disagreement is worth naming precisely, because it is not the cancellation catastrophe — it is the ordinary, tiny fp32 rounding you'd expect. The strategy fires on the crossover of two Hull averages, h vs h3. A relative error of ~1×10⁻⁶ on an indicator at price ~30,000 is an absolute wobble of ~0.03. On the vast majority of bars the two curves are further apart than that and the crossover is unambiguous. But on a borderline bar — where h − h3 is itself within 0.03 of zero — that wobble can flip the sign of the comparison, moving one crossover by a single bar, adding or removing one trade.

This is why "fraction of combos that differ" is a worthless health metric, and our first parity check embarrassed itself by using it. At 150k bars each combo has thousands of crossovers, so at least one borderline bar shows up on almost every combo — 37 of 80 combos "differed," which sounds alarming and means nothing. The metric that matters is by how much:

• PnL delta across all 80 combos: max |Δ| = 1.843 percentage points, max relative = 1.25×10⁻²; crash threshold 5 p.p.
• Trade-count drift per combo: max |Δn| = 4 trades out of thousands, max relative = 2.5×10⁻³; crash threshold 1%.
• In aggregate: 90 shifted trades out of 479,016 — 0.019%.

Ninety trades out of nearly half a million, every one of them a borderline crossover nudged by a rounding wobble smaller than a price tick, none of it anywhere near the crash thresholds. That is the signature of an fp32 method that is correct — small, bounded, explainable disagreement — and it is a completely different animal from a relative error of 211. The thresholds exist to catch a broken formulation masquerading as "well, it's just fp32"; the real deltas come in an order of magnitude under them. The trade count is the oracle the equity curve refused to be.

The payoff, and where the GPU stops helping

With correctness established, the speed is worth stating — and then honestly qualified, because the GPU's advantage is not uniform across the pipeline.

On the pure WMA convolutions in isolation — the operation the whole method exists to accelerate — the fp32 conv1d batch runs 55.9× faster than single-thread numba, at that 8.2×10⁻⁷ relative error. That is the clean, apples-to-apples GPU-vs-CPU number: same math, one thread of compiled CPU code against the Metal GPU.

But a sweep is not only convolutions. Once the HMA/HMA3 matrices are computed on the GPU, the trades still have to be extracted — an O(n) walk over each combo's crossovers — and that we do on the CPU in fp64, reusing the exact trade semantics of the other methods rather than re-implementing them on the GPU. The end-to-end timed() figure includes everything: kernel warm-up excluded (symmetric with excluding numba's compile), but GPU→CPU transfer and CPU trade extraction included. On 150k bars × 80 combos, best-of-3, M2 Max:

*M0 extrapolated from a uniform 5-combo sample.

The full-engine M5 does 779 combos/s — 2,796× over the pandas baseline, 19.6× over serial numba, and 3.2× over the entire 12-core CPU pool running numba (M4). One GPU beats every CPU core the machine has, threefold.

Now the honest qualifier: notice that the end-to-end GPU advantage (19.6× over M2) is smaller than the convolution-only advantage (55.9× over numba). That gap is Amdahl's law arriving on schedule. The GPU annihilates the convolutions so completely that they stop being the bottleneck; what is left — the O(n) CPU trade extraction that the GPU sped up not at all — now dominates M5's wall time. This is the same lesson the speed-ladder and the IPC-tax pieces in this series keep landing on: past a point, the win is not "make the fast part faster," it is orchestration — where the data lives, which stage is now serial, what you are paying to move between device and host. Chasing a hypothetical M6 that pushes trade extraction into a custom Metal kernel would claw back only the shrinking CPU slice, which is why we didn't build it.

The general lesson: silent numerical garbage is the default

Step back from Hull averages and MLX, because the trap generalizes far past this one indicator.

The seductive pitch of GPU backtesting is "one big matrix": stack every parameter combination into a tensor, run the whole sweep as a handful of dense array ops, let the hardware eat it. That pitch is real — the speedups above are real. But it quietly changes the numerical regime underneath you, and the change is invisible in the code. On the CPU your defaults protected you: fp64, NaN propagation, a cumsum that could run to 10¹⁴ without noticing. Move the same array expression to Metal and you are in fp32 with a hard 1.6×10⁷ integer ceiling, and the identical line of code — cumsum(j * price) — goes from exact to garbage. Nothing in the syntax warns you. The compiler is happy. The output is finite and plausible. fp32 doesn't fail loudly; it fails politely, with numbers.

The three habits that actually protect you are cheap:

1. Know where your intermediates live, not just your inputs and outputs. The inputs (prices ~10⁴) and outputs (a WMA ~10⁴) were both comfortably inside fp32's exact range. The disaster was entirely in a hidden intermediate (S2 ~10¹⁴) that neither the API nor the types made visible. Before trusting any fp32 reduction, ask what the largest partial sum reaches — and if it crosses ~10⁷, change the formulation.
2. Prefer formulations that keep magnitudes bounded. Direct convolution over prefix sums; local windows over global scans; centering/differencing before summing rather than after. Big-then-cancel is the anti-pattern. The correct algorithm is often the one that looks asymptotically worse on paper but never manufactures a quantity it has to cancel away.
3. Validate against an fp64 oracle through a discrete invariant. Do not compare curves; compare something quantized and downstream — trade count, number of crossovers, position-change events. A discrete invariant turns a silent 211× error into a screaming assertion failure, and turns an acceptable 8×10⁻⁷ error into a small, bounded, explainable delta. This is the same discipline as the one-bar shift test for look-ahead bias: a cheap diagnostic that converts an invisible failure into a visible one.

None of this is exotic numerical analysis. It is the ordinary hygiene of not trusting a fast backtest until a slow, trusted one has ratified it — extended to the one place where the language stops warning you that precision has silently dropped by twelve orders of magnitude.

Takeaways

1. Apple's GPU has no float64 — every GPU number in your backtest is fp32. Integers are exact only to ~1.6×10⁷ and precision is ~1.2×10⁻⁷. There is no flag, no fallback. Most of a backtest survives this; exactly one operation usually doesn't.
2. The prefix-sum WMA is the trap. cumsum(j · price) climbs to ~10¹⁴, seven orders past fp32's exact ceiling, and recovering a window forces you to subtract two such numbers whose rounding error (±8×10⁶) already dwarfs the answer. Measured max relative error: 211×. It never crashes — it returns plausible garbage.
3. The fix is a different sum, not more bits. A direct windowed convolution (mx.conv1d) keeps every partial sum near price scale, so fp32 holds seven honest digits: 8.2×10⁻⁷ relative error, and 55.9× faster than single-thread numba. You cannot buy fp64 on Metal, and you don't need to.
4. Verify with a discrete invariant, never the curve. Trade-count parity caught it: fp32 conv1d disagreed with fp64 on 90 of 479,016 trades (0.019%), all borderline crossovers, all far under the crash thresholds — the signature of a correct method, unmistakably different from a 211× error. "Fraction of combos that differ" is a decoy metric; measure by how much.
5. The full sweep does 779 combos/s — 2,796× over the pandas baseline, 3.2× over the entire 12-core CPU pool — but the end-to-end win (19.6× over serial numba) is smaller than the convolution-only win (55.9×) because the CPU trade extraction is now the bottleneck. Past a point, speed is orchestration, not arithmetic.

The GPU port was 2,796× faster and, on its first working version, completely wrong — and the two facts had nothing to do with each other. The speed was real. The garbage was a hidden 10¹⁴ intermediate that fp32 could not hold and no error message would mention. If a backtest gets dramatically faster and the numbers still look fine, that is not confirmation. On Metal, "looks fine" is exactly what a relative error of 211 looks like.

This is the GPU rung of the ladder this series has been climbing: the backtest-engine speed ladder, the IPC tax of multiprocessing, the look-ahead taxonomy of leaks, and the objective-function design that decides what "good" even means. Speed is worthless if it is fast at computing the wrong number.