Technology · Discipline

Streaming statistics in CUDA

Welford and Chan moments, histogram quantiles and compare-and-swap float atomics, written once in a device library and tested on their own, so every metric kernel is numerically stable in single precision.

Inference FabricCUDA device librarySingle precision
Welford per thread · Chan merge per level
n, μ, M₂n, μ, M₂n, μ, M₂n, μ, M₂n, μ, M₂n, μ, M₂n, μ, M₂n, μ, M₂
warp mergewarp merge
block merge, broadcast to every thread
  • 256bins in the shared-memory block histogram
  • 10 KBof shared memory for a median of millions of distances
  • 0.4%resolution of that median, about

Three ways single precision goes wrong

Statistics over a hundred thousand vectors in single precision across thousands of threads go wrong in familiar ways: a naive variance cancels catastrophically, a sum depends on the order of the adds, and the minimum of a float cannot be taken atomically because the hardware has no such instruction. The device library under my CUDA engine solves each of those once, and every metric kernel in the Inference Fabric reads from it.

One sum, two orders of arrival
Naive, order oneNaive, order twoChan merge, either order
lengths illustrative: the naive totals drift apart, the merged ones agree

Every primitive has its own CUDA test, including the block Welford and its broadcast form, the reductions and the tile machinery.

The kernels that call tested primitives inherit tested ground.

Stable moments, finished in one launch

Moments are accumulated with Welford's method, each thread carrying a count, a running mean and a sum of squared deviations, and merged with Chan's combination formula at the warp and block levels, so the moments come out stable whatever the order of arrival. The convention is population statistics throughout, written into the reference so the GPU and CPU agree on the denominator.

Block reductions for sums, minima and maxima are templated with broadcast variants that leave the result in every thread, which the streaming kernels use to finish a cluster without a second launch.

Plain reduction
Broadcast reduction
the result sits in one thread; a second launch reads itthe result sits in every thread; the kernel carries on

Quantiles without keeping the data

A block histogram lives in shared memory, every value is binned as it is produced, and the median or any percentile is read from the cumulative counts by one thread, with a documented resolution instead of a sort. That is how the characteristic scale in Physics as a lens on embeddings finds the median of millions of pairwise distances.

every value binned as it is producedmedian: the bin where the running count passes half

Atomics the hardware does not provide

Floating-point minimum and maximum atomics are built from compare-and-swap on the integer bit pattern, looping until the swap succeeds, and are issued once per tile, never per element. They are what let the diversity kernels find a cluster's distance range in one pass before the histogram pass that needs it.

Attempt 1
Attempt 2
Attempt 3
read the bits, keep the smaller value, compare-and-swap; another thread got there first twice, the third swap holds
Once per value
Once per tile
Each lit cell is one atomic on global memory; the tile reduces first and contends once.