0005: 二分图的最大(权)匹配简述

· · 算法·理论

0005: 二分图的最大(权)匹配简述

二分图的最大(权)匹配问题是图匹配问题的的一个特例, 一般图的算法需要带花树算法, 先咕了.

本文依赖于 0004: 图匹配和增广路的定义与性质简述, 0002: 最大流与最小割问题简述 中的 Dinic 算法以及 0006: 费用流简述.

0. 相关定义

为了方便, 我们给出一些定义: 我们设要求解的二分图为 G=(X,Y,E), 即左部点集合为 X, 右部点集合为 Y, 边集为 E. 此外, 总点数为 n, 总边数为 m, 若一个点 u 为匹配点, 则与它匹配的点为 m_u, 若 u 为非匹配点, 则 m_u=0.

G 是一个带权图, 则边 \langle u,v\rangle 的权值为 w(u,v).

1. 二分图的最大匹配问题

1.1 增广路算法

\text{Augmenting Path Algorithm}

1.1.1 算法思路

此算法本质为第 4 篇(本文依赖) 中的 3.3.

我们考虑从每个非匹配点开始寻找增广路, 由于增广路的长度定为奇数, 所以其两个端点必然分别处于二分图的左右部, 所以我们不妨从左部点开始寻找增广路, 在右部点结束. 因此, 增广路上左部点到右部点的边为非匹配边, 反之为匹配边.

我们枚举左部点, 设当前枚举的点为 x\in X, 如果 x 为匹配点, 直接跳过即可, 因为每个点只要变为匹配点就不会变回非匹配点. 否则, 我们从 x 开始 DFS, 枚举 x 的邻接点 y\in Y, 如果 y 是非匹配点, 则找到增广路, 否则, 我们从 m_y 开始继续 DFS. 于是, 我们可以写出以下的代码.

int mx[maxn],my[maxn],vis[maxn];
//mx 表示左部点的匹配点, my 表示右部点的匹配点, vis 表示 DFS 中是否访问. 
bool dfs(int cur)
{
    vis[cur]=true;
    for(int v1=0;v1<G[cur].size();v1++)
    {
        int v=G[cur][v1];
        if(vis[my[v]])continue;
        if(!my[v]||dfs(my[v]))
        {
            mx[cur]=v;
            my[v]=cur;
            return true;
        }
    }
    return false;
}

但这样会有一个问题, 那就是每对一个左部点开始进行 DFS 时, 就需要将 vis 数组清空, 这个时间复杂度是不可接受的, 其实, 我们只要将每次 DFS 起始点的编号作为 vis 的值即可, 具体代码见 1.1.3.

1.1.2 时间复杂度分析

我们需要枚举每个左部点, 复杂度为 O(n), 每次 DFS 需要 O(m) 的时间, 所以总时间复杂度为 O(nm).

这个复杂度劣于接下来的两种算法, 但它是二分图最大权匹配中 KM 算法的基础.

1.1.3 代码实现

模板题: P3386 【模板】二分图最大匹配.

#include<iostream>
#include<vector>
using namespace std;
const int maxn=5e2+5;
vector<int> G[maxn];
int mx[maxn],my[maxn],vis[maxn];
//mx 表示左部点的匹配点, my 表示右部点的匹配点, vis 表示 DFS 中是否访问. 
int in,im,ie;
bool dfs(int cur,int tag)
{
    vis[cur]=tag;
    for(int v1=0;v1<G[cur].size();v1++)
    {
        int v=G[cur][v1];
        if(vis[my[v]]==tag)continue;
        if(!my[v]||dfs(my[v],tag))
        {
            mx[cur]=v;
            my[v]=cur;
            return true;
        }
    }
    return false;
}
int main()
{
    cin>>in>>im>>ie;
    for(int v1=1;v1<=ie;v1++)
    {
        int iu,iv;
        scanf("%d %d",&iu,&iv);
        G[iu].push_back(iv);
    }
    int ans=0;
    for(int v1=1;v1<=in;v1++)if(!mx[v1])ans+=dfs(v1,v1);
    cout<<ans<<endl;
    return 0;
}

1.2 朴素的基于网络流的算法

1.2.1 算法思路

我们考虑进行建立一个网络 G' 以进行求解. 我们建立超级源汇点 ST, 然后由 SX 中的每个点连一条容量为 1 的边, Y 中的每个点也向 T 连一条容量为 1 的边, 代表任意一个点只能与一条边相匹配. 对于原图 G 中的每条边 \langle u,v\rangle, 建立一条 u\rightarrow v 的容量的边.

随后, 我们对这张图运行 ST 的最大流 (一般采用 Dinic 算法), 其流值即为最大匹配大小, 原图中的边对应在 G' 中的边若满流则说明其属于最大匹配.

1.2.2 时间复杂度分析

见第二篇中的分析, 该情况属于其中最优的性质, 复杂度可以降低到 O(|E|\sqrt{|V|}), 在此处的时间复杂度即为 O(m\sqrt{n}).

1.3 基于网络流的改进算法: HK 算法

\text{Hopcroft-Karp Algorithm}

1.3.1 算法思路

HK 算法与基于网络流的算法本质相同, 但对其进行了常数优化. 我们设分层网络上每个点 u 的层次为 d_u. 特别地, 令 d_T=\min_{v\in Y}\{d_y+1\}, 其中, T 点指虚拟网络中的汇点, 虽然它在 HK 算法中并不存在.

首先, 我们通过 BFS 对网络进行分层. 我们先将所有左部非匹配点加入队列, 并令它们的 d 值为 0. 随后, 当我们取出队头元素 u 时, 枚举它的所有邻接点 v, 若 v 已访问过, 则跳过, 否则令 d_v=d_u+1. 若 v 点为匹配点, 则令 d_{m_v}=d_v+1(因为 m_v 是匹配点, 所以 v 第一次访问说明 m_v 也是第一次访问), 如果 d_v<d_T 则将 m_v 入队, 因为一个 d 值大于 d_T 的点显然无法推流到 T.

接下来, 我们通过 DFS 寻找增广路, 这里, 网络流里的 "增广路" 与图匹配里的 "增广路" 达成了统一, 我们只需要用与增广路算法同样的代码即可, 但我们只走 d_u+1=d_y 的边.

1.3.2 时间复杂度分析

由于 HK 算法本质与网络流算法相同, 所以其复杂度也为 O(m\sqrt n), 但常数有所减小.

1.3.3 代码实现

#include<iostream>
#include<vector>
#include<queue>
using namespace std;
const int maxn=505,inf=0x3f3f3f3f;
int mx[maxn],my[maxn],vis[maxn],dx[maxn],dy[maxn],tag;
//dx[i] 表示左部第 i 个点距 S 点的距离, dy 同理
vector<int> G[maxn];
int in,im,ie;
bool bfs()
{
    queue<int> q;
    for(int v1=1;v1<=in;v1++)
    {
        dx[v1]=inf;
        if(!mx[v1])
        {
            q.push(v1);
            dx[v1]=0;
        }
    }
    for(int v1=1;v1<=im;v1++)dy[v1]=inf;
    int dt=inf;
    while(!q.empty())
    {
        int cur=q.front();
        q.pop();
        for(int v1=0;v1<G[cur].size();v1++)
        {
            int v=G[cur][v1];
            if(dy[v]!=inf)continue;
            dy[v]=dx[cur]+1;
            if(!my[v])dt=dy[v]+1;
            else
            {
                dx[my[v]]=dy[v]+1;
                q.push(my[v]);
            }
        }
    }
    return dt<inf; 
}
bool dfs(int cur)
{
    vis[cur]=tag;
    for(int v1=0;v1<G[cur].size();v1++)
    {
        int v=G[cur][v1];
        if(vis[my[v]]==tag||dx[cur]+1!=dy[v])continue;
        if(!my[v]||dfs(my[v]))
        {
            mx[cur]=v;
            my[v]=cur;
            return true;
        }
    }
    return false;
}
int main()
{
    scanf("%d %d %d",&in,&im,&ie);
    for(int v1=1;v1<=ie;v1++)
    {
        int iu,iv;
        scanf("%d %d",&iu,&iv);
        G[iu].push_back(iv);
    }
    int ans=0;
    while(bfs())
    {
        tag++;
        for(int v1=1;v1<=in;v1++)if(!mx[v1])ans+=dfs(v1);
    }
    cout<<ans<<endl;
    return 0;
}

2. 二分图的最大权匹配问题

2.1 KM算法(匈牙利算法)

\text{Kuhn–Munkres Algorithm}(\text{Hungarian Algorithm})

KM 算法又称为匈牙利算法(NOI大纲 称这两个算法并不相同, 这里依据OI wiki), 可以在 O(n^3) 时间内求出二分图的最大权完美匹配.

2.1.1 应用 KM 算法前的转换

首先, 需要将两个集合中点数比较少的补点, 使得两边点数相同, 此时, 就可以使用邻接矩阵存储该图, 且任意一组完美匹配定然包含所有点.

由于最大权完美匹配与最大权匹配不一定相同, 我们可以按照 第 4 篇中的 2.1 进行转换: 将不存在的边权或者负边权重设为 0, 这种情况下, 问题就转换成求最大权完美匹配问题, 从而能应用 KM 算法求解.

2.1.2 可行顶标&相等子图

我们给每个点 u 赋一个值 l(u), 称为顶标. 若一组顶标是可行的, 则对于任意边 \langle u,v\rangle, 有 l(u)+l(v)\ge w(u,v).

我们称相等子图是一个包含原图 G 的所有点, 但只包含 l(u)+l(v)=w(u,v) 的边 \langle u,v\rangle 的导出子图.

定理 5.1

对于一组可行顶标, 若其相等子图上存在完美匹配, 则该完美匹配就是原图的最大权完美匹配.

证明

设任意一个原图上的完美匹配 M 的边权和为:

val(M)=\sum_{\langle u,v\rangle\in M}w(u,v)\le\sum_{\langle u,v\rangle\in M}l(u)+l(v)=\sum_{u\in V} l(u)

设任意一个相等子图上的完美匹配 M' 的边权和为:

val(M')=\sum_{\langle u,v\rangle\in M'}w(u,v)=\sum_{\langle u,v\rangle\in M'}l(u)+l(v)=\sum_{u\in V} l(u)

因此, 始终有 val(M)\le val(M'), 从而得证.

因此, 我们求解的方法就是不断调整可行顶标, 使得相等子图上存在完美匹配.

由于我们进行了补点, 设 n 表示左右两部各自的点数, 即 |X|=|Y|=n, 设 lx(i) 表示左部第 i 个点的顶标, ly(i) 表示右部第 i 个点的顶标.

2.1.3 寻找增广路&修改顶标

与增广路算法一样, 我们依次枚举左部点, 若当前点 u 为非匹配点, 则查找增广路, 否则跳过.

查找增广路的过程与增广路算法的 DFS 不同, 使用的是 BFS (使用 DFS 会导致时间复杂度达到 O(n^4), 但多数情况下也可通过, 模板不行 ), 但注意我们只能走相等子图上的边.

如果找到增广路就进行增广, 否则就需要修改顶标. 我们设 BFS 访问过的左部点集为 S, 没访问过的为 S'; 访问过的右部点集为 T, 没访问过的为 T'. 则其划分应当类似下图:

我们可以发现:

从而, 我们考虑给 S 中的点的顶标 -a, 给 T 中的点的顶标 +a, 于是会产生以下影响:

因此, 我们每次选择的 a 值应当为 S-T' 中的最小值:

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

我们发现, 每次修改顶标后, 定然有一点 v\in T' 加入新的 T, 而这个 v 要么是非匹配点从而找到增广路, 要么会使与它匹配的 S' 中的点加入 S, 所以至多 n 次修改顶标后就可以找到增广路.

由于 S 中的点不会进入 S', T 中的点不会进入 T' (实际上, S\cup T 构成了一棵 交错树), 所以我们维护每个右部点相邻的所有 S 中的点可以产生的最小的 a, 设为 slack 值:

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

从而我们可以在 O(n) 时间内求出 a 值:

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

2.1.4 时间复杂度分析

我们需要 O(n) 时间枚举左部点, O(n) 次修改顶标, O(n) 时间求出 a 值, 所以总时间复杂度为 O(n^3).

2.1.5 代码实现

模板题: P6577 【模板】二分图最大权完美匹配.

其中, 我们使用 pre 数组记录每个右部点在交错树中的父节点.

#include<iostream>
#include<queue>
#include<cstring>
using namespace std;
const int maxn=505;
typedef long long ll;
const ll inf=0x3f3f3f3f3f3f3f3f;
ll lx[maxn],ly[maxn],G[maxn][maxn],slack[maxn];//使用邻接矩阵存储二分图
bool visx[maxn],visy[maxn];//左右部点在本轮增广过程中是否访问
int matx[maxn],maty[maxn],pre[maxn],in,im;//matx 表示左部点的匹配点, maty 表示右部点的匹配点
queue<int> q;
bool check(int cur)
{
    visy[cur]=true;
    if(maty[cur])
    {
        q.push(maty[cur]);
        visx[maty[cur]]=true;
        return false;
    }
    while(cur)
    {
        maty[cur]=pre[cur];
        swap(cur,matx[pre[cur]]);
    }
    return true;
}
void bfs(int st)
{
    while(!q.empty())q.pop();
    visx[st]=true;
    q.push(st);
    while(true)
    {
        while(!q.empty())
        {
            int cur=q.front();
            q.pop();
            visx[cur]=true;
            for(int v1=1;v1<=in;v1++)if(!visy[v1])
            {
                ll delta=lx[cur]+ly[v1]-G[cur][v1];
                if(slack[v1]>=delta)
                {
                    pre[v1]=cur;
                    if(delta)slack[v1]=delta;
                    else if(check(v1))return;
                }
            }
        }
        ll a=inf;
        for(int v1=1;v1<=in;v1++)if(!visy[v1])a=min(a,slack[v1]);
        for(int v1=1;v1<=in;v1++)
        {
            if(visx[v1])lx[v1]-=a;
            if(visy[v1])ly[v1]+=a;
            else slack[v1]-=a;
        }
        for(int v1=1;v1<=in;v1++)if(!visy[v1]&&!slack[v1]&&check(v1))return;
        //这行要放在上个循环外, 因为如果提前 return 可能会没修改完后边的顶标,我就被这个坑了很多次qwq.
    }
}
int main()
{
    cin>>in>>im;
    for(int v1=1;v1<=in;v1++)
    {
        for(int v2=1;v2<=in;v2++)G[v1][v2]=-inf;
        //由于是最大权完美匹配, 不存在的边要置为 -inf, 只有在最大权匹配中才置为 0.
        lx[v1]=-inf;
    }
    for(int v1=1;v1<=im;v1++)
    {
        int iu,iv;
        ll iw;
        scanf("%d %d %lld",&iu,&iv,&iw);
        G[iu][iv]=max(G[iu][iv],iw);
    }
    for(int v1=1;v1<=in;v1++)for(int v2=1;v2<=in;v2++)lx[v1]=max(lx[v1],G[v1][v2]);
    //初始化顶标
    for(int v1=1;v1<=in;v1++)if(!matx[v1])
    {
        for(int v2=1;v2<=in;v2++)
        {
            visx[v2]=visy[v2]=false;
            slack[v2]=inf;
        }
        bfs(v1);
    }
    ll ans=0;
    for(int v1=1;v1<=in;v1++)ans+=lx[v1]+ly[v1];
    cout<<ans<<endl;
    for(int v1=1;v1<=in;v1++)printf("%d ",maty[v1]);
    return 0;
}

2.2 基于网络费用流的算法

只需要将 1.2 中的基于网络最大流的二分图最大匹配算法中的最大流算法改为费用流算法即可. 与 ST 相邻的边权为 0, 原图上的边\langle u,v\rangle 所对应在网络上的边权为 -w(u,v), 这是因为费用流算法只能求解最小费用最大流问题, 所以我们要取原边权的相反数, 最后再将总费用取反即为答案.

2.2.1 时间复杂度分析

如果使用 SSP 算法, 复杂度即为 O(n^2m), 因为最大流(即匹配数)至多为 n.

3. 习题

见洛谷题单: 二分图最大(权)匹配.

4. 参考资料

  1. 二分图最大匹配 OI wiki
  2. 二分图最大权匹配 OI wiki
  3. 《算法竞赛入门经典, 训练指南》刘汝佳 陈锋
  4. 《算法导论》 第三版
  5. 网络流,二分图与图的匹配 cnblogs