这篇文章所针对的是熟悉网络流理论的基本知识的读者。如果你对网络流理论的基本知识不是很了解的话,我会建议你先看参考文献[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)的流,那么我们就能得到最大流问题优化的模型:
如下图:
![最大流问题:增广路径算法的比较_网络](http://my.csdn.net/uploads/201205/25/1337912983_1321.jpg)
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。最大流的值从上文中得知其是有界的,根据假设三。而这种推理表明了该方法的有限性。
有了以上的这些准备,我们开始讨论算法。
第三节 最短增广路径算法,O(n*n*m)
Edmonds and Karp在1972年,以及Dinic在1970年都独立的证明了如果每步增广路径都是最短的话,那么整个算法将会执行O(n*m)步。之所以能实现这个最短路径算法(每条边的长度等于一)是利用了广度优先搜索算法BFS的,参考文献[2]、参考文献[6].增广路径算法被广泛的讨论和研究在许多书籍和文章中,包括参考文献[5]。我们回顾一下此算法:
如下图:
![最大流问题:增广路径算法的比较_网络_02](http://my.csdn.net/uploads/201205/25/1337926168_2449.jpg)
该算法通过执行O(n*m)步以找出一条增广路径。由于在广度优先搜索时最坏情况下需O(m)次操作,所有算法的总复杂度应是O(n*m*m)。我会在下面举个简单的例子。
第四节 改进的最短增广路径算法,O(n*n*m)
如我们早先提到的,找到任意最短增广路径的方法就是在残留网络中,通过执行广度优先搜索来找到这些路径。在最坏情况下BFS需要O(m)次操作以及规定最大流的时间复杂度为O(n*n*m)。于是在1987年,Ahuja和Orlin改进了最短增广路径算法,参考文献[1]。他们利用这一事实:在所有增广中,从顶点i到汇点t的最小距离是单调递增的,并且将每次增广的平均时间减少到O(n)。改进后的增广路径算法,运行的时间仍然是O(n*n*m)。现在我们就可以根据参考文献[1]来讨论它了。
定义一:距离函数d
残留容量Rij表示的是一个从节点集到非负整数的函数。如果距离函数满足一下的几个条件,那么我们就说它是有效的。
● d(t) = 0;
● d(i) <= d(j) + 1,且Rij > 0
很容易证明,在残留网络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。如下图所示:
![最大流问题:增广路径算法的比较_distance_03](http://my.csdn.net/uploads/201205/25/1337926848_9095.jpg)
该算法保留部分可容许的路径。比如,从源点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开发了一些基于运行时间的例子。利用他的想法,我们组成一个较为简单的网,这个网络不依赖下一条的选择。
如下图:
![最大流问题:增广路径算法的比较_网络_04](http://my.csdn.net/uploads/201205/25/1337927967_6202.jpg)
首先,最短增广路径算法沿着路径(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倍。
如下图所示:
![最大流问题:增广路径算法的比较_算法_05](http://my.csdn.net/uploads/201205/25/1337928493_1405.jpg)
第六节 最大容量路径算法,O(n*n*mlognU) / O(m*m lognU logn) / O(m*m lognU logU)
1972年,Edmonds and Karp发明了另一种找到增广路径的方法。在每一步,他们试图用尽可能最大的数来增加这个流。这个算法的另一个叫法是:Ford-Fulkerson方法梯度修正。这个修正算法代替了用BFS寻找最短路径,而改为利用Dijkstra算法来建立最大可能容量的路径。在增广之后,算法会在残留网络中找到另一条这样的路径,并沿着这条路径增加流,一直重复这几步直到找到最大流。
![最大流问题:增广路径算法的比较_distance_06](http://my.csdn.net/uploads/201205/25/1337928744_4942.jpg)
我们可以根据某一引理而得到算法的运行时间限制。为了理解这一证明,我们应该记住在网络中,任一流的值小于或者等于割的容量,或者阅读参考文献[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的边相交。所以,很显然有:
定理三:最大容量路径算法执行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)。不过这种方法在实际中的表现去不怎么样。
第七节 容量调整算法,O(m*mlogU)
1985年,Gabow描述了所谓的“位缩放”算法,由于Ahuja and Orlin在本节中描述的是类似容量调整算法。
非正式的,该算法的主要思想是增加沿路径有足够大容量的流,而不是沿着最大容量增加。正式的,我们引入一个参数Δ。首先,Δ是个很大的数,例如,令Δ = U。此算法试图找出一条增广路径,且其容量不小于Δ,当在残留网络中存在这样的Δ-路径时,那么沿着这条路径增加流,并重复此过程。
该算法可建立一个最大流或者令Δ/2,并且用新的Δ继续寻找路径和增加流量。沿着路径增加流(容量至少是Δ)的阶段被称为“Δ缩进阶段”或者“Δ阶段”。Δ是一个整数值,算法将会执行O(logU)次“Δ阶段”。当Δ等于1的时候,容量调整算法和Edmonds-Karp算法将没有任何区别。
![最大流问题:增广路径算法的比较_算法_07](http://my.csdn.net/uploads/201205/25/1337929800_4537.jpg)
引理四:对于每个“Δ-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(如果少于这些,算法运行的太快)。所有的运行时间都是以毫秒为单位。
![最大流问题:增广路径算法的比较_network_08](http://my.csdn.net/uploads/201205/25/1337930233_1263.jpg)
现在,我们来看看第二组测试实例。总共做了184次测试,所有网络的顶点都限制在400个。
![最大流问题:增广路径算法的比较_algorithm_09](http://my.csdn.net/uploads/201205/25/1337930374_7156.jpg)
我们很有兴趣的想知道这些算法在密集网络中是怎样运行的。我们来看看第三组测试:有200个密集网络,且其顶点限制在400个。
![最大流问题:增广路径算法的比较_distance_10](http://my.csdn.net/uploads/201205/25/1337930501_5327.jpg)
毫无疑问,经过改进后实现的Edmonds-Karp算法赢得了这场游戏。第二名则是被改进的调整算法拿下,使用BFS的调整容量算法拿下了第三名。
至于最大容量路径,最好使用一种用堆实现的变种方法;在稀疏网络中,它能收到很好的效果。而对于其他算法,它们只是适用于理论研究和兴趣爱好。
正如你看到的,复杂度为O(n*mlogU)的算法并不是那么快的,它甚至比复杂度为O(n*n*m)的算法还要慢。而我们最常用的却是复杂度为O(n*m*m)的算法,虽然此算法有更糟糕的时间范围,但是它的运行速度比一般算法都要快。
我的建议:始终使用BFS的容量调整路径算法,因为它很容易实现。改进的最短增广路径算法也是很相当容易实现的,但是你必须要非常小心,正确编写程序。在比赛中,它是很容易错过的一个bug。
在结束本文之前,我给出了改进的虽短增广路径算法的完整实现。我用邻接矩阵表示这个网络,这样能更好的理解算法。在实际分析中我们用的是不一样的实现,邻接矩阵比邻接表实现起来相对慢一些。不过,最终还是由读者选择最适合自己的数据结构。
/******improved shortest augmenting path algorithm******/ #include <stdio.h> #define N 2007 // Number of nodes #define oo 1000000000 // Infinity // Nodes, Arcs, the source node and the sink node int n, m, source, sink; // Matrixes for maintaining // Graph and Flow int G[N][N], F[N][N]; int pi[N]; // predecessor list int CurrentNode[N]; // Current edge for each node int queue[N]; // Queue for reverse BFS int d[N]; // Distance function int numbs[N]; // numbs[k] is the number of nodes i with d[i]==k // Reverse breadth-first search // to establish distance function d int rev_BFS() { int i, j, head(0), tail(0); // Initially, all d[i]=n for(i = 1; i <= n; i++) numbs[ d[i] = n ] ++; // Start from the sink numbs[n]--; d[sink] = 0; numbs[0]++; queue[ ++tail ] = sink; // While queue is not empty while( head != tail ) { i = queue[++head]; // Get the next node // 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]]++; } } return 0; } // Augmenting the flow using predecessor list pi[] int Augment() { int i, j, tmp, width(oo); // Find the capacity of the path for(i = sink, j = pi[i]; i != source; i = j, j = pi[j]) { tmp = G[j][i]; if(tmp < width) width = tmp; } // Augmentation itself for(i = sink, j = pi[i]; i != source; i = j, j = pi[j]) { G[j][i] -= width; F[j][i] += width; G[i][j] += width; F[i][j] -= width; } return width; } // Relabel and backtrack int Retreat(int &i) { int tmp; int j, mind(n-1); // Check all adjacent edges // to find nearest for(j=1; j <= n; j++) // If there is an arc // and j is "nearer" if(G[i][j] > 0 && d[j] < mind) mind = d[j]; tmp = d[i]; // Save previous distance // Relabel procedure itself numbs[d[i]]--; d[i] = 1 + mind; numbs[d[i]]++; // Backtrack, if possible (i is not a local variable! ) if( i != source ) i = pi[i]; // If numbs[ tmp ] is zero, algorithm will stop return numbs[ tmp ]; } // Main procedure int find_max_flow() { int flow(0), i, j; rev_BFS(); // Establish exact distance function // For each node current arc is the first arc for(i=1; i<=n; i++) CurrentNode[i] = 1; // Begin searching from the source i = source; // The main cycle (while the source is not "far" from the sink) for( ; d[source] < n ; ) { // 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; } } // End of the main cycle // We return flow value return flow; } // The main function // Graph is represented in input as triples <from, to, capacity> // No comments here int main() { int i, p, q, r; scanf("%d %d %d %d", &n, &m, &source, &sink); for(i = 0; i < m; i++) { scanf("%d %d %d", &p, &q, &r); G[p][q] += r; } printf("%d", find_max_flow()); return 0; }参考文献:
[1] Ravindra K. Ahuja, Thomas L. Magnanti, and James B. Orlin. Network Flows: Theory, Algorithms, and Applications.
[2] Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest. Introduction to Algorithms.
[3] Ford, L. R., and D. R. Fulkerson. Maximal flow through a network.
[4] Norman Zadeh. Theoretical Efficiency of the Edmonds-Karp Algorithm for Computing Maximal Flows.
[5] _efer_. Algorithm Tutorial: MaximumFlow.
[6] gladius. Algorithm Tutorial: Introduction to graphs and their data structures: Section 1.
[7] gladius. Algorithm Tutorial: Introduction to graphs and their data structures: Section 3.
[8] http://elib.zib.de/pub/mp-testdata/generators/index.html -- A number of generators for network flow problems.