bug: Bidirectional search returns a sub-optimal cost path

The bug is in include/cpp_common/bidirectional.hpp. Start by reading the found() function and the explore_forward/backward methods. Understand how best_cost is updated and the stopping condition. Run the provided SQL examples to reproduce the issue, then examine the algorithm's literature to verify the proposed fix. The fix involves updating best_cost for tentative costs and adjusting the stopping condition. Test changes with the reproduction script to ensure bdDijkstra returns optimal paths.

My apologies if the bug report is a bit verbose, I wanted to really clearly walk through my understanding of what happens currently vs what I think should happen based on my understanding of bidirectional search, in the hopes it more clearly justifies what I'm suggesting.

Platform/versions

SELECT version();
PostgreSQL 17.6 on aarch64-unknown-linux-gnu, compiled by gcc (GCC) 15.2.0, 64-bit

SELECT postgis_full_version();
POSTGIS="3.3.7 a0c7967" [EXTENSION] PGSQL="170" GEOS="3.14.1-CAPI-1.20.5" PROJ="9.7.1" LIBXML="2.15.1" LIBJSON="0.18" LIBPROTOBUF="1.5.2"WAGYU="0.5.0 (Internal)"

SELECT pgr_version();
3.4.1

The bug

pgr_bdDijkstra in about 10% of cases returns a path with higher cost than pgr_dijkstra does, given the same edge query and the same start/end pair. A simple reprex, which I'll give code for replicating at scale on real data at the end:

-- This returns 11
SELECT max(agg_cost) AS dijkstra_cost FROM pgr_dijkstra(
  $$SELECT id, source, target, cost, -1 AS reverse_cost FROM (VALUES
      (1, 1, 2,  6.0),  
      (2, 2, 5,  6.0),   
      (3, 1, 3,  7.0),  
      (4, 3, 5,  6.0),   
      (5, 1, 4, 10.0),   
      (6, 4, 5,  1.0)   
    ) AS e(id, source, target, cost)$$,
  1, 5, directed := true);

-- This returns 12
SELECT max(agg_cost) AS bd_cost FROM pgr_bdDijkstra(
  $$SELECT id, source, target, cost, -1 AS reverse_cost FROM (VALUES
      (1, 1, 2,  6.0),
      (2, 2, 5,  6.0),
      (3, 1, 3,  7.0),
      (4, 3, 5,  6.0),
      (5, 1, 4, 10.0),
      (6, 4, 5,  1.0)
    ) AS e(id, source, target, cost)$$,
  1, 5, directed := true);

Hypothesized cause

best_cost is updated on too narrow a set of nodes, and the stopping condition tests it incorrectly. In include/cpp_common/bidirectional.hpp:

bool found(const V &node) {
        /*
         * Update common node
         */
        if (forward_finished[node] && backward_finished[node]) {
            if (best_cost >= forward_cost[node] + backward_cost[node]) {
                v_min_node = node;
                best_cost =  forward_cost[node] + backward_cost[node];
                return false;
            } else {
                return true;
            }
        }
        return false;
    }

The code above only updates best_cost on doubly-settled nodes. But the literature on bidirectional search agrees the shortest path can meet at a node that is only on the frontier, i.e. still sitting in the priority queue. See
these notes - "The key result in proving correctness is that once a vertex u enters S, its d-value d[u] is equal to the true shortest path distance δ(s,u) - this is CLRS Theorem 24.6."

or, this StackOverflow thread where someone makes this same mistake then fixes it themselves.

Walking through the reprex, take nodes s, t, p, q, and r, with source node s and target t:

s --6--> q, q --6--> t    cost 12
s --7--> r, r --6--> t    cost 13
s --10-> p, p --1--> t    cost 11   <- the optimum

#	direction	pops	after exploring	best_cost
1	forward	(0, s)	s=0 is settled, q=6, r=7, p=10 sit in the queue	—
2	backward	(0, t)	t=0 settled, q=6, r=6, p=1 sit in the queue	—
3	backward	(1, p)	p=1 settled, s=11 (1 + 10) is added to the queue	—
4	forward	(6, q)	q=6 settled, t=12 (6 + 6) added to queue	—
5	backward	(6, q)	q=6 settled, found(q) sees both settled, 6+6=12	12
6	backward	(6, r)	r=6 settled	12
7	forward	(7, r)	r=7 settled, found(r) sees both settled, 7+6=13, so 12 >= 13 is false and we break	12

The optimum is sitting in the arrays for p from iteration 2 onward: forward_cost[p] = 10 and backward_cost[p] = 1, summing to 11. It is never added up, because found() requires both finished flags first, and the forward search terminates before it settles p.

The termination is arbitrary, too. The search only stops at iteration 7 because r with cost (7 + 6) happened to be the next node settled from both sides. That is just a coincidence of how the nodes pop for a given combination of weights, not a proof of the shortest path being found.

A possible fix

Disclaimer, I am not a C++ guy nor do I have a full picture of the pgRouting code. But I think, generally, best_cost needs to be updated from forward_cost[node] and backward_cost[node] directly, which hold the tentative cost as soon as we pop the prior node.

1. A new function update_mu

/* Update best cost if both searches have a cost for this node, track it in v_min_node */
void update_mu(const V &node) {
        if (forward_cost[node] == INF || backward_cost[node] == INF) return;
        if (forward_cost[node] + backward_cost[node] < best_cost) {
            best_cost = forward_cost[node] + backward_cost[node];
            v_min_node = node;
        }
    }

Then call it in explore_forward/backward when we relax each node's edges.

void explore_forward(const Cost_Vertex_pair &node) override {
   ...
         if (edge_cost + current_cost < forward_cost[next_node]) {
            forward_cost[next_node] = edge_cost + current_cost;
            forward_predecessor[next_node] = current_node;
            forward_edge[next_node] = graph[*out].id;
            forward_queue.push({forward_cost[next_node], next_node});
            /* update mu for tentative cost if applicable */
            update_mu(next_node);
         }
   }
   forward_finished[current_node] = true;
   /* update mu for settled cost if applicable */
   update_mu(current_node);
}

Everything cheaper than forward/backward_queue.top() is already settled, so any path we haven't already priced costs at least the sum of the two - this is what the literature around bidirecional search stopping conditions says (Pohl's stopping condition). Once that gets higher than best_cost, we know no further nodes will find a better path. Stop the while loop when this happens:

while (!forward_queue.empty() && !backward_queue.empty()) {
   if (forward_queue.top().first
      + backward_queue.top().first >= best_cost) break;

v_min_node is stored by update_mu, and forward_predecessor[p] is written when the edge is relaxed, so the path out of p in both directions exists when we exit. This means the existing reconstruction at bidirectional.hpp:163-185 should be unchanged and the bug fix hopefully is lean.

Going back to the reprex, best_cost gets updated to 11 on the second iteration:

#	forward_queue.top().first + backward_queue.top().first >= best_cost?	direction	pop	relax	best_cost
1	0 + 0 < INF	forward	s=0	q=6, r=7, p=10; no backward costs yet	—
2	6 + 0 < INF	backward	t=0	q: 6+6=12, r: 7+6=13, p: 10+1=11	11 at p
3	6 + 1 = 7 < 11	backward	p=1	s gives 1+10=11, no improvement on best_cost (it's the same path we already found just unidirectional)	11
4	6 + 6 = 12 >= 11	backward	break	11

Reproduction at scale

Perhaps optional for review, but with an edges table with real network data, a reviewer can run this through psql, altering the target query in edges to match their data:

\timing on

\if :{?edges}
\else
\set edges 'SELECT id, source, target, cost, reverse_cost FROM ways'
\endif
\if :{?n}
\else
\set n 500
\endif
\if :{?seed}
\else
\set seed .50
\endif

SELECT setseed(:seed) AS _seed \gset

DROP TABLE IF EXISTS public.bd_vs_dijkstra;

CREATE UNLOGGED TABLE public.bd_vs_dijkstra AS
WITH e AS (:edges),
starts AS (
  SELECT row_number() OVER () AS rn, vid
  FROM (SELECT source AS vid FROM e ORDER BY random() LIMIT :n) x
),
ends AS (
  SELECT row_number() OVER () AS rn, vid
  FROM (SELECT target AS vid FROM e ORDER BY random() LIMIT :n) y
)
SELECT
  s.vid AS source_vid,
  t.vid AS target_vid,
  (SELECT max(agg_cost) FROM pgr_dijkstra(:'edges',   s.vid, t.vid, directed := true)) AS dijkstra_cost,
  (SELECT max(agg_cost) FROM pgr_bdDijkstra(:'edges', s.vid, t.vid, directed := true)) AS bd_cost
FROM starts s
JOIN ends t ON t.rn = s.rn
WHERE s.vid <> t.vid;

\echo ''
\echo '=== how often does pgr_bdDijkstra return a worse path? ==='
SELECT
  count(*)                                                       AS pairs,
  -- Total n
  count(*) FILTER (WHERE dijkstra_cost IS NOT NULL
                     AND bd_cost IS NOT NULL)                     AS routed,
  -- bdDijkstra has high cost path than regular
  count(*) FILTER (WHERE bd_cost > dijkstra_cost)                 AS bd_worse,
  -- Sanity check regular dijkstra never is worse than bdDijkstra
  count(*) FILTER (WHERE bd_cost < dijkstra_cost)                 AS bd_cheaper_impossible,
  -- bd_worse / routed
  ROUND((100.0 * count(*) FILTER (WHERE bd_cost > dijkstra_cost))
        / (count(*) FILTER (WHERE dijkstra_cost IS NOT NULL
                                    AND bd_cost IS NOT NULL)), 1) AS pct_worse
FROM public.bd_vs_dijkstra;

(in my data)

pairs	routed	bd_worse	bd_cheaper_impossible	pct_worse
500	500	55	0	11.0

\echo ''
\echo '=== and by how much? ==='
SELECT
  ROUND(MIN(bd_cost / dijkstra_cost)::numeric, 3) AS min_ratio,
  ROUND(AVG(bd_cost / dijkstra_cost)::numeric, 3) AS avg_ratio,
  ROUND(MAX(bd_cost / dijkstra_cost)::numeric, 3) AS max_ratio
FROM public.bd_vs_dijkstra
WHERE bd_cost > dijkstra_cost;

min_ratio	avg_ratio	max_ratio
1.000	1.005	1.033