KM算法

原文地址：http://www.cnblogs.com/crazyac/articles/1932808.html

HDU2255

引用1：
KM算法是通过给每个顶点一个标号（叫做顶标）来把求最大权匹配的问题转化为求完备匹配的问题的。设顶点Xi的顶标为A[i]，顶点Yi的顶标为B [i]，顶点Xi与Yj之间的边权为w[i,j]。在算法执行过程中的任一时刻，对于任一条边(i,j)，A[i]+B[j]>=w[i,j]始终 成立。KM算法的正确性基于以下定理： 
　　若由二分图中所有满足A[i]+B[j]=w[i,j]的边(i,j)构成的子图（称做相等子图）有完备匹配，那么这个完备匹配就是二分图的最大权匹配。 
　　这个定理是显然的。因为对于二分图的任意一个匹配，如果它包含于相等子图，那么它的边权和等于所有顶点的顶标和；如果它有的边不包含于相等子图，那么它的边权和小于所有顶点的顶标和。所以相等子图的完备匹配一定是二分图的最大权匹配。 
　　初始时为了使A[i]+B[j]>=w[i,j]恒成立，令A[i]为所有与顶点Xi关联的边的最大权，B[j]=0。如果当前的相等子图没有完备匹配，就按下面的方法修改顶标以使扩大相等子图，直到相等子图具有完备匹配为止。 
　　我们求当前相等子图的完备匹配失败了，是因为对于某个X顶点，我们找不到一条从它出发的交错路。这时我们获得了一棵交错树，它的叶子结点全部是X顶点。现在我们把交错树中X顶点的顶标全都减小某个值d，Y顶点的顶标全都增加同一个值d，那么我们会发现： 
两端都在交错树中的边(i,j)，A[i]+B[j]的值没有变化。也就是说，它原来属于相等子图，现在仍属于相等子图。 
两端都不在交错树中的边(i,j)，A[i]和B[j]都没有变化。也就是说，它原来属于（或不属于）相等子图，现在仍属于（或不属于）相等子图。 
X端不在交错树中，Y端在交错树中的边(i,j)，它的A[i]+B[j]的值有所增大。它原来不属于相等子图，现在仍不属于相等子图。 
X端在交错树中，Y端不在交错树中的边(i,j)，它的A[i]+B[j]的值有所减小。也就说，它原来不属于相等子图，现在可能进入了相等子图，因而使相等子图得到了扩大。 
　　现在的问题就是求d值了。为了使A[i]+B[j]>=w[i,j]始终成立，且至少有一条边进入相等子图，d应该等于min{A[i]+B[j]-w[i,j]|Xi在交错树中，Yi不在交错树中}。 
　　以上就是KM算法的基本思路。但是朴素的实现方法，时间复杂度为O(n4)——需要找O(n)次增广路，每次增广最多需要修改O(n)次顶 标，每次修改顶标时由于要枚举边来求d值，复杂度为O(n2)。实际上KM算法的复杂度是可以做到O(n3)的。我们给每个Y顶点一个“松弛量”函数 slack，每次开始找增广路时初始化为无穷大。在寻找增广路的过程中，检查边(i,j)时，如果它不在相等子图中，则让slack[j]变成原值与A [i]+B[j]-w[i,j]的较小值。这样，在修改顶标时，取所有不在交错树中的Y顶点的slack值中的最小值作为d值即可。但还要注意一点：修改 顶标后，要把所有的slack值都减去d。

 

引用2：

【最优完备匹配】
对于二分图的每条边都有一个权（非负），要求一种完备匹配方案，使得所有匹配边的权和最大，记做最优完备匹配。（特殊的，当所有边的权为1时，就是最大完备匹配问题）
KM算法：（全称是Kuhn-Munkras，是这两个人在1957年提出的，有趣的是，匈牙利算法是在1965年提出的）
为每个点设立一个顶标Li，先不要去管它的意义。
设vi,j��为(i,j)边的权，如果可以求得一个完备匹配，使得每条匹配边vi,j=Li+Lj，其余边vi,j≤Li+Lj。
此时的解就是最优的，因为匹配边的权和=∑Li，其余任意解的权和都不可能比这个大

定理：二分图中所有vi,j=Li+Lj的边构成一个子图G，用匈牙利算法求G中的最大匹配，如果该匹配是完备匹配，则是最优完备匹配。
（不知道怎么证明）
问题是，现在连Li的意义还不清楚。
其实，我们现在要求的就是L的值，使得在该L值下达到最优完备匹配。

L初始化：
Li=max{wi,j}(i∈x,j∈y)
Lj=0

建立子图G，用匈牙利算法求G的最大匹配，如果在某点i (i∈x)找不到增广轨，则得不到完备匹配。
此时需要对L做一些调整：
设S为寻找从i出发的增广轨时访问的x中的点的集合，T为访问的y中的点的集合。
找到一个改进量dx，dx=min{Li+Lj-wi,j}(i∈S,j不∈T)
Li=Li-dx (i∈S)
Li=Li+dx (i∈T)

重复以上过程，不断的调整L，直到求出完备匹配为止。
从调整过程中可以看出：
每次调整后新子图中在包含原子图中所有的边的基础上添加了一些新边。
每次调整后∑Li会减少dx，由于每次dx取最小，所以保证了解的最优性。

复杂度分析：
设n为点数，m为边数，从每个点出发寻找增广轨的复杂度是O(m)，如果找不到增广轨，对L做调整的复杂度也是O(m)，而一次调整或者找到一条增广轨，或者将两个连通分量合成一个，而这两种情况最多都只进行O(n)次，所以总的复杂度是O(nm)

扩展：
根据KM算法的实质，可以求出使得所有匹配边的权和最小的匹配方案。
L初始化：
Li=min{wi,j}(i∈x,j∈y)
Lj=0
dx=min{wi,j-Li-Lj}(i∈S,j不∈T)
Li=Li+dx (i∈S)
Li=Li-dx (i∈T)

#include <iostream>#include <cstdio>#include <algorithm>#include <cstring>#include <cstdlib>#include <queue>#include <stack>#include <set>#include <map>using namespace std;#define REP(i,n)   for(int i=0;i<(n);++i)#define FOR(i,l,r) for(int i=(l);i<=(r);++i)#define DSC(i,r,l) for(int i=(r);i>=(l);--i)#define NX 310#define NY 310#define INF 1e9int n,nx,ny;int lx[NX];int ly[NY],link[NY],slack[NY];bool visitx[NX],visity[NY];int g[NX][NY];bool dfs(int x){    visitx[x]=1;    FOR(y,1,ny)    {        if(visity[y]) continue;        int temp=lx[x]+ly[y]-g[x][y];        if(temp==0)        {            visity[y]=1;            if(link[y]==-1 || dfs(link[y]) )            {                link[y]=x;                return 1;            }        }        else if(slack[y]>temp)            slack[y]=temp;            //不在相等子图中slack 取最小的    }    return 0;}int KM(){    memset(link,-1,sizeof(link));    memset(ly,0,sizeof(ly));    FOR(i,1,nx)    {        lx[i]=-INF;        FOR(j,1,ny) lx[i]=max(lx[i],g[i][j]);        //lx初始化为与它关联边中最大的    }    FOR(x,1,nx)    {        FOR(i,1,ny) slack[i]=INF;        while(1)        {            memset(visitx,0,sizeof(visitx));            memset(visity,0,sizeof(visity));            if(dfs(x)) break;            //若成功（找到了增广轨），则该点增广完成，进入下一个点的增广            //若失败（没有找到增广轨），则需要改变一些点的标号，使得图中可行边的数量增加。            //方法为：将所有在增广轨中（就是在增广过程中遍历到）的X方点的标号全部减去一个常数d，            //所有在增广轨中的Y方点的标号全部加上一个常数temp            int temp=INF;            FOR(i,1,ny)                if(!visity[i])                    temp=min(temp,slack[i]);            FOR(i,1,nx)                if(visitx[i])                    lx[i]-=temp;            FOR(i,1,ny)            if(visity[i]) ly[i]+=temp;            else slack[i]-=temp;            //修改顶标后，要把所有不在交错树中的Y顶点的slack值都减去temp        }    }    int ret=0;    FOR(i,1,ny)        if(link[i]>-1)            ret+=g[ link[i] ][i];    return ret;}int main(){    while(cin>>n)    {        memset(g,0,sizeof(g));        nx=ny=n;        FOR(i,1,n)            FOR(j,1,n)                scanf("%d",&g[i][j]);        cout<<KM()<<endl;    }    return 0;}