Deciding When to Offload R-tree Queries to the GPU

The break-even batch size $Q^{*}$ 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 $Q^{*}$ 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 $Q^{*} = (t_{u p} + t_{d o w n}) / (t_{c p u} - t_{g p u})$ — is derived in compute shader vs CPU spatial indexing. Here we make it concrete and runnable.

Runnable reference implementation

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 $Q^{*}$ .

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

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 $Q^{*}$ .
Query box size	~2% of extent	Result count $k$ drives $t_{c p u}$ ; match it to your real query selectivity.
Feature count $N$	your dataset	Larger $N$ raises CPU $lo g N$ slightly but mostly raises $k$ ; 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 $t_{d o w n}$ minimal; full id lists inflate readback and move $Q^{*}$ right.
`dataInVram` shortcut	boolean	If queries are GPU-generated, skip the sweep — offload wins; $t_{u p}$ is zero.

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

Failure modes

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 $Q^{*}$ 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 $Q^{*}$ 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 $Q^{*}$ , 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.

Compute Shader vs CPU Spatial Indexing — the break-even model, decision matrix, and grid-query kernel this harness measures.
Performance Tuning & Profiling for WebGPU Spatial Pipelines — the profiling context these timings feed into.
WGSL Spatial Algorithms on the GPU — building the GPU-resident grid the offload path queries.
Benchmarking WebGPU and WebGL 2.0 Point Rendering — the sibling harness that applies the same fair-timing discipline to rendering.

Up: Compute Shader vs CPU Spatial Indexing

Deciding When to Offload R-tree Queries to the GPU

Runnable reference implementation

Parameter reference

Failure modes

Backend and Python interop note

Related