[学习笔记]二分图最大权完美匹配-KM算法
壹、模板测试链接
贰、说明与概念
我们现在要解决的问题：在一个二分图中 \(G=\lang V_1,V_2,E\rang\) 中，\(V_1,V_2\) 是两个点集且 \(|V_1|=|V_2|\)，现在所有的 \(e_i\in E\) 都有一个边权 \(w_i\in R\)，求在 \(G\) 中的一种完美匹配方案，使得所有匹配边的权和最大，记做最大权完美匹配。
特殊地，如果 \(|V_1|\neq |V_2|\)，那么我们需要将两个集合中点数比较少的补点，使得两边点数相同，再将不存在的边权重设为 \(0\).
这个问题对于网络流来说是一种费用流情况，用其思想时间复杂度不是十分优秀，我们考虑使用一种新的算法 —— \(\tt KM\) 算法。
\(\tt KM\) 算法使用的一些概念需要进行定义。
2.1.可行顶标
对于原图中的任意一个结点，给定一个函数 \(L(node)\) 求出结点的顶标值。我们用数组 \(L_1\) 记录集合 \(V_1\) 中的结点顶标值，用数组 \(L_2\) 记录集合 \(V_2\) 中的结点顶标值。顶标值并非没有任何要求，在图 \(G\) 中，若有 \(e=\lang u,v,w\rang\in E\)，那么一定有
2.2.相等子图
相等子图的定义是基于可行顶标的，相等子图是 \(G\) 的一个子图。
定义可行边满足 \(e=\lang u,v,w\rang\in E\) 且 \(L_1(u)+L_2(v)=w\). 那么相等子图即为所有可行边构成的图.
叁、一些定理
说明算法流程前，首先引出一些定理。
3.1.定理一
定理内容：
如果原图的一个相等子图中包含完备匹配（完美匹配针对原图而言），那么这个匹配就是原图的最大权完美匹配。
证明：
由于算法中一直保持顶标的可行性，所以任意一个匹配的权值之和肯定小于等于所有结点的顶标之和，则相等子图中的完美匹配（完美匹配针对原图而言）肯定是最大权完美匹配。
3.2.定理二
对于原图 \(G\) 的任何一个相等子图 \(G'=\lang V_1',V_2',E'\rang\)，一定满足
并且，如果 \(G'\) 中存在一种完美匹配 \(M\)（完美匹配针对原图而言）有
那么匹配 \(M\) 就是 \(G\) 的最大权完美匹配。
肆、算法流程
根据定理一与定理二，\(\tt EK\) 算法其实就是不断调整 \(L_1,L_2\) 使得相等子图有完备匹配。我们的目的就是如何实现调整的过程。
现在假设点数为 \(n\)，\(w(i,j)\) 表示由 \(i\in V_1\) 连向 \(j\in V_2\) 的边的边权，并且如果 \(\lang i,j\rang \notin E\)，那么 \(w(i,j)=-\infty\)，首先进行初始化：
不难发现，初始化之后每条边 \(e=\lang u,v,w\rang\in E\) 都满足 \(L_1(u)+L_2(v)\ge w\).
初始化之后，考虑对 \(V_1\) 中尚未匹配的点在相等子图中进行增广，如果找到了增广路就增广，如果找不到，我们就会得到一个交错树。
交错树 ：增广路径形成的树。
假设 \(S,T\) 表示 \(V_1,V_2\) 中在交错树中的点的点集，\(S',T'\) 表示 \(V_1,V_2\) 中不在在交错树中的点的点集，我们现在取一个 \(d\)，并对于所有 \(i\in S\)，让 \(L_1(i)\) 减去 \(d\)，并让所有 \(j\in T\) 的点的 \(L_2(j)\) 加上 \(d\)，对于不同的边会出现什么情况：
- \(S-T\) 或 \(S'-T'\)，没有影响；
- \(S-T'\)，这种边 \(L_1(u)+L_2(v)-w(u,v)\) 值会变小，如果这个值变成 \(0\) 就被加进相等子图中去了，并且可以证明所有相等子图中的边都不可能是这种边，不然交错树会变大；
- \(S'-T\)，这种边 \(L_1(u)+L_2(v)-w(u,v)\) 的值会变大，如果它原来是相等子图中的边，那么它之后就不是了，如果它原来就不是，那它在修改之后就更不可能是了；
现在考虑我们选择的 \(d\) 是什么，由于我们目的是让至少一条边加入到相等子图中去，并且始终满足 \(\sum_{u\in V_1' \\v\in V_2'}L_1(u)+L_2(v)\ge \sum_{e\in E'}w_e\)，那么我们选择的 \(d\) 一定来源于第 \(2\) 种情况，因为只有它能让某条边加入到相等子图中去，并且这个 \(d\) 一定是
之后我们重复这个过程，直到每个点都完成一次增广，这个时候就找到了最大权完美匹配。
时间复杂度：每个点都会被交错树中一次，每拉进去一次就会执行一次尝试增广，所以这里的复杂度是 \(\mathcal O(n(n+m))\)，但是每个点都要求完成一次增广，所以总复杂度是 \(\mathcal O(n\times n(n+m))\)，当 \(m=n^2\) 时，总体复杂度是 \(\mathcal O(n^4)\) 的。
伍、代码（优化前）
using namespace Elaina;
const int maxn=500;
const ll inf=(1<<30)-1;
int e[maxn+5][maxn+5];
int n,m;
inline void input(){
n=readin(1),m=readin(1);
int u,v;ll w;
rep(i,1,m){
u=readin(1),v=readin(1),w=readin(1ll);
e[u][v]=w;
}
}
int vall[maxn+5],valr[maxn+5];
int visl[maxn+5],visr[maxn+5];
int matchl[maxn+5],matchr[maxn+5];
// record the value of expression : vall[u]+valr[v]-weight_edge
int down[maxn+5];
int findAugRoad(const int u){
visl[u]=1;
for(int v=1;v<=n;++v)if(e[u][v] && !visr[v]){
if(vall[u]+valr[v]==e[u][v]){
if(!matchr[v] || findAugRoad(matchr[v])){
matchl[u]=v,matchr[v]=u;
return 1;
}
}
else down[v]=Min(down[v],vall[u]+valr[v]-e[u][v]);
}
return 0;
}
inline void KM(){
// initialize value array
memset(valr+1,0,n<<2);
for(int u=1;u<=n;++u){
vall[u]=-inf;
for(int v=1;v<=n;++v)vall[u]=Max(vall[u],e[u][v]);
}
// enumerate each node as start
for(int s=1;s<=n;++s){
while(1){
memset(visl+1,0,n<<2);
memset(visr+1,0,n<<2);
memset(down,0x3f,n<<2);
// if an augmented road has existed in the equal graph
// the current node has finished augmenting.
if(findAugRoad(s))break;
int delta=inf;
for(int i=1;i<=n;++i)if(!visr[i])
delta=Min(delta,down[i]);
for(int i=1;i<=n;++i){
if(visl[i])vall[i]-=delta;
if(visr[i])valr[i]+=delta;
}
}
}
ll ans=0;
rep(i,1,n)ans+=e[matchr[i]][i];
writc(ans,'\n');
rep(i,1,n)writc(matchr[i],' ');Endl;
}
signed main(){
input();
KM();
return 0;
}
陆、优化
这个算法是 \(\mathcal O(n^4)\) 的，不是十分理想，能不能进行一些优化？
原来的算法中，每次将顶标修改之后，需要使用一次匈牙利进行拓展或者找到当前最小的 \(\min_{u\in S \\v\in T'}\{L_1(u)+L_2(v)-w(u,v)\}\)，匈牙利给我们带来了无穷无尽的复杂度，考虑能否简化这个过程？
对于 \(V_2\) 中的点，我们定义一个数组 \(\tt def[i]\) 表示最少需要花费 \(def[i]\) 的修改将其拉进交错树中，即
定义松弛操作为将 \(V_1\) 中的某个点去取出，遍历其连接的所有的 \(V_2\) 中的点，执行 \(def[i]=\min_{1\le j\le n}\{L_1(j)+L_2(i)-w(i,j)\}\quad (i\in T',j\in S)\).
最开始，\(S\) 中只有起点，然后将起点拿去松弛，然后在所有的 \(def\) 中找到一个最小的，这个值就是 \(d\)，然后进行修改操作，在修改完之后一定会有一个 \(V_2\) 中的点被加入到 \(T\)，如果这个点还没有被匹配，那么我们就找到了一条增广路，只需要跑一遍匈牙利增广这条路即可，如果这个点 \(v\) 在之前被匹配了，即连向 \(v\) 的点想要 \(v\)，但是可惜它已经名花有主，这个时候我们想办法让它的 “主” 找另外的点，这样我们就将它的主人拿去松弛，松弛完又去找 \(d\)，然后重复这个过程...直到某个时刻找到的 \(v\) 没有主人，那么我们就成功找到了一条增广路，然后跑一边匈牙利增广就可以了.
由于我们松弛的过程不修改匹配，所以每个匹配的主人最多松弛一次，松弛复杂度是 \(\mathcal O(n)\) 的，每个点作为起点，找到增广路就执行一次匈牙利并到下一个点，所以整体复杂度是 \(\mathcal O(n((m+n)+n^2))=\mathcal O(nm+n^3)\) 的。
柒、真正的代码
using namespace Elaina;
const int maxn=500;
const int inf=(1<<30)-1;
struct edge{int to,nxt,w;}e[maxn*maxn+5];
int tail[maxn+5],ecnt;
inline void add_edge(const int u,const int v,const int w){
e[++ecnt]=edge{v,tail[u],w};tail[u]=ecnt;
}
int n,m;
inline void input(){
n=readin(1),m=readin(1);
int u,v,w;
rep(i,1,m){
u=readin(1),v=readin(1),w=readin(1);
add_edge(u,v,w);
}
}
int def[maxn+5];
int l1[maxn+5],l2[maxn+5];
/** @brief record the match of the node set V_2 */
int match[maxn+5];
/** @brief whether we have used this pair of match */
int used[maxn+5];
inline void slack(const int u){
for(int i=tail[u],v;i;i=e[i].nxt){
v=e[i].to;
def[v]=Min(def[v],l1[u]+l2[v]-e[i].w);
}
}
int hungary(const int u){
for(int i=tail[u],v;i;i=e[i].nxt)if(!used[v=e[i].to]){
if(l1[u]+l2[v]==e[i].w){
used[v]=1;
if(!match[v] || hungary(match[v])){
match[v]=u;
return 1;
}
}
}
return 0;
}
inline void KM(){
/* initial */
rep(u,1,n){
l1[u]=-inf;
for(int i=tail[u];i;i=e[i].nxt)
l1[u]=Max(l1[u],e[i].w);
}
/* numerate the start point */ ;
rep(s,1,n){
/* initial */ ;
rep(i,1,n)used[i]=0,def[i]=inf;
/* first you should slack point @p s */ ;
slack(s);
while(true){
/* find the point with the least def */ ;
int mini=inf,p;
rep(i,1,n)if(!used[i] && def[i]<mini)
mini=def[i],p=i;
// printf("mini == %d, p == %d\n",mini,p);
if(mini==0){
/* we've found the augmented road which meets our demand */ ;
if(!match[p]){
/* initial the array used as the vis array of hungary */ ;
memset(used+1,0,n<<2);
hungary(s);
break;
}
else{
/* or we're going to try to extend the staggered tree */ ;
/* mark this pair of match */ ;
used[p]=1;
slack(match[p]);
}
}
else{
/* if there's no such point whose def equals to 0 */ ;
/* we're going to adapt l1 and l2 */ ;
rep(i,1,n){
/* a node has been used, which means it has already matched, so we need to update l1 and l2 */ ;
if(used[i])l2[i]+=mini,l1[match[i]]-=mini;
/* pay attention to this if, if you do not want to use this if, please set INF a larger number*/
else if(def[i]!=inf)def[i]-=mini;
}
/* pay attention !! point s is simple, but it exists in the staggered tree */ ;
l1[s]-=mini;
}
}
}
ll ans=0;
rep(i,1,n)ans+=l1[i]+l2[i];
writc(ans,'\n');
rep(i,1,n)writc(match[i],' ');Endl;
}
signed main(){
input();
KM();
return 0;
}
捌、个人理解
以上的算法都是在干一件事：
在相等子图中找到一条增广路，这条增广路的开头由我们指定，结尾是一个无主的节点。如果在当前的相等子图中没有这样的增广路，那么我们就会得到一个交错树，然后我们需要干一些更改可行顶标的事，让我们的交错树变大到在其之内有一条我们想要的增广路，然后就增广.
时间复杂度不同是出现在我们如何更改可行顶标，如果每次都暴力用匈牙利改，那么时间复杂度就是 \(\mathcal O(n^2m)\)，如果使用优化之后的思路，复杂度就是 \(\mathcal O(nm+n^3)\).
其实，优化之后的思路更像是用 \(\tt bfs\) 实现匈牙利，难道不是吗？
