週刊 spaghetti_source

Matrix Chain Multiplication / Optimal triangulation (approximation)

17:13

いつか理解して簡単化したいと思っていながら全然進まないMCMの話題です．

連鎖行列積問題（matric chain multiplication）とは，適当なサイズの長方行列 M1, M2, ... が与えられたとき，M1 x M2 x ... x Mn を計算するための最良の方法（カッコ付け）を求める問題です．ここで p x q 行列と q x r 行列の積は pqr で計算できるとします．たとえば M1 = 1x100, M2 = 100x100, M3 = 100x10 だった場合，(M1 M2) M3 だと 10000 + 1000 の計算コスト，M1 (M2 M3) だと 100000 + 1000 の計算コストになります．この場合は前者が最適解で，次善の解と比べて約10倍の差になります．

この問題は動的計画法の説明で用いられ，補足として「実は O(n log n) で解ける[Hu-Shing 1982 & 1984」とコメントされるのが様式美です．O(n log n) は比較モデルのもとで最善です [Bradford-Choppella-Rawlins 1993]．

本当はここで Hu-Shing の O(n log n) を紹介したいのですが，原論文がややこしくて，現状積ん読状態です．（誰か一緒に読みませんか？）……私見ですが，現在 O(n log n) を厳密に理解している人はほとんど居ないのではないか，と思っています（結果に触れる論文は多いが，内容について触れる論文が全く無い）．

そこで今回は近似アルゴリズムを紹介します．

連鎖行列積問題のままでは扱いにくいので，これを最適三角形分割（optimal triangulation）とよばれる次の問題に帰着します：正n角形の各頂点に値 w(1), w(2), ..., w(n) が対応しているものが与えられたとき，これの三角形分割であってコストが最小のものを求めよ．ここで三角形分割のコストは現れる三角形 (i,j,k) すべてについて f(i,j,k) := w(i) w(j) w(k) を足しあわせたものとして定義する．

連鎖行列積問題と最適三角形分割問題が同値であることを冒頭の例を変形することで確認しておきます．1x100, 100x100, 100x10 行列に対して重みが 1, 100, 100, 10 なる四角形を対応させます．四角形の三角形分割は２種類あり，それらが連鎖行列積の (M1 M2) M3 と M1 (M2 M3) に対応します．

さて，連鎖行列積問題の近似アルゴリズムは次のアイデアに基づきます：

• (1) とても横長の行列があれば，先に計算する．
  • 例： 3x4, 4x6, 6x10000, 10000x3, 3x2 なら最初に (6x10000) x (10000x3) を計算する．
• (2) そのようなものがなければ小さい方から貪欲に計算する．
  • 例： 2x4, 4x3, 3x5, 5x4 なら，2 が一番多く評価に現れるように (((2x4) x (4x3)) x (3x5)) x (5x4) とする．

これを最適三角形分割の言葉で言い直すと次のようになります：

• (1) とても重い点があれば，その両側を結ぶ．
• (2) そのようなものがなければ一番軽い点から扇状に分割する．

どの程度重ければ左右を直結してよいかが問題となりますが，次の定理を使います．

Theorem. 任意の添字 * について 1/w(*) + 1/w(i) ≦ 1/w(i-1) + 1/w(i+1) が成り立つとき (i-1) と (i+1) を直結する最適解が存在する．

証明：三角形分割を１つ固定し，その中の四角形 (i-1), i, (i+1), p を考える．三角形分割なので辺は (i-1)-(i+1) もしくは i-p である．前者であるための条件は，前者の分割コスト ≦ 後者の分割コストであり，整理すると 1/w(p) + 1/w(i) < 1/w(i-1) + 1/w(i+1) を得る．最適三角形分割ではどの頂点 p が i の対角点になるかわからないが，すべての頂点でこの条件が成り立っていれば十分．

上の不等式は，* として w(s) 最小の点を取るのが一番厳しいことは容易にわかります．これをあわせて，次のアルゴリズムを得ます[Chin 78, Hu-Shing 81]

s := argmin w(*)
while 不等式を満たす頂点 i が存在
  i-1, i+1 を結んで三角形 (i-1,i,i+1) を切り落とす
s を軸とする扇を出力

不等式を満たす i を探す部分の計算量がかかるように思えますが，実のところ，前から順番に見ていくだけで十分です．簡単のため多角形を回転して s = 1 にとります．S を「最終的に扇の羽になる頂点」の集合をあらわすスタックとすると次の手続きが得られます．

S = []
for i = 1, 2, ...
  while 1, Sの先頭2つ, i が不等式を満たす
    S の先頭要素を削除
  S に i を追加
1 を軸とする扇を出力

この手続きの誤差率は1/(2√3+3)≒0.1547...となります．これは，この手続きに対する最悪インスタンスを構成して最悪の値を取ったものです．計算量は同じ点が 2 回スタックに入らないので O(n)です．

※歴史的には，この手続きは Chin が（近似率の悪い）バージョンを提案し，それを見たHu and Shingがすぐに改善できることに気づいた，という流れになっています．

ランダム入力などでやってみると，かなりうまく合うことがわかります．以下の実装を用いた自分の実験では，n = 100, w[i] \in [100,150] uniform のランダム列に対して平均で近似率 1.005 程度の解が出ています．

これだけ簡単な手続きでこれだけ良い近似が出るのだから，厳密アルゴリズムでこれを使って枝刈りするとうまくいかないかなあ，と一時期思っていましたが，あまりうまく行かないモノですねえ．

実装例

typedef long long Int;
bool cond(Int a, Int b, Int c, Int d) {
  // 1/a + 1/c <= 1/b + 1/d  iff b-d separation
  return b*c*d + a*b*d <= a*c*d + a*b*c;
}
Int approxTriangulation(vector<Int> w) {
  int n = w.size(), k = 0, S[n];
  rotate(w.begin(), min_element(w.begin(), w.end()), w.end());
  Int cost = 0;
  for (int i = 1; i < n; ++i) {
    while (k >= 2 && cond(w[0], w[S[k-2]], w[S[k-1]], w[i])) {
      cost += w[S[k-2]] * w[S[k-1]] * w[i];
      --k;
    }
    S[k++] = i;
  }
  for (int i = 0; i+1 < k; ++i) 
    cost += w[0] * w[S[i]] * w[S[i+1]];
  return cost;
}

参考文献

• L. E. Deimel, Jr and T. A. Lampe (1979): "An invariance theorem concerning optimal computation of matrix chain products", North Carolina State University Technical Report #TR79-14.
• F. Y. Chin (1978): "An O(n) algorithm for determining a near optimal computation order of matrix chain products", Communications of the ACM, vol.21, no.7, vol.544-549.
• T. C. Hu and M. T. Shing (1981): "An O(n) algorithm to find a near optimal partition of a conve polygon", Journal of Algorithms, vol.2, pp.122-138.
• T. C. Hu and M. T. Shing (1982): "Computation of matrix chain products, Part I" SIAM Journal on Computing, vol.11, no.2, pp.362-373.
• T. C. Hu and M. T. Shing (1984): "Computation of matrix chain products, Part II" SIAM Journal on Computing, vol.13, no.2, pp.241-261.
• P. G. Bradford, V. Choppell, and G. J. E. Rawlins: "On Lower Bounds for the Matrix Chain Ordering Problem", Indiana University, Technical Report TR391.

Tractability of Packing/Covering/Partitioning Problems

17:15

前回は全域木を複数個取る問題を題材に，マトロイド分割の説明をしました．マトロイドに関する話はこれで一段落にして，今回は複数全域木問題の周辺問題について，解ける・解けないを論じた[Bernath-Kiraly 2011]を紹介します．今回はコードなしで，軽めの話です．

＃前回の，コード紛失したものは見つかりました → https://topcoder-g-hatena-ne-jp.jag-icpc.org/spaghetti_source/20121124/1354351349 ．いつもどおり自信が無いのでプロの実装を待っています．

グラフ G = (V, E) の辺に関する判定問題を考えます．辺集合の性質 A, B について，次の３つのクラスの問題を考えます．

• 詰込問題（Packing） A ∩ B：詰込 X, Y であって（i.e., X ∩ Y = φ），X が A, Y が B を満たすものがあるか．
• 被覆問題（Covering） A ∪ B：被覆 X, Y であって（i.e., X ∪ Y = E），X が A, Y が B を満たすものがあるか．
• 分割問題（Partitioning） A + B：分割 X, Y であって（i.e., 詰込かつ分割），X が A, Y が B を満たすものがあるか．

たとえば，前回の全域木を k（=2）個取る問題は，全域木と全域木の詰め込みで，Spanning Tree ∩ Spanning Tree となります．重み（位数）の最適化と，存在判定にはギャップがありますが，少なくとも存在判定が解けない問題は最適化問題もNP-hardです．

結果

以下の７個の条件の組合せを考えます．

• P: パス
• Pst: s-t パス
• C: サイクル
• T: 木
• SpT: 全域木
• F: 森
• Cut: カット（i.e., 取り除くと非連結）

分割問題

Theorem. F+F, F+SpT, SpT+SpT, Cut+SpT, Cut+Cut 以外はNP-Complete

F/SpT系はマトロイド分割です．Cut+SpTは取れません（Cutを除くとSpTが存在しない）．Cut+Cutは，こう分割できることがbipartiteであることと同値であることから従います（可能側は自明．不可能側はbipartiteでないとき奇サイクルが存在することを使う）．

注目すべきは SpT+SpTがPで，T+T がNP-Completeになることでしょうか．SpTは極大Fとして特徴づけられるので，Tの難しさ（連結性）を正面から考えなくて済んでいる，ということだと思います．

被覆問題

被覆では，T, F, SpT は，重なっている部分を取り除くことで F との分割問題に帰着します．つまり P, Pst, C, Cut の４つの組合せが考慮対象になります．

Theorem. すべての組合せがNP-Complete

カバーは，難しい．

詰込問題

詰込では，P, T, F が絡むものは自明にすればよいので，そうでない Pst, C, SpT, Cut の４つの組合せが考慮対象になります．

Theorem. Pst∩SpT, C∩SpT 以外は P

Cut絡みは，基本的に自明です．Cut∩Cは小さくサイクルを取ってそれをカットの片側にするように頂点を分割します．Cut∩PstとCut∩Cutも同じノリです．SpT∩SpTはマトロイド分割です．Cut∩SpTが取れないのは分割と同じです．

Pst絡みですが，Pst+Ps't'は辺素パス問題のk=2の場合なので，グラフマイナータイプのアルゴリズムがあります [Robertson-Seymour 1995]．k=2の特殊ケースでは，グラフマイナー理論の確立よりも前に[Seymour 1980; Thomassen 1980; Shiloach 1980]でも論じられています．

Pst∩C, C∩Cは論文でアルゴリズムが提示されており，どちらも幅優先探索に基づきます．それほど実装は重くなさそうなので，気が向いたら実装するかもしれません．

参考文献

• [1] A. Bernath and Z. Kiraly (2011): "On the tractability of some natural packing, covering and partitioning problems", Technical Report.
• [2] N. Robertson and P. D. Seymour (1995): "Graph Minors XIII. The Disjoint Paths Problem", Journal of Combinatorial Theory, Series B, vol.63, no.1, pp.65-110.
• [3] P.D. Seymour (1980): "Disjoint paths in graphs", Discrete Mathematics, vol. 29, pp. 293-309.
• [4] Y. Shiloach (1980): "A polynomial solution to the undirected two paths problem", Journal on ACM, vol. 27, pp. 445-456.
• [5] C. Thomassen (1980): "2-linked graph", European Journal of Combinatorics, vol. 1, pp. 371-378.

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;
  }
};

Matroid Intersection

15:17

今回はマトロイド交差です．複雑といわれますが，フレームワーク自体は簡単です．

定義

マトロイド (V, I1) と (V, I2) に対して，I1 と I2 両方で独立な最大サイズの X を求める問題をマトロイド交差（matroid intersection）といいます．重みがある場合には，区別して，重み付きマトロイド交差（weighted matroid intersection）といいますが，ここでは特に区別しません．

重要な性質として「両方のマトロイドで独立な集合」がマトロイドの公理を満たさないことと，それにもかかわらず，マトロイド交差問題そのものは多項式時間で解けることがあります．

マトロイド交差では，より一般化した「独立割り当て問題（independent assignment problem）」を考えたほうがわかりやすいと思うので，そちらを考えることにします[Iri-Tomizawa 1976]．この問題は次のように定義されます（重み付きver）．

• 入力：マトロイド (V1,I1,w1), (V2,I2,w2)，これらを左右の頂点集合とする重み付き二部グラフ (W,A,w)
• 出力：二部グラフマッチング M であって，左端点の重み + 右端点の重み + 二部グラフの重みを最大化するもの．i.e., max. w1(M1) + w(M2) + w(M) ... （M1 は左端点，M2 は右端点）

通常のマトロイド交差は左右の頂点集合がイコールで，二部グラフが左右を一対一に結ぶケースです．独立割り当てに一般化しても，問題の難度は変わりません．

独立割り当て問題で，左右のマトロイドが自由マトロイドのときは，普通の二部グラフマッチングになります．つまり，独立割り当ては「二部グラフマッチングに毛が生えた」程度の問題です．実はアルゴリズムも「二部グラフマッチングに毛が生えた」程度です．

※３つ以上のマトロイド交差も対応する割り当て問題が作れて「『３次元マッチング』に毛が生えた」問題と対応します．無毛の３次元マッチングがNP-Hardなので，３マトロイド交差もNP-Hardとなります．逆に，Pで解ける３次元マッチングの部分クラスに対応するような３マトロイド交差は，解ける見込みがあります．

応用例

適当に２つマトロイドをもってきて交差をとれば，いくらでも気持ち悪い格好良い問題が生成できます．

• 有向全域木＝グラフ的＋横断：全域木であって各頂点に入ってくる枝が高々１本であるようなもの．
• カラフル全域木＝グラフ的＋分割：枝に色がついているグラフにおいて，「同じ色の辺が高々 c 個までしか無い」的な全域木．
• 非巡回向き付け＝横断＋分割：無向グラフの辺を「前後向きの二倍に増やし，同じ色に塗る（全く別の枝とは別の色にする）」と，「同じ色の枝は高々１本」＋「入次数は高々１」でグラフの非巡回向き付けになる．

あたりは別名が有名な問題です．他にも

• 分割＋分割：たとえば：人の集合を「人種で組み分け」したものと「宗教で組み分け」したものが与えられたときに「同じ人種の人は２人まで」かつ「同じ宗教の人は２人まで」を満たす最大部分集合を求める問題．
• グラフ的＋補グラフ的：全域木であって，取り除いたグラフが連結であるようなもの．
• グラフ的＋ガンモ：グラフ G の部分集合 H と特別な頂点 s が指定されたとき，H 内の森で，s から各点に辺素パスが引けるもの
• 分割＋スケジューリング：単位時間仕事で納期・報酬のほかに色がついていて，同じ色の仕事は k 個まで的問題．
• グラフ的＋スケジューリング：グラフの枝に納期が設定されていて・・・（略）
• etc

などが作れます．

分割＋分割がマトロイド交差，つまり二部グラフマッチング，つまり最小費用流に帰着することは，マトロイドということを意識せずに使ったことのある人も多いかと思います．

アルゴリズム

マトロイド交差のアルゴリズムは以下のようになります．

while true:
  入力と現在の実行可能解 M をもとに二部グラフマッチングと同じようなグラフ G を作る．
  s-t 最短路を求め，それにそって実行可能解を更新する．最短路がない場合は現在の解が最適．

これも二部グラフの最小重みマッチングと同じフレームワークです．違いは作るグラフ（と，更新手続き）．概要は以下のとおり：

• 頂点集合は左から {s}, V1, V2, {t}
• 枝集合は 5 種類（詳細後述）：
  • s-V1：左マトロイドの増加公理に対応
  • V1-V1：左マトロイドの交換公理に対応
  • V1-V2：入力される二部グラフ
  • V2-V2：右マトロイドの交換公理に対応
  • V2-t：右マトロイドの増加公理に対応

二部グラフマッチングとの差異は，V1-V1枝と，V2-V2枝があること．これらはジャンプアーク（jump arc）や交換枝（exchange arc）といわれます．口頭では「ヒゲ」ということもあります．各枝の具体的な定義は

• s-V1：x ∈ V1-M1 について M1 + {x} が実行可能なら s → x を引く
• V1-V2：マッチングに入っていないものは左から右に，入っているものは右から左に引く
• V2-V2 枝： x ∈ V2-M2 かつ y ∈ M1 について M1+{y} が実行不能かつ M1-{x}+{y} が実行可能なら x → y を引く

です．書いていないところは対称的なので省略します．このグラフで s-t パスにそって増加すると，サイズの1つ大きな交差が得られます．これは例えば

• s → v1 → v2 → u2 → u1 → w2 → t というパスを考えると（数字は左右に対応），対応する操作は
• 左のマトロイド条件に適合している要素 v1 を追加して，マッチングで対応する v2 も追加する．v2 を追加することで右のマトロイド条件が崩れるから u2 を抜いて，マッチングで対応していた u1 を新たに w2 につなぎ替えて，ゴール

となることから，そうかな，という気がしてきます．本気で証明しようとすると，実行可能性が不安なケースが１つ発生するのですが，最短パスであるということから，大丈夫だと証明できます（正確には：優先度付きキューによる最短路なら大丈夫．任意パスではダメ）．

悲しいお知らせ

計算量を見積もりましょう．

作られるグラフの枝数は，交換枝が全部で O(n^2) なので，普通にやると O(n^2) 回のオラクルコールと最短路アルゴリズムの時間 O(n^2 log n) になります．augment の総数をかけて，全体で O(n^3) オラクルコール＋O(n^3 log n)時間．※最短路でDijkstraが使えることは，二部グラフマッチングと同じくポテンシャルが取れることに基づきます）．

実のところ，一般論ではこのあたりが限界だと思われています．現行最速のマトロイド交差は [Cunningham 1986] の O(n r^{1.5}) = O(n^{2.5}) 回のオラクルコールですが，長い間更新がありません．※線形マトロイドに限ると [Harvey 2009] の O(n r^{ω-1}) （ωは行列積の指数）です．

ところが，多くのマトロイド交差に落ちる問題は，その問題独自のずっと効率的な解法が知られています．たとえば有向全域木問題は Chu-Liu/Edmonds の O(E log V) がありますが，マトロイド交差の一般論では作るグラフがそもそも O(E^2) サイズなので，O(E^3) とかそんな計算量になります．

要するに：競技プログラミングなどの「応用問題」では，マトロイド交差一般形をそのまま使うのは極めて稀です．多項式性の基準にしたり，アルゴリズムの着想点にするくらいに思っておいて，普通は問題ありません．多分．

交差好きなのになあ・・・

実装例

次の問題を考えます：

• 単位時間・締め切り付きの仕事が与えられる．
• 各仕事には「発注者」のタグがついていて，同じ発注者からの仕事は2つまでしかできない．
• 可能な最大仕事数を求めよ．

分割マトロイドとスケジューリングマトロイドの重みなし交差なので，多項式時間で解けます．最適化の余地はいろいろありますが，以下では凝った最適化をせずに実装します（多少の最適化は，最短路の O(n^2 log n) に吸われる）．

左に分割マトロイドをおき，右にスケジューリングマトロイドを置くことにしますが，実装の簡単のため，前回述べた，マトロイドのフローへの一般化をしておきます．すなわち，左頂点集合を発注者の集合（容量２），右頂点集合をジョブの集合にします．

左-右の枝は，各ジョブと対応する発注者の間に枝をおきます．左-左の枝はありません．

右-右の枝がポイントです．「x を入れるために抜くべき y」の列挙が必要になりますが，前回でスケジューリングマトロイドに親しんだ方なら「自分よりも締め切りが早い y」すべてが候補になることがわかるはずです（前回このマトロイドを取り上げたのは，交換枝が非自明だけど列挙が簡単な例だからでした）．

以下が実装例ですが，いつも以上に実装に自信がないので，現役勢は頑張ってください．

struct JobSchInt {
  vector<int> job;
  vector<vector<int>> JobToOwner, OwnerToJob;
  JobSchInt(int owner) : OwnerToJob(owner) { }

  void addJob(int time, int owner) {
    int j = job.size();
    job.push_back(time);
    JobToOwner.push_back( vector<int>(1, j) );
    OwnerToJob[owner].push_back(j);
  }
  int J, O;
  vector<int> capacity;
  vector<int> assign;

  enum { LEFT, RIGHT };
  struct Node { int side, index; };
  vector<int> X;
  vector<vector<Node>> prev;

  bool build() {
    J = JobToOwner.size();
    O = OwnerToJob.size();
    capacity.assign(O, 2); // at most 2 jobs per owner
    assign.assign(J, -1);
  }

  bool append(int u) {
    int i = X.size();
    for (; i > 0; --i) {
      if (X[i-1] <= job[u]) break;
      if (X[i-1] <= i) return false;
    }
    if (i >= job[u]) return false;

    int side = RIGHT;
    while (1) {
      int pside = prev[side][u].side;
      int pu = prev[side][u].index;
      if (pu == -1) {
        --capacity[u];
        break;
      }
      if (side == RIGHT && pside == LEFT) assign[u] = pu;
      if (side == LEFT && pside == RIGHT) assign[u] = -1;
      u = pu;
      side = pside;
    }
    return true;
  }

  bool augment() {
    X.clear();
    for (int j = 0; j < J; ++j) 
      if (assign[j] >= 0) X.push_back(job[j]);
    sort(X.begin(), X.end());

    prev.clear();
    prev.push_back(vector<Node>(O, {-1,-1}));
    prev.push_back(vector<Node>(J, {-1,-1}));

    queue<Node> Q;
    for (int o = 0; o < O; ++o) {
      if (capacity[o] == 0) continue;
      Q.push({LEFT, o});
      prev[LEFT][o] = {0, -1};
    }

    while (!Q.empty()) {
      Node n = Q.front(); Q.pop();
      int u = n.index;
      
      if (n.side == LEFT) {
        for (int v : OwnerToJob[u]) {
          if (prev[RIGHT][v].side != -1) continue;
          if (assign[v] != u) {
            prev[RIGHT][v] = {LEFT, u};
            Q.push({RIGHT, v});
          }
        }
      } else {
        if (assign[u] != -1) {
          int v = assign[u];
          if (prev[LEFT][v].side == -1) {
            prev[LEFT][v] = {RIGHT, u};
            Q.push({LEFT, v});
          }
        } else {
          if (append(u)) return true;
          for (int v = 0; v < J; ++v) {
            if (prev[RIGHT][v].side != -1) continue;
            if (assign[v] != -1 && job[v] <= job[u]) {
              prev[RIGHT][v] = {RIGHT, u};
              Q.push({RIGHT, v});
            }
          }
        }
      }
    }
    return false;
  }
  int maximumProfit() {
    build();
    int size = 0;
    while (augment()) {
      ++size;
    }
    return size;
  }
};

void test()  {
  JobSchInt JSI(2);
  JSI.addJob(1, 0);
  JSI.addJob(1, 1);
  JSI.addJob(2, 0);
  JSI.addJob(3, 0);

  int a = JSI.maximumProfit();
  int b = 3;
  if (a != b) {
    cout << a << " != " << b << endl;
    exit(-1);
  }
}

参考文献

• M. Iri and N. Tomizawa (1976): "An algorithm for finding an optimal independent assignment", Journal of the Operations Research Society of Japan, vol.19, no.1, pp.32-57.
• W. H. Cunningham (1986): "Improved bounds for matroid partition and intersection algorithms", SIAM Journal on Computing, vol.15, no.4, pp.948-957.
• N. J. Harvey (2009): "Algebraic structures and algorithms for matching and matroid Problems", SIAM Journal on Computing, vol.39, no.2, pp.679-702.

Matroid maximization

12:00

今回から数回連続でマトロイドの基本的な話題を扱います．（想定１ヶ月程度）

割と教科書的な内容となりますが，「『マトロイド』が競技プログラミング界隈でよく使われている」といった印象があまり無いので（※事実誤認でしたらすみません．全域木としては頻出なのですが……），今回は布教活動のつもりの紹介です．一応定義から書きますが，きちんとした解説を書くのはキリがないし，目的でもないので，「昔どこかで勉強したけどよく覚えてない」程度の知識は仮定します．

なお，余談として述べていることは，普通のマトロイドの節には載っていない内容ですが，これはこれで別の布教活動のつもりなので，ふーんと思ってくれると嬉しいです．興味のある方は [室田 2007] が読みやすい入門書です．

定義

マトロイドというのは行列（matrix）のようなもの（-oid）から来ている言葉で，行列の線形独立性の抽象化として [Whitney 1935] に提案されました．集合 V とその部分集合の族 I の組 (V, I) がマトロイドであるとは，以下の３条件を満たすことをいいます．

• 空集合は I に入る
• X ∈ I のとき，任意の X の部分集合 Z もまた Z ∈ I
• X, Y ∈ I, |X| < |Y| のとき，ある y ∈ Y-X が存在して X+{x} ∈ X

上２つの条件をみたすものを独立性システムとよび，３つ目の条件；増加公理（augmentation property）が加わったものをマトロイドといいます．X ∈ I なる X のことを独立集合といいます．

極大独立集合を基といいます．上では独立集合から基を定義しましたが，逆に基から独立集合を定義することができます（基の部分集合を独立集合にとる）．基全体を B であらわします．このとき基の満たすべき公理は以下のものです．

• X, Y ∈ B のとき，任意の x ∈ X-Y に対してある y ∈ Y-X が存在して X-{x}+{y}, Y-{y}+{x} ∈ B

これは基交換公理（base exchange property）と呼ばれています．マトロイドの基はすべて同サイズとなることが証明できますが，これは簡単な非マトロイド性の基準として使えます（i.e., 基のサイズがバラバラならマトロイドでない，すなわち貪欲解法の望みは薄い）．

※余談：基交換公理は「離散凸性」を表現しているものと見ることができます．すなわち，X, Y を {0,1}^n 上の点と同一視すると，X-{x}+{y} は X から Y に近づく点で，Y-{y}+{x} は Y から X に近づく点なので，凸集合の性質「2点を結ぶ線分の内側に進むとまた集合に入る」に対応していると考えられます．これを Z^n に拡張したものは M凸 集合とよばれます（Mはmatroidの頭文字）．

覚えておくべきマトロイド

後述するように，マトロイド上の最適化問題は貪欲法で解けます．したがって，多くのマトロイドを知っておけば，問題を見た瞬間に「あっこれマトロイドだ！」といって解法を思いつくことができます．「典型DP」を覚えておくのと同じノリで「典型マトロイド」を覚えておくのが重要だと思います．

具体的に抑えておくべきこととしては，

1. 有名なマトロイド：自由マトロイド・線型マトロイド・グラフ的マトロイド・マッチングマトロイド，ガンモイドなど
2. 既知のマトロイドから新しいマトロイドを作る操作： 制限，縮約，マイナー，ユニオン，誘導

の2種類があります．有名なマトロイドは上の名前で検索すれば出るのでここでは繰り返しません．

既知のマトロイドから新しいマトロイドを作る操作ですが，誘導以外は Wikipedia を参照すれば出ています．誘導（induction）はあまり説明されているものを見ないのですが，とても便利です．これは

• 既知のマトロイドを左におき，新しい台集合を右におき，任意に二部グラフを描く
• 二部グラフマッチングの左端点集合が左で独立のとき，右端点集合を右で独立と定義する

というものです（いわゆる横断マトロイドの一般化バージョン）．これを知っているだけで，多くのマトロイドを構成することができます．特に自由マトロイドからの誘導が重要で，以下のマトロイドはそのようなレシピで構成できます．

• 一様マトロイド：すべてのサイズ k 以下の部分集合を独立とするマトロイド．左にサイズ k の自由マトロイドをおいて，完全二部グラフにして誘導．
• 分割マトロイド：各要素に色（赤，青，…）がついているとして，|赤い要素|≦c1，|青い要素|≦c2，... なるものを独立とするマトロイド．左側に c1+c2+... 個の点をおいて，c1から赤要素に完全，c2から青要素に完全，... というふうに枝を引いて自由マトロイドを誘導．
• 横断マトロイド：有向グラフの頂点を左，辺を右におき，点 u から辺 (v, u) (for all v) に枝を引いて自由マトロイドを誘導すると，各頂点への入次数が1であるような辺集合を独立としてマトロイドになる．マトロイド交差の片側としてよく出てくる例．

※余談：本当は，マッチングよりもフローに一般化したほうが綺麗です．例えば一様マトロイドは「容量 k の 1点」を左においてフローを流すと考える．分割マトロイドは「容量 c1, ..., cm の m 点」からフローを流すと考える．この構成のためにはマトロイド（{0,1}^n上）をM凸集合（Z^n上）に一般化しておく必要があります．

マトロイドと貪欲法

マトロイドで（最も）重要なことは「マトロイド上の最適化問題は貪欲法で解ける」という点です．V の各要素 x に重み w(x) が対応しているとし，w(X) = Σ[x∈X] w(x) と書きます（重みがあるマトロイドを重み付きマトロイドと呼びます）．このとき

  maximize w(X) subject to X ∈ B

という問題は，重みの値によらず>，以下の貪欲法で解けます．

for x ∈ X : w(x) decreasing order
  if X + {x} ∈ I
    X = X + {x}

逆に「任意の重みで貪欲法が最適解を出す」ならば，マトロイドであることが示せます．このことを用いると，「重みが入力として与えられる問題」で「貪欲で解けそう」ならばマトロイド性を疑ってみる手があります．

最適化問題のパターンとしては，独立集合の中の最大を求める問題

  maximize w(X) subject to X ∈ I

というものもあり，先の問題と解は不一致です．こちらの場合は，貪欲法の途中で得られる解の max を覚えておけば大丈夫です．

※余談：V の各部分集合 X に対して独立に w(X) の値が決まっているものを付値（valuation）といいます．マトロイドに「凸っぽい」付値が備わったものを付値マトロイドとよび，その場合も貪欲法で最適解が出ます．これは凸関数が降下法で最小値が求まることと対応します．

具体例：単位時間ジョブスケジューリング

マトロイドの典型例ですが，次回に備えてやっておきます．[Cormen-Leiserson-Rivest-Stein 2009] にも載っていたはずです．

仕事 V = { v1, ..., vn } が与えられる．各仕事は単位時間で終わらせることができる．取り掛かった仕事を途中で中断することはできない．各仕事には締め切り t(v1), ..., t(vn) と報酬 p(v1), ..., p(vn) が決められていて，締め切りより前に終わらせた場合はその仕事の報酬を得ることができる．

最も多くの報酬を得る方法を求めよ（報酬のみ）．

なんとなく貪欲で解ける気がするので，マトロイドを疑います．「実行可能な仕事の集合」の性質を考えると，以下が同値であることがわかります：

• 集合 X の仕事をすべて締め切り内に終わらせるスケジューリングが存在する
• 任意の x ∈ X について，t(x) までに終わらせるべき仕事の数は t(x) 個以下．すなわち |{ y | t(y) ≦ t(x) }| ≦ t(x)

上から下は自明で，下から上は X を時刻でソートしたものがスケジューリングになります．このことを踏まえると増加公理の成立がチェックできるため※，締め切り内に終わらせられることを独立とするとマトロイドになることが従います．（※：サイズの大きい仕事集合の一番最後の仕事を挿入する．同じ仕事がある場合は後詰めして次へ．）

さて，マトロイドとわかったので，貪欲で解けます．問題は「X ∈ I のとき X+{x} ∈ I かどうかを判定する」ルーチン（増加オラクル）の実装部分．自分の体感では「マトロイド上の最適化」に落ちる問題の実装でハードなのは，だいたい具体的なオラクルの実装です．問題（マトロイド）依存の部分なので解説書もなく，とにかく数をこなすしか無いのが現状だと思います．

この問題の場合，上の説明にあるように X を時刻でソートしたものが実行可能スケジュールになるという点に注目します．X に対して後ろから挿入ソートしつつ要素数制約をチェックすれば，1回の判定が O(n) で実行できてトータルで O(n^2) になります．

struct JobSchedulingN {
  struct Job { int t, p; };
  vector<Job> V;
  bool insert(vector<int>& X, int t) {
    int i = X.size();
    for (; i > 0; --i) {
      if (X[i-1] <= t) break;
      if (X[i-1] <= i) return false;
    }
    if (i >= t) return false;
    X.insert(X.begin()+i, t);
    return true;
  }
  void addJob(int t, int p) { V.push_back({t, p}); }
  int maximumProfit() {
    sort(V.begin(), V.end(), [](Job a, Job b) { return a.p > b.p; });
    vector<int> X;
    int profit = 0;
    for (Job x : V) 
      if (insert(X, x.t)) profit += x.p;
    return profit;
  }
};

もうひといき

もうすこし頑張って計算量を落とすことができます．ソートしたものを使うのだから平衡二分探索木を使うのが自然です．要素数に関する不等式を効率的にチェックできればいいので，平衡二分探索木の各ノード x に margin(x) := t(x) - |{ y | t(y) ≦ t(x) }| を覚えさせて，これの min を評価してやれば十分です．

平衡二分探索木の選定ですが，要素を挿入する際に「その要素より上のすべてのノードのmarginを１減らす」という範囲操作が発生するため，「その要素より上」の切り出し，すなわちsplitを実装しやすいものを使うと実装が簡単です．具体的にはTreapやSplay Treeをおすすめします．

以下が実装例です．merge-splitベースのTreapで，範囲操作を遅延評価で実装します．計算量 O(n log n)．

struct JobScheduling {
  struct Job { int t, p; };
  vector<Job> V;

  struct Treap {
    struct Node { 
      int key, pr, size, margin, dmargin, mmargin;
      Node *ch[2];
    } *root;
    Treap() : root(0) { }

    int MMARGIN(Node *n) { return n ? n->mmargin+n->dmargin : 99999999; }
    int SIZE(Node *n) { return n ? n->size : 0; }
    Node *update(Node *n) {
      if (!n) return n;
      if (n->ch[0]) n->ch[0]->dmargin += n->dmargin;
      if (n->ch[1]) n->ch[1]->dmargin += n->dmargin;
      n->margin += n->dmargin; n->dmargin = 0;
      n->mmargin = min({n->margin, MMARGIN(n->ch[0]), MMARGIN(n->ch[1])});
      n->size = 1 + SIZE(n->ch[0]) + SIZE(n->ch[1]);
      return n;
    }
    Node *link(Node *n, Node *c, int b) { n->ch[b] = update(c); return update(n); }
    Node *merge(Node *l, Node *r) {
      l = update(l); r = update(r);
      if (!l || !r) return l ? l : r;
      if (l->pr < r->pr) return link(l, merge(l->ch[1], r), 1);
      else               return link(r, merge(l, r->ch[0]), 0);
    }
    pair<Node*, Node*> split(Node *n, int key) {
      if (!update(n)) return make_pair(n, n);
      if (n->key < key) {
        pair<Node*, Node*> p = split(n->ch[1], key);
        return make_pair(link(n, p.first, 1), update(p.second));
      } else {
        pair<Node*, Node*> p = split(n->ch[0], key);
        return make_pair(update(p.first), link(n, p.second, 0));
      }
    }
    bool insert(int t) {
      pair<Node*, Node*> p = split(root, t);
      bool augmentable = MMARGIN(p.second) > 0 && SIZE(p.first) < t;
      if (augmentable) {
        if (p.second) p.second->dmargin -= 1;
        Node *s = update(new Node({t, rand(), 1, t-SIZE(p.first)-1}));
        p.second = merge(s, p.second);
      }
      root = merge(p.first, p.second);
      return augmentable;
    }
  };
  void addJob(int t, int p) { V.push_back({t, p}); }
  int maximumProfit() {
    sort(V.begin(), V.end(), [](Job a, Job b) { return a.p > b.p; });
    int profit = 0;
    Treap X;
    for (Job x : V) 
      if (X.insert(x.t)) profit += x.p;
    return profit;
  }
};

参考文献

• H. Whitney (1935): "On the abstract properties of linear dependence", American Journal of Mathematics, vol. 57, pp.509-533.
• 室田一雄 (2007): "離散凸解析の考えかた 最適化における離散と連続の数理", 共立出版.
• T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein (2009): "Introduction to Algorithms", The MIT Press; 3rd ed.