週刊 spaghetti_source

matroid partition

16:12

今回はマトロイド分割，特に，k全域森問題を扱います．

一般論

台集合が共通しているマトロイド (V, I1), ..., (V, Ik) が与えられたとき，台集合 V の分割 {X1, ..., Xk}であって，各 Xj が Ij で独立になるようなものを求める問題をマトロイド分割（matroid partition）といいます．重み付きバージョンも自然に定義されます．マトロイド分割はマトロイド和（matroid sum）と呼ぶことも多いのですが，個人的には『分割』のほうが明確なので好きです．

マトロイド分割はマトロイド交差の特殊ケースです．つまり M1 と M2 への分割は，M1 と M2^* の交差になります．ここで M^* は M の双対マトロイドで，補集合が M で基のものを M^* で基と定義したものです．

マトロイド交差は2つのマトロイドまでしか効率的に解けなかったことに対し，マトロイド分割は任意個のマトロイドで効率的に解けます．これはマトロイド分割を次の独立割り当て問題とみなすことで理解できます（独立割り当ては前回を参照）：

• 左マトロイド：台集合 V の自由マトロイド
• 右マトロイド：与えられたマトロイドの直和．つまり V の K 個のコピーで，各コピー上で独立性が Ik で定義されたもの
• 枝：左のノード x から，対応する右の K 個の x のコピーに枝を引く

コピーの中から１つしか選べない（＝分割になっている）という条件を，二部グラフを使って表現している，ということです．

独立割り当ての特殊ケースなので，独立割り当てのアルゴリズムが適用できます．上から定義されるグラフを G := (L+R,A) とします．

(1) G の左から右への最短路を見つける（無ければ終了）
(2) 道にそってアップデートし，(1) に戻る

これをそのまま実装すると，グラフ G のサイズは L+R = O(kn)，A = O(kn^2) となり，重みなしで O(kn^3) 回の独立判定になります．

...

ここまでを前提として，次の論文を紹介します．非常に読みやすい論文です．

• J. Roskind and R. E. Tarjan (1985): "A Note on Finding Minimum-Cost Edge-Disjoint Spanning Trees", Mathematics of Operations Research, vol.10, pp.701-708;

この論文はk全域森問題（k disjoinit forest problem）：

重み付き単純グラフが与えられる．互いに素な森 k 個を選んで最小コストにせよ

という問題を扱ったものです．この問題はマトロイド分割の特殊ケース（グラフィックマトロイド k 個への分割）です．

論文ではマトロイド分割をマトロイド交差の特殊ケースと認識した上で，さらに，効率的に実装するにはどうしたらよいかを，この問題を例に説明しています．

1. matroid greedyとしての定式化

k全域森問題では，重みはマトロイドで共通しています．つまり，独立割り当てでいえば，重みは超ソースと左側のノードの間にのみあることになります．このとき一般の独立割り当てアルゴリズムのかわりに，matroid greedy的な実装を行うことができます：

F = (F1, ..., Fk) をすべて空集合で初期化
for e \in E: 重みの軽い順
  e を現在の解に追加できるなら追加

もちろん，この「追加できるなら」という部分でマトロイド交差的な交互道の問題を解く必要があるのですが，だいぶ見通しが良くなっています．

「追加できるなら」の部分を掘り下げると次のようになります．

1. e =: e(1) をどこかの Fj =: F(1) に入れる．代わりに抜かないといけない辺を e(2) とする
2. e(2) をどこかの Fj =: F(2) に入れる．代わりに抜かないといけない辺を e(3) とする
...
m. 何も抜かずにどこかの Fj =: F(m) に e(m) を追加できたら追加可能，できなければ追加不能

新しい辺 e で，玉突き式に，追い出していくイメージです．増加道探索（最短路）なので，同じ辺を２回チェックする必要が無いことに注意してください（∵独立割り当てのグラフでいうところの左側のノードを考える）．

2. cyclic scanning

augmentの手続きを簡単化します．上の実装では e(j) を入れるべき森 F(j) を毎回全部調べる必要がありましたが，実は F(j) = F_j にとれます．つまり F_1 から抜いた辺を F_2 に入れられるか試し，代わりに抜いた辺を F_3 に入れられるか試し，…… としても構いません．この手続きを cyclic scanning といいます．

この手続きで得られるパスは最短パスではないので，マトロイド交差の一般論だけでは実行可能性が保証されないパスになるのですが，実はこれでも大丈夫です（∵異なるマトロイドを跨る部分では，実行可能性が問題にならないため）．

cyclic scanning自体は計算量にも実装量にもたいして寄与しないのですが，次に述べる辺列挙の技術と組み合わさることで，大きな効果となります．

3. enumerating edges

実装で面倒くさいのは，辺 e(j) を森 F(j) に追加するために抜かないといけない辺 e(j+1) の列挙の部分です．

理屈の上では F(j) を e(j) の一方の端点で rootify し，残りの端点から根へのパスをたどることで実現できるのですが．このとおりに実装するには rootify 付きの Link-Cut Tree を作る必要があって，憂鬱な気持ちになります．

ところが，この操作は，実のところ，難しいデータ構造を使うこと無く，比較的簡単に実装できます．

辺 e = e(1) を新たに追加するステージを考えます（つまり，matroid greedyの i ステップ目）．このとき，すべての F_1, ..., F_m を e(1) の片側の頂点で rootify しておき，各 F(j) の根に「チェック済み」のラベルを張ります．この操作は普通のグラフ探索によって O(|F|) で実行できます．

e(1) を F(1) に追加する際には，片側は必ず根なので，根でない側を見ます．その頂点が F(1) に入っていなければ F(1) に追加できるので終了．F(1) に入っていれば，サイクルができたといことで，e(1) の代わりに抜く辺をすべてキューに突っ込みます．このとき「根に近い順に」突っ込みつつ，突っ込んだ辺の両端点を F(1) で「チェック済み」にしておきます．突っ込む順序が重要です．

次のステップで e(2) を F(2) に追加しようと試みます．キューに突っ込む順番から，e(2) の片側の頂点は必ず根です．なので F(1) の場合と同様の手続きを踏みます．

さらにその次のステップ；e(2)' を F(2) に追加できるかを見るステップが最重要ポイントです．e(2)' は両端点とも根ではありませんが，端点の少なくとも一方はチェック済みです（∵e(2) の根でなかった側がチェック済みになる）．よって未チェック側から根にさかのぼりつつ，チェック済み頂点にぶつかったら終了し，根に近い方から辺を追加する，という手続きでまだキューに突っ込まれていない交換枝をすべて追加することができます．

この議論を一般化すると，次の手続きで交換辺候補を効率的に列挙できることがわかります．

1. e(j) = (u(j), v(j)) のうち，少なくとも一方はチェック済みなのでそれを u(j) とする
2. v(j) から木 F(j) を根方向にさかのぼりつつチェックを付ける．
   チェック済み頂点が出たら終了
3. さかのぼる際に現れた辺を根方向から順にキューに追加する

このようにすると，同じ辺を2回チェックすることがなくなるため，増加道探索のトータル計算量が O(|F|) = O(nk) でバウンドできます．増加道探索の最初に rootify で O(|F|) 費やしましたが，それとバランスしていることに注意してください．

4. clump

現段階で計算量は O(mnk) になっています．すなわち O(m) 個の辺について，O(nk) で増加道チェックする計算量です．この，m の部分を高速化する概念として clump というものがあります（clumpは「茂み」の意）．

概念や正当性をきちんと説明するのが面倒なので，高速化した後の手続きのみを与えます．

n 頂点からなる union find を初期化
F = (F1, ..., Fk) をすべて空集合で初期化
for e = (u,v) \in E: 重みの軽い順
  if find(u,v) then e は追加不能
  e を現在の解に追加できるなら追加
  追加できなかった場合 unite(u,v)

union find で管理されている互いに素な集合 1 つ 1 つを clump といいます．同じ clump 同士を結ぶ辺が存在しないことを証明することで上記の手続きの正当性が証明されます．

※正当性は[Roskind-Tarjan 1985] を参照．また，clump は [Gabow-Westermann 1992] によって一般のマトロイド分割に一般化されているので，理解はこちらのほうが読みやすいかもしれません．

辺 e の augment をチェックするたびに，必ず |F| が増加するか，もしくは union find 内の互いに素な集合の数が１つ減ります．よって kn + n で augment が走る回数がバウンドでき，トータルで O(k^2 n^2) の計算量となります．

なお，現行最速はおそらく [Gabow-Westermann 1992] の O(k n^2 log k) なので，k を定数と思えば，このアルゴリズムが best possible です．

実装例

コード紛失 ＼(^o^)／

見つかったら貼ります；；

参考文献

• J. Roskind and R. E. Tarjan (1985): "A Note on Finding Minimum-Cost Edge-Disjoint Spanning Trees", Mathematics of Operations Research, vol.10, pp.701-708;
• H. N. Gabow, H. H. Westermann (1992): Forests, frames, and games: Algorithms for matroid sums and applications, Algorithmica, vol.7, no.1-6, pp.465-497.

追記

17:42

追記（2012/12/01）：みつかりました．k = 2，n ≦ 8 に対しては，全探索で動作確認していますが，自信はありません．

#define REP(i,n) for(int i=0;i<n;++i)
#define FOR(i,c) for(__typeof(c.begin())i=c.begin();i!=c.end();++i)
#define ALL(c) (c).begin(),(c).end()

struct UnionFind {
  vector<int> data;
  UnionFind(int size) : data(size, -1) { }
  bool unite(int x, int y) {
    x = root(x); y = root(y);
    if (x != y) {
      if (data[y] < data[x]) swap(x, y);
      data[x] += data[y]; data[y] = x;
    }
    return x != y;
  }
  bool find(int x, int y) { return root(x) == root(y); }
  int root(int x) { return data[x] < 0 ? x : data[x] = root(data[x]); }
  int size(int x) { return -data[root(x)]; }
};

struct Graph {
  int n, m;
  struct Edge {
    int s, t;
    int c;
    Edge(int s, int t, int c) : s(s), t(t), c(c) { };

    int id; // e in F_{e->id}
    Edge *prev; // e is an alternative for e->prev
  } NIL;
  vector<Edge> edge;
  void addEdge(int s, int t, int c) {
    edge.push_back(Edge(s, t, c));
  }
  Graph(int n) : n(n), m(0), NIL(-1,-1,0) { NIL.prev = &NIL; }
  typedef vector< set<Edge*> > Forest; // complexity increased by log n factor
  vector<Forest> forest;

  // pick connected component which contains x, and orient edges to x
  void orientTree(int u, Forest &F, vector<Edge *> &p) {
    FOR(i, F[u]) {
      Edge *e = (*i);
      if (e == p[u]) continue;
      if (e->s == u) swap(e->s, e->t);
      e->prev = 0; 
      p[e->s] = e;
      orientTree(e->s, F, p);
    }
  }
  bool augment(Edge *e0, int k) {
    int x = e0->s, y = e0->t;
    vector< vector<Edge *> > parent(k, vector<Edge*>(n));
    REP(i, k) { 
      parent[i][x] = &NIL;
      orientTree(x, forest[i], parent[i]);
    }
    queue<Edge*> Q; Q.push(e0);
    while (!Q.empty()) { // find augment path 
      Edge *e = Q.front(); Q.pop();
      int u = e->s, v = e->t, j = (e->id+1) % k;
      if (!parent[j][u] || !parent[j][v]) { // e can be added to F_j
        forest[j][u].insert(e);
        forest[j][v].insert(e);
        e->id = j;
        for (Edge *f; f = e->prev; e = f) { 
          f->id = (e->id+k-1)%k;
          forest[f->id][e->s].erase(e);
          forest[f->id][e->t].erase(e);
          forest[f->id][f->s].insert(f);
          forest[f->id][f->t].insert(f);
        }
        return true;
      } else { // e makes cycle -> enqueue circuit edges
        if (parent[j][u]->prev) swap(u, v);
        stack<Edge *> S; // near x edges first
        for (; !parent[j][u]->prev; u = parent[j][u]->t) S.push(parent[j][u]);
        for (; !S.empty(); S.pop()) { S.top()->prev = e; Q.push(S.top()); }
      }
    }
    return false;
  }
  static bool Compare(Edge i, Edge j) { return i.c > j.c; }
  int disjointForest(int k) {
    FOR(i, edge) { i->id = -1; i->prev = 0; }
    sort(ALL(edge), Compare);
    forest = vector<Forest>(k, Forest(n));
    int cost = 0;
    UnionFind clump(n); // clumps form disjoint sets
    FOR(i, edge) { 
      Edge &e = *i;
      if (clump.find(i->s, i->t)) continue;
      if (augment(&(*i), k)) cost += i->c;
      else clump.unite(i->s, i->t);
    }
    return cost;
  }
};