跳转至

一般圖最大匹配

帶花樹算法（Blossom Algorithm）

開花算法（Blossom Algorithm，也被稱做帶花樹）可以解決一般圖最大匹配問題（maximum cardinality matchings）。此算法由 Jack Edmonds 在 1961 年提出。 經過一些修改後也可以解決一般圖最大權匹配問題。 此算法是第一個給出證明説最大匹配有多項式複雜度。

一般圖匹配和二分圖匹配（bipartite matching）不同的是，圖可能存在奇環。

general-matching-1

以此圖為例，若直接取反（匹配邊和未匹配邊對調），會使得取反後的 \(M\) 不合法，某些點會出現在兩條匹配上，而問題就出在奇環。

下面考慮一般圖的增廣算法。 從二分圖的角度出發，每次枚舉一個未匹配點，設出發點為根，標記為 「o」，接下來交錯標記 「o」 和 「i」，不難發現 「i」 到 「o」 這段邊是匹配邊。

假設當前點是 \(v\)，相鄰點為 \(u\)，可以分為以下兩種情況：

  1. \(u\) 未拜訪過，當 \(u\) 是未匹配點，則找到增廣路徑，否則從 \(u\) 的配偶找增廣路。
  2. \(u\) 已拜訪過，遇到標記「o」代表需要 縮花，否則代表遇到偶環，跳過。

遇到偶環的情況，將他視為二分圖解決，故可忽略。縮花 後，再新圖中繼續找增廣路。

general-matching-2

設原圖為 \(G\)，縮花 後的圖為 \(G'\)，我們只需要證明：

  1. 若 \(G\) 存在增廣路，\(G'\) 也存在。
  2. 若 \(G'\) 存在增廣路，\(G\) 也存在。

general-matching-3

設非樹邊（形成環的那條邊）為 \((u,v)\)，定義花根 \(h=LCA(u,v)\)。 奇環是交替的，有且僅有 \(h\) 的兩條鄰邊類型相同，都是非匹配邊。 那麼進入 \(h\) 的樹邊肯定是匹配邊，環上除了 \(h\) 以外其他點往環外的邊都是非匹配邊。

觀察可知，從環外的邊出去有兩種情況，順時針或逆時針。

general-matching-4

於是 縮花 與 不縮花 都不影響正確性。

實作上找到 花 以後我們不需要真的 縮花，可以用數組紀錄每個點在以哪個點為根的那朵花中。

複雜度分析 Complexity Analysis

每次找增廣路，遍歷所有邊，遇到 花 會維護 花 上的點，\(O(|E|^2)\)。

枚舉所有未匹配點做增廣路，總共 \(O(|V||E|^2)\)。

參考代碼

參考代碼 
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
// graph
template <typename T>
class graph {
 public:
  struct edge {
    int from;
    int to;
    T cost;
  };

  vector<edge> edges;
  vector<vector<int> > g;
  int n;

  graph(int _n) : n(_n) { g.resize(n); }

  virtual int add(int from, int to, T cost) = 0;
};

// undirectedgraph
template <typename T>
class undirectedgraph : public graph<T> {
 public:
  using graph<T>::edges;
  using graph<T>::g;
  using graph<T>::n;

  undirectedgraph(int _n) : graph<T>(_n) {}

  int add(int from, int to, T cost = 1) {
    assert(0 <= from && from < n && 0 <= to && to < n);
    int id = (int)edges.size();
    g[from].push_back(id);
    g[to].push_back(id);
    edges.push_back({from, to, cost});
    return id;
  }
};

// blossom / find_max_unweighted_matching
template <typename T>
vector<int> find_max_unweighted_matching(const undirectedgraph<T> &g) {
  std::mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
  vector<int> match(g.n, -1);   // 匹配
  vector<int> aux(g.n, -1);     // 時間戳記
  vector<int> label(g.n);       // 「o」或「i」
  vector<int> orig(g.n);        // 花根
  vector<int> parent(g.n, -1);  // 父節點
  queue<int> q;
  int aux_time = -1;

  auto lca = [&](int v, int u) {
    aux_time++;
    while (true) {
      if (v != -1) {
        if (aux[v] == aux_time) {  // 找到拜訪過的點 也就是LCA
          return v;
        }
        aux[v] = aux_time;
        if (match[v] == -1) {
          v = -1;
        } else {
          v = orig[parent[match[v]]];  // 以匹配點的父節點繼續尋找
        }
      }
      swap(v, u);
    }
  };  // lca

  auto blossom = [&](int v, int u, int a) {
    while (orig[v] != a) {
      parent[v] = u;
      u = match[v];
      if (label[u] == 1) {  // 初始點設為「o」找增廣路
        label[u] = 0;
        q.push(u);
      }
      orig[v] = orig[u] = a;  // 縮花
      v = parent[u];
    }
  };  // blossom

  auto augment = [&](int v) {
    while (v != -1) {
      int pv = parent[v];
      int next_v = match[pv];
      match[v] = pv;
      match[pv] = v;
      v = next_v;
    }
  };  // augment

  auto bfs = [&](int root) {
    fill(label.begin(), label.end(), -1);
    iota(orig.begin(), orig.end(), 0);
    while (!q.empty()) {
      q.pop();
    }
    q.push(root);
    // 初始點設為「o」，這裏以「0」代替「o」，「1」代替「i」
    label[root] = 0;
    while (!q.empty()) {
      int v = q.front();
      q.pop();
      for (int id : g.g[v]) {
        auto &e = g.edges[id];
        int u = e.from ^ e.to ^ v;
        if (label[u] == -1) {  // 找到未拜訪點
          label[u] = 1;        // 標記「i」
          parent[u] = v;
          if (match[u] == -1) {  // 找到未匹配點
            augment(u);          // 尋找增廣路徑
            return true;
          }
          // 找到已匹配點 將與她匹配的點丟入queue 延伸交錯樹
          label[match[u]] = 0;
          q.push(match[u]);
          continue;
        } else if (label[u] == 0 && orig[v] != orig[u]) {
          // 找到已拜訪點 且標記同為「o」代表找到「花」
          int a = lca(orig[v], orig[u]);
          // 找LCA 然後縮花
          blossom(u, v, a);
          blossom(v, u, a);
        }
      }
    }
    return false;
  };  // bfs

  auto greedy = [&]() {
    vector<int> order(g.n);
    // 隨機打亂 order
    iota(order.begin(), order.end(), 0);
    shuffle(order.begin(), order.end(), rng);

    // 將可以匹配的點匹配
    for (int i : order) {
      if (match[i] == -1) {
        for (auto id : g.g[i]) {
          auto &e = g.edges[id];
          int to = e.from ^ e.to ^ i;
          if (match[to] == -1) {
            match[i] = to;
            match[to] = i;
            break;
          }
        }
      }
    }
  };  // greedy

  // 一開始先隨機匹配
  greedy();
  // 對未匹配點找增廣路
  for (int i = 0; i < g.n; i++) {
    if (match[i] == -1) {
      bfs(i);
    }
  }
  return match;
}
UOJ #79. 一般圖最大匹配 
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
#include <bits/stdc++.h>
using namespace std;

// graph
template <typename T>
class graph {
 public:
  struct edge {
    int from;
    int to;
    T cost;
  };

  vector<edge> edges;
  vector<vector<int> > g;
  int n;

  graph(int _n) : n(_n) { g.resize(n); }

  virtual int add(int from, int to, T cost) = 0;
};

// undirectedgraph
template <typename T>
class undirectedgraph : public graph<T> {
 public:
  using graph<T>::edges;
  using graph<T>::g;
  using graph<T>::n;

  undirectedgraph(int _n) : graph<T>(_n) {}

  int add(int from, int to, T cost = 1) {
    assert(0 <= from && from < n && 0 <= to && to < n);
    int id = (int)edges.size();
    g[from].push_back(id);
    g[to].push_back(id);
    edges.push_back({from, to, cost});
    return id;
  }
};

// blossom / find_max_unweighted_matching
template <typename T>
vector<int> find_max_unweighted_matching(const undirectedgraph<T> &g) {
  std::mt19937 rng(114514);  // 这里随机种子是无关紧要的
  // 也可以用 chrono::steady_clock::now().time_since_epoch().count()
  // 获取当前时间
  vector<int> match(g.n, -1);   // 匹配
  vector<int> aux(g.n, -1);     // 时间戳记
  vector<int> label(g.n);       // "o" or "i"
  vector<int> orig(g.n);        // 花根
  vector<int> parent(g.n, -1);  // 父节点
  queue<int> q;
  int aux_time = -1;

  auto lca = [&](int v, int u) {
    aux_time++;
    while (true) {
      if (v != -1) {
        if (aux[v] == aux_time) {  // 找到拜访过的点 也就是LCA
          return v;
        }
        aux[v] = aux_time;
        if (match[v] == -1) {
          v = -1;
        } else {
          v = orig[parent[match[v]]];  // 以匹配点的父节点继续寻找
        }
      }
      swap(v, u);
    }
  };  // lca

  auto blossom = [&](int v, int u, int a) {
    while (orig[v] != a) {
      parent[v] = u;
      u = match[v];
      if (label[u] == 1) {  // 初始点设为"o" 找增广路
        label[u] = 0;
        q.push(u);
      }
      orig[v] = orig[u] = a;  // 缩花
      v = parent[u];
    }
  };  // blossom

  auto augment = [&](int v) {
    while (v != -1) {
      int pv = parent[v];
      int next_v = match[pv];
      match[v] = pv;
      match[pv] = v;
      v = next_v;
    }
  };  // augment

  auto bfs = [&](int root) {
    fill(label.begin(), label.end(), -1);
    iota(orig.begin(), orig.end(), 0);
    while (!q.empty()) {
      q.pop();
    }
    q.push(root);
    // 初始点设为 "o", 这里以"0"代替"o", "1"代替"i"
    label[root] = 0;
    while (!q.empty()) {
      int v = q.front();
      q.pop();
      for (int id : g.g[v]) {
        auto &e = g.edges[id];
        int u = e.from ^ e.to ^ v;
        if (label[u] == -1) {  // 找到未拜访点
          label[u] = 1;        // 标记 "i"
          parent[u] = v;
          if (match[u] == -1) {  // 找到未匹配点
            augment(u);          // 寻找增广路径
            return true;
          }
          // 找到已匹配点 将与她匹配的点丢入queue 延伸交错树
          label[match[u]] = 0;
          q.push(match[u]);
          continue;
        } else if (label[u] == 0 && orig[v] != orig[u]) {
          // 找到已拜访点 且标记同为"o" 代表找到"花"
          int a = lca(orig[v], orig[u]);
          // 找LCA 然后缩花
          blossom(u, v, a);
          blossom(v, u, a);
        }
      }
    }
    return false;
  };  // bfs

  auto greedy = [&]() {
    vector<int> order(g.n);
    // 随机打乱 order
    iota(order.begin(), order.end(), 0);
    shuffle(order.begin(), order.end(), rng);

    // 将可以匹配的点匹配
    for (int i : order) {
      if (match[i] == -1) {
        for (auto id : g.g[i]) {
          auto &e = g.edges[id];
          int to = e.from ^ e.to ^ i;
          if (match[to] == -1) {
            match[i] = to;
            match[to] = i;
            break;
          }
        }
      }
    }
  };  // greedy

  // 一开始先随机匹配
  greedy();
  // 对未匹配点找增广路
  for (int i = 0; i < g.n; i++) {
    if (match[i] == -1) {
      bfs(i);
    }
  }
  return match;
}

int main() {
  ios::sync_with_stdio(false);

  int n, m;
  cin >> n >> m;

  undirectedgraph<int> g(n);

  while (m--) {
    int u, v;
    cin >> u >> v;
    g.add(u - 1, v - 1);  // 0-based
  }

  auto match = find_max_unweighted_matching(g);

  cout << count_if(match.begin(), match.end(), [](int x) { return x != -1; }) /
              2
       << endl;
  for (int i = 0; i < n; i++) cout << match[i] + 1 << " \n"[i == n - 1];

  return 0;
}

基於高斯消元的一般圖匹配算法

提示

在閲讀以下內容前，你可能需要先閲讀「線性代數」部分中關於矩陣的內容：

這一部分將介紹一種基於高斯消元的一般圖匹配算法。與傳統的帶花樹算法相比，它的優勢在於更易於理解與編寫，同時便於解決「最大匹配中的必須點」等問題；缺點在於常數比較大，因為高斯消元的 \(O(n^3)\) 基本是跑滿的，而帶花樹一般跑不滿。

前置知識：Tutte 矩陣

定義：對於一張 \(n\) 個點的無向圖 \(G = (V, E)\)，其 Tutte 矩陣 \(\tilde{A}(G)\) 為一個 \(n \times n\) 的矩陣，其中：

\[ \tilde{A}(G)_{i,j} = \begin{cases} x_{i,j}, & i<j,\; (v_i, v_j)\in E \\ -x_{i,j}, & i > j,\; (v_i, v_j) \in E \\ 0, & \text{otherwise} \end{cases} \]

其中 \(x_{i, j}\) 是一個變量，因此 \(\tilde{A}(G)\) 中共有 \(|E|\) 個變量。

在無歧義的情況下，以下將 \(\tilde{A}(G)\) 簡寫為 \(\tilde{A}\)。

定理（Tutte 定理）：\(G\) 存在完美匹配當且僅當 \(\det \tilde{A} \ne 0\)。

證明

這裏引入「偶環覆蓋」的概念：一個無向圖 \(G\) 的偶環覆蓋指用若干偶環（包括二元環）不重不漏地覆蓋所有的點。

易證 \(G\) 存在完美匹配當且僅當 \(G\) 存在偶環覆蓋。

  • 如果 \(G\) 存在偶環覆蓋，我們只需要在每個環都隔一條取一條邊，就可以得到一個完美匹配。
  • 如果 \(G\) 存在完美匹配，我們只需要將匹配邊對應的二元環取出，就可以得到一個偶環覆蓋。

然後證明 \(G\) 存在偶環覆蓋當且僅當 \(\tilde{A} \ne 0\)。

考慮行列式的定義

\[ \det A = \sum_{\pi} (-1)^{\pi} \prod_{i} A_{i, \pi_i} \]

其中 \(\pi\) 是任意排列，\((-1)^{\pi}\) 表示若 \(\pi\) 中的逆序對數為奇數，則取 \(-1\)，否則取 \(1\)。

不難看出每個排列都可以被看作 \(G\) 的一個環覆蓋。如果這個環覆蓋中存在奇環，則將這個環翻轉後的和一定為 \(0\)，因此只有偶環覆蓋才能使行列式不為 \(0\)，證畢。

定理：\(\operatorname{rank}\tilde{A}\) 一定為偶數，並且 \(G\) 的最大匹配的大小等於 \(\operatorname{rank}\tilde{A}\) 的一半。

證明

反對稱矩陣的秩只能是偶數；後者請讀者自行思考。

實際應用中不可能帶着 \(|E|\) 個變量進行計算，不過可以取一個數域，例如取某個素數 \(p\) 的剩餘系 \(\mathcal{Z}_p\)，將變量分別隨機替換為 \(\mathcal{Z}_p\) 中的數，再進行計算。方便起見，在無歧義的情況下，以下用 \(\tilde{A}\) 直接指代替換後的矩陣。

定理：\(\operatorname{rank}\tilde{A}\) 至多為 \(G\) 的最大匹配大小的兩倍，並且二者相等的概率至少為 \(1 - \frac n p\)。

考慮到一般圖最大匹配中 \(n\) 基本不會超過 \(10^3\)，實際中 \(p\) 取 \(10^9\) 數量級的素數就足夠了。

由定理可知，如果只需要求最大匹配數，而無需匹配方案，那麼只需要用一次高斯消元求出 \(\operatorname{rank}\tilde{A}\) 即可，遠比帶花樹簡潔。不過如果需要輸出方案，會稍微複雜一些，需要用到下面介紹的算法。

構造完美匹配

由 Tutte 定理和上面的定理可知，如果 \(G\) 存在完美匹配，那麼 \(\tilde{A}\) 有很大概率滿秩。方便起見，以下敍述中均省略「有很大概率」。

記 \(G\) 中標號為 \(i\) 的點為 \(v_i\)，進一步地我們有如下定理：

定理：\(\tilde{A}^{-1}_{j,i} \ne 0 \iff G - \{v_i, v_j\}\) 有完美匹配。

逆矩陣與伴隨矩陣

對任意 \(n\) 階方陣 \(A\)，定義其伴隨矩陣為 \(A^*_{i, j} = (-1)^{i + j} M_{j, i}\)，其中 \(M_{j, i}\) 為刪去第 \(j\) 行第 \(i\) 列的餘子式。換言之，設 \(A\) 的代數餘子式矩陣為 \(M\)，則 \(A^* = M^T\)。

定理：如果 \(A\) 可逆，那麼 \(A^{-1} = \frac 1 {\det A} A^*\)。

所以這裏的 \(A^{-1}_{j, i} \ne 0 \iff M_{i, j} \ne 0\)，也就是 \(A\) 刪去第 \(i\) 行第 \(j\) 列後的部分滿秩。

換言之，如果 \((v_i, v_j) \in E\)，並且 \(\tilde{A}^{-1}_{j, i} \ne 0\)，就表明存在一個完美匹配方案包含 \((v_i, v_j)\) 這條邊。以下將這種邊稱為 可行邊。

由如上定理，對於一個有完美匹配的無向圖 \(G\)，我們可以得到一個比較顯然的暴力算法來尋找一組完美匹配：每次枚舉 \(i, j\)，如果 \((v_i, v_j)\) 是一條可行邊（連邊存在，並且 \(\tilde{A}^{-1}_{j, i} \ne 0\)），就將 \((v_i, v_j)\) 加入匹配方案，並在 \(G\) 中都刪掉這兩個點，再重新計算新的 \(\tilde{A}^{-1}\)。

總共要做 \(\frac n 2\) 輪，每輪都是 \(O(n^3)\) 的，總的複雜度是 \(O(n ^ 4)\)，有點慢了。實際上我們在重新計算 \(\tilde{A}^{-1}\) 時，不必每次都重新用高斯消元求逆矩陣，而是可以利用如下定理：

定理（消去定理）：令

\[ A = \begin{bmatrix} a_{1, 1} & v^T \\ u & B \end{bmatrix} \quad A^{-1} = \begin{bmatrix} \hat a^{1, 1} & \hat v^T \\ \hat u & \hat B \end{bmatrix} \]

並且 \(\hat a_{1, 1} \ne 0\), 那麼就有

\[ B^{-1} = \hat B - \frac {\hat u \hat v^T} {\hat a_{1, 1}} \]

定理中描述的是消去第一行第一列的情況。實際上，它可以非常顯然地推廣到消去任意一行一列的情況，因此我們只需在算法最開始計算一次 \(\tilde{A}^{-1}\)，後面每次刪除兩個點時，只需執行兩次 \(O(n^2)\) 的消去過程即可。

描述有些抽象，可以參考 C++ 代碼 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
void eliminate(int A[][MAXN], int r, int c) {  // 消去第 r 行第 c 列
  row_marked[r] = col_marked[c] = true;        // 已經被消掉

  int inv = quick_power(A[r][c], p - 2);  // 逆元

  for (int i = 1; i <= n; i++)
    if (!row_marked[i] && A[i][c]) {
      int tmp = (long long)A[i][c] * inv % p;

      for (int j = 1; j <= n; j++)
        if (!col_marked[j] && A[r][j])
          A[i][j] = (A[i][j] - (long long)tmp * A[r][j]) % p;
    }
}

總共要做 \(\frac n 2\) 輪，每輪複雜度為 \(O(n^2)\)，因此上述算法可以在 \(O(n^3)\) 的時間內找到一組完美匹配。

構造最大匹配

我們剛剛已經解決了構造一組完美匹配的問題，但是求解問題時一般需要最大匹配。

前面已經提到，\(G\) 的最大匹配大小等於 \(\operatorname{rank}\tilde{A}\) 的一半。如果我們能找到 \(\tilde{A}\) 的一個極大滿秩子矩陣，那麼對子矩陣對應的導出子圖求出一組完美匹配，即可找到 \(G\) 的一組完美匹配。

換一個角度考慮，如果 \(G\) 有完美匹配，那麼 \(\tilde{A}\) 滿秩，換言之，\(\tilde{A}\) 是線性無關的。那麼如果 \(\tilde{A}\) 不是滿秩的，我們可以求出 \(\tilde{A}\) 的一組線性基，然後只保留線性基對應的行列，就可以得到 \(\tilde{A}\) 的一個極大滿秩子矩陣。

求出極大滿秩子矩陣之後，再用上面的算法找出導出子圖的一組完美匹配，即可得到原圖的一組最大匹配。注意由於高斯消元中可能會有行的交換，因此實現時要注意維護好點的編號。

UOJ #79. 一般圖最大匹配 
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
#include <bits/stdc++.h>
using namespace std;

const int maxn = 505, p = (int)1e9 + 7;

int qpow(int a, int b) {
  int ans = 1;
  while (b) {
    if (b & 1) ans = (long long)ans * a % p;
    a = (long long)a * a % p;
    b >>= 1;
  }
  return ans;
}

int A[maxn][maxn], B[maxn][maxn], t[maxn][maxn], id[maxn];

// 高斯消元 O(n^3)
// 在传入 B 时表示计算逆矩阵, 传入 nullptr 则只需计算矩阵的秩
void Gauss(int A[][maxn], int B[][maxn], int n) {
  if (B) {
    memset(B, 0, sizeof(t));
    for (int i = 1; i <= n; i++) B[i][i] = 1;
  }

  for (int i = 1; i <= n; i++) {
    if (!A[i][i]) {
      for (int j = i + 1; j <= n; j++)
        if (A[j][i]) {
          swap(id[i], id[j]);
          for (int k = i; k <= n; k++) swap(A[i][k], A[j][k]);

          if (B)
            for (int k = 1; k <= n; k++) swap(B[i][k], B[j][k]);
          break;
        }

      if (!A[i][i]) continue;
    }

    int inv = qpow(A[i][i], p - 2);

    for (int j = 1; j <= n; j++)
      if (i != j && A[j][i]) {
        int t = (long long)A[j][i] * inv % p;

        for (int k = i; k <= n; k++)
          if (A[i][k]) A[j][k] = (A[j][k] - (long long)t * A[i][k]) % p;

        if (B) {
          for (int k = 1; k <= n; k++)
            if (B[i][k]) B[j][k] = (B[j][k] - (long long)t * B[i][k]) % p;
        }
      }
  }

  if (B)
    for (int i = 1; i <= n; i++) {
      int inv = qpow(A[i][i], p - 2);

      for (int j = 1; j <= n; j++)
        if (B[i][j]) B[i][j] = (long long)B[i][j] * inv % p;
    }
}

bool row_marked[maxn] = {false}, col_marked[maxn] = {false};

int sub_n;  // 极大满秩子矩阵的大小

// 消去一行一列 O(n^2)
void eliminate(int r, int c) {
  row_marked[r] = col_marked[c] = true;  // 已经被消掉

  int inv = qpow(B[r][c], p - 2);

  for (int i = 1; i <= sub_n; i++)
    if (!row_marked[i] && B[i][c]) {
      int t = (long long)B[i][c] * inv % p;

      for (int j = 1; j <= sub_n; j++)
        if (!col_marked[j] && B[r][j])
          B[i][j] = (B[i][j] - (long long)t * B[r][j]) % p;
    }
}

int vertices[maxn], girl[maxn];  // girl 是匹配点, 用来输出方案

int main() {
  auto rng = mt19937(chrono::steady_clock::now().time_since_epoch().count());

  int n, m;
  scanf("%d%d", &n, &m);  // 点数和边数

  while (m--) {
    int x, y;
    scanf("%d%d", &x, &y);
    A[x][y] = rng() % p;
    A[y][x] = -A[x][y];  // Tutte 矩阵
  }

  for (int i = 1; i <= n; i++)
    id[i] = i;  // 输出方案用的，因为高斯消元的时候会交换列
  memcpy(t, A, sizeof(t));

  Gauss(A, nullptr, n);

  for (int i = 1; i <= n; i++)
    if (A[id[i]][id[i]]) vertices[++sub_n] = i;  // 找出一个极大满秩子矩阵

  for (int i = 1; i <= sub_n; i++)
    for (int j = 1; j <= sub_n; j++) A[i][j] = t[vertices[i]][vertices[j]];

  Gauss(A, B, sub_n);

  for (int i = 1; i <= sub_n; i++)
    if (!girl[vertices[i]])
      for (int j = i + 1; j <= sub_n; j++)
        if (!girl[vertices[j]] && t[vertices[i]][vertices[j]] && B[j][i]) {
          // 注意上面那句 if 的写法, 现在 t 是邻接矩阵的备份，
          // 逆矩阵 j 行 i 列不为 0 当且仅当这条边可行
          girl[vertices[i]] = vertices[j];
          girl[vertices[j]] = vertices[i];

          eliminate(i, j);
          eliminate(j, i);
          break;
        }

  printf("%d\n", sub_n / 2);
  for (int i = 1; i <= n; i++) printf("%d ", girl[i]);

  return 0;
}

習題

參考資料

  1. Mucha M, Sankowski P.Maximum matchings via Gaussian elimination
  2. 周子鑫，楊家齊《基於線性代數的一般圖匹配》
  3. ZYQN 《基於線性代數的一般圖匹配算法》
本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：H-J-Granger, accelsao, Ir1d, Early0v0, Henry-ZHR, HeliumOI, AntiLeaf
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用