KM算法小记

这个人踩死了，连 KM 都不会。

之前一直以为费用流一定优于 KM，前几天做一道题，费用流过不了，非要用 KM 才能过，后来才知道 KM 是能做到 \(O(n^3)\) 的。

二分图最大权完美匹配

给定一张二分图，保证有完美匹配。每条边有权值，求权值和最大的完美匹配。

顶标

顶标是 KM 算法的核心工具。

我们给左右两边每个点分配一个数值，称为“顶标”。不妨设 \(lx[i]\) 为左边第 \(i\) 个点的顶标，\(ly[i]\) 则为右边。

这里，我们分配的顶标需要满足 \(\forall (u,v) \in E,\space lx[u]+ly[v]\ge w(u,v)\)，也就是说一条边两端的顶标之和大于等于该边权值。

• 定义 - 相等子图：二分图保留 \(E\) 中满足 \(lx[u]+ly[v]=w(u,v)\) 的边后的图。

我们需要知道一个事实：如果当前的相等子图中存在完美匹配，那么这个匹配一定是原图的最大权完美匹配。

注意到一点，因为 \(lx[u]+ly[v]\ge w(u,v)\)，所以任何一个匹配的权值和 \(\le \sum\limits_u lx[u]+\sum\limits_u ly[u]\)。

而相等子图中的完美匹配正好取到上届 \(\sum\limits_u lx[u]+\sum\limits_u ly[u]\)，不存在匹配的权值之和比他大，所以一定是我们要的答案。

于是问题转化为如何合理分配顶标，使得相等子图存在完美匹配。

分配顶标 —— 增量+调整

我们考虑调整法，初始化顶标，可以把所有顶标赋值为无穷大，或者所连边的权值最大值。

以 \(i=1...n\) 的顺序为每个 \(i\) 找一个匹配。使用匈牙利算法，在相等子图中找一条增广路。

当然可能找不到增广路，此时我们需要调整顶标。设 \(slack[v]\) 表示不在搜索树中的右部点 \(v\) 与有连边的在搜索树中的左部点 \(u\) 的 \(lx[u]+ly[v]-w(u,v)\) 的最小值。

显然一个合法的 \(v\) 的 \(slack[v]\) 一定不为 \(0\)。令 \(d=\min\limits_v \{slack[v]\}\)，我们令所有在搜索树中的左部点 \(u\) 修改 \(lx[u]\gets lx[u]-d\)，令所有在搜索树中的右部点 \(v\) 修改 \(ly[v]\gets ly[v]-d\)。显然，一条在搜索树中的点一定还在相等子图中；一条 \(u\) 在搜索树、\(v\) 不在搜索树的边，修改后有可能出现在相等子图中，容易发现至少会有一条出现。

由于原图保证了一定有完美匹配，所以这样更新下去，一定可以找到合法的增广路。

一次顶标的修改会使最多一个点进入搜索树，每次进行一个增广，每个点都要做，时间为 \(O(n^2m)\)。当 \(m=n^2\) 时，为 \(O(n^4)\)。

点击查看代码

#include<bits/stdc++.h>
#define ll long long
#define pir pair<ll,ll>
#define fi first
#define se second
#define mkp make_pair
#define pb push_back
using namespace std;
const ll maxn=510, mod=998244353;
ll n,m; ll a[maxn][maxn];
vector<pir>to[maxn];
ll mch[maxn],lx[maxn],ly[maxn];
ll visx[maxn], visy[maxn], d[maxn];
ll dfs(ll u){ visx[u]=1;
	for(ll v=1;v<=n;v++)
		if(!visy[v]){
			if(lx[u]+ly[v]==a[u][v]){
				visy[v]=1;
				if(!mch[v]||dfs(mch[v])){
					mch[v]=u; return 1;
				}
			} d[v]=min(d[v],lx[u]+ly[v]-a[u][v]);
		} return 0;
}
void KM(){
	memset(lx,0xcf,sizeof lx);
	for(ll u=1;u<=n;u++)
		for(ll v=1;v<=n;v++) lx[u]=max(lx[u],a[u][v]);
	for(ll i=1;i<=n;i++){
		while(1){
			memset(visx,0,sizeof visx);
			memset(visy,0,sizeof visy);
			memset(d,0x3f,sizeof d);
			if(dfs(i)) break;
			ll dif=1e17;
			for(ll j=1;j<=n;j++)
				if(!visy[j]) dif=min(dif,d[j]);
			for(ll j=1;j<=n;j++)
				if(visx[j]) lx[j]-=dif;
			for(ll j=1;j<=n;j++)
				if(visy[j]) ly[j]+=dif;
		}
	}
	ll res=0;
	for(ll i=1;i<=n;i++) res+=a[mch[i]][i];
	printf("%lld\n",res);
	for(ll i=1;i<=n;i++) printf("%lld ",mch[i]);
}
int main(){
	scanf("%lld%lld",&n,&m);
	for(ll i=1;i<=n;i++)
		for(ll j=1;j<=n;j++) a[i][j]=-1e16;
	for(ll i=1;i<=m;i++){
		ll u,v,w; scanf("%lld%lld%lld",&u,&v,&w);
		a[u][v]=max(a[u][v],w);
	}
	KM();
	return 0;
}

BFS 的优化

考虑每次修改完顶标后，搜索树中的边本质不变，唯一变的是会向外拓展一些新的点。我们每次都做一遍增广，浪费了大量时间。

考虑 DFS 转 BFS，我们每次修改顶标后只需要检查没搜过的所有右部点，其是否可以加入搜索树。

这样一个点只会被搜一次，BFS 部分总时间为 \(O(nm)\)。同时，最多修改 \(n\) 次顶标，每次是 \(O(n)\)，这部分是 \(O(n^3)\)。所以总时间为 \(O(nm+n^3)\)，当 \(m=n^2\) 时，为 \(O(n^3)\)。

点击查看代码

#include<bits/stdc++.h>
#define ll long long
#define pir pair<ll,ll>
#define fi first
#define se second
#define mkp make_pair
#define pb push_back
using namespace std;
const ll maxn=510, mod=998244353;
ll n,m; ll a[maxn][maxn];
vector<pir>to[maxn];
ll mchx[maxn],mchy[maxn], lx[maxn],ly[maxn];
ll visx[maxn],visy[maxn], sla[maxn], pre[maxn]; ll q[maxn], l, r;
void aug(ll u){
	while(u){
		mchy[u]=pre[u];
		swap(u,mchx[pre[u]]);
	}
}
void bfs(ll s){
	memset(visx,0,sizeof visx);
	memset(visy,0,sizeof visy);
	memset(pre,0,sizeof pre);
	memset(sla,0x3f,sizeof sla);
	q[l=r=1]=s;
	while(1){
		while(l<=r){
			ll u=q[l++]; visx[u]=1;
			for(ll v=1;v<=n;v++)
				if(!visy[v]&&sla[v]>lx[u]+ly[v]-a[u][v]){
					sla[v]=lx[u]+ly[v]-a[u][v], pre[v]=u;
					if(!sla[v]){ visy[v]=1;
						if(!mchy[v]){
							aug(v); return;
						} else q[++r]=mchy[v];
					}
				}
		}
		ll dif=1e17;
		for(ll i=1;i<=n;i++)
			if(!visy[i]) dif=min(dif,sla[i]);
		for(ll i=1;i<=n;i++)
			if(visx[i]) lx[i]-=dif;
		for(ll i=1;i<=n;i++)
			if(visy[i]) ly[i]+=dif;
			else sla[i]-=dif;
		for(ll i=1;i<=n;i++)
			if(!visy[i]&&!sla[i]){
				visy[i]=1;
				if(!mchy[i]){
					aug(i); return;
				} q[++r]=mchy[i];
			}
	}
}
void KM(){
	for(ll u=1;u<=n;u++)
		for(ll v=1;v<=n;v++) lx[u]=max(lx[u],a[u][v]);
	for(ll i=1;i<=n;i++){
		bfs(i);
	}
	ll res=0;
	for(ll i=1;i<=n;i++) res+=a[mchy[i]][i];
	printf("%lld\n",res);
	for(ll i=1;i<=n;i++) printf("%lld ",mchy[i]);
}
int main(){
	scanf("%lld%lld",&n,&m);
	for(ll i=1;i<=n;i++)
		for(ll j=1;j<=n;j++) a[i][j]=-1e16;
	for(ll i=1;i<=m;i++){
		ll u,v,w; scanf("%lld%lld%lld",&u,&v,&w);
		a[u][v]=max(a[u][v],w);
	}
	KM();
	return 0;
}