TL;DR: tsbootstrap's fast path used to build every bootstrap replicate into one giant batched tensor: on a routine job, a 160 MB tensor pushed out to main memory and hauled back to produce one summary number per replicate. We rebuilt the hot path around two rules: never materialize an array whose only consumer is a reduction, and never carry state you can derive. Every head-to-head benchmark cell now favors the fused path; the settled benchmark grid sustains 3.1x to 20x on the longer series. The last stubborn cell, where the runtime was mostly Python building random-number seeds, dropped from 13.1 ms to 0.55 ms in the redesign's validation run once we deleted the seed objects.
The full replicate tensor is still available when you want it, and on that materializing path our default backend remains slower than the incumbent. That row is printed below with the rest.
Early this summer, in the middle of a performance push I thought was going well, we ran tsbootstrap's block-bootstrap engine head-to-head against arch, a mature econometrics library whose bootstrap loop has had years of polish, for the first time. We lost. At a realistic workload, a few thousand replicates of a few-thousand-point series, the whole call came back several times slower end-to-end, on the machine we developed on, before the compiled engine this article describes existed. Every release we had ever shipped carried that loss; users hitting exactly that workload had been paying it the whole time, and nothing in our tooling would ever have told us.
{/* The opening's historical loss is deliberately number-free: the underlying dev-machine measurements are not yet in the number ledger, so no multiplier may be printed here until rows are added. */}
I had not seen it coming: our benchmark suite had only ever compared the library against its own history. A regression suite tells you when you get slower than yesterday; it cannot tell you that you ship slower than the other library, at the exact workload where users would notice. The day a head-to-head existed, the loss appeared. It had been there all along.
The profiler's first answer was the ordinary one. Most of the measured deficit was removable interpreter tax: a per-replicate Python index loop, thousands of seed objects spawned per call, stacked object wrapping. But beneath it sat a second, structural problem, the reason the reflexive fix could not win either. The reflexive fix for Python overhead is to batch: replace the loop with a single vectorized NumPy operation. Batching removes the interpreter from the inner loop and swaps it for B-fold traffic through main memory (B being the number of bootstrap replicates, the count every cost in this story scales with) because the batched design materializes every replicate at once. The incumbent's replicate loop is ordinary Python too. It wins by keeping its working set (the data its loop actively touches) small enough to live in cache, consuming each resample the moment it exists. We were not going to beat that by batching harder.
Here is the machine's side of the argument. A CPU core computes only on what it holds in its registers. Think of that as your hands. Behind them sit the caches: L1, a desk you can reach without looking up; L2, the shelf behind you; L3, shared with the other cores, a cabinet down the hall. Behind all of that is DRAM: main memory, the warehouse across the street, bigger and slower at every step, orders of magnitude in latency top to bottom. The hardware hides none of that from your runtime; it hides it only from your source code.
Now put a bootstrap on that machine. A bootstrap, one time for anyone who has not run one: it re-draws your data thousands of times to see how much a statistic would wobble across plausible alternate samples, and each redrawn copy is a replicate. Take a modest, entirely realistic job: a series of n = 10,000 observations, B = 2,000 bootstrap replicates, float64 (8 bytes per number). Build all the replicates at once, the way a batched NumPy design wants to, and you have asked for a B x n tensor: 10,000 x 2,000 x 8 bytes = 160 MB. That is far beyond the tens of megabytes of L3 a core complex sees. Every byte of that tensor is written out to main memory, then hauled back in so the statistic can read it once, and if the statistic is a mean, the answer you keep is 2,000 numbers. Sixteen kilobytes, bought with 160 MB of round trips.
{/* Asset BUILT: TensorScrub wraps this beat. Readers with scroll-timeline support and no reduced-motion preference get the scroll-scrubbed stage (assemble / flood / collapse, byte counter 160 MB -> 16 kB, rest caption 20 MB where 1.94 GB stood); everyone else (reduced-motion, print, browsers without animation-timeline: scroll()) instantly gets this static three-panel comic (SVG twin at src/figures/charts/tsb-03-tensor-never-fits.svg, class prefix tnf-) with the caption and table twin. The switch is pure CSS (@supports + @media), no JS involved. */}
The batched design
The tsbootstrap that ships today never builds that tensor on its fast path. The story of how it stopped: The first wall is memory traffic: the bytes themselves, moving through the hierarchy. The second took longer to see, because it hides inside the language runtime and the library's own design, not the hardware: the cost of state (allocated objects and seeded generators) that exists only to make randomness happen in the right order. We stopped materializing arrays, and we stopped carrying state.
Wall one: memory traffic
Start with what the batched design actually does. To bootstrap a series with blocks (contiguous runs resampled whole, so the short-range dependence that makes a series a series survives inside each block; the first piece in this series is about what happens when it doesn't), you draw random block starts, expand them into resample indices, gather the data through those indices, and compute your statistic on each replicate. The batched version does each of those as one array operation over all B replicates: build a (B, n) index matrix, gather into a (B, n, d) values tensor (d is just the number of columns in your series), then reduce along the middle axis.
Count the work in each cell of that tensor. Eight bytes written on the way out. Eight bytes read back for the reduction. Roughly one floating-point addition, one addition per sixteen bytes moved, where a compute-bound kernel like a matrix multiply does hundreds of operations for every byte it pulls. That ratio, the arithmetic intensity (how much computing you do per byte you move), is the number that decides whether your workload is limited by the processor or by the memory system. Ours is about as low as arithmetic intensity gets. The cores were never the constraint; they spent the benchmark waiting on memory. The gather is slow not because the code is interpreted but because its arithmetic intensity is far too low to keep the processor busy: it is a bytes problem being graded as a FLOPs problem.
The incumbent's loop, seen this way, stops looking old-fashioned and starts looking correct. One replicate at a time means one (n, d) resample at a time: a working set that fits in cache, is consumed by the statistic immediately, and is gone before the next one arrives. That working set never grows with B, so DRAM barely enters its story. The batched design's traffic grows linearly with B, and once the tensor outgrows the last cache level, every additional replicate's bytes take the round trip through main memory. NumPy's programming model quietly couples two decisions that are really separate: batching the dispatch and batching the bytes. The loop was winning because it never built the thing we were building.
Here is the whole argument as accounting, on the worked example: n = 10,000, B = 2,000, float64, one mean per replicate. Redo it on a napkin for your own workload; the ratio is the point.
| quantity | materialize path | fused, streaming path |
|---|---|---|
| intermediate allocated | 160 MB, the full replicate tensor | one (n,) index row + one (n, d) scratch per thread (~80 kB each) |
| bytes through DRAM | at least twice the tensor: written out, read back | consumed in cache, per work item |
| useful output kept | 16 kB, one mean per replicate | 16 kB, the same answer |
| bytes moved per byte of answer | ~20,000 : 1 | ~1 : 1 |
{/* This table is the bytes-per-FLOP device. If clean hardware-counter measurements land, ONE measured row (last-level-cache misses per iteration, empty-loop baseline subtracted) may be added; if the counters come back muddy the table ships exactly as-is, never a sketched row. */}
None of this is a new insight: compilers have fused loops to avoid intermediate arrays for decades, and every numerical programmer eventually meets the memory wall. What is worth knowing is that the same wall now dominates transformer training and serving. FlashAttention, the kernel that reshaped both, describes its own core move in the caption of its first figure: it uses "tiling to prevent materialization of the large N × N attention matrix", and that one refusal, plus fusing the surrounding operations into a single kernel, yielded a 7.6x speedup on GPT-2's attention. The hierarchy it plays against is the GPU's version of ours: an A100's main memory streams at up to 2.0 TB/s, while the tiny on-chip SRAM (192 KB per streaming multiprocessor) runs roughly ten times faster. Same shape, same conclusion: keep the working set in fast memory; never write the big intermediate at all. FlashAttention is the proof that the old compiler principle still binds at datacenter scale.
One pass, and nothing left behind
The fix is a single compiled kernel that does everything the batched pipeline did, per replicate and in one pass, without ever building the B-wide tensor. Per-worker scratch still exists, but it is n-sized and transient, not B times the data. The standard name for this is kernel fusion: instead of separate passes that hand each other arrays, one loop body builds, gathers, and reduces for one replicate and emits only the statistic.
Each replicate is one unit of parallel work. The kernel gives it a private (n,) index buffer and a private (n, d) scratch, builds that replicate's resample indices into the buffer, gathers values through them column by column, applies the reducer (mean, variance, standard deviation, or a quantile) in the same loop body, and writes one row of the (B, d) output. The 160 MB tensor from the opening never exists; peak memory is the source data, per-worker scratch, and the 16 kB we actually wanted. The kernel is compiled to machine code the first time it runs, a just-in-time (JIT) compilation, and that first-call cost is paid off the hot path by an explicit warm-up sweep when the compiled backend loads.
The same byte-mindedness shows up in smaller knobs: resample index tensors are
int32rather than int64, halving their footprint (with an explicit raise the moment an index would exceed what 32 bits can address, rather than a silent overflow), and the returned simulation tensor can be requested in float32 while every model fit and every reduction staysfloat64. Storage-dtype choices, both of them, made for traffic.
From the outside it is one argument. The materializing call and the fused call take the same method spec and the same seed:
from tsbootstrap import MovingBlock, bootstrap, bootstrap_reduce
x = np.random.default_rng(0).standard_normal(10_000)
# the materializing path: every replicate, as arrays
res = bootstrap(
x, method=MovingBlock(block_length="auto"),
n_bootstraps=2_000, random_state=0)
# the fused path: same resampling, one pass, only the statistic survives
red = bootstrap_reduce(
x, method=MovingBlock(block_length="auto"), statistic="mean",
n_bootstraps=2_000, random_state=0, backend="compiled")
lo, hi = red.quantile(0.05), red.quantile(0.95)
The memory consequence is not subtle. In a peak-memory sweep over B (measured 2026-07-04 on version 0.4.0; the streaming path is unchanged since), holding the full replicate tensor at B = 50,000 on an n = 2,000 series cost 1.94 GB; the reduce path over the same fifty thousand replicates peaked at 20 MB, roughly 97x lighter. The reduce path's peak grows with the number of statistics you keep, not with the length of the paths (B rows of answers instead of B full series, which is the whole point), so fifty thousand replicates stops being a capacity-planning question.
Fusing the pass won the big workloads. It did not win everything. The cell it refused to win is where the second wall shows itself.
Wall two: the state you carry
After the fused kernel shipped, one benchmark cell refused to move: the smallest one, an i.i.d. resample (independent draws, no blocks, the simplest bootstrap there is) of a short series: 200 observations, 999 replicates. Still losing to the incumbent, run after run, while its neighbors fell by multiples. The cell's total time held at 13.1 milliseconds like it was nailed there. A kernel this small finishes in microseconds; something else owned the clock. I finally profiled the entire call, not the part I was proud of.
The kernel was innocent. It was barely running.
About 95% of that cell's wall time was spent before the kernel launched, in Python, setting up random-number seeds. And this is not a fixed dispatch cost: it is an O(B) scaling tax, a serial object-allocation loop that dominates the clock on small jobs and grows linearly to tax the biggest ones before the kernel even starts.
{/* Figure: the "Derive, don't carry" split-screen loop (HyperFrames, 40 s = two ~20 s passes; pass B fires the fan-out in a fixed shuffled permutation while the checksum fingerprint stays identical, the thread-count-invariance beat). Figure.astro serves the .gif reference as a muted-loop
The same workload, twice, in the state redesign
Reproducibility in tsbootstrap is a per-replicate contract: replicate i is bound to its own random stream, stable no matter how the work is chunked or how many threads run it. The obvious implementation is to materialize the streams: spawn B seed objects up front, one per replicate, and hand them out. Correct, and O(B) Python objects allocated on every call. For big workloads the kernel dwarfed it; for small ones the setup was the runtime. We had stopped materializing the data tensor and were still materializing a tensor of state.
Terminology, one time: a PRNG, short for pseudorandom number generator, is the deterministic algorithm behind every "random" draw in scientific computing. Give it the same seed and it replays the same sequence, which is what makes a bootstrap reproducible at all.
Wall two has two faces. The first is ordinary Python overhead: allocating B objects, looping B times, on every call. The second is deeper than Python, and it is the reason GPU people solved this a decade before we needed it. Conventional PRNGs are, in the words of the Random123 paper that changed that corner of computing, "designed as sequentially dependent state transformations": each draw mutates a state block that the next draw depends on. State means memory: a Mersenne Twister instance drags around about 2.5 kB of it. Worse than the bytes is the sequencing. State that must evolve in order resists being handed to ten thousand independent workers.
The alternative the Random123 authors proposed, counter-based generators that "require little or no memory for state," inverts the design. Stop treating randomness as a machine you advance; treat it as a pure function: hand it a key and a counter, get back bits, no state anywhere. Any worker can compute draw k of stream b without computing any other draw first.
So we deleted the seed objects, and the nailed-down cell fell: 13.1 ms before, 0.55 ms after, on the same box in the redesign's validation run, with identical kernel arithmetic. Nothing about the math changed. The setup work simply no longer exists.
Here is what replaced it. The executor now takes a single 128-bit root key and the replicate count, and each backend derives its streams from there; inside the compiled kernel, each replicate's key is computed from the root and the replicate index by a few integer mixing rounds: nothing allocated, nothing carried, no Python in sight. Our disease was Pythonic, object-allocation cost rather than per-thread kilobytes, but it was exactly the disease Random123 diagnosed: state carried in sequence where a pure function should be. The key for replicate b is now a pure function of (root, b); any specific draw is a pure function of that key and a counter. That sentence is the whole design, and its consequences are the guarantees: run the replicates on one thread or eight, in any order, and the output is bitwise identical, because there is no shared state for a schedule to touch. That guarantee is not pedantry: it is what makes a CI failure mean a real regression rather than thread-schedule noise, and what lets a published interval replay byte for byte months later on a machine with a different core count.
The generator doing the drawing deserves its precise name. The compiled kernels run Philox-4x32-10, the counter-based generator from that same Random123 paper (Salmon et al., SC 2011): ten rounds, the canonical multiplier constants. It is the identical algorithm, with the identical constants, that NVIDIA ships as cuRAND's Philox device generator, and the same splittable-key philosophy JAX built its entire random API around. The 32-bit variant was chosen deliberately so the integer index streams can one day be made bit-identical between this CPU kernel and an array-backend GPU port. A statistics library on a laptop and a datacenter training run, drawing their randomness from the same ten rounds of the same function, for the same reason: state does not parallelize; functions do. None of this touches the default: the shipping default backend is unchanged and byte-for-byte compatible with previous releases, and the Philox path is opt-in via backend="compiled": equal to the default in distribution, not bit-identical to it.
A hand-rolled generator inside a compiled kernel is also exactly the kind of code that fails silently: a wrong-but-plausible RNG sails through smoke tests and quietly skews every interval you ever publish. So the kernel is never trusted on its own word. The key-derivation math lives twice (once inside the kernel, once in a deliberately slow, dependency-free reference module that takes a different arithmetic path), so a bug in one implementation cannot hide behind agreement with itself; the round function is pinned against the published Random123 known-answer vectors (the paper's own input-output pairs, which the algorithm must reproduce exactly); and the fused, materializing, and index-only paths all share one index-generation code path, so they cannot drift apart per seed.
The default backend's own machinery: plain NumPy PCG64 streams, one per replicate; even its chunking uses a fixed chunk size (never derived from available RAM), so accumulation order and results cannot vary across machines. And the compiled path raises loudly if the compiled extra is missing rather than falling back silently.
Two optimizations we measured and deleted
The rebuild was not a straight line. Two other ideas were measured during development, on the development machine, and both were reverted; the verdicts, not the multipliers, are what matters.
| hypothesis | verdict | why |
|---|---|---|
| shared scratch: pre-allocate one buffer block per chunk of replicates instead of per-work-item temporaries | implemented, proven bitwise-identical, net regression on every cell, reverted | the allocator was already pooling per-iteration buffers; the shared scratch defeated SIMD auto-vectorization |
| cache-blocked gather: sort or tile the resample reads for locality | every locality fix measured as a net loss, disqualified | the stall is bandwidth-bound, and any reorder breaks bitwise identity |
The shared scratch was classic advice (allocation is overhead, pooling is virtue), and it lost even single-threaded, because the JIT runtime already pools small per-iteration allocations, and rows of a shared block are strided views (rows spaced apart in memory rather than contiguous), which the auto-vectorizer, the compiler pass that processes several array elements per instruction, refuses to touch. The cache-blocked gather's diagnosis even measured out (the gather really is memory-stall-bound), but every locality-improving reorder cost more in bookkeeping than it recovered in bandwidth, and because floating-point addition is not associative, reordering the accumulation changes the last bits of the result, which the bitwise-reproducibility contract forbids.
Both hypotheses were plausible enough that a design debate would have shipped them. The measurements killed them in an afternoon, and both verdicts are pinned so neither gets re-attempted. Burying an idea is far cheaper than maintaining one.
The scoreboard, losses included
Everything above was measured on dedicated hardware so the numbers mean something: 8-core AMD EPYC Milan cloud instances with dedicated vCPUs, 30 GB of RAM, Ubuntu 24.04 (kernel 6.8), exact package versions pinned in each run's committed manifest. Sixteen cells: the observation-resampling methods crossed with a range of series lengths and replicate counts, our fused reduce path against the incumbent's equivalent loop. Every cell's time is the median of 15 settled timed runs after an explicit warm-up, on the shipping versions of both libraries (tsbootstrap 0.6.1 against arch 8.0.0), one receipt, measured 2026-07-11.
All sixteen cells favor the fused path, and the concession below holds beside them. On the sustained workloads (the longer series where the engines are past their fixed overheads) the settled grid holds 3.1x to 20x. The settled timing exists because single-draw timing failed us first: when we re-timed the grid with one draw per cell, the incumbent's own control numbers moved 22 to 74 percent between runs of identical code on identical instance types (the box noise exceeded the effect), so nothing measured single-draw is quoted here. The smallest cells post far larger multiples still, and we do not quote them either: at that size the comparison measures the incumbent's per-replicate dispatch overhead as much as it measures our kernel.
One more disclosure that belongs next to the range: the fused path runs parallel across all eight cores, and the incumbent's loop is single-threaded by construction. That concurrency is part of what the fused design buys, not a handicap we imposed on the comparison, but you should know the multiple contains it.
And one path still loses. If you want the full (B, n, d) replicate tensor as plain NumPy arrays (the materializing path, the thing the whole article is about not doing), the incumbent's loop is faster than our default path in every cell of the grid, by 1.2x to 3.1x on the settled 2026-07-11 grid, depending on the cell. That row is a design position, printed with the rest. Materializing is the incumbent's home game, and it plays it well; we spent our budget on the path where the statistic is the product, because for the statistic-shaped questions a bootstrap exists to answer, nobody needed the tensor in the first place.
{/* Chart asset regenerated 2026-07-11 (settled-grid adoption): SVG twin at src/figures/charts/tsb-03-results-grid.svg, regenerated from data by docs/media-src/tsb-03-results-grid/gen.py (JSON copies in its data/ dir). ALL THREE series come from the settled 2026-07-11 grid (vs_arch_ccx33_2026-07-11_settled.json, median of 15 settled samples per cell); the earlier single-draw grids are superseded (their arch control drifted 22-74% between runs). Two small-workload cells (moving and circular block, n=200, B=10,000) render as capped arrows past the 40x clip, no on-axis multiple; badge on the smallest i.i.d. cell is hoisted above the row line and scoped to the validation run. Alt, caption, and table twin verified against the regenerated SVG. */}
The grid, wins and losses on one chart, every series from the settled 2026-07-11 grid (median of 15 timed runs per cell). Fused reduce path: all sixteen cells favor tsbootstrap, 3.1x to 20x on the sustained workloads, a multiple that includes the fused path
{/* Chart asset BUILT 2026-07-11: SVG twin at src/figures/charts/tsb-03-cache-counters.svg, regenerated from data by docs/media-src/tsb-03-cache-counters/gen.py (counter JSON copy in its data/ dir). All values computed from the JSON; render matches the alt text and table twin. Screenshot-QA'd at 928px and 390px. */}
What the hardware counted while both paths ran the same workload (a stationary-block resample at the worked example
The whole rebuild, in three moves:
| the move | the reason |
|---|---|
| never materialize an array whose only consumer is a reduction | the tensor's bytes cost more than its arithmetic; the answer was 16 kB, the intermediate was 160 MB |
| fuse the pass: build, gather, reduce in one loop body | a working set that fits in cache is the difference between cache latency and DRAM round trips |
| derive state, never carry it | B seed objects were 95% of a workload; a key derived from (root, replicate) is free, order-independent, and bitwise reproducible |
{ q: "Why would a bootstrap be limited by memory rather than by computation?", a: "Because its arithmetic intensity is very low: resampling moves eight-byte values around and then does roughly one floating-point operation per value for a statistic like the mean. When a batched implementation materializes all *B* replicates at once, the resulting tensor is far larger than the caches, so every byte takes a round trip through main memory. The processor spends the workload waiting on that traffic, not computing. The fix is to process one replicate at a time inside a fused kernel so the working set stays cache-resident." },
{ q: "What is a counter-based random number generator?", a: "A conventional pseudorandom generator is a state machine: each draw updates internal state that the next draw depends on, which forces sequential execution and per-stream memory. A counter-based generator such as Philox is a pure function instead: it maps a (key, counter) pair straight to random bits, with no carried state. That means any parallel worker can compute any draw of any stream independently and in any order. tsbootstrap's compiled backend derives a Philox-4x32-10 key for each replicate from a single 128-bit root key inside the kernel." },
{ q: "Does the compiled backend change my results?", a: "The default backend is unchanged: NumPy PCG64 streams, one per replicate, reproducing previous releases byte for byte. The compiled backend is opt-in and uses its own counter-based streams, so its replicates are equal in distribution to the default path but not bit-identical to it. Within the compiled backend, results are bitwise reproducible for a given seed regardless of thread count or execution order, and the library raises an explicit error if the compiled extra is unavailable rather than silently falling back." },
{ q: "When is the fused compiled path the wrong choice?", a: "When you actually need the full materialized tensor of resampled series (for a custom diagnostic that inspects whole paths, for example), the materializing path exists for exactly that, and the incumbent's per-replicate loop is quicker at producing the fully materialized array than our default path. The fused path is built for the common case where the bootstrap's output is a statistic per replicate, and there it wins the benchmark grid in every cell." },
]} />
The fastest code in this library is the code that no longer runs: the tensor nobody builds, the seed objects nobody allocates. Performance work sounds like adding cleverness; the two biggest wins here were subtractions. Somewhere in your own hot path there is probably a 160 MB array whose only job is to be read once and averaged, and next to it a few thousand small state objects built on every call, because carrying state was the design that came to mind first, and main memory bills you by the byte for both, whether or not you ever look. You find these things the way we did: profile the whole call, count the bytes. The FLOPs were never the bill.
The first piece in this series was about a resampler quietly deleting the structure in your data. This one was about a design quietly deleting your performance, byte by byte, in arrays you never needed and state you never had to carry. What happens when the data refuses to be a rectangle? tsbootstrap already packs panels of unequal-length series into one flat array behind a row-offset table (compressed-sparse-row bookkeeping, the layout sparse linear algebra has leaned on for decades), and how a thousand ragged series share one fused pass without padding a single one is where this track goes next.
Top comments (0)