Abstract
We present original average-case results on the performance of the Ford–Fulkerson maxflow algorithm on grid graphs (sparse) and random geometric graphs (dense). The analysis technique combines experiments with probability generating functions, stochastic context free grammars and an application of the maximum likelihood principle enabling us to make statements about the performance, where a purely theoretical approach has little chance of success. The methods lends itself to automation allowing us to study more variations of the Ford–Fulkerson maxflow algorithm with different graph search strategies and several elementary operations. A simple depth-first search enhanced with random iterators provides the best performance on grid graphs. For random geometric graphs a simple priority-first search with a maximum-capacity heuristic provides the best performance. Notable is the observation that randomization improves the performance even when the inputs are created from a random process.
Similar content being viewed by others
Notes
The tool can be downloaded at http://wwwagak.cs.uni-kl.de/malijan.html.
\([a..b]\) with \(a, b\in \mathbb {N}\) and \(a\le b\) is the integer interval, thus \([a..b]:=\{a, a+1, \ldots ,b\}\).
References
Aldrich, J.: R. A. Fisher and the making of maximum likelihood 1912–1922. Stat. Sci. 12(3), 162–176 (1997)
Chandran, B., Hochbaum, D.: A computational study of the pseudoflow and push-relabel algorithms for the maximum flow problem. Oper. Res. 57, 358–376 (2009)
Chi, T., Geman, S.: Estimation of probabilistic context-free grammars. Comput. Linguist. 24(2), 299–308 (1998)
Chomsky, N., Schützenberger, M.P.: The algebraic theory of context-free languages. In: Braffort, P., Hirschberg, D. (eds.) Computer Programming and Formal Languages, pp. 118–161. North Holland, New York (1963)
David, H.A., Nagaraja, H.N.: Order Statistics, 3rd edn. Wiley, New York (2003)
Dinic, E.A.: An algorithm for the solution of the max-flow problem with the polynomial estimation. Sov. Math. Dokl. 11, 1277–1280 (1970). Dokl. Akad. Nauk SSSR 194, 1970, no. 4 (in Russian)
Durstenfeld, R.: Algorithm 235: random permutation. Commun. ACM 7(7), 420 (1964)
Edmonds, J., Karp, R.M.: Theoretical improvements in algorithmic efficiency for network flow problems. J. ACM 19(2), 248–264 (1972)
Erdos, P., Rényi, A.: On random graphs I. Publ. Math. Debr. 6, 290–297 (1959)
Flajolet, P., Sedgewick, R.: Analytic Combinatorics. Cambridge University Press, Cambridge (2009)
Ford, L.R., Fulkerson, D.R.: Flows in Networks. Princeton University Press, Princeton (1962)
Gilbert, E.N.: Random graphs. Ann. Math. Stat. 30, 1141–1144 (1959)
Goldberg, A.V., Rao, S.: Beyond the flow decomposition barrier. J. Assoc. Comput. Mach. 45, 753–782 (1998)
Goldfab, D., Grigoriads, M.D.: A computational comparison of the dinic and network simplex methods for maximum flow. Ann. Oper. Res. 13, 83–123 (1988)
Knuth, D.E.: The Art of Computer Programming, Volume 1: Fundamental Algorithms, 3rd edn. Addison Wesley, Reading (1997)
Knuth, D.E.: The Art of Computer Programming, Volume 2: Seminumerical Algorithms, 3rd edn. Addison Wesley, Reading (1998)
Laube, U., Nebel, M.E.: Maximum likelihood analysis of algorithms and data structures. Theor. Comput. Sci. 411(1), 188–212 (2010)
Motwani, R.: Average-case analysis of algorithms for matchings and related problems. J. ACM 41(6), 1329–1356 (1994)
Orlin, J.B.: Max flows in \({O}(nm)\) time, or better. In: STOC ’13: Proceedings of the 45th Annual ACM Symposium on Symposium on Theory of Computing, pp. 765–774. ACM, New York, NY, USA (2013)
Penrose, M.: Random Geometric Graphs. Oxford University Press, Oxford (2003)
Sedgewick, R.: Putting the science back into computer science. www.cs.princeton.edu/rs/talks/ScienceCS10.pdf (2010). Accessed 25 Nov 2014
Sedgewick, R.: The role of the scientific method in programming. www.cs.princeton.edu/rs/talks/ScienceCS.pdf (2010). Accessed 25 Nov 2014
Sedgewick, R., Wayne, K.: Algorithms, 4th edn. Addison Wesley, Reading (2011)
Vizing, V.G.: The cartesian product of graphs. Vyčisl. Sistemy 9, 30–43 (1963)
Wild, S., Nebel, M., Reitzig, R., Laube, U.: Engineering java 7’s dual pivot quicksort using malijan. In: ALENEX 2013: Proceedings of the Meeting on Algorithm Engineering and Experiments, pp. 55–70. New Orleans, USA (2013)
Acknowledgments
This work was inspired by the slides of two talks of R. Sedgewick on “The Role of the Scientific Method in Programing” [22] and “Putting the Science back into Computer Science” [21]. The authors would like to thank Sebastian Wild, Raphael Reitzig, Michael Holzhauser, Vasil Tenev and Florian Furbach for their ongoing effort in the development of MaLiJAn, the tool to do the maximum likelihood average case analysis of Java bytecode programs semi-automatically.
Author information
Authors and Affiliations
Corresponding author
Additional information
This work is supported by DFG Grants NE 1379/2-1 and NE 1379/3-1.
Appendix
Appendix
1.1 Grid Graph Properties: Proof of Lemma 2
The grid graph \(G_{j, k}\) has the following properties:
there are 4 vertices of degree 2, \(2k+2j-8\) vertices of degree 3 and \((k-2)(j-2)\) vertices of degree 4.
Proof
Assume without loss of generality, that the underlying path graphs \(P_j\) and \(P_k\) of the grid graph \(G_{j,k}=P_j\mathrel {\square }P_k\) are not only defined on \([1..j]\) and \([1..k]\) respectively, but the elements of the edge sets have the form \((i, i+1)\). According to the definition of the grid graph its edges then have the form \(\bigl \{(u, i),(u,i+1)\bigr \}\) or \(\bigl \{(i, v), (i+1,v)\bigr \}\), \(u\in [1..k], v\in [1..j]\), intuitively the vertical and horizontal edges of the grid graph.
The \(2\cdot (k-1+j-1)\) vertices on the “edges” of the grid are of the form \((1, \bullet ),(k,\bullet ), (\bullet ,j)\) and \((\bullet , 1)\). Four of them are on two “edges” of the grid: \((1, 1), (1, j), (k, j)\) and \((k, 1)\) thus in the “corners” of the grid and have degree 2, as there are only two ways to increase or decrease each element of the pairs and stay within the ranges \([1..k]\) and \([1..j]\), the others have 3 degrees of freedom. The remaining \((k-2)(j-2)\) inner vertices have 4 degrees of freedom.
Thus \({|}{V(G_{j, k})}{|}=4+2\cdot (k-2+j-2)+(k-2)(j-2)=4+2k+2j-8+jk-2j-2k+4=j\cdot k\) and counting each edge from both ends yields
1.2 Random Geometric Graph Properties: Proof of Theorem 2
The average degree of an arbitrary node in a random geometric graph \(v\in V\bigl (G_r(\fancyscript{X}_n)\bigr )\) on point set \(\fancyscript{X}_n\) uniformly distributed over the unit square is:
The expected degree of a vertex \(v\) depends on its position in the unit square. In the center \([r,1-r]^2\) square we have with \(r=15/100\)
as this fraction of all \(n\) nodes in the disk around node \(v\) are connected to \(v\). This is less near the boundary as no full disk with radius \(r=15/100\) fits in there. On the boundary we have
as there is only half a disk possible. In the corner only a quarter disk fits, thus
Thus to proof the theorem we have to determine the average size of the neighborhood including the boundary effects when the vertex \(v\) is near the boundaries of the unit square.
Proof
To do this we integrate over the unit square to determine the average size of the neighborhood \(p(r)\). Exploiting the symmetry of the unit square, integrating over a single quadrant \([0,1/2]^2\) is enough:
The contribution \(A(x, y, r)\) obviously varies with the location \((x, y)\) and the radius \(r\) of the neighborhood. We can distinguish four rectangles in the quadrant by splitting it horizontally and vertically in the point \((r,r)\). The large square \([r, 1/2]^2\) allows for a full disk at all positions. In the two rectangles \([0, r]\times [r, 1/2]\) and \([r, 1/2]\times [0, r]\) a segment of the disk is cut off by the boundary. And on the small square \([0, r]^2\) this happens on both boundaries. Thus we arrive at
The first integral is easily solved as the contribution does not depend on the location \((x, y)\).
Thus the contribution \(\pi r^2\) is weighted with the area of the larger square \((1/2-r)^2\).
The two rectangles are easy too, as the contribution only depends on the change of the location in one dimension.
The contribution \(A(r, y)\) is just the area of a disk minus the segment that protrudes over the boundary by \(r-y\).
The area of the segment \(A_{SG}(r, y)\) is just the area of the sector of the disk less the equal sided triangle formed by the chord and the radii whose height is \(y\) when the segment protrudes by \(r-y\). Thus
The length of half the chord ist just
by Pythagoras’ theorem. Accordingly
The area of the sector \(A_{SK}(r, y)\) depends on its apex angle \(\alpha \):
We have to express this angle by \(y\) and \(r\). Another application of Pythagoras’ theorem shows that half the angle is just \(\frac{\alpha }{2}=\arccos \left( 1-\frac{r-y}{r}\right) \) yielding
The final contribution is
and integrating over it yields
This is again of the form where the contribution \(r^2(\pi -2/3)\) is weighted with area of the two rectangles \(r(1/2-r)\). The antiderivatives used above appear again in the next cases and are derived afterwards.
The small square is more difficult to solve. When the relation \(\sqrt{x^2+y^2}\ge r\) holds the easier case is encountered. The segments of the disk protruding over two boundaries do not overlap and they can simply be subtracted as done before. However when \(\sqrt{x^2+y^2}< r\) holds the overlapping part would be subtracted twice. So we add a correction term. Separating the two cases gives
The integrands are
To determine \(A_\ge (r, x, y)\) we already know all the necessary expressions and get:
Thus it remains to determine the correction term \(K(r, x, y)\). It describes a part of a disk that looks like a segment offset from the center. We start the derivation with a full disk. In Cartesian coordinates we have
because when moving along the diameter \([-r, r]\) we can only move as far as \([-\sqrt{r^2-x^2}, \sqrt{r^2-x^2}]\) in the perpendicular direction to stay with in the disk. Halving that interval gives a semi disk
Halving the other direction too yields a quarter disk
Limiting the integral further by cutting of stripes with width \(0<a<r\) in one and \(0<b<r\) in the other direction we get:
Obviously \(r>\sqrt{a^2+b^2}\) has to be true otherwise the area would be empty. Thus
Solving the remaining integrals
and
and combining the above, we get
Thus once more the contribution \(r^2(\pi -29/24)\) is weighted by the area \(r^2\). With all integrals solved we can combine the results in (5) and arrive at
for the average size of the neighborhood. As expected with \(r=15/100\) it is \(p(r)\approx 0.062\), which is a little less than a full disk \(\pi r^2\approx 0.071\) due to the boundary effects. \(\square \)
The following five integrals have been used in the previous proof.
Lemma 6
Proof
Starting with integration by parts
with
yields an easier integral. We must first determine \(v'\). Let \(y=\arcsin (x)\) then \(x=\sin (y)\) and use implicit differentiation. Then solve for \(dy/dx\) and apply Pythagoras’ Theorem and substitute \(x=\sin (y)\).
Hence we must now solve the integral on the right-hand side of
Subsituting with \(k=1-x^2\) thus
Backsubstitution yields
Inserting this into (8) gives the result. \(\square \)
Lemma 7
Proof
Substituting with \(k=r^2-y^2\) hence
\(\square \)
Lemma 8
Proof
Starting with integration by parts
with
yields an easier integral. We must first determine \(v'\). Let \(y=\arcsin (x/r)\) then \(x/r=\sin (y)\) and use implicit differentiation. Then solve for \(dy/dx\) and apply Pythagoras’ Theorem and substitute \(x/r=\sin (y)\).
Or with Lemma 6 from before the chain rule gives us
Hence we must now solve the integral on the right-hand side of
Subsituting with \(k=r^2-x^2\) thus
Backsubstitution yields
Inserting this into (9) gives the result. \(\square \)
Lemma 9
Proof
Starting with integration by parts
with
yields an easier integral. We must first determine \(v'\). Using the chain rule we find
Hence we must now solve the integral on the right-hand side of
Substituting \(k=r^2-x^2\) thus
Backsubstitution yields
Inserting this in (10) gives the result. \(\square \)
Lemma 10
Proof
By linearity we get
And we can continue with the standard trigonometric substitution \(x/r=\sin (y)\), thus \(y=\arcsin (x/r)\) and we continue with:
Using the cosine double angle formula
we remove the squared cosine
Using the sine double angle formula \(\sin (2y)=2\sin (y)\cos (y)\), and Pythagoras’ Theorem to replace the cosine
we can rewrite the right-hand side of (12)
Finally reverting the substitution with \(y=\arcsin (x/r)\)
Inserting this in (12) and its result in (11) completes the proof.\(\square \)
1.3 Average Edge Capacities: Proof of Lemma 3
The average capacity of an edge \(e\in E(G_k)\) is
with \(\epsilon \ll 1\) and \({\varPhi }(x)\) denoting the cdf of the standard normal distribution.
Proof
The capacity of an edge is chosen in the interval \([0..\infty ]\) according to a truncated normal distribution. \(250-\lfloor 50\cdot X\rfloor \) where \(X\sim \fancyscript{N}(0, 1)\). Truncation is due to the fact, that negative capacities are pointless. Whenever \(250-\lfloor 50\cdot X\rfloor <0\) is true we discard the result and draw again. Due to
we accept only draws within the half-open interval \(X\in [-\infty , 251/50)\), consequently the normal distribution is truncated to this interval. The pdf of the truncated normal distribution is
With \(\phi (x)\) being the pdf and \({\varPhi }(x)\) denoting the cdf of the standard normal distribution \(\fancyscript{N}(0, 1)\)
The pdf and cdf of a random variable \(X\) normal distributed with mean \(\mu \) and variance \(\sigma ^2\), that is \(X\sim \fancyscript{N}(\mu , \sigma ^2)\), can be written as
in terms of the standard normal distribution. This is due to the standard normalization \(Z=\frac{X-\mu }{\sigma }\), that turns an RV \(X\) with mean \(\mu \) and variance \(\sigma ^2\) into a RV \(Z\) with mean 0 and variance 1. Conversely \(X=\mu \pm \sigma Z\) can be used to transform a standard normal random variable into one with mean \(\mu \) and variance \(\sigma ^2\).
Integrating
confirms that we have a probability distribution again.
Now for the expectation
we deal with the floor function by treating it as a piecewise constant function. This step function allows us to rewrite the integral as a sum over the steps.
Thus we integrate
Checking that we have a distribution
clearly telescopes to
Thus the expectation now simply is
\(\square \)
1.4 Average Maxflow in a Square Grid Graph: Proof of Lemma 5
The average maxflows in a square grid graph \(G_k\) are bounded by:
Proof
There is no more flow possible than the minimum of the total capacities of source and sink. Due to the random edge weights there may not be enough edges in between with a high enough total capacity. As these bottlenecks may only lower the maxflow we ignore them and proceed with:
with \(c=\deg (s)=\deg (t)\). This expectation is the first moment of the first order statistics of two iid discrete random variables, it can be expressed via the “tail” of the cdf of the discrete random variables, see [5]:
In case of the discrete uniform distribution \(\fancyscript{U}\) on \([0..100]\) we find:
In case of the gauss distribution \(X\sim \fancyscript{N}(0,1)\) we redraw whenever \(250-\lfloor 50\cdot X\rfloor <0\) is true. Due to
we accept only draws within the half-open interval \(X\in [-\infty , 251/50)\). The truncated pdf \(f_X(x)\) of the standard normal distribution is then:
therefore its cdf is
Transforming with \(Z=250-50X\) yields:
Thus
where \(\phi (x)\) and \({\varPhi }(x)\) are the pdf and cdf of the standard normal distribution respectively. \(\square \)
1.5 DFS, BFS and PFS Implementations of the Graph Search
The annotation produceCost("EdgeFlowChanged", 1) inside the method addflowRto() is used to keep track of the updates made to the graph data structure. The next code fragments show the similarity of the dfs() and bfs() routines, just the stack is replaced by a queue. The part where the iterator visits the incident edges of a vertex is shown in both versions. First we present the code for depth-first-search:
Next the code for breadth-first-search is given:
The methods push()/put() and pop()/get() are annotated to be able to track them.
The maximum capacity augmenting path heuristic is implemented via a priority queue that uses the leftover capacity as the priority.
Rights and permissions
About this article
Cite this article
Laube, U., Nebel, M.E. Maximum Likelihood Analysis of the Ford–Fulkerson Method on Special Graphs. Algorithmica 74, 1224–1266 (2016). https://doi.org/10.1007/s00453-015-9998-5
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s00453-015-9998-5