- 网上讲的对KM算法的优化都是加个slack,看似减掉了 O ( n 2 ) O(n^2) O(n2)的找最小值过程,但仔细想想理论上是没啥用的。因为,众所周知,每次重新建立交替树的过程理论最坏情况就是 O ( m ) O(m) O(m)的。而一般用KM算法的图是完备二部图, O ( m ) = O ( n 2 ) O(m)=O(n^2) O(m)=O(n2)。所以网上的传统优化法实际上复杂度还是 O ( n 4 ) O(n^4) O(n4)的。这篇博文实验了一种更快的,最坏情况真正做到O(n^3)的算法。
- 我用了两个题测了一下,发现果然自己脑补的算法跑的快很多。一道hdu的题,一道ICPC2019南京区赛的题(这道题用网上的假 O ( n 3 ) O(n^3) O(n3)优化算法会超时),坑了不知多少英雄好汉。
二分图最优完美匹配的KM算法极其优化
KM算法是干啥的
KM算法力图求出二分图完美匹配中边权和最大的那个。设二部图分为X点集和Y点集。同点集内任意两点无边相连。
KM算法基本原理和正确性证明
- 首先你得会匈牙利算法,并理解它。考虑到这是一个大家都会的东西这里就不说了,高校很多专业课程都讲过。
- 可行子图引理
- KM算想法理很简单:如果能找到一个顶标的填数方案,使得可行子图有完美匹配,那么这个匹配就是答案。
- 顶标填数方案:每个X里的点i和Y里的点j都填上数,记作 l x [ i ] , l x [ j ] lx[i],lx[j] lx[i],lx[j],满足 l x [ i ] + l x [ j ] ≥ w [ i , j ] lx[i]+lx[j]\ge w[i,j] lx[i]+lx[j]≥w[i,j]。这么一种天数方案就是一个顶标。
- 把所有点和满足 l x [ i ] + l x [ j ] = = w [ i , j ] lx[i]+lx[j]== w[i,j] lx[i]+lx[j]==w[i,j]的边都拎出来,形成一个子图,就是可行子图。
- 若可行子图有完美匹配那这个匹配就是答案。
- 正确性证明:任何一种顶标方案下,所有完美匹配边权和必然小于等于 ∑ l x [ i ] + ∑ l y [ j ] \sum lx[i]+\sum ly[j] ∑lx[i]+∑ly[j]。而根据相等子图的定义,若它有完备匹配,这个匹配的边权和等于顶标和。所以不可能有别的匹配的边权和大于这个匹配。
- 现在就要想办法找到一个有完备匹配顶标填数方案。KM算法过程就是用逐步逼近的思想找这个顶标方案的过程。
- 初始时可以令 l x [ i ] = m a x ( w ( i , j ) ) , l y [ j ] = 0 lx[i]=max(w(i,j)),ly[j]=0 lx[i]=max(w(i,j)),ly[j]=0显然这是满足顶标合法性的。
- 然后按着匈牙利算法那样扫描X集合里的点,每个点都在相等子图里尝试找增广路。
- 然而很可能因为能用的边太少了,导致当前的点i搜不到增广路。所以目前的结果是,以当前的X集中的点i为根(第一层),形成了一棵奇数层数的交替树,奇数层是X集里的点,偶数层是Y里的点。如果不太懂,就请务必再理解理解匈牙利算法。
- 然后满足 l x [ i ] + l y [ j ] = = w [ i , j ] lx[i]+ly[j]== w[i,j] lx[i]+ly[j]==w[i,j]边太少了,导致交替树扩展不动了,就没找到增广路。
- 所以现在要想办法调整顶标,让更多的相等边(满足 l x [ i ] + l x [ j ] = = w [ i , j ] lx[i]+lx[j]== w[i,j] lx[i]+lx[j]==w[i,j]的边)加入到树上,并且保持顶标的合法性。
- 最好的办法是,给树中的所有属于X集的点每个顶标减一个小量delta,同时给在树中的属于Y集的点加一个delta。
- 这样保证了原来树里的边还在树里,不在树里的边现在树里的可能性变大。
- 如何保证顶标合法性呢?只要 d e l t a = m i n ( l x [ i ] + l y [ j ] − w ( i , j ) ) , i ∈ X ∧ i ∈ t r e e , j ∈ Y ∧ j ∉ t r e e delta=min(lx[i]+ly[j]-w(i,j)),i\in X\wedge i\in tree, j\in Y\wedge j\notin tree delta=min(lx[i]+ly[j]−w(i,j)),i∈X∧i∈tree,j∈Y∧j∈/tree
- 这样原树中的边两端点一个+delta一个-delta,和不变,还在树中;X在树中Y不在树中边只有X减了delta,并且delta是这种情况里面能减量的最小值,保证剪完之后还满足顶标合法性。
- 根据上面delta=min…一定有新的边加入树中(就是 l x [ i ] + k x [ j ] − w ( i , j ) = = d e l t a lx[i]+kx[j]-w(i,j)==delta lx[i]+kx[j]−w(i,j)==delta那些边),所以接着搜索增广路搜到的几率就变大了。而顶标的合法性始终满足。
- 正确性:在满足顶标合法性的前提下,每次必然加入新边,若一直没有增广路就会一直加入所有边,这时整个原图都进来了,一定能找到增广路。所以算法是能在有限次调整后收敛到合理的顶标方案的。
KM算法的传统优化代码
这是我发现的网上写的最好懂的板子,结果是正确的。
bool findpath(int x)
{
int tempDelta;
visx[x] = true;
for(int y = 0 ; y < ny ; ++y){
if(visy[y]) continue;
tempDelta = lx[x] + ly[y] - G[x][y];
if(tempDelta == 0){//(x,y)在相等子图中
visy[y] = true;
if(match[y] == -1 || findpath(match[y])){
match[y] = x;
return true;
}
}
else if(slack[y] > tempDelta)
slack[y] = tempDelta;//(x,y)不在相等子图中且y不在交错树中
}
return false;
}
void KM()
{
for(int x = 0 ; x < nx ; ++x){
for(int j = 0 ; j < ny ; ++j) slack[j] = INF;//这里不要忘了,每次换新的x结点都要初始化slack
while(true){
memset(visx,false,sizeof(visx));
memset(visy,false,sizeof(visy));//这两个初始化必须放在这里,因此每次findpath()都要更新
if(findpath(x)) break;
else{
int delta = INF;
for(int j = 0 ; j < ny ; ++j)//因为dfs(x)失败了所以x一定在交错树中,y不在交错树中,第二类边
if(!visy[j] && delta > slack[j])
delta = slack[j];
for(int i = 0 ; i < nx ; ++i)
if(visx[i]) lx[i] -= delta;
for(int j = 0 ; j < ny ; ++j){
if(visy[j])
ly[j] += delta;
else
slack[j] -= delta;
//修改顶标后,要把所有的slack值都减去delta
//这是因为lx[i] 减小了delta
//slack[j] = min(lx[i] + ly[j] -w[i][j]) --j不属于交错树--也需要减少delta,第二类边
}
}
}
}
}
- 我们给每个Y顶点一个“松弛量”函数slack,每次开 始找增广路时初始化为无穷大。在寻找增广路的过程中,检查边(i,j)时,如果它不在相等子图中,则让slack[j]变成原值与 A[i]+B[j]-w[i,j]的较小值。这样,在修改顶标时,取所有不在交错树中的Y顶点的slack值中的最小值作为d值即可。但还要注意一点:修 改顶标后,要把所有不在交错树中的Y顶点的slack值都减去d。
复杂度分析
- 我发现了上面代码一个不爽的问题,就是每次都会把vis清空,重新重头开始搜索建立交替树,这岂不是浪费了很多时间?
- 在原来树的基础上接着只搜索新加入的边岂不美哉?这样每条边至多搜索一次,每个找增广路点的点最坏是 O ( m ) = O ( n 2 ) O(m)=O(n^2) O(m)=O(n2),要搜n个点的增广路,总复杂度就真的是 O ( n 3 ) O(n^3) O(n3)了。
- 否则按上面算法每次重新建树,每次 O ( m ) O(m) O(m),所有点一共最多加m次边,总复杂度最坏是 O ( m 2 ) = O ( n 4 ) O(m^2)=O(n^4) O(m2)=O(n4)。
- 我就是据此脑补出自己的优化的。
不重复构造交替树的更快算法
为了改变上面重复建立交替树的浪费,我改了算法,主要思想就是每次接着原树接着跑增广路,只跑新加入的边。
- 代码:
int findpath(int x)
{
int tempDelta;
visx[x]=true;
for(int y=1;y<=ny;y++){
if(visy[y])continue;
tempDelta =lx[x]+ly[y]-G[x][y];
if(tempDelta == 0){
visy[y] = true;
fa[y+nx]=x;
if(match[y] == -1){
return y+nx;
}
fa[match[y]]=y+nx;//记录交替树的父亲信息(为了区别X,Y集合,Y的点都映射成n+y)
int res=findpath(match[y]);
if(res>0)return res;//返回增广路的末端叶子节点
}
else if(slack[x] > tempDelta)//统计以x为准的slack值。
slack[x] = tempDelta;
}
return -1;
}
void KM()
{
for(int x = 1 ; x <= nx ; ++x){
for(int i = 1 ; i <= nx ; ++i) slack[i] =oo;
for(int i=1;i<=nx+ny;i++)fa[i]=-1;
memset(visx,false,sizeof(visx));
memset(visy,false,sizeof(visy));//换到外面,可以保留原树
int fir=1;int leaf=-1;
while(true){
if(fir==1){
leaf=findpath(x);
fir=0;
}
else{
for(int i=1;i<=nx;i++){
if(slack[i]==0){//只接着搜有新边加入的X点
slack[i]=oo;//slack要重新清空,方以后接着用
leaf=findpath(i);
if(leaf>0)break;
}
}
}
if(leaf>0){
int p=leaf;
while(p>0){
match[p-nx]=fa[p];
p=fa[fa[p]];//顺着记录一路找找上去
}
break;
}
else{
int delta =oo;
for(int i = 1 ; i <= nx ; ++i)
if(visx[i] && delta > slack[i])
delta = slack[i];
for(int i = 1 ; i <= nx ; ++i)
if(visx[i]) {lx[i] -= delta;slack[i]-=delta;}//X点的slack要响应改变,slack变0说明有新边加入
for(int j = 1 ; j <= ny ; ++j){
if(visy[j])
ly[j] += delta;
}
}
}
}
}
- 为了方便,我把slack定义改成slack[i]表示X集中树里的i点到Y集中不在树里的所有j点的 m i n ( l x [ i ] + l y [ j ] − w ( i , j ) ) min(lx[i]+ly[j]-w(i,j)) min(lx[i]+ly[j]−w(i,j))。每次减掉delta后,slack变为0的点下肯定有新加入的边,否则肯定没有。所以不重新搜原来的树,从这些点开始直接接着原树跑,直到交替树扩展到存在增广路为止。
- 所以原来while循环里的vis清空是万万不能要的,要了就是推翻一切重新建树,不行的。要放在外面。
- 这样一条增广路可能是经过很多次交替树的扩大得到的,没法直接统计match的答案。所以我开了个fa数组记录了增广树的信息,搜到增广路的结尾后顺着fa一路向上在填match的答案。
- 注意slack以X集点为准的话和Y为准有不同,每次变0后要重新变为INF,因为这个X集的点可能还要和剩的Y点做delta的最小值选择。而传统算法以Y为准不用这样,因为归零之后这个Y的点就进树了,以后就不会再用到它做delta选择了。
实验比较证据
-
hdu2255最优匹配板题
用网上的板子:
我脑补的板子:
快的挺多的。 -
ICPC2019南京现场赛(我校摘得一块金牌。。。)
-
题目链接:https://nanti.jisuanke.com/t/42404
这个更猛,网上的板子TLE了,我脑补的跑贼快。。。 -
不幸的是现在网上的板子都是 O ( n 4 ) O(n^4) O(n4)或假 O ( n 3 ) O(n^3) O(n3)的,无语了。
-
有人被这道区域赛题坑了,然后评论别人的博客时简单提了自己的想法被我看到了,然后我才脑补了这种方法。