Spaghetti Source logo

最小シュタイナー木 (Dreyfus-Wagner)

説明

グラフ G と,部分頂点集合 T が与えられたとき,T をすべて張る木のことを T のシュタイナー木という.辺に重みのあるグラフにおいて,T のシュタイナー木のうちもっとも小さいものを求める問題を最小シュタイナー木問題,または単にシュタイナー木問題という.最小シュタイナー木問題は NP 困難である.

シュタイナー木問題のためには,理論的にはすべての木を試せばよいが,グラフのすべての木の数は最悪 2^{E} のオーダーになるため,全探索すら現実的でない.Dreyfus-Wagner は,動的計画法に基づくアルゴリズムで,うまく実装すれば時間計算量 O(n 3^t + n^2 2^t + n^3) を達成する.

方針は次のとおり.まず,Floyd Warshall などで全点対間最短路行列 d を計算しておく.次に,配列 OPT[S][p] を管理する.これは S のシュタイナー木で,p を端点として持つもののコストである.このとき,次の漸化式が成立する.

第一の式(ベースケース)は自明であろう.第二の式は,S を二つに分割し,それぞれの q を端とするシュタイナー木を計算し,それに p -- q という道をくっつけて p を端とする S のシュタイナー木を計算している.この配列が得られた後,全体 T のシュタイナー木は T を二つに分割し,それぞれの q を端とするシュタイナー木を求め,q を根として張り合わせることで得られる.以上で計算量は O(n^2 3^t + n 2^t + n^3) になる.

ところで,注意深く2ステップ目の計算を追うと,計算順序を変更することにより計算量を下げることができる.具体的には,これをさらに2ステップに分解する.

こうすると,ステップ 2 の計算はステップ 2-1 の計算をやった後にステップ 2-2 の計算を行っても同じ結果を与えることがわかる.これによって計算量は O(n^2 2^t + n 3^t + n^3) に下がる.

なお,Dreyfus-Wagner を改良したアルゴリズムとして Molle-Richter-Rossmanith がある.これは計算量を O( (2+δ)^t poly(n) ) に改良するが,ここではそれは実装していない.

使い方

const vector<int> &T
部分頂点集合.同じものが複数入っていても平気.
const matrix &g
隣接行列.無向性(対称性)を仮定する.つながっていない場合は inf.
戻り値
シュタイナー木の重み.ただし ≧ inf のときは存在しない(連結でない)ことをあらわす.シュタイナー木自体を求めたい場合は,Floyd Warshall の逆再生のように,中継する点を覚える配列を用意すればよい(はず).

計算量

|V| = n, |T| = t としたとき,

ただし以下の実装では E \subset S の生成で手を抜いているため,3^t が 4^t になっている.もちろん適当に実装すると 3^t にできるのだが,簡潔に実装すれば定数の小ささによって十分ペイすることが実測によりわかったため,こうしてある.

なお,この計算量は,大雑把に見て t = 11 程度が限界だと思ってよい.

ToDo

部分的な Verify は終わったが,完全な Verify は未だなのでどこかにあれば.

ソースコード

weight minimum_steiner_tree(const vector<int>& T, const matrix &g) {
  const int n = g.size();
  const int numT = T.size();
  if (numT <= 1) return 0;

  matrix d(g); // all-pair shortest
  for (int k = 0; k < n; ++k)
    for (int i = 0; i < n; ++i)
      for (int j = 0; j < n; ++j)
        d[i][j] = min( d[i][j], d[i][k] + d[k][j] );

  weight OPT[(1 << numT)][n];
  for (int S = 0; S < (1 << numT); ++S)
    for (int x = 0; x < n; ++x)
      OPT[S][x] = inf;

  for (int p = 0; p < numT; ++p) // trivial case
    for (int q = 0; q < n; ++q)
      OPT[1 << p][q] = d[T[p]][q];

  for (int S = 1; S < (1 << numT); ++S) { // DP step
    if (!(S & (S-1))) continue;
    for (int p = 0; p < n; ++p)
      for (int E = 0; E < S; ++E)
        if ((E | S) == S)
          OPT[S][p] = min( OPT[S][p], OPT[E][p] + OPT[S-E][p] );
    for (int p = 0; p < n; ++p)
      for (int q = 0; q < n; ++q)
        OPT[S][p] = min( OPT[S][p], OPT[S][q] + d[p][q] );
  }
  weight ans = inf;
  for (int S = 0; S < (1 << numT); ++S)
    for (int q = 0; q < n; ++q)
      ans = min(ans, OPT[S][q] + OPT[((1 << numT)-1)-S][q]);
}
/* 以下は PKU 3123 用
weight ans = inf;
for (int S = 0; S < (1 << numT); ++S) {
  weight sub = 0;
  for (int P = 0; P < 4; ++P) {
    int E = 0, p = -1;
    for (int k = 0; k < 8; k += 2) 
      if (((S >> k) & 3) == P)
        E += 3 << k, p = T[k];
    sub += (!!E) * OPT[E][p];
  }
  ans = min(ans, sub);
}
return ans;
*/

Verified

参考文献

前原 貴憲(maehara@prefield.com).

Last Modified: 2007.05.21 11:45:45.