Steiner Tree, in depth: the metric-closure 2-approximation, exact Dreyfus-Wagner dynamic programming, a worked example and geometric variants

The problem and why it is hard

Formally: given an undirected graph G = (V, E) with nonnegative edge weights and a terminal set T of size k, find a tree in G containing every terminal with minimum total weight. The tree may include any non-terminal vertices. Because weights are nonnegative an optimal connected subgraph can always be reduced to a tree, and every leaf of an optimal tree is a terminal, since a non-terminal leaf could be removed for free.

The two extremes show where the difficulty lives. With k = 2 the answer is a shortest path, solved by Dijkstra's algorithm. With k = |V| the answer is a minimum spanning tree, solved greedily by Prim's algorithm. In between, the hard part is choosing which Steiner points to use: once you fix the vertex set, the best tree on it is just its MST. Karp listed the decision version among his 1972 NP-complete problems, and it is also hard to approximate below a constant: no polynomial algorithm can guarantee a ratio better than 96/95 unless P = NP.

Two facts make it tractable in practice. The difficulty grows with the number of terminals, not the number of vertices, so instances with a dozen terminals in a large graph are solvable exactly. And a simple 2-approximation is fast and often much better than its guarantee.

A 2-approximation via the metric closure

The idea is to forget the Steiner points at first. Build the metric closure on the terminals: a complete graph whose edge between terminals a and b has the weight of the shortest a-to-b path in G. Take its MST, then turn each closure edge back into the path it stands for.

def steiner_approx(adj, terminals):
    sp = {t: dijkstra(adj, t) for t in terminals}             # (dist, prev) per terminal
    closure = mst_edges(terminals, lambda a, b: sp[a][0][b])  # 1-2: MST of the closure
    used = set()
    for a, b, _ in closure:                                   # 3: expand to real paths
        p = path(sp[a][1], b)
        used |= {frozenset(e) for e in zip(p, p[1:])}
    sub = subgraph(adj, used)
    tree = mst_edges(list(sub), lambda a, b: edge_weight(sub, a, b))  # remove cycles
    term = set(terminals)
    while True:                                               # 4: prune
        deg = degrees(tree)
        leaves = {x for x, d in deg.items() if d == 1 and x not in term}
        if not leaves:
            return sum(w for *_, w in tree), tree
        tree = [e for e in tree if e[0] not in leaves and e[1] not in leaves]

The helpers dijkstra, path, mst_edges, subgraph, edge_weight and degrees are the textbook routines their names suggest and are omitted for space.

Why at most twice optimal: take an optimal Steiner tree and walk around it, visiting each edge twice. That walk costs twice the optimum and passes through every terminal. Shortcutting it to visit terminals in order gives a cycle in the metric closure no more expensive than the walk, and dropping one edge of that cycle gives a spanning tree of the closure. The closure MST is no more expensive than that tree. A careful count gives the bound 2(1 - 1/l), where l is the number of leaves of the optimal tree. This is the algorithm of Kou, Markowsky and Berman (1981). The same doubling argument appears in the metric TSP approximation in approximation algorithms.

Exact: Dreyfus-Wagner dynamic programming

Dreyfus and Wagner's dynamic program computes the optimum by building trees over subsets of terminals. Let dp[S][v] be the cost of the cheapest tree connecting the terminals in S together with vertex v. Two moves build every such tree. Merge: a tree for S rooted at v may split at v into a tree for a subset S1 and a tree for S minus S1, both containing v. Grow: the meeting point can be moved along a path, so after merging, run Dijkstra on dp[S] seeded with every vertex's current value.

def dreyfus_wagner(adj, terminals):
    nodes, k = list(adj), len(terminals)
    full = (1 << k) - 1
    dp = [{v: INF for v in nodes} for _ in range(1 << k)]
    for i, t in enumerate(terminals):
        dp[1 << i][t] = 0
    for S in range(1, full + 1):                  # subsets in increasing order
        cur = dp[S]
        sub = (S - 1) & S
        while sub:                                # merge: every split of S, once
            if sub < (S ^ sub):
                left, right = dp[sub], dp[S ^ sub]
                for v in nodes:
                    cur[v] = min(cur[v], left[v] + right[v])
            sub = (sub - 1) & S
        pq = [(d, v) for v, d in cur.items() if d < INF]
        heapq.heapify(pq)                         # grow: multi-source Dijkstra
        while pq:
            d, u = heapq.heappop(pq)
            if d > cur[u]:
                continue
            for v, w in adj[u]:
                if d + w < cur[v]:
                    cur[v] = d + w
                    heapq.heappush(pq, (d + w, v))
    return dp[full][terminals[0]]

Enumerating every subset of every subset costs 3^k per vertex, and the Dijkstra pass costs O(m + n log n) per subset, so the total is O(3^k n + 2^k (m + n log n)) time and O(2^k n) memory. That is exponential only in k: twelve terminals means 531,441 times n merge operations, which is fine for graphs of thousands of vertices. The subset-iteration pattern is the same one used by Held-Karp in TSP via dynamic programming. To recover the tree as well as the cost, store for each dp[S][v] whether it came from a merge (and which split) or from a grow (and which predecessor), then unwind from dp[full][t].

Both functions were checked against brute force (enumerate every set of Steiner vertices and take the MST of the induced subgraph) on 300 random graphs with 6 to 14 vertices and 2 to 6 terminals: the dynamic program matched brute force every time, and the approximation always landed between the optimum and twice it, at worst 1.21 times.

Worked example: where the approximation loses

Four terminals A, B, C and D sit on the corners of a square whose sides weigh 5. A non-terminal hub s connects to each corner with weight 3.

Approximation. In the metric closure, adjacent corners are 5 apart (the side beats 3 + 3 through s) and opposite corners are 6 apart (through s). The MST of the closure picks three sides of weight 5, total 15, and expanding them changes nothing because each closure edge is a single graph edge. The hub is never used, because no single pair of terminals benefits from it enough.

Exact. The dynamic program finds dp[{A, C}] = 6, through s. Adding B costs only one more spoke: dp[{A, B, C}] = 9. All four give 12, the star through s. The ratio is 15 / 12 = 1.25, inside the bound 2(1 - 1/4) = 1.5 for an optimal tree with four leaves.

The lesson generalises: the approximation is blind to points that are only worth visiting for three or more terminals at once. Every improved algorithm works by considering such small groups explicitly.

Better approximations and practical solvers

In practice, large real instances are often solved exactly by combining graph reductions (removing vertices and edges provably absent from some optimal tree) with an integer programme over a directed cut formulation. When an exact answer is out of reach, a local search that starts from the 2-approximation, tries inserting and deleting Steiner points, and re-runs the MST on the chosen vertex set, usually closes most of the gap.

Euclidean and rectilinear variants

In the plane, Steiner points may be placed anywhere, not just at given vertices. In the Euclidean version every Steiner point has exactly three edges meeting at 120 degrees, and an optimal tree for k terminals has at most k - 2 Steiner points. The ratio between the MST and the Steiner tree is at most 2/sqrt(3), about 1.155, if the Gilbert-Pollak conjecture (Steiner ratio sqrt(3)/2) holds; a published proof was later found to have a gap, so treat the conjecture as open.

With rectilinear (Manhattan) distances, the version that matters for chip wiring, Hanan's theorem says an optimal tree can be found on the grid formed by drawing horizontal and vertical lines through every terminal. That turns a continuous problem into a graph Steiner problem on at most k squared grid points. Hwang showed the rectilinear MST is at most 3/2 times the rectilinear Steiner tree, which is why routers start from an MST and then improve it.

Failure modes

• Skipping the re-MST and pruning. Expanded shortest paths can overlap; summing closure weights double counts shared edges and keeps useless dangling paths.
• Disconnected terminals. If a terminal is unreachable, distances stay infinite and the code returns infinity or crashes. Check reachability first.
• Memory blow-up. The DP table is 2^k times n entries; twenty terminals in a million-vertex graph will not fit. Bound k before choosing the exact method.
• Negative weights. Dijkstra-based steps assume nonnegative weights; with negative edges the problem changes (adding edges can help) and these algorithms are wrong.
• Directed graphs. Multicast with asymmetric links is the directed Steiner tree problem, which is much harder to approximate; undirected algorithms give invalid trees.
• Ties and float sums. Comparing floating-point costs for equality in tests causes flaky results; compare with a tolerance or use integer weights.

Trade-offs: choosing a method

Choose by the number of terminals and how much the last few percent are worth. Under about 15 terminals, Dreyfus-Wagner gives the optimum directly. For larger terminal sets where cost matters, such as fibre build-outs, use reductions plus an integer programming solver. When speed matters more, as in recomputing multicast trees as group membership changes, use the metric-closure approximation, optionally followed by local search, and measure its gap against the exact method on sampled instances. For dynamic membership, incremental heuristics that graft a new terminal onto the existing tree by its shortest path avoid rebuilding at the cost of drifting from optimal; rebuild periodically. Kruskal's union-find approach, shown in the Kruskal MST walkthrough, is an equally valid choice for the MST steps.

What to do next

1. Implement the metric-closure approximation including expansion, re-MST and pruning.
2. Implement Dreyfus-Wagner with tree recovery, and reproduce the square example: 15 versus 12.
3. Cross-check both against brute force on small random graphs, as described above.
4. Measure the approximation's gap on samples of your real instances before trusting it.
5. If k is large and cost matters, try an integer programming solver with reductions.
6. For geometric problems, build the Hanan grid (rectilinear) or use a dedicated Euclidean solver rather than discretising by hand.
7. If links are asymmetric, treat it as a directed Steiner problem from the start.