Bidirectional Dijkstra, in depth: the mu stop rule, a counterexample to the first-meeting rule, tested Python and path recovery

Why two searches settle fewer nodes

Why does meeting in the middle help? Picture a road network as roughly two-dimensional with uniform density. Dijkstra from s to a target at distance d settles everything within distance d, a disk of area proportional to d squared. Two searches that each reach about d/2 settle two disks of radius d/2, total area proportional to 2 x (d/2) squared, which is half. So the gain on road-like graphs is a constant factor of about two, not an asymptotic change.

The famous exponential gain belongs to a different setting. In a graph where the number of nodes within k hops grows like b to the power k, such as a tree with branching factor b or a puzzle state space, one search explores about b^d nodes while two explore about 2 b^(d/2). That is the argument behind meet in the middle. Keep the two cases apart when you estimate the speedup for your graph.

The algorithm, precisely

Keep two tentative distance maps, d_f from the source and d_b to the target, two priority queues, and a number mu, the length of the best complete s-t path seen so far, initially infinity. The backward search runs Dijkstra from t on the reverse graph: for a directed edge (x, y) of weight w, the backward search relaxes y to x. On an undirected graph both directions use the same adjacency.

Each step settles the minimum node of one queue and relaxes its edges as usual. Two rules make it correct:

1. Update mu on every relaxation. When the forward search relaxes edge (u, v) with weight w and v already has a backward label, set mu = min(mu, d_f[u] + w + d_b[v]), and symmetrically for the backward search. Updating mu only when a node is settled by both sides is the classic bug, because the best path may cross between the two balls on an edge whose endpoints are each settled by only one side.
2. Stop when top_f + top_b is at least mu, where top_f and top_b are the smallest keys in the two queues. Then return mu and the edge that produced it.

Which side to advance is a free choice that does not affect correctness. Strict alternation is simplest. Advancing the side with fewer queued entries tends to balance the work when one end sits in a dense city and the other on a sparse rural road.

Worked example: where the first meeting lies

Run the algorithm with strict alternation, forward first. The table records each settle.

After step 4 the smallest forward key is u at 4 and the smallest backward key is u at 4, so the sum 8 is at least mu = 7 and the search stops, returning 7 via edge a-b. No node has been settled by both sides. The intuitive rule would keep going: forward settles u at 4, backward settles u at 4, u is now settled by both, and the rule returns d_f[u] + d_b[u] = 8 through s-u-t, which is wrong. Step 2 also shows why relaxation-time updates matter: mu became finite before either search had settled u.

Implementation with path recovery

This implementation handles directed graphs given an adjacency list and its reverse, uses lazy deletion in binary heaps, checks mu on every relaxed edge, and rebuilds the path from the meeting edge. Heap costs are discussed in priority queues.

import heapq, math

def bidirectional_dijkstra(adj, radj, s, t):
    """adj[x] = [(y, w), ...]; radj is adj with every edge reversed. w >= 0."""
    if s == t:
        return 0, [s]
    dist = [{s: 0}, {t: 0}]
    parent = [{s: None}, {t: None}]
    heaps = [[(0, s)], [(0, t)]]
    done = [set(), set()]
    graph = [adj, radj]
    mu, meet = math.inf, None

    def top(i):
        h = heaps[i]
        while h and (h[0][1] in done[i] or h[0][0] > dist[i][h[0][1]]):
            heapq.heappop(h)            # discard stale entries
        return h[0][0] if h else math.inf

    side = 0
    while True:
        if top(0) + top(1) >= mu:
            break
        i = side if heaps[side] else 1 - side
        d, u = heapq.heappop(heaps[i])
        done[i].add(u)
        for v, w in graph[i].get(u, ()):
            nd = d + w
            if nd < dist[i].get(v, math.inf):
                dist[i][v] = nd
                parent[i][v] = u
                heapq.heappush(heaps[i], (nd, v))
            other = dist[1 - i].get(v)
            if other is not None and nd + other < mu:
                mu = nd + other
                meet = (u, v) if i == 0 else (v, u)   # original edge direction
        side = 1 - side
    if meet is None:
        return math.inf, None
    x, y = meet
    left = []
    while x is not None:
        left.append(x)
        x = parent[0][x]
    right = []
    while y is not None:
        right.append(y)
        y = parent[1][y]
    return mu, left[::-1] + right

Since top() cleans stale entries before the stop check, a non-empty heap always has a valid minimum when popped. If either queue empties, its top is infinity and the loop stops, which is correct: that side has settled everything it can reach, so every crossing edge has been checked.

Why the stop rule is correct

Why is mu optimal when top_f + top_b is at least mu? Suppose a shorter path P of length L below mu exists. Every node on P with d_f below top_f has been settled forward with its exact distance, and every node with d_b below top_b has been settled backward. Because L is below top_f + top_b, each node on P has either its true distance from s below top_f or its true distance to t below top_b; otherwise its two distances would sum to at least top_f + top_b, more than L. Walking along P from s, there is a last forward-settled node x followed by a node y that is backward-settled. Whichever of the two was settled later relaxed the edge between them and found the other's label already final, so mu was set to at most d_f[x] + w + d_b[y] = L. Thus mu is at most L, a contradiction.

That argument uses only non-negative weights, the same requirement as plain Dijkstra. With negative edges neither direction's settled labels are final, and the proof collapses.

Testing against plain Dijkstra

Test against plain Dijkstra on many small random graphs, including directed ones, unreachable pairs and zero weights. Random graphs find the meeting bug quickly, which is why the five-node example above is also worth keeping as a named unit test.

import random

def dijkstra(adj, s, t):
    dist, h = {s: 0}, [(0, s)]
    while h:
        d, u = heapq.heappop(h)
        if u == t: return d
        if d > dist[u]: continue
        for v, w in adj.get(u, ()):
            if d + w < dist.get(v, math.inf):
                dist[v] = d + w; heapq.heappush(h, (d + w, v))
    return math.inf

for trial in range(5000):
    n = random.randint(2, 12)
    adj, radj = {}, {}
    for _ in range(random.randint(0, 3 * n)):
        a, b, w = random.randrange(n), random.randrange(n), random.randint(0, 9)
        adj.setdefault(a, []).append((b, w)); radj.setdefault(b, []).append((a, w))
    s, t = random.randrange(n), random.randrange(n)
    got, path = bidirectional_dijkstra(adj, radj, s, t)
    assert got == dijkstra(adj, s, t)
    if path:
        assert path[0] == s and path[-1] == t

Also assert that the returned path's summed weight equals the distance, which catches reconstruction bugs that leave the number right and the route wrong.

Engineering and extensions

For production graphs with millions of nodes, replace the dictionaries with flat arrays indexed by node id and a per-query timestamp array, so a query does not pay to clear state. Store the reverse graph once in compressed sparse row form beside the forward graph. Bidirectional search also combines with goal direction: bidirectional A* needs consistent heuristics and its own stop rule, covered in A* variants. In contraction hierarchies both searches only move upward in the node ordering, and the stop rule changes again: each direction stops independently once its own smallest key is at least mu, because the two searches no longer meet in the middle of the path's length.

When you need more than one route, bidirectional search does not directly give alternatives; Yen's algorithm can use it as the inner shortest path call.

Trade-offs

Memory roughly doubles, since each query keeps two distance maps, two parent maps and two queues, and the graph needs a reverse adjacency for directed edges. In exchange the number of settled nodes falls, and settled nodes dominate run time because each one means heap operations and cache misses on the adjacency arrays.

The two searches are independent until they meet, so it is tempting to run them on two threads. This can help on large graphs, but the shared mu and the stop check need synchronisation, and the extra coordination often eats the gain for queries that finish in a millisecond. Measure before you parallelise. Bidirectional search also needs an explicit target: when you want the nearest of many facilities, either run a single search from the user or add a virtual super target connected to every facility with zero-weight edges and search backwards from it.

Failure modes

• Stopping at the first doubly settled node. Returns a non-optimal path, as in the worked example.
• Updating mu only at settle time. Misses paths that cross on an edge between the balls.
• Forgetting the reverse graph. On directed graphs the backward search must follow edges backwards.
• Stale heap tops. Reading the heap minimum without discarding stale entries makes the stop check use an outdated key and stop late or, with some bookkeeping bugs, too early.
• Wrong path join. Mixing up the edge orientation when the backward search found the meeting edge.

What to do next

1. Implement the code above and run the random test against plain Dijkstra for at least a few thousand graphs.
2. Add the five-node graph as a named regression test that would catch the first-meeting bug.
3. Count settled nodes for both versions on your real graph to measure the actual speedup.
4. Move to array-based state and a stored reverse graph before benchmarking at scale.
5. If you need more speed, try bidirectional A* with landmarks or a contraction hierarchy next.