The CPU backend
The correctness oracle. Deliberately naive, deliberately slow, and the reason the Vulkan backend can be trusted at all.
It is a reference, not a fallback
Correctness here is a chain: the CPU backend is checked against PyTorch for semantics, then the Vulkan backend against the CPU for kernel bugs. That second link only means anything because the CPU backend shares vkML's exact semantics — so a mismatch is unambiguously a kernel defect rather than a difference of convention.
Two consequences follow, and both are enforced rather than hoped for:
- CPU support must be a superset of Vulkan support. A test asserts it
directly, because it broke once: widening the Vulkan
supports()gates for float16 made the GPU accept operators the CPU still rejected, and every one of those had a GPU result no oracle could check while the suite stayed green. - Performance comes last here. The matmul is a naive triple loop with no blocking and no vectorisation, and the source says so where someone might otherwise "fix" it.
Measured on this machine, the CPU backend is roughly 116× slower than PyTorch on the same work. That is the cost of the trade, not a defect — and it is why the README's CPU-only install path is for checking answers, not for getting them.
Pairwise summation, and why it is not optional
This is the single most important numerical decision in the backend, and it is a correctness requirement rather than a refinement.
Sequential summation of n values in float32 has a worst-case relative error of
n·ε, with ε = 2⁻²³ ≈ 1.19e-7. The project's gate against PyTorch is
atol = rtol = 1e-5. So:
| n | What it is | Sequential error | Against the 1e-5 gate |
|---|---|---|---|
| 784 | MNIST input features | ≈ 9.3e-5 | fails by ~9× |
| 4096 | transformer hidden | ≈ 4.9e-4 | fails by ~49× |
A naive accumulator does not merely lose a little precision — it misses the acceptance
criterion outright, and does so in a way that looks like a kernel bug. Pairwise
summation splits the range recursively, giving a bound of about
(B + log₂(n/B))·ε for a sequential base case of size B:
| B = 32, n | Pairwise error | Margin |
|---|---|---|
| 784 | ≈ 4.3e-6 | 2.3× |
| 4096 | ≈ 4.6e-6 | 2.2× |
| 16384 | ≈ 4.9e-6 | 2.0× |
kPairwiseBlock = 32 is chosen to keep at least 2× margin out to n = 16384, which
covers every reduction length these models produce. NumPy uses 128 and would still pass in
practice, because real rounding errors random-walk rather than aligning — but designing to the
worst-case bound costs nothing measurable and removes a whole class of "why is this test flaky
at large K" investigation later.
It also mirrors what the GPU has to do anyway. That device has no global
float atomicAdd, so its reductions must be tree-shaped regardless — subgroup
reduction, then shared memory, then a deterministic second pass. Matching the CPU reference to
that structure is what keeps the two comparable, and it is why the GEMM shaders use
BK = 32: the same block size, so both backends fold K identically.
float16 is a storage format, never an accumulator
Every kernel computes in float whatever it stores. The widening lives in exactly
two overloads in iterate.h, which is what keeps the f32 and f16 paths from drifting
apart.
It matters most in matmul: an f16 accumulator over K = 784 would lose roughly three decimal digits, well outside the 1e-3 the f16 tolerance allows. The dot product accumulates in float regardless of the storage type.
Strided iteration
iterate.h converts a flat logical index into a byte offset per operand.
Broadcasting is handled for free: a stride of 0 contributes nothing to the
offset, so every index along that axis reads the same element — no special case, no branch.
That is why a broadcast operand costs no memory and no separate code path anywhere in the backend.
Randomness is counter-based
rand and dropout use Philox4x32-10, a counter-based
generator rather than a stateful one. The value at each index is a pure function of
(seed, offset, index).
That property is what lets the GPU produce the same draw as the CPU without any sequencing between invocations — there is no stream to keep in step, because there is no stream. It is also why the signatures take a seed and an offset rather than carrying hidden state.
Each value takes the top 24 bits of a 32-bit output, which is the float significand's width, so every result is exactly representable. The round keys are the fractional parts of the golden ratio and of √3, the paper's own "nothing up my sleeve" construction.
Allocation
CPU allocations are aligned to 64 bytes — one cache line on every CPU this will run on, and the alignment AVX-512 wants. The backend will not be hand-vectorised, but aligning costs nothing and removes a variable if the compiler ever vectorises part of it.