Maximum Likelihood Analysis of the Ford–Fulkerson Method on Special Graphs

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.

Appendix

1.1 Grid Graph Properties: Proof of Lemma 2

The grid graph $G_{j, k}$ has the following properties:

$$\begin{aligned} {|}{V(G_{j, k})}{|}=j \cdot k, \quad {|}{E(G_{j, k})}{|}=(j-1)k+j(k-1)=2j \cdot k-j-k, \end{aligned}$$

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

$$\begin{aligned} 2\cdot {|}{E(G_{j, k})}{|}= & {} 4\cdot 2+3\cdot 2\cdot (k-2+j-2)+4\cdot (k-2)(j-2)\\= & {} 4j\cdot k-2j-2k. \end{aligned}$$

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:

$$\begin{aligned} \mathbb {E} \bigl [\hbox {deg}(v)\bigr ]=n\cdot p(r) \quad \hbox {with} \quad p(r) = r^2\left( \pi +\frac{r(3r-16)}{6}\right) . \end{aligned}$$

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$

$$\begin{aligned} \mathbb {E}[\hbox {deg}(v)]=n\frac{\pi r^2}{1}\approx 0.071n \end{aligned}$$

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

$$\begin{aligned} \mathbb {E}[\hbox {deg}(v)]=n\frac{\pi r^2}{2} \approx 0.035 n \end{aligned}$$

as there is only half a disk possible. In the corner only a quarter disk fits, thus

$$\begin{aligned} \mathbb {E}[\hbox {deg}(v)]=n\frac{\pi r^2}{4}\approx 0.018n. \end{aligned}$$

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:

$$\begin{aligned} p(r)=4 \cdot \int _0^{\frac{1}{2}}\int _0^{\frac{1}{2}}A(r, x, y)\,dxdy. \end{aligned}$$

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

$$\begin{aligned} p(r)\!=\!4\cdot \left( \int _r^{\frac{1}{2}}\!\int _r^{\frac{1}{2}}A(r)\,dxdy \!+\! 2\cdot \int _0^r\!\int _r^{\frac{1}{2}}A(r, y)\,dxdy \!+\! \int _0^r\int _0^r\!A(r, x, y)\,dxdy\!\right) \!. \end{aligned}$$

(5)

The first integral is easily solved as the contribution does not depend on the location $(x, y)$.

$$\begin{aligned} \int _r^{\frac{1}{2}}\int _r^{\frac{1}{2}}A(r)\,dxdy&=A(r)\int _r^{\frac{1}{2}}\int _r^{\frac{1}{2}}1\,dxdy=\pi r^2\int _r^{\frac{1}{2}}\int _r^{\frac{1}{2}}1\,dxdy\\ {}&=\pi r^2\cdot \left( \frac{1}{2}-r\right) ^2. \end{aligned}$$

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.

$$\begin{aligned} 2\cdot \int _0^r\int _r^{\frac{1}{2}}A(r, y)\,dxdy= & {} 2\cdot \int _0^rA(r, y)\int _r^{\frac{1}{2}}1\,dxdy\\= & {} 2\cdot \int _0^rA(r,y)\left( \frac{1}{2}-r\right) \,dy=(1-2r)\int _0^rA(r, y)\,dy. \end{aligned}$$

The contribution $A(r, y)$ is just the area of a disk minus the segment that protrudes over the boundary by $r-y$.

$$\begin{aligned} A(r, y)=A_\circ (r)-A_{SG}(r, y)=\pi r^2-A_{SG}(r,y). \end{aligned}$$

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

$$\begin{aligned} A_{SG}(r, y)= & {} A_{SK}(r, y)-A_\varDelta (r, y),\\ A_\varDelta (r, y)= & {} \frac{sy}{2}. \end{aligned}$$

The length of half the chord ist just

$$\begin{aligned} \frac{s}{2}=\sqrt{r^2-y^2} \end{aligned}$$

by Pythagoras’ theorem. Accordingly

$$\begin{aligned} A_\varDelta (r, y)=y \cdot \sqrt{r^2-y^2}. \end{aligned}$$

The area of the sector $A_{SK}(r, y)$ depends on its apex angle $\alpha$:

$$\begin{aligned} A_{SK}(r,y) = \pi r^2\frac{\alpha }{2\pi }=\frac{\alpha }{2}r^2. \end{aligned}$$

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

$$\begin{aligned} A_{SG}(r, y)&=r^2\arccos \left( 1-\frac{r-y}{r}\right) -y\cdot \sqrt{r^2-y^2}\\&=r^2\arccos \left( \frac{y}{r}\right) -y\cdot \sqrt{r^2-y^2}. \end{aligned}$$

The final contribution is

$$\begin{aligned} A(r, y)&=\pi r^2-\left( r^2\arccos \left( \frac{y}{r}\right) -y\cdot \sqrt{r^2-y^2}\,\right) \\&= r^2\left( \pi -\arccos \left( \frac{y}{r}\right) \right) +y\cdot \sqrt{r^2-y^2}, \end{aligned}$$

and integrating over it yields

$$\begin{aligned} (1-2r)\int _0^rA(r, y)\,dy&=(1-2r)\left[ \int _0^rr^2\left( \pi -\arccos \left( \frac{y}{r}\right) \right) +y\cdot \sqrt{r^2-y^2}\,dy\right] \\&=(1-2r)\left[ r^2\!\int _0^r\!\pi \!-\!\arccos \left( \frac{y}{r}\right) \,dy\!+\!\int _0^r\!y\cdot \sqrt{r^2-y^2}\,dy\!\right] \\&=(1-2r)\left[ r^2\bigl (\pi r-r\bigr )+\frac{r^3}{3}\right] \\&=(1-2r)\left[ r^3(\pi -1)+\frac{r^3}{3}\right] \\&=(1-2r)r^3\left( (\pi -1)+\frac{1}{3}\right) \\&=(1-2r)r^3\left( \pi -\frac{2}{3}\right) \\&=2\left( \frac{1}{2}-r\right) r\cdot r^2\left( \pi -\frac{2}{3}\right) . \end{aligned}$$

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

$$\begin{aligned}&\int _0^r \int _0^rA (r, x, y)\,dydx\nonumber \\&\quad = \int _0^r \left[ \int _0^{\sqrt{r^2-x^2}}A_<(r, x, y)\,dy+ \int _{\sqrt{r^2-x^2}}^rA_\ge (r, x, y)\,dy\right] dx. \end{aligned}$$

(6)

The integrands are

$$\begin{aligned} A_\ge (r, x, y)&=\pi r^2-A_{SG}(r, x) -A_{SG}(r, y)\nonumber \\ A_<(r, x, y)&=A_\ge (r, x, y)+K(r, x, y)\nonumber \\&=\pi r^2-A_{SG}(r, x)-A_{SG}(r, y)+K(r, x, y). \end{aligned}$$

(7)

We can use (7) to rewrite (6)

$$\begin{aligned} \int _0^r\int _0^rA(r, x, y)\,dydx&=\int _0^r\left[ \int _0^{\sqrt{r^2-x^2}}A_\ge (r, x, y)+K(r, x, y)\,dy\right. \\&\quad \ \left. +\int _{\sqrt{r^2-x^2}}^rA_\ge (r, x, y)\,dy\right] dx\\&=\int _0^r\left[ \int _0^{\sqrt{r^2-x^2}}K(r, x, y)\,dy+\int _0^rA_\ge (r, x, y)\,dy\right] dx. \end{aligned}$$

To determine $A_\ge (r, x, y)$ we already know all the necessary expressions and get:

$$\begin{aligned} A_\ge (r, x, y)&=\pi r^2-r^2\arccos \left( \frac{y}{r}\right) +y\cdot \sqrt{r^2\!-\!y^2}\!-\!r^2\arccos \left( \frac{x}{r}\right) +x\cdot \sqrt{r^2-x^2}\\&=r^2\left( \pi -\arccos \left( \frac{y}{r}\right) -\arccos \left( \frac{x}{r}\right) \right) +y\cdot \sqrt{r^2-y^2}+x\cdot \sqrt{r^2-x^2}\\&=r^2\left( \arcsin \left( \frac{y}{r}\right) +\arcsin \left( \frac{x}{r}\right) \right) +y\cdot \sqrt{r^2-y^2}+x\cdot \sqrt{r^2-x^2}. \end{aligned}$$

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

$$\begin{aligned} \int _{-r}^{r}\int _{-\sqrt{r^2-x^2}}^{\sqrt{r^2-x^2}}1\,dydx=\pi r^2, \end{aligned}$$

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

$$\begin{aligned} \int _{-r}^{r}\int _0^{\sqrt{r^2-x^2}}1\,dydx=\frac{\pi r^2}{2}. \end{aligned}$$

Halving the other direction too yields a quarter disk

$$\begin{aligned} \int _0^{r}\int _0^{\sqrt{r^2-x^2}}1\,dydx=\frac{\pi r^2}{4}. \end{aligned}$$

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:

$$\begin{aligned} K(r, a, b)=\int _a^{\sqrt{r^2-b^2}}\int _b^{\sqrt{r^2-x^2}}1\,dydx. \end{aligned}$$

Obviously $r>\sqrt{a^2+b^2}$ has to be true otherwise the area would be empty. Thus

$$\begin{aligned} K(r, a, b)&=\int _a^{\sqrt{r^2-b^2}}\sqrt{r^2-x^2}-b\,dx\\&=\left[ \frac{r^2}{2}\arcsin (x/r)+\frac{x}{2}\sqrt{r^2-x^2}-bx\right] _a^{\sqrt{r^2-b^2}}\\&=\left( \frac{r^2}{2}\arcsin \left( \sqrt{r^2-b^2}/r\right) +\frac{\sqrt{r^2\!-\!b^2}}{2}\sqrt{r^2\!-\!\sqrt{r^2-b^2}^2}\!-\!b\sqrt{r^2\!-\!b^2}\right) \\&\quad \ -\left( \frac{r^2}{2}\arcsin (a/r)+\frac{a}{2}\sqrt{r^2-a^2}-ba\right) \\&=\frac{r^2}{2}\arcsin \left( \frac{\sqrt{r^2-b^2}}{r}\right) +\frac{1}{2}\sqrt{r^2-b^2}\sqrt{r^2-r^2+b^2}-b\sqrt{r^2-b^2}\\&-\frac{r^2}{2}\arcsin \left( \frac{a}{r}\right) -\frac{a}{2}\sqrt{r^2-a^2}+ba\\&=ba-\frac{a}{2}\sqrt{r^2-a^2}+\sqrt{r^2-b^2}\frac{\sqrt{b^2}}{2}-b\sqrt{r^2-b^2}\\&\quad \ +\frac{r^2}{2}\arcsin \left( \frac{\sqrt{r^2-b^2}}{r}\right) -\frac{r^2}{2}\arcsin \left( \frac{a}{r}\right) \\&=ba-\frac{a}{2}\sqrt{r^2-a^2}+\frac{b}{2}\sqrt{r^2-b^2}-b\sqrt{r^2-b^2}\\&\quad \ +\frac{r^2}{2}\left( \arcsin \left( \frac{\sqrt{r^2-b^2}}{r}\right) -\arcsin \left( \frac{a}{r}\right) \right) \\&=ba-\frac{a}{2}\sqrt{r^2-a^2}-\frac{b}{2}\sqrt{r^2-b^2}\\&\quad \ +\frac{r^2}{2}\left( \arcsin \left( \frac{\sqrt{r^2-b^2}}{r}\right) -\arcsin \left( \frac{a}{r}\right) \right) . \end{aligned}$$

Solving the remaining integrals

$$\begin{aligned} \int _0^rA_\ge (r, x, y)\,dy&=\int _0^rr^2\left( \arcsin \left( \frac{y}{r}\right) +\arcsin \left( \frac{x}{r}\right) \right) +y\cdot \sqrt{r^2-y^2}\\&\quad \ +x\cdot \sqrt{r^2-x^2}\,dy\\&=r^2\arcsin \left( \frac{x}{r}\right) \left[ y\right] _0^r+r^2\int _0^r\arcsin \left( \frac{y}{r}\right) \,dy+x\sqrt{r^2-x^2}[y]_0^r\\&\quad \ +\int _0^ry\sqrt{r^2-y^2}\,dy\\&=r^3\arcsin \left( \frac{x}{r}\right) \!+\!r^2\left[ y\arcsin \left( \frac{y}{r}\right) +\sqrt{r^2\!-\!y^2}\right] _0^r+rx\sqrt{r^2\!-\!x^2}\\&\quad \ +\left[ -\frac{1}{3}(r^2-y^2)^{3/2}\right] _0^r\\&=r^3\arcsin \left( \frac{x}{r}\right) +\frac{r^3\pi }{2}-r^3+rx\sqrt{r^2-x^2}+\frac{r^3}{3}\\&=r^3\arcsin \left( \frac{x}{r}\right) +r^3\left( \frac{\pi }{2}-\frac{2}{3}\right) +rx\sqrt{r^2-x^2} \end{aligned}$$

and

$$\begin{aligned}&\int _0^{\sqrt{r^2-x^2}}K(r, x, y)\,dy\\&\quad =\int _0^{\sqrt{r^2-x^2}}xy-\frac{x}{2}\sqrt{r^2-x^2}-\frac{y}{2}\sqrt{r^2-y^2}\\&\qquad \ +\frac{r^2}{2}\left( \arcsin \left( \frac{\sqrt{r^2-y^2}}{r}\right) -\arcsin \left( \frac{x}{r}\right) \right) \,dy\\&\quad =x\bigl [y^2/2\bigr ]_0^{\sqrt{r^2-x^2}}-\frac{x}{2}\sqrt{r^2-x^2}[y]_0^{\sqrt{r^2-x^2}} -\frac{1}{2}\int _0^{\sqrt{r^2-x^2}}y\sqrt{r^2-y^2}\,dy\\&\qquad \ +\frac{r^2}{2}\int _0^{\sqrt{r^2-x^2}}\arcsin \left( \frac{\sqrt{r^2-y^2}}{r}\right) \,dy-\frac{r^2}{2}\arcsin \left( \frac{x}{r}\right) [y]_0^{\sqrt{r^2-x^2}}\\&\quad =\frac{x}{2}(r^2-x^2)-\frac{x}{2}(r^2-x^2)-\frac{1}{2}\left[ -\frac{1}{3}(r^2-y^2)^{3/2}\right] _0^{\sqrt{r^2-x^2}}\\&\qquad \ -\frac{r^2}{2}\sqrt{r^2-x^2}\arcsin \left( \frac{x}{r}\right) +\frac{r^2}{2}\left[ y\arcsin \left( \sqrt{1-\frac{y^2}{r^2}}\,\right) +\sqrt{r^2-y^2}\right] _0^{\sqrt{r^2-x^2}}\\&\quad =\frac{x^3}{6}-\frac{r^3}{6}-\frac{r^2}{2}\sqrt{r^2-x^2}\arcsin \left( \frac{x}{r}\right) \\&\qquad \ +\frac{r^2}{2}\sqrt{r^2-x^2}\arcsin \left( \sqrt{1-\frac{\sqrt{r^2-x^2}^2}{r^2}}\,\right) -\sqrt{r^2-\sqrt{r^2-x^2}^2}+\frac{r^2}{2}\sqrt{r^2}\\&\quad =\frac{x^3}{6}-\frac{r^3}{6}-\frac{r^2}{2}\sqrt{r^2-x^2}\arcsin \left( \frac{x}{r}\right) +\frac{r^2}{2}\sqrt{r^2-x^2}\arcsin \left( \frac{x}{r}\right) -\frac{xr^2}{2}+\frac{r^3}{2}\\&\quad =\frac{x^3}{6}+\frac{r^3}{3}-\frac{xr^2}{2} \end{aligned}$$

and combining the above, we get

$$\begin{aligned} \int _0^r\int _0^rA(r, x, y)\,dydx&=\int _0^rr^3\arcsin \left( \frac{x}{r}\right) +r^3\left( \frac{\pi }{2}-\frac{2}{3}\right) \\&\quad \ +rx\sqrt{r^2-x^2}+\frac{x^3}{6}+\frac{r^3}{3}-\frac{xr^2}{2}\,dx\\&=r^3\int _0^r\arcsin \left( \frac{x}{r}\right) \,dx+r^3\left( \frac{\pi }{2}-\frac{2}{3}\right) [x]_0^r\\&\quad \ +r\int _0^rx\sqrt{r^2-x^2}\,dx+\frac{1}{6}\left[ \frac{x^4}{4}\right] _0^r+\frac{r^3}{3}[x]_0^r-\frac{r^2}{2}\left[ \frac{x^2}{2}\right] _0^r\\&=r^3\left[ x\arcsin \left( \frac{x}{r}\right) +\sqrt{r^2-x^2}\right] _0^r+r^4\left( \frac{\pi }{2}-\frac{2}{3}\right) \\&\quad \ +r\left[ -\frac{1}{3}(r^2-x^2)^{3/2}\right] _0^r+\frac{r^4}{24}+\frac{r^4}{3}-\frac{r^4}{4}\\&=\frac{\pi }{2}r^4-r^4+r^4\left( \frac{\pi }{2}-\frac{1}{3}\right) +\frac{r^4}{3}+\frac{r^4}{24}-\frac{r^4}{4}\\&=\pi r^4-\frac{24}{24}r^4+\frac{1}{24}r^4-\frac{6}{24}r^4=r^2\cdot r^2\left( \pi -\frac{29}{24}\right) . \end{aligned}$$

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

$$\begin{aligned} p(r)&=4\cdot \left( \pi r^2\cdot \left( \frac{1}{2}-r\right) ^2+(1-2r)r^3\left( \pi -\frac{2}{3}\right) + r^4\left( \pi -\frac{29}{24}\right) \right) \\&=\pi r^2\left( 1-2r\right) ^2+(1-2r)r^3\left( 4\pi -\frac{8}{3}\right) + r^4\left( 4\pi -\frac{29}{6}\right) \\&=\pi r^2-4\pi r^3+4\pi r^4+r^3\left( 4\pi -\frac{8}{3}\right) -r^4\left( 8\pi -\frac{16}{3}\right) +r^4\left( 4\pi -\frac{29}{6}\right) \\&=\pi r^2-\frac{16}{6}r^3+\frac{32}{6}r^4-\frac{29}{6}r^4\\&=r^2\left( \pi +\frac{3r^2-16r}{6}\right) =r^2\left( \pi +\frac{r(3r-16)}{6}\right) \end{aligned}$$

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

$$\begin{aligned} \int \arcsin (x)\,dx=x\arcsin (x)+\sqrt{1-x^2}+C \end{aligned}$$

Proof

Starting with integration by parts

$$\begin{aligned} \int u'v=uv-\int uv' \end{aligned}$$

with

$$\begin{aligned} u'=1\quad \leadsto \quad u=x\quad \hbox {and}\quad v=\arcsin (x)\quad \leadsto \quad v'=\frac{d}{dx}\arcsin (x) \end{aligned}$$

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)$.

$$\begin{aligned} \frac{d}{dx}x&=\frac{d}{dx}\sin (y),\\ 1&=\frac{d}{dx}y\cos (y),\\ \frac{d}{dx}\arcsin (x)=\frac{d}{dx}y&=\frac{1}{\cos (y)}=\frac{1}{\sqrt{1-\sin ^2(y)}}=\frac{1}{\sqrt{1-x^2}}. \end{aligned}$$

Hence we must now solve the integral on the right-hand side of

$$\begin{aligned} \int \arcsin (x)\,dx=x\arcsin (x)-\int \frac{x}{\sqrt{1-x^2}}\,dx. \end{aligned}$$

(8)

Subsituting with $k=1-x^2$ thus

$$\begin{aligned} \frac{d}{dx}k= & {} \frac{d}{dx}1-x^2=-2x \quad \leadsto \quad \frac{dk}{-2x}=dx,\\ \int \frac{x}{\sqrt{1-x^2}}\,dx= & {} \int \frac{x}{\sqrt{k}}\frac{dk}{-2x}=-\frac{1}{2}\int \frac{1}{\sqrt{k}}\,dk=-\frac{1}{2}\int k^{-1/2} \,dk\\ {}= & {} -\frac{1}{2}\frac{k^{1/2}}{1/2}=-\sqrt{k}. \end{aligned}$$

Backsubstitution yields

$$\begin{aligned} \int \frac{x}{\sqrt{1-x^2}}\,dx=-\sqrt{1-x^2}. \end{aligned}$$

Inserting this into (8) gives the result. $\square$

Lemma 7

$$\begin{aligned} \int y\sqrt{r^2-y^2}\,dy=-\frac{1}{3}\bigl (r^2-y^2\bigr )^{3/2}+C \end{aligned}$$

Proof

Substituting with $k=r^2-y^2$ hence

$$\begin{aligned} \frac{d}{dy} k= & {} \frac{d}{dy}r^2-y^2=-2y \quad \leadsto \quad \frac{dk}{-2y }=dy,\\ \int y\sqrt{r^2-y^2}\,dy= & {} \int y\sqrt{k}\frac{dk}{-2y}=-\frac{1}{2}\int \sqrt{k}\,dk=-\frac{1}{2}\frac{k^{3/2}}{3/2}\\= & {} -\frac{1}{3}(r^2-y^2)^{3/2}+C. \end{aligned}$$

$\square$

Lemma 8

$$\begin{aligned} \int \arcsin (x/r)\,dx=x\arcsin (x/r)+\sqrt{r^2-x^2}+C \end{aligned}$$

Proof

Starting with integration by parts

$$\begin{aligned} \int u'v=uv - \int uv' \end{aligned}$$

with

$$\begin{aligned} u'=1 \quad \leadsto \quad u=x \quad \hbox {and}\quad v=\arcsin (x/r)\quad \leadsto \quad v'=\frac{d}{dx}\arcsin (x/r) \end{aligned}$$

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)$.

$$\begin{aligned} \frac{d}{dx} \frac{x}{r}&=\frac{d}{dx}\sin (y),\\ \frac{1}{r}&=\frac{d}{dx}y\cos (y),\\ \frac{d}{dx}\arcsin (x/r)=\frac{d}{dx}y&=\frac{1}{r\cos (y)}=\frac{1}{r\sqrt{1-\sin ^2(y)}}\\&=\frac{1}{r\sqrt{1-(x/r)^2}}\mathop {=}\limits ^{r>0}\frac{1}{\sqrt{r^2-x^2}}. \end{aligned}$$

Or with Lemma 6 from before the chain rule gives us

$$\begin{aligned} \frac{d}{dx}\arcsin (x/r)=\frac{1}{\sqrt{1-\frac{x^2}{r^2}}}\cdot \left( \frac{x}{r}\right) '\mathop {=}\limits ^{r>0}\frac{1}{\sqrt{r^2-x^2}}. \end{aligned}$$

Hence we must now solve the integral on the right-hand side of

$$\begin{aligned} \int \arcsin (x/r)\,dx=x\arcsin (x/r)-\int \frac{x}{\sqrt{r^2-x^2}}\,dx. \end{aligned}$$

(9)

Subsituting with $k=r^2-x^2$ thus

$$\begin{aligned} \frac{d}{dx} k= & {} \frac{d}{dx}r^2-x^2=-2x\quad \leadsto \quad \frac{dk}{-2x}=dx,\\ \int \frac{x}{\sqrt{r^2-x^2}}\,dx= & {} \int \frac{x}{\sqrt{k}}\frac{dk}{-2x}=-\frac{1}{2}\int \frac{1}{\sqrt{k}}\,dk=-\frac{1}{2}\int k^{-1/2}\,dk\\= & {} -\frac{1}{2}\frac{k^{1/2}}{1/2}=-\sqrt{k}. \end{aligned}$$

Backsubstitution yields

$$\begin{aligned} \int \frac{x}{\sqrt{r^2-x^2}}\,dx=-\sqrt{r^2-x^2}. \end{aligned}$$

Inserting this into (9) gives the result. $\square$

Lemma 9

$$\begin{aligned} \int \arcsin \left( \sqrt{1-\frac{y^2}{r^2}}\right) \,dy=y\arcsin \left( \sqrt{1-\frac{y^2}{r^2}}\right) -\sqrt{r^2-x^2}+C \end{aligned}$$

Proof

Starting with integration by parts

$$\begin{aligned} \int u'v=uv-\int uv' \end{aligned}$$

with

$$\begin{aligned} u'\!=\!1\quad \!\leadsto \!\quad u\!=\!x\quad \hbox {and} \quad v\!=\!\arcsin \left( \sqrt{1-\frac{y^2}{r^2}}\right) \quad \!\leadsto \quad \! v'\!=\!\frac{d}{dy}\arcsin \left( \sqrt{1\!-\!\frac{y^2}{r^2}}\right) \end{aligned}$$

yields an easier integral. We must first determine $v'$. Using the chain rule we find

$$\begin{aligned} v'&=\frac{1}{\sqrt{1-\sqrt{1-\frac{y^2}{r^2}}^2}}\cdot \left( \sqrt{1-\frac{y^2}{r^2}}\right) '\\&=\frac{1}{\sqrt{1-1+\frac{y^2}{r^2}}}\cdot \frac{1}{2}\left( 1-\frac{y^2}{r^2}\right) ^{-\frac{1}{2}}\cdot \left( \frac{-2}{r^2}y\right) \\&=\frac{1}{y/r}\cdot \frac{-1}{\sqrt{1-\frac{y^2}{r^2}}}\frac{y}{r^2}\\&=\frac{-1}{\sqrt{r^2-y^2}}. \end{aligned}$$

Hence we must now solve the integral on the right-hand side of

$$\begin{aligned} \int \arcsin \left( \sqrt{1-\frac{y^2}{r^2}}\right) \,dy=y\arcsin \left( \sqrt{1-\frac{y^2}{r^2}}\right) -\int y\frac{-1}{\sqrt{r^2-y^2}}\,dy. \end{aligned}$$

(10)

Substituting $k=r^2-x^2$ thus

$$\begin{aligned}&\frac{d}{dx}k=\frac{d}{dx}r^2-x^2=-2x\quad \leadsto \quad \frac{dk}{-2x}=dx,\\&\int \frac{-y}{\sqrt{r^2-y^2}}\,dy=\int \frac{-y}{\sqrt{k}}\frac{1}{-2y}\,dk=\frac{1}{2}\int k^{-1/2}\,dk=\frac{1}{2}2k^{1/2}=\sqrt{k}. \end{aligned}$$

Backsubstitution yields

$$\begin{aligned} \int \frac{-y}{\sqrt{r^2-y^2}}\,dy=\sqrt{r^2-x^2}. \end{aligned}$$

Inserting this in (10) gives the result. $\square$

Lemma 10

$$\begin{aligned} \int \sqrt{r^2-x^2}-b\,dx=\frac{r^2}{2}\arcsin \left( \frac{x}{r}\right) +\frac{x}{2}\sqrt{r^2-x^2}-bx+C \end{aligned}$$

Proof

By linearity we get

$$\begin{aligned} \int \sqrt{r^2-x^2}-b\,dx=\int \sqrt{r^2-x^2}\,dx-b\int 1\,dx=r\int \sqrt{1-x^2/r^2}\,dx-bx. \end{aligned}$$

(11)

And we can continue with the standard trigonometric substitution $x/r=\sin (y)$, thus $y=\arcsin (x/r)$ and we continue with:

$$\begin{aligned}&\frac{d}{dx}\frac{x}{r}=\frac{d}{dx}\sin (y)\quad \leadsto \quad \frac{1}{r}=\frac{d}{dx}y\cos (y)\quad \leadsto \quad dx=dy r\cos (y),\\&r\int \sqrt{1-\sin ^2(y)} r\cos (y)\,dy=r^2\int \sqrt{\cos ^2(y)}\cos (y)\,dy=r^2\int \cos ^2(y)\,dy. \end{aligned}$$

Using the cosine double angle formula

$$\begin{aligned} \cos (2x)=2\cos ^2(x)-1\quad \leadsto \quad \frac{1+\cos (2x)}{2}=\cos ^2(x), \end{aligned}$$

we remove the squared cosine

$$\begin{aligned} r^2\int \frac{1+\cos (2y)}{2}\,dy= & {} \frac{r^2}{2}\int 1+\cos (2y)\,dy=\frac{r^2}{2}y+\frac{r^2}{2}\int \cos (2y)\,dy\nonumber \\= & {} \frac{r^2}{2}y+\frac{r^2}{4}\sin (2y), \end{aligned}$$

(12)

Using the sine double angle formula $\sin (2y)=2\sin (y)\cos (y)$, and Pythagoras’ Theorem to replace the cosine

$$\begin{aligned} \left( \frac{x}{r}\right) ^2&=\sin ^2(y)=1-\cos ^2(y),\\ \frac{r^2-x^2}{r^2}&=\cos ^2(y),\\ \frac{\sqrt{r^2-x^2}}{r}&=\cos (y) \end{aligned}$$

we can rewrite the right-hand side of (12)

$$\begin{aligned} \frac{r^2}{2}y+\frac{r^2}{4}\sin (2y)&=\frac{r^2}{2}y+\frac{r^2}{4}2\sin (y)\cos (y)\\&=\frac{r^2}{2}y+\frac{r^2}{2}\sin (y)\frac{\sqrt{r^2-x^2}}{r}. \end{aligned}$$

Finally reverting the substitution with $y=\arcsin (x/r)$

$$\begin{aligned}&=\frac{r^2}{2}\arcsin (x/r)+\frac{r^2}{2}\sin (\arcsin (x/r))\frac{\sqrt{r^2-x^2}}{r}\\&=\frac{r^2}{2}\arcsin (x/r)+\frac{r^2}{2}\frac{x}{r}\frac{\sqrt{r^2-x^2}}{r}\\&=\frac{r^2}{2}\arcsin (x/r)+\frac{x}{2}\sqrt{r^2-x^2}. \end{aligned}$$

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

$$\begin{aligned} \mathbb {E}\bigl [\hbox {cap}_\mathrm{gauss}(e)\bigr ]=\frac{1}{{\varPhi }(251/50)}\sum _{a=0}^\infty {\varPhi }\left( \frac{250-a}{50}\right) =\frac{501}{2}+\epsilon \end{aligned}$$

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

$$\begin{aligned} 250<\lfloor 50 \cdot \frac{251}{50}\rfloor =\lfloor 251\rfloor =251 \end{aligned}$$

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

$$\begin{aligned} f(x)=\frac{g(x)}{{\varPhi }(251/50)}\qquad g(x)={\left\{ \begin{array}{ll}0&{}x\ge \frac{251}{50}\\ \phi (x)&{}x<\frac{251}{50}\end{array}\right. } \end{aligned}$$

With $\phi (x)$ being the pdf and ${\varPhi }(x)$ denoting the cdf of the standard normal distribution $\fancyscript{N}(0, 1)$

$$\begin{aligned} \phi (x)=\frac{1}{\sqrt{2\pi }}e^{-\frac{1}{2}x^2}\quad {\varPhi }(z)=\int _{-\infty }^z\phi (t)dt\quad \frac{d}{dz}{\varPhi }(z)=\phi (z). \end{aligned}$$

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

$$\begin{aligned} f_X(x)=\frac{1}{\sigma }\phi \left( \frac{x-\mu }{\sigma }\right) \quad F_X(x)={\varPhi }\left( \frac{x-\mu }{\sigma }\right) \end{aligned}$$

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

$$\begin{aligned} \int _{-\infty }^\infty f(x)dx= & {} \frac{1}{{\varPhi }(251/50)}\int _{-\infty }^\infty g(x)dx=\frac{1}{{\varPhi }(251/50)}\int _{-\infty }^{\frac{251}{50}}\phi (x)dx\nonumber \\= & {} \frac{{\varPhi }(251/50)}{{\varPhi }(251/50)}=1 \end{aligned}$$

confirms that we have a probability distribution again.

Now for the expectation

$$\begin{aligned} \mathbb {E}\bigl [250-\lfloor 50\cdot X\rfloor \bigr ]=250-\int _{-\infty }^\infty \lfloor 50\cdot x\rfloor f(x)dx \end{aligned}$$

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.

$$\begin{aligned} 250-\lfloor 50 \cdot x\rfloor&=a\in \mathbb {N}_0\\ 250-a&=\lfloor 50 \cdot x\rfloor \Leftrightarrow x\in \left[ \frac{250-a}{50},\frac{250-a+1}{50}\right) \end{aligned}$$

Thus we integrate

$$\begin{aligned} h(a)=\int _{\frac{250-a}{50}}^{\frac{250-a+1}{50}}f(x)dx= & {} \frac{1}{{\varPhi }(251/50)}\\ \int _{\frac{250-a}{50}}^{\frac{250-a+1}{50}}g(x)dx= & {} \frac{{\varPhi }\left( \frac{250-a+1}{50}\right) -{\varPhi }\left( \frac{250-a}{50}\right) }{{\varPhi }(251/50)} \end{aligned}$$

Checking that we have a distribution

$$\begin{aligned} \sum _{a=0}^{\infty }h(a)=\sum _{a=0}^{\infty }\frac{{\varPhi }\left( \frac{250-a+1}{50}\right) -{\varPhi }\left( \frac{250-a}{50}\right) }{{\varPhi }(251/50)}=\frac{\sum _{a=0}^{\infty }{\varPhi }\left( \frac{250-a+1}{50}\right) -{\varPhi }\left( \frac{250-a}{50}\right) }{{\varPhi }(251/50)} \end{aligned}$$

clearly telescopes to

$$\begin{aligned} \sum _{a=0}^{\infty }h(a)=\frac{{\varPhi }\left( \frac{250-0+1}{50}\right) }{{\varPhi }(251/50)}=1 \end{aligned}$$

Thus the expectation now simply is

$$\begin{aligned} \sum _{a=0}^\infty a\cdot h(a)&=\sum _{a=0}^\infty a\cdot \frac{{\varPhi }\left( \frac{250-a+1}{50}\right) -{\varPhi }\left( \frac{250-a}{50}\right) }{{\varPhi }(251/50)}\\&=\frac{1}{{\varPhi }(251/50)}\left( \sum _{a=0}^\infty a\cdot {\varPhi }\left( \frac{251-a}{50}\right) -\sum _{a=0}^\infty a\cdot {\varPhi }\left( \frac{250-a}{50}\right) \right) \\&=\frac{1}{{\varPhi }(251/50)}\left( \sum _{a=1}^\infty a\cdot {\varPhi }\left( \frac{251-a}{50}\right) -\sum _{a=1}^\infty a\cdot {\varPhi }\left( \frac{250-a}{50}\right) \right) \\&=\frac{1}{{\varPhi }(251/50)}\left( \sum _{a=0}^\infty (a+1)\cdot {\varPhi }\left( \frac{251-(a+1)}{50}\right) \right. \\&\quad \ \left. -\,\sum _{a=1}^\infty a\cdot {\varPhi }\left( \frac{250-a}{50}\right) \right) \\&=\frac{1}{{\varPhi }(251/50)}\left( {\varPhi }\left( \frac{250}{50}\right) +\sum _{a=1}^\infty (a+1)\cdot {\varPhi }\left( \frac{250-a}{50}\right) \right. \\&\quad \ \left. -\,a\cdot {\varPhi }\left( \frac{250-a}{50}\right) \right) \\&=\frac{1}{{\varPhi }(251/50)}\left( {\varPhi }\left( \frac{250}{50}\right) +\sum _{a=1}^\infty {\varPhi }\left( \frac{250-a}{50}\right) \right) \\&=\frac{\sum _{a=0}^\infty {\varPhi }\left( \frac{250-a}{50}\right) }{{\varPhi }(251/50)}\approx 250.500067\ldots \end{aligned}$$

$\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:

$$\begin{aligned} \begin{array}{r|ccc} \mathbb {E}\bigl [\hbox {mf}(\cdot )\bigr ]&{}\hbox {unit}&{}\hbox {unif}&{}\hbox {gauss}\\ \hline llur&{}2\cdot 1&{}\quad 2\cdot \frac{3350}{101}&{}\quad 2\cdot 222.3\\ lmum&{}3\cdot 1&{}\quad 3\cdot \frac{3350}{101}&{}\quad 3\cdot 222.3\\ \end{array} \end{aligned}$$

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:

$$\begin{aligned} \mathbb {E}[\hbox {maxflow(s, t)}]\le \mathbb {E}\bigl [\min \bigl \{\hbox {cap}(s), \hbox {cap}(t)\bigr \}\bigr ]=c\cdot \mathbb {E}\bigl [\min \bigl \{\hbox {cap}(e),\hbox {cap}(e)\bigr \}\bigr ] \end{aligned}$$

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]:

$$\begin{aligned} \mathbb {E}\bigl [X_{(1)}\bigr ]=\mathbb {E}\bigl [\min \{X_1, X_2 \}\bigr ]=\sum _{x=0}^\infty \bigl (1-F_X(x)\bigr )^2. \end{aligned}$$

In case of the discrete uniform distribution $\fancyscript{U}$ on $[0..100]$ we find:

$$\begin{aligned} \mathbb {E}\bigl [U_{(1)}\bigr ]&=\sum _{u=0}^{100}\left( 1-\frac{1+u}{101}\right) ^2=\sum _{u=0}^{100}\left( \frac{100-u}{101}\right) ^2=\sum _{u=1}^{100}\left( \frac{u}{101}\right) ^2=\frac{1}{101^2}\sum _{u=1}^{100}u^2\\&=\frac{1}{101^2}\frac{(200+1)(100+1)100}{6}=\frac{3350}{101}=33.1683. \end{aligned}$$

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

$$\begin{aligned} 250<\left\lfloor 50 \cdot \frac{251}{50}\right\rfloor =251 \end{aligned}$$

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:

$$\begin{aligned} f_X(x)=\frac{g(x)}{{\varPhi }(251/50)}\quad g(x)={\left\{ \begin{array}{ll}0&{}x\ge 251/50\\ \phi (x)&{} x<251/50 \end{array}\right. }, \end{aligned}$$

therefore its cdf is

$$\begin{aligned} F_X(x)={\left\{ \begin{array}{ll} 1 &{} x>251/50\\ {\varPhi }(x)/{\varPhi }(215/50) &{} x\le 251/50 \end{array}\right. }. \end{aligned}$$

Transforming with $Z=250-50X$ yields:

$$\begin{aligned} F_Z(z)={\left\{ \begin{array}{ll}0&{}z<-1\\ 1-{\varPhi }\bigl ((250-z)/50\bigr )/{\varPhi }(251/50)&{}z\ge -1 \end{array}\right. }. \end{aligned}$$

Thus

$$\begin{aligned} \mathbb {E}\bigl [Z_{(1)}\bigr ]&=\sum _{z=0}^\infty \bigl (1-F_Z(z)\bigr )^2=\sum _{z=0}^\infty \left( \frac{{\varPhi }\left( \frac{250-z}{50}\right) }{{\varPhi }(251/50)}\right) ^2\approx 222.3 \end{aligned}$$

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.