找到图中最小有向圈的一个时间复杂度为O(mn)的算法

找到图中最小有向圈的一个时间复杂度为O(mn)的算法
本文介绍一种在n个点 m条边并且没有负长度有向圈的有向图中找到最小有向圈的一个O(mn)时间复杂度的算法，这一结果是对于之前的最好的时间界O(mn+n^2 loglogn)的一个改进。
此算法首先在O(mn)时间复杂度内找到了具有最小平均长度λ^*的圈，因此它选择了点势函数来使得节约成本是λ^*或者更大。之后，该算法解决了所有点对之间的最短路问题，但是接下来将注意力集中于长度最多为nλ^*的最短路，将每一个源点的最短路时间的计算加速到O(m)，这样所有的点时间的计算就到了O(mn)。
同时，本文还进行了该种算法与现有的在有向图中找到最小有向圈的算法的实验并加以比较。实验表明，该算法在解决最小有向圈问题时比解决最短路问题更有优势。分析表明解决最小有向圈问题的运行时间比平均的O(n^2 )时间要快得多。
一：介绍
我们在有n个点 m条边并且没有负长度有向圈的有向图中来找最小有向圈(MLDC)。Floyd[8]与Warshall[19]在O(n^3 )时间复杂度内解决了该问题。一个可以选择的办法是找到所有点对之间的最短路。如果图中有负向弧，用修正标号算法来解决第一个最短路问题，因此，可以用节约成本来将此问题转化为非负长度的问题，之后的n-1个最短路问题可以用Dijkstra’s 算法来解决。用Fredman和Tarjan[10]算法来解决最短路问题，其运行时间是O(mn+n^2 logn)。Pettis[17]提供了一种在O(mn+n^2 loglogn)时间内针对于稀疏图的改善算法。基于Itai和Rodeh[13]的结果，Roditty和Williams[18]用快速矩阵乘积，对于有向图的权重为[-M,M]之间的整数时，引入了一种在O(mn^2.376 )时间复杂度内解决MLDC问题的算法，但是不是很实用。
现在阐述O(mn)解决MLDC问题的算法。该算法改进了其他的组合算法，并且是最好的强多项式时间界，其运行时间达到了Bellman[2]、Ford[9]和Moore[15]算法中探测有向负圈的时间界。该算法最先用Karp[14]的O(mn)算法来计算最小平均圈(MMC)问题。在用点势来修正花费之后，其实现了n个单源点的最短路问题的计算问题，其中长度最多为nλ^*，λ^*为圈的最小平均长度。n个点，每一个单源点的最短路的计算用Dijkstra’s 算法执行中的Dial[6]的变形可在O(m)时间内解决，因此，可以得到MLDC问题的O(mn)时间复杂度的算法。
有一个启发性的证据证明比O(mn)更快地解决MLDC问题是困难的，假设存在MLDC的最优圈W，要想证明W的最优性，需要证明任一点j∈W，包含j的最小圈的长度至少为C(W)，这需要Ω(n)个计算，每一个有需要Ω(m)步，所以共需要Ω(mn)步。对于MLDC的这一情形与MMC问题做一个比较。当前对于MMC最好的强多项式时间复杂度是O(mn)，然而，可以快速证明W是MMC问题的最优解。λ^*为圈W的平均长度，p为引理1中的点势函数，则W的每条边的节约成本是λ^*，图G的其他边的节约成本至少为λ^*。这在O(m)的时间内说明了W 的最优性，同时也为在比O(mn)更快地时间内解决MMC问题提供了可能。
关于MMC算法的研究已经由Dasdan et al. [4]，Dasdan[3]，和Georgiadis et al. [11]进行了研究。这篇文章主要是对于MLDC算法的实验性研究，截止到目前，这是唯一一篇对MLDC算法进行的实验性研究。实验表明，该算法的平均情形的运行情况都比最坏情况下的O(mn)时间复杂度快得多。
在本文中，我们依赖于Howard[12]算法来找到MMC。Howard 算法在实际中更加有效。在用点势来修正花费之后，其实现了n个单源点的最短路问题的计算问题，其中长度最多为nλ^*，λ^*为圈的最小平均长度。n个点，每一个单源点的最短路的计算用Dijkstra’s 算法执行中的Dial[6]的变形可在O(m)时间内解决，因此，可以得到MLDC问题的O(mn)时间复杂度的算法。对于我们的大多数实验，先解决MMC问题确实能够加速整体的运行时间，但是解决MMC问题的运行时间是一大瓶颈，其运行时间比解决所有的最短路问题时间还要长。
二、关于最小有向圈的初期结果。
给出有向图G=(V,A)，V={1、2⋯n}，A为m条弧，对于∀i∈V，Γ_i^+={j∈V│(i,j)∈A}即为i点的后继点， Γ_i^-={j∈V│(j,i)∈A}即为i点的前继点，对于∀(i,j)∈A，C_ij∈R表示其长度，有时将其称为弧(i,j)的费用。从点s到点t的有向路P_st 是序列〈i_1,(i_1,i_2 ),i_2⋯(i_(l-1),i_l ),i_l 〉，满足i_1=s，i_l=t，∀1≤ω≤l-1，(i_(ω-1),i_ω )∈A，有向路的长度是所有弧长之和，点i_1 到i_r 的有向路加上弧(i_r,i_1 )构成有向圈，C(W)=∑_((i,j)∈W)▒C_ij 表示有向圈W的长度，MLDC即为判别最小长度的有向圈问题，若G不含负长度且至少含有一个有向圈则称为MLDC问题可行。对于每个圈W，λ(W)=C(W)/|W| 表示圈W的平均长度，MMC表示确定最小平均长度的有向圈，W^*表示MMC的最优解，λ^*=λ(W^* )，有向图是强连通的若是图中任何两点都可以相互到达。图G的强连通分支是G的最大的强连通子图。
引理1众所周知，例如可从Ahuja et al.[1]知道。
引理1：W^*表示强连通图G=(V,A)的MMC的最优解，λ^*=λ(W^* )，通过令C_ij^'=C_ij-λ^*，(i,j)∈A，由G得到G^' ，P(i)是G^'中从点1到点i的最短路的长度，∀(i,j)∈A，节约成本为C_ij^P=C_ij+P(i)-P(j)=C_ij^'+P(i)-P(j)+λ^*，则有：(1)：对G中的任何有向圈W有C(W)=C^P (W)。(2)：∀(i,j)∈A，C_ij^P≥λ^*，∀(i,j)∈W，C_ij^P=C_ij+P(i)-P(j)=C_ij+λ^* 〖-C_ij=λ〗^*。
可用Karp[14]的算法在O(mn)时间复杂度内确定λ^*和P，由引理1知道，MLDC带有花费C⟺在G中带有花费为C^P的MLDC问题。对于每一对顶点对i，j，计算从点i到点j且经过弧C^P的最短路的长度d^P (i,j)。最小有向圈可以由min{d^P (i,j)+C_ji^P:(j,i)∈A}来找到。通过解决所有点对的最短路问题来解决MLDC问题是众所周知的。因为d^P (i,j)≤C^P (W^* )=C(W^* )=λ^* |W^* |，可以通过求d^P (i,j)来加速计算P的过程。
三、在O(mn)时间复杂度内确定n条最短路来解决MLDC问题。
假设C=C^P，d=d^P，考虑图G，其中每条弧长C_ij≥λ^*，用s点来表示Dijkstra’s 算法的源点，则从s点到j点的最短路d^P (s,j)<(|W^* |-1) λ^*，则每一个包含弧(k,s)的圈的长度至少为d(s,k)+C_ks，从而有C(W^* )=λ^* |W^* |。
我们采用Dijkstra’s [7]算法的Dial[6]变式，该算法至多占用|W^* |-1个储存桶，采取由Denardo 和Fox[5]发展的Dial算法来加速计算。假设在某些迭代中，已经得到从s点出发的最短路。下面描述接下来的实现过程：d_s (j)表示点j的当前距离编号，所以d_s (j)是在G中从s到点j的某些路长，所以有d_s (j)≥d(s,j)，d(s,j)表示点s到点j的最短路的路长。pred(j)表示点s到点j，长度为d_s (j)的路上点j的前继点。设初值d_s (s)=0，d_s (j)=∞(j≠s)。该算法保持一列|W^* |-1个储存桶，其中储存桶k包含使得〖kλ〗^*≤d_s (j)<(k+1)λ^*的任意点j。d_s (s)=0所以点s不在任何一个储存桶中。当从点s出发找到最短路时，将点j嵌入到适当地储存桶中，这一过程成为点对(s,j)的储存桶的更新，如下程序所示：
Procedure bucket-update(s,j) 点对(s,j)的储存桶的更新
Dial算法迭代地选择点i，第一个被选择的点是点s，随后，该算法从最小的非空储存桶中任意选择一点i，然后从这个桶里面删除点i 。并不需要选择具有最小距离编号的点i。若所有的储存桶都是空的，则停止寻找包含点s 的最小圈。选择了顶点i之后，算法判断从点s到点i 的路加弧(i,s)构成的圈是否为当前最小圈的改善，随后检查每一条从点i出发的弧，对于每一条被浏览的弧(i,s)，令d_s (j)=min{d_s (j)，d_s (i)+C_ij }，并适当地更新pred(j)。若d_s (j)值下降，且d_s (j)<λ^* |W^* |，则该算法检查点j是否在适当的储存桶，若不在，则将其挪到适当地储存桶中。
更新距离标号的程序如下图2所示。图3是寻找包含点s的最小有向圈的算法。在这个算法中，mldc表示目前找到的最小有向圈的长度。图4是找最小有向圈的全部算法。
在Dial最短路算法中的通常实施中，从储存桶k中选择一个点i时，就选择了一个最小的距离标号。在包含点s的最小有向圈中，从储存桶k中任意选择一个点都是可以的，因为在Denardo 和Fox[5]中有如下结果：
引理2：若点j是非空储存桶里面的任一点，则d_s (j)=d(s,j)。这一引理可以归纳性地证明，依赖于下列的观察。设点j是储存在储存桶k中，i是从点s到点j路上j的前继点，则d(s,i)≤d_s (j)-λ^*，所以点i可以从一些b<k的储存桶b中选择
定理1：包含点s的最小有向圈算法在O(m)时间内可以准确决定是否存在长度小于λ^* |W^* |且包含点s的有向圈。证明思路为：由引理1，该算法恰好与从点s出发的，路长最多为(|W^* |-1) λ^*的最短路问题一致，因此，它决定了长度最多为λ^* |W^* |且包含点s的最小有向圈。而浏览一条弧并且更新距离标号和储存桶需要O(1)时间，所以一共需要O(m)时间。而确认最小非空储存桶的时间共为O(n)这是因为最小非空储存桶的标号在迭代过程中非降，一次最多浏览一个空的储存桶，在选择操作中，平均每点产生顶点的时间是O(1)，所以总共的时间为O(n)。因此，当从点s出发浏览时算法运行时间是O(m)。
定理2：在图G中，算法MLDC能在O(mn)时间复杂度内计算最小有向圈。用Karp[14]算法来计算MMC问题需要O(mn)时间。距离p也可以在O(mn)时间内获得。从而，由定理1， n条独立的最短路问题中的每一条都可以在O(m)内获得，从而，MLDC算法可以在O(mn)时间复杂度内计算完。
在试验中加速时可以将一系列实践性的改善包含在算法中。这些改善对于最坏情况的运行时间没有影响，除了解决MMC问题的Howard算法，它确实能减小最坏运行时间。(1)在实践中用Howard算法算法来解决MMC问题。我们也实施了Ahuja和 Orlin的缩放单纯性算法，其运行时间为O(mnlognC)，且和Howard算法思想相似，但Howard算法更快，所以对实验实施Howard算法来加速运算。(2)用整数运算。在计算MLDC问题时，对于问题实例用距离标号时令C_ij-λ^*代替C_ij，设试验中最初的弧长都是整数的，则λ^*是有理数，且分母数小于n，为了避免出现分数弧长，用λ^*来乘以最初的弧长，这样可将原来问题转化为等价的整数问题。(3)在计算MMC之后停止计算。若最小平均圈W^*只有两条弧，则是一个最小圈。若W^*中只有三条弧，则它为最小弧或者有只有两条弧的最小圈，后者时只需要检查所有的两条弧的圈即可，这需要O(m)时间。(4)将算法限制到强连通分支上可以加速时间，因为强连通分支可以周期性地计算。在实际中，考虑下面的加速建议：1：删除下标小于s的点。对于s∈[1,n]，在实施包含点s的最小有向圈时永久性地删除点s和点s的入射弧。2：删除源点和汇点。算法记录每个点的出度和入度，永久性地删除出度入度为0的点。3：删除大长度弧。用mldc表示某次迭代中出现的最小圈，永久性地删除长度至少为mldc-λ^*的弧，因为该弧不可能出现在最短路中。4：初始化距离标号的向量。在计算第一个源点的最短路前，这一向量在每个位置都被设置成mldc-λ^*，当找到从一个给定点的最短路后，修正的动态距离标号包含每一个距离最多为mldc-λ^*的标号。在选择下一个源点之前，将该点的距离标号重新设置为mldc-λ^*，并且将动态标号变空。
四、实验结果
4.1算法
(1)Floyd与Warshall算法(FW算法)。d^k [i,j]表示点{i,j}⋃{1,2⋯k-1}上的点i到点j的最短路长度，有以下的递归关系：d^k [i,j]=min{d^k [i,j]，d^k [i,k]+d^k [k,j] }，对于固定的k，i，j变化时，可以不再计算d^k [i,j]≥mldc的部分加速计算。
(2)BHEAP(二叉堆算法)。当图包含负长度时，初始最优距离用Pollotino[16]计算，剩下来的最短路树用Dijkstra’s 算法的二叉堆来计算，所以得堆操作需要log(size(heap) )时间，计算时限制到一个强连通分支中比点s下标大的，距离标号小于mldc的弧。另外，考虑Dijkstra’s 算法的同时执行。当一个点的两个方向的最小距离下标之和比mldc更小时，将BHEAP算法编码为BBHEAP算法。
(3)在被称为OSN的算法中，有6个变式，通过乘以和除以|W^* |将所有的运算操作都等价于在整点上运算。
(4)针对被称为OSNB的Dial的漏桶算法，第一步，用双向关联列表来储存桶，可在O(1)时间内完成删除。当储存桶k是最小指标储存桶时，删除储存桶k的第一个元素；当点j从储存桶k移向一个更低指标的储存桶中时，也删除第一个元素。第二种情况需要双向列表。
(5)松弛的漏桶算法(OSNBL)。在松弛的漏桶算法中，用单向列表来实施储存桶，采用松弛删除法：若当点j从储存桶k移向一个更低指标的储存桶中时，将点j添加到更低指标的储存桶中，但是不从储存桶k中删除点j。
(6)在一个被称为OSNH的执行中，用二叉堆来执行最短路算法。OSNH与BHEAP的区别在于，前者首先使用了Howard算法来解决MCC且修正了距离，BHEAP并没有解决MCC问题。
最后三部操作与之前的三部一样，除了Dijkstra’s 算法的双向版本，分别称为BOSNB、BOSNBL、BOSNH。
4.2测试实例
第一组实例由NETGEN生成器生成，生成随机n∈{8000,16000⋯256000}，m∈{2n,4n⋯128n}，下面这些参数固定，源点为1，汇点为1，T源点=T汇点=0，最小花费为1，最大花费为10000，每个n和m的组合，产生10个复制，共产生420个不同的实例。第二组实例由GRIDGEN生成，N从8000到256000变化时，图的广度属于集合{2,4，,128}，下面这些参数固定，源点为1，汇点为1，平均度等于广度，最小花费为1，最大花费为10000，每个n和m的组合，产生11个复制，共产生462个不同的实例。第三组实例由HPGEN生成，首先构造一个长度在[1,10]的任选的哈密顿路，剩下的m-n+1条弧由随机选择的点对决定，这些弧的弧长在[1,10000]区间内任选，产生随机图，n和m与第一组中相同，每对组合产生10个复制，共产生420个不同的实例。试验中所有弧长非负。
4.3结果
1：Floyd与Warshall算法(FW算法)没有其他算法有竞争力。
2：Dijkstra’s 算法的双向变形并不能实质性地提高运行时间。
3：BOSNB、BOSNBL是OSN的6步执行中最快的，同时也比BHEAP算法的两步执行都快。
BOSNB与BOSNBL运行时间差不多，BOSNB比OSNH的运行时间快10%--30%。BOSNB比BHEAP运行更快，BHEAP是渐进更快地。
4：在实验中运行Howard算法会使得解决MLDC实例的时间缩短，尽管如此，Howard算法仍然是操作的阻碍。
5：当图的密度d=m/n值较小时，BHEAP比BOSNB快，而当d值较大时，情况相反。d的增大汇导致运行时间的增大，但是，d的加倍并不能导致任何算法的CPU运行时间的加倍，并且，Howard算法远远没有BHEAP算法那样对d值敏感。
6：回归分析提供了Howard算法和BOSNB-H的运行时间的一个很好的估计，并且给出了BHEAP运行时间的一个极好的估计。
7：BOSNB和BHEAP将比呈现出的最坏情况明显性地快速解决MLDC问题。回顾对MLDC问题最坏情况分析是O(mn)， BHEAP的运行时间是O(mnlogn)，若BHEAP回归分析对于n，m很大时是准确的，用d值来代替m/n，则BHEAP在实际中的运行时间将如下：对于NETGEN 实例是O(n^63 m^73 )GRIDGEN实例是O(n^5 m^76 )，对HPGEN实例是O(n^46 m^87 )，在每一种情况下都比最坏运行时间快得多。
:五、总结
我们初步的理论结果是对于解决有向图中的最小圈问题的一个更快的算法，我们在原有的最好时间O(mn+n^2 loglogn)的基础上将之改善为O(mn)时间复杂度内解决该问题。
在我们的计算中，依赖于Howard算法来解决MMC问题。一般的，在解决最小圈问题之前解决平均圈问题确实会导致在NETGEN、GRIDGEN、HPGEN三种生成器生成的实例的运行时间的加快。这种加快需要浏览确定每一条最短路计算的点的个数来实现。然而，由实验的结果，当n>1000000，d很小时，似乎可以略过MMC的结果，这似乎是更好地结果。
同时，可以观察到，在执行运算Dijkstra’s 算法的双向版本来解决MLDC问题可导致运行时间的略微提高。一个令人吃惊而又意义重大的结果是：Howard算法是解决MLDC问题的瓶颈操作尽管Howard算法在实际云中中比O(mn)时间复杂度快得多。对于MMC运行时间的加快会导致解决MLDC问题的运行时间的相应的提高。
六、参考文献
1. Ahuja, R., Magnanti, T., & Orlin, J.B., (1993). Network Flows, Prentice-Hall, Inc.
2. Bellman, R., (1958). On a Route Problem. Quart. Appl. Math. 16, 87-90.
3. Dasdan, A., (2004). Experimental Analysis of the Fastest Optimum Cycle Ratio and Mean Algorithms. ACM Transactions on Design Automation of Electronic Systems 9 (4),385-418.
4. Dasdan, A., Irani, S., & Gupta, R.K., (1999). Efficient Algorithms for Optimum Cycle Mean and Optimum Cost to Time Ratio Problem. In Proceedings of the 36th Design Automation Conference.
5. Denardo, E.V. and Fox, B.L. (1979). Shortest-route methods: 1. Reaching, pruning, and buckets. Operations Research 27-1, 161-186.
6. Dial, R.B., (1969). Algorithm 360: Shortest Path Forest with Topological Ordering. Communications of the ACM 12, 632-633.
7. Dijkstra., E.W., (1959). A Note on Two Problems in Connection with Graphs. Number . Math. 1, 269-271.
8. Floyd, R.W., (1962). Algorithm 97: shortest path. Comm. ACM 5, 345.
9. Ford, L.R., (1956). Network Flow Theory, The rand Corporation Report P-923, Santa Monica, Calif.
10. Fredman, M.L. & Tarjan, R.E., (1987). Fibonacci heaps and their uses in improved network optimization algorithms. Journal of the ACM (JACM) 34-3, 596-615.
11. Georgiadis, L., Goldberg, A.V., Tarjan, R.E., & Werneck, R. F., (2009). An experimental study of minimum mean cycle algorithms. In Proceedings of the 11th Workshop on Algorithm Engineering and Experiments (ALENEX). SIAM, Philadelphia, 1-13.
12. Howard, R.A. (1960). Dynamic Programming and Markov Processes. The M.I.T. Press, Cambridge, Mass.
13. Itai, A., & Rodeh, M., (1978). Finding a minimum circuit in a graph. SIAM J. Computing 7 (4), 413 - 423.
14. Karp, R.M., (1978). A characterization of the minimum cycle mean in a digraph. Discrete Mathematics 23, 309-11.
15. Moore, M.H., (1976). On the fastest route for convoy-type traffic in flow rate-constrained networks. Transportation Science,10 (2), 113 - 124.
16. Pallottino, s., (1984). Shortest-Path Methods: Complexity, In-terrelations and New Propositions. Networks 14, 257-267.
17. Pettie, S., (2004). A new approach to all-pairs shortest paths on real-weighted graphs. Theoretical Computer Science 312 (1), 47-74.
18. Roditty, L. & Williams, V.V. (2011). Minimum weight cycles and triangles: equivalences and algorithms. arXiv:1104.2882v1 14 Apr 2011.
19. Warshall, S., (1962). A theorem on Boolean matrices. Journal of the ACM 9 (1), 11-12.