zoukankan      html  css  js  c++  java
  • 从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)只能处理完备匹配,万能还是费用流啊。

  • 相关阅读:
    use paramiko to connect remote server and execute command
    protect golang source code
    adjust jedi vim to python2 and python3
    install vim plugin local file offline
    add swap file if you only have 1G RAM
    datatables hyperlink in td
    django rest framework custom json format
    【JAVA基础】网络编程
    【JAVA基础】多线程
    【JAVA基础】String类的概述和使用
  • 原文地址:https://www.cnblogs.com/zhangjianjunab/p/13812944.html
Copyright © 2011-2022 走看看