- できるだけ毎週土曜日更新．次回更新は2月16日．
- ネタ切れ深刻なので記事リクエスト募集中．twitterないしはコメント欄まで．
- 貼っているコードはあんまり検証していないので，信頼性はアレです．アレ．
- twitter: @tmaehara フォローはお気軽に．
2012-11-03
Matroid maximization
今回から数回連続でマトロイドの基本的な話題を扱います．（想定１ヶ月程度）
割と教科書的な内容となりますが，「『マトロイド』が競技プログラミング界隈でよく使われている」といった印象があまり無いので（※事実誤認でしたらすみません．全域木としては頻出なのですが……），今回は布教活動のつもりの紹介です．一応定義から書きますが，きちんとした解説を書くのはキリがないし，目的でもないので，「昔どこかで勉強したけどよく覚えてない」程度の知識は仮定します．
なお，余談として述べていることは，普通のマトロイドの節には載っていない内容ですが，これはこれで別の布教活動のつもりなので，ふーんと思ってくれると嬉しいです．興味のある方は [室田 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」を覚えておくのと同じノリで「典型マトロイド」を覚えておくのが重要だと思います．
具体的に抑えておくべきこととしては，
- 有名なマトロイド：自由マトロイド・線型マトロイド・グラフ的マトロイド・マッチングマトロイド，ガンモイドなど
- 既知のマトロイドから新しいマトロイドを作る操作： 制限，縮約，マイナー，ユニオン，誘導
の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; } };
参考文献