跳转至

二分图最大权匹配

二分图的最大权匹配是指二分图中边权和最大的匹配。

Hungarian Algorithm（Kuhn–Munkres Algorithm）

匈牙利算法又称为 KM 算法，可以在 时间内求出二分图的 最大权完美匹配。

考虑到二分图中两个集合中的点并不总是相同，为了能应用 KM 算法解决二分图的最大权匹配，需要先作如下处理：将两个集合中点数比较少的补点，使得两边点数相同，再将不存在的边权重设为 ，这种情况下，问题就转换成求 最大权完美匹配问题，从而能应用 KM 算法求解。

可行顶标

给每个节点 分配一个权值 ，对于所有边 满足 。

相等子图

在一组可行顶标下原图的生成子图，包含所有点但只包含满足 的边 。

定理 1 : 对于某组可行顶标，如果其相等子图存在完美匹配，那么，该匹配就是原二分图的最大权完美匹配。

证明 1.

考虑原二分图任意一组完美匹配 ，其边权和为

任意一组可行顶标的相等子图的完美匹配 的边权和

即任意一组完美匹配的边权和都不会大于 ，那个 就是最大权匹配。

有了定理 1，我们的目标就是透过不断的调整可行顶标，使得相等子图是完美匹配。

因为两边点数相等，假设点数为 ， 表示左边第 个点的顶标， 表示右边第 个点的顶标， 表示左边第 个点和右边第 个点之间的权重。

首先初始化一组可行顶标，例如

然后选一个未匹配点，如同最大匹配一样求增广路。找到增广路就增广，否则，会得到一个交错树。

令 ， 表示二分图左边右边在交错树中的点，， 表示不在交错树中的点。

bigraph-weight-match-1

在相等子图中：

  • 的边不存在，否则交错树会增长。
  • 一定是非匹配边，否则他就属于 。

假设给 中的顶标 ，给 中的顶标 ，可以发现

  • 边依然存在相等子图中。
  • 没变化。
  • 中的 有所减少，可能加入相等子图。
  • 中的 会增加，所以不可能加入相等子图。

所以这个 值的选择，显然得是 当中最小的边权，

。

当一条新的边 加入相等子图后有两种情况

  • 是未匹配点，则找到增广路
  • 和 中的点已经匹配

这样至多修改 次顶标后，就可以找到增广路。

每次修改顶标的时候，交错树中的边不会离开相等子图，那么我们直接维护这棵树。

我们对 中的每个点 维护

。

所以可以在 算出顶标修改值

交错树新增一个点进入 的时候需要 更新 。修改顶标需要 给每个 减去 。只要交错树找到一个未匹配点，就找到增广路。

一开始枚举 个点找增广路，为了找增广路需要延伸 次交错树，每次延伸需要 次维护，共 。

参考代码 
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. 修改单点 和所有 之间的权重，即权重矩阵中的一行
    • 修改顶标
    • 删除 相关的匹配
  2. 修改所有 和单点 之间的权重，即权重矩阵中的一列
    • 修改顶标
    • 删除 相关的匹配
  3. 修改单点 和单点 之间的权重，即权重矩阵中的单个元素
    • 做 1 或 2 两种操作之一即可
  4. 添加某一单点 ，或者某一单点 ，即在权重矩阵中添加或者删除一行或者一列
    • 对应地做 1 或 2 即可，注意此处加点操作仅为加点，不额外设定权重值，新加点与其他点的权重为 0.
算法证明
  • 设原图为 G，左右两边的顶标为 和 ，可行顶标为 l，那 是 G 的一个子图，包含图 G 中满足 的点和边。
  • 在上面匈牙利算法的部分，定理一证明了：对于某组可行顶标，如果其相等子图存在完美匹配，那么，该匹配就是原二分图的最大权完美匹配。
  • 假设原来的最优匹配是 , 当一个修改发生的时候，我们会根据规则更新可行顶标，更新后的顶标设为 , 或者 ，会出现以下情况：
    1. 权重矩阵的一整行被修改了，设被修改的行为 行，即 的所有边被修改了，所以 原来的顶标可能不满足条件，因为我们需要 ，但对于其他的 来说，除了 相关的边，他们的边权是不变的，因此他们的顶标都是合法的，所以算法中修改了 相关的顶标使得这组顶标是一组可行顶标。
    2. 权重矩阵的一整列被修改了，同理可得算法修改顶标使得这组顶标是一组可行顶标。
    3. 修改权重矩阵某一元素，任意修改其中一个顶标即可满足顶标条件
  • 每一次权重矩阵被修改，都关系到一个特定节点，这个节点可能是左边的也可能是右边的，因此我们直接记为 , 这个节点和某个节点 在原来的最优匹配中匹配上了。每一次修改操作，最多让这一对节点 unpair，因此我们只要跑一轮匈牙利算法中的搜索我们就能得到一个新的 match，而根据定理一，新跑出来的 match 是最优的。

以下代码应该为论文 2 作者提交的代码（以下代码为最大化权重版本，原始论文中为最小化 cost）

动态匈牙利算法参考代码 
#include <algorithm>
#include <cstring>
#include <iostream>
#include <list>
using namespace std;
using LL = long long;
constexpr LL INF = (LL)1e15;
constexpr 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() {
  cin >> N;
  for (int i = 0; i < N; i++)
    for (int j = 0; j < N; j++) cin >> 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 && 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 && 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);
  }

  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;

  cin >> type;
  if (type[0] == 'C') {
    cin >> 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') {
    cin >> i;
    for (int c = 0; c < N; c++) cin >> 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') {
    cin >> j;
    for (int r = 0; r < N; r++) cin >> 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')
    cout << hungarian() << '\n';
}

int main() {
  cin.tie(nullptr)->sync_with_stdio(false);
  readMatrix();
  initHungarian();
  LL ans = hungarian();
  int M;
  cin >> M;
  while (M--) dynamicHungarian();
  return 0;
}

转化为费用流模型

与 二分图最大匹配 类似，二分图的最大权匹配也可以转化为网络流问题来求解。

首先，在图中新增一个源点和一个汇点。

从源点向二分图的每个左部点连一条流量为 ，费用为 的边，从二分图的每个右部点向汇点连一条流量为 ，费用为 的边。

接下来对于二分图中每一条连接左部点 和右部点 ，边权为 的边，则连一条从 到 ，流量为 ，费用为 的边。

另外，考虑到最大权匹配下，匹配边的数量不一定与最大匹配的匹配边数量相等，因此对于每个左部点，还需向汇点连一条流量为 ，费用为 的边。

求这个网络的 最大费用最大流 即可得到答案。此时，该网络的最大流量一定为左部点的数量，而最大流量下的最大费用即对应一个最大权匹配方案。

习题

UOJ #80. 二分图最大权匹配

模板题

#include <iostream>
#include <limits>
#include <queue>
#include <vector>
using namespace std;

template <typename T>
struct hungarian {  // km
  int n;
  vector<int> matchx, matchy, pre;
  vector<bool> visx, visy;
  vector<T> lx, ly;
  vector<vector<T>> g;
  vector<T> slack;
  T inf, res;
  queue<int> q;
  int org_n, 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(false), cin.tie(nullptr);
  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();
}
最后更新: 2024年10月9日
创建日期: 2020年2月20日
回到页面顶部