【魔术】最大流相关

phigy · 2023-03-18 09:06:51 · 算法·理论

0.前言

本文总长度约 6000 字（较短）

你需要熟练掌握 Dinic 及预流推进算法及其特点。（OI-wiki，洛谷日报#56，洛谷日报#127）

同时理解预流推进中 BFS，GAP等优化的价值。（本文代码若包含则不再复述）

由于相关内容缺乏中文资料，可能存在一定错误。

对于过于繁琐的证明，本文给出论文链接。

本文将对最大流进行一个有序的简单探究，并更正许多普遍的认知错误，所有内容最终指向 Enhanced LMES 算法。

1.常见的一些预流推进算法

[1] A New Approach to the Maximum-Flow Problem 1988

1-1.通用预流推进 Generic Push Relabel

极大部分的博客将FIFO预流推进当作了通用预流推进算法。

实际上通用预流推进算法是更加“通用”的。

你只需要做到每次选一个溢出节点然后推流即可。

由于非饱和推流次数不会超过 2n^2m 次。

因此通用预流推进的时间复杂度为 O(n^2m)。

1-2.先进先出预流推进 FIFO Preflow Push

这是我们所熟知的最简单的预流推进算法。

即将有超额流的推送目标用队列存储作为之后的被推送点以免去遍历时间并让每次推送选择的点离源点较近。

论文第四部分详细讲了复杂度的证明方式。

简而言之由于通过队列的次数不会超过 4n^2 次。

因此先进先出预流推进的时间复杂度为 O(n^3)。

评测记录

1-3.后进先出预流推进 LIFO Preflow Push

这可能是本文最实用的内容了。

顾名思义即将先进先出预流推进中的队列存储改为栈存储。

虽然时间复杂度仍然是 O(n^3)，但在绝大部分数据下，因其每次推送选择的点离汇点较近，而拥有良好的实际运行效率，能够超过使用优先队列的HLPP，且写起来简单，推荐使用。

评测记录

1-4 动态树优化的预流推进 Dynamic Trees optimized Preflow Push

论文第五部分 Tarjan 介绍了一种使用动态树让平均每次非饱和流推送操作复杂度低于 O(1) 的办法。

通过使用适当的树操作，我们可以沿着树中的整个路径推流，要么导致饱和，要么将超额流从树中的某个点一直移动到树根。

我们能够证明一个顶点变成溢出点被存储的次数是 O(nm+\frac{n^3}{k})。在每个树操作的成本为 O(\log k) 的情况下，算法的总运行时间为 O((nm+\frac{n^3}{k})\log k)，对于选择 k=\frac{n^2}{m}，其被最小化到 O(nm\log (\frac{n^2}{m})) 的常数因子内。

1-5 最高标号预流推进 Highest Label Preflow Push

[2]Analysis of preflow push algorithms for maximum network flow 1988

相信大家已经对 HLPP 很熟悉了。

但是几乎所有人的实现方法都是将先进先出预流推进中的队列存储改为优先队列存储，从而每次能够得到标号最高的节点推流。

但是实际上我们可以通过对每个高度开一个链表存储节点实现 O(1) 寻找溢出最高点。（具体实现相关文章）

Cheriyan and Maheshwari (1988) 证明了在每次用最高标号的节点推流时，非饱和推流可以降至 n^2\sqrt{m} 次

因此时间复杂度为 O(n^2\sqrt{m})。

评测记录

2.Binary Blocking Flow

二进制阻塞流算法？

The Binary Blocking Flow Algorithm

[8]Beyond the flow decomposition barrier 1998

阻塞流指一条使得不存在一条 S 到 T 的路径使得每条边都没有满流的流。

Dinic 算法的全称为 Dinic 阻塞流算法，它与 EK 算法最大的区别就是它是阻塞流算法，它利用每次在分层图上进行阻塞增广，即会使分层图的源点与汇点不再连通的增广，这样重新构建分层图时最短路会增大从而限制时间复杂度。

关于 LCT 优化的 Dinic

在 Binary Blocking Flow 这个算法中我们通过二进制拆位实现缩点继而在无环图上做阻塞流而达到了一个优秀的复杂度。

我们首先定义 d(i) 为 i 到 t 的距离。

我们认为 (i,j) 是特殊边当 2\Delta\leq u_f(i,j)<3\Delta,d(i)=d(j),u_f(j,i)\ge 3\Delta。

我们定义另一个函数 -\ell 在特殊边上为 0，在其它边上等于 \ell。

注意这并不会改变距离函数 d_\ell=d_{-\ell}。

我们每一轮开始先计算 d 与 -\ell，注意二进制长度函数为 0 的边可能会形成 d 相同点组成的环，我们把 0 边组成的强联通分量缩掉，这样就可以得到一个无环图，只要我们对缩掉的点的流量改动不超过 \Delta，那么扩展时初始图可以重新维持流量守恒。

我们在无环图上增广，每次要么计算出阻塞流，要么计算出流量为 \Delta 的流。

和 Dinic 一样二进制阻塞流算法的时间复杂度也是基于每轮迭代源点到汇点的最短路会增大性质。

并且每一轮中在 O(\sqrt m) 次找到流量为 \Delta 的增广后 F 变小，在 O(\sqrt m) 阻塞增广后 F 变小。

那么我们就有了完整的算法流程：

• 计算 \ell,d。
• 计算 -\ell。
• 缩点。
• 计算阻塞流或计算流量为 \Delta 的流。
• 扩展回初始图。

对于其中第四条，我们有 O(m\log(n^2/m)) 阻塞流算法与 O(\min(n^{2/3},m^{1/2})m) 计算单位容量最大流的方法（就是 Dinic）。

而对于轮数，实质上我们只有 O(\sqrt m\log U) 次,最终得到 O(\min(n^{2/3},m^{1/2})m\log(n^2/m)\log U)。

3.Orlin 的 O(nm) 改进

好像没起名字。

[7]Max flows in o(nm) time, or better 2012

我们沿用 Binary Blcking Flow。

在这个算法中我们通过超大边缩点与对边种类的划分带来的紧密网络而达到了一个优秀的复杂度。

3-1. 超大边缩点 abundance graph

定义 \Delta 是从 s 到 t 的最大残余流量的上限。

我们认为一个边是超大边当它 \ge2\Delta

如果存在超大边 (i,j),(j,i)，那么可以将节点 i 和 j 缩为一个节点，并在缩小后的图中寻找最大流。

在缩小后的图中找到最大流之后，可以将缩点扩展为它们原来的边 (i,j),(j,i)。在扩展后的图中，可以通过在 (i,j),(j,i) 上推送流量来平衡节点 i 和 j 中的流量，从而使得流量的分配合法。

我们在下图中以 (5,6),(6,5) 为例来说明这种缩点。在缩点后，每个边中流量的总变化都小于 2\Delta。

在最终展开标记为 5-6 的节点时，节点 5,6 的流量守恒可能会被破坏，但破坏的程度小于 2\Delta。通过从 5 到 6 或从 6 到 5 发送流量，可以重新保证流量守恒。

另外，也可以缩掉超大的外部弧（指连着超级源点或超级汇点）。我们在下图中以边 (s, 1) 和 (3, t) 为例来说明这种缩点。

在最终展开标记为 s-1 和 3-t 的节点时，节点 1 或 3 可能不满足流量守恒。可以通过在 (s, 1) 和 (3, t) 中推送流量来重新保证流量守恒。

缩点的总时间复杂度为 O(m)。 缩点后的扩展操作时间复杂度也是 O(m)。

3-2 牢不可破的网络 compact network

非常紧密网络有两部分组成。

第一部分是入边残量为正的非超大边的节点的子集。

第二部分为又两部分的并，一部分是第一部分的点间的属于网络边，称为原始边，容量为初始容量；另一部分是第一部分的点对间不属于网络的边，称为伪边，容量为 2\Delta。

接下来我们来得到紧密网络。

我们这个算法能够做到 O(nm) 就是利用在紧密网络中寻找近似最大流比在常规网络中找最大流快来实现。

外部弧是反超大边若它不是超大边。

我们认为一个点是可紧的当它在缩点后的网络中的入流量和比出流量和大的不超过 \Delta/(16mn)且不出中边。

可紧点的定义帮助我们获得紧密网络，若一个点是可紧的且每个入反超大边的容量小于 \Delta/(16mn)，那么它就是特别可紧的。

我们建立一个网络存储那些不是特别可紧的点。

这些不是特别可紧的点被称为临界点，紧密网络中都是临界点。

紧密网络的最大流最多损失 \Delta/(8m) 流量。

3-3.时间复杂度

设 c 为紧密网络点个数。

若 c>m^{9/16} 每轮的时间为 O(m^{3/2}\log^2 n) 乘以 c/m^{9/16} 为 O(cm^{15/16}\log^2n)。

若 m^{1/3}\leq c\leq m^{9/16} ，紧密网络的边为 O(c^2)，在紧密网络中使用上述方法求近似最大流时间为 O(c^{8/3}\log n) ，即 O(cm^{15/16}\log^2n)。

若 c<m^{1/3} ，O(c^3)=O(cm^{2/3})。

接下来我们考虑以下的其它过程：

• 缩点扩点 O(m^{\frac{5}{3}})。
• 建立紧密网络的伪边 O(nm)
• 将伪边流转化为初始边流 O(nm)
• 在紧密网络中将流量从反超大边推送到伪边 O(m^{\frac{5}{3}}\log n)
• 将紧密网络的反超大伪边流转换为初始边流。 O(m^{\frac{5}{3}}\log n)

其中 c 是 O(m) 的，因此总时间复杂度 O(nm+m^{\frac{31}{16}}\log^2n)。

4.超额缩放算法 Excess Scaling

[3]A FAST AND SIMPLE ALGORITHM FOR THE MAXIMUM FLOW PROBLEM 1988

回到了我们的预流推进。

4-1.方法

Ahuja and Orlin 在 1988 年提出使用名为 Excess Scaling 的技巧把通用预流推进非饱和推流的次数从 O(n^2m) 降至 O(n^2\log U)。

类似 Binary Blcking Flow 算法对值域的缩放，我们进行 K=\lceil \log U\rceil +1 轮 Scaling，在每一轮中，我们定义 \Delta 为比所有超额流都大的最小 2 的幂，每一轮 Scaling 结束时 \Delta 减半。

在一次 Scaling 中，我们只考虑超额流大于\Delta/2的节点进行推流，并选择其中标号最小的节点，确保推到一个拥有较小超额流的节点去。

在 8n^2 次非饱和推流之后，不再有超额流大于 \Delta/2，一次 Scaling 结束。

在最多 K 次 Scaling 操作后，所有的节点超额流都会降至 0，至此我们可以得到最大流。

具体的，我们使用链表对每个标号的超额流足够大的点进行存储，每次去有点的最小标号中任意点推一次流，并更新当前弧。

当一个点无法推流时重贴标签并重置当前弧。

本算法复杂度正确性基于无重边和自环，因此预处理合并重边，去掉自环。

4-2 伪代码

4-3 实现相关

该算法饱和推流一共 O(nm) 次，非饱和推流一共 O(n^2\log U) 次。

时间复杂度为 O(nm+n^2\log U)，并可以调整 \Delta/ \beta 的 \beta\ge 2 这个参数变为 O(nm+\beta n^2\log_{\beta}U)。

该算法按 OI 的简单化实现方式常数巨大无比，千万不要考场使用。

评测记录(重测变慢了:( )

[7.5评测记录]

在 n=1000,m=1000000 的随机数据下 Excess Scaling 为 15s HLPP 为 600s

#include <iostream>
#include <queue>
#include <cstring>
#include <cmath>
#include <cctype>
#include <cstdio>
#include <cstring>
#include <map>
#include <ctime>

#define int long long

using namespace std;

int res=0;

int K,level;
const int inf=(long long)1<<40;
int n,m,s,t;
struct edge
{
    int to,val,next;
}e[1200005];
int cnt=1,head[1200005],now[1200005];
void add(int x,int y,int z)
{
    cnt++;
    e[cnt].to=y;
    e[cnt].val=z;
    e[cnt].next=head[x];
    now[x]=head[x]=cnt;
}

int h[120005],w[120005],cnth[120005];
queue<int>list[120005];
queue<int>qq;
int book[120005];
void bfs()
{
    memset(h,0x3f,sizeof(h));
    h[t]=0;
    qq.push(t);
    while(!qq.empty())
    {
        int u=qq.front();
        qq.pop();
        for(int i=head[u];i;i=e[i].next)
        {
            int v=e[i].to;
            if(e[i^1].val&&h[u]+1<h[v])
            {
                h[v]=h[u]+1;
                if(book[v]==0)
                {
                    qq.push(v);
                    book[v]=1;
                }
            }
        }
    }
}

void push_relabel(int x,int Delta)
{
    int i;
    int found=0;
    for(i=now[x];i;i=e[i].next)
    {
        int v=e[i].to;
        if(h[x]==h[v]+1&&e[i].val)
        {
            found=1;
            break;
        }
        else
        {
            now[x]=e[i].next;
        }
    }
    if(found==1)
    {
        int v=e[i].to;
        int flow=min(e[i].val,w[x]);
        if(v!=t)flow=min(flow,Delta-w[v]);
        e[i].val-=flow;
        e[i^1].val+=flow;
        w[x]-=flow;
        w[v]+=flow;
        if(w[x]<=Delta/2)list[level].pop();
        if(w[v]>Delta/2&&v!=s&&v!=t)
        {
            list[h[v]].push(v);
            level--;
        }
    }
    else
    {
        res++;
        list[level].pop();
        cnth[level]--;
        if(cnth[level]==0)
        {
            for(int i=1;i<=n;i++)
            {
                if(h[i]>level&&i!=s&&i!=t&&h[i]<n+1)
                {
                    h[i]=n+1;
                }
            }
        }
        int minn=inf;
        for(int i=head[x];i;i=e[i].next)
        {
            int to=e[i].to;
            if(e[i].val)
            {
                minn=min(minn,h[to]);
            }
        }
        h[x]=minn+1;
        if(minn==inf)h[x]=n+1;
        if(h[x]<=1200000)cnth[h[x]]++;
        if(h[x]<=1200000)list[h[x]].push(x);
        now[x]=head[x];
    }
}

int Excess_Scaling()
{
    bfs();
    if(h[s]==0x3f3f3f)return 0;
    h[s]=n;
    for(int i=1;i<=n;i++)
    {
        if(h[i]<0x3f3f3f)
        {
            cnth[h[i]]++;
        }
    }
    for(int i=head[s];i;i=e[i].next)
    {
        int to=e[i].to,val=e[i].val;
        if(val)
        {
            e[i].val-=val;
            e[i^1].val+=val;
            w[s]-=val;
            w[to]+=val;
        }
    }
    for(int k=1;k<=K;k++)
    {
        int Delta=(long long)1<<(K-k);
        for(int i=1;i<=n;i++)
        {
            if(w[i]>Delta/2)
            {
                if(h[i]<=1200000)list[h[i]].push(i);
            }
        }
        level=1;
        while(level<=n)
        {
            if(list[level].empty())level++;
            else
            {
                int x=list[level].front();
                push_relabel(x,Delta);
            }
        }
        for(int i=1;i<=32;i++)
        {
            while(!list[i].empty())
            {
                list[i].pop();
            }
        }
    }
    return w[t];
}

signed main()
{
    int i,j,k;
    ios::sync_with_stdio(0);
    cin>>n>>m>>s>>t;
    for(i=1;i<=m;i++)
    {
        int x,y,z;
        cin>>x>>y>>z;
        if(x==y)continue;
        add(x,y,z);
        add(y,x,0);
    }
    K=31;
    cout<<Excess_Scaling();
    return 0;
}

4-4 其它相关

文末作者提到使用另一种超额缩放算法可优化至 O(nm+n^2\sqrt{\log U})，但这不关键。

而使用动态树可以优化至 O(nm\log(n\sqrt{\log U}/m+2))。

[5]Improved time bounds for the maximum flow problem. 1989

最后，使用一种充分地大缩放可优化至 O(nm+n^2\frac{\log U}{\log\log U})，这也就是下文的内容了。

4-5.Large Medium Excess Scaling

[6]A FAST MAX FLOW ALGORITHM 2019

大中超额缩放算法？

基于超额缩放算法，在 LMES 中我们取一个 2 的幂 k，每轮 Scaling 当没有超额流大于 \Delta /k 时结束并将 \Delta 除以 k，总共缩放 O(\log_kU) 次。

对于每次 Scaling 我们将超额流较大即 e_i\geq \Delta/2 的点名为大溢出点存起来，将 \Delta/k\leq e_i<\Delta/2 的点称为中溢出点存起来。

若存在大溢出点取标号最小的推流，否则取标号最大的中溢出点推流，推不了流就重贴，没法推就结束这轮。

该算法的时间复杂度为 O(nm+n^2(k+\log_k U))， Ahuja 等人为了平衡时间，选择 k 为超过 2+\frac{\log U}{\log\log U} 的最小 2 的幂，其时间复杂度最终为 O(nm+n^2\frac{\log U}{\log\log U})。

4-6 随机化预流推进

[9]A Randmized Maximum-Flow Algorithm 1989

Cheriyan 和 Hagerup 提出对每个点的边随机一个编号每次取编号最小的推流等其它随机化的优化方法再基于 LMES 可以实现 O(nm+n^2\log^2n) 的最大流，后续也有一系列的随机化预流推进算法的进一步研究。

具体有空再说吧。

5. Enhanced LMES

我们考虑将弱多项式复杂度 O(nm+n^2(k+\log_k U)) 优化为强多项式 O(n^2k+nmlog_k n)。

对于 O(n^2k) 这项是由重贴标号导致的，在增强型 LMES 中保留，故不深入考虑。

相反我们考虑 O(n^2\log_k U) 这项，它出现在两个部分的复杂度限制中。

我们有两个部分的复杂度限制（详见论文第9页），对这两个地方的改进将使大中溢出点的推流数量降低至 O(kn^2+nm\log_k n)。

参考2012年Orlin的做法即上方 O(nm) 算法，我们的措施有：

• 将边划分为大中小三个分类。
• 对超大边进行缩点
• 特殊点
• 还流森林

接下来我们从这四个方面详细分析沿着这些思路是如何创造出一个强多项式的时间复杂度。

5-0 参数

我们先定义一些参数：Q,\epsilon ,M

• \frac{1}{4nk}\leq \epsilon\leq \frac{1}{n}
• 最小的整数 Q 满足 k^{-Q}\leq \frac{1}{4n}
• M=\epsilon^{-2},16n^2\leq M\leq 16k^2M^2
• Q$ 次缩放后 $\Delta$ 变为 $\epsilon \Delta

5-1 小，中，大边

我们设边 (i,j) 双向流量之和为 val_{(i,j)}。

若 val_{(i,j)}<\epsilon^5\Delta ，则 (i,j) 为小边。

若 \epsilon^5\Delta\leq val_{(i,j)}<2M\Delta ，则 (i,j) 为中边。

若 2M\Delta\leq val_{(i,j)} ，则 (i,j) 为大边。

考虑中等边，每个中边在 O(\log_k n) 个缩放阶段都是中等边，因此在整个缩放过程中,与中等边相连的节点会有 O(nm\log_k n) 次大，中超额流推送。

在考虑不与中等边相连的节点前，我们考虑超大边缩点。

5-2 超大边缩点

当一次 Scaling 开始前我们认为一条边可用流量大于 M\Delta 则其为超大边，超大边会一直都是超大边。

我们可以通过缩点超大边环来将问题变得更简单。

当发现超大边环时我们将它合并为一个节点并继续进行增强型 LMES 最后反向扩展合并节点。

但该算法因为存在没有推流且还流森林为空的无用缩放阶段而无法达到强多项式复杂度，我们要通过修改下次缩放大小来达到 O(nm\log_k n) 的复杂度。

5-3 特殊点

特殊点在 Scaling 开始前 e_i<1.5\epsilon^4 ，当最小的中超额流点的标号减小时 e_i<1.5\epsilon^4，由于有 m\log_k n 轮缩放，特殊点有 O(nm\log_{k} n) 次大超额流和中超额流的推送。

非特殊点，则会在 O(\log_k n)个缩放阶段内与中边相交或者成为被合并的超大边的节点。

这意味着在所有缩放阶段中，非特殊节点的总数为 O(m\log_k n)。因此，非特殊节点将也有 O(nm\log_k n) 个大超额流和中超额流推送的数量。

综上大超额流和中超额流推送的数量为 O(nm\log_k n)，但实际上为确保每个特殊点满足以上的性质，我们需要允许特殊点的超额流为负数，前提是 e_i\ge -1.5\epsilon^4。

这种差异在分析中产生了巨大的影响。通用预流推进算法通过确保每个超量在终止时都是非负的而终止。我们需要改进我们的算法，以使每个超额流在终止时都为 0。这需要新的方法将流推送到超额流为负数的节点。

我们将描述一种新的数据结构，该数据结构将流返回到超量为负的节点。使用此数据结构，超额流始终至少为 -\Delta，并且在终止时所有超额流都为非负数。

但是首先，我们定义一个节点成为特殊节点的含义。我们还解释了为什么允许特殊节点的超量为负值是有意义的。

5-4 还流森林

如果一个点不是特殊点且 e_i<\epsilon \Delta 我们认为他是不正常的，当一轮 Scaling 开始时，我们将不正常的点放入还流森林。

还流森林具有以下要求：

• 没环

• 每个连通分量都有根

• 根到每个点都有有向的超大边路径。

• 森林边 (i,j) 是超大边且 h_j\leq h_i+1 。

• 若某一轮中不正常点 v 被插入森林钟有欠流 \Delta \times\epsilon^2/k 。设 \Delta' = \epsilon ^ 2 \times \Delta，若 \Delta' 里 v 还是反常点，那我们把 \epsilon ^ 2 / k 单位的流从 v 的根寄给 v。

• 每个根节点都是一个正常的点。每个根 i 维护一个额外的变量 Reserve_i 表示所有后代的欠流和，故能够推出 e_i\ge\epsilon\times\Delta+Reserve_i

• 每个叶节点都是不正常的点，正常点的叶节点要被删掉

• 某轮中若 v 在还流森林中，则e_v\ge-\epsilon\times\Delta/2k

• 还流森林中 v 是一个不正常的点且 v 的欠流有 \Delta/k 个单位，那么 v 的根要把 \Delta/k 单位的流送给 v，送流完成后e_v\ge\Delta/2k

5-5 伪代码

5-6 时间复杂度