Skip to content

Maximum Weight Bipartite Matching

The maximum weight matching in a bipartite graph refers to the matching with the maximum sum of edge weights.

Hungarian Algorithm (Kuhn–Munkres Algorithm)

The Hungarian algorithm, also known as the KM algorithm, can find the maximum weight perfect matching of a bipartite graph in \(O(n^3)\) time.

Considering that the numbers of vertices in the two sets of a bipartite graph are not always equal, to apply the KM algorithm to solve the maximum weight matching problem, the following preprocessing is required: add dummy vertices to the smaller set to make both sides have the same number of vertices, and set the weight of non-existent edges to \(0\). In this case, the problem is transformed into finding a maximum weight perfect matching, which can then be solved using the KM algorithm.

Feasible Labeling

Assign a weight \(l(i)\) to each vertex \(i\), such that for all edges \((u,v)\), \(w(u,v) \leq l(u) + l(v)\).

Equality Subgraph

The spanning subgraph of the original graph under a feasible labeling, containing all vertices but only including edges \((u,v)\) that satisfy \(w(u,v) = l(u) + l(v)\).

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

Proof 1.

Consider any perfect matching \(M\) of the original bipartite graph. Its edge weight sum 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)\)

The edge weight sum of any perfect matching \(M'\) in the equality subgraph of a feasible labeling is

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

Thus, the edge weight sum of any perfect matching will not exceed \(val(M')\), meaning \(M'\) is the maximum weight matching.

With Theorem 1, our goal is to continuously adjust the feasible labeling so that the equality subgraph contains a perfect matching.

Since both sides have the same number of vertices, assume the number is \(n\). Let \(lx(i)\) denote the label of the \(i\)-th vertex on the left, \(ly(i)\) denote the label of the \(i\)-th vertex on the right, and \(w(u,v)\) denote the weight between the \(u\)-th vertex on the left and the \(v\)-th vertex on the right.

First, initialize a feasible labeling, for example:

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

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

Let \(S\) and \(T\) denote the vertices on the left and right sides of the bipartite graph that are in the alternating tree, and \(S'\) and \(T'\) denote those not in the alternating tree.

bigraph-weight-match-1

In the equality subgraph:

  • Edges \(S-T'\) do not exist; otherwise, the alternating tree would grow.
  • Edges \(S'-T\) must be non-matching edges; otherwise, they would belong to \(S\).

Suppose we subtract \(a\) from the labels of vertices in \(S\) and add \(a\) to the labels of vertices in \(T\). We observe:

  • Edges \(S-T\) remain in the equality subgraph.
  • Edges \(S'-T'\) remain unchanged.
  • For edges \(S-T'\), \(lx + ly\) decreases, possibly allowing them to join the equality subgraph.
  • For edges \(S'-T\), \(lx + ly\) increases, so they cannot join the equality subgraph.

Thus, the value \(a\) must be the minimum edge weight among \(S-T'\):

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

When a new edge \((u,v)\) joins the equality subgraph, two cases arise:

  • \(v\) is an unmatched vertex, then an augmenting path is found.
  • \(v\) is already matched to a vertex in \(S'\).

After at most \(n\) label adjustments, an augmenting path can be found.

Each time the labels are adjusted, edges in the alternating tree do not leave the equality subgraph, so we can maintain this tree directly.

For each vertex \(v\) in \(T\), maintain:

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

Thus, the adjustment value \(a\) can be computed in \(O(n)\):

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

When a new vertex joins \(S\) in the alternating tree, \(slack(v)\) needs to be updated in \(O(n)\). Adjusting labels requires subtracting \(a\) from each \(slack(v)\) in \(O(n)\). Once the alternating tree finds an unmatched vertex, an augmenting path is found.

Initially, enumerate \(n\) vertices to find augmenting paths. To find an augmenting path, the alternating tree needs to be extended \(n\) times, each extension requiring \(n\) maintenance operations, resulting in a total of \(O(n^3)\).

Reference Code 
  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;  // matching point for left set
  vector<int> matchy;  // matching point for right set
  vector<int> pre;     // left point connecting to right set
  vector<bool> visx;   // visited array for left
  vector<bool> visy;   // visited array for 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);  // negative values are worse than no match, so set to 0
  }

  bool check(int v) {
    visy[v] = true;
    if (matchy[v] != -1) {
      q.push(matchy[v]);
      visx[matchy[v]] = true;  // in S
      return false;
    }
    // found new unmatched point, update matching points; pre records the connected point on the "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 chance to join equality subgraph, find augmenting path
                                      // return if found, rebuild alternating tree
                return;
              }
            }
          }
        }
      }
      // no augmenting path, adjust 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 labeling
    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

Original paper: The Dynamic Hungarian Algorithm for the Assignment Problem with Changing Costs

Paper with clearer pseudocode: A Fast Dynamic Assignment Algorithm for Solving Resource Allocation Problems

Related OJ problem: DAP

Algorithm Idea
  1. Modify the weights between a single vertex \(u_i\) and all \(v_j\), i.e., a row in the weight matrix.
    • Update label \(lx(u_i) = max(w_{ij} - v_{j}), \forall j\)
    • Delete matches related to \(u_i\)
  2. Modify the weights between all \(u_i\) and a single vertex \(v_j\), i.e., a column in the weight matrix.
    • Update label \(ly(v_j) = max(w_{ij} - u_{i}), \forall i\)
    • Delete matches related to \(v_j\)
  3. Modify the weight between a single vertex \(u_i\) and a single vertex \(v_j\), i.e., a single element in the weight matrix.
    • Perform either operation 1 or 2.
  4. Add a single vertex \(u_i\) or a single vertex \(v_j\), i.e., add or delete a row or column in the weight matrix.
    • Correspondingly perform operation 1 or 2. Note that adding a vertex here only adds the vertex without setting additional weight values; the weight between the new vertex and others is 0.
Algorithm Proof
  • Let the original graph be G, with left and right labels \(\alpha^{i}\) and \(\beta^{j}\), and feasible labeling l. Then \(G_l\) is a subgraph of G, containing vertices and edges from G that satisfy \(w_{ij} = alpha_{i}+beta_{j}\).
  • In the Hungarian algorithm section, Theorem 1 proves: For a feasible labeling, if its equality subgraph has a perfect matching, then that matching is the maximum weight perfect matching of the original bipartite graph.
  • Assume the original optimal matching is \(M^*\). When a modification occurs, we update the feasible labeling according to the rules. Let the updated labels be \(\alpha^{i^*}\) or \(\beta^{j^*}\). The following situations may occur:
    1. An entire row of the weight matrix is modified. Let the modified row be \(i^*\), meaning all edges of \(v_{i^*}\) are modified. Thus, the original label of \(v_{i^*}\) may no longer satisfy the condition because we require \(w_{i^{*}j} \leq alpha_{i^*}+beta_{j}\). However, for other \(u_j\), except for edges related to \(i^*\), their edge weights remain unchanged, so their labels are still valid. Therefore, the algorithm modifies the labels related to \(v_{i^*}\) to make this set of labels a feasible labeling.
    2. An entire column of the weight matrix is modified. Similarly, the algorithm modifies the labels to make this set a feasible labeling.
    3. A single element of the weight matrix is modified. Arbitrarily modify one of the labels to satisfy the labeling condition.
  • Each time the weight matrix is modified, it relates to a specific node, which may be on the left or right. Denote this node as \(x\). This node is matched with some node \(y\) in the original optimal matching. Each modification operation at most unpairs this pair of nodes. Therefore, we only need to run one round of the search in the Hungarian algorithm to obtain a new match, which, according to Theorem 1, is optimal.

The following code should be the code submitted by the second author of paper 2 (the following code is for maximizing weight; the original paper minimizes cost).

Dynamic Hungarian Algorithm Reference Code 
  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
#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;
}

Transformation to Min-Cost Max-Flow Model

Similar to Bipartite Maximum Matching, the maximum weight matching in a bipartite graph can also be transformed into a network flow problem.

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

Connect the source node to each left vertex in the bipartite graph with an edge of capacity \(1\) and cost \(0\). Connect each right vertex in the bipartite graph to the sink node with an edge of capacity \(1\) and cost \(0\).

Next, for each edge in the bipartite graph connecting left vertex \(u\) and right vertex \(v\) with weight \(w\), add an edge from \(u\) to \(v\) with capacity \(1\) and cost \(w\).

Additionally, considering that in maximum weight matching, the number of matched edges may not equal the number of edges in a maximum matching, for each left vertex, also add an edge to the sink node with capacity \(1\) and cost \(0\).

Finding the maximum cost maximum flow of this network yields the answer. At this point, the maximum flow of the network must equal the number of left vertices, and the maximum cost under maximum flow corresponds to a maximum weight matching scheme.

Exercises

UOJ #80. Bipartite Maximum Weight Matching

Template problem

  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
#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);  // Negative values are worse than no match, so setting to 0 has no effect
  }

  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 found, modify 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();
}
本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用