匈牙利算法以及优化

匈牙利算法可以在 $O(nm)$ 的复杂度内解决 $O(n)$ 个点 $O(m)$ 条遍的二分图最大匹配问题。

若左部点 $X$ 和右部点 $Y$ 匹配，我们称 $X, Y$ 都为匹配点，边 $(X, Y)$ 为一条匹配边。此外的点为非匹配点，此外的边为非匹配边。

简述匈牙利算法的基本流程：

• 枚举未匹配左部点：

  1. 从当前左部点开始，搜索出一条以一个未匹配点结尾的交错路，称为增广路。交错路是指，匹配边与非匹配边交替出现的路径。
  2. 如果找到一条增广路，将这条增广路上所有非匹配边变为匹配边，所有匹配边变为非匹配边。这个过程让匹配数加一。

该算法的正确性这篇博客做的很好，可以参考。

操作 1 和 3 都是简单的。如何实现操作 2？通常我们用 dfs 搞定：

const int N = 501;

int n, m, q; // 左部点编号 1 ~ n，右部点编号 1 ~ m，总共有 q 条边

vector<int> e[N]; // 左部点 i 连向 e[i] 的所有点

bool vis[N]; // 本轮 dfs 是否访问过这个点
int pr[N]; // pr[i] 表示右部点 i 与哪个左部点匹配

bool dfs(int x)
{
    for (int i : e[x]) if (!vis[i]) {
        vis[i] = 1;
        // 右部点 i 是非匹配点，那么路径 s 到 i 是一条增广路
        if (!pr[i]) { 
            pr[i] = x;
            return true;
        }
        // (i, pr[i]) 是一条匹配边，把它作为交错路的一部分，然后从 pr[i] 开始继续搜索
        if (dfs(pr[i])) { 
            pr[i] = x;
            return true;
        }
    }
    return false;
}

int match()
{
    int ans = 0;
    for (int i = 1; i <= n; i++) {
        memset(vis, 0, sizeof(vis));
        ans += dfs(i);
    }
    return ans;
}

可以发现每个点最多被访问一遍，那自然每条边最多被访问一遍。故 dfs 搜索增广路的复杂度为 $O(m + n)$。我们会枚举 $n$ 次，每次搜索一遍，所以总复杂度为 $O(nm)$。

bfs 实现

显然我们可以将上文的 dfs 换成 bfs 来搜索一条增广路。但是 bfs 要多开点数组，因为在找到增广路后我们不能回溯，所以要用一个数组存一下当前点在增广路上的前一个点。

const int N = 501;

int n, m, q;
vector<int> e[N];
bool vis[N];
int pre[N]; // 增广路上，右部点 i 的上一个点是左部点 pre[i]
int pl[N]; // 左部点 i 匹配右部点 pl[i]
int pr[N]; // 右部点 i 匹配左部点 pr[i]

// 切换路径匹配状态
void flip(int x)
{
    // x 一定是右部点
    while (x) {
        // x 现在应该和 pre[x] 匹配了
        pr[x] = pre[x];  
        // 将 x 恢复到 pre[x] 原来匹配的点继续切换路径状态。
        int tmp = pl[pre[x]];
        pl[pre[x]] = x;
        x = tmp;
        // 以上 3 行代码可以简写成 swap(x, pl[pre[x]]);
    }
}

bool bfs(int s)
{
    queue<int> Q;
    Q.push(s);
    pre[s] = 0; // 初始化
    while (!Q.empty()) {
        int x = Q.front(); Q.pop();
        for (int i : e[x]) if (!vis[i]) {
            vis[i] = 1;
            pre[i] = x;
            if (!pr[i]) {
                flip(i);
                return true;
            }
            // 入队继续去找增广路，道理和 dfs 一样的。
            Q.push(pr[i]); 
        }
    }
    return false;
}

int match()
{
    int ans = 0;
    for (int i = 1; i <= n; i++) {
        memset(vis, 0, sizeof(vis));
        ans += bfs(i);
    }
    return ans;
}

和 dfs 一样的方法分析复杂度，可得复杂度仍为 $O(nm)$。

在稀疏图上，由于增广路长度通常不大，数量也少，所以 bfs 通常会快一些（dfs 搜索失败代价大）。

bitset 优化

我们并不满足于 $O(nm)$ 的匈牙利算法。实际上有题会给出 $n = 2000, m = O(n^2)$ 规模的二分图匹配问题。此时不管是 dfs 还是 bfs 都会达到 $O(nm) = O(n^3)$ 的复杂度，然后被卡掉。但事实是 $O(nm)$ 已经是匈牙利的极限了，没有什么简单的优化该算法复杂度的方法。所以我们需要优化常数。

注意到，找增广路复杂度 $O(m)$ 大部分贡献都是遍历出边时，访问了已经访问过的点(vis[i] == true)。剩下的过程都是 $O(n)$ 的（增广路长度不超过 $2n$）。要是我们能跳过已经访问过的点就好了。

具体来说，我们用 bitset 存图。若 $g_{i, j} = 1$ 那么左部点 $i$ 和右部点 $j$ 之间有一条边。

我们再开一个 bitset $v$，若 $v_i = 1$ 表示右部点 $i$ 还未被处理过。那么 $g_x \& v$ 就是未被处理过，并且与左部点连通的右部点。使用 _Find_first() 和 _Find_next(i) 遍历即可。设 $v$ 中为 $1$ 的位有 $c$ 个，这样遍历的复杂度为 $O(\frac nw + c)$。那么整个程序的复杂度就是 $O(\frac {n^3}w)$。

代码实现如下：

bitset<N> nw, g[N], vis;

bool bfs(int s)
{
    queue<int> Q;
    Q.push(s);
    pre[s] = 0;
    vis.set();
    while (!Q.empty()) {
        int x = Q.front(); Q.pop();
        nw = vis & g[x];
        for (int i = nw._Find_first(); i <= n; i = nw._Find_next(i)) {
            vis.reset(i);
            pre[i] = x;
            if (!pr[i]) {
                while (i) {
                    pr[i] = pre[i];
                    swap(i, pl[pre[i]]);
                }
                return true;
            }
            Q.push(pr[i]);
        }
    }
    return false;
}

其它部分和原版 bfs 写法一致。这样就可以跑点数较多的稠密图了。