Partial Implementation of Max Flow and Min Cost Flow in Almost-Linear Time

Abraham and Neiman give two algorithms in [1] for create-petal which both take in a graph, a starting point $x_{0},$ a target $t$ and a radius $r$ (where $r$ is expressed as a range in the first algorithm) and which both return a cluster (which is the petal) as well as a center of that petal. The algorithm given on page 235 is not fast enough for the overall desired bound in this paper: it requires that the ratio between the largest and smallest distances is $O(n^{k})$ for some $k$ which does not necessarily hold in general, and the create-petal algorithm has a for loop which could potentially run for many iterations, and there is no clear bound on how many iterations this would take. For this reason, we implement the faster algorithm on page 245 in section 6 of [1], which the authors themselves state is necessary for the best asymptotic efficiency.

In the code, we begin with an ordinary implementation of Dijkstra’s Algorithm using Fibonacci Heaps along with a corresponding shortest_path function. Everything in this portion is standard. After that, we write a function to generate the directed graph $\tilde{G}$ as specified on page 245 of [1]. As Abraham and Neiman point out, all edges on the directed graph have nonnegative length.

With that out of the way, we can write the create_petal function. In essence, this algorithm iteratively carves petals off of the graph and then carves petals off of the remaining part of the graph until all petals have been carved, when any remaining vertices are made part of the stigma of the flower. This function takes in the following inputs: a graph $Y$, a target $t$, a starting point $x_{0}$ and a radius $r.$ Note that this is slightly different from the slower function that Abraham and Neiman propose. In particular, the create_petal function does not take in the original graph as input except at the beginning; rather, as described in the following subsection, each time a petal is carved, create_petal is iteratively called on the remaining subgraph. We also do not have any hi and lo parameters, as the faster algorithm described on page 245 of [1] does not include any.

To create a petal, we first run Dijkstra on the given graph with the given starting point $x_{0}.$ We then generate the corresponding directed graph based on the shortest distances from $x_{0}$ to all other points in the graph and calculate the shortest path to the target $t.$ Next, we need to find the special vertex $r^{\prime}$, defined on page 245 of [1] as the furthest point from $t$ on the shortest path from $t$ to $x_{0}$ in the original graph $G$ that is within distance $r$ of $t.$ To do this, we could perform binary search on the sorted distances on the shortest path from $t$ to $x_{0}$ in order to find the largest number less than or equal to $r.$ However, in our implementation we simply do a linear search as this does not impact the overall asymptotic time complexity, which we empirically check later.

Next, the algorithm stipulates that the petal is the ball of radius $r/2$ around $t$ in the directed graph where each edge on the path from $t$ to $p_{r^{\prime}}$ is set to half of its original weight in $G.$ At first, this may seem to be equivalent to the following: include every vertex within the ball of radius $r/2$ and then also include every vertex on the path from $t$ to $p_{r^{\prime}}$ including $p_{r^{\prime}}.$ This would make intuitive sense because $p_{r^{\prime}}$ is defined as the furthest point along the path from $t$ to $x_{0}$ within distance $r$ of $t$. If all of the closest points to $t$ are further than $r$ away, then $p_{r^{\prime}}=t.$ An implementation which does this would not be quite right for the following reason: while edges in the directed graph all have nonnegative length by the triangle inequality, they can have length 0, and this fact influences the implementation of create_petal. Therefore, in the directed graph, $p_{r^{\prime}}$ may not be the furthest point within distance $r.$ Edges of length 0 mean that some points “further” along the path may be of equal distance.

Therefore, it is essential to first set the weight of each edge on the path from $t$ to $p_{r^{\prime}}$ to half its original weight in $G.$ After that, we run Dijkstra on the directed graph $\tilde{G}$ and select all vertices within distance $r/2$ of $t.$

The function then returns the petal and the point $p_{r^{\prime}}$ which is the center of the petal and plays a key role in the petal decomposition. In addition, we also choose to return the shortest path from $x$ to $t$, since this will come up in the petal_decomposition function described in the next subsection.

In our implementation of this function, we begin with some sanity checks that the center $x$ and target $t$ are both in the graph. While this may seem unnecessary, we will see in the next subsection that there are many recursive calls which involve modifying the graph. These assert and raise statements helped catch many errors earlier on in the implementation and ensured everything worked properly once those errors were corrected.

We begin by computing the radius of graph, defined as the maximum length across all vertices $v$ of the shortest path from $x_{0}$ to $v$ found using Dijkstra. We then initialize the following lists with their defined purposes:

• •

  Ys: $Y_{j}$ is defined as the original graph minus the vertices and edges contained in petals $X_{1},\ldots,X_{j}$. In particular, $Y_{0}$ is initialized to the entire graph before any petals have been removed.

• •

  Xs: $X_{j}$ is the petal found on iteration $j$ and is carved from $Y_{j-1}$, except for $X_{0}$ which is the remaining part of the graph after all petals have been carved (what the authors in [1] call the “stigma” of the flower)

• •

  xs: $x_{j}$ is the center of petal $X_{j},$ or the point $p_{r^{\prime}}$ when computing the petal $X_{j}$, except for $x_{0}$ which is the original input

• •

  ts: $t_{j}$ is the target when creating petal $X_{j}$

• •

  ys: $y_{j}$ is the neigbhor of $x_{j}$ on the path from $x_{0}$ to $t_{j}$ which is closer to $x_{0}$

The rest of the implementation follows the pseudocode. The first petal that is carved may be a special petal if the distance from $x_{0}$ to the target $t$ is greater than $\frac{5}{8}$ of the radius from $x_{0}.$ Abraham and Neiman describe why they choose the $\frac{5}{8}$ constant in [1]. The purpose of this petal is to preserve the shortest path from $x_{0}$ to $t.$

To carve all of the other petals, we iteratively check for points whose distance from $x_{0}$ is at least $\frac{3r}{4}.$ If no other points exist, we have found all petals. Otherwise, we arbitrarily choose one of those points to be the target $t_{j}$ for petal $X_{j}.$ In our implementation, we simply choose the first key in the dictionary of key value pairs where the keys are vertices and the values are distances from the center $x_{0}$, but this is arbitrary because the keys are only ordered as they were in the original input graph, which is arbitrary.

We then call create_petal on the remaining part of the original graph (denoted by the element most recently appended to the Ys list) along with our center $x_{0}$ and target $t_{j}.$ For the radius, we use $r/8$. This choice comes from the pseudocode for petal_decomposition in [1], where the call to create_petal sets lo = 0 and hi = $r/8.$ However, as we recall from the previous subsection, the faster implementation of create_petal does not have any lo or hi parameters but rather just a single radius parameter. For this reason, it was not obvious whether to use $r/8$ or something like $r/16.$ However, we decided to use $r/8$ because this parameter is a radius and all points within this radius will have distance less than or equal to $r/8$ from the center. For this reason, having a max possible value of $r/8$ (the hi parameter) vs a single value of $r/8$ seemed more analogous. However, that is not a fully rigorous argument, and the testing described in Section 4.1.4 is what provides empirical validation for this choice.

Once we create the petal, we subtract it off from the remaining part of the graph. Using the path from $x_{0}$ to $t_{j}$, we compute $y_{j}$, which will be used in the following subsection. We then halve all of the edges on the shortest path to $t$ (but beginning at $x_{j}$ rather than $x_{0}$) returned by the create_petal function as given in the pseudocode.

Finally, we let $X_{0}$ be the last $Y_{j}$, namely the remaining part of the graph after all petals have been carved.

This is the third and final function necessary to generate a low-stretch spanning tree. For inputs, it takes in the graph, a starting point $x_{0}$ and a target $t.$ This function is recursive with a base case that the tree on a graph with one vertex is just the single vertex itself.

On each recursive call, we compute the petal decomposition of the graph. As explained in the previous subsection, this returns a set of petals $X_{j}$, their centers $x_{j}$, a set of neighbors $y_{j}$ and a set of targets $t_{j}$. We then run hierarchical-petal-decomposition recursively on each petal $X_{j}$ where the starting point is $x_{j}$ and the target is $t_{j}.$ Each of these recursive calls generates a tree $T_{j}$ for each petal $X_{j}.$ These trees are all disjoint, so we create a new graph which combines all of these $T_{j}$ through the edges $(x_{j},y_{j}).$ This gives us a final low-stretch spanning tree $T.$

Lastly, to compute the tree itself, we call hierarchical-petal-decomposition on the graph with arbitrary starting point $x_{0}$ but critically also with target $x_{0}.$ The justification for this decision is explained in [1].

The function hierarchical_petal_decomposition was the last function described in [1]. In order to test the LSST implementation, we generate many Erdős-Rényi random graphs. For each random graph, we pick a probability $p$ and let there be an edge between every pair of distinct vertices with probability $p$ (so $p=1$ would correspond to a complete graph).

Then, for each graph, we compute an LSST using hierarchical_petal_decomposition, measuring its stretch and the duration to compute it. We repeat this for 10 trials and take the average for both the stretch and the duration. There are two primary things to verify: that the stretch is $O(m\log n\log\log n)$ and that the time complexity is $O(m\log n\log\log n).$ We also separately checked that the output was always a spanning tree (it reaches every vertex of the original graph and has $n-1$ edges). However, when timing the code, we did not include this extra check since it is not a necessary part of the algorithm once we are confident that the algorithm is running correctly.

To check the asymptotic nature of the stretch and duration to compute the LSST, for different values of $m,n$ we graph the stretch and duration vs $m\log n\log\log n.$ We expect that there are some finite constants $c_{1},c_{2}\in\mathbb{R}$ such that the stretch with $m$ edges and $n$ vertices is bounded above by $c_{1}m\log n\log\log n$ and likewise the duration is bounded above by $c_{2}m\log n\log\log n.$ We check this graphically below.

In both cases, we can see that for sufficiently large values of $m$ and $n,$ both the duration and stretch are bounded above by a linear constant multiplied by $m\log n\log\log n.$ Of course, this does not definitively verify that the implementation is correct (for example, there could hypothetically be a $\log\log\log n$ factor that would not show up graphically in the experiments run above) but given that the authors in [1] prove their algorithm and stretch are both $O(m\log n\log\log n)$ and the implementation followed their algorithm, it is reasonable to conclude that the graphs above provide empirical validation that our implementation works correctly.

The heavy_light_decomposition function takes in one input: the tree. Within this function, we define two inner functions: dfs_size and decompose_chain, which we describe below.

For the dfs_size function which takes in a vertex $v$ and its parent as inputs. Then, the function uses Depth First Search to recursively visit every vertex and compute the size of the subtree (number of vertices) rooted at that vertex. If the vertex is a leaf node, we say its subtree has size 1 (just itself). Otherwise, we initialize the size of the subtree to 1 and compute $\operatorname{size}(v)=1+\sum_{w\in c(v)}\operatorname{size}(w)$ where $c(v)$ is the set of all children of $v$.

Next, we have the decompose_chain function. This function takes in a vertex $v$, the head of the chain of $v$ and the parent of $v$ as inputs. In this function, we first identify the heavy child for that vertex, which corresponds to condition 2 in Definition 4.5 and is the one child of $u$ within $P_{i}$, the heavy chain containing $u$. We find this heavy child by iterating over the children of $u$ and picking the child of $u$ which has the maximum subtree size.

The decompose_chain function then recursively calls itself as long as we keep finding new “heavy children.” For each recursive call, we take in the previous heavy_child as our new vertex $v,$ the same chain_head as before and the vertex $v$ as the parent (since the heavy_child is a child of vertex $v$).

Once the recursive calls are complete, we iterate over the children of $v$, and for every child that isn’t the heavy child, we initialize a new chain and recursively call decompose_chain where our new vertex $v$ is the child we’re iterating over, the chain_head is the index of the current chain and the parent is our old $v.$ By setting the code up in this way, we ensure that we create a chain for every vertex $v$ and its descendants and that we satisfy the requirements of Lemma 4.5.

With these two functions implemented, we can finish our heavy_light_decomposition implementation. We initialize variables for the length of the tree, the subtree sizes, and the chains. We then call dfs_size on the root vertex $0$ and its “parent” (which is not a real vertex) but which we denote $-1.$ This allows us to figure out the sizes of all of the subtrees in the tree. Then, we call decompose_chain on our root vertex $0$, the chain_head which we initialize to $0$ and the parent $-1$. Finally, the heavy_light decomposition function returns the list of chains.

In order to check our implementation, we check the three conditions of Lemma 4.5 empirically. We produce a function that generates random trees. These trees are implemented as a list of lists, with nodes numbered $0$ to $n$ where $0$ is the root. Further, the list at index $i$ consists of the children of node $i$, and the list is empty if and only if node $i$ is a leaf node. Based on this format, such random trees can be generated as follows: for vertex $n$, randomly choose one of the vertices $i=0,1,\ldots,n-1$ to be its parent by appending vertex $n$ to the list at index $i$. This also has the convenient property that the parent’s number is always strictly less than the child’s number.

We then call our heavy light decomposition function on 1000 randomly generated trees. We iterate over all possible vertices for that tree and check that vertex $u$ is represented in one and only one heavy chain. We can then easily check the children of $u$ and ensure that at most one child of $u$ is in that same heavy chain.

For the third condition, we calculate $u^{\uparrow T}$ for each $u$ and count the number of heavy chains that intersect with it. This is the trickiest condition to check because Lemma 4.5 only promises an $O(\log n)$ bound on the number of intersecting heavy chains. To check this, we consider increasing values of $n$ and graph the average number of intersecting heavy chains vs $\log n$, where the average is calculated by summing up the number of intersections across all vertices and then dividing by both the number of vertices and the number of randomly generated trees, in this case 1000. Below, we plot the average number of intersections vs $\log n$, and we see that a line with slope 1.15 dominates the growth. Thus, empirically it appears that the number of intersections is bounded by $1.15\log n.$ In addition, we also computed the quotient of the number of intersections divided by $\log n$, which gave us the following list: [1.276, 1.196, 1.144, 1.134, 1.126, 1.103, 1.092, 1.090, 1.084, 1.081, 1.079, 1.072]. We can see that this list is monotonically decreasing, which provides further assurance that the number of intersections is $O(\log n).$

In addition, we also empirically check the time complexity for our implementation just like we did for low-stretch spanning trees. Lemma 4.4 states that the algorithm should be $O(n)$ time, which indeed appears to be the case below.