<script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/3.2.2/es5/tex-mml-chtml.min.js">
</script>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
tex2jax: {
inlineMath: [['$', '$'] ],
displayMath: [ ['$$','$$'], ["\\[","\\]"] ]
}
});
</script>

# Hungarian Algorithm

## 概要

https://judge.yosupo.jp/problem/assignment
$n \times n$ 行列 $A$ が与えられる。 $Opt=\sum A[i][p_i]$ が最小値をとるような順列 $p$ を1つ求める。

https://cp-algorithms.com/graph/hungarian-algorithm.html　の和訳みたいなものです。

## $O(N^4)$ 解法
適当に定数を足して $A[i][j] \geq 0$ としてよい。

**potential** と呼ばれる量を定義する。 $u,v$ を長さ $n$ の列とし、

$$u[i]+v[j] \leq A[i][j]$$

を満たすものとする。また、 $f=\sum u[i]+\sum v[j]$ とする。

#### 補題1　$f \leq Opt$

#### 証明
任意の順列 $p$ に対し、

$$ f=\sum u[i]+v[p[i]] \leq \sum A[i][p[i]] \leq Opt $$

が成り立つ。 $\square$

この $u,v$ を調整することで $f=Opt$ が成り立つようにすることを目指す。

ここで、左右 $n$ 頂点の二部グラフ $G$ を用意し、 $u[i]+v[j]=A[i][j]$ であるとき $(i,j)$ に辺を張る（この条件を **tight** と呼ぶ）。
すると、補題中の不等式の等号成立は $G$ が完全マッチングを持つことと等価になる。そこで、$G$ 内の最大マッチング $M$ も一緒に管理することを考える。

アルゴリズムの概要は以下。

#### Step 1.
$u[i]=v[i]=0,M=\emptyset$ で初期化する。

#### Step 2.
$(G,M)$ に対する増加パス、つまり、 $G$ の辺で $M$ に属するなら右から左、そうでないなら左から右に有向辺を張り、左側の $M$ に含まれない頂点から始まって右側の $M$ に含まれない頂点で終わるパスが存在するかをDFSないしBFSで探す。

もし存在するならそのパスでマッチングを反転すれば $M$ のサイズを増やせる。そうでないなら $M$ が最大マッチングであることが知られている(https://manabitimes.jp/math/1147)。

#### Step 3.
逆に $u,v$ をチューニングして $M$ を増やすことが出来るかを確認する。

$Z_1$ （resp. $Z_2$ ）を、左側(右側)の頂点でStep 2の探索中に到達したものの集合とし、

$$ \Delta = \min_{i \in Z_1,j \notin Z_2} A[i][j]-u[i]-v[j] $$

と置く。

ここで、$\forall i \in Z_1$ について $u[i]$ を $\Delta$ 増やし、$\forall j \in Z_2$ について $v[j]$ を $\Delta$ 減らすという操作を行い、Step 2に戻る。

#### 補題2 $\Delta>0$

#### 証明
$\Delta=0$ とすると、tightな $i \in Z_1,j \notin Z_2$ が存在する。もし $i \to j$ の有向辺が存在すれば探索によって $j \in Z_2$ となるはずなので、有向グラフの作り方から $(i,j) \in M$ でなければならない。しかし、$i \in Z_1$ となるためには辺 $j \to i$ を通る必要があるため $j \in Z_2$ でなければならない。よって矛盾が言えた。 $\square$

#### 補題3
$u,v$ は依然potentialの条件を満たしている。
#### 証明
- $i \in Z_1,j \in Z_2$ のとき、 $i \notin Z_1,j \notin Z_2$ のとき $u[i]+v[j]$ は不変。
- $i \notin Z_1,j \in Z_2$ のとき、 $u[i]+v[j]$ は $\Delta$ 減少する。
- $i \in Z_1,j \notin Z_2$ のとき、 $\Delta$ の定義から不等式は保たれる。　$\square$

#### 補題4
$M$ 内の辺はtightなままである。

#### 証明
補題3の証明から、等号が<に変わるのは $i \notin Z_1,j \in Z_2$ のケースのみだが、もし $(i,j) \in M$ なら辺 $j \to i$ が存在するので $i \notin Z_1$ ではおかしい。$\square$

#### 補題5
$u,v$ の変更によって $|Z_1|+|Z_2|$ のサイズは真に増加する。

#### 証明
変更前に到達できていた頂点は保存されることに注意する。実際、用いる辺は $i \in Z_1,j \in Z_2$ に属し、補題3からtightであることは保存されている。
また、$\Delta$ の作り方から、少なくとも1つの $i \in Z_1,j \notin Z_2$ なるペアがtightになっており、変更後は $j \in Z_2$ となることが分かる。

このサイクルを、Step 2で $M$ が完全マッチングになるまで行う。
補題5から、$u,v$ の更新は高々 $N$ 回行われ、マッチングの更新は高々 $N$ 回であることと、$\Delta$ の計算などで $O(N^2)$ かかることから、$O(N^4)$ のアルゴリズムが得られる。

## $O(N^3)$ 解法

元のアルゴリズムでは有向グラフの探索を全ての頂点から行っていたが、逆に各頂点1つずつ増加パスの探索していくことを考える。

また、到達可能な頂点はincrementalであったことから、途中経過を保持しておけば $u,v$ の更新のあとで探索を再開することができる。

さらに、$\Delta$ の計算についての補助変数を $minv[j]=\min_{i \in Z_1} A[i][j]-u[i]-v[j]$ とすると、 $Z_1$ の要素が追加されたときの更新と $\Delta$ の計算は $O(N)$ でできる。

以上の改善を施すと、各列について
- 増加パスが見つかるまでにpotentialの更新は高々 $O(N)$ 回起こり、更新自体は $O(N)$ で出来る
- $minv$ の管理、有向グラフの探索に $O(N^2)$

よって全体の計算量は $O(N^3)$ となる。

## 実装
- $A$ は 0-indexed　
- 最適解をとる $p$ が返ってくる

template <typename T, T MX> vector<int> Hungarian(int n, vector<vector<T>> &A) {
    vector<int> u(n + 1), v(n + 1), q(n + 1, n), way(n + 1, n);
    for (int i = 0; i < n; ++i) {
        q[n] = i;
        int j0 = n;
        vector<T> minv(n + 1, MX);
        vector<bool> used(n + 1, false);
        do {
            used[j0] = true;
            int i0 = q[j0], j1 = n;
            T delta = MX;
            for (int j = 0; j < n; ++j)
                if (!used[j]) {
                    T cur = A[i0][j] - u[i0] - v[j];
                    if (cur < minv[j])
                        minv[j] = cur, way[j] = j0;
                    if (minv[j] < delta)
                        delta = minv[j], j1 = j;
                }
            for (int j = 0; j <= n; ++j)
                if (used[j])
                    u[q[j]] += delta, v[j] -= delta;
                else
                    minv[j] -= delta;
            j0 = j1;
        } while (q[j0] != n);
        do {
            int j1 = way[j0];
            q[j0] = q[j1];
            j0 = j1;
        } while (j0 != n);
    }
    vector<int> p(n);
    rep(i, 0, n) p[q[i]] = i;
    return p;
}