Measuring the Crossover Point for GPU Spatial Joins

“Should this spatial join run on the GPU” has no general answer, because the crossover depends on four things that vary per application: how many queries arrive together, how selective they are, whether the results are needed on the CPU, and how large the indexed set is. What does generalise is the method for finding out. A harness that varies those four parameters over a real dataset and reports both implementations’ timings answers the question in an afternoon and keeps answering it as the data changes. This page is that harness, the parameters worth sweeping, and the reading of the resulting surface. It is one stage of compute shader vs CPU spatial indexing.

The four parameters, ranked by how much they move the answer A table of four sweep parameters with the range worth covering and how strongly each moves the crossover. Query batch size, swept from one to a thousand, moves it most, because the GPU fixed cost is amortised across the batch. Whether the result must be read back moves it almost as much, because a round trip often costs more than the join. Selectivity, the fraction of features a query matches, moves it moderately through atomic contention and result size. Feature count, the axis most people sweep first, moves it least because both sides scale with it. SWEEP PARAMETERS · BY INFLUENCE Range Influence Query batch size 1 … 1000 strongest Readback needed true / false near-decisive Selectivity 0.001 … 0.1 moderate Feature count 10 K … 10 M weakest Sweep readback as two separate surfaces rather than as a parameter — the cases barely compare.
The ordering is the useful finding. Feature count is the axis that gets swept first and explains the least, because both implementations scale with it in much the same way.

Runnable reference implementation

The harness runs both implementations over the same fixture and reports timings that are directly comparable, which means controlling for the things that otherwise dominate.

typescript
interface JoinCase {
  featureCount: number;     // size of the indexed set
  queryCount: number;       // how many arrive in one batch
  selectivity: number;      // fraction of features a query matches
  readback: boolean;        // does the CPU need the result this frame?
}

async function measure(c: JoinCase, gpu: GpuJoin, cpu: CpuJoin): Promise<Result> {
  const features = makeFixture(c.featureCount);
  const queries = makeQueries(c.queryCount, c.selectivity);

  // Warm both: JIT for the CPU, pipeline compilation and clocks for the GPU.
  for (let i = 0; i < 5; i++) { cpu.run(features, queries); }
  await gpu.warm(features, queries);

  const cpuMs = median(times(9, () => timed(() => cpu.run(features, queries))));
  const gpuMs = median(await timesAsync(9, () => gpu.timed(features, queries, c.readback)));

  return { case: c, cpuMs, gpuMs, winner: gpuMs < cpuMs ? "gpu" : "cpu" };
}

Two controls carry most of the weight. Warming matters on both sides — the CPU path needs the JIT to have optimised the hot loop, and the GPU path needs pipelines compiled and clocks at steady state. And the median of nine runs rather than a mean of three, because both sides have occasional long runs for unrelated reasons and a mean lets one of them decide the answer.

Parameter reference

Parameter Sweep over Why it matters
featureCount 10 K … 10 M, ×10 The axis everyone expects to matter, and the weakest of the four.
queryCount 1, 10, 100, 1000 The strongest: the GPU’s fixed dispatch cost is amortised across the batch.
selectivity 0.001, 0.01, 0.1 High selectivity means many survivors, which means atomics contend and readback grows.
readback true / false A result the CPU reads costs a round trip that usually decides the answer on its own.
Repetitions 9, take the median Both sides have outliers; a mean lets one of them dominate.
Where the crossover sits at four batch sizes A bar chart in thousands of features showing the crossover point — the indexed-set size at which the GPU join becomes faster — for four query batch sizes. With a single query the crossover is beyond ten million features, meaning the CPU wins across the whole practical range. At ten queries it is about 3 200 thousand. At a hundred queries it is about 640 thousand. At a thousand queries it is about 180 thousand, so the GPU wins on quite modest datasets. CROSSOVER FEATURE COUNT · THOUSANDS 1 query CPU wins throughout 10 queries ≈3.2 M 100 queries ≈640 K 1000 queries ≈180 K 0 3500 7000 10500 K never very large sets only mid-sized sets modest sets Modelled proportions with results staying resident. Add a readback and every bar moves right.
Two orders of magnitude of movement from one parameter. That is why a threshold quoted in feature count alone, with no batch size attached, is not a usable rule.

Reading the surface

Sweeping four parameters produces a surface rather than a number, and three features of it are worth naming because they recur across applications.

Batch size dominates. Moving from one query to a thousand typically shifts the crossover by an order of magnitude in feature count, because the GPU’s fixed cost — the dispatch, the bind group, the pass — is paid once per batch and the CPU’s cost is paid once per query. An application with genuinely single queries will rarely find a GPU win at any size.

Readback is close to decisive. A result that stays on the GPU to feed a draw call costs nothing extra; a result the CPU reads in the same frame costs a synchronisation that is frequently larger than the join itself. The two cases are different enough that they should be swept as separate surfaces rather than as one parameter.

Feature count matters least. It is the axis people reach for first and it moves the answer least, because both implementations scale roughly linearly with it. What changes with size is not which side wins but by how much.

The practical output is not a single threshold but a rule with two or three conditions — “GPU when the batch exceeds fifty queries and the result stays resident” — which is a form the application can actually evaluate at runtime.

Turning the surface into a runtime rule

The harness produces a table; what the application needs is a predicate it can evaluate before dispatching. Getting from one to the other is a modelling step worth doing explicitly rather than by eyeballing.

The form that works is a small decision function with two or three conditions, derived from the surface rather than fitted to it. “Use the GPU when the batch exceeds fifty queries and the result stays resident” is a rule an application can evaluate in a line, and it captures the two parameters that actually move the answer.

Where the surface suggests a threshold that varies by device — and it will, by a factor of several — the rule should read the threshold from the capability record rather than hard-coding it. A conservative default, adjusted upward on devices that report large limits, gets most of the benefit without a per-device calibration run.

The last piece is a fallback for being wrong. A rule that picks the GPU on a device where it happens to lose costs a few milliseconds per frame and nothing else, which is a tolerable failure mode — but it is worth measuring the chosen path with a timestamp query occasionally and logging when the prediction disagrees with the measurement. That turns the rule into something that can be improved with data from real users rather than only from the harness.

Failure modes

  • The GPU always loses. Check for per-run pipeline creation in the harness. Compiling inside the timed region measures the compiler.
  • Results vary by more than 20 per cent between runs. Clocks are not at steady state. Discard the first run entirely and take the median of the rest.
  • The crossover moves between machines by a factor of five. Expected. The output of this exercise is a rule, evaluated per device against a capability record, not a constant baked into the source.
  • The GPU wins in the harness and loses in the application. The harness omitted the readback, or the application’s join is not batched the way the harness assumed.
  • The harness says GPU and the frame says otherwise. The join is winning in isolation and losing in context, because it now competes with the rest of the frame for the same device. Measure it with a timestamp query inside a real frame rather than alone.
  • Both sides return slightly different result sets. A boundary condition — inclusive against exclusive comparison at the edge of a query rectangle. Fix it before trusting any timing, because they are not measuring the same operation.
The three controls that make the numbers comparable Three controls applied before any timing is recorded. Both implementations are warmed — five runs for the CPU so the just-in-time compiler has optimised the hot loop, and a full pipeline build for the GPU so compilation is outside the measurement. Nine runs are taken and the median reported, because both sides produce occasional outliers that a mean would let decide the result. And both implementations are checked to return identical result sets first, because a boundary difference means they are not performing the same operation and the comparison is meaningless. BEFORE TIMING · THREE CONTROLS 1 Warm both sides JIT for CPU, compile for GPU or you are timing start-up 2 Median of nine discard the first run a mean lets one outlier decide 3 Assert equal results same set, exactly otherwise they differ in what they compute Half-open intervals on both sides is the convention worth standardising on.
The third control is the one that gets skipped and the one that invalidates everything else. An implementation that is faster because it silently drops boundary cases is not faster.

Backend / Python interop note

There is a third option the comparison usually omits, and on large static datasets it beats both: do the join on the server.

A spatial join between a fixed set of polygons and a large point set — points in administrative areas, sensors in catchments — has a result that does not depend on the camera. Computing it once in duckdb or geopandas and shipping the answer as a column removes the join from the client entirely.

python
import duckdb

def precompute_join(points_path: str, areas_path: str) -> None:
    duckdb.sql("INSTALL spatial; LOAD spatial;")
    duckdb.sql(f"""
        COPY (
          SELECT p.*, a.area_id
          FROM read_parquet('{points_path}') p
          LEFT JOIN read_parquet('{areas_path}') a
            ON ST_Within(p.geom, a.geom)
        ) TO 'joined.parquet' (FORMAT PARQUET)
    """)

The same reasoning applies to partially-static joins. A join whose polygon set is fixed and whose point set changes — sensors reporting into fixed catchments — can be precomputed per polygon and updated incrementally, so the client never performs a join at all and instead looks up a precomputed area identifier per point. That is a different shape of solution from either row of the comparison and it is frequently the right one, which is the argument for asking whether the join has to be dynamic before measuring how fast it can be made.

That option belongs in the harness as a third row, because a comparison that only weighs CPU against GPU will pick one of them for a problem that should not have been on the client at all. The client-side join earns its place when the query changes with the interaction — a lasso selection, a moving buffer distance, a filter the user is dragging — and those are exactly the cases where the batch is large and the result stays on the GPU.