每日算法 — 使用java实现Steiner树:最小斯坦纳树与启发式近似算法

斯坦纳树(Steiner Tree)问题是图论中经典的NP-hard组合优化问题:给定无向带权图和终端节点集合,寻找连接所有终端节点的最小权重子图。与最小生成树(MST)不同,斯坦纳树允许引入额外的Steiner点来降低总代价。该问题在VLSI布线、网络设计、生物信息学等领域有广泛应用。本文用Java实现斯坦纳树的精确DP算法与启发式近似解法。

核心原理

问题定义

给定无向图G=(V, E, w)和终端节点集合T⊆V,斯坦纳树问题要求找到一棵包含T中所有节点的最小权重树。当|T|=2时退化为最短路径;当T=V时退化为最小生成树。

关键性质

  • Steiner点度数:最优解中每个Steiner点的度数至少为3。
  • 边数上界:n个终端的斯坦纳树最多包含2n-2个节点和2n-3条边。
  • NP-hard性:一般图上的斯坦纳树问题是NP-hard,但小规模实例可用动态规划精确求解。

Dreyfus-Wagner算法

对于终端数k较小的场景,使用状态压缩动态规划:
– 状态dp[mask][v]表示以mask描述的终端集合连通、且树以v为根的最小代价
– 转移:将mask拆分为两个非空子集sub和mask\sub,dp[mask][v] = min(dp[sub][v] + dp[mask\sub][v])
– 对每个状态运行Dijkstra在图上松弛边权

近似算法(Kou-Markowsky-Berman)

  1. 计算所有终端点对之间的最短路径
  2. 构造完全图( metric closure),边权为最短路径距离
  3. 求该完全图的最小生成树
  4. 将MST边展开为原图中的最短路径
  5. 在展开子图上求MST,即为2-近似解

完整Java实现

import java.util.*;

/**
 * 斯坦纳树(Steiner Tree)实现
 * 包含:Dreyfus-Wagner精确DP、KMB 2-近似算法、最短路径计算
 */
public class SteinerTree {

    private final int n;                  // 节点数
    private final List<Edge>[] graph;     // 邻接表
    private final int INF = Integer.MAX_VALUE / 4;

    public SteinerTree(int n) {
        this.n = n;
        this.graph = new ArrayList[n];
        for (int i = 0; i < n; i++) graph[i] = new ArrayList<>();
    }

    /**
     * 添加无向边
     */
    public void addEdge(int u, int v, int w) {
        graph[u].add(new Edge(v, w));
        graph[v].add(new Edge(u, w));
    }

    /**
     * 边类
     */
    static class Edge {
        int to, weight;
        Edge(int to, int weight) { this.to = to; this.weight = weight; }
    }

    /**
     * Dijkstra单源最短路径
     */
    public int[] dijkstra(int src) {
        int[] dist = new int[n];
        Arrays.fill(dist, INF);
        dist[src] = 0;
        PriorityQueue<int[]> pq = new PriorityQueue<>(Comparator.comparingInt(a -> a[1]));
        pq.offer(new int[]{src, 0});

        while (!pq.isEmpty()) {
            int[] cur = pq.poll();
            int u = cur[0], d = cur[1];
            if (d > dist[u]) continue;
            for (Edge e : graph[u]) {
                int v = e.to;
                if (dist[u] + e.weight < dist[v]) {
                    dist[v] = dist[u] + e.weight;
                    pq.offer(new int[]{v, dist[v]});
                }
            }
        }
        return dist;
    }

    /**
     * Floyd-Warshall全源最短路径
     */
    public int[][] allPairsShortestPath() {
        int[][] dist = new int[n][n];
        for (int i = 0; i < n; i++) {
            Arrays.fill(dist[i], INF);
            dist[i][i] = 0;
        }
        for (int u = 0; u < n; u++) {
            for (Edge e : graph[u]) {
                dist[u][e.to] = Math.min(dist[u][e.to], e.weight);
            }
        }
        for (int k = 0; k < n; k++) {
            for (int i = 0; i < n; i++) {
                for (int j = 0; j < n; j++) {
                    if (dist[i][k] + dist[k][j] < dist[i][j]) {
                        dist[i][j] = dist[i][k] + dist[k][j];
                    }
                }
            }
        }
        return dist;
    }

    /**
     * Dreyfus-Wagner精确算法:状态压缩DP求解斯坦纳树
     * 时间复杂度: O(3^k * n + 2^k * n * log n)
     * 适用于终端数k <= 10-12的情况
     */
    public Result dreyfusWagner(int[] terminals) {
        int k = terminals.length;
        int fullMask = (1 << k) - 1;

        // dp[mask][v]: 连接mask中所有终端、且以v为根的最小代价
        int[][] dp = new int[1 << k][n];
        for (int[] row : dp) Arrays.fill(row, INF);

        // 初始化:单个终端到自身的代价为0
        for (int i = 0; i < k; i++) {
            dp[1 << i][terminals[i]] = 0;
        }

        // 按mask大小递增处理
        for (int mask = 1; mask <= fullMask; mask++) {
            // 子集划分转移:mask = sub | (mask ^ sub)
            for (int sub = (mask - 1) & mask; sub > 0; sub = (sub - 1) & mask) {
                int other = mask ^ sub;
                if (other == 0) continue;
                // 确保不重复划分(sub < other)
                if (sub > other) continue;

                for (int v = 0; v < n; v++) {
                    if (dp[sub][v] + dp[other][v] < dp[mask][v]) {
                        dp[mask][v] = dp[sub][v] + dp[other][v];
                    }
                }
            }

            // 对当前mask运行Dijkstra-like松弛(SPFA风格)
            // 找到所有dp[mask][v] < INF的节点,在图上传播
            boolean[] inQueue = new boolean[n];
            Queue<Integer> queue = new LinkedList<>();
            for (int v = 0; v < n; v++) {
                if (dp[mask][v] < INF) {
                    queue.offer(v);
                    inQueue[v] = true;
                }
            }

            while (!queue.isEmpty()) {
                int u = queue.poll();
                inQueue[u] = false;
                for (Edge e : graph[u]) {
                    int v = e.to;
                    if (dp[mask][u] + e.weight < dp[mask][v]) {
                        dp[mask][v] = dp[mask][u] + e.weight;
                        if (!inQueue[v]) {
                            queue.offer(v);
                            inQueue[v] = true;
                        }
                    }
                }
            }
        }

        // 找到最优根节点
        int minCost = INF, bestRoot = -1;
        for (int v = 0; v < n; v++) {
            if (dp[fullMask][v] < minCost) {
                minCost = dp[fullMask][v];
                bestRoot = v;
            }
        }

        // 回溯重建树边
        Set<int[]> treeEdges = new HashSet<>();
        if (bestRoot != -1) {
            reconstructTree(dp, fullMask, bestRoot, terminals, treeEdges, new boolean[1 << k][n]);
        }

        return new Result(minCost, bestRoot, treeEdges);
    }

    /**
     * 递归回溯重建斯坦纳树的边
     */
    private void reconstructTree(int[][] dp, int mask, int root, int[] terminals,
                                  Set<int[]> edges, boolean[][] visited) {
        if (visited[mask][root]) return;
        visited[mask][root] = true;

        int k = terminals.length;

        // 检查是否是子集合并点
        boolean foundSplit = false;
        for (int sub = (mask - 1) & mask; sub > 0; sub = (sub - 1) & mask) {
            int other = mask ^ sub;
            if (other == 0 || sub > other) continue;

            if (dp[sub][root] + dp[other][root] == dp[mask][root]) {
                reconstructTree(dp, sub, root, terminals, edges, visited);
                reconstructTree(dp, other, root, terminals, edges, visited);
                foundSplit = true;
                break;
            }
        }

        if (!foundSplit) {
            // 是边上移:找到前驱节点
            for (Edge e : graph[root]) {
                int prev = e.to;
                if (dp[mask][prev] + e.weight == dp[mask][root]) {
                    int a = Math.min(root, prev);
                    int b = Math.max(root, prev);
                    edges.add(new int[]{a, b, e.weight});
                    reconstructTree(dp, mask, prev, terminals, edges, visited);
                    break;
                }
            }
        }
    }

    /**
     * Kou-Markowsky-Berman (KMB) 2-近似算法
     * 时间复杂度: O(k * n^2 + k^2 * log k)
     * 适用于大规模图,保证结果不超过最优解的2倍
     */
    public Result kmbApproximation(int[] terminals) {
        int k = terminals.length;

        // 步骤1:计算所有终端之间的最短路径
        int[][] allDist = allPairsShortestPath();

        // 步骤2:构造metric closure(完全图,边权为最短路径距离)
        // 步骤3:求完全图的最小生成树
        List<int[]> metricEdges = new ArrayList<>();
        for (int i = 0; i < k; i++) {
            for (int j = i + 1; j < k; j++) {
                metricEdges.add(new int[]{i, j, allDist[terminals[i]][terminals[j]]});
            }
        }
        metricEdges.sort(Comparator.comparingInt(e -> e[2]));

        // Kruskal求MST
        int[] parent = new int[k];
        for (int i = 0; i < k; i++) parent[i] = i;

        List<int[]> mstEdges = new ArrayList<>();
        for (int[] e : metricEdges) {
            int ri = find(parent, e[0]);
            int rj = find(parent, e[1]);
            if (ri != rj) {
                parent[ri] = rj;
                mstEdges.add(e);
            }
        }

        // 步骤4:展开MST边为原图中的最短路径,收集所有节点和边
        Set<Integer> steinerNodes = new HashSet<>();
        Set<String> steinerEdges = new HashSet<>();

        for (int[] me : mstEdges) {
            int u = terminals[me[0]];
            int v = terminals[me[1]];
            // 展开u到v的最短路径
            List<Integer> path = getShortestPath(u, v, allDist);
            for (int node : path) steinerNodes.add(node);
            for (int i = 0; i < path.size() - 1; i++) {
                int a = path.get(i), b = path.get(i + 1);
                int w = getEdgeWeight(a, b);
                steinerEdges.add(a + "#" + b + "#" + w);
            }
        }

        // 步骤5:在展开子图上求MST
        // 构建展开子图的邻接表
        Map<Integer, List<Edge>> subGraph = new HashMap<>();
        for (String es : steinerEdges) {
            String[] parts = es.split("#");
            int a = Integer.parseInt(parts[0]);
            int b = Integer.parseInt(parts[1]);
            int w = Integer.parseInt(parts[2]);
            subGraph.computeIfAbsent(a, x -> new ArrayList<>()).add(new Edge(b, w));
            subGraph.computeIfAbsent(b, x -> new ArrayList<>()).add(new Edge(a, w));
        }

        // Prim算法求MST
        int totalCost = 0;
        Set<Integer> mstNodes = new HashSet<>();
        Set<String> finalEdges = new HashSet<>();
        PriorityQueue<int[]> pq = new PriorityQueue<>(Comparator.comparingInt(e -> e[1]));

        int start = terminals[0];
        mstNodes.add(start);
        for (Edge e : subGraph.getOrDefault(start, Collections.emptyList())) {
            pq.offer(new int[]{e.to, e.weight, start});
        }

        while (!pq.isEmpty() && mstNodes.size() < steinerNodes.size()) {
            int[] cur = pq.poll();
            int v = cur[0], w = cur[1], u = cur[2];
            if (mstNodes.contains(v)) continue;
            mstNodes.add(v);
            totalCost += w;
            int a = Math.min(u, v), b = Math.max(u, v);
            finalEdges.add(a + "#" + b + "#" + w);

            for (Edge e : subGraph.getOrDefault(v, Collections.emptyList())) {
                if (!mstNodes.contains(e.to)) {
                    pq.offer(new int[]{e.to, e.weight, v});
                }
            }
        }

        // 转换边格式
        Set<int[]> edgeSet = new HashSet<>();
        for (String es : finalEdges) {
            String[] parts = es.split("#");
            edgeSet.add(new int[]{Integer.parseInt(parts[0]), 
                                  Integer.parseInt(parts[1]), 
                                  Integer.parseInt(parts[2])});
        }

        return new Result(totalCost, start, edgeSet);
    }

    /**
     * 获取a到b的最短路径(基于Floyd结果回溯)
     */
    private List<Integer> getShortestPath(int a, int b, int[][] dist) {
        // 简化为Dijkstra路径回溯
        int[] prev = new int[n];
        Arrays.fill(prev, -1);
        int[] d = new int[n];
        Arrays.fill(d, INF);
        d[a] = 0;
        PriorityQueue<int[]> pq = new PriorityQueue<>(Comparator.comparingInt(x -> x[1]));
        pq.offer(new int[]{a, 0});

        while (!pq.isEmpty()) {
            int[] cur = pq.poll();
            int u = cur[0];
            if (u == b) break;
            if (cur[1] > d[u]) continue;
            for (Edge e : graph[u]) {
                if (d[u] + e.weight < d[e.to]) {
                    d[e.to] = d[u] + e.weight;
                    prev[e.to] = u;
                    pq.offer(new int[]{e.to, d[e.to]});
                }
            }
        }

        List<Integer> path = new ArrayList<>();
        int cur = b;
        while (cur != -1) {
            path.add(cur);
            cur = prev[cur];
        }
        Collections.reverse(path);
        return path;
    }

    private int getEdgeWeight(int u, int v) {
        for (Edge e : graph[u]) {
            if (e.to == v) return e.weight;
        }
        return INF;
    }

    private int find(int[] parent, int x) {
        if (parent[x] != x) parent[x] = find(parent, parent[x]);
        return parent[x];
    }

    /**
     * 最小生成树(Prim算法)——用于对比
     */
    public int mstWeight() {
        boolean[] used = new boolean[n];
        int[] minEdge = new int[n];
        Arrays.fill(minEdge, INF);
        minEdge[0] = 0;
        int total = 0;

        for (int i = 0; i < n; i++) {
            int v = -1;
            for (int j = 0; j < n; j++) {
                if (!used[j] && (v == -1 || minEdge[j] < minEdge[v])) {
                    v = j;
                }
            }
            used[v] = true;
            total += minEdge[v];
            for (Edge e : graph[v]) {
                if (!used[e.to] && e.weight < minEdge[e.to]) {
                    minEdge[e.to] = e.weight;
                }
            }
        }
        return total;
    }

    /**
     * 结果封装
     */
    static class Result {
        int cost;                   // 总代价
        int root;                   // 根节点(参考)
        Set<int[]> edges;           // 树的边集

        Result(int cost, int root, Set<int[]> edges) {
            this.cost = cost;
            this.root = root;
            this.edges = edges;
        }

        @Override
        public String toString() {
            StringBuilder sb = new StringBuilder();
            sb.append("Steiner Tree Cost: ").append(cost).append("\n");
            sb.append("Edges (").append(edges.size()).append("):\n");
            for (int[] e : edges) {
                sb.append("  ").append(e[0]).append(" - ").append(e[1])
                  .append(" (w=").append(e[2]).append(")\n");
            }
            return sb.toString();
        }
    }

    /**
     * 主函数:演示斯坦纳树的精确解与近似解
     */
    public static void main(String[] args) {
        // 构造测试图:网格状拓扑
        SteinerTree st = new SteinerTree(12);

        // 构建图(0-1-2-3 横向,0-4-8 纵向等)
        int[][] edges = {
            {0,1,2}, {1,2,3}, {2,3,2},
            {0,4,4}, {1,5,1}, {2,6,3}, {3,7,2},
            {4,5,2}, {5,6,1}, {6,7,3},
            {4,8,2}, {5,9,4}, {6,10,1}, {7,11,2},
            {8,9,3}, {9,10,2}, {10,11,1}
        };
        for (int[] e : edges) st.addEdge(e[0], e[1], e[2]);

        // 终端节点:0, 3, 8, 11(四个角落)
        int[] terminals = {0, 3, 8, 11};

        System.out.println("=== 斯坦纳树问题演示 ===");
        System.out.println("图节点数: 12, 终端节点: " + Arrays.toString(terminals));
        System.out.println("全图MST权重: " + st.mstWeight());

        // 精确算法
        System.out.println("\n--- Dreyfus-Wagner 精确算法 ---");
        long start1 = System.currentTimeMillis();
        Result exact = st.dreyfusWagner(terminals);
        long time1 = System.currentTimeMillis() - start1;
        System.out.println("耗时: " + time1 + "ms");
        System.out.println(exact);

        // 2-近似算法
        System.out.println("\n--- KMB 2-近似算法 ---");
        long start2 = System.currentTimeMillis();
        Result approx = st.kmbApproximation(terminals);
        long time2 = System.currentTimeMillis() - start2;
        System.out.println("耗时: " + time2 + "ms");
        System.out.println(approx);

        // 对比
        System.out.println("\n=== 结果对比 ===");
        System.out.printf("精确解代价: %d%n", exact.cost);
        System.out.printf("近似解代价: %d%n", approx.cost);
        System.out.printf("近似比: %.2f%n", (double) approx.cost / exact.cost);

        // 测试更多终端的情况
        System.out.println("\n=== 更多终端测试 (6个终端) ===");
        int[] moreTerminals = {0, 2, 5, 7, 9, 11};
        Result exact2 = st.dreyfusWagner(moreTerminals);
        Result approx2 = st.kmbApproximation(moreTerminals);
        System.out.printf("6终端精确解: %d, 近似解: %d, 近似比: %.2f%n",
            exact2.cost, approx2.cost, (double) approx2.cost / exact2.cost);
    }
}

关键设计解析

状态压缩DP的子集枚举

Dreyfus-Wagner算法的核心在于for (int sub = (mask-1) & mask; sub > 0; sub = (sub-1) & mask)这一枚举技巧。它按降序枚举mask的所有非空真子集,每个mask被划分为两个互补子集sub和mask^sub。通过sub > other的剪枝避免重复划分,将时间复杂度控制在O(3^k)。

SPFA式图松弛

在DP状态转移后,需要通过图上的边松弛来传播代价改善。不同于对每次状态运行完整Dijkstra,这里使用SPFA(Shortest Path Faster Algorithm)风格:将所有已更新节点入队,逐层松弛邻接边。这是因为状态间的依赖具有局部性,SPFA在小规模图上效率更高。

Metric Closure的构造

KMB算法的精髓在于将原图”度量闭包化”:终端间的最短路径距离满足三角不等式,构成一个度量空间。在这个度量空间上的MST,展开回原始图后代价不会超过最优斯坦纳树的2倍。这一证明基于树度量到图度量的Lipschitz嵌入理论。

复杂度分析

算法 时间复杂度 空间复杂度 适用场景
Dreyfus-Wagner O(3^k · n + 2^k · m) O(2^k · n) k ≤ 12的精确解
KMB近似 O(k·n² + k²·log k) O(n²) 大规模图的快速近似
全对最短路径 O(n³)或O(k·(n+m)log n) O(n²) KMB前置步骤

其中k为终端数,n为节点数,m为边数。

扩展方向

  1. 终端数大于12的精确算法:结合整数线性规划(ILP)求解器处理中等规模实例。
  2. Steiner森林:多个不相交终端组分别连通,推广为Steiner森林问题。
  3. 有向斯坦纳树:处理有向图上的连通问题,近似比更高。
  4. ** prize-collecting变体**:允许不连接某些终端,但需支付惩罚代价。
  5. GPU并行:使用CUDA加速全对最短路径和DP状态计算。