网络单纯形算法学习笔记

参考链接

网络单纯形算法解说【翻译】 - zhiyin123123 - 博客园

前言

阅读本篇，你不需要了解单纯形算法，但是你需要了解费用流相关概念。

如果你会使用 SPFA + EK 算法（一种 SSP 算法）求解费用流问题，并想学会更高效的算法，那么你适合阅读本篇。

约定

• 如无特殊说明，使用 \(n,m\) 代表图中的点数和边数，用 \(S,T\) 代表源点和汇点。

介绍

网络单纯形算法可以直接解决如下问题：

• 有源汇最小费用最大流问题。图中可能有负圈。

它能做到期望 \(O(nm)\)（常数极小，类似于 SPFA）或者 \(O(m\log n)\)（要写 LCT）的复杂度。通常我们实现 \(O(nm)\) 的版本。不过，据说该算法能被卡到指数，如果实现错误还可能引起死循环。

介绍以下大致思路。首先，连接一条 \(S\to T\) 的边，容量为 \(\inf\)，费用为 \(\inf\)，并钦定让其流满。在接下来的讨论中，我们认为，图中所有边都有虚拟反边，对于边 \(u\to v\)，容量为 \(\text{cap}\)，费用为 \(\text{cost}\)，有虚拟反边 \(v\to u\)，容量为 \(\text{cap}\)，费用为 \(-\text{cost}\)；每当一条边流量增加 \(f\) 时，其对应反边流量自动增加 \(-f\)。然后，重复进行如下操作：

• 在图中找到一个负圈，使得其所有边均能增加流量，然后找到最大的流量 \(f\)，让负圈上的所有边的流量增加 \(f\)，并均不超出限制。

直到找不到负圈为止。最后，把一开始添加的 \(S\to T\) 边和其虚拟反边删掉，得到的图就是答案。

该思路中，添加的辅助边 \(S\to T\) 将费用和流量联系在了一起，若要流这条边，费用将会昂贵，不流这条边，其容量又太大。

该算法的正确性基于：费用流的局部最优解就是全局最优解。

正文

维护生成树解

为了方便讨论，先做如下定义。

• 对于边 \(u\to v\)，定义它和其对应的虚拟反边构成一个边对。
• 对于边对 \((u\to v,v\to u)\)，若 \(u\to v\) 的流量可以增加（而不超过限制），\(v\to u\) 的流量也可以增加，则称其为自由边对。
• 对于一个边对组成的圈，如果组成它的所有边对都是自由边对，则称其为自由圈。

在算法的过程中，我们会不让自由圈出现，因为总有一种推流方式，让该圈的费用变得更小。

我们要维护一个该图的生成树（由边对构成）（不妨设图联通），满足

• 生成树外的边对均不是自由边对。

这样，就一定没有自由圈啦。（因为生成树无圈）

为了方便，我们强制认为，生成树上的边对 \(e\) 可以向两个方向推流，即使不能，我们也可以认为它能向某个方向推流 \(0\)。我们还认为，生成树以外的边对，只有一个推流方向。

然后，我们就要不断在图中找负圈了。设生成树为 \(\text{Tree}\)，则对于任意边对 \(e\)，我们很容易判断 \(\text{Tree}\cap e\) 中是否有负圈。至于如何判断，后面再讲。

设图为 \(G\)，现在我们想证明一个引理：

• 图 \(G\) 中无负圈，当且仅当对于所有边对 \(e\)，图 \(\text{Tree}\cup e\) 中无负圈。

如果它是正确的，我们就可以不断的消除和树有关的负圈来正确求解了。

证明其实很简单，画画图写几个不等式就能搞定。

假如在算法的执行过程中，我们找到了边对 \(e\)，使得 \(\text{Tree}\cup e\) 中有负圈，那么，我们直接将这个圈流满，然后将 \(e\) 加入 \(\text{Tree}\) 中，再在该圈中从 \(\text{Tree}\) 删除一条非自由边对，然后就有能维持 \(\text{Tree}\) 是一颗生成树了。

我们不断找到这样的边 \(e\)，然后执行如此推流操作，直到找不到 \(e\)，算法的流程也就如此了。

维护生成树的具体实现

在维护生成树时，首先要解决的问题就是，对于任意一条边对 \(e\)，如何判断 \(\text{Tree}\cap e\) 中是否有负圈。这可以树上差分。

不妨先给 \(\text{Tree}\) 钦定一个根节点 \(\text{root}\)，然后，我们维护每个点 \(u\) 在树 \(\text{Tree}\) 中到 \(\text{root}\) 的距离 \(d_u\)（一条边贡献的距离为其费用）。需要注意的是，这个距离是有向的，\(\text{root}\) 到 \(u\) 的距离是 \(-d_u\)。然后，对于边 \(u\to v\)，设其费用为 \(\text{cost}\)，它和 \(\text{Tree}\) 组成的圈的费用和为

\[\text{cost}-d_u+d_v \]

然后，我们还要维护从树中加边和删边，用 LCT？太难写了，直接单次 \(O(n)\) 暴力维护。

还有，对于一条边，我们要找出它和生成树构成的圈包含哪些边，这样才能推流。可以使用单次 \(O(n)\) 的暴力 LCA 维护。

如何找到边使得 \(\text{Tree}\cap e\) 有负圈

实际上，有如下可行方法：

1. 循环遍历所有边对 \(e\)，如果发现 \(\text{Tree}\cap e\) 中有负圈，那么就立刻选择它推流。
2. 暴力遍历所有边对 \(e\)，找到 \(\text{Tree}\cap e\) 中负圈费用最小的，用它推流。
3. 将边对分为 \(B\) 组（\(B\) 可能取 \(\sqrt m\)），循环遍历每个组，在组中选择产生负圈权值最小的边，用它推流。

选第 \(1\) 个即可，下面两个方法是🤡，跑的更慢。

如何避免死循环

容易发现，这个算法有可能会产生死循环，具体来讲，有如下两个原因：

1. 有的时候，推流减少的费用可能为 \(0\).
2. 用 \(e\) 推流并加入 \(\text{Tree}\) 中后，有可能有多条圈中的边对是非自由边对，不知道删哪一个。

下面是解决方案，不过我不会证明其正确性：

• 用 \(e\)（设其为 \(u\to v\)）推流并加入 \(\text{Tree}\) 中后，设 \(u\) 和 \(v\) 在 \(\text{Tree}\) 中的 LCA 为 \(\text{lca}\)，则想象你从 \(\text{lca}\) 沿着圈上的边走到了 \(u\)，再走到 \(v\)，再走到了 \(\text{lca}\)，如此绕一圈，然后，删掉你遇到的第一个非自由边对。

该算法能爆切的题目

P7173 【模板】有负圈的费用流 - 洛谷

暴力建图就能过，还比正解跑得快：P5331 [SNOI2019] 通信 - 洛谷

代码实现

#include<bits/stdc++.h>
template<int Node_cnt,int Edge_cnt,typename cost_t,typename flow_t>
class min_cost_max_flow_t{
    template<typename Tp>
    static bool tomin(Tp &x,const Tp &y){if(y<x){x=y; return 1;} return 0;}
    static constexpr int NSI=Node_cnt,ESI=(Edge_cnt+1)*2;
    struct ed_t{int nxt,to; flow_t f; cost_t c;};
    int n,m,S,T;
    std::vector<char> vis;
    std::vector<int> hd,Fa,Fe,dep,buf;
    std::vector<cost_t> dis;
    std::vector<ed_t> E;
    void build(int it,int fa,int fe){
        vis[it]=1;
        Fa[it]=fa; Fe[it]=fe;
        for(int j=hd[it];j!=0;j=E[j].nxt)if(!vis[E[j].to]){
            build(E[j].to,it,j^1);
        } 
        return ;
    }
    void upd(int it){
        if(it==0||vis[it]) return ;
        vis[it]=1;
        upd(Fa[it]);
        dep[it]=dep[Fa[it]]+1;
        dis[it]=dis[Fa[it]]+E[Fe[it]].c;
        return ;
    }
    int getlca(int x,int y){
        while(x!=y) (dep[y]>=dep[x])?y=Fa[y]:x=Fa[x];
        return x;
    }
    void pushflow(int j){
        int l=NSI,r=NSI;
        buf[l]=j;
        int u=E[j^1].to,v=E[j].to;
        int lca=getlca(u,v);
        for(int i=u;i!=lca;i=Fa[i]) buf[--l]=Fe[i]^1;
        for(int i=v;i!=lca;i=Fa[i]) buf[++r]=Fe[i];
        int del=NSI; cost_t f=E[j].f;
        for(int i=l;i<=r;i++) if(tomin(f,E[buf[i]].f)) del=i;
        for(int i=l;i<=r;i++){int s=buf[i]; E[s].f-=f; E[s^1].f+=f;}
        if(del==NSI) return ;  
        if(del<NSI){
            for(int i=del+1;i<=NSI;i++){
                int s=buf[i];
                Fa[E[s^1].to]=E[s].to;
                Fe[E[s^1].to]=s;
            }
        }else{
            for(int i=NSI;i<=del-1;i++){
                int s=buf[i];
                Fa[E[s].to]=E[s^1].to;
                Fe[E[s].to]=s^1;
            }
        }
        std::fill(vis.begin()+1,vis.begin()+n+1,0);
        return ;
    }
public:
    min_cost_max_flow_t():n(0),m(1),S(0),T(0),
        vis(NSI+5,0),hd(NSI+5,0),Fa(NSI+5,0),Fe(NSI+5,0),dep(NSI+5,0),buf(NSI*2+5,0),
        dis(NSI+5,0),E(ESI+5){}
    void connect(int u,int v,flow_t f,cost_t c){
        E[++m]={hd[u],v,f,c}; hd[u]=m;
        E[++m]={hd[v],u,0,-c}; hd[v]=m;
        return ;
    }
    std::pair<flow_t,cost_t> solve(int _n,int _s,int _t){
        n=_n; S=_s, T=_t;
        flow_t sflow=0; cost_t scost=0;
        for(int j=2;j<=m;j++){
            auto &e=E[j];
            sflow+=e.f; if(e.c>0) scost+=e.c;
        }
        connect(S,T,sflow,scost+1); 
        std::swap(E[m-1].f,E[m].f);
        flow_t flow=0; cost_t cost=0;
        build(S,0,0); std::fill(vis.begin()+1,vis.begin()+n+1,0);
        while(1){
            bool ok=0;
            for(int j=2;j<=m;j++){
                auto &e=E[j];
                if(e.f>0){
                    int u=E[j^1].to,v=e.to;
                    upd(u); upd(v);
                    if(e.c-dis[u]+dis[v]<0){
                        pushflow(j);
                        ok=1;
                    }
                }
            }
            if(!ok) break;
        }
        for(int j=3;j<m;j+=2){
            auto &e=E[j];
            if(e.to==S) flow+=e.f;
            if(E[j^1].to==S) flow-=e.f;
            cost+=-e.c*e.f;
        }
        return {flow,cost};
    }
};