マトロイド交叉 (src/Optimization/matroid_intersection.hpp)

• View this file on GitHub
• View document part on GitHub
• Last update: 2023-08-10 15:00:06+09:00
• Include: #include "src/Optimization/matroid_intersection.hpp"

概要

2つのマトロイド $M_1 = (S, \mathcal{I}_1), M_2 = (S, \mathcal{I}_2)$ が与えられたとき、共通独立集合 $I \in \mathcal{I}_1 \cap \mathcal{I}_2$ のうち、要素数が最大のもの（最大共通独立集合）や、要素に重みがある場合に重みが最大／最小となるものを求める問題。

このライブラリは、重みなし・重み付きのマトロイド交叉問題を解くための関数を提供します。

---

マトロイドインターフェース

マトロイドを表すクラスは、以下のインターフェースを満たす必要があります。

struct YourMatroid {
  // 現在の独立集合 I に基づいて内部状態を構築する
  void build(const std::vector<int>& I);

  // 要素 e が I に追加可能かを判定する
  // I は build で渡された集合
  bool oracle(int e) const;

  // I に e を追加したとき、I から f を除けば独立性を保てるかを判定する
  // (I \ {f}) U {e} が独立であるかを判定する
  bool oracle(int e, int f) const;
};

---

提供されるマトロイド

GraphicMatroid

グラフ理論における森（閉路を含まない辺集合）を独立集合とするマトロイドです。

• GraphicMatroid(int n): n頂点のグラフマトロイドを構築する。
• int add_edge(int u, int v): 辺(u, v)を追加し、その辺のIDを返す。

PartitionMatroid

台集合 $S$ を disjoint な部分集合 $S_1, S_2, \dots, S_k$ に分割し、各部分集合から高々 $R_i$ 個の要素を選ぶ集合を独立集合とするマトロイドです。

• PartitionMatroid(int m, const std::vector<std::vector<int>>& parts, const std::vector<int>& R):
  • m: 台集合の要素数
  • parts: 台集合の分割。parts[i]は $S_i$ の要素のリスト。
  • R: 各部分集合から選べる要素数の上限。R[i]は $S_i$ から選べる上限数。

---

関数

重みなし

関数名	概要	計算量
matroid_intersection(n,M1,M2)	2つのマトロイドM1, M2の最大共通独立集合を求める。 n: 台集合のサイズ。	$O(nr^{1.5}\tau)$

重みあり

関数名	概要	計算量
weighted_matroid_intersection<sgn>(n,M1,M2,w)	2つのマトロイドM1, M2の重み付き最大／最小共通独立集合を求める。 n: 台集合のサイズ。 w: 各要素の重み。 sgn: MAXIMIZEまたはMINIMIZEを指定。 返り値: std::vector<std::vector<int>>。ret[k]に、サイズkの共通独立集合のうち重みが最大／最小のものが格納される。	$O(nr^2\tau)$

補足:

• $n$は台集合のサイズ, $r$は最大共通独立集合のサイズ, $\tau$は独立性オラクル呼び出しにかかる時間。

---

コード例

色付きの辺を持つグラフで、各色の辺の数を制限しながら全域森を求める問題を考えます。これはグラフマトロイドと分割マトロイドの交叉問題として定式化できます。

#include <iostream>
#include "Optimization/matroid_intersection.hpp"

int main() {
    // 頂点数6、辺数8のグラフ
    int n_vertices = 6;
    int n_edges = 8;

    GraphicMatroid m1(n_vertices);
    // edges[i] = {u, v, color}
    std::vector<std::tuple<int, int, int>> edges = {
        {0, 1, 0}, {0, 2, 1}, {1, 2, 0}, {1, 3, 2},
        {2, 4, 1}, {3, 4, 0}, {3, 5, 2}, {4, 5, 1}
    };
    for (const auto& [u, v, c] : edges) {
        m1.add_edge(u, v);
    }

    // 分割マトロイド: 各色の辺は1本まで
    // part 0: color 0, part 1: color 1, part 2: color 2
    std::vector<std::vector<int>> parts(3);
    for (int i = 0; i < n_edges; ++i) {
        parts[std::get<2>(edges[i])].push_back(i);
    }
    // 各色の上限は1
    std::vector<int> R = {1, 1, 1};
    PartitionMatroid m2(n_edges, parts, R);

    // 最大共通独立集合（= 各色1本以下の辺からなる森）を求める
    std::vector<int> result = matroid_intersection(n_edges, m1, m2);

    // 結果: {0, 4, 6} (辺0, 4, 6)
    // 辺0: (0,1) color 0
    // 辺4: (2,4) color 1
    // 辺6: (3,5) color 2
    // これらは森をなし、各色1本ずつという制約を満たす。
    std::cout << "Maximum common independent set (edge IDs):" << std::endl;
    for (int i = 0; i < result.size(); ++i) {
        std::cout << result[i] << (i == result.size() - 1 ? "" : " ");
    }
    std::cout << std::endl;

    return 0;
}

参考

• https://hitonanode.github.io/cplib-cpp/combinatorial_opt/matroid_intersection.hpp
• William H. Cunningham. Improved bounds for matroid partition and intersection algorithms. SIAMJournal on Computing (SICOMP), 15(4):948–957, 1986.

Verify

• [beecrowd

  2128 Demonstration of Honesty!](https://www.beecrowd.com.br/judge/en/problems/view/2128) (グラフ+分割)

• Hello 2020 G. Seollal - Codeforces (グラフ+分割)
• 2019 Petrozavodsk Winter Camp, Yandex Cup D. Pick Your Own Nim - Codeforces (二値+分割)
• 幾何コンテスト2013 B - 玉座の間 (重み二部マッチング, 重みが実数)

• ACM-ICPC Japan Alumni Group Summer Camp 2012 Day 2 まるかいて (重み二部マッチング,

  台集合

  <= 1e4, 構築, sp judge)

• AtCoder Library Practice Contest E - MinCostFlow (重み付き, 分割+分割)
• Deltix Round, Summer 2021 H. DIY Tree - Codeforces (重み付き, グラフ+分割)
• NAIPC 2018 G Rainbow Graph(重み付き, あるマトロイド (グラフマトロイドの双対マトロイドと自由マトロイドの直和?) 同士の交差)

Depends on

• 最大最小を指定するための列挙型 (src/Optimization/MinMaxEnum.hpp)

Verified with

• test/aoj/1163.matroid_intersection.test.cpp
• test/aoj/1605.weighted_matroid_intersection.test.cpp
• test/aoj/GRL_2_B.weighted_matroid_intersection.test.cpp
• test/aoj/GRL_7_A.matroid_intersection.test.cpp
• test/atcoder/abc231_h.weighted_matroid_intersection.test.cpp
• test/misc/cf_gym103373_i.test.cpp
• test/yukicoder/421.matroid_intersection.test.cpp
• test/yukicoder/421.weighted_matroid_intersection.test.cpp

Code

#pragma once
#include <vector>
#include <algorithm>
#include <limits>
#include <array>
#include <queue>
#include <cassert>
#include "src/Optimization/MinMaxEnum.hpp"
template <typename Matroid1, typename Matroid2> std::vector<int> matroid_intersection(int n, Matroid1 M1, Matroid2 M2) {
 std::vector<int> b(n, false), pre(n), I[2];
 for (int e= 0; e < n; e++) I[0].push_back(e);
 M1.build(I[1]), M2.build(I[1]);
 for (bool converged= false; !converged;) {
  pre.assign(n, false);
  std::vector L(1, std::vector<int>());
  for (int u: I[0])
   if (M1.oracle(u)) pre[u]= true, L[0].push_back(u);
  int m= 0;
  for (; L.back().size(); m+= 2) {
   L.push_back({});
   for (int e: L[m]) {
    if (converged= M2.oracle(e)) break;
    for (int f: I[1])
     if (!pre[f] && M2.oracle(f, e)) L[m + 1].push_back(f), pre[f]= true;
   }
   if (converged) break;
   L.push_back({});
   for (int e: L[m + 1])
    for (int f: I[0])
     if (!pre[f] && M1.oracle(e, f)) L[m + 2].push_back(f), pre[f]= true;
  }
  if (!converged) break;
  std::vector<std::vector<int>> L2(m + 1);
  for (int e: L[m])
   if (M2.oracle(e)) L2[m].push_back(e);
  for (int i= m; i; i-= 2) {
   for (int e: L[i - 1])
    for (int f: L2[i])
     if (M1.oracle(e, f)) {
      L2[i - 1].push_back(e);
      break;
     }
   for (int e: L[i - 2])
    for (int f: L2[i - 1])
     if (M2.oracle(f, e)) {
      L2[i - 2].push_back(e);
      break;
     }
  }
  pre.assign(n, -1);
  for (int e: L2[0])
   if (M1.oracle(e)) {
    bool isok= false;
    if (m) {
     std::vector<size_t> ei(m);
     for (int i= 0; e != -1;) {
      if (ei[i] == L2[i + 1].size()) e= pre[e], ei[i--]= 0;
      else if (int f= L2[i + 1][ei[i]++]; pre[f] == -1 && (i & 1 ? M1.oracle(e, f) : M2.oracle(f, e)))
       if (pre[f]= e, e= f; ++i == m) {
        if (M2.oracle(e))
         for (isok= true; e != -1; e= pre[e]) b[e]= !b[e];
        else e= pre[e], --i;
       }
     }
    } else if (M2.oracle(e)) isok= true, b[e]= 1;
    if (isok) {
     converged= false, I[0].clear(), I[1].clear();
     for (int u= 0; u < n; u++) I[b[u]].push_back(u);
     M1.build(I[1]), M2.build(I[1]);
    }
   }
 }
 return I[1];
}
template <MinMaxEnum sgn, class Matroid1, class Matroid2, class cost_t> std::vector<std::vector<int>> weighted_matroid_intersection(int n, Matroid1 M1, Matroid2 M2, std::vector<cost_t> c) {
 assert(n == (int)c.size());
 bool b[n];
 std::fill_n(b, n, false);
 std::vector<int> I[2], p;
 std::vector<std::vector<int>> ret(1);
 for (int u= 0; u < n; u++) I[0].push_back(u);
 if constexpr (sgn == MAXIMIZE) {
  auto cmx= *std::max_element(c.begin(), c.end());
  for (auto &x: c) x-= cmx;
 } else {
  auto cmi= *std::min_element(c.begin(), c.end());
  for (auto &x: c) x-= cmi;
 }
 for (auto &x: c) x*= -sgn * (n + 1);
 for (bool converged= false; !converged;) {
  converged= true, M1.build(I[1]), M2.build(I[1]);
  std::priority_queue<std::pair<cost_t, int>> pq;
  std::vector<cost_t> dist(n, std::numeric_limits<cost_t>::lowest());
  for (int u: I[0])
   if (M1.oracle(u)) pq.emplace(dist[u]= c[u] - 1, u);
  for (p.assign(n, -1); pq.size();) {
   auto [d, u]= pq.top();
   if (pq.pop(); d != dist[u]) continue;
   if (b[u]) {
    for (int v: I[0])
     if (M1.oracle(u, v))
      if (cost_t cost= d + c[v] - 1; dist[v] < cost) pq.emplace(dist[v]= cost, v), p[v]= u;
   } else {
    if (M2.oracle(u)) {
     for (int v= u; v != -1; v= p[v]) b[v]= !b[v];
     I[0].clear(), I[1].clear(), converged= false;
     for (int u= 0; u < n; u++) I[b[u]].push_back(u);
     ret.emplace_back(I[1]);
     break;
    }
    for (int v: I[1])
     if (M2.oracle(v, u))
      if (cost_t cost= d - c[v] - 1; dist[v] < cost) pq.emplace(dist[v]= cost, v), p[v]= u;
   }
  }
 }
 return ret;
}
class GraphicMatroid {
 int n;
 std::vector<std::array<int, 2>> es;
 std::vector<int> g, pos, comp, in, out;
 inline bool is_ancestor(int u, int v) const { return in[u] <= in[v] && in[v] < out[u]; }
public:
 GraphicMatroid(int n_): n(n_), comp(n), in(n), out(n) {}
 int add_edge(int u, int v) { return es.push_back({u, v}), es.size() - 1; }
 void build(const std::vector<int> &I) {
  in.assign(n, -1), g.resize(I.size() * 2), pos.assign(n + 1, 0);
  for (int e: I) {
   auto [u, v]= es[e];
   ++pos[u], ++pos[v];
  }
  for (int i= 0; i < n; ++i) pos[i + 1]+= pos[i];
  for (int e: I) {
   auto [u, v]= es[e];
   g[--pos[u]]= v, g[--pos[v]]= u;
  }
  std::vector<int> ei(pos.begin(), pos.begin() + n), pre(n, -1);
  for (int u= 0, t= 0, p; u < n; u++)
   if (in[u] == -1)
    for (in [comp[u]= p= u]= t++; p >= 0;) {
     if (ei[p] == pos[p + 1]) out[p]= t, p= pre[p];
     else if (int v= g[ei[p]++]; in[v] == -1) comp[v]= comp[u], pre[v]= p, in[p= v]= t++;
    }
 }
 inline bool oracle(int e) const { return comp[es[e][0]] != comp[es[e][1]]; }
 inline bool oracle(int e, int f) const {
  if (oracle(f)) return true;
  return e= es[e][in[es[e][0]] < in[es[e][1]]], is_ancestor(e, es[f][0]) != is_ancestor(e, es[f][1]);
 }
};
struct PartitionMatroid {
 std::vector<int> belong, R, cnt;
 PartitionMatroid(int m_, const std::vector<std::vector<int>> &parts, const std::vector<int> &R_): belong(m_, -1), R(R_) {
  assert(parts.size() == R.size());
  for (int i= parts.size(); i--;)
   for (int e: parts[i]) belong[e]= i;
 }
 void build(const std::vector<int> &I) {
  cnt= R;
  for (int e: I)
   if (belong[e] != -1) cnt[belong[e]]--;
 }
 inline bool oracle(int e) const { return belong[e] == -1 || cnt[belong[e]] > 0; }
 inline bool oracle(int e, int f) const { return oracle(f) || belong[e] == belong[f]; }
};

#line 2 "src/Optimization/matroid_intersection.hpp"
#include <vector>
#include <algorithm>
#include <limits>
#include <array>
#include <queue>
#include <cassert>
#line 2 "src/Optimization/MinMaxEnum.hpp"
enum MinMaxEnum { MAXIMIZE= -1, MINIMIZE= 1 };
#line 9 "src/Optimization/matroid_intersection.hpp"
template <typename Matroid1, typename Matroid2> std::vector<int> matroid_intersection(int n, Matroid1 M1, Matroid2 M2) {
 std::vector<int> b(n, false), pre(n), I[2];
 for (int e= 0; e < n; e++) I[0].push_back(e);
 M1.build(I[1]), M2.build(I[1]);
 for (bool converged= false; !converged;) {
  pre.assign(n, false);
  std::vector L(1, std::vector<int>());
  for (int u: I[0])
   if (M1.oracle(u)) pre[u]= true, L[0].push_back(u);
  int m= 0;
  for (; L.back().size(); m+= 2) {
   L.push_back({});
   for (int e: L[m]) {
    if (converged= M2.oracle(e)) break;
    for (int f: I[1])
     if (!pre[f] && M2.oracle(f, e)) L[m + 1].push_back(f), pre[f]= true;
   }
   if (converged) break;
   L.push_back({});
   for (int e: L[m + 1])
    for (int f: I[0])
     if (!pre[f] && M1.oracle(e, f)) L[m + 2].push_back(f), pre[f]= true;
  }
  if (!converged) break;
  std::vector<std::vector<int>> L2(m + 1);
  for (int e: L[m])
   if (M2.oracle(e)) L2[m].push_back(e);
  for (int i= m; i; i-= 2) {
   for (int e: L[i - 1])
    for (int f: L2[i])
     if (M1.oracle(e, f)) {
      L2[i - 1].push_back(e);
      break;
     }
   for (int e: L[i - 2])
    for (int f: L2[i - 1])
     if (M2.oracle(f, e)) {
      L2[i - 2].push_back(e);
      break;
     }
  }
  pre.assign(n, -1);
  for (int e: L2[0])
   if (M1.oracle(e)) {
    bool isok= false;
    if (m) {
     std::vector<size_t> ei(m);
     for (int i= 0; e != -1;) {
      if (ei[i] == L2[i + 1].size()) e= pre[e], ei[i--]= 0;
      else if (int f= L2[i + 1][ei[i]++]; pre[f] == -1 && (i & 1 ? M1.oracle(e, f) : M2.oracle(f, e)))
       if (pre[f]= e, e= f; ++i == m) {
        if (M2.oracle(e))
         for (isok= true; e != -1; e= pre[e]) b[e]= !b[e];
        else e= pre[e], --i;
       }
     }
    } else if (M2.oracle(e)) isok= true, b[e]= 1;
    if (isok) {
     converged= false, I[0].clear(), I[1].clear();
     for (int u= 0; u < n; u++) I[b[u]].push_back(u);
     M1.build(I[1]), M2.build(I[1]);
    }
   }
 }
 return I[1];
}
template <MinMaxEnum sgn, class Matroid1, class Matroid2, class cost_t> std::vector<std::vector<int>> weighted_matroid_intersection(int n, Matroid1 M1, Matroid2 M2, std::vector<cost_t> c) {
 assert(n == (int)c.size());
 bool b[n];
 std::fill_n(b, n, false);
 std::vector<int> I[2], p;
 std::vector<std::vector<int>> ret(1);
 for (int u= 0; u < n; u++) I[0].push_back(u);
 if constexpr (sgn == MAXIMIZE) {
  auto cmx= *std::max_element(c.begin(), c.end());
  for (auto &x: c) x-= cmx;
 } else {
  auto cmi= *std::min_element(c.begin(), c.end());
  for (auto &x: c) x-= cmi;
 }
 for (auto &x: c) x*= -sgn * (n + 1);
 for (bool converged= false; !converged;) {
  converged= true, M1.build(I[1]), M2.build(I[1]);
  std::priority_queue<std::pair<cost_t, int>> pq;
  std::vector<cost_t> dist(n, std::numeric_limits<cost_t>::lowest());
  for (int u: I[0])
   if (M1.oracle(u)) pq.emplace(dist[u]= c[u] - 1, u);
  for (p.assign(n, -1); pq.size();) {
   auto [d, u]= pq.top();
   if (pq.pop(); d != dist[u]) continue;
   if (b[u]) {
    for (int v: I[0])
     if (M1.oracle(u, v))
      if (cost_t cost= d + c[v] - 1; dist[v] < cost) pq.emplace(dist[v]= cost, v), p[v]= u;
   } else {
    if (M2.oracle(u)) {
     for (int v= u; v != -1; v= p[v]) b[v]= !b[v];
     I[0].clear(), I[1].clear(), converged= false;
     for (int u= 0; u < n; u++) I[b[u]].push_back(u);
     ret.emplace_back(I[1]);
     break;
    }
    for (int v: I[1])
     if (M2.oracle(v, u))
      if (cost_t cost= d - c[v] - 1; dist[v] < cost) pq.emplace(dist[v]= cost, v), p[v]= u;
   }
  }
 }
 return ret;
}
class GraphicMatroid {
 int n;
 std::vector<std::array<int, 2>> es;
 std::vector<int> g, pos, comp, in, out;
 inline bool is_ancestor(int u, int v) const { return in[u] <= in[v] && in[v] < out[u]; }
public:
 GraphicMatroid(int n_): n(n_), comp(n), in(n), out(n) {}
 int add_edge(int u, int v) { return es.push_back({u, v}), es.size() - 1; }
 void build(const std::vector<int> &I) {
  in.assign(n, -1), g.resize(I.size() * 2), pos.assign(n + 1, 0);
  for (int e: I) {
   auto [u, v]= es[e];
   ++pos[u], ++pos[v];
  }
  for (int i= 0; i < n; ++i) pos[i + 1]+= pos[i];
  for (int e: I) {
   auto [u, v]= es[e];
   g[--pos[u]]= v, g[--pos[v]]= u;
  }
  std::vector<int> ei(pos.begin(), pos.begin() + n), pre(n, -1);
  for (int u= 0, t= 0, p; u < n; u++)
   if (in[u] == -1)
    for (in [comp[u]= p= u]= t++; p >= 0;) {
     if (ei[p] == pos[p + 1]) out[p]= t, p= pre[p];
     else if (int v= g[ei[p]++]; in[v] == -1) comp[v]= comp[u], pre[v]= p, in[p= v]= t++;
    }
 }
 inline bool oracle(int e) const { return comp[es[e][0]] != comp[es[e][1]]; }
 inline bool oracle(int e, int f) const {
  if (oracle(f)) return true;
  return e= es[e][in[es[e][0]] < in[es[e][1]]], is_ancestor(e, es[f][0]) != is_ancestor(e, es[f][1]);
 }
};
struct PartitionMatroid {
 std::vector<int> belong, R, cnt;
 PartitionMatroid(int m_, const std::vector<std::vector<int>> &parts, const std::vector<int> &R_): belong(m_, -1), R(R_) {
  assert(parts.size() == R.size());
  for (int i= parts.size(); i--;)
   for (int e: parts[i]) belong[e]= i;
 }
 void build(const std::vector<int> &I) {
  cnt= R;
  for (int e: I)
   if (belong[e] != -1) cnt[belong[e]]--;
 }
 inline bool oracle(int e) const { return belong[e] == -1 || cnt[belong[e]] > 0; }
 inline bool oracle(int e, int f) const { return oracle(f) || belong[e] == belong[f]; }
};

Back to top page