Skip to content

9. Parallelism & determinism

The idea

Histogram GBT parallelizes naturally (features are independent during histogram building, rows are independent during predict), but the naive version silently sacrifices something valuable: reproducibility. Floating-point addition is not associative: \((a+b)+c \ne a+(b+c)\) in the last bit, so any parallel scheme that lets thread scheduling change the order of additions produces models that differ run to run. The reference libraries mostly accept this (fixed thread count → same model, different count → close-but-different). bonsai's contract is the same, held deliberately rather than by accident: when training on the CPU, the model is bit-identical across runs at a fixed thread count, and for everything except the u8 histogram fill, bit-identical to a serial run at any count. The one exception was bought, measured, and is this chapter's closing story.

The "on the CPU" is not a hedge, it is the boundary. GPU training accumulates histogram cells with device atomics whose landing order is the scheduler's, not the configuration's; the cells are int64 fixed point, so the sum is one integer in any order and the same GPU fit produces the same model bytes every run on one device and one build. What the device does not promise is the rest of the CPU rule: identity across devices and toolkits, where the gradient kernels' transcendentals and fused multiply-adds differ, and identity with the host fit, which accumulates float cells and agrees to 1e-4 (chapter 10).

The math (of non-determinism)

A histogram cell is a sum of \(\sim 10^5\) doubles. Change the addition order and the result moves in the last ulp; a split whose gain ties within an ulp can then flip; one flipped split changes a subtree; 200 trees later the models disagree visibly. Determinism is not about the average case: it's about ties, and gradient data produces exact ties routinely (symmetric gradients, duplicated rows).

The design rule that buys determinism: no unordered cross-thread floating-point reductions. Every accumulator is either written by exactly one thread in serial order, or assembled from per-thread partials merged in a fixed order, so the addition order is a pure function of configuration, never of scheduling.

In bonsai

  • The seam: parallel::for_each_index(n, f) in include/bonsai/parallel.hpp: an OpenMP parallel for (dynamic schedule, size-scaled chunks so asymmetric P/E cores stay busy) with a serial fallback. Worker count: [parallel] n_threads, 0 = auto (hardware threads, capped at 16 and at any cgroup CPU quota, so a quota-limited container sizes to what it is allowed to burn rather than to the cores it is shown). An explicit count over that quota is left alone and warns once on stderr, because the resolved count is part of the model's identity and a clamp would change the model bytes a container produces. Design notes in the invariants.
  • Histogram fill (src/grower.cpp): the routing is per node, not per bin width. u16 bins keep the feature-parallel shape (fill_columns): each feature's histogram owned by one thread, filled in row order, determinism by ownership. u8 bins take that same fill whenever the node holds at least a quarter of the rows (k_col_fill_den = 4), and only sparser u8 nodes take the row-wise fill (populate_many/fill_sparse, decisions 49 and 105): row chunks stream a row-major mirror, tiled into 2048-feature blocks, and a node touched by more than one thread range owns a partial per extra range, reduced in fixed range order; the partition derives from the configured thread count, which is exactly where the fixed-count contract comes from. So the dense nodes near the root stay bit-identical at any thread count, and the contract narrows only for the sparse nodes deeper in the tree.
  • Split scan (src/split.cpp): the levelwise finder is feature-parallel, each feature finding its own best independently and the winners merging serially in feature order, so gain ties break exactly as a serial scan would; the per-node finder walks its features serially and inherits that tie-break outright.
  • Node-parallel split scan (src/level_step.hpp): a level runs one worker per frontier node (LevelStep::host_find), so each node's serial scan is independent and nothing is shared between them.
  • Row-parallel predict / gradients / scatter: one row, one thread, no shared accumulator.

Relaxing the contract was a deliberate trade, not an accident (decision 49): ground-truth instrumentation showed deep sparse nodes filling at a fifth of the dense rate, a cache problem no feature-parallel scan could fix, and the row-wise fill bought a measured 1.6–1.7× on real cells for the narrowing above. The models are still exactly reproducible; you just have to hold n_threads fixed, like every reference library. What did NOT survive contact with measurement: a per-parity two-way split promising 1.6× in a microbenchmark delivered nothing in the real loop (the microbenchmark's arrays were cache-resident; the streaming loop's are not), and an LLC-size-based automatic layout choice was rejected because it would have made models hardware-dependent.

Try it

import numpy as np
import bonsai

rng = np.random.default_rng(0)
X = rng.normal(size=(5000, 10)).astype(np.float32)
y = (X[:, 0] * 2.0 + X[:, 1] + rng.normal(0, 0.1, 5000)).astype(np.float32)


def fit(n_threads):
    # u16 bins (max_bin > 255) are bit-identical at any thread count.
    return bonsai.BonsaiRegressor(
        n_iters=50, max_bin=511, n_threads=n_threads).fit(X, y)


p1 = np.asarray(fit(1).predict(X))
p8 = np.asarray(fit(8).predict(X))
print("bit-identical:", np.array_equal(p1, p8))

With u16 bins (max_bin > 255) this compares equal at any thread pair; with u8 bins, re-run either side at the same n_threads twice and compare those: run-to-run identity is the contract. Try the same with XGBoost's nthread and diff the dumped models.

Gotchas & war stories

  • thread_local inside a parallel region is per-worker. The ordered grad/hess gather buffer is a main-thread thread_local; naming it inside the OpenMP lambda resolved to each worker's own empty vector and out-of-bounds writes segfaulted 74 tests at once. Fix: capture the data pointers before entering the region. If you use OpenMP with any TLS, this one is waiting for you.
  • Two OpenMP runtimes in one process deadlock. The Python extension originally linked Homebrew's libomp.dylib; the moment XGBoost (which bundles its own libomp) built a DMatrix in the same process, one OpenMP call stack spanned two different libomp images and parked forever at a join barrier. Fix, standard for wheels: link libomp statically into the module and export only the module-init symbol (BONSAI_OPENMP_STATIC=ON; decision 36). bonsai, XGBoost, and LightGBM now interleave in one process.
  • NUMA first-touch can halve your fill rate silently. The row-wise fill's partial histograms were first zeroed by the main thread, which homed every page on one socket of a dual-socket EPYC; workers on the other socket then paid remote-memory latency for every add (2× fill penalty, 3.1× thread scaling). Per-block zeroing into make_unique_for_overwrite storage fixed both (5.8× scaling). Where memory is touched first decides where it lives.
  • User time lies under OpenMP. Idle workers spin-wait (blocktime), so user CPU time looks saturated even when threads are starved at barriers. Profile with a sampler and read the stacks, not the totals: that's how the join-barrier bottlenecks in the v0.2.0 round were found.