週刊 spaghetti_source

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.