Deciding When to Offload R-tree Queries to the GPU

The break-even batch size is not a constant you can copy from a blog post — it is a property of your device, your data, and your result size, and it can differ by two orders of magnitude between a phone and a workstation. Guessing it wrong is expensive in both directions: offload too early and you pay upload and readback latency for savings that never materialize; offload too late and a single CPU thread grinds through a query flood that thousands of GPU lanes would have absorbed. The reliable move is to measure on the hardware the code actually runs on, then let the router pick the path from that number. This page gives you a self-contained harness that sweeps batch sizes, times the CPU flatbush path and the GPU compute path fairly, and reports the crossover so the offload decision is data, not a hunch.

The model behind the measurement — the fixed upload-and-readback overhead, the per-query costs, and the inequality — is derived in compute shader vs CPU spatial indexing. Here we make it concrete and runnable.

Runnable reference implementation

A pointer-based R-tree node flattened for the GPU Two strips show the same R-tree node. The pointer layout scatters the node across memory: a bounding box of four f32 values, a child pointer of eight bytes, a child count, and a parent pointer, with the children themselves at unrelated addresses. The flattened layout packs the same node into a fixed 32-byte record — four f32 bounds, a u32 first-child index into the same array, a u32 child count, and eight bytes of padding to keep the stride aligned — so a workgroup reads a run of sibling nodes in one coalesced transaction. BYTE OFFSET 0 4 8 12 16 20 24 28 32 36 bounds (4 × f32) child ptr count parent ptr Pointer layout scattered each pointer dereference is an uncoalesced read at an unrelated address bounds (4 × f32) first count pad Flattened 32 B stride siblings are contiguous — one workgroup reads a whole level in a few transactions Build the flattened array in breadth-first order so a level is a contiguous range.
Flattening is what makes a tree usable on a GPU at all. The traversal logic barely changes; what changes is that following a child index reads memory the wavefront has probably already fetched, instead of jumping somewhere the cache has never seen.

The harness builds a flatbush R-tree and an equivalent GPU grid over the same feature boxes, then times a batch of range queries on each path across a geometric sweep of batch sizes. The CPU path is timed with performance.now() around the tree walk; the GPU path is timed end to end — upload, dispatch, and readback — because that whole span is what the CPU alternative avoids. It returns the smallest batch size at which the GPU total beats the CPU total: your measured .

typescript
import Flatbush from "flatbush";

interface BreakEven { qStar: number | null; samples: { Q: number; cpuMs: number; gpuMs: number }[]; }

// --- CPU baseline: flatbush range queries, timed with the wall clock ---
function timeCpu(index: Flatbush, queries: Float32Array, Q: number): number {
  const t0 = performance.now();
  let sink = 0;
  for (let q = 0; q < Q; q++) {
    const o = q * 4;
    sink += index.search(queries[o], queries[o + 1], queries[o + 2], queries[o + 3]).length;
  }
  const t1 = performance.now();
  if (sink < 0) throw new Error("unreachable");   // defeat dead-code elimination
  return t1 - t0;                                  // ~ Q * t_cpu, single thread
}

// --- GPU path: upload + dispatch + readback, timed end to end ---
async function timeGpu(
  device: GPUDevice, pipeline: GPUComputePipeline,
  queryBuf: GPUBuffer, hitsBuf: GPUBuffer, readback: GPUBuffer,
  bind: GPUBindGroup, queries: Float32Array, Q: number,
): Promise<number> {
  const view = queries.subarray(0, Q * 4);
  const t0 = performance.now();
  device.queue.writeBuffer(queryBuf, 0, view);           // t_up
  const enc = device.createCommandEncoder();
  const pass = enc.beginComputePass();
  pass.setPipeline(pipeline);
  pass.setBindGroup(0, bind);
  pass.dispatchWorkgroups(Math.ceil(Q / 64));            // one lane per query
  pass.end();
  enc.copyBufferToBuffer(hitsBuf, 0, readback, 0, Q * 4);
  device.queue.submit([enc.finish()]);
  await readback.mapAsync(GPUMapMode.READ);              // t_down: forces completion
  readback.getMappedRange();                             // touch the result
  readback.unmap();
  const t1 = performance.now();
  return t1 - t0;                                        // t_up + dispatch + t_down
}

// --- Sweep batch sizes and find the crossover ---
export async function findBreakEven(
  device: GPUDevice, pipeline: GPUComputePipeline,
  queryBuf: GPUBuffer, hitsBuf: GPUBuffer, readback: GPUBuffer, bind: GPUBindGroup,
  index: Flatbush, queries: Float32Array,
  batchSizes: number[], reps: number,
): Promise<BreakEven> {
  const samples: BreakEven["samples"] = [];
  let qStar: number | null = null;

  // Warm both paths so JIT and pipeline compilation stay out of the medians.
  timeCpu(index, queries, Math.min(1024, queries.length / 4));
  await timeGpu(device, pipeline, queryBuf, hitsBuf, readback, bind, queries, 1024);

  for (const Q of batchSizes) {
    const cpu: number[] = [];
    const gpu: number[] = [];
    for (let r = 0; r < reps; r++) {
      cpu.push(timeCpu(index, queries, Q));
      gpu.push(await timeGpu(device, pipeline, queryBuf, hitsBuf, readback, bind, queries, Q));
    }
    const cpuMs = median(cpu);
    const gpuMs = median(gpu);
    samples.push({ Q, cpuMs, gpuMs });
    if (qStar === null && gpuMs < cpuMs) qStar = Q;       // first Q where offload wins
  }
  return { qStar, samples };
}

function median(xs: number[]): number {
  const s = [...xs].sort((a, b) => a - b);
  const m = s.length >> 1;
  return s.length % 2 ? s[m] : (s[m - 1] + s[m]) / 2;
}

// --- Build the CPU index and a matching random query set for the sweep ---
export function buildIndexAndQueries(featBoxes: Float32Array, maxQ: number, extent: Float32Array): {
  index: Flatbush; queries: Float32Array;
} {
  const n = featBoxes.length / 4;
  const index = new Flatbush(n);
  for (let i = 0; i < n; i++) {
    const o = i * 4;
    index.add(featBoxes[o], featBoxes[o + 1], featBoxes[o + 2], featBoxes[o + 3]);
  }
  index.finish();

  // Random small query boxes across the extent — same set feeds both paths.
  const [minX, minY, maxX, maxY] = extent;
  const w = (maxX - minX) * 0.02, h = (maxY - minY) * 0.02;
  const queries = new Float32Array(maxQ * 4);
  for (let q = 0; q < maxQ; q++) {
    const cx = minX + Math.random() * (maxX - minX);
    const cy = minY + Math.random() * (maxY - minY);
    const o = q * 4;
    queries[o] = cx - w; queries[o + 1] = cy - h;
    queries[o + 2] = cx + w; queries[o + 3] = cy + h;
  }
  return { index, queries };
}

Drive it with a geometric sweep so the crossover is bracketed, not stepped past:

typescript
const { index, queries } = buildIndexAndQueries(featBoxes, 500_000, extent);
// gridPipeline, queryBuf, hitsBuf, readback, and bind come from the grid setup
// in the compute-vs-CPU-indexing reference; the grid queries the same featBoxes.
const result = await findBreakEven(
  device, gridPipeline, queryBuf, hitsBuf, readback, bind,
  index, queries,
  [500, 1_000, 5_000, 10_000, 50_000, 100_000, 250_000, 500_000], // batch sweep
  15,                                                              // reps per size
);
console.log(`break-even Q* = ${result.qStar ?? "GPU never won in range"}`);

The WGSL grid-query kernel this drives is the one from the compute shader vs CPU spatial indexing reference; reuse it verbatim so the only variable across the sweep is the batch size.

Parameter reference

What each traversal strategy costs on a GPU A comparison of three traversal strategies across three properties. A recursive descent has no GPU equivalent because WGSL has no recursion, needs unbounded stack, and is not implementable. An explicit stack per invocation is implementable, costs registers proportional to tree depth, and diverges heavily because each lane descends a different path. A level-by-level breadth-first sweep is implementable, uses no per-lane stack, and keeps lanes uniform because every lane processes the same level at the same time. TRAVERSAL STRATEGY ON A GPU Implementable? Cost Divergence Recursive descent no — — Per-lane stack yes registers × depth heavy Level sweep yes one pass per level none Depth is bounded and small, so the number of sweep passes is a compile-time-ish constant.
The level sweep costs more total work — it visits nodes a depth-first descent would have pruned — and still usually wins, because uniform lanes at full occupancy beat a smarter algorithm that leaves most of the wavefront idle.
Parameter Typical value Guidance
batchSizes 500 → 500k, geometric Space multiplicatively (×2 or ×5); a linear sweep wastes reps far from the crossover.
reps 10–20 Median over enough runs to reject GC and scheduler spikes; fewer is noisy near .
Query box size ~2% of extent Result count drives ; match it to your real query selectivity.
Feature count your dataset Larger raises CPU slightly but mostly raises ; keep it realistic.
Warm-up batch 1024 Runs both paths once before timing so JIT and pipeline compile stay out of medians.
@workgroup_size 64 Query lanes diverge on cell counts; 64 keeps partially-idle workgroups cheap.
Result payload count (u32) per query A count keeps minimal; full id lists inflate readback and move right.
dataInVram shortcut boolean If queries are GPU-generated, skip the sweep — offload wins; is zero.

The query box size is the parameter that most often makes a measured mislead: sweep with tiny boxes that return few candidates and the CPU looks unbeatable; sweep with large boxes that return thousands each and balloons, pulling down. Match the selectivity to production, or measure a range of box sizes and store a per selectivity band.

Failure modes

Offload failures and what they really indicate A table of three failure signatures with cause and response. A GPU path slower than the CPU at every size usually means the batch is too small for the fixed dispatch cost to amortise, and the response is to batch more queries or stay on the CPU. A GPU path that is fast but produces a stuttering map means results are being read back inside the frame, and the response is to move the readback to a ring read two frames later. A GPU path returning slightly different results from the CPU means the bounds comparison used a different epsilon or a different inclusivity at the edges, and the response is to define the predicate once and port it exactly. OFFLOAD FAILURE · CAUSE · RESPONSE Cause Response slower at every size batch too small batch more, or stay on CPU fast but stuttering readback inside the frame read from a ring, later results differ slightly edge inclusivity differs one predicate, ported exactly Half-open intervals on both sides is the convention worth standardising on — pick one and write it down.
The third row deserves a test rather than a fix. Run both implementations over the same fixture in CI and compare the result sets exactly; a boundary disagreement that shows up as three features in a million is impossible to find by inspection and trivial to find by assertion.
  • Timing that hides the readback. Timing only dispatchWorkgroups and not the mapAsync readback measures a GPU cost the CPU path never has to beat, making offload look free. Detection: measured is implausibly small and offload underperforms in production. Fix: time the full upload-dispatch-readback span as the harness does; that span is exactly what the CPU alternative avoids.
  • Dead-code-eliminated CPU baseline. A JIT that proves the flatbush.search results are unused can delete the loop, reporting a near-zero CPU time and a false of zero. Detection: cpuMs is orders of magnitude below a hand check. Fix: accumulate a value derived from every result — the sink sum here — and reference it after the loop so the work cannot be elided.
  • Warm-up folded into the median. Without the warm-up call, the first sweep entries carry JIT and pipeline-compilation cost and distort the crossover. Detection: the smallest batch sizes are anomalously slow on one path. Fix: run both paths once before the timed sweep and take medians over reps, never a single sample.
  • Crossover outside the swept range. If the sweep tops out below the real , the harness reports null and you wrongly conclude offload never helps. Detection: qStar is null while GPU time is still falling relative to CPU at the largest batch. Fix: extend batchSizes upward, or recognize that for this selectivity the CPU genuinely wins across your operating range and keep queries on the CPU.

Backend and Python interop note

When the feature boxes originate server-side, the break-even shifts in your favour if the backend ships them in a layout the GPU grid can consume without repacking. Emit bounding boxes as a contiguous, C-ordered float32 array of [minX, minY, maxX, maxY] rows so a single writeBuffer lands them in the grid-build buffer, and pre-sort features by a coarse cell key so the on-GPU grid build does less scatter work:

python
import numpy as np
import geopandas as gpd

gdf = gpd.read_parquet("features_z12.parquet")          # pre-projected, Web Mercator metres
bounds = gdf.geometry.bounds                            # minx, miny, maxx, maxy columns

boxes = np.empty((len(gdf), 4), dtype=np.float32)
boxes[:, 0] = bounds["minx"].to_numpy()
boxes[:, 1] = bounds["miny"].to_numpy()
boxes[:, 2] = bounds["maxx"].to_numpy()
boxes[:, 3] = bounds["maxy"].to_numpy()

# Coarse cell key for locality; a pre-sorted feed cuts GPU grid-build scatter.
cell = (np.floor(boxes[:, 1] / 1000.0).astype(np.int64) << 20) | \
        np.floor(boxes[:, 0] / 1000.0).astype(np.int64)
order = np.argsort(cell, kind="stable")
payload = np.ascontiguousarray(boxes[order]).tobytes()  # ships straight into the grid buffer

Projecting on the server matters as much for offload as it does for aggregation: a grid cell index is only stable when every box is in one linear space, so mixing geographic degrees with metres corrupts the index on both the CPU and the GPU regardless of how the break-even lands.

Up: Compute Shader vs CPU Spatial Indexing