AIOJ

Blossom algorithm

The Blossom algorithm, also called the blossom tree algorithm, solves the maximum-cardinality matching problem in general graphs. It was proposed by Jack Edmonds in 1961. With some modifications, it can also solve the maximum-weight matching problem. It was the first algorithm to prove that maximum matching has polynomial complexity.

Unlike bipartite matching, matching in a general graph may involve odd cycles.

general-matching-1

For this graph, directly taking the complement (exchanging matched and unmatched edges) can make the resulting $M$ invalid: some vertices may belong to two matched edges. The problem is caused by the odd cycle.

Now consider an augmenting algorithm for general graphs. From the perspective of bipartite graphs, enumerate an unmatched vertex each time, make it the root, mark it “o”, and then mark vertices alternately “o” and “i”. It is easy to see that an edge from “i” to “o” is a matched edge.

Suppose the current vertex is $v$ and an adjacent vertex is $u$ . There are two cases:

$u$ has not been visited. If $u$ is unmatched, an augmenting path has been found; otherwise, continue searching from the mate of $u$ .
$u$ has already been visited. Encountering a vertex marked “o” means that we need to shrink a blossom; otherwise, we have encountered an even cycle and skip it.

An even cycle can be handled as in a bipartite graph, so it may be ignored. After shrinking the blossom, continue looking for an augmenting path in the new graph.

general-matching-2

Let the original graph be $G$ and the graph after shrinking the blossom be $G'$ . We only need to prove:

If $G$ has an augmenting path, then $G'$ also has one.
If $G'$ has an augmenting path, then $G$ also has one.

general-matching-3

Let the non-tree edge (the edge that forms the cycle) be $(u,v)$ , and define the blossom base as $h=LCA(u,v)$ . An odd cycle is alternating, and the two edges incident to $h$ have the same type—both are unmatched—if and only if this is the case. Therefore, the tree edge entering $h$ must be matched, and every vertex other than $h$ on the cycle has an unmatched edge leading outside the cycle.

There are two ways to leave the cycle through an outside edge: clockwise and counterclockwise.

general-matching-4

Thus, shrinking or not shrinking the blossom does not affect correctness.

In an implementation, after finding a blossom we do not need to actually shrink it. We can use an array to record the blossom rooted at which vertex contains each vertex.

Complexity analysis

For each augmenting-path search, traverse all edges. When a blossom is encountered, maintaining its vertices costs $O(|E|^2)$ .

Enumerating every unmatched vertex as the start of an augmenting-path search gives a total complexity of $O(|V||E|^2)$ .

Reference code

// 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(std::random_device{}());
  vector<int> match(g.n, -1);   // matching
  vector<int> aux(g.n, -1);     // timestamp
  vector<int> label(g.n);       // "o" or "i"
  vector<int> orig(g.n);        // blossom base
  vector<int> parent(g.n, -1);  // parent
  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) {  // found a visited vertex, i.e. the LCA
          return v;
        }
        aux[v] = aux_time;
        if (match[v] == -1) {
          v = -1;
        } else {
          v = orig[parent[match[v]]];  // continue from the parent of the matched vertex
        }
      }
      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) {  // mark the root as "o" and find an augmenting path
        label[u] = 0;
        q.push(u);
      }
      orig[v] = orig[u] = a;  // shrink the blossom
      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);
    // mark the root as "o"; use "0" for "o" and "1" for "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) {  // found an unvisited vertex
          label[u] = 1;        // mark as "i"
          parent[u] = v;
          if (match[u] == -1) {  // found an unmatched vertex
            augment(u);          // find an augmenting path
            return true;
          }
          // found a matched vertex; put its mate in the queue to extend the alternating tree
          label[match[u]] = 0;
          q.push(match[u]);
          continue;
        } else if (label[u] == 0 && orig[v] != orig[u]) {
          // found a visited vertex with the same "o" label: a blossom
          int a = lca(orig[v], orig[u]);
          // find the LCA, then shrink the blossom
          blossom(u, v, a);
          blossom(v, u, a);
        }
      }
    }
    return false;
  };  // bfs

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

    // match vertices whenever possible
    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

  // make a random matching first
  greedy();
  // find augmenting paths from unmatched vertices
  for (int i = 0; i < g.n; i++) {
    if (match[i] == -1) {
      bfs(i);
    }
  }
  return match;
}

UOJ #79. General graph maximum matching

#include <algorithm>
#include <cassert>
#include <iostream>
#include <numeric>
#include <queue>
#include <random>
#include <vector>
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);  // the random seed is irrelevant here
  // chrono::steady_clock::now().time_since_epoch().count() can also be used
  // to obtain the current time
  vector<int> match(g.n, -1);   // matching
  vector<int> aux(g.n, -1);     // timestamp
  vector<int> label(g.n);       // "o" or "i"
  vector<int> orig(g.n);        // blossom base
  vector<int> parent(g.n, -1);  // parent
  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) {  // found a visited vertex, i.e. the LCA
          return v;
        }
        aux[v] = aux_time;
        if (match[v] == -1) {
          v = -1;
        } else {
          v = orig[parent[match[v]]];  // continue from the parent of the matched vertex
        }
      }
      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) {  // mark the root as "o" and find an augmenting path
        label[u] = 0;
        q.push(u);
      }
      orig[v] = orig[u] = a;  // shrink the blossom
      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);
    // mark the root as "o"; use "0" for "o" and "1" for "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) {  // found an unvisited vertex
          label[u] = 1;        // mark as "i"
          parent[u] = v;
          if (match[u] == -1) {  // found an unmatched vertex
            augment(u);          // find an augmenting path
            return true;
          }
          // found a matched vertex; put its mate in the queue to extend the alternating tree
          label[match[u]] = 0;
          q.push(match[u]);
          continue;
        } else if (label[u] == 0 && orig[v] != orig[u]) {
          // found a visited vertex with the same "o" label: a "blossom"
          int a = lca(orig[v], orig[u]);
          // find the LCA, then shrink the blossom
          blossom(u, v, a);
          blossom(v, u, a);
        }
      }
    }
    return false;
  };  // bfs

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

    // match vertices whenever possible
    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

  // make a random matching first
  greedy();
  // find augmenting paths from unmatched vertices
  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;
}

General-graph matching based on Gaussian elimination

Tip

Before reading the following material, you may need to read the sections on matrices in [linear algebra]:

This section introduces a general-graph matching algorithm based on Gaussian elimination. Compared with the traditional blossom algorithm, it is easier to understand and implement, and is convenient for problems such as finding vertices that must occur in a maximum matching. Its disadvantage is a larger constant factor: Gaussian elimination essentially runs at its full $O(n^3)$ complexity, whereas the blossom algorithm usually does not.

Prerequisite: the Tutte matrix

Definition: For an undirected graph $n$ with $G = (V, E)$ vertices, its Tutte matrix $\tilde{A}(G)$ is an $n \times n$ matrix defined by:

\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}

Here $x_{i, j}$ is a variable, so $\tilde{A}(G)$ contains $|E|$ variables.

When there is no ambiguity, we abbreviate $\tilde{A}(G)$ as $\tilde{A}$ below.

Theorem (Tutte's theorem): $G$ has a perfect matching if and only if $\det \tilde{A} \ne 0$ .

Proof

Introduce the notion of an even-cycle cover: an even-cycle cover of an undirected graph $G$ is a collection of even cycles (including 2-cycles) that covers every vertex exactly once.

It is easy to prove that $G$ has a perfect matching if and only if $G$ has an even-cycle cover.

If $G$ has an even-cycle cover, take every other edge on each cycle to obtain a perfect matching.
If $G$ has a perfect matching, turn each matched edge into a 2-cycle to obtain an even-cycle cover.

We now prove that $G$ has an even-cycle cover if and only if $\tilde{A} \ne 0$ .

Consider the definition of the determinant:

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

Here $\pi$ is an arbitrary permutation, and $(-1)^{\pi}$ is $\pi$ when $-1$ has an odd number of inversions and $1$ otherwise.

It is not difficult to see that every permutation can be viewed as an even-cycle cover of $G$ . If this cover contains an odd cycle, the sum of the term and the term obtained by reversing that cycle is $0$ . Therefore, only even-cycle covers can make the determinant different from $0$ . This proves the theorem.

Theorem: $\operatorname{rank}\tilde{A}$ is always even, and the size of a maximum matching in $G$ equals half of $\operatorname{rank}\tilde{A}$ .

Proof

The rank of a skew-symmetric matrix is even; the second statement is left to the reader.

In practice, it is impossible to compute with $|E|$ variables. Instead, choose a field, such as the residue field $p$ of a prime $\mathcal{Z}_p$ , randomly replace the variables with values in $\mathcal{Z}_p$ , and then compute. For convenience, when there is no ambiguity, $\tilde{A}$ below refers directly to the resulting matrix.

Theorem: $\operatorname{rank}\tilde{A}$ is at most twice the size of a maximum matching in $G$ , and the two are equal with probability at least $1 - \frac n p$ .

Since $n$ in general-graph maximum matching is usually no greater than $10^3$ , a prime $p$ on the order of $10^9$ is sufficient in practice.

The theorem shows that if we only need the size of a maximum matching and not the matching itself, a single Gaussian elimination to obtain $\operatorname{rank}\tilde{A}$ is enough, which is much simpler than the blossom algorithm. If we need to output the matching, the problem is slightly more complicated and uses the algorithm below.

Constructing a perfect matching

By Tutte's theorem and the theorem above, if $G$ has a perfect matching, then $\tilde{A}$ is full rank with high probability. For convenience, the following discussion omits “with high probability.”

Let $G$ denote the vertex numbered $i$ in $v_i$ . We have the following theorem:

Theorem: $\tilde{A}^{-1}_{j,i} \ne 0 \iff G - \{v_i, v_j\}$ has a perfect matching.

Inverse and adjugate matrices

For any $n$ -order square matrix $A$ , define its adjugate matrix by $A^*_{i, j} = (-1)^{i + j} M_{j, i}$ , where $M_{j, i}$ is the minor obtained by deleting row $j$ and column $i$ . Equivalently, if $A$ is the matrix of cofactors of $M$ , then $A^* = M^T$ .

Theorem: If $A$ is invertible, then $A^{-1} = \frac 1 {\det A} A^*$ .

Therefore, $A^{-1}_{j, i} \ne 0 \iff M_{i, j} \ne 0$ , which means that the submatrix left after deleting row $A$ and column $i$ of $j$ has full rank.

In other words, if $(v_i, v_j) \in E$ and $\tilde{A}^{-1}_{j, i} \ne 0$ , then there is a perfect matching containing the edge $(v_i, v_j)$ . We call such an edge feasible below.

The theorem gives an obvious brute-force algorithm for finding a perfect matching in an undirected graph $G$ that has one: enumerate $i, j$ each time; if $(v_i, v_j)$ is a feasible edge (the edge exists and $\tilde{A}^{-1}_{j, i} \ne 0$ ), add $(v_i, v_j)$ to the matching, delete both vertices from $G$ , and recompute the new $\tilde{A}^{-1}$ .

There are $\frac n 2$ rounds, each costing $O(n^3)$ , for a total complexity of $O(n ^ 4)$ , which is somewhat slow. When recomputing $\tilde{A}^{-1}$ , however, we need not calculate the inverse from scratch with Gaussian elimination each time; we can use the following theorem:

Theorem (elimination theorem): Let

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}

and suppose $\hat a_{1, 1} \ne 0$ . Then

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

The theorem describes eliminating the first row and first column. It clearly generalizes to eliminating any row and column, so we only need to compute $\tilde{A}^{-1}$ once at the start; deleting two vertices thereafter requires only two $O(n^2)$ elimination operations.

The description is somewhat abstract; see the C++ code

void eliminate(int A[][MAXN], int r, int c) {  // eliminate row r and column c
  row_marked[r] = col_marked[c] = true;        // already eliminated

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

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

There are $\frac n 2$ rounds, each with complexity $O(n^2)$ , so the above algorithm finds a perfect matching in $O(n^3)$ time.

Constructing a maximum matching

We have just solved the problem of constructing a perfect matching, but most problems ask for a maximum matching.

As mentioned above, the size of a maximum matching in $G$ equals half of $\operatorname{rank}\tilde{A}$ . If we can find a maximum-rank nonsingular submatrix of $\tilde{A}$ , then finding a perfect matching in the subgraph induced by that submatrix gives a maximum matching of $G$ .

Equivalently, if $G$ has a perfect matching, then $\tilde{A}$ is full rank; in other words, $\tilde{A}$ has linearly independent rows. If $\tilde{A}$ is not full rank, we can find a basis of $\tilde{A}$ 's rows and retain the corresponding rows and columns to obtain a maximum-rank nonsingular submatrix of $\tilde{A}$ .

After finding the maximum-rank nonsingular submatrix, use the algorithm above to find a perfect matching in the induced subgraph, obtaining a maximum matching of the original graph. Since Gaussian elimination may exchange rows, take care to maintain the vertex numbering during implementation.