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)只能处理完备匹配,万能还是费用流啊。

  • 相关阅读:
    Nginx proxy开启cache缓存
    Nginx Server 配置
    nginx 全局配置
    Nginx 配置文件解析
    Python-Scrapy框架
    Python-shutil模块
    os模块3
    os模块
    python os模块atime ,ctime,mtime意义
    Python-正则表达式
  • 原文地址:https://www.cnblogs.com/zhangjianjunab/p/13812944.html
Copyright © 2011-2022 走看看