Edmond-Karp 算法 网络流 增广路

比如说这个，EK算法首先由俄罗斯科学家Dinic在1970年提出，没错，就是dinic算法的创始人，实际上他提出的也正是dinic算法，在EK的基础上加入了层次优化，这个我们以后再说，1972年Jack Edmonds和Richard Karp发表了没有层次优化的EK算法。但实际上他们是比1790年更早的时候就独立弄出来了。



你看，研究一下历史也是很有趣的。



扯远了，首先来看一下基本的网络流最大流模型。



有n个点，有m条有向边，有一个点很特殊，只出不进，叫做源点，通常规定为1号点。另一个点也很特殊，只进不出，叫做汇点，通常规定为n号点。每条有向边上有两个量，容量和流量，从i到j的容量通常用c[I,j]表示,流量则通常是f[I,j]。通常可以把这些边想象成道路，流量就是这条道路的车流量，容量就是道路可承受的最大的车流量。很显然的，流量<=容量。而对于每个不是源点和汇点的点来说，可以类比的想象成没有存储功能的货物的中转站，所有”进入”他们的流量和等于所有从他本身”出去”的流量。



把源点比作工厂的话，问题就是求从工厂最大可以发出多少货物，是不至于超过道路的容量限制，也就是，最大流。



比如这个图。每条边旁边的数字表示它的容量。





下面我们来考虑如何求最大流。



首先，假如所有边上的流量都没有超过容量(不大于容量)，那么就把这一组流量，或者说，这个流，称为一个可行流。一个最简单的例子就是，零流，即所有的流量都是0的流。

这样我们就得到了一个更大的流，他的流量是之前的流量+delta，而这条路就叫做增广路。



我们不断地从起点开始寻找增广路，每次都对其进行增广，直到源点和汇点不连通，也就是找不到增广路为止。当找不到增广路的时候，当前的流量就是最大流，这个结论非常重要。

   这篇文章我们将重温最大流问题，实现一些最有名的增广路径算法的实际分析的目标。我们将讨论的这几种算法的复杂度在O(n*m*m)到O(n*mlogU)之间，并且从讨论的结果中得到在实践中最有效的一种。正如我们所想的，理论上的复杂度并不能揭示该算法在实际中的价值。



   这篇文章所针对的是熟悉网络流理论的基本知识的读者。如果你对网络流理论的基本知识不是很了解的话，我会建议你先看参考文献[1]、参考文献[2]以及参考文献[5]---算法教程之最大流问题。



   在第一节，将涉及到最大流理论的一些必要的定义和声明。在中间的章节，我们将会着重讨论增光路径算法。在最后一节，将会展示时间分析的结果，并强调在实践中的最佳算法，同时我也会给出该算法的简单实现。

   假设一有向网G = (V, E)，其中V表示顶点集，E表示边集。相关联的两个节点i和j所组成的弧arc(i,j)均有非负的容量Uij。同时我们在有向网G中定义了两个特殊的顶点，即一个源点s和一个汇点t。

   对于V集中的节点i，我们用E(i)表示从节点i产生的所有边。

   令U = max Uij。

   令n表示顶点的个数，m表示边的个数。



   我们期望从源点s到汇点t之间找到这个最大流，并且在所有节点中满足:从一个顶点到另一个顶点的流不能超过设定的容量。用Xij代表边arc(i,j)的流，那么我们就能得到最大流问题优化的模型： 如下图：

 



   Xij被称为一个可行的解决方案或者可行流，并且它满足所有的约束条件。给一个流x，我们能根据以下的想法来构造出残留网络。假设边(i,j)是流中的单元Xij，那么我们定义边(i,j)的残留容量Rij = Uij - Xij。这就意味着我们可以从顶点i到顶点j压入额外的单位流量Rij。如果我们从j到i的弧(i,j)压入Xij单位流量，就可以抵消从i到j的流Xij。



   因此，给定一个可行流x,我们定义流的残留网络x如下：假设一网络G = (V, E)，一个可行的解决方案x可产生一个新的残留网络，我们用Gx = (V, Ex)来定义这个残留网络，其中，Ex是一个可行的解决方案对应的残留边的集合x。 附注：残留的容量 + 反向平衡的流量共同构成了残留网络。



   那什么是Ex呢？我们用边(i,j)、边(j,i)来代表弧(i,j)：边(i,j)的残留容量Rij = Uij - Xij，并且边(j,i)的残留容量Rij = Xij。然后我们就能从一个正的残留容量的新边集中构造集合Ex。

   在这节中，我将描述一种构造所有增广路径算法的方法，这种方法是由Ford and Fulkerson在1956年发明的。

   增广路径是找出在残留网络中从源点到汇点的有向路径。增广路径的残留容量是路径中任意边所形成的最小残留容量。显然，我们可以沿着增广路径从源点到汇点发送额外的流。



   假如有这么一条路，这条路从源点开始一直一段一段的连到了汇点，并且，这条路上的每一段都满足流量 < 容量。那么，我们一定能找到这条路上的每一段的(容量-流量)的值当中的最小值delta。我们把这条路上每一段的流量都加上这个delta，一定可以保证这个流依然是可行流。这样我们就得到了一个更大的流，他的流量是之前的流量+delta，而这条路就叫做增广路径。



   所有的增广路径算法的构造是基于增广路径定理的：



   定理一(增广路径定理)：流x是最大流当且仅当这个残留网络不包含其他增广路经。



   由这个定理我们得到一种找到最大流的方法。这种方法通过在所有路径中不断地找出增广路径和增广流，直到网络中不在包含这样的路径。我们要讨论的一些算法，它们所不不同的只是寻找增广路径的方法。



   我们认为最大流问题基于以下假设：



          假设一：这个流网络是一个有向网。



          假设二：网络中的所有容量都是非负整数。

          附注：这个假设对于某些算法不是必须的，这些算法的复杂边界涉及到数据的完整性。



          假设三：这个问题有一个最佳解决方案，且这个方案是有界的。

          附注：这个特定的假设意味着从源点到汇点是有容量限制的路径。



          假设四：这个网中不包含平行的弧。

          附注：这个假设的规定不失一般性，因为我们可以总结出所有平行弧的容量。



   至于这些假设为什么是正确的，我将证明留给读者自己。



   其实我们很容易确定上述办法的正确性。根据假设二，对于每一增广步骤，我们增加流值至少一个单位流量，通常开始的流值是0。最大流的值从上文中得知其是有界的，根据假设三。而这种推理表明了该方法的有限性。



   有了以上的这些准备，我们开始讨论算法。

   Edmonds and Karp在1972年，以及Dinic在1970年都独立的证明了如果每步增广路径都是最短的话，那么整个算法将会执行O(n*m)步。之所以能实现这个最短路径算法（每条边的长度等于一）是利用了广度优先搜索算法BFS的，参考文献[2]、参考文献[6].增广路径算法被广泛的讨论和研究在许多书籍和文章中，包括参考文献[5]。我们回顾一下此算法： 如下图：

   如我们早先提到的，找到任意最短增广路径的方法就是在残留网络中，通过执行广度优先搜索来找到这些路径。在最坏情况下BFS需要O(m)次操作以及规定最大流的时间复杂度为O(n*n*m)。于是在1987年，Ahuja和Orlin改进了最短增广路径算法，参考文献[1]。他们利用这一事实：在所有增广中，从顶点i到汇点t的最小距离是单调递增的，并且将每次增广的平均时间减少到O(n)。改进后的增广路径算法，运行的时间仍然是O(n*n*m)。现在我们就可以根据参考文献[1]来讨论它了。

   很容易证明，在残留网络Gx中，从某一顶点i到汇点t，节点i的有效距离标号（由d(i)表示）是最短路径长度的下界。若残量网络中任意一点的距离标号正好等于该顶点至汇点t的最短路路长，则称距离函数是精确的。同时我们也很容易的证明，如果d(s) >= n，那么残留网络中不再有从源点到汇点的路径。



   如果满足d(i) = d(j) + 1的边，我们称这条边(i,j)是可容许的，反之，其他的边是不容许的。如果一条路径包含了从源点s到汇点t包含了可容许的边，那么这条路径是可容许的。显然，一条了容许的路径是从源点到汇点的最短路径。对于可容许路径中的每一条边需满足条件Rij > 0，它是一条增广路径。



   因此，改进后的最短增广路径算法包括四步：main cycle, advance, retreat and augment。如下图所示：

   该算法保留部分可容许的路径。比如，从源点s到某一顶点i包含了容许边。算法从部分可容许路径的末节点（这类节点也称当前节点）开始执行advance或者retreat步骤。如果从当前节点始发有一些可容许的边，那么将会执行算法的advance步骤，并会将这条边添加到部分可容许的边中。否则，算法将会执行retreat步骤。



   如果部分容许路径到达了汇点，我们就执行一次增广。当d(s) >= n是算法结束。另外，Ex(i)的正规表达式：Ex(i) = { (i,j) in E(i): Rij > 0 }。



   现在我大概证明一下该算法的运行时间O(n*n*m)。



   引理一：算法每步都会保存距离标号。此外，每次重标号都要严格地增加一个节点的距离标号。



   证明描述：对一些重标号操作和增广执行归纳法。



   引理二：每个节点的距离标号之多增添n次。连续的，重标号操作至多执行n*n次。



   证明：引理二是引理一的延伸，如果d(s) >= n，那么残留网络中不在包含增广路径。



   因为改进的最短增广路径算法产生增广是沿着最短路径（和没有改进的算法一样），所以增广的总数都是相同的O(n*m)。执行一次retreat步就重标一个节点，这就是为什么retreat steps需要O(n*n)（根据引理二）。执行retreat/relabel步骤的时间是O( n ∑i in V |E(i)| ) = O(nm)。由于一次增广需要时间O(n)，所以总的增广时间应是O(n*n*m)。advance步骤执行的总时间是增广时间加上retreat/relabe的时间，也是O(n*n*m)。于是我们得到以下定理：



   定理二：改进后的最短增广路径算法的运行时间为O(n*n*m)。



   Ahuja and Orlin认为这是对该算法非常实用的一次改进。因为当最大流被找到的时候，算法执行了许多无用的重标号操作，解决无效操作更好的办法就是添加一个终止的条件。我们引入一个(n+1)维的数组numbs，下标从0到n。numbs(k)代表的值是节点的个数，它的参数k等于距离标号。当算法利用BFS计算初始距离标号的同时，初始化这个数组numbs。



   当算法从节点x到节点y增加节点的距离标号时，将会从numbs(x)减1，而numbs(y)加1，同时检查numbs(x)是否等于0。如果等于0，算法终止。



   这种方法是一种启发式，但是它在实际中确实很好用。证明留给读者。（提示：当节点i有d(i) > x，以及节点j有d(j) < x时产生割，此时就要利用最大流最小割定理。）

   本节，我们将在最坏情况下来比较两种最短增广路径算法的运行时间。

   最坏情况下，不论是改进的还是没有改进的算法都会执行O(n*n*n)次。如果m = n*n，Norman Zade开发了一些基于运行时间的例子。利用他的想法，我们组成一个较为简单的网，这个网络不依赖下一条的选择。 如下图：



   除了源点s和汇点t之外，其他的节点被分为四个子集：S={s1,...,sk},T={t1,...,tk},U={u1,...,u2p},V={v1,...,v2p}。集合S和集合T包含k个节点，而集合U和集合V包含2p个节点。k和p都是定整数。上图中用粗体线连接的边（连接S和T）表示单位容量，用虚线连接的边表示无穷大容量，其他的边表示的容量为k。



   首先，最短增广路径算法沿着路径(s,S,T,t)增加流k*k次，此时流的长度等于3，这些路径的容量都是单位容量。之后，残留网络中将包含反向弧(T,S)，并且算法将会选择另外k*k个长度为7的增广路径(s,u1,u2,T,S,v2,v1)。接下来，算法会继续选择长度为11的增广路径(s,u1,u2,u3,u4,S,T,v4,v3,v2,v1,t)。如此这般，这般如此，一直执行下去。



   这时候，让我们来计算一下网络中的一些参数。顶点的个数n = 2*k + 4*p + 2，边的个数m = k*k + 2*p*k + 2*k +4*p。那么就很容易的得到增广的次数a = k*k*(p+1)。



   我们在最坏情况下做了五次测试，每次测试的顶点分别为：100个、148个、202个、250个、298个，并比较了改进后的算法和没改进的算法的运行时间。从下图中我们得知，改进的算法更快一些。对于有298个顶点的网络，改进的算法比没改进的算法快23倍。通过实践分析后我们得知：一般情况下，改进后的算法比没改进的算法快14倍。 如下图所示：



   然而，我们的比较结果并不是最可靠的，因为我们只是用了其中的一种网络。我们只想证明改进的算法比没改进的算法的运行速度快，并且快的数量级是线性的。我将在文章的末尾讲解一个更为准确的比较。

   1972年，Edmonds and Karp发明了另一种找到增广路径的方法。在每一步，他们试图用尽可能最大的数来增加这个流。这个算法的另一个叫法是：Ford-Fulkerson方法梯度修正。这个修正算法代替了用BFS寻找最短路径，而改为利用Dijkstra算法来建立最大可能容量的路径。在增广之后，算法会在残留网络中找到另一条这样的路径，并沿着这条路径增加流，一直重复这几步直到找到最大流。



   毫无疑问，算法在整数容量的条件下是正确的。然而，对于非整数边容量，经过测试，算法很有可能由于失败而终止。

   我们可以根据某一引理而得到算法的运行时间限制。为了理解这一证明，我们应该记住在网络中，任一流的值小于或者等于割的容量，或者阅读参考文献[1]，参考文献[2]。我们用c(S,T)来表示割(S,T)的容量。



   引理三：让F表示最大流的值，那么G包含了容量不小于F/m的增广路径。



   证明：假设G不包含这样的路径。我们构造一个集合E' = { (i,j) in E: Uij ≥ F/m }。令网络G' = (V, E')，且网络中没有从源点s到汇点t的路径。S是从G中和T = V \ S中获得的节点的集合。很明显，(S,T)是一个割并且有c(S,T) >= F。但是割(S,T)只和Uij < F/m的边相交。所以，很显然有： c(S,T) < (F/m)_m = F， 于是，这与事实c(S,T) ≥ F相矛盾。



   定理三：最大容量路径算法执行O(mlog(nU))次增广。



   证明：假设算法经过k次增广后终止。让f1表示第一次发现增广路径的容量，f2表示第二次，依此类推，fk表示第k次增广路径的容量。此时，令Fi = f1 + f2 +...+ fi，让F*表示最大流的值。根据定理三，就可以证明：

   为了找到路径的最大容量，我们用Dijkstra算法，该算法在每次迭代时会带来额外的开销。因为Dijkstras算法的简单实现的复杂度为O(n*n)，最大容量路径算法总的运行时间是O(n2mlog(nU))。



   对于稀疏网络，Dijkstra算法利用堆实现的运行时间是O(mlogn)，对于最大流则需O(m2 logn log(nU))。看起来这比改进后的Edmonds-Karp算法更好一些，然而，这个估计是极具欺骗性的。



   还有另一种的变种方法来找到最大容量路径，可以利用二分查找来建立最大容量路径。设找最大容量路径的区间为[0,U]，如果一些路径的容量等于U/2，那么我们继续在区间[U/2,U]上找这条路径；否则，我们将在区间[0，U/2-1]上找这条路径。这种方法需要额外的O(mlogU)开销，并给出了最大流算法的时间约束O(m*mlog(nU)logU)。不过这种方法在实际中的表现去不怎么样。

   1985年，Gabow描述了所谓的“位缩放”算法，由于Ahuja and Orlin在本节中描述的是类似容量调整算法。



   非正式的，该算法的主要思想是增加沿路径有足够大容量的流，而不是沿着最大容量增加。正式的，我们引入一个参数Δ。首先，Δ是个很大的数，例如，令Δ = U。此算法试图找出一条增广路径，且其容量不小于Δ，当在残留网络中存在这样的Δ-路径时，那么沿着这条路径增加流，并重复此过程。



   该算法可建立一个最大流或者令Δ/2，并且用新的Δ继续寻找路径和增加流量。沿着路径增加流（容量至少是Δ）的阶段被称为“Δ缩进阶段”或者“Δ阶段”。Δ是一个整数值，算法将会执行O(logU)次“Δ阶段”。当Δ等于1的时候，容量调整算法和Edmonds-Karp算法将没有任何区别。



   我们可以很容易得到一条容量至少是Δ的路径---在O(m)时间内（用BFS算法）。开始，我们令Δ的值可以是U或者是Δ的二次方但不能超过U。



   引理四：对于每个“Δ-phase”，算法的最坏情况是执行O(m)次增广。

   定理四：容量调整算法的运行时间是O(m2logU)。



   请记住，此时当寻找一条增广路径时，使用BFS和DFS是没有任何区别的。但是，在实践中却是截然相反的。

   在上一节，我们介绍了一种运行时间为O(m*mlogU)的寻找最大流的算法。本节我们将改进此算法，将其运行时间提高至O(n*mlogU)。



   现在我们独立的看看每个“Δ-phase”。回想上一节，每个“Δ-scaling phase”都包含了O(m)次增广。当描述最短增广路径算法的改进型时，我们会将相似的技术应用到“Δ-phase”中。在每个阶段，我们通过仅使用的路径（容量至少等于Δ）来寻找最大流。对改进的最短增广路径算法的复杂度分析意味着：如果算法保证执行O(m)次增广，那么它将运行O(nm)的时间内，这是因为增广的时间从O(n*n*m)减少到O(n*m)以及其他的一些操作，就像前面，需要O(n*m)的时间。这些原因立即对改进的容量调整算法的运行时间形成了O(nmlogU) 的约束。



   不幸的是，这种改进在实践中几乎对运行时间的降低起不了作用。

   现在，我们来做一些有意思的事情。在这节，我将会以实际应用的观点来比较前面所有介绍的算法。为了实现这一目标，在超链8的帮助下，我做了一些测试案例，并将它们借助密度分成三组。第一组测试的网络满足：m ≤ n1.4---一些稀疏的网络；第二组测试的网络--中等密度网络满足：n1.6 ≤ m ≤ n1.7；第三组测试的网络--几乎是完全图（包括完整的非循环网络）满足：m ≥ n1.85。



   我在前面已经讲过所有算法的一些简单实现。所有的实现都是用邻接表来表示网络。

   我们先来对第一组做些测试。有564个稀疏网络，且它们的顶点数都现在2000（如果少于这些，算法运行的太快）。所有的运行时间都是以毫秒为单位。



   从图表得知，在稀疏网络中，试图不用堆实现的Dijkstra的最大容量路径算法确实是一个严重的错误。因为用堆实现的运行速度确实比期望的要快。大约在同一时间执行容量调整算法（使用DFS和BFS），然而改进后的实现时间几乎是原来的两倍快。然而令人不解的是，在稀疏网络中，改进的最短路径算法被证明是最快的。



   现在，我们来看看第二组测试实例。总共做了184次测试，所有网络的顶点都限制在400个。

 

   在中等密度网络中，通过二分查找实现的最大容量路径算法留下了许多不足之处，但是用堆实现仍然比没有用堆实现的要快。用BFS实现容量调整算法要比用DFS实现快。改进后的调整算法和改进后的最短增广路径算法在这次测试中都是很优秀的。



   我们很有兴趣的想知道这些算法在密集网络中是怎样运行的。我们来看看第三组测试：有200个密集网络，且其顶点限制在400个。



   现在，我们看看容量调整算法的BFS和DFS的版本之间的差异。出乎预料的是，用堆实现的Dijkstra的最大容量路径算法被证明快于没有用堆实现的算法。



   毫无疑问，经过改进后实现的Edmonds-Karp算法赢得了这场游戏。第二名则是被改进的调整算法拿下，使用BFS的调整容量算法拿下了第三名。



   至于最大容量路径，最好使用一种用堆实现的变种方法；在稀疏网络中，它能收到很好的效果。而对于其他算法，它们只是适用于理论研究和兴趣爱好。



   正如你看到的，复杂度为O(n*mlogU)的算法并不是那么快的，它甚至比复杂度为O(n*n*m)的算法还要慢。而我们最常用的却是复杂度为O(n*m*m)的算法，虽然此算法有更糟糕的时间范围，但是它的运行速度比一般算法都要快。



   我的建议：始终使用BFS的容量调整路径算法，因为它很容易实现。改进的最短增广路径算法也是很相当容易实现的，但是你必须要非常小心，正确编写程序。在比赛中，它是很容易错过的一个bug。



   在结束本文之前，我给出了改进的虽短增广路径算法的完整实现。我用邻接矩阵表示这个网络，这样能更好的理解算法。在实际分析中我们用的是不一样的实现，邻接矩阵比邻接表实现起来相对慢一些。不过，最终还是由读者选择最适合自己的数据结构。 /******improved shortest augmenting path algorithm******/ #include <stdio.h>

// Check all adjacent nodes

for(j = 1; j <= n; j++)  {			

		

  // If it was reached before or there is no edge

  // then continue

  if(d[j] < n || G[j][i] == 0)  continue;

 

  // j is reached first time

  // put it into queue

  queue[ ++tail ] = j;

 

  // Update distance function

  numbs[n]--;

  d[j] = d[i] + 1;

  numbs[d[j]]++;

		

}   }

// Start searching an admissible arc from the current arc

for(j = CurrentNode[i]; j <= n; j++)

  // If the arc exists in the residual network

  // and if it is an admissible

  if( G[i][j] > 0 && d[i] == d[j] + 1 )

    // Then finish searhing

    break;

 

// If the admissible arc is found

if( j <= n )  {			

  CurrentNode[i] = j;  // Mark the arc as "current"

  pi[j] = i;  //  j is reachable from i			

  i = j;  // Go forward

 

  // If we found an augmenting path

  if( i == sink )  {

    flow += Augment();  // Augment the flow

    i = source;  // Begin from the source again

  }

}

// If no an admissible arc found

else  {

  CurrentNode[i] = 1;  // Current arc is the first arc again

 

  // If numbs[ d[i] ] == 0  then the flow is the maximal

  if( Retreat(i) == 0 )

    break;	

 

}