Contraction hierarchies (shortcut-based shortest paths)

Deep dive

How it works

Parts:

• Graph — static, positive weights (directed: in and out lists).
• Importance queue — min-heap of node priorities with lazy re-evaluation.
• Witness search — bounded local Dijkstra deciding whether a shortcut is needed.
• Upward graph — per node, edges (original and shortcut) to higher-ranked nodes; directed graphs add a reverse-upward graph.
• Middle node — stored per shortcut, for path unpacking.

1. Prioritise. Priority = edge difference (shortcuts a contraction adds minus edges it removes), plus terms for already contracted neighbours or the node's level (Bast et al. 2016).
2. Pop and re-check. Recompute the popped node's priority; if it grew, push it back (lazy update).
3. Contract v. For each in-neighbour u and out-neighbour w: c = c(u,v) + c(v,w). Search from u without v, stopping at distance c or a settled-node cap. If no path of cost ≤ c reaches w, add u→w with cost c and middle v. A missed witness only adds a needless shortcut.
4. Rank. v takes the next rank; its edges to remaining nodes become its upward edges.
5. Query. Dijkstra from s on upward edges, from t on reverse-upward edges. Stall-on-demand drops a node when a neighbour above it reaches it more cheaply.
6. Meet. Answer = min d_s(m) + d_t(m) over nodes both searches reached. They cannot stop at the first meeting: the path's top node is on both sides.
7. Unpack. Replace each shortcut u→w by u→v→w, recursively.

contract v:    u ──3── v ──4── w      u ────7──── w   (shortcut, only if no
 (rank low)                           u→w path of cost ≤ 7 avoids v)

query s → t:   rank 4        m ◉   best meeting node: min d_s(m) + d_t(m)
                            ↗   ↖
               rank 3     a ◉     ◉ b
                          ↗         ↖
               rank 1   s ◉           ◉ t     (only upward edges, both sides)

Code

Undirected grid with random weights; order by edge difference with lazy updates; bounded witness search; upward query with stalling. It compares 1,000 random pairs with plain Dijkstra. Distances only; unpacking needs a middle node per shortcut.

use std::{cmp::Reverse, collections::BinaryHeap, time::Instant};
type Adj = Vec<Vec<(usize, u32)>>; // undirected: each edge stored at both ends
const INF: u32 = u32::MAX / 2;

/// Dijkstra from s, stopping at `target`. With `stall` (upward graph only) a node is skipped if a
/// neighbour above it already offers a shorter way in. `d` must be all INF; returns (reached, popped).
fn search(g: &Adj, s: usize, target: usize, stall: bool, d: &mut [u32]) -> (Vec<usize>, usize) {
    let (mut seen, mut pops) = (vec![s], 0);
    let mut heap = BinaryHeap::from([Reverse((0u32, s))]);
    d[s] = 0;
    while let Some(Reverse((k, x))) = heap.pop() {
        if k > d[x] { continue; }
        pops += 1;
        if x == target { break; }
        if stall && g[x].iter().any(|&(y, w)| d[y] < INF && d[y] + w < k) { continue; }
        for &(y, w) in &g[x] {
            if k + w < d[y] {
                if d[y] == INF { seen.push(y); }
                d[y] = k + w;
                heap.push(Reverse((k + w, y)));
            }
        }
    }
    (seen, pops)
}

/// Shortcuts needed if v were contracted now: (u, w, cost) with u < w.
fn needed(adj: &Adj, done: &[bool], v: usize, d: &mut [u32]) -> Vec<(usize, usize, u32)> {
    let nb: Vec<(usize, u32)> = adj[v].iter().copied().filter(|&(x, _)| !done[x]).collect();
    let mut out = Vec::new();
    for &(u, cu) in &nb {
        let limit = nb.iter().filter(|&&(w, _)| w > u).map(|&(_, cw)| cu + cw).max().unwrap_or(0);
        let mut touched = vec![u];
        let mut heap = BinaryHeap::from([Reverse((0u32, u))]);
        d[u] = 0;
        while let Some(Reverse((k, x))) = heap.pop() { // witness search: avoids v, cost- and size-bounded
            if k > d[x] { continue; }
            if k > limit || touched.len() > 500 { break; }
            for &(y, c) in &adj[x] {
                if y == v || done[y] || k + c >= d[y] { continue; }
                if d[y] == INF { touched.push(y); }
                d[y] = k + c;
                heap.push(Reverse((k + c, y)));
            }
        }
        out.extend(nb.iter().filter(|&&(w, cw)| w > u && d[w] > cu + cw).map(|&(w, cw)| (u, w, cu + cw)));
        for x in touched { d[x] = INF; }
    }
    out
}

/// Contracts all nodes; returns the upward graph (edges to higher-ranked nodes) and the shortcut count.
fn build(mut adj: Adj) -> (Adj, usize) {
    let n = adj.len();
    let (mut done, mut d, mut nd) = (vec![false; n], vec![INF; n], vec![0i64; n]);
    let prio = |adj: &Adj, done: &[bool], d: &mut [u32], nd: &[i64], v: usize| {
        let deg = adj[v].iter().filter(|&&(x, _)| !done[x]).count() as i64;
        2 * (needed(adj, done, v, d).len() as i64 - deg) + nd[v] // edge difference + contracted neighbours
    };
    let mut heap: BinaryHeap<_> = (0..n).map(|v| Reverse((prio(&adj, &done, &mut d, &nd, v), v))).collect();
    let (mut up, mut count) = (vec![Vec::new(); n], 0);
    while let Some(Reverse((p, v))) = heap.pop() {
        let fresh = prio(&adj, &done, &mut d, &nd, v); // lazy update: re-check before contracting
        if fresh > p { heap.push(Reverse((fresh, v))); continue; }
        for (u, w, c) in needed(&adj, &done, v, &mut d) {
            count += 1;
            for (a, b) in [(u, w), (w, u)] {
                match adj[a].iter_mut().find(|e| e.0 == b) {
                    Some(e) => e.1 = c, // cheaper parallel edge
                    None => adj[a].push((b, c)),
                }
            }
        }
        up[v] = adj[v].iter().copied().filter(|&(x, _)| !done[x]).collect(); // all stay above v
        for &(x, _) in &up[v] { nd[x] += 1; }
        done[v] = true;
    }
    (up, count)
}
/// CH distance: two upward searches, best meeting node. `df`/`db` are scratch arrays of INF.
fn query(up: &Adj, s: usize, t: usize, df: &mut [u32], db: &mut [u32]) -> (u32, usize) {
    let ((sf, pf), (sb, pb)) = (search(up, s, usize::MAX, true, df), search(up, t, usize::MAX, true, db));
    let best = sf.iter().filter(|&&x| db[x] < INF).map(|&x| df[x] + db[x]).min().unwrap_or(INF);
    sf.iter().for_each(|&x| df[x] = INF);
    sb.iter().for_each(|&x| db[x] = INF);
    (best, pf + pb)
}

fn main() {
    let (w, n) = (150, 150 * 150); // grid with random weights 1..=100
    let mut x = 88172645463325252u64;
    let mut rnd = move || { x ^= x << 13; x ^= x >> 7; x ^= x << 17; x as usize };
    let mut adj: Adj = vec![Vec::new(); n];
    for v in 0..n {
        for u in [(v % w + 1 < w).then(|| v + 1), (v + w < n).then(|| v + w)].into_iter().flatten() {
            let c = 1 + (rnd() % 100) as u32;
            adj[v].push((u, c)); adj[u].push((v, c));
        }
    }
    let t0 = Instant::now();
    let (up, shortcuts) = build(adj.clone());
    println!("{n} nodes, {shortcuts} shortcuts, preprocessing {:?}", t0.elapsed());
    let pairs: Vec<_> = (0..1000).map(|_| (rnd() % n, rnd() % n)).collect();
    let (mut df, mut db) = (vec![INF; n], vec![INF; n]);
    let run = |name: &str, f: &mut dyn FnMut(usize, usize) -> (u32, usize)| {
        let (t0, mut pops) = (Instant::now(), 0);
        let r: Vec<u32> = pairs.iter().map(|&(s, t)| { let (d, p) = f(s, t); pops += p; d }).collect();
        println!("{name}: {:?}/query, {} popped", t0.elapsed() / 1000, pops / 1000);
        r
    };
    let fast = run("CH", &mut |s, t| query(&up, s, t, &mut df, &mut db));
    let slow = run("Dijkstra", &mut |s, t| {
        let (seen, p) = search(&adj, s, t, false, &mut df);
        let d = df[t];
        seen.iter().for_each(|&x| df[x] = INF);
        (d, p)
    });
    assert_eq!(fast, slow); // same distances
}

Complexity and constants

• Preprocessing: one witness search per neighbour pair per contraction, plus priority re-evaluations; the order sets its cost. Optimal ordering is NP-hard (Dibbelt et al. 2016), so real systems use heuristics.
• Query: the upward search space is the cost driver; 280 settled vertices on average in Western Europe (Bast et al. 2016, Table 1).
• With a distance metric and 12 cores the survey reports 12 min preprocessing, 858 scans and 0.87 ms; with compact turn expansion 2.27 ms.

Efficiency

• Cost: measured on this Mac (Apple M4 Pro), one thread, the code above on a 150 × 150 grid with weights 1-100 (22,500 nodes, 44,700 edges): preprocessing 0.65 s, 51,061 shortcuts; 14 µs and 215 popped nodes per query against 610-690 µs and 11,272 for Dijkstra with early stop. On 220 × 220: 1.48 s, 19 µs, 281 popped against 1.4 ms, 24,510. Bound: preprocessing is witness-search compute; queries are heap and cache-miss bound.
• Levers:
  • stall-on-demand — 2.2× here (measured: 493 → 215 pops, 31 → 14 µs) — every query
  • capped witness search — faster preprocessing for some extra shortcuts (Dibbelt et al. 2016 cap at 50 settled nodes) — large graphs
  • bucket-based many-to-many — about |S|+|T| searches instead of |S|·|T| (Bast et al. 2016) — distance matrices
  • PHAST — over 10× faster than Dijkstra for one-to-all on roads (Bast et al. 2016) — full distance trees
• Naive → standard → state of the art:
  • Dijkstra per query · 2.2 s on 18 M vertices (survey)
  • contraction hierarchies · 110 µs after 5 min
  • hub labels · 0.56 µs with 18.8 GiB (survey); customizable variants · new weights in under 1 s

Implementation notes

• Renumber nodes by rank and store the upward graph as CSR: both searches then read near-contiguous memory.
• Use a binary or 4-ary heap (d-ary heap) and u32 distances; reset only touched entries between queries.
• Store each shortcut's middle node; unpacking is linear in path length (Bast et al. 2016).
• Directed graphs: contract with in and out edges separately.

Pitfalls

• Stopping at the first meeting node gives wrong distances.
• The witness search must skip v, or it finds paths that no longer exist.
• Keep the cheaper edge when a shortcut duplicates an existing one.
• Uniform-cost grids have no obvious hierarchy (Storandt 2013); measure before adopting.
• Any weight change invalidates the hierarchy.

Variants and alternatives

In production

Road routing: OSRM (CH and a multi-level algorithm; its README recommends the latter by default and CH for very large distance matrices), GraphHopper (speed mode), RoutingKit. Studied elsewhere: game grid maps (Storandt 2013), pedestrian routing with indoor and outdoor sections (Li et al. 2021), public transport (adaptation described in the survey).

Implementations

Sources

• R. Geisberger, P. Sanders, D. Schultes, D. Delling, Contraction Hierarchies: Faster and Simpler Hierarchical Routing in Road Networks, WEA 2008. paper
• R. Geisberger, P. Sanders, D. Schultes, C. Vetter, Exact Routing in Large Road Networks Using Contraction Hierarchies, Transportation Science 46(3), 2012. paper
• H. Bast et al., Route Planning in Transportation Networks, in Algorithm Engineering (LNCS 9220), 2016. paper · arXiv
• J. Dibbelt, B. Strasser, D. Wagner, Customizable Contraction Hierarchies, ACM JEA 21, 2016. paper · arXiv
• D. Delling, A. Goldberg, A. Nowatzyk, R. Werneck, PHAST: Hardware-accelerated shortest path trees, JPDC 73(7), 2013. paper
• S. Storandt, Contraction Hierarchies on Grid Graphs, KI 2013. paper
• R. Li, S. Storandt, U. Müller, D. Weber, Barrier-Free Pedestrian Routing with Contraction Hierarchies, ACM SIGSPATIAL 2021. paper