Matrix Decomposition

The /sources_to_targets endpoint computes a many-to-many time/distance matrix between N sources and M targets. calculon-matrix does this by factoring the matrix so the expensive work is shared and the rest runs in parallel.

The algebra

The matrix is the Cartesian product of source set \(S=\{s_1,\dots,s_N\}\) and target set \(T=\{t_1,\dots,t_M\}\). It is separable along either axis:

Each cell is the shortest-path distance \(d(s_i, t_j)\). Computing it does not require an independent search per cell: a single forward search from \(s_i\) reaches every target, and a single reverse search from \(t_j\) reaches every source. So the whole matrix costs \(N + M\) searches, not \(N \cdot M\) - and the two factorings above cost the same. Factoring only decides which set is computed once and shared (frozen) versus iterated (looped).

A forward search from \(s\) and a reverse search from \(t\) that meet on a directed edge \(e\) compose to a candidate path:

where \(\overrightarrow{g}\) and \(\overleftarrow{g}\) are the forward/reverse shortest costs to \(e\). The sum is independent of evaluation order, so the searches may run in any sequence (or in parallel) without changing the result - this is what makes the decoupling safe.

Freeze the smaller side

if N ≤ M:  freeze SOURCES (forward), loop TARGETS (reverse probe the index)
else:      freeze TARGETS (reverse), loop SOURCES (forward probe the index)

The common factor (the frozen side) is computed once into a shared, read-only index; the looped side runs against it. Freeze the smaller set because:

• Parallelism - the parallel loop runs over the non-frozen set, so loop over the larger one for task count ≥ cores.
• Memory - the frozen index size grows with the frozen-set size.
• Degenerate case is optimal for free - for 1 × M, freezing the single source builds one forward one-to-many tree and the M targets become plain lookups into it (no M reverse searches). Symmetric for N × 1.

Two phases

graph LR
    A[Snap all unique points<br/>parallel] --> B[Phase 1: freeze min side<br/>build shared index, parallel]
    B --> C[Phase 2: loop max side<br/>probe index, parallel]
    C --> D[Assemble matrix]

Each search is a self-contained DirSearch (its own labels, queue, edge-status, heuristic). A search settling an edge reports it; the frozen phase records edge → (elem, cost, dist), the probe phase looks up its opposing edge in the frozen index and keeps the minimum connection per counterpart.

Because the meeting cost is order-independent, each cell takes the minimum over all meeting edges (not the first meeting by interleave order, as the old round-robin did). This is deterministic and slightly more accurate.

DirSearch

One directional search with self-contained state, parametrised by direction:

Forward expands a node's edges directly; Reverse expands opposing edges (using allowed_reverse and the opposing edge's cost). Both seed from snapped candidates and run to a distance threshold.

Parallelism & determinism

• Snapping (edge projection) for unique points runs with rayon.
• Phase 1 seeds + runs each frozen search in parallel, then merges their settled edges into the index.
• Phase 2 seeds + runs each probe search in parallel. Each probe element writes a disjoint set of result cells (its own row/column), so threads never touch the same cell.

Determinism: per-cell value is a min, which is commutative; within a single search it is computed sequentially. The result is therefore identical regardless of thread count - verified by a golden snapshot test (tests/golden_matrix.rs).

Results

100 × 100 grid (Montpellier metro, Monaco-scale tiles), Apple Silicon:

Roughly 10×, deterministic, and marginally closer to the routing engine.

Accuracy vs. the route engine

Matrix cells diverge from the point-to-point /route result by a modest margin in both directions: on some pairs the matrix is cheaper, on others the route is. Comparing extracted paths shows the cause is not a cost-accounting bug (the sequences differ in their middle edges, and distance - not just cost - changes), but that the matrix and the route engine are two independent approximate searches with different hierarchy limits, A heuristics, and pruning. Each is suboptimal on different inputs; the matrix frequently finds the shorter* path of the two.

Tightening this needs a ground-truth oracle (an exhaustive, heuristic-free Dijkstra) to decide which search is wrong on a given pair, then aligning the two on hierarchy/threshold behaviour - a separate accuracy effort, not a localized fix.