AIOJ

The maximum-weight matching of a bipartite graph is a matching whose sum of edge weights is as large as possible.

Hungarian Algorithm (Kuhn–Munkres Algorithm)

The Hungarian algorithm, also called the KM algorithm, can find a maximum-weight perfect matching in a bipartite graph in $O(n^3)$ time.

Because the two vertex sets of a bipartite graph do not always have the same size, to apply the KM algorithm to maximum-weight matching, first add dummy vertices to the smaller set until the two sides have equal size, and set the weights of nonexistent edges to $0$ . The problem is then converted into a maximum-weight perfect matching problem, which can be solved with the KM algorithm.

Feasible Vertex Labels

Assign every vertex $i$ a weight $l(i)$ such that every edge $(u,v)$ satisfies $w(u,v) \leq l(u) + l(v)$ .

Equality Subgraph

Under a set of feasible vertex labels, the spanning subgraph of the original graph contains all vertices but only the edges satisfying $w(u,v) = l(u) + l(v)$ , namely the edges $(u,v)$ .

Theorem 1: For a set of feasible vertex labels, if its equality subgraph has a perfect matching, that matching is a maximum-weight perfect matching of the original bipartite graph.

Proof 1.

Consider any perfect matching $M$ of the original bipartite graph. Its sum of edge weights is

$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)$

For any perfect matching $M'$ in the equality subgraph of a set of feasible vertex labels, the sum of edge weights is

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

Thus, the weight of any perfect matching is no greater than $val(M')$ , so $M'$ is a maximum-weight matching.

Given Theorem 1, our goal is to repeatedly adjust the feasible vertex labels until the equality subgraph has a perfect matching.

Since the two sides have equal size, let their common size be $n$ . Let $lx(i)$ denote the label of the $i$ -th vertex on the left, $ly(i)$ the label of the $i$ -th vertex on the right, and $w(u,v)$ the weight of the edge between the $u$ -th left vertex and the $v$ -th right vertex.

First initialize a set of feasible vertex labels, for example:

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

Then choose an unmatched vertex and find an augmenting path as in maximum matching. If an augmenting path is found, augment the matching; otherwise, an alternating tree is obtained.

Let $S$ and $T$ denote the left-side and right-side vertices in the alternating tree, and let $S'$ and $T'$ denote the vertices not in the alternating tree.

bigraph-weight-match-1

In the equality subgraph:

The edge $S-T'$ does not exist; otherwise the alternating tree would grow.
$S'-T$ must be a non-matching edge; otherwise its left endpoint would belong to $S$ .

Suppose that the label change gives $S$ the amount $-a$ and $T$ the amount $+a$ . We can see that

Edges $S-T$ remain in the equality subgraph.
$S'-T'$ is unchanged.
On $S-T'$ , $lx + ly$ decreases, so some of these edges may enter the equality subgraph.
On $S'-T$ , $lx + ly$ increases, so none of these edges can enter the equality subgraph.

Therefore, the value of $a$ should clearly be the smallest slack among $S-T'$ :

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

When a new edge $(u,v)$ enters the equality subgraph, there are two cases:

If $v$ is unmatched, an augmenting path has been found.
If $v$ is matched, it is matched with a vertex in $S'$ .

Thus, an augmenting path can be found after changing the vertex labels at most $n$ times.

Whenever the vertex labels are changed, the edges in the alternating tree do not leave the equality subgraph, so we can maintain this tree directly.

For each vertex in $T$ , denoted $v$ , maintain

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

$O(n)$ time therefore suffices to compute the value $a$ used to change the vertex labels:

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

When a new vertex enters $S$ as the alternating tree grows, it takes $O(n)$ to update $slack(v)$ . Changing the vertex labels costs $O(n)$ to subtract from every $slack(v)$ the value $a$ . As soon as the alternating tree reaches an unmatched vertex, an augmenting path is found.

At the beginning, we enumerate the $n$ vertices to find an augmenting path. Finding one augmenting path may require extending the alternating tree $n$ times, and each extension requires $n$ maintenance operations, for a total of $O(n^3)$ .

Reference Code

template <typename T>
struct hungarian {  // km
  int n;
  vector<int> matchx;  // the matching vertex in the left set
  vector<int> matchy;  // the matching vertex in the right set
  vector<int> pre;     // the left vertex connected to a right vertex
  vector<bool> visx;   // visited array, left
  vector<bool> visy;   // visited array, right
  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);  // a negative value is no better than not matching, so setting it to 0 does not affect the result
  }

  bool check(int v) {
    visy[v] = true;
    if (matchy[v] != -1) {
      q.push(matchy[v]);
      visx[matchy[v]] = true;  // in S
      return false;
    }
    // found a new unmatched vertex; update the matching, and let the pre array record the vertex connected to it by a non-matching edge
    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 means the edge may enter the equality subgraph; find an augmenting path
                                      // found one, so return and rebuild the alternating tree
                return;
              }
            }
          }
        }
      }
      // no augmenting path; update the vertex labels
      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() {
    // initial vertex labels
    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 original paper is The Dynamic Hungarian Algorithm for the Assignment Problem with Changing Costs.

A paper with clearer pseudocode is A Fast Dynamic Assignment Algorithm for Solving Resource Allocation Problems.

The related OJ problem is DAP.

Algorithm Idea

Change the weights between a single vertex $u_i$ $u_{i}$ and all $v_j$ $v_{j}$ , that is, a row of the weight matrix.
- Update the vertex label $lx(u_i) = max(w_{ij} - v_{j}), \forall j$ .
- Delete the matching associated with $u_i$ .
Change the weights between all $u_i$ $u_{i}$ and a single vertex $v_j$ $v_{j}$ , that is, a column of the weight matrix.
- Update the vertex label $ly(v_j) = max(w_{ij} - u_{i}), \forall i$ .
- Delete the matching associated with $v_j$ .
Change the weight between a single vertex $u_i$ $u_{i}$ and a single vertex $v_j$ $v_{j}$ , that is, one element of the weight matrix.
- Either operation 1 or operation 2 is sufficient.
Add a single vertex $u_i$ $u_{i}$ or a single vertex $v_j$ $v_{j}$ , that is, add or delete a row or column in the weight matrix.
- Perform operation 1 or operation 2 accordingly. Note that adding a vertex only adds the vertex; no additional weight is specified, and the weights between the new vertex and all other vertices are 0.

Algorithm Proof

Let the original graph be G, let the labels on the two sides be $\alpha^{i}$ and $\beta^{j}$ , and let l be a feasible vertex labeling. Then $G_l$ is a subgraph of G containing the vertices and edges of G that satisfy $w_{ij} = alpha_{i}+beta_{j}$ .
Theorem 1 in the Hungarian algorithm section above proves that, for a set of feasible vertex labels, if its equality subgraph has a perfect matching, that matching is a maximum-weight perfect matching of the original bipartite graph.
Suppose the original optimal matching is $M^*$ $M^{*}$ . When a change occurs, we update the feasible vertex labels according to the rules, and the updated label is denoted by $\alpha^{i^*}$ $α^{i^{*}}$ or $\beta^{j^*}$ $β^{j^{*}}$ . The following cases may occur:
1. An entire row of the weight matrix is changed. Let the changed row be row $i^*$ , meaning that all edges incident to $v_{i^*}$ are changed. The original label of $v_{i^*}$ may no longer satisfy the condition because we need $w_{i^{*}j} \leq alpha_{i^*}+beta_{j}$ . For the other $u_j$ , however, the edge weights are unchanged except for edges related to $i^*$ , so their labels remain valid. The algorithm therefore changes the label related to $v_{i^*}$ so that the labels are feasible.
2. An entire column of the weight matrix is changed. The same reasoning shows that changing the labels makes them feasible.
3. One element of the weight matrix is changed. Changing either one of the two associated labels is enough to satisfy the labeling condition.
Every change to the weight matrix involves a particular vertex, which may be on the left or the right. Denote it by $x$ , and suppose this vertex is matched with some vertex $y$ in the original optimal matching. Each change can unpair at most this pair of vertices, so one search round of the Hungarian algorithm is enough to obtain a new match. By Theorem 1, the new match is optimal.

The following code is presumably the code submitted by the authors of paper 2 (the code below is the maximum-weight version, while the original paper uses minimum cost).

Reference Code for the Dynamic Hungarian Algorithm

#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;
}

Reduction to a Min-Cost Flow Model

Similar to maximum bipartite matching, maximum-weight matching in a bipartite graph can also be solved as a network flow problem.

First, add a source and a sink to the graph.

Add an edge of capacity $1$ and cost $0$ from the source to every left-side vertex of the bipartite graph, and an edge of capacity $1$ and cost $0$ from every right-side vertex to the sink.

Next, for every edge connecting a left-side vertex $u$ and a right-side vertex $v$ with weight $w$ , add an edge from $u$ to $v$ with capacity $1$ and cost $w$ .

In addition, because the number of matching edges in a maximum-weight matching need not equal the number of matching edges in a maximum-cardinality matching, add an edge of capacity $1$ and cost $0$ from every left-side vertex directly to the sink.

The answer is obtained by finding the maximum-cost maximum flow of this network. The maximum flow is necessarily equal to the number of left-side vertices, and the maximum cost among flows of that value corresponds to a maximum-weight matching.

Problems

UOJ #80. Maximum-Weight Bipartite Matching

Template problem

#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);  // a negative value is no better than not matching, so setting it to 0 does not affect the result
  }

  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;
              }
            }
          }
        }
      }
      // no augmenting path; update the vertex labels
      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() {
    // initial vertex labels
    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();
}