从KM算法看二分图带权匹配(雾)

参考文献

BFS的KM写法:https://blog.csdn.net/c20182030/article/details/73330556
正确性证明:算法竞赛进阶指南
DFS和BFS的写法的时间复杂度:https://www.zhihu.com/question/53200316

算法讲解

例题

题目

讲解

转换思路

很明显,每个小人和每个房子的曼哈顿距离就是小人走到房子的代价,于是对于每个小人,向每个房子连边,边权为曼哈顿距离,那么就是求权值最小的完备匹配(每个人都有房子的匹配)。

KM算法讲解

事实上KM求的是最大权匹配

取个负不就行了

下文的证明默认是求最大权匹配,且边权都是整数(当然,后面会讲小数做法)。

事实上,如果没有的边会默认边权为(-∞)

交错树

以下摘自算法进阶(因为这种专业术语我是真的不会啊)。
在匈牙利算法的(DFS)中,如果从某个左部借点出发寻找匹配失败,那么所有访问过的节点,以及为了访问这些节点而经过的边,共同构成的一棵树。
这棵树的根是一个左部节点,所有叶子节点也都是左部节点(因为最终匹配失败了),并且树上第(1,3,5,...)层都是非匹配边,第(2,4,6...)层的边都是匹配边。因此,这棵树被称为交错树。

顶标

全称“顶点标记点”。

在图中,我们会给每个点一个权值,但是在初始化的时候,左边的点和右边的点是不一样的(左边的点集为(|X|),右边的点集为(|Y|)),右边的点权一开始为(0),而右边的点(x_i)的初始化为(max(w(x_{i},y_{j}))),其中(w)表示边权,(y_{j}∈Y)

相等子图

二分图的所有节点,以及满足(d_{x_{i}}+d_{y_{j}}=w(x_{i},y_{j}))的边所构成的子图叫相等子图。

过程

首先,在相等图上跑匈牙利算法。

假如目前是(dfs(x_{i})),如果没有找到增广路,说明相等子图不够大,那么对于访问过的(x_{i})(y_{j})(d_{x_{i}}--,d_{y_{j}}++)

然后再跑匈牙利,反复如此,知道(x_{i})找到增广路。

那么,为什么这样子干,最终一定会找到增广路呢?

如果(x_{i},y_{j})都被访问过了,那么(x_{i}+y_{j})并不会变,所以((x_{i},y_{j}))这条边是否在子图中并没有变化。

但是如果(x_{i})访问过,(y_{j})没有被访问过,那么(x_{i}+y_{j})会比原来小(1),有可能被加入到子图中,这样每次顶标的修改,就意味着相等子图可能在扩大,实际上,如果这条边加入进去,则意味着能多访问到一个点了。

但是如果(y_{j})被访问过了,而(x_{i})没有被访问,那么(x_{i}+y_{j})会比原来大(1),直接被删除出去,但是这条边只有在(x_{i})要访问(y_{j})才会用到,但是(y_{j})我原本就能访问到,还需要这条边干嘛?所以不会导致能访问到的点减少。

所以修改顶标只会导致能访问到的点增加,所以最后一定能匹配上。

最后全部匹配完了,答案就是顶标之和。

事实上,改变顶标还有一个性质,因为左边的点除了(DFS)进入的点,都是通过匹配边到达的,而右边的点被访问之后,因为不存在增广路,则一定会有匹配边,便一定会访问一个没有被访问过的左部点,所以左部点的个数是右部点个数(+1),所以这样修改会使顶标之和(-1)

正确性证明

首先,完备匹配的权值和一定小于等于顶标之和。

我们先证明(x_{i}+y_{j}≥w(x_{i},y_{j})),刚开始的初始化肯定满足,且如果这条边被加入了相等子图之后,(x_{i})被访问,(y_{j})一定也会被访问到,所以(x_{i}+y_{j})不变,而(y_{j})单独被访问,则(x_{i}+y_{j})只会增加,也不会小于(w(x_{i},y_{j}))(还会被踢出相等子图),如果不在相等子图,减一最多从('>')变成('='),然后加入相等子图,如果还是('>'),那也还是大于等于(w(x_{i},y_{j})),得证。(简单来说就是上述过程一定一直满足(x_{i}+y_{j}≥w(x_{i},y_{j}))

根据上述结论,每种完备匹配的匹配边权值和小于等于顶标之和。

而在找到完备匹配方案时,这种方案恰好等于顶标之和,又根据方案都小于等于顶标之和,所以这种方案是最大权的,等于现在的顶标之和。

事实上用此方法应该可能可以证明用此方法也可以做(n)个人,(m)个房子,每个人都要有个房子的最大权匹配。也许

优化

我们发现,这种方法是在是太慢了,尤其是当边权很大的时候,时间简直太慢了,我们称其为(1)号做法,而且这种做法还很难处理小数(难道小数用(1e-4)什么的?)。

我们发现,有时候(+1,-1)还不一定有边会加入子图当中,因为(min(d_{x_{i}}+d_{y_{j}}-w(x_{i},y_{j})))可能(>1),因此可以直接用(slack)存储不在相等子图的边的(d_{x_{i}}+d_{y_{j}}-w(x_{i},y_{j})),顶标一次(+slack)(-slack),那么每次加减都可以让一条边加入到相等子图中。

且这样子做,小数边权的问题也就迎刃而解了,因为(slack)也可以是小数吗。

我们称其为(2)号做法。

以下代码的边权就是小数的。

double  vwh[N],vbl[N],slack;//边权
int  match[N];
bool  st1[N],st2[N];
bool  check(int  x)
{
	st1[x]=1;
	for(int  k=last[x];k;k=a[k].next)//实际上n^2直接用邻接矩阵即可
	{
		int  y=a[k].y;
		if(st2[y])continue;
		if(vwh[x]+vbl[y]-a[k].c>1e-4)slack=mymin(slack,(vwh[x]+vbl[y])-a[k].c);
		else
		{
			if(!match[y]){match[y]=x;return  1;}
			else  if(!st1[match[y]])
			{
				st2[y]=1;
				if(check(match[y])==1){match[y]=x;return  1;}
			}
		}
	}
	return  0;
}

int  main()
	for(int  i=1;i<=n;i++)
	{
		while(1)
		{
			slack=999999999.0;
			memset(st1,0,sizeof(st1));memset(st2,0,sizeof(st2));
			if(check(i)==1)break;
			for(int  j=1;j<=n;j++)
			{
				if(st1[j])vwh[j]-=slack;
				if(st2[j])vbl[j]+=slack;
			}
		}
	}

但是,这种做法并不是最快的。

因为我们发现,(slack)如果在(DFS)中判断,加入的边的两个顶点可能都是被访问的点,所以这条边加入进去也不会增加能被访问点的个数(换而言之这条边很废)

因此,考虑(DFS)后再外面找(slack),只找左边被访问过,右边没被访问过的边。

我们称其为(3)号做法。

double  ma[N][N];//邻接矩阵
double  vwh[N],vbl[N],slack;//让白点找黑点
int  match[N],n;
bool  st1[N],st2[N];
bool  check(int  x)
{
	st1[x]=1;
	for(int  y=1;y<=n;y++)
	{
		if(st2[y])continue;
		if(vwh[x]+vbl[y]-ma[x][y]<=1e-4)
		{
			if(!match[y]){match[y]=x;return  1;}
			else  if(!st1[match[y]])
			{
				st2[y]=1;
				if(check(match[y])==1){match[y]=x;return  1;}
			}
		}
	}
	return  0;
}

int  main()
	for(int  i=1;i<=n;i++)
	{
		while(1)
		{
			slack=999999999.0;
			memset(st1,0,sizeof(st1));memset(st2,0,sizeof(st2));
			if(check(i)==1)break;
			for(int  j=1;j<=n;j++)
			{
				if(st1[j])
				{
					for(int  k=1;k<=n;k++)if(!st2[k])slack=mymin((vwh[j]+vbl[k])-ma[j][k],slack);
				}
			}
			for(int  j=1;j<=n;j++)
			{
				if(st1[j])vwh[j]-=slack;
				if(st2[j])vbl[j]+=slack;
			}
		}
	}

但是,这样还不是最优秀的,因为我们发现,每一次跑完(DFS)并修改顶标,能访问到的点集不会缩小,只会扩大,且原本就在点集中的元素不会消失,所以并不用每次(DFS)都去跑一遍,可以用(BFS)代替修改顶标扩展子图的过程。

此代码中的(slack)是个数组,表示如果下一次要把这个点加入进去,(slack)要等于多少(在下文代码中,真正的(slack)(d)变量)

我们称其为(4)号做法。

int  pre[N],match[N];
double  slack[N],val[N][N],d1[N],d2[N];
bool  v[N];//判断右边的点是否被访问过
void  bfs(int  st)
{
	double  d;
	for(int  i=1;i<=n;i++)slack[i]=999999999.0,v[i]=0;//初始化
	match[0]=st;//为了方便后面的while设置的
	int  py=0,px,yy;
	while(match[py]!=0)//如果现在找到的右边的点没有配偶,那么找到增广路
	{
		px=match[py];v[py]=1;
		d=999999999.0;
		for(int  i=1;i<=n;i++)
		{
			if(!v[i])
			{
				if(slack[i]>d1[px]+d2[i]-val[px][i])slack[i]=d1[px]+d2[i]-val[px][i],pre[i]=py;//更新这个点的slack
				if(slack[i]<d)d=slack[i],yy=i;//判断下一个要加入的点是谁
			}
		}
		for(int  i=0/*注意这个地方为0*/;i<=n;i++)
		{
			if(v[i])d1[match[i]]-=d,d2[i]+=d;
			else  slack[i]-=d;//边权被改动了,所以slack[i]-=d
		}
		py=yy;
	}
	while(py)match[py]=match[pre[py]],py=pre[py];//匹配
}
void  KM()
{
	for(int  i=1;i<=n;i++)bfs(i);//对每个左边的点都这么干
}

但实际上,在最后一道练习题我才用了这个实现方式,其他的都是用3号做法的。

接下来分析没种做法的时间复杂度问题。

(1)号做法,理论上等于(找一个点的对象时扩展(n)个点的顶标变化量(×n^3)),但是由于顶标改动会给后面带来影响(后面找对象更少的改动顶标),而且这种本质上也不是定值,以及一些更加神奇的事情,分析难度很大(因为我原本就不是很会分析时间复杂度,只是乱搞罢了),但是其实还有一种比这个更加准确一点应该的时间复杂度,为(O((开始顶标和-答案)n^2))。(这个我是真的不会。。。)

(2)号做法,每次修改顶标一定会导致一条边加入到相等子图中,所以是(O(mn^2))(n^2)是匈牙利的时间),但是发现可能在这次(DFS)中,会把一些边从相等子图删除,但是每条边对于一个点只会删除一次,加入一次,所以应该是(O(mn^3))的,也就是(O(n^5))的。但是他们都说是n^4的,果然是我错了吗

(3)号做法,由于每次修改顶标一定会有一个点加入进来,所以为(O(n*n*n^2)=O(n^4))

(4)号做法,每次修改顶标只会有一个点加入进来,加入进来只会跑其的(n)条边,所以是(O(n*n*n)=O(n^3))的。(实际上就是把(3)号做法(n)次修改顶标重复跑的部分合成了一次跑),正好对应了(KM)的时间复杂度是(O(n^3))的。

所以,(KM)算法本质是(O(n^3))的,和匈牙利一个样子,真的是优美的算法啊。

代码真的是很神奇的东西,稍微加上一点优化,就少了一个(n)当然也有可能是我分析错了。

代码

本代码采用(3)号做法。

#include<cstdio>
#include<cstring>
#include<cmath>
#define  N  110
#define  M  11000
using  namespace  std;
template<class  T>
inline  T  mymin(T  x,T  y){return  x<y?x:y;}
template<class  T>
inline  T  mymax(T  x,T  y){return  x>y?x:y;}
template<class  T>
inline  T  zabs(T  x){return  x<0?-x:x;}
struct  ZUOBIAO
{
	int  x,y;
}a1[N],a2[N];
inline  double  dis(ZUOBIAO  x,ZUOBIAO  y){return  zabs(x.x-y.x)+zabs(x.y-y.y);}
int  n,m,n1,n2;
char  st[N];
int  ma[N][N];

bool  st1[N],st2[N];
int  match[N];
int  slack,d1[N],d2[N];
bool  dfs(int  x)
{
	st1[x]=1;
	for(int  i=1;i<=n1;i++)
	{
		if(!st2[i]  &&  d1[x]+d2[i]==ma[x][i])
		{
			st2[i]=1;
			if(!match[i]  ||  dfs(match[i])==1){match[i]=x;return  1;}
		}
	}
	return  0;
}
int  KM()
{
	for(int  i=1;i<=n1;i++)
	{
		d1[i]=-999999999;d2[i]=0;match[i]=0;
		for(int  j=1;j<=n1;j++)
		{
			ma[i][j]=-dis(a1[i],a2[j]);//取负号,找最大
			d1[i]=mymax(d1[i],ma[i][j]);
		}
	}
	for(int  t=1;t<=n1;t++)
	{
		while(1)
		{
			slack=999999999;
			memset(st1,0,sizeof(st1));
			memset(st2,0,sizeof(st2));
			if(dfs(t)==1)break;
			for(int  i=1;i<=n1;i++)
			{
				if(!st1[i])continue;
				for(int  j=1;j<=n1;j++)
				{
					if(!st2[j])slack=mymin(slack,d1[i]+d2[j]-ma[i][j]);
				}
			}
			for(int  i=1;i<=n1;i++)if(st1[i])d1[i]-=slack;
			for(int  i=1;i<=n1;i++)if(st2[i])d2[i]+=slack;
		}
	}
	int  ans=0;
	for(int  i=1;i<=n1;i++)ans+=d1[i]+d2[i];//顶标之和就是答案
	return  -ans;
} 
int  main()
{
	while(1)
	{
		n1=n2=0;
		scanf("%d%d",&n,&m);
		if(!n)break;
		for(int  i=1;i<=n;i++)
		{
			scanf("%s",st+1);
			for(int  j=1;j<=m;j++)
			{
				if(st[j]=='H')a2[++n2].x=i,a2[n2].y=j;
				else  if(st[j]=='m')a1[++n1].x=i,a1[n1].y=j;
			}
		}
		printf("%d
",KM());
	}
	return  0;
}

练习

蚂蚁

题目

这道题目是很巧妙的证明加上套模板。

考虑黑边点匹配的边权为他们的距离,权值和最小的匹配方案就是没有交点的方案。

在这里插入图片描述
如果有交点,也就是上图直线方案,换成虚线方案距离和一定会变小,根据三角形两边之和大于第三边就可以证明了。

由于方案的距离和不可能无限变小,因为方案总数也就(n!),最终会变成没有交点的方案,证毕。

然后建模跑模板即可,本题仍然使用(3)号做法。

#include<cstdio>
#include<cstring>
#include<cmath>
#define  N  110
#define  M  11000
using  namespace  std;
template<class  T>
inline  T  mymin(T  x,T  y){return  x<y?x:y;}
template<class  T>
inline  T  mymax(T  x,T  y){return  x>y?x:y;}
struct  DIAN
{
	int  x,y;
}bl[N],wh[N];
inline  double  dis(DIAN  x,DIAN  y){return  sqrt((x.x-y.x)*(x.x-y.x)+(x.y-y.y)*(x.y-y.y));}
double  ma[N][N];
double  vwh[N],vbl[N],slack;//让白点找黑点
int  match[N],n;
bool  st1[N],st2[N];
bool  check(int  x)
{
	st1[x]=1;
	for(int  y=1;y<=n;y++)
	{
		if(st2[y])continue;
		if(vwh[x]+vbl[y]-ma[x][y]<=1e-4)
		{
			if(!match[y]){match[y]=x;return  1;}
			else  if(!st1[match[y]])
			{
				st2[y]=1;
				if(check(match[y])==1){match[y]=x;return  1;}
			}
		}
	}
	return  0;
}
int  main()
{
	scanf("%d",&n);
	for(int  i=1;i<=n;i++)scanf("%d%d",&bl[i].x,&bl[i].y);
	for(int  i=1;i<=n;i++)scanf("%d%d",&wh[i].x,&wh[i].y);
	for(int  i=1;i<=n;i++)
	{
		vwh[i]=-999999999.0;
		for(int  j=1;j<=n;j++)
		{
			ma[i][j]=-dis(wh[i],bl[j]);
			vwh[i]=mymax(vwh[i],ma[i][j]);
		}
	}
	for(int  i=1;i<=n;i++)
	{
		while(1)
		{
			slack=999999999.0;
			memset(st1,0,sizeof(st1));memset(st2,0,sizeof(st2));
			if(check(i)==1)break;
			for(int  j=1;j<=n;j++)
			{
				if(st1[j])
				{
					for(int  k=1;k<=n;k++)if(!st2[k])slack=mymin((vwh[j]+vbl[k])-ma[j][k],slack);
				}
			}
			for(int  j=1;j<=n;j++)
			{
				if(st1[j])vwh[j]-=slack;
				if(st2[j])vbl[j]+=slack;
			}
		}
	}
	for(int  i=1;i<=n;i++)printf("%d
",match[i]);
	return  0;
}

[SDOI2018]新生舞会

题目

很明显的(01)分数规划+KM算法。

这道题目放在最后只是单纯的因为我打的是(n^3)打法。只是单纯的因为我是第二快

#include<cstdio>
#include<cstring>
#define  N  110
using  namespace  std;
//都是n^2的了,用什么边目录啊,正经人谁用边目录啊。
template<class  T>
inline  T  mymax(T  x,T  y){return  x>y?x:y;}
int  n;
int  pre[N],match[N];
double  slack[N],val[N][N],d1[N],d2[N];
bool  v[N];
void  bfs(int  st)
{
	double  d;
	for(int  i=1;i<=n;i++)slack[i]=999999999.0,v[i]=0;
	match[0]=st;
	int  py=0,px,yy;
	while(match[py]!=0)
	{
		px=match[py];v[py]=1;
		d=999999999.0;
		for(int  i=1;i<=n;i++)
		{
			if(!v[i])
			{
				if(slack[i]>d1[px]+d2[i]-val[px][i])slack[i]=d1[px]+d2[i]-val[px][i],pre[i]=py;
				if(slack[i]<d)d=slack[i],yy=i;
			}
		}
		for(int  i=0;i<=n;i++)
		{
			if(v[i])d1[match[i]]-=d,d2[i]+=d;
			else  slack[i]-=d;
		}
		py=yy;
	}
	while(py)match[py]=match[pre[py]],py=pre[py];
}
int  a[N][N],b[N][N];
bool  pd(double  shit)
{
	for(int  i=1;i<=n;i++)
	{
		d2[i]=0;d1[i]=-999999999.0;match[i]=0;
		for(int  j=1;j<=n;j++)val[i][j]=a[i][j]-shit*b[i][j],d1[i]=mymax(val[i][j],d1[i]);
	}
	for(int  i=1;i<=n;i++)bfs(i);
	double  ans=0;
	for(int  i=1;i<=n;i++)ans+=d1[i]+d2[i];
	return  ans>=0;
}
int  main()
{
	scanf("%d",&n);
	for(int  i=1;i<=n;i++)
	{
		for(int  j=1;j<=n;j++)scanf("%d",&a[i][j]);
	}
	for(int  i=1;i<=n;i++)
	{
		for(int  j=1;j<=n;j++)scanf("%d",&b[i][j]);
	}
	double  l=0,r=10000,mid,ans=0;
	while(r-l>=1e-9)
	{
		mid=(l+r)/2;
		if(pd(mid)==1)ans=mid,l=mid;
		else  r=mid;
	}
	printf("%.6lf
",ans);
	return  0;
}

竟然还比费用流跑得快。

事实上我也不知道费用流时间复杂度是多少。

最快的貌似也是BFS KM

小结

但是实际上(KM)只能处理完备匹配,万能还是费用流啊。

【推广】 免费学中医,健康全家人
原文地址:https://www.cnblogs.com/zhangjianjunab/p/13812944.html