跳转至

二分圖最大權匹配

二分圖的最大權匹配是指二分圖中邊權和最大的匹配。

Hungarian Algorithm（Kuhn–Munkres Algorithm）

匈牙利算法又稱為 KM 算法，可以在 \(O(n^3)\) 時間內求出二分圖的 最大權完美匹配。

考慮到二分圖中兩個集合中的點並不總是相同，為了能應用 KM 算法解決二分圖的最大權匹配，需要先作如下處理：將兩個集合中點數比較少的補點，使得兩邊點數相同，再將不存在的邊權重設為 \(0\)，這種情況下，問題就轉換成求 最大權完美匹配問題，從而能應用 KM 算法求解。

可行頂標

給每個節點 \(i\) 分配一個權值 \(l(i)\)，對於所有邊 \((u,v)\) 滿足 \(w(u,v) \leq l(u) + l(v)\)。

相等子圖

在一組可行頂標下原圖的生成子圖，包含所有點但只包含滿足 \(w(u,v) = l(u) + l(v)\) 的邊 \((u,v)\)。

定理 1 : 對於某組可行頂標，如果其相等子圖存在完美匹配，那麼，該匹配就是原二分圖的最大權完美匹配。

證明 1.

考慮原二分圖任意一組完美匹配 \(M\)，其邊權和為

\(val(M) = \sum_{(u,v)\in M} {w(u,v)} \leq \sum_{(u,v)\in M} {l(u) + l(v)} \leq \sum_{i=1}^{n} l(i)\)

任意一組可行頂標的相等子圖的完美匹配 \(M'\) 的邊權和

\(val(M') = \sum_{(u,v)\in M} {l(u) + l(v)} = \sum_{i=1}^{n} l(i)\)

即任意一組完美匹配的邊權和都不會大於 \(val(M')\)，那個 \(M'\) 就是最大權匹配。

有了定理 1，我們的目標就是透過不斷的調整可行頂標，使得相等子圖是完美匹配。

因為兩邊點數相等，假設點數為 \(n\)，\(lx(i)\) 表示左邊第 \(i\) 個點的頂標，\(ly(i)\) 表示右邊第 \(i\) 個點的頂標，\(w(u,v)\) 表示左邊第 \(u\) 個點和右邊第 \(v\) 個點之間的權重。

首先初始化一組可行頂標，例如

\(lx(i) = \max_{1\leq j\leq n} \{ w(i, j)\},\, ly(i) = 0\)

然後選一個未匹配點，如同最大匹配一樣求增廣路。找到增廣路就增廣，否則，會得到一個交錯樹。

令 \(S\)，\(T\) 表示二分圖左邊右邊在交錯樹中的點，\(S'\)，\(T'\) 表示不在交錯樹中的點。

bigraph-weight-match-1

在相等子圖中：

  • \(S-T'\) 的邊不存在，否則交錯樹會增長。
  • \(S'-T\) 一定是非匹配邊，否則他就屬於 \(S\)。

假設給 \(S\) 中的頂標 \(-a\)，給 \(T\) 中的頂標 \(+a\)，可以發現

  • \(S-T\) 邊依然存在相等子圖中。
  • \(S'-T'\) 沒變化。
  • \(S-T'\) 中的 \(lx + ly\) 有所減少，可能加入相等子圖。
  • \(S'-T\) 中的 \(lx + ly\) 會增加，所以不可能加入相等子圖。

所以這個 \(a\) 值的選擇，顯然得是 \(S-T'\) 當中最小的邊權，

\(a = \min \{ lx(u) + ly(v) - w(u,v) | u\in{S} , v\in{T'} \}\)。

當一條新的邊 \((u,v)\) 加入相等子圖後有兩種情況

  • \(v\) 是未匹配點，則找到增廣路
  • \(v\) 和 \(S'\) 中的點已經匹配

這樣至多修改 \(n\) 次頂標後，就可以找到增廣路。

每次修改頂標的時候，交錯樹中的邊不會離開相等子圖，那麼我們直接維護這棵樹。

我們對 \(T\) 中的每個點 \(v\) 維護

\(slack(v) = \min \{ lx(u) + ly(v) - w(u,v) | u\in{S} \}\)。

所以可以在 \(O(n)\) 算出頂標修改值 \(a\)

\(a = \min \{ slack(v) | v\in{T'} \}\)

交錯樹新增一個點進入 \(S\) 的時候需要 \(O(n)\) 更新 \(slack(v)\)。修改頂標需要 \(O(n)\) 給每個 \(slack(v)\) 減去 \(a\)。只要交錯樹找到一個未匹配點，就找到增廣路。

一開始枚舉 \(n\) 個點找增廣路，為了找增廣路需要延伸 \(n\) 次交錯樹，每次延伸需要 \(n\) 次維護，共 \(O(n^3)\)。

參考代碼 
  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
template <typename T>
struct hungarian {  // km
  int n;
  vector<int> matchx;  // 左集合對應的匹配點
  vector<int> matchy;  // 右集合對應的匹配點
  vector<int> pre;     // 連接右集合的左點
  vector<bool> visx;   // 拜訪數組 左
  vector<bool> visy;   // 拜訪數組 右
  vector<T> lx;
  vector<T> ly;
  vector<vector<T> > g;
  vector<T> slack;
  T inf;
  T res;
  queue<int> q;
  int org_n;
  int org_m;

  hungarian(int _n, int _m) {
    org_n = _n;
    org_m = _m;
    n = max(_n, _m);
    inf = numeric_limits<T>::max();
    res = 0;
    g = vector<vector<T> >(n, vector<T>(n));
    matchx = vector<int>(n, -1);
    matchy = vector<int>(n, -1);
    pre = vector<int>(n);
    visx = vector<bool>(n);
    visy = vector<bool>(n);
    lx = vector<T>(n, -inf);
    ly = vector<T>(n);
    slack = vector<T>(n);
  }

  void addEdge(int u, int v, int w) {
    g[u][v] = max(w, 0);  // 負值還不如不匹配 因此設為0不影響
  }

  bool check(int v) {
    visy[v] = true;
    if (matchy[v] != -1) {
      q.push(matchy[v]);
      visx[matchy[v]] = true;  // in S
      return false;
    }
    // 找到新的未匹配點 更新匹配點 pre 數組記錄着"非匹配邊"上與之相連的點
    while (v != -1) {
      matchy[v] = pre[v];
      swap(v, matchx[pre[v]]);
    }
    return true;
  }

  void bfs(int i) {
    while (!q.empty()) {
      q.pop();
    }
    q.push(i);
    visx[i] = true;
    while (true) {
      while (!q.empty()) {
        int u = q.front();
        q.pop();
        for (int v = 0; v < n; v++) {
          if (!visy[v]) {
            T delta = lx[u] + ly[v] - g[u][v];
            if (slack[v] >= delta) {
              pre[v] = u;
              if (delta) {
                slack[v] = delta;
              } else if (check(v)) {  // delta=0 代表有機會加入相等子圖 找增廣路
                                      // 找到就return 重建交錯樹
                return;
              }
            }
          }
        }
      }
      // 沒有增廣路 修改頂標
      T a = inf;
      for (int j = 0; j < n; j++) {
        if (!visy[j]) {
          a = min(a, slack[j]);
        }
      }
      for (int j = 0; j < n; j++) {
        if (visx[j]) {  // S
          lx[j] -= a;
        }
        if (visy[j]) {  // T
          ly[j] += a;
        } else {  // T'
          slack[j] -= a;
        }
      }
      for (int j = 0; j < n; j++) {
        if (!visy[j] && slack[j] == 0 && check(j)) {
          return;
        }
      }
    }
  }

  void solve() {
    // 初始頂標
    for (int i = 0; i < n; i++) {
      for (int j = 0; j < n; j++) {
        lx[i] = max(lx[i], g[i][j]);
      }
    }

    for (int i = 0; i < n; i++) {
      fill(slack.begin(), slack.end(), inf);
      fill(visx.begin(), visx.end(), false);
      fill(visy.begin(), visy.end(), false);
      bfs(i);
    }

    // custom
    for (int i = 0; i < n; i++) {
      if (g[i][matchx[i]] > 0) {
        res += g[i][matchx[i]];
      } else {
        matchx[i] = -1;
      }
    }
    cout << res << "\n";
    for (int i = 0; i < org_n; i++) {
      cout << matchx[i] + 1 << " ";
    }
    cout << "\n";
  }
};

Dynamic Hungarian Algorithm

原論文 The Dynamic Hungarian Algorithm for the Assignment Problem with Changing Costs

偽代碼更清晰的論文 A Fast Dynamic Assignment Algorithm for Solving Resource Allocation Problems

相關 OJ 問題 DAP

算法思路
  1. 修改單點 \(u_i\) 和所有 \(v_j\) 之間的權重，即權重矩陣中的一行
    • 修改頂標 \(lx(u_i) = max(w_{ij} - v_{j}), \forall j\)
    • 刪除 \(u_i\) 相關的匹配
  2. 修改所有 \(u_i\) 和單點 \(v_j\) 之間的權重，即權重矩陣中的一列
    • 修改頂標 \(ly(v_j) = max(w_{ij} - u_{i}), \forall i\)
    • 刪除 \(v_j\) 相關的匹配
  3. 修改單點 \(u_i\) 和單點 \(v_j\) 之間的權重，即權重矩陣中的單個元素
    • 做 1 或 2 兩種操作之一即可
  4. 添加某一單點 \(u_i\)，或者某一單點 \(v_j\)，即在權重矩陣中添加或者刪除一行或者一列
    • 對應地做 1 或 2 即可，注意此處加點操作僅為加點，不額外設定權重值，新加點與其他點的權重為 0.
算法證明
  • 設原圖為 G，左右兩邊的頂標為 \(\alpha^{i}\) 和 \(\beta^{j}\)，可行頂標為 l，那 \(G_l\) 是 G 的一個子圖，包含圖 G 中滿足 \(w_{ij} = alpha_{i}+beta_{j}\) 的點和邊。
  • 在上面匈牙利算法的部分，定理一證明了：對於某組可行頂標，如果其相等子圖存在完美匹配，那麼，該匹配就是原二分圖的最大權完美匹配。
  • 假設原來的最優匹配是 \(M^*\), 當一個修改發生的時候，我們會根據規則更新可行頂標，更新後的頂標設為 \(\alpha^{i^*}\), 或者 \(\beta^{j^*}\)，會出現以下情況：
    1. 權重矩陣的一整行被修改了，設被修改的行為 \(i^*\) 行，即 \(v_{i^*}\) 的所有邊被修改了，所以 \(v_{i^*}\) 原來的頂標可能不滿足條件，因為我們需要 \(w_{i^{*}j} \leq alpha_{i^*}+beta_{j}\)，但對於其他的 \(u_j\) 來説，除了 \(i^*\) 相關的邊，他們的邊權是不變的，因此他們的頂標都是合法的，所以算法中修改了 \(v_{i^*}\) 相關的頂標使得這組頂標是一組可行頂標。
    2. 權重矩陣的一整列被修改了，同理可得算法修改頂標使得這組頂標是一組可行頂標。
    3. 修改權重矩陣某一元素，任意修改其中一個頂標即可滿足頂標條件
  • 每一次權重矩陣被修改，都關係到一個特定節點，這個節點可能是左邊的也可能是右邊的，因此我們直接記為 \(x\), 這個節點和某個節點 \(y\) 在原來的最優匹配中匹配上了。每一次修改操作，最多讓這一對節點 unpair，因此我們只要跑一輪匈牙利算法中的搜索我們就能得到一個新的 match，而根據定理一，新跑出來的 match 是最優的。

以下代碼應該為論文 2 作者提交的代碼（以下代碼為最大化權重版本，原始論文中為最小化 cost）

動態匈牙利算法參考代碼 
  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
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <list>
using namespace std;
typedef long long LL;
const LL INF = (LL)1e15;
const int MAXV = 105;

int N, mateS[MAXV], mateT[MAXV], p[MAXV];
LL u[MAXV], v[MAXV], slack[MAXV];
LL W[MAXV][MAXV];
bool m[MAXV];
list<int> Q;

void readMatrix() {
  scanf("%d", &N);
  for (int i = 0; i < N; i++)
    for (int j = 0; j < N; j++) scanf("%lld", &W[i][j]);
}

void initHungarian() {
  memset(mateS, -1, sizeof(mateS));
  memset(mateT, -1, sizeof(mateT));
  for (int i = 0; i < N; i++) {
    u[i] = -INF;
    for (int j = 0; j < N; j++) u[i] = max(u[i], W[i][j]);
    v[i] = 0;
  }
}

void augment(int j) {
  int i, next;
  do {
    i = p[j];
    mateT[j] = i;
    next = mateS[i];
    mateS[i] = j;
    if (next != -1) j = next;
  } while (next != -1);
}

LL hungarian() {
  int nres = 0;
  for (int i = 0; i < N; i++)
    if (mateS[i] == -1) nres++;

  while (nres > 0) {
    for (int i = 0; i < N; i++) {
      m[i] = false;
      p[i] = -1;
      slack[i] = INF;
    }
    bool aug = false;
    Q.clear();
    for (int i = 0; i < N; i++)
      if (mateS[i] == -1) {
        Q.push_back(i);
        break;
      }

    do {
      int i, j;
      i = Q.front();
      Q.pop_front();
      m[i] = true;
      j = 0;

      while (aug == false && j < N) {
        if (mateS[i] != j) {
          LL minSlack = u[i] + v[j] - W[i][j];
          if (minSlack < slack[j]) {
            slack[j] = minSlack;
            p[j] = i;
            if (slack[j] == 0) {
              if (mateT[j] == -1) {
                augment(j);
                aug = true;
                nres--;
              } else
                Q.push_back(mateT[j]);
            }
          }
        }
        j++;
      }

      if (aug == false && Q.size() == 0) {
        LL minSlack = INF;
        for (int k = 0; k < N; k++)
          if (slack[k] > 0) minSlack = min(minSlack, slack[k]);
        for (int k = 0; k < N; k++)
          if (m[k]) u[k] -= minSlack;

        int x = -1;
        bool X[MAXV];
        for (int k = 0; k < N; k++)
          if (slack[k] == 0)
            v[k] += minSlack;
          else {
            slack[k] -= minSlack;
            if (slack[k] == 0 && mateT[k] == -1) x = k;
            if (slack[k] == 0)
              X[k] = true;
            else
              X[k] = false;
          }

        if (x == -1) {
          for (int k = 0; k < N; k++)
            if (X[k]) Q.push_back(mateT[k]);
        } else {
          augment(x);
          aug = true;
          nres--;
        }
      }
    } while (aug == false);
  }

  LL ans = 0;
  for (int i = 0; i < N; i++) ans += (u[i] + v[i]);
  return ans;
}

void dynamicHungarian() {
  char type[2];
  LL w;
  int i, j;

  scanf("%s", type);
  if (type[0] == 'C') {
    scanf("%d%d%lld", &i, &j, &w);
    if ((w < W[i][j]) && (mateS[i] == j)) {
      W[i][j] = w;
      if (mateS[i] != -1) {
        mateT[mateS[i]] = -1;
        mateS[i] = -1;
      }
    } else if ((w > W[i][j]) && (u[i] + v[j] < w)) {
      W[i][j] = w;
      u[i] = -INF;
      for (int c = 0; c < N; c++) u[i] = max(u[i], W[i][c] - v[c]);
      if (mateS[i] != j) {
        mateT[mateS[i]] = -1;
        mateS[i] = -1;
      }
    } else
      W[i][j] = w;
  } else if (type[0] == 'X') {
    scanf("%d", &i);
    for (int c = 0; c < N; c++) scanf("%lld", &W[i][c]);
    if (mateS[i] != -1) {
      mateT[mateS[i]] = -1;
      mateS[i] = -1;
    }
    u[i] = -INF;
    for (int c = 0; c < N; c++) u[i] = max(u[i], W[i][c] - v[c]);
  } else if (type[0] == 'Y') {
    scanf("%d", &j);
    for (int r = 0; r < N; r++) scanf("%lld", &W[r][j]);
    if (mateT[j] != -1) {
      mateS[mateT[j]] = -1;
      mateT[j] = -1;
    }
    v[j] = -INF;
    for (int r = 0; r < N; r++) v[j] = max(v[j], W[r][j] - u[r]);
  } else if (type[0] == 'A') {
    i = j = N++;
    u[i] = -INF;
    for (int c = 0; c < N; c++) u[i] = max(u[i], W[i][c] - v[c]);
    v[j] = -INF;
    for (int r = 0; r < N; r++) v[j] = max(v[j], W[r][j] - u[r]);
  } else if (type[0] == 'Q')
    printf("%lld\n", hungarian());
}

int main() {
  readMatrix();
  initHungarian();
  LL ans = hungarian();
  int M;
  scanf("%d", &M);
  while (M--) dynamicHungarian();
  return 0;
}

轉化為費用流模型

在圖中新增一個源點和一個匯點。

從源點向二分圖的每個左部點連一條流量為 \(1\)，費用為 \(0\) 的邊，從二分圖的每個右部點向匯點連一條流量為 \(1\)，費用為 \(0\) 的邊。

接下來對於二分圖中每一條連接左部點 \(u\) 和右部點 \(v\)，邊權為 \(w\) 的邊，則連一條從 \(u\) 到 \(v\)，流量為 \(1\)，費用為 \(w\) 的邊。

求這個網絡的 最大費用最大流 即可得到答案。

習題

UOJ #80. 二分圖最大權匹配

模板題

  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
#include <bits/stdc++.h>
using namespace std;

template <typename T>
struct hungarian {  // km
  int n;
  vector<int> matchx;
  vector<int> matchy;
  vector<int> pre;
  vector<bool> visx;
  vector<bool> visy;
  vector<T> lx;
  vector<T> ly;
  vector<vector<T> > g;
  vector<T> slack;
  T inf;
  T res;
  queue<int> q;
  int org_n;
  int org_m;

  hungarian(int _n, int _m) {
    org_n = _n;
    org_m = _m;
    n = max(_n, _m);
    inf = numeric_limits<T>::max();
    res = 0;
    g = vector<vector<T> >(n, vector<T>(n));
    matchx = vector<int>(n, -1);
    matchy = vector<int>(n, -1);
    pre = vector<int>(n);
    visx = vector<bool>(n);
    visy = vector<bool>(n);
    lx = vector<T>(n, -inf);
    ly = vector<T>(n);
    slack = vector<T>(n);
  }

  void addEdge(int u, int v, int w) {
    g[u][v] = max(w, 0);  // 负值还不如不匹配 因此设为0不影响
  }

  bool check(int v) {
    visy[v] = true;
    if (matchy[v] != -1) {
      q.push(matchy[v]);
      visx[matchy[v]] = true;
      return false;
    }
    while (v != -1) {
      matchy[v] = pre[v];
      swap(v, matchx[pre[v]]);
    }
    return true;
  }

  void bfs(int i) {
    while (!q.empty()) {
      q.pop();
    }
    q.push(i);
    visx[i] = true;
    while (true) {
      while (!q.empty()) {
        int u = q.front();
        q.pop();
        for (int v = 0; v < n; v++) {
          if (!visy[v]) {
            T delta = lx[u] + ly[v] - g[u][v];
            if (slack[v] >= delta) {
              pre[v] = u;
              if (delta) {
                slack[v] = delta;
              } else if (check(v)) {
                return;
              }
            }
          }
        }
      }
      // 没有增广路 修改顶标
      T a = inf;
      for (int j = 0; j < n; j++) {
        if (!visy[j]) {
          a = min(a, slack[j]);
        }
      }
      for (int j = 0; j < n; j++) {
        if (visx[j]) {  // S
          lx[j] -= a;
        }
        if (visy[j]) {  // T
          ly[j] += a;
        } else {  // T'
          slack[j] -= a;
        }
      }
      for (int j = 0; j < n; j++) {
        if (!visy[j] && slack[j] == 0 && check(j)) {
          return;
        }
      }
    }
  }

  void solve() {
    // 初始顶标
    for (int i = 0; i < n; i++) {
      for (int j = 0; j < n; j++) {
        lx[i] = max(lx[i], g[i][j]);
      }
    }

    for (int i = 0; i < n; i++) {
      fill(slack.begin(), slack.end(), inf);
      fill(visx.begin(), visx.end(), false);
      fill(visy.begin(), visy.end(), false);
      bfs(i);
    }

    // custom
    for (int i = 0; i < n; i++) {
      if (g[i][matchx[i]] > 0) {
        res += g[i][matchx[i]];
      } else {
        matchx[i] = -1;
      }
    }
    cout << res << "\n";
    for (int i = 0; i < org_n; i++) {
      cout << matchx[i] + 1 << " ";
    }
    cout << "\n";
  }
};

int main() {
  ios::sync_with_stdio(0), cin.tie(0);
  int n, m, e;
  cin >> n >> m >> e;

  hungarian<long long> solver(n, m);

  int u, v, w;
  for (int i = 0; i < e; i++) {
    cin >> u >> v >> w;
    u--, v--;
    solver.addEdge(u, v, w);
  }
  solver.solve();
}
本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：accelsao, Enter-tainer, guodong2005, StudyingFather, Backl1ght, Chrogeek, H-J-Granger, Henry-ZHR
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用