从匈牙利算法到带权带花树——详解对偶问题在图匹配上的应用

Potassium · 2021-09-10 21:48:05 · 题解

友情提示：前往 github blog 可以获得更好的阅读体验 :)

引入

前些天看到了一个比较有趣的题目，需要用到一般图最大权匹配。可是我只会二分图最大匹配，甚至不会 KM 和带花树的原理，于是就进行了一个资料的搜，顺便增长一下板子库。然而——

对一般图最大权匹配，网上现成高质量资料较少。

苦于资料杂乱、且无论正确性还是代码具体实现都没有一个清晰的全流程教程，在吸取众长、奋斗三天后，终于通过了 UOJ 上的模板题。

这篇博客主要目的在于，梳理整个流程（包含代码实现）中应用的原理和使用的技巧。一般图最大权匹配主要通过下面四个问题逐步转化：

1. 二分图最大匹配
2. 一般图最大匹配
3. 二分图最大权匹配
4. 一般图最大权匹配

本文主要描述，如何通过解决前三个问题的思路，解决第四个更为复杂的问题。本文的前置知识只有匈牙利算法。

基本定义

匹配：对于图 G=(V,E) 的一个匹配 M 是边集 E 的子集，其中两两不共用顶点。

最大匹配：边数最多的匹配。

完美匹配：包含原图所有点的匹配。

最大权匹配：对于图 G(V,E)，边 e\in E 的权重为 w(e)，最大权匹配 M=(V',E') 为 \sum_e w(e) 最大的匹配。

最大权最大匹配：最大权最大匹配 M(V',E') 为 \sum_e w(e) 最大的最大匹配。

匹配点、未匹配点、匹配边、非匹配边：字面含义，对应于某个匹配的状态。

交错路径：从一个未匹配点出发，依次交错经过非匹配边、匹配边形成的路径。

增广路：一个末尾是未匹配点的交错路径。这样的路径可以通过翻转整条路上边的选择状态获得一个长度恰好增加 1 的匹配。

下列各图展示了一个无向二分图（箭头只是表示增广路的顺序）某时刻对增广路增广的情况，分别为一个二分图、一条交错路、一条增广路和对增广路增广的结果。（怎么不支持 html 啊，不会横着放，再次传送门获得更好的阅读体验）

交错树：从未匹配点 r 开始寻找增广路时，交错路径组成的树。其中，我们称向根方向匹配的点为 S 点（黑点，偶点），背向根方向匹配的点为 T 点（白点，奇点）。下图展示了二分图和以 1 为根的交错树。

二分图最大匹配

二分图匹配可以用网络流做，Dinic 在二分图上的复杂度和 Hopcroft-Karp 相同，均为 O(E\sqrt{V}) ，不赘述。

匈牙利算法流传较广的是 DFS 做法，即枚举左部每个未匹配点，记录右部节点匹配的左部标号 match_q，检查交错路是否为增广路，复杂度 O(VE)。正确性是通过 Berge's Theorm （如果找不到某点为端点的增广路，则最大匹配可以不包括这个点）证明的，比较显然。

int dfs(int p,int t){ // p: 左半边编号 
    for(auto q:G[p])
        if(t!=vis[q]){
            vis[q]=t;
            if(match[q]==-1||dfs(mt[q],t))
                return match[q]=p,1;
        }
    return 0;
}
int maxmatch(){
    int ans=0;
    memset(match,-1,sizeof(match));
    for(int i=1;i<=n;i++)
        if(match[i]==-1&&dfs(i,i))ans++;
    return ans;
}

同样地，可以进行 BFS，即初始时所有点都未匹配，将未匹配点入队。逐个取出队首元素作为交错树的根，对交错树进行 BFS 匹配。可以看出，增广的本质是在以每个点为根的交错树上寻找增广路。而二分图左部永远是黑点，右部永远是白点，因此增广路的首末必分别在左部和右部。

一般图最大匹配

对于一般图，就没有这么好的性质了。但唯一的不同是，一般图中可能有奇数长度的环，下简称奇环。在奇环上的每一个点连出去的非匹配边，都能够成为某个增广路上的边，这导致在增广的过程中难以判断如何翻转。下图是增广路经过奇环的两种形式。

我们发现，奇环具有了所有黑点需要具备的性质：任意点连出去的非匹配边都可以在从根开始的交错路上；奇环的根向根方向匹配。于是，将奇环缩成一个点，称为花。奇环的根，也就是整个环中唯一连了两个环内非匹配边的点，称为花托（上图中的 root）。

将花缩成一个点后，环上所有点变成黑点，并可以在交错树上向外延伸。

具体实现中，不需要真正地记录花，只需要记录每个点对应花根位置 fa（本质是个并查集）；同时为了方便上述两种情况的判断，对于所有黑点，记录交错树上从根开始连向当前点的非匹配边 pre（对于非花点 i，pre_i 就是父节点）。

每当找到黑-黑边（S-S 边）时进行缩花（shrink）操作，花托就是两点的 LCA ，可以暴力求出；将花内所有白点变黑入队，并更新这些点的 pre 。由于实现比较容易，下面不解释地给出参考代码（后来还是写了一些注释便于理解）。

int mat[N],color[N],n; // color 1: black 2: white
int pre[N]; // walk through a unmatched edge
int fa[N]; // union-find set
vector<int>G[N];
queue<int>Q; 
int find(int x){return fa[x]==x?x:fa[x]=find(fa[x]);}
void rev(int p){ // reverse augment path from rt to p
    // !!!!!!!!!!!! graph needs to be [1,n]
    if(p){
        rev(mat[pre[p]]);
        mat[p]=pre[p];mat[pre[p]]=p;
    }
}
void shrink(int u,int v,int r){
    // shrink odd cycle into 1 point: color all along to black ([x])
    // pre: -> [x] ->
    //             <- [y] <-
    while(find(u)!=r){
        // [pre[mat[u]] <- mat[u] -- [u] -- [v]
        // turn into (-> stands for pre)
        //         [u]  <- [v]    -- []  -> []
        pre[u]=v;
        fa[mat[u]]=fa[u]=r;
        if(color[mat[u]]==2)color[mat[u]]=1,Q.push(mat[u]);
        v=mat[u];u=pre[v];
    }
}
int vis[N];
int lca(int x,int y){
    static int t=0;
    ++t;
    while(1){
        if(vis[x=find(x)]==t)return x;
        else if(x)vis[x]=t;
        x=pre[mat[x]];
        swap(x,y);
    }
}
int augment(int rt){
    int i;
    for(i=0;i<=n;i++)fa[i]=i,color[i]=pre[i]=0;
    color[rt]=1;
    Q=queue<int>({rt});
    while(!Q.empty()){
        p=Q.front();Q.pop();
        assert(color[p]==1); // black
        for(auto q:G[p]){
            if(color[q]==1){
                int r=lca(p,q); // root of flower
                shrink(p,q,r);
                shrink(q,p,r);
            }
            else if(!color[q]){
                color[q]=2; // white
                pre[q]=p;
                if(!mat[q]){
                    rev(q);
                    return 1;
                }else{
                    color[mat[q]]=1;
                    Q.push(mat[q]);
                }
            }
        }
    }
    return 0;
}
int maxmatch(){
    int i,m;
    scanf("%d%d",&n,&m);
    int ans=0;
    for(i=0;i<m;i++){
        int x,y;
        scanf("%d%d",&x,&y);
        G[x].pb(y);G[y].pb(x);
        if(!mat[x]&&!mat[y])mat[x]=y,mat[y]=x,ans++;
    }
    for(i=1;i<=n;i++)if(!mat[i])ans+=augment(i);
    return ans;
}

二分图最大权匹配

和二分图最大匹配一样，最大权匹配也可以通过网络流那一套求解：将每条边的边权看做一个费用，然后求一下最大费用最大流即可。但为了扩展到一般图，考虑 KM 算法。

学之前，我对二分图最大权匹配的印象仅停留在 KM 算法的板子上，没有真正去理解 KM 算法。下面简单介绍一下 KM 算法（如果你是第一次学习 KM 算法，推荐一个比较有趣的引入：KM算法详解+模板，仅看找对象部分大致了解一下即可）。

最大权完美匹配

我们先考虑最大权完美匹配（最大权匹配可以通过增加零边变成完美匹配）。设对于边 e，选/不选该边表示为 x_e=0 或 x_e=1；设 \delta(v) 表示 v 所连的边集，则我们要求的最大权匹配即为

这是一个整数规划问题，通过原始对偶（Primal-Dual）思想转化为对偶问题：

同时，根据互补松弛条件， x_e>0\Rightarrow z_e=0。即，对于选中的边 e，必有 z_e=0。

这里的 z 就是 KM 算法中的顶标（vertex labeling）。定义边 e(u,v) 为”等边“当且仅当 z_u+z_v=w(e)，也即 z_e=0，根据上述互补松弛条件，匹配只能在等边构成的子图上建。初始时，左部（黑点）每个点的顶标均为其能够连到的最大边权，右部（白点）每个点的顶标均为 0。与匈牙利算法类似，每次增广先在“等边”构成的一张二分子图上进行尝试匹配，如果匹配不到，则松弛（将标准放低一些）后继续尝试；直至所有边都成为等边仍未找到增广，算法结束。

松弛的目的是使一些新边成为等边，加入二分子图，让左部点能尝试匹配的点更多；同时保证已有匹配边 z_e=0 。对于所有未匹配边，求出最小的 delta=\min_e z_e=z_u+z_v-w(e)，然后将交错树中黑点的顶标降低 delta，交错树中白点的顶标增加 delta。这样暴力做的复杂度是 O(n^4)，下面给出参考代码。

暴力实现

const ll inf=0x3f3f3f3f3f3f3f3f;
ll w[N][N],lx[N],ly[N];
int mat[N],visx[N],visy[N],n;
int dfs(int x){
    visx[x]=1;
    int y;
    for(y=1;y<=n;y++){
        if(!visy[y]&&lx[x]+ly[y]==w[x][y]){
            visy[y]=1;
            if(!mat[y]||dfs(mat[y]))return mat[y]=x,1;
        }
    }
    return false;
}
ll km(){
    int i,x,y;
    memset(lx,0,sizeof(lx));
    memset(ly,0,sizeof(ly));
    memset(mat,0,sizeof(mat));
    for(x=1;x<=n;x++)
        for(y=1;y<=n;y++)
            lx[x]=max(lx[x],w[x][y]);
    for(i=1;i<=n;i++){
        while(1){
            memset(visx,0,sizeof(visx));
            memset(visy,0,sizeof(visy)); 
            if(dfs(i))break;
            ll d=inf;
            for(x=1;x<=n;x++)
                for(y=1;y<=n;y++)
                    if(visx[x]&&!visy[y])
                        d=min(d,lx[x]+ly[y]-w[x][y]);
            if(d==inf)break;
            for(x=1;x<=n;x++)if(visx[x])lx[x]-=d;
            for(y=1;y<=n;y++)if(visy[y])ly[y]+=d;
        }
    }
    ll ans=0;
    for(i=1;i<=n;i++)ans+=w[mat[i]][i];
    return ans;
}

不太行，考虑优化。

优化方案

一种可行的方案是，对于所有右部点 y，记录 slv_y=\text{argmin}_{x\in T}(z_x+z_y-w(x,y)) ，即固定了 y 后， z_e 最小值对应 x ，其中 T 为交错树。每次扩展的时候，都是取 slv_y 计算出结果最小的 y 进行扩展，这样就能够保留已有的交错树进而继续增广，每次增广的复杂度变成 O(n^2)，总复杂度 O(n^3)。下面给出示例代码，其中为方便增广，使用 pre 记录黑点在交错上父亲的父亲（最近黑点祖先）。

const ll inf=0x3f3f3f3f3f3f3f3f;
ll w[N][N],lx[N],ly[N];
int mx[N],my[N],visx[N],visy[N],n;
int slv[N],pre[N];
queue<int>Q;
ll calc_slv(int u,int v){return lx[u]+ly[v]-w[u][v];}
void add_to_tree(int x){ // add x to tree, update all slv[y]
    Q.push(x);
    visx[x]=1;
    int y;
    for(y=1;y<=n;y++)
        if(!slv[y]||calc_slv(x,y)<calc_slv(slv[y],y))slv[y]=x;
}
void link(int x,int y){
    int tx,ty;
    while(x){
        tx=pre[x];ty=mx[x];
        my[y]=x;mx[x]=y;
        x=tx,y=ty;
    }
}
int augment(int p){
    int x,y;
    while(1){
        while(!Q.empty()){
            int x=Q.front();Q.pop();
            for(y=1;y<=n;y++)if(!visy[y]&&!calc_slv(x,y)){
                visy[y]=1;
                if(!my[y])return link(x,y),1;
                pre[my[y]]=x;
                add_to_tree(my[y]);
            }
        } 
        ll d=inf;
        for(y=1;y<=n;y++)
            if(!visy[y])
                d=min(d,calc_slv(slv[y],y));
        if(d==inf)break;
        for(x=1;x<=n;x++)if(visx[x])lx[x]-=d;
        for(y=1;y<=n;y++)if(visy[y])ly[y]+=d;
        // add new edges
        for(y=1;y<=n;y++)if(!visy[y]&&!calc_slv(slv[y],y)){
            visy[y]=1;
            if(!my[y])return link(slv[y],y),1;
            pre[my[y]]=slv[y]; // tree: slv -> y == my[y], == stands for match
            add_to_tree(my[y]);
        }
    }
    return 0;
}
ll km(){
    int i,x,y;
    memset(lx,0,sizeof(lx));
    memset(ly,0,sizeof(ly));
    memset(mx,0,sizeof(mx));
    memset(my,0,sizeof(my));
    for(x=1;x<=n;x++)
        for(y=1;y<=n;y++)
            lx[x]=max(lx[x],w[x][y]);
    for(i=1;i<=n;i++){
        memset(visx,0,sizeof(visx));
        memset(visy,0,sizeof(visy)); 
        memset(slv,0,sizeof(slv)); 
        Q=queue<int>();
        add_to_tree(i);
        augment(i);
    }
    ll ans=0;
    for(i=1;i<=n;i++)ans+=w[my[i]][i];
    return ans;
}

Bonus

还有一种使用同样思想，利用局部最优获得全局最优的实现方式。对于每一个黑点 x，找到某个白点加入交错树，这个过程中始终是选择 z_e 较小的边进行扩展的，直到找到 x 所对应的增广路。

假设当前已经处理到 x，1\dots x-1 都已经有对应的增广路了。设对于未匹配白点 y，与所有已匹配黑点 x 中的 z_e 最小值为 slack_y，有 slack_y=\min_{i\in[1,x)}z_i+z_y-w_{i,y}。那么加入交错树的 yy=\text{argmin}_yslack_y；同时，根据 yy 的 z_{(i,yy)} 去更新黑白已匹配点的顶标和下一轮的 slack，同时将 mat_{yy} 这个黑点加入交错树（如果 yy 没有匹配那就说明增广路找着了，退出循环）。

这样的方法对于顶标的初始值没有要求（可以不初始化），每次得到的是 x 的局部最优解，但由于可以进行增广，因此最后得到的也是全局最优解，非常巧妙。下面给出代码。其中 pre 记录的是白点在交错树上父亲的父亲（最近白点祖先），visy 记录白点 i 是否在交错树中。这个实现常数略小一点。

const ll inf=0x3f3f3f3f3f3f3f3f;
ll a[N],p[N],b[N],c[N];
ll w[N][N],slack[N],lx[N],ly[N];
// for right points
//   pre: last right point in augment path
//   visy: if i is in augment path
//   left: left match
int pre[N],visy[N],left[N]; 
int n;
void augment(int p){
    int i,x=p,y=0,yy;
    memset(pre,0,sizeof(pre));
    memset(slack,0x3f,sizeof(slack));
    left[0]=p;
    while(1){
        x=left[y],visy[y]=1;
        ll delta=inf;
        for(i=1;i<=n;i++){
            if(visy[i])continue;
            ll d=lx[x]+ly[i]-w[x][i];
            if(slack[i]>d)slack[i]=d,pre[i]=y;
            if(slack[i]<delta)delta=slack[i],yy=i;
        }
        for(i=0;i<=n;i++){
            if(visy[i])lx[left[i]]-=delta,ly[i]+=delta;
            else slack[i]-=delta;
        }
        y=yy;
        if(left[y]==-1)break;
    }
    // reverse augment path
    while(y)left[y]=left[pre[y]],y=pre[y]; 
}
ll km(){
    int i,j;
    memset(left,-1,sizeof(left));
    memset(lx,0,sizeof(lx));
    memset(ly,0,sizeof(ly));
    for(i=1;i<=n;i++){
        memset(visy,0,sizeof(visy));
        augment(i);
    }
    ll res=0;
    for(i=1;i<=n;i++)if(left[i]>0)res+=w[left[i]][i];
    return res;
}

其他匹配的形式

上述讨论的都是左右部大小相等的完全二分图上的最大权完美匹配；对于 n_1\neq n_2 的情况，可以补齐不足元素进行计算；对于不要求完美匹配的最大权匹配，可以直接应用本模板；对于非完全二分图，增广过程中对不存在的边进行特判即可。

一般图最大权匹配

仍然考虑 KM 算法在带花树上跑的过程。当进行松弛操作时，如果仍然按照上面的处理方式，白点 +delta，黑点 -delta，那么对于一个缩成点的花来说，其中有很多交错树上的边 (u,v)，u 和 v 的顶标都被加或者减掉了 delta，导致这条边不满足 z_e\ge 0 的条件。于是需要对于整个花 B，设置额外的花标 z_B，在松弛时额外对所有黑花 +2delta，对白花 -2delta，才能保持平衡。

对偶问题

设 O 为大小为 \ge 3 奇数的集合的集合（包含所有花），\gamma(S) 表示 S 集合中的边，这个过程仍然可通过线性规划转化为对偶问题：

对偶问题

与二分图最大权匹配相同，根据互补松弛条件， x_e>0\Rightarrow z_e=0。即，对于选中的边 e，必有 z_e=0；除此以外，还有 z_B>0\Rightarrow x(\gamma(B))=\lfloor\frac{|B|}2\rfloor，即所有 z_B>0 的集合 B，都被选了集合大小一半的边，也即 B 是一朵花。同时，我们加入一个条件： z_B>0\Rightarrow x(\delta(B))=1，即只有花 B 向外连了一条边的时候， z_B>0 才是有意义的。这很容易根据 z_B 的用途理解。

具体实现

因此，我们需要维护的就是：

1. 任意边 z_e\ge 0，取等时是交错树里（包括形成花的边）的边；

通过类似 KM 的办法，我们就可以容易地维护带花树了！本文完

难点就在它的实现步骤。这里我分享一下通过借鉴一派 UOJ 神仙提交代码，提纯结晶后得到的模板实现思路与细节。min25 的板子真心看不懂，有懂哥欢迎补充！

先摆一下是如何存储这张带花树的：使用 edge 类型存边，大有讲究！对于有花参与的连边，edge 存储的 u,v 不是花本身，而是花内真实节点的边。

#define N M*2+1
struct edge{
    int u,v;T w;
    edge(){}
    edge(int u,int v,T w):u(u),v(v),w(w){}
};
// Graph
int n,n_x; // [1, n]: point; [n+1, n_x]: flower
edge g[N][N]; // adjacent matrix
// flower
vector<int>flower[N]; // nodes in flower i (outer flower)
int root[N]; // flower root, root<=n root=i: normal nodes
int flower_from[N][N]; // flower_from[b][x]: outermost flower in b that contains x

其中，与普通带花树不同，这里的花是额外进行存储的（存储在 n+1 到 n_x 内），因此需要开两倍空间。花可以嵌套，而且花需要在某些时刻进行展开（后面会详细讲述），因此需要保存额外的信息。对于一个花 B，存储其根 root，对于花内任意点有 root_{i\in B}=B；存储花内所有节点 flower_B，是以花托为起始的环；同时，花内每一点 i 保存“花内最大父花” flower\_from_{B,i} 表示最大的包含 i 的 B 的子花。

对于上面这朵嵌套花：\{1,2,3\}\in B1,\{B1,4,5,6\}\in B2，存储为：

flower[B1]={1,2,3} // 1 是花托，按照某个顺序转圈，也可能是 {1,3,2}
flower[B2]={6,B1,4,5} // 6 是花托，按照某个顺序转圈，也可能是 {6,5,4,B1}
root[1]=root[6]=root[B1]=root[B2]=B2 // i 最外层的花托
flower_from[B2][B1]=B1
flower_from[B2][1]=B1
flower_from[B2][4]=4
flower_from[B1][1]=1
root[7]=7 // 非花点 root 为自身

然后是 slack 和匹配相关：

// slack
T label[N]; // node label, [1, n] point label, [n+1, n_x] flower label
int col[N]; // color saved at flower root
int slv[N]; // slack node of NON LEAF NODES, slv[y]=x z(x,y) min_x
// match
int mat[N]; // match, mat[x]=y (x,y)\in E
int fa[N]; // fa in cross tree
int vis[N]; // if in path

queue<int>Q; // bfs queue