Calculating Distance Matrices with SQL

Building an N×N pairwise distance matrix in DuckDB is an inherently $O(N^2)$ operation whose production failures come not from the algorithm but from the unbounded Cartesian product it materializes before any distance predicate is applied. This walkthrough sits under Vectorized Aggregations and isolates the exact bottleneck — a naive CROSS JOIN over a coordinate column that spills, thrashes NVMe, and OOM-kills the session — then gives a minimal session setup, a candidate-pruning two-stage rewrite, plan-validation queries, and fallback routing for invalid geometry and out-of-memory (OOM) conditions.

When DuckDB Spatial encounters an unfiltered self-join on a coordinate array, the planner expands the full product first and only then evaluates ST_Distance. For a million points that is $10^{12}$ pairs and roughly a terabyte of intermediate state — orders of magnitude past any sane memory_limit. The objective is to never construct rows that the distance threshold will immediately discard.

Root-Cause Analysis of Matrix Blowup #

Distance-matrix queries fail along four distinct axes, and each has a different fix. Treating an OOM crash as “needs more RAM” usually masks one of the first three.

• No spatial pre-filter. A bare CROSS JOIN forces an exact distance kernel on every pair. The fix is a cheap coordinate-range (bounding-box) test that runs before ST_Distance, so the planner discards distant pairs during the scan rather than after materializing them. This is the same envelope-pruning principle that drives R-tree-backed spatial indexing for join workloads.
• Predicate placement that defeats pruning. Wrapping the coordinate columns in ST_Point(...) and filtering only on the constructed geometry hides the cheap numeric comparison from the optimizer. The degree-space ABS(a.x - b.x) test must operate on the bare numeric columns so it can be pushed down ahead of geometry construction.
• CRS / unit mismatch. ST_Distance on EPSG:4326 inputs returns degrees, not metres, so a <= 50000 threshold silently means nothing. Both inputs must be projected to a metric CRS first; the rules for that live in the CRS mapping & transformations reference.
• Unbounded intermediate set. DuckDB decodes geometry — stored internally as WKB, as covered in ST_Geometry vs WKB — before applying predicates, so an unpruned product spikes the heap and spills. Whether the working set fits at all is governed by the in-memory vs disk storage boundaries for your dataset.

Deterministic Configuration #

Enforce strict resource boundaries before generating any matrix. DuckDB’s out-of-core executor handles spilling, but distance workloads need temp storage on a low-latency block device and join reordering enabled. Confirm the prerequisites first:

• spatial extension installed and loaded in the session
• Both coordinate sets share one CRS (a metric CRS if you need metre results)
• memory_limit set below total RAM with headroom for the candidate set
• A fast local temp_directory available for graceful spill

INSTALL spatial;
LOAD spatial;

-- Cap resident memory below total RAM: the candidate set, decoded geometry,
-- and result region coexist, so leave headroom to avoid an OS OOM-kill.
SET memory_limit = '6GB';

-- Match physical cores. Distance kernels are SIMD-bound; oversubscribing
-- hyperthread siblings contends for the same units and slows the scan.
SET threads = 4;

-- Spill to a fast local NVMe path, never a network mount — distance spills
-- are write-heavy and a slow device turns a spill into a stall.
SET temp_directory = '/mnt/fast_nvme/duckdb_spill';

-- Allow hash-based partitioning and join reordering; the matrix does not
-- depend on input order, so preserving it only blocks the optimizer.
SET preserve_insertion_order = false;

Boundary check: validate CRS alignment before execution. Mismatched projections silently corrupt distance output rather than raising an error. Project both inputs with ST_Transform to a metric CRS (EPSG:3857 or a local UTM zone) so ST_Distance returns metres.

Optimized Execution Pattern #

The fix replaces full Cartesian expansion with a two-stage filter: a cheap degree-space bounding-box pre-filter that prunes the overwhelming majority of pairs, then an exact ST_Distance only on the surviving candidates. Pushing the spatial predicate into the scan phase is the core technique shared across Modern Spatial SQL Query Patterns.

Before — the naive form materializes every pair before measuring:

-- ANTI-PATTERN: full N×N product built before any distance is evaluated.
SELECT a.id AS src_id, b.id AS tgt_id,
       ST_Distance(a.geom, b.geom) AS dist_meters
FROM points a
CROSS JOIN points b
WHERE a.id <> b.id
  AND ST_Distance(a.geom, b.geom) <= 50000;  -- runs on EVERY pair

After — the numeric pre-filter collapses the candidate set before geometry is touched:

WITH source_points AS (
    SELECT id, x, y, ST_Point(x, y) AS geom
    FROM staging.coordinates
    WHERE x BETWEEN -180 AND 180
      AND y BETWEEN -90 AND 90      -- drop out-of-range coords up front
),
filtered_candidates AS (
    SELECT
        a.id AS src_id,
        b.id AS tgt_id,
        ST_Distance(a.geom, b.geom) AS dist_meters
    FROM source_points a
    CROSS JOIN source_points b
    WHERE a.id <> b.id
      -- Cheap degree-space window on bare numeric columns; the optimizer
      -- pushes this down so ST_Distance never sees the far pairs.
      AND ABS(a.x - b.x) < 0.45     -- ±0.45° ≈ 50 km at mid-latitudes
      AND ABS(a.y - b.y) < 0.45
)
SELECT src_id, tgt_id, dist_meters
FROM filtered_candidates
WHERE dist_meters <= 50000;          -- exact threshold on survivors only

The behavioral change is entirely in when the comparison happens. The degree-space ABS() test is a coarse, fast gate on indexed-friendly numeric columns; the precise dist_meters <= 50000 check only ever runs on the few pairs that survive it. On clustered point sets this prunes 90%+ of pairs before any geometry kernel executes.

Trade-off: the 0.45° window is latitude-dependent — a degree of longitude shrinks toward the poles. For high-latitude data, project to a metric CRS first and pre-filter on projected metre coordinates, or tighten the window per latitude band rather than trusting one global constant.

Bounded aggregation over the candidate set #

Once the matrix is bounded, summarize it with columnar batch aggregation instead of holding the raw pairs. Chaining the threshold scan straight into a GROUP BY keeps the work inside vectorized aggregations and avoids re-materializing the full set:

SELECT
    CASE WHEN dist_meters <= 1000  THEN 'local'
         WHEN dist_meters <= 10000 THEN 'regional'
         ELSE 'distant' END AS proximity_tier,
    count(*)               AS pair_count,
    avg(dist_meters)       AS mean_dist,
    percentile_cont(0.5) WITHIN GROUP (ORDER BY dist_meters) AS median_dist
FROM filtered_candidates
WHERE dist_meters <= 50000
GROUP BY proximity_tier
ORDER BY proximity_tier;

For nearest-neighbour or top-K extraction, rank inside a partition instead of materializing the whole matrix — a pattern detailed under window functions for geospatial. ROW_NUMBER() over a per-source partition lets the planner stop after K rows rather than full-sorting every pair:

WITH ranked AS (
    SELECT src_id, tgt_id, dist_meters,
           ROW_NUMBER() OVER (PARTITION BY src_id ORDER BY dist_meters) AS nn_rank
    FROM filtered_candidates
)
SELECT src_id, tgt_id, dist_meters
FROM ranked
WHERE nn_rank <= 5;   -- 5 nearest neighbours per source point

Diagnostic Queries & Plan Validation #

Confirm the pre-filter is actually pushed down rather than applied after the join. Capture the plan and read the row estimate at each node:

EXPLAIN ANALYZE
SELECT src_id, tgt_id, dist_meters
FROM filtered_candidates
WHERE dist_meters <= 50000;

Diagnostic thresholds: in the EXPLAIN ANALYZE output, the row count entering the ST_Distance projection should be a small fraction of N². If the join node reports a cardinality near the full product, the ABS() predicate did not push down — check that the filter references bare numeric columns and not the constructed geometry. A healthy plan shows the FILTER node below the distance projection; an anti-pattern shows the distance computed first.

Track memory pressure and spill volume during a run with the duckdb_memory() table function:

-- One-liner monitor: which operators hold memory and how much spilled.
SELECT tag, memory_usage_bytes, temporary_storage_bytes
FROM duckdb_memory()
ORDER BY memory_usage_bytes DESC;

If temporary_storage_bytes climbs into the gigabytes for a query you expected to stay in memory, the candidate set is not being pruned — return to the plan and tighten the bounding-box window. If runtime degrades by more than ~15% across identical workloads after an engine upgrade or a data-distribution shift, re-capture the plan and re-run on a fresh connection to discard any session-level state before concluding it is a regression.

Geometry Validation & Fallback Routing #

Null or nonsensical distances trace back to malformed coordinates or invalid geometry entering the kernel. Gate construction with a range check and route failures to a quarantine table rather than letting them halt the pipeline:

WITH validated_points AS (
    SELECT id, x, y,
        CASE
            WHEN x BETWEEN -180 AND 180 AND y BETWEEN -90 AND 90
                 AND ST_IsValid(ST_Point(x, y))
                THEN ST_Point(x, y)
            ELSE NULL                 -- quarantine out-of-range / invalid
        END AS geom
    FROM staging.coordinates
)
SELECT a.id AS src_id, b.id AS tgt_id,
       ST_Distance(a.geom, b.geom) AS dist_meters
FROM validated_points a
CROSS JOIN validated_points b
WHERE a.id <> b.id
  AND a.geom IS NOT NULL
  AND b.geom IS NOT NULL
  AND ABS(a.x - b.x) < 0.45
  AND ABS(a.y - b.y) < 0.45;

For inputs that are polygons rather than points, repair self-intersections with ST_MakeValid before measuring, since ST_Distance on an invalid ring can return spurious zeros. Never allow an unvalidated cast into a production distance kernel.

When the candidate set still breaches memory_limit despite pruning — common with dense urban clusters where the bounding box no longer prunes — fall back to chunked execution. Process one batch of source points at a time and append, so the heap only ever holds one chunk’s worth of pairs:

-- Chunked matrix: bound the working set to one source batch at a time.
INSERT INTO distance_results
SELECT a.id AS src_id, b.id AS tgt_id,
       ST_Distance(a.geom, b.geom) AS dist_meters
FROM source_points a
CROSS JOIN source_points b
WHERE a.id BETWEEN ? AND ?          -- driver loops over id ranges
  AND a.id <> b.id
  AND ABS(a.x - b.x) < 0.45
  AND ABS(a.y - b.y) < 0.45
  AND ST_Distance(a.geom, b.geom) <= 50000;

Driving that loop from Python — sizing chunks against memory_limit and streaming each batch out — is the batch processing pipeline pattern; collecting the bounded result straight into a GeoDataFrame follows the DuckDB-to-GeoPandas sync workflow.

External Reference Standards. For strict geometry-type compliance, see the Open Geospatial Consortium Simple Feature Access specification. Kernel-specific optimization flags are documented in the DuckDB Spatial Extension overview.

Related #

See also

• Optimizing Point-in-Polygon Queries in DuckDB — the same two-stage pruning applied to containment joins.
• Using ST_ClusterDBSCAN for Spatial Grouping — distance-driven grouping once the matrix is bounded.

Up: Vectorized Aggregations · Modern Spatial SQL Query Patterns