ブログ名

競技プログラミングやお歌のお話をする高菜です。

Hopcroft–Karp の Rust による短くて易しい実装

Hopcroft–Karp のアルゴリズムは、二部グラフ $G = (V (= L \cup R), E)$ の最大マッチングを $O (E \sqrt V)$ で計算します。(ただし、$E = \Omega(V)$ を仮定します。)詳しくは Wikipedia へどうぞです…!

この記事の方法は行数でいうと全部で 40 行以下です。相当詰めたつもりですから、多分これが一番易しいと思いますの気持ちなのですが、さらに易しくできる場合はコッソリ教えていただきたいかもです。

実装だけ欲しい方は一番下までスクロールしていただいてですね。

仕様

いずれにしてもアルゴリズムが最悪 $\omega(V + E)$ 時間ですから、入出力は実行速度を気にせず使いやすいようにしてしまいましょう!

fn hopcroft_karp(edges: &[(usize, usize)]) -> Vec<(usize, usize)>
  • edges: 辺のリスト
  • return value: マッチングの辺のリスト

$|L|, |R|$ は受け取らず、中の実装で何とかすることにします。

また添え字は左右ともに $0$ 始まりとしましょう。その理由は、仮に $L = \lbrace 0, \dots, K - 1 \rbrace$, $R = \lbrace K, \dots, V - 1 \rbrace$ のような形であったとしても、両側 $\lbrace 0, \dots V - 1 \rbrace$ と思って呼び出してもバグりませんし、実行速度もそこまで大きくは変わらないはずだからです。(一応、$E = O ( V )$ ならば定数部分に影響はしますが。)一方、逆に頂点番号が被っていない前提で実装してしまうと、被っているときにバグる可能性がありますから、ここは利用時の安心を買うことにしましょう。

アルゴリズムと実装

アルゴリズムを通して残余グラフ $\tilde G$ に相当する情報を管理ししょう。この記事では

  • 変数 g: もとのグラフ $G$ に $L$ 側から $R$ 側に向きづけたもの(注:実装の都合上、マッチ辺は消しません)
  • 変数 h: マッチングに $R$ 側から $L$ 側に向きづけたもの

を管理します。

Hopcroft–Karp は $O ( \sqrt V )$ 回、各頂点について$L$ の未マッチ頂点からの $\tilde G$ における距離に相当する情報を計算します。 しかし $\tilde G$ において $R$ の各頂点の出次数は $1$ 以下ですから、$L$ か $R$ 片方だけに絞っても そこで、$L$ の各頂点について $L$ の未マッチ頂点からの $\tilde G$ における距離の $1 / 2$ 倍を変数 dist に格納しておきます。これは普通にキュー(queue)を用いた BFS で計算します。

アルゴリズム中で用いる、可変サイズの値を持つ変数は上に述べた g, h, dist, queue の 4 つですべてです。変数名は少ないほうがよいというのが、Rust の考え方ですね。

グラフの構築

まず、グラフ g, h を作ってしまいましょう。h は最初は空グラフですから、これでよいです。

  let left = edges.iter().map(|&(i, _)| i).max().map_or(0, |i| i + 1);
  let right = edges.iter().map(|&(_, j)| j).max().map_or(0, |j| j + 1);
  let mut g = vec![Vec::new(); left];
  edges.iter().for_each(|&(i, j)| g[i].push(j));
  let mut h = vec![None::<usize>; right];

増分路の極大集合の発見 Step 1: 距離計算

大まかに言うと、こういう流れです。

なお while { a; b; c } { } は、C++ でいうと do { a; b; } while (c) に相当します。Rust の場合は見た目ちょっとテクニカルになってしまいますね。

  while {
      let dist = /*距離計算*/;
      /* 増分路の極大集合を DFS で計算 */;
      /* 更新があれば true */
  } {}

距離計算です。みなさま無限回書いたコードなので余裕かと思いますが、コードの分量だけでいうと自体はここが一番重いですね。

  • 始点集合は $L$ の未マッチ頂点全体です。これは変数 h を読むとわかります。
  • 隣接する頂点は、未マッチ辺、マッチ辺の順にたどれるものです。未マッチ辺は g から h を除けばわかります。マッチ辺は h をみるとわかります。
    let mut dist = vec![0_usize; left];
    h.iter().flatten().for_each(|&i| dist[i] = !0);
    let mut queue = (0..left).filter(|&i| dist[i] == 0).collect::<VecDeque<_>>();
    while let Some(i) = queue.pop_front() {
      for i1 in g[i].iter().flat_map(|&j| h[j]) {
        if dist[i1] == !0 {
          dist[i1] = dist[i] + 1;
          queue.push_back(i1);
        }
      }
    }

増分路の極大集合の発見 Step 2: DFS による増分路の列挙

最短路全体のなす $\tilde G$ の部分グラフを DFS しながら、増分路を発見し次第帰りにマッチング h を更新していきます。帰りに処理のある DFS ですから、残念ながらスタックで書いてもあまり簡単にならないため、泣く泣く関数を定義します。ラムダ再帰はまだですか?? Rust さん????

  • 一度訪問した頂点は、増分路に使用されようとされまいと再び訪問する必要がありませんから、DFS で使う訪問済みフラグは共用です。
  • 一度訪問した頂点の dist には興味がないため、そこに $+ \infty$ (!0, std::usize::MAX) を代入することで変数を節約します。
  • map(..).fold(false, bitor) のところ、「要するに .any() でしょう」と思うかもしれないのですが、実はそうではなく true を発見しても最後まで探索を続ける必要があります。さもなくば増分路が高々 1 つしか見つかりません。
  • .any(|&j| ..) のところ、読みづらければみなさまのほうで match などを用いて書き下してみていただいてですね。
        (0..left)
            .map(|i| dist[i] == 0 && path(&g, &mut h, &mut dist, i))
            .fold(false, std::ops::BitOr::bitor)
/* 中略 */

fn path(g: &[Vec<usize>], h: &mut [Option<usize>], dist: &mut [usize], i: usize) -> bool {
    let d = std::mem::replace(&mut dist[i], !0);
    g[i].iter()
        .copied()
        .find(|&j| h[j].map_or(true, |i1| d + 1 == dist[i1] && path(g, h, dist, i1)))
        .into_iter()
        .inspect(|&j| h[j] = Some(i))
        .next()
        .is_some()
}

補遺

アホォ~イ

Ford–Fulkerson 型の遅い版アルゴリズム

Hopcroft–Karp は Dinic 法の特殊な場合なわけですが、同様にして Ford–Fulkerson を特殊化することで二部グラフの最大マッチングのアルゴリズムを作ることができますね。つまり、距離計算をせず増分路を 1 つずつ見つけていく最悪 $\Theta(EV)$ 時間の方法です。

短いのでコードを全部乗せてしまいますが、実質ほぼ空の while 式と path 関数だけになっております。

  • 距離配列を訪問済みフラグの代わりにするのができないため、訪問済みフラグは別で管理しておく必要があることに注意です。
  • 今回は map(..).fold(false, bitor) の代わりに .any() を使ってよいです。
fn slow_bipartite_matching(edges: &[(usize, usize)]) -> Vec<(usize, usize)> {
    let left = edges.iter().map(|&(i, _)| i).max().map_or(0, |i| i + 1);
    let right = edges.iter().map(|&(_, j)| j).max().map_or(0, |j| j + 1);
    let mut g = vec![Vec::new(); left];
    let mut h = vec![None::<usize>; right];
    edges.iter().for_each(|&(i, j)| g[i].push(j));
    while {
        let mut used = vec![false; left];
        (0..left).any(|i| !used[i] && path(&g, &mut h, &mut used, i))
    } {}
    h.iter()
        .enumerate()
        .filter_map(|(j, i)| i.map(|i| (i, j)))
        .collect::<Vec<_>>()
}
fn path(g: &[Vec<usize>], h: &mut [Option<usize>], used: &mut [bool], i: usize) -> bool {
    used[i] = true;
    g[i].iter()
        .copied()
        .find(|&j| h[j].map_or(true, |i1| !used[i1] && path(g, h, used, i1)))
        .into_iter()
        .inspect(|&j| h[j] = Some(i))
        .next()
        .is_some()
}

昔ツイートしたこれはこれのことでした。

さらに余談ですが nightly API の Option::inspect を使うとメソッドチェーンからメソッドが 2 個(.into_iter().next())減ります。

距離配列など使わなくとも BFS 木を計算してそれに沿った増分路を探せば良いのではないでしょうか。

だめです。左側、右側の部集合ともに $0, \dots, N - 1$ であるとして、$i \le j$ のときに辺 $(i, j)$ を張るようなグラフが反例です。

頂点数の小さい順に探索して BFS 木を作っていくとすると、毎回最初の頂点が右側の全頂点を制覇して回ってしまい、増分路が 1 つづつしか見つからず結局 $N$ 回イテレートする羽目になってしまいます。また毎回の BFS では当然全ての辺をチェックする必要がありますから、計算量は $\Theta ( E V )$ になってしまいます。

これは密グラフのときの反例ですが、疎グラフでも反例があるのかどうかはわかりません。実は $O ( E \sqrt E )$ くらいで抑えられていたりしませんか?(泣) なお Library Checker さんにジャッジをお願いしたところ、普通に通ってしまうなど、真実は闇の中へです…… どなたか!(→ 提出リンク

付録: 完全なコード

結局こういう感じになります。Rust 1.42.0 で動きますので、コピペしたい方はどうぞです。もちろん責任は取りません! Libray Checker さんに提出もしてきましたので、心配な方はそちらからぜひです。(→ 提出リンク

use std::collections::VecDeque;

fn hopcroft_karp(edges: &[(usize, usize)]) -> Vec<(usize, usize)> {
    let left = edges.iter().map(|&(i, _)| i).max().map_or(0, |i| i + 1);
    let right = edges.iter().map(|&(_, j)| j).max().map_or(0, |j| j + 1);
    let mut g = vec![Vec::new(); left];
    let mut h = vec![None::<usize>; right];
    edges.iter().for_each(|&(i, j)| g[i].push(j));
    while {
        let mut dist = vec![0_usize; left];
        h.iter().flatten().for_each(|&i| dist[i] = !0);
        let mut queue = (0..left).filter(|&i| dist[i] == 0).collect::<VecDeque<_>>();
        while let Some(i) = queue.pop_front() {
            for i1 in g[i].iter().flat_map(|&j| h[j]) {
                if dist[i1] == !0 {
                    dist[i1] = dist[i] + 1;
                    queue.push_back(i1);
                }
            }
        }
        (0..left)
            .map(|i| dist[i] == 0 && path(&g, &mut h, &mut dist, i))
            .fold(false, std::ops::BitOr::bitor)
    } {}
    h.iter()
        .enumerate()
        .filter_map(|(j, i)| i.map(|i| (i, j)))
        .collect::<Vec<_>>()
}
fn path(g: &[Vec<usize>], h: &mut [Option<usize>], dist: &mut [usize], i: usize) -> bool {
    let d = std::mem::replace(&mut dist[i], !0);
    g[i].iter()
        .copied()
        .find(|&j| h[j].map_or(true, |i1| d + 1 == dist[i1] && path(g, h, dist, i1)))
        .into_iter()
        .inspect(|&j| h[j] = Some(i))
        .next()
        .is_some()
}

ご指摘

[追記]Option を使わない方法

Option のメソッドが使えたほうが短いと思っていましたが、見た感じあまり変わらなさそうですし、こちらでもよさそうです。むしろ map_or とかするくらいなら !0 のほうがマシまであります。

誤字

https://twitter.com/0x3b8001f5/status/1627999134912516097