2023.3.28 【图论】KM算法 | 二分图最大权完美匹配
2023.3.28 【模板】KM算法 | 二分图最大权完美匹配
题目概述
给定一张二分图,左右部均有 \(n\) 个点,共有 \(m\) 条带权边,且保证有完美匹配。
求一种完美匹配的方案,使得最终匹配边的边权之和最大。
数据规模与约定
- 对于 \(100\%\) 的数据,满足 \(1\leq n\leq 500\),\(1\leq m\leq n^2\),\(-19980731\leq h \leq 19980731\) 。且保证没有重边。
我们看到,在一般情况下,这道题可以用普通的最大费用最大流解决,但是dinic的上限复杂度是\(O(n^2m)\)的,在这道题中毒瘤出题人就会将其卡成\(O(n^4)\),难以通过这题,所以我们要用一种新的算法——KM算法
km算法(Kuhn-Munkres算法),是一种在二分图上求解最大权完美匹配的算法,用邻接矩阵存图即可。相比费用流较高的复杂度,km算法有着更为优秀的效率,但局限性在于只能做匹配,而不像费用流那样灵活可变。
前置芝士:匈牙利算法
我们记一个点的顶标为\(lx[i]\)(左部点),\(ly[i]\)(右部点),满足性质:对于任意一条边(u,v),\(lx[u] + ly[v] >= e[u][v]\),当取等号时,我们把这条边叫做相等边。
所有点和所有相等边所组成的子图称为相等子图。
核心算法:贪心地将増广所需的边中,边权最大的那些边变成相等边,即逐渐扩大相等子图。
核心性质:扩大相等子图至其刚好有完美匹配时,该匹配即为原图的最大权完美匹配(很好理解,因为扩大相等子图的过程是贪心的)
因为所有点都要被匹配,所以最终的顶标和满足\(\Sigma{lx[u]} + \Sigma{ly[v]} >= \Sigma e_i\)
\(e_i\)即匹配选中的边
所以当\(\Sigma e_i = \Sigma{lx[u]} + \Sigma{ly[v]}\)时,边权和最大,同时我们尽量让顶标和最大,就能得到答案,这是基于贪心思想得出的。
参考Luogu上第一篇题解的定义,设:
\(lx[i]\)左部点的顶标
\(ly[i]\)右部点的顶标
\(vx[i]\)左部点遍历标记
\(vy[i]\)右部点遍历标记
\(px[i]\)左部点匹配
\(px[i]\)右部点匹配
\(pre[i]\)对于一个右部点,当前匹配到的左部点(暂未确定)
\(slack[i]\)对于指向右部点i的所有边,\(min(lx[u]+ly[i]−e[u][i])\)的值,即松弛量(初始化为inf)
为什么\(slack[i]\)如此定义呢?
可以简单地理解为,当无法匹配的时候,我们要适当修改顶标,让其匹配,修改量就应该是所有点可以修改的量中最小的,而为了满足顶标的定义,即\(lx[u] + ly[i] - e[u][i] >= 0\),我们取\(lx[u] + ly[i] - e[u][i]\)中最小的,这样就算是最小的顶标,减去后仍然满足\(lx + ly - e == 0\),满足以上条件。
所以当\(slack[i] == 0\)时,代表i进入了一个相等子图。
为了满足要求,初始的左部点顶标\(lx[i]\)设为与i连接的边中最大的边权,右部点\(ly[i]\)全部设为0。
算法流程
首先,对于每个点,我们都尝试用匈牙利算法来将其匹配一次,为了减小时间复杂度,我们使用bfs搜索(详见Luogu P6577 题解第一篇)。队列里的点即试图匹配的点,我们将当前待匹配点s推进去。
对于当前队头x,我们试图在它上面扩展相等子图,顺便匹配。对于一个右部点i,如果它在当前轮没有被搜索过,即\(vy[i] == 0\)(不代表它没有匹配),我们尝试与其匹配,并且扩展相等子图,即当且仅当\(lx[x] + ly[i] - e[u][i] < slack[i]\)时进行更新(因为我们要尽量使\(slack[i] == 0\),才能添加进相等子图),将\(pre[i]\)记为x,如果\(slack[i] == 0\)的话,就将i与x正式匹配,即记\(vy[i] = 1\),因为队列中待匹配的点是由s引申出来的,所以如果\(py[i] == 0\),即i没有被匹配过,就意味着找到了一条增广路,将增广路上的点的\(px\)和\(py\)
更新。然后退出当前轮匹配,不然就让x与i匹配,让原来匹配的左部点(\(py[i]\))尝试匹配,推进队列。
以上过程更新结束过后,代表没有找到一个增广路,需要修改顶标,记录一个变量
\(d = Min_{vy[i] == 0} {slack[i]}\),我们要将一个暂未匹配的右部点,通过修改顶标,将其变成相等子图中的点。
找好d以后,对于每一条边(u,v):
1.若\(vx[u] == 1 \ \& \ vy[i] == 1\):两个点都已经是相等子图中的点,\(slack\)不变
2.若\(vx[u] == 0\ \& \ vy[i] == 1\):\(slack[i] += d\),因为不知道这条边是不是相等子图 中的边,如果是,则会拉出相等子图,但是\(vy[i] == 1\),下次增广不会用到这条边,如果不是,加后仍然不是,所以没有影响
3.若\(vx[u] == 1\ \& \ vy[i] == 0\),代表这条边不是相等边,但是修改后可能使其成为相等子图中的边,所以\(slack[i] -= d\)
4.若\(vx[u] == 0 \ \& \ vy[i] == 0\),没有遍历到,\(slack[i]\)不变
以上四条可以巧妙地用以下四行解决
for(int i = 1;i <= n;i++)
{
if(vx[i]) lx[i] -= d;
if(vy[i]) ly[i] += d; else slack[i] -= d;
}
更新以后,再扩大一遍相等子图即可。
所有的扩大相等子图操作都在\(v[i] == 0\)的情况下操作,因为\(v[i] == 1\)时代表这个右部点已经加入了相等子图。
因为每扩大一次相等子图就要找一次在相等子图上找一次增广路,所以以上操作是循环的,通过函数实现,找到增广路并更新后直接return即可。
关于更新操作,我们在之前记了一个\(pre[i]\)数组,这样就可以记下当前搜索状态,用bfs替代dfs,更新时沿着pre数组倒回去,每次更新一个右部点(x)和一个左部点(\(pre[x]\))即可。
完结撒花
Code
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll N = 505,inf = 0x3f3f3f3f3f3f3f3f;
ll e[N][N],vx[N],vy[N],px[N],py[N],lx[N],ly[N],slack[N],d,pre[N];
int n,m;
queue <int> q;
inline void upd(int x)
{
int t;
while(x)
{
t = px[pre[x]];
px[pre[x]] = x;
py[x] = pre[x];
x = t;
}
}
inline void bfs(int s)
{
memset(vx,0,sizeof(vx));
memset(vy,0,sizeof(vy));
memset(slack,0x3f,sizeof(slack));
while(!q.empty()) q.pop();
q.push(s);
while(1)
{
while(!q.empty())
{
int x = q.front();
q.pop();
vx[x] = 1;
for(int i = 1;i <= n;i++)//右部点
{
if(!vy[i])
{
if(lx[x] + ly[i] - e[x][i] < slack[i])
{
slack[i] = lx[x] + ly[i] - e[x][i];
pre[i] = x;//pre[i]代表“可能成为增广路的预选值”
if(slack[i] == 0)
{
vy[i] = 1;
if(!py[i])
{
upd(i);
return;
}
else q.push(py[i]);
}
}
}
}
}
d = inf;
for(int i = 1;i <= n;i++)
if(!vy[i])
d = min(d,slack[i]);
for(int i = 1;i <= n;i++)
{
if(vx[i]) lx[i] -= d;
if(vy[i]) ly[i] += d;
else slack[i] -= d;
}
for(int i = 1;i <= n;i++)
if(!vy[i])
if(slack[i] == 0)
{
vy[i] = 1;
if(!py[i])
{
upd(i);
return;
}
else q.push(py[i]);
}
}
}
int main()
{
ll x,y,z;
cin>>n>>m;
for(int i = 1;i <= n;i++)
for(int j = 1;j <= n;j++)
e[i][j] = -inf;
memset(ly,0,sizeof(ly));
for(int i = 1;i <= m;i++)
{
cin>>x>>y>>z;
e[x][y] = z;
lx[x] = max(lx[x],e[x][y]);
}
for(int i = 1;i <= n;i++) bfs(i);//左部点
ll ans = 0;
for(int i = 1;i <= n;i++) ans += e[py[i]][i];
cout<<ans<<endl;
for(int i = 1;i <= n;i++)
cout<<py[i]<<" ";
return 0;
}