NEAREST is quadratic (correlated LATERAL); use ASOF joins + a fast overlap candidate instead

[conradbzura]

Description

NEAREST scales quadratically. It transpiles to a correlated CROSS JOIN LATERAL (SELECT ... FROM b WHERE a.chrom = b.chrom ORDER BY ABS(distance) LIMIT k) — for every left row it scans and sorts all of b on that chromosome. Timings (1 thread), found by the benchmark suite:

N	GIQL NEAREST k=1	bioframe closest k=1
1e3	0.086 s	—
1e4	0.87 s	—
1e5	hangs (> 2 min)	0.083 s
1e6	hangs	~1 s

Expected Behavior

NEAREST should scale roughly log-linearly like a sorted sweep (as bioframe/bedtools do), not quadratically.

Root Cause

The per-left-row correlated ORDER BY distance LIMIT k subquery is an O(n · m · log m) nested loop.

Proposed fix (validated for k = 1)

DuckDB's ASOF JOIN finds the single nearest match on one side in log-linear time. A correct k = 1 nearest-interval combines three candidate sources, then ranks:

1. left: a ASOF JOIN b ON a.chrom = b.chrom AND a.start >= b.end (nearest interval ending at/before the left start),
2. right: a ASOF JOIN b ON a.chrom = b.chrom AND a.end <= b.start (nearest starting at/after the left end),
3. overlap (distance 0): a column-to-column INTERSECTS join — which already reaches DuckDB's IE_JOIN fast path via DuckDB INTERSECTS SEMI/ANTI joins fall back to BLOCKWISE_NL_JOIN instead of IE_JOIN (~80x slowdown) #208 (per-chromosome partitioning),

UNION ALL-ed and reduced with ROW_NUMBER() OVER (PARTITION BY <left key> ORDER BY ABS(distance), b.start, b.end) keeping rn = 1 (the same distance formula and tie-break the current LATERAL uses).

Validated end-to-end: exact row-for-row match with the current NEAREST at 1e4, and 1.51 s at 1e6 (runtime per-chromosome getvariable overlap form) vs the current quadratic hang. The ASOF core alone is 0.98 s at 1e6; the overlap candidate is the only piece that needs the #208 IEJoin treatment.

Scope

k = 1, unstranded, unsigned, no max_distance is the flagship case (bedtools closest / bioframe closest defaults) and is what this fix targets first. k > 1 needs the k-nearest-per-side (a windowed scan, not a single ASOF); strand/signed/max_distance add partitioning/filtering. Those shapes should decline to the current LATERAL plan until extended, mirroring how the INTERSECTS override declines shapes it does not yet handle.

[conradbzura]

Investigation write-up — overlap-materialization is an open R&D problem

Deferring the DuckDB NEAREST fast path to a future R&D task. The algorithm is validated and correct; the blocker is purely how to materialize the overlap intermediate without either (a) a bad DuckDB plan or (b) emitting DDL as a side effect of a transpile.

Note: the prototype used for the measurements below was not committed — no branch or PR captures it. This comment is the sole record. The design, SQL shapes, measurements, and open questions here are intended to be complete enough to re-derive the implementation when we pick this up.

Environment for all numbers below: Apple M4 Pro, DuckDB 1.5.4, synthetic hg38-shaped data, n = 1e6 per side, single-thread unless a thread count is stated, median of ≥3 reps.

The algorithm (works, exact-match correct)

For the flagship shape (k = 1, unstranded, unsigned, no max_distance, correlated, canonical registered base tables), k=1 nearest = the min-distance winner among three candidate sources, ranked with the same ROW_NUMBER() OVER (PARTITION BY <ref key> ORDER BY ABS(distance), start, end) tie-break the portable LATERAL/DataFusion-fallback forms use:

1. left: ref ASOF JOIN target ON ref.chrom = target.chrom AND ref.start >= target.end — nearest interval ending at/before the left start,
2. right: ref ASOF JOIN target ON ref.chrom = target.chrom AND ref.end <= target.start — nearest starting at/after the left end,
3. overlap (distance 0): a column-to-column range-overlap join.

This structurally mirrors the existing DataFusion _fallback_form decorrelation (a SELECT DISTINCT reference-key relation, per-key ranking, the surrounding CROSS JOIN LATERAL rewritten to a plain join on the reference key, and the reserved-column finalizer that hides the synthetic rk_*/rn columns from a surfacing SELECT */b.*). Verified exact row-for-row match (reference key, target row, and the (start,end) tie-break) against both the current LATERAL plan and bedtools closest on real DuckDB. The two ASOF branches are cheap (~0.25s each @1e6) and parallelize fine — they read ordinary base tables. The entire difficulty is the overlap branch.

Why the overlap is hard

The overlap join carries a chrom equality key (a.chrom = b.chrom AND a.start < b.end AND b.start < a.end). DuckDB only selects its fast IE_JOIN for INNER pure-inequality joins with no equality key (this is the #208 finding); with the chrom equality it plans HASH_JOIN + inequality residual, which is quadratic per chromosome:

overlap formulation	plan	time @1e6
plain chrom-equality range join (materialized inputs)	HASH_JOIN	82.8s
genome-linearized global coords (cid*span + pos, pure inequality)	PIECEWISE_MERGE_JOIN	>60s (killed)
integer chrom-id double-inequality (cid<=cid AND cid>=cid AND …)	range op	>60s (killed)
per-chromosome IE_JOIN (one WHERE chrom='chrN' subquery per chrom)	IE_JOIN	~0.4s ✅

Only per-chromosome partitioning reaches IE_JOIN. Linearization does remove the equality key and get a range operator, but the two inequalities aren't selective in 1-D global space (a.gstart < b.gend matches ~half the rows), so PIECEWISE_MERGE_JOIN degrades. GIQL doesn't know a table's chromosomes at transpile time, so per-chromosome partitioning requires runtime chromosome discovery — the #208 dynamic-SQL pattern: SET VARIABLE ov_sql = (string_agg of one WHERE chrom='<c>' subquery per chromosome); … query(getvariable('ov_sql')).

The DuckDB planner pathology

The per-chrom query(getvariable(...)) overlap on its own is fast (~0.4s to scan). But it's an opaque table function with no cardinality statistics, and the moment it sits in the candidate UNION ALL that feeds the ranking window, DuckDB mis-estimates it and re-executes the entire per-chromosome IEJoin per probe:

overlap read as…	time @1e6
bare SELECT COUNT(*) FROM query(getvariable(...)) (sole source)	~0.4s
window over query(getvariable(...)) alone (sole source)	~0.47s
inside the candidate UNION ALL feeding the ranking window	~105s

Everything that keeps the overlap inline re-triggers this:

• WITH ov AS MATERIALIZED (SELECT * FROM query(getvariable(...))) → still 101s (the MATERIALIZED hint does not give the optimizer a real cardinality or prevent re-execution),
• pre-reducing the overlap through its own window first, then unioning → still >90s.

The re-execution is of the SQL, so the only fix is to materialize the overlap into a concrete value the planner can't re-run.

The two materializations that work — and their trade-off

approach	@1e6 (1 thread)	@1e6 (8 threads)	side effect
temp table (CREATE OR REPLACE TEMP TABLE)	1.28s	0.41s	emits DDL from a transpiled SELECT
list variable (SET VARIABLE ov = (SELECT list(row) …) + UNNEST)	2.66s	1.58s	only SET VARIABLE (same mechanism as the INTERSECTS override)
bioframe closest (reference)	~1s	—	—

Both are exact-match correct. The tension:

• The temp table is fast and parallelizes well (3.1× from 1→8 threads), but a transpiler emitting CREATE … TABLE as a side effect of a SELECT is a real API-cleanliness problem (fails on read-only connections, surprising in a transaction). Note it is an internal per-query scratch for the computed overlap — not input-data materialization (which is the user's job and belongs in ingestion docs).
• The list variable avoids all DDL, but building one monolithic value (list() aggregates every row into a single value) and reading it (UNNEST of one list) are serialization points: it scales only 1.7× from 1→8 threads, so the gap to the temp table widens with parallelism (2.1× → 3.8×) and it never beats bioframe.

Multithread scaling detail:

threads	temp table	list var
1	1.28s	2.66s
2	0.82s	2.10s
4	0.55s	1.80s
8	0.41s	1.58s

Open R&D directions (not yet tried / not yet working)

• Make the list-var parallelizable — e.g. split into per-chromosome list values so UNNEST can fan out, or store parallel flat lists per column instead of one list-of-structs.
• Replace the ranking window with a GROUP BY … arg_min over the reference key, which may plan the getvariable differently than a window does (untested).
• Reduce the materialized overlap first — collapse to one best (min start,end) overlapping target per reference inside the getvariable, shrinking whatever we materialize (helps only when overlaps-per-reference > 1).
• A DuckDB cardinality hint / stats API for query() / table functions, if/when available, would let the inline getvariable overlap plan correctly with no materialization at all — the cleanest possible fix.
• A GIQL-level "prepared/registered" NEAREST where the overlap intermediate is built once via an explicit, documented setup step rather than smuggled into every transpiled query.

Related data-layout findings (context for ingestion docs)

Separately confirmed that the input Parquet layout matters a lot for the per-chromosome WHERE chrom='chrN' scans (INTERSECTS join, n=1e6, 1 thread), because DuckDB reads Parquet footer stats but our files currently have a single row group (nothing to skip):

input layout	time
VIEW over single-row-group Parquet (current E2E)	0.659s
VIEW, small row groups (50k rows)	0.389s
hive PARTITION_BY (chrom), VIEW	0.437s
materialized CREATE TABLE (native zonemaps)	0.306s
on-disk DuckDB file, 200 MB memory cap (out-of-core, keeps stats)	0.342s

Takeaways for a future "ingesting data for GIQL" doc: DuckDB already reads Parquet footer stats, but a single row group can't be skipped; smaller row_group_size or PARTITION_BY (chrom) unlocks row-group / partition pruning; and a persistent DuckDB file gives the full stats benefit out-of-core (doesn't require pulling everything into RAM). These help the E2E/view regime; the compute-only regime already materializes and benefits from the existing (chrom,start,end) sort.

Status

Algorithm and correctness approach are settled (exact-match verified against bedtools closest and the portable LATERAL plan); the only open question is the overlap-materialization strategy above. The prototype was reverted to keep the branch pristine, so this write-up is the only surviving artifact — nothing is captured on a branch or in a PR. Parking as R&D.