On Finding Optimal (Dynamic) Arborescences

Let G = (V, E) be a directed and weighted graph with vertex set V of size n and edge set E of size m, such that each edge (u, v) \in E has a real-valued weight w(u, c). An arborescence in G is a subgraph T = (V, E') such that for a vertex u \in V, the root, there is a unique path in T from u to any other vertex v \in V. The weight of T is the sum of the weights of its edges. In this paper, given G, we are interested in finding an arborescence in G with minimum weight, i.e., an optimal arborescence. Furthermore, when G is subject to changes, namely edge insertions and deletions, we are interested in efficiently maintaining a dynamic arborescence in G. This is a well known problem with applications in several domains such as network design optimization and in phylogenetic inference. In this paper we revisit algorithmic ideas proposed by several authors for this problem, we provide detailed pseudo-code as well as implementation details, and we present experimental results on large scale-free networks and on phylogenetic inference. Our implementation is publicly available at \url{https://gitlab.com/espadas/optimal-arborescences}.


Introduction
The problem of finding an optimal arborescence in directed and weighted graphs is one of the fundamental problems in graph theory with several practical applications.It has been found in modeling broadcasting [1], network design optimization [2], subroutines to approximate other problems, such as the traveling salesman problem [3], and it is also closely related to the Steiner problem [4].Arborescences are also found in multiple clustering problems, from taxonomy to handwriting recognition and image segmentation [5].In phylogenetics, optimal arborescences are useful representations of probable phylogenetic trees [6,7].
Chu and Liu [8], Edmonds [9], and Bock [10] proposed independently a polynomial time algorithm for the static version of this problem.The algorithm by Edmonds relies on a contraction phase followed by an expansion phase.A faster version of Edmonds algorithm was proposed by Tarjan [11], running in O(m log n) time.Camerini et al. [12] corrected the algorithm proposed by Tarjan, namely the expansion procedure.The fastest known algorithm was proposed later by Gabow et al. [13], with improvements in the contraction phase and running in O(n log n + m) time.Fischetti and Toth [14] also address this problem restricted to complete directed graphs, relying on the Edmonds algorithm.The algorithms proposed by Tarjan, Camerini et al., and Gabow et al., rely on elaborated constructions and advanced data structures, namely for efficiently keeping mergeable heaps and disjoint sets.
As stated by Aho et al. [15], "efforts must be made to ensure that promising algorithms discovered by the theory community are implemented, tested and refined to the point where they can be usefully applied in practice."The transference of algorithmic ideas and results from algorithm theory to practical applications can be however considerable, in particular when dealing with elaborated constructions and data structures, a well known challenge in algorithm engineering [16].
Although there are practical implementations of the Edmonds algorithm, such as the implementation by Tofigh and Sjölund [17] or the implementation in NetworkX [18], most of them neglect these elaborated constructions.Even though the Tarjan version is mentioned in the implementation by Tofigh and Sjölund, they state in the source code that arXiv:2311.03262v1[cs.DS] 6 Nov 2023

Optimal arborescences
Both Edmonds and Tarjan algorithms proceed in two phases: a contraction phase followed by an expansion phase.The contraction phase then maintains a set of candidate edges for the optimal arborescence under construction.This set is empty in the beginning.As this phase proceeds, edges selected may form cycles which are contracted to form super-vertices.The contraction phase ends when no contraction is possible and all vertices have been processed.In the expansion phase super-vertices are expended in reverse order of their contraction, and one edge is discarded per cycle to form the arborescence of the original graph.The main difference between both algorithms is on the contraction phase.

Edmonds algorithm
Let G = (V, E) be a directed and weighted graph with vertex set V of size n and edge set E of size m, such that each edge (u, v) ∈ E has a real-valued weight w(u, v).Let each contraction mark the end of an iteration of the algorithm and, for iteration i, let G i denote the graph at that iteration, D i be the set of selected vertices in iteration i, E i be the set of selected edges incident on the selected vertices, and Q i be the cycle formed by the edges in E i , if any.
The algorithm starts with G 0 = G, D i and E i initialized as empty sets, and Q i initialized as an empty graph, for all iterations i.
The algorithm proceeds by selecting vertices in G i which are not yet in D i .If such a vertex exists, then it is added to D i and the minimum weight incident edge on it is added to E i .The algorithm stops if either a cycle is formed in E i or if all vertices of G i are in D i .
If E i holds a cycle, then we add the edges forming the cycle to Q i and we build a new graph G i+1 from G i , where the vertices in the cycle are contracted into a single super-vertex v i+1 .Edges (u, v) ∈ G i are added to G i+1 and updated as follows.

3.
Edges originating from the new vertex: Edges incident to the new vertex: Here w(u, v i+1 ) is the weight of the edge (u, v i+1 ) in G i+1 , w(u, v) is the weight of (u, v) in G i , σ Q i denotes the maximum edge weight in the cycle Q i , and w(u ′ , v) is the weight of the edge in the cycle incident to vertex v.After the weights are updated, the algorithm continues the contraction phase with the next iteration i + 1.
The contraction phase ends when there are no more vertices to be selected in G i , for some iteration i.

Expansion phase
The final content of E i holds an arborescence for the graph G i .Let H denote a subgraph formed by the edges of E i .For every contracted cycle Q i , we add to H all cycle edges except one.If the contracted cycle is a root of H, we discard the cycle edge with the maximum weight.If the contracted cycle if not a root of H, we discard the cycle edge that shares the same destination as an edge currently in H.The algorithm proceeds in reverse with respect to the contraction phase, examining graph G i−1 and the cycle Q i−1 .This process continues until all contractions are undone, and H is the final arborescence.
The contracted version of G 1 , named G 2 , is presented in Figure 3a.Minimum weight incident edges in every vertex of G 2 is marked in red in Figure 3b, and they are added to 3c.This cycle must be also contracted.The maximum weight edge of Q 2 is edge (0, C ′ ) and σ Q 2 = 15.Since E 2 contains a cycle, a contraction is required and we obtain the graph G 3 in Figure 4a.The minimum incident edges in G 3 are marked in red in Figure 4b, and added to   The expansion phase expands the cycles formed in reverse order kicking out one edge per cycle.The removed edges are presented as dashed edges.Let H = E 4 and decrement i. Vertex C ′′′ is a root of H since there is no edge directed towards C ′′′ .In this case every edge of Q 3 is added to H except the maximum weight edge of the cycle as shown in Figure 5a.In iteration i = 2, note that H = {(C ′ , 1)}, and vertex C ′ is a root, therefore every edge from Q 2 must be added to H except the maximum weight edge (C ′ , 0).In the next iteration, i = 1, vertex C is not a root of H since edge (0, C) ∈ H.In this case, we add every edge from Q 1 except the ones that share the destination with edges in H, as illustrated in Figure 5c.Regarding the final expansion, H = {(0, C), (C, 1), (C, 6), (6, 1)} implies that C is not a root.Every edge in Q 0 = {(3, 2), (2,4), (4,5), (5,3)} except edge (3,2) is added to H.The optimal arborescence of G 0 is shown in Figure 5d.

Tarjan algorithm
The algorithm proposed by Tarjan [11] is built on Edmonds algorithm, but it relies on advanced data structures to become more efficient, namely in the contraction phase.The algorithm builds in this phase a subgraph H = (V, E ′ ) of G = (V, E) such that H contains the selected edges.The optimum arborescence could then be extracted from H through a depth-first search, taking into account Lemma 2 in Tarjan paper [11].This lemma states that there is always a simple path in H from any vertex u in a root strongly connected component S to any vertex v in the weakly connected component containing S. Camerini et al. [12] provided however a counter-example for this construction, and they proposed a correction that relies on a auxiliary forest F and that we discuss below.
The algorithm by Tarjan keeps track of weakly and strongly connected components in G, as well as non examined edges entering each strongly connected component.The bookkeeping mechanism used the union-find data structure [24] to maintain disjoint sets.Let SFIND, SUNION, SMAKE-SET denote operations on strongly connected components and WMAKE-SET, WFIND, WUNION on weakly connected components.Find operations allow to find the component where a given vertex lies in, union operations allow to merge two components together, and make-set operations allow to initialize singleton components for each vertex.Non examined edges are kept through a collection of priority queues, implemented as mergeable heaps.Let MELD, EXTRACT-MIN, and INIT denote the operations on heaps, where the meld operation allow to merge two heaps, the extract-min allows to get and remove the minimum weight element, and the initialization operation allows to initialize a heap from a list of elements.We consider also the SADD-WEIGHT operation that allows to add a constant weight to all edges incident on a given strongly connected component in constant time.Note that edges incident on given strongly connected component are maintained in a priority queue where they are compared taking into account its weight and the constant weight added to that strongly connected component.
The correction proposed by Camarini et al. requires us then to maintain a forest F and a set rset that holds the roots of the optimal arborescence, i.e., the vertices without incident edges.Each node of forest F has associated an edge of G, a parent node, and a list of children.

Initialization
Data structures are initialized as follows.queues is an array of heaps, initialized with an heap for each vertex v containing incident edges on v. roots is the list of vertices to be processed, initialized as V.The forest F is initialized as empty as well as the set rset.Four auxiliary arrays are also needed to build F and the optimal arborescence, namely inEdgeNode that for each vertex v stores a node of F associated with the minimum weight edge incident in v, π that stores the leaf nodes of F, cycleEdgeNode that stores for each cycle representative vertex v the list of cycle edge nodes in F, and max that stores for each strongly connected component the target of the maximum weight edge.These arrays are initialized as follows.
roots ← ∅ ▷ Set of vertices to process.

Contraction phase
The contraction phase proceeds while roots ̸ = ∅ as follows.It pops a vertex r from roots and it verifies if there are incident edges in r such that they do not belong to a contracted strongly connected component.If there are such edges, then it extracts the one with minimum weight; otherwise it stops and it continues with another vertex in roots.The pseudo-code is as follows.
Once an incident edge on r is found that does not lie within a strongly connected component, i.e., that is incident on a contracted strongly connected component, we must update forest F. Hence we create a new node minNodeF in forest F associated with edge (u, r).If r is not part of a strongly connected component, i.e., r is not part of a cycle, then minNodeF becomes a leaf of F. Otherwise, F must be updated by making minNodeF a parent of the trees of F that are part of the strongly connected component.The following pseudo-code details this updating of forest F.
Create the node minNodeF in forest F associated with the edge (u, r) The next step is to verify if (u, r) forms a cycle with minimum weight edges formerly selected.It is enough to check if (u, r) connects vertices in the same weakly connected components.Note that (u, r) is incident on a root and, if u lies in the same weakly connected component as r, then adding (u, r) forms necessarily a cycle.Assuming that adding (u, r) does not form a cycle, we perform the union of the sets representing the two weakly connected components to which u and r belong, i.e., WUNION(u, r).We update also the inEdgeNode[r] array as r now has an incident edge selected.
If adding (u, r) forms a cycle a contraction is performed.The contraction procedure starts firstly by finding the edges involved in the cycle, using a backward depth-first search.During this process, a map is initialized where the edge is associated to its F node (the map key).Then the maximum weight edge in the cycle is found, the reduced costs are computed and the weight of the edges is updated.Note that the min-heap property is always maintained when reducing the costs without running any kind of procedure to ensure it, since the constant reduced is added to every edge in a given priority queue.Arrays inEdgeNode and cycleEdgeNode are updated, and heaps involved in the cycle are merged.The pseudo-code is as follows.

Expansion phase
We obtain the optimal arborescence from the forest F, which is decomposed to break the cycles of G.Note that the nodes of F will represent the edges of H seen in Edmonds algorithm.The expansion phase is as follows.We first take care of the super-nodes of F which are roots of the optimal arborescence, represented by the set rset.Each vertex u in rset is the representative element of a cycle, i.e. the destination of the maximum edge of a cycle.Hence u becomes a root of the optimal arborescence, and every edge incident to u in F must be deleted.The tree F is decomposed by deleting the node π[u] and all its ancestors.For the other cycles, which corresponding super-vertices are not optimal arborescence roots, the incident edge (u, v) represented by a root in F is added to H, and the other incident edges represented in F by π [v] and its ancestors are deleted.The procedure ends when there are no more nodes in F. The optimal arborescence is given by H.The pseudo-code is as follows.

Illustrative example
Let us consider the graph G = (V, E) in Figure 6.At the beginning of the contraction phase the forest F is empty.There is a priority queue associated to each vertex and their content is (1,2,10)}, and We have also roots = {0, 1, 2, 3}, rset = ∅, and max We start popping vertices, denoted by r, from the set roots and finding the minimum weighted edge incident to each r.We can safely pop 0, 1 and 2 from roots, and the respective minimum weight incident edges (3, 0), (0, 1), and (3, 2), with weights 1, 6 and 8, respectively, without forming a cycle.These edges are added to forest F as nodes, leading to the state seen in Figure 7. Since each vertex in {0, 1, 2} forms a strongly connected component with a single vertex, we have π[0] = (3, 0), π [1] = (0, 1), and π [2] = (3, 2).
Since R = ∅, the expansion phase proceeds with evaluation of nodes in set N. Set N is processed similarly to set R with two minor changes: the elements of N when removed, are added to H; since N contains edges (u, v) as nodes, then the path P v is traced from the leaf node stored in π [v].This process terminates when N = ∅ and H holds the optimal arborescence.The final arborescence H = {(3, 2), (3, 0), (0, 1)} for our example is depicted in Figure 12.

Optimal dynamic arborescences
Pollatos, Telelis and Zissimopoulos [21,22] proposed two variations of an intermediary tree data structure which is built during the execution of the Edmonds algorithm on G and that is then updated when G changes.We will present the data structure by Pollatos et al. [22], named augmented tree data structure (ATree), that encodes the set of edges H introduced in the previous section, along with all vertices (simple and contracted) processed during the contraction phase of Edmonds algorithm.When G is modified, the ATree is decomposed and processed, yielding a partially contracted graph G ′ = (V ′ , E ′ ).Then the Edmonds algorithm is executed for G ′ .Note that only G and the ATree are kept in memory.
Let us assume that the graph G = (V, E) is strongly connected and that w(u, v) > 0, for all (u, v) ∈ E. If G is not strongly connected, we can add a vertex v ∞ and 2n edges such that w(v

ATree
A simple node of the ATree, represented as N s v encodes an edge with target v ∈ V.A complex-node, represented as N c i...j encodes an edge which target a super-vertex that represents a contraction of the vertices i . . .j ∈ V.In what follows, whenever the type of an ATree node is not known or relevant in the context, we just use N to represent it.The parent of an ATree node is the complex-node which edge targets the super-vertex into which the child edge target is contracted.Since G is strongly connected, all vertices will eventually be contracted into a single super-vertex and the ATree will have a single root.A null edge is encoded in the ATree root node.See Figure 13.The ATree takes O(n) space and its construction can be embedded into the Edmonds algorithm implementation without affecting its complexity.Let us detail how an update in G affects the ATree F, namely edge insertions and deletions.Edge weight updates are easily achieved by deleting the edge and adding it again with the new weight.Vertex deletions are solved by deleting all related edges, and vertex insertions are trivial solved by considering G ′ with the existing super-vertex and the new vertex (and related edges).null The corresponding ATree.Figure 13.A graph and its ATree.The root represents the edge incident to the contraction of all graph nodes (null).

Edge deletion
Let (u, v) ∈ E be the edge we want to delete from G. If (u, v) / ∈ F, we just remove it from G. If (u, v) ∈ F, we remove (u, v) from G and we decompose the ATree: we delete the node N which represents the edge (u, v) and, as we broke the cycles containing (u, v), we also delete every ancestor node of N in F. Each child of a removed node becomes the root of its sub-tree.Then, we create a partially contracted graph G ′ with the remaining nodes in the ATree and we run the Edmonds algorithm for G ′ to rebuild the full ATree F, and find the new optimal arborescence.
The graph G ′ = (V ′ , E ′ ) is obtained from the decomposed ATree as follows.Note that if a complex-node N c i,...,j is a root of F, the super-vertex representing the contraction of vertices i, ..., j belongs to V ′ .Let then {N x 1 , . . ., N x ℓ } be the roots of the ATree F, where x k is the representant of the contraction when N x k is a complex-node in the ATree.Then V ′ = {x 1 , . . ., x ℓ }.E ′ is the set of the incident edges in V ′ .

Edge insertion
Inserting a new edge (u, v) is handled by reducing the problem to an edge deletion.We first add (u, v) in G. Then we check if (u, v) should replace an edge present in ATree F.
Starting from the leaf N s v of the ATree F representing an edge incident to v, and then following its ancestors, we check if w(u, v) is smaller than the weight of the edge represented by each N (see Figure 14).We can replace an edge if the previous condition holds and if N s u is not present in its sub-tree, i.e. (u, v) should not be an edge connecting two nodes of the current cycle.We then engage a virtual deletion of the candidate node (the edge is not deleted but the ATree is decomposed), we build the graph G ′ = (V ′ , E ′ ∪ {(u, v)}) from the decomposed ATree, and we execute the Edmonds algorithm for G ′ to rebuild the full ATree F. (a) null ) Figure 14.We add the edge (2, 0) with weight 2 to the graph ((a) edge represented in red).The process starts with the analysis of node N s 0 , which represents edge (3, 0) with weight 1. (2, 0) cannot replace N s 0 edge as it is heavier.Then N s 0 parent is examined, N c 0,1,3 .The corresponding edge is heavier than (2, 0) and N s 2 is not present in its sub-tree.Then, (2, 0) should replace this node, and N c 0,1,3 and its ancestors are virtually deleted ((b), nodes represented in white).The Edmonds algorithm is executed on the remaining nodes (represented in grey).

ATree data structure
ATree is an extension of the forest F data structure presented in Section 2.2.The nodes of the ATree maintain the following records: the edge of G the node N represents, EDGE(N), the cost of the edge at the time it was selected, w N , its parent PARENT(N), the list of its children CHILDREN(N), its kind (simple or contracted), the list of contracted-edges during the creation of the super-node, the edge of maximum weight in the cycle e max and its weight w max .
In an edge deletion, the edge is removed from the contracted-edges list into which it belongs.In the process of decomposing and reconstructing the ATree, the set of edges E ′ corresponds to the concatenation of the lists contracted-edges associated to each deleted ATree nodes.And we need to update the weight of every edge (u, v) of E ′ .Let N s v be the simple node whose contracted-edges contains (u, v).The new weight w ′ (u, v) is w ′ (u, v) = w(u, v) − ∑ N∈P w N where P is the set of ancestors of N s v , w(u, v) is the original weight, and w N the weight of the edge represented by N at the time it was selected.We run a BFS on each tree to find the subtracted sum r i of each simple node N s i in O(n) time.Then, we scan the edges e to assign the reduced cost w ′ (e) = w(e) − r i .
While adding an edge (u, v), we look for a candidate node to replace in the ATree.The process starts with N s v and it checks every ancestor until the root is inspected or if a node N verifies w ′ (EDGE(N)) > w ′ (u, v), where w ′ (EDGE(N)) denotes the reduced cost of the edge presented by N. If the root is reached without verifying the condition, we insert (u, v) in the contracted-edges list of the lowest common ancestor of N s v and N s u .If the condition is met, we found a candidate node N where (u, v) could be added, but we must determine if (u, v) is safe to be added.We check if N s u is already present with a BFS in the sub-tree of root N. If we find N s u , we insert (u, v) in the contracted-edges list of the lowest common ancestor of N s v and N s u .Otherwise, we engage a virtual deletion of EDGE(N) (the edge is not deleted but the ATree is decomposed), then we build the graph G ′ = (V ′ , E ′ ∪ {(u, v)}) and execute the Edmonds algorithm for G ′ as mentioned before.The pseudo-code for finding a candidate is as follows, where (u, v) is edge to be inserted.nodeF = N s v mechanisms proposed by Tarjan, and the forest data structure introduced by Camerini et al., and further extended as the ATree data structure.Recall that the bookkeeping mechanism adjusted to maintain the forest data structure relies on the following data structures: for every node v, a list L(v) stores each edge incident to v; disjoint sets keep track of the strongly and weakly connected components; a collection of queues keeps track of the edges entering each vertex; and a forest or, in the dynamic case, an ATree F.

Incidence lists
Since edges of G are processed by incidence and not by origin, G is represented as an array of edges sorted with respect to target vertices.This is beneficial since it takes advantage of memory locality bringing improvements to the overall performance.

Disjoint sets
Two implementations of the union-find data structure for managing disjoint sets are used, with both supporting the standard operations.One is used to represent weakly connected components, while the other is employed for strongly connected components.The latter is an augmented version.In the case of the first implementation, the following common operations are supported: WFIND(x) that returns a pointer to the representative element of the unique set containing x; WUNION(x) which unites the sets that contain x and y; and WMAKE-SET(x) that creates a new set whose only element and representative is x.For the augmented implementation, the same operations are supported, but named SFIND, SUNION and SMAKE-SET; two extra operations are also supported, namely SADD-WEIGHT and SFIND-WEIGHT detailed below.
Our implementations of union-find data structure rely on the conventional heuristics, namely union by rank and path compression, achieving nearly constant time per operation in practice; m operations over n elements take O(mα(n)) amortized time, where α is the inverse of the Ackermann [25].Both implementations use two arrays of integers, namely the rank and the parent array instead of pointer based trees.Even though operations computational complexity is theoretically speaking the same, using arrays instead of pointers promotes again memory locality since arrays are allocated in contiguous memory.
The purpose of having a different implementation for strongly connected components is to bring a constant time solution for the computation of the reduced costs, exploiting the path compression and union by rank heuristics.While finding the minimum weight incident edges in every vertex in the contraction phase, cycles may arise.Then the maximum weight edge in the cycle is found, the reduced costs are computed, and the weight of incident edges is updated by summing the reduced costs.In this context, the augmented version of the union-find data structure supports then the following operations as mentioned above: SADD-WEIGHT(x, k) which adds a constant k to the weight of all elements of the set containing x, and SFIND-WEIGHT(x) that returns the accumulated weight for the set containing x. Supporting these operations requires an additional attribute weight, represented internally as a third array to store the weights.The weights are initialized with 0. The SADD-WEIGHT(x, k) operation adds value k to the root or representative element of the set containing x in constant time.The operation SFIND(x) has been rewritten for updating the weights whenever the underlying union-find tree structure changes due to the path compression heuristic; this change does not change the complexity of this operation.The operation SFIND-WEIGHT(x) performs the sum of all values stored in field weight on the path from x until we meet the root of disjoint-set containing x; the cost of this operation is identical to the cost of operation SFIND.A constant time solution is obtained then for updating the weight of all elements in a given set, which allows us to update the weight of all edges incident on a given vertex also in constant time.

Queues
Heaps are used to implement the priority queues which track the edges incident in each vertex.In this context, three types of heaps were implemented and tested, namely binary heaps [26], binomial heaps [27] and pairing heaps [28,29].The pairing heaps is the alternative that has simultaneously better theoretical and expected experimental results; although binary heaps are faster than all other heap implementations when the decreasekey operation is not needed, pairing heaps are often faster than d-ary heaps (like binary heaps) and almost always faster than other pointer-based heaps [30].Our experimental results consider also this comparison (see Section 5).With respect to theoretical results, using pairing heaps to implement priority queues, and assuming that n is the size of a heap, the common heap operations are as follows: INIT(L) creates a heap with elements in list L in O(n) time; INSERT(h, e) inserts an element e in the heap h in Θ(1) time; GET-MIN(h) obtains the element with minimum weight in Θ(1) time; EXTRACT-MIN(h) returns and removes from the heap h the element with minimum weight in O(log n) amortized time; DECREASE-KEY(h, e) decreases the weight of element e in o(log n) amortized time; and MELD(h 1 , h 2 ) merges two heaps h 1 and h 2 in Θ(1) time.
Our implementation does not rely on the DECREASE-KEY operation, but it relies heavily on the MELD and EXTRACT-MIN operations.In this context it is important to note that the MELD operation takes O(n) time for binary heaps and O(log n) time for binomial heaps.The EXTRACT-MIN(h) runs in O(log n) time for both binary and binomial heaps.On the other hand both pairing and binomial heaps are pointer based data structures, while binary heaps are array based.Hence, it is not clear a priori which heap implementation would be better in practice and, hence, it is a topic of analysis in our experimental evaluation as mentioned.

Forest
Several data structures were introduced to manage F and the cycles in G.A set rset holds the roots of the optimal arborescence, i.e. the vertices which do not have any incident edge.A table max holds the destination of the maximum edges in a strongly connected component.A table π points to the leaves of F, where π[v] = (u, v) means that the node (u, v) of F was created during the evaluation of vertex v.The table inEdgeNode holds for each v, the unique node of F entering the strongly connected component represented by v. Finally, the list cycleEdgeNode holds the lists of nodes in a cycle, where cycleEdgeNode[rep] holds the nodes of the cycle represented by rep.
These data structures allow us to construct and maintain the forest F within the contraction phase without burdening the overall complexity of the algorithm.They allow also to extract an optimal arborescence in linear time during the expansion phase.Detailed pseudo-code has been presented in Section 2.2.
This representation is extended for implementing the ATree taking into account the data structure description and the pseudo-code presented in Section 3.4.

Complexity
Let us discuss the complexity of our implementation for finding a (static) optimal arborescence in a graph G with n vertices and m edges.
In the initialization phase we mainly have the n INIT operations for the priority queues, the n SMAKE-SET operations on augmented disjoint sets, the n WMAKE-SET operations on disjoint sets, and O(n) operations on other data structures.All these operations take constant time each, thus the initialization takes O(n) time.
In the contraction phase, only the operations on priority queues and disjoint sets may not take constant time.The operations on priority queues are at most m EXTRACT-MIN operations and n MELD operations.Since EXTRACT-MIN takes O(log n) time and (for pairing heaps) the MELD operation takes constant time, then it takes O(m log n) total time for maintaining priority queues.The operations on disjoint-sets are m WFIND and SFIND operations, n WUNION and SUNION operations, and n SADD-WEIGHT operations.The disjoint set operations take O(mα(n)) total time where α is the inverse of the Ackermann function [25].The other operations run in O(m + n) time.Therefore, the contraction phase takes O(m log n) time.
In the expansion phase, F contains no more than 2n − 2 nodes and each node of F is visited exactly once, so the procedure takes O(n) time.The total time required to find an optimal arborescence is therefore dominated by the priority queue operations yielding a final time complexity of O(m log n).
Let us analyse now the cost of maintaining dynamically the optimal arborescence.Let ρ be the set of affected vertices and edges, |ρ| the number of affected vertices, and ||ρ|| the number of affected edges.A vertex is affected if it is included in a different contraction in the new output after an edge insertion or removal.Note that |ρ| < n, that all operations in an addition or deletion of an edge occur in O(n) time and that a re-execution of Edmonds algorithm processes only the affected vertices.The update of an optimal arborescence, using the implementation presented in Section 2.2, can then be achieved in O(n + ||ρ|| log |ρ|) time per edge insertion or removal.

Experimental evaluation
We implemented the original Edmonds algorithm as described in Section 2.1, and Tarjan algorithm as described in Section 2.2.The implementation of Tarjan algorithm has three variants which differ only on the heap implementation.As discussed before, we considered binary heaps, binomial heaps and pairing heaps in our experiments.Algorithms were implemented in Java 11, and binaries were compiled with javac 11.0.20.Experiments were performed on a computer with the following hardware: Intel(R) Xeon(R) Silver 4214 CPU @ 2.20GHz and 16 GB of RAM.
The aim of this experimental evaluation is to compare the performance of Edmonds original algorithm with Tarjan algorithm, to evaluate the use of different heap implementations, and to investigate the practicality of the dynamic algorithm for dense and sparse graphs.As datasets we used randomly generated graphs, both dense and sparse, and real phylogenetic data.

Datasets
Graphs datasets are comprised by sparse and dense graphs generated accordingly to well known random models.To generate sparse graphs we considered three different models.One of them was the Erdos-Rényi (ER) model [31], with p = c log n n , c ≥ 1, where p denotes the probability of linking a node u with a node v and n is the number of nodes in the network.Whenever p has the previously defined value, the network has one giant component and some isolated nodes.Moreover, these graphs were generated using fast_gnp_random_graph generator of the NetworkX library [32], with p = 0.02.
Sparse scale-free directed graphs were also generated using the model by Bollobás et al. [33] (identified as scale-free in our experiments) and a variant of the duplication model by Chung et al. [34].The first were generated using the scale_free_graph function of the NetworkX library, with all parameters set with their omission value except the number of nodes.The later were generated using our own implementation, where given 0 ≤ p ≤ 1, the partial duplication model builds a graph G = (V, E) by partial duplication as follows: start with a single vertex at time t = 1 and, at time t > 1, perform a duplication step: uniformly select a random vertex u of G; add a new vertex v and edges (u, v) and (v, u) with (different) random weights; for each neighbor w of u, add edges (v, w) and/or (w, v) with probability p, and random integer weights chosen uniformly from [0, 1000].
Dense graphs were generated using the complete_graph generator of the NetworkX library, that creates a complete graph, i.e. all pairs of distinct nodes have an edge connecting them.Edge weights were assigned randomly.
Running time and memory is averaged over five runs and for five different graphs of each size, for all models.
We used also real phylogenetic data in the dynamic updating evaluation, namely real dense graphs using phylogenetic datasets available on EnteroBase [35], respective details are shown in Table 1.Graphs were built based on the pairwise distance among genetic pro-files, as usual in distance based phylogenetic inference [7].The experiments on these data were carried out by considering increasing volumes of data, namely [10%, 20%, 30%, ..., 100%].

Edmonds versus Tarjan
We compared both Edmonds and Tarjan algorithms for complete and sparse graphs using generated graph datasets.This comparison is presented in Figure 15.As expected, Tarjan algorithm is faster and the experimental running time follows the expected theoretical bound of O(m log n).The memory requirements are also lower for the Tarjan algorithm, growing linearly with the size of the graph, as expected.
Given these results, we omit Edmonds algorithm from the remaining evaluation.

Different heap implementations
Results for scale-free graphs are presented in Figure 16.The running time and memory requirements are according to expectations and to the analysis for complete and sparse graphs in the previous section.The somewhat strange behavior in memory plots for a lower number of vertices and edges is due to Java's garbage collector and it can be ignored.
The focus in this section is the performance of different heap implementations together with Tarjan algorithm.The improved theoretical performance of binomial and pairing heaps is not supported by our experiments and in fact fared no better than binary heaps.Pairing heaps obtained a similar time performance to binary heaps in the duplication models while simultaneously using less space.This is particularly interesting since the meld operation is more efficient for pairing heaps.However the memory locality exploited by binary heaps plays here an important role.

Dynamic optimal arborescences
Let us compare the performance of maintaining dynamic optimal arborescences versus ab initio computation on edge updates.Both implementations rely on the algorithm by Tarjan described in this paper.Our experiments consist on evaluating the running time and required memory for adding and deleting edges.Results are averaged over a sequence of 10 independent DELETE operations, and also over a sequence of 10 independent ADD operations.The sequences of edges subject to deletion or insertion were randomly selected.
Figure 17a provides the results for the DELETE operations.We observe that updating the arborescence is twice as fast compared with its ab initio computation.Note that these results are aligned with the results presented by Pollatos, Telelis and Zissimopoulos [21,22].Figures 18a to 22a provide the results for the DELETE operation over phylogenetic data described above.As the size of the dataset grows, and the inferred graph becomes larger, the dynamic updating becomes also more competitive, being twice as fast when compared with the ab initio computation.
The results for the ADD operations are presented in Figure 17b for complete graphs, and in Figures 18b to 22b.It is clearly perceived that the ab initio computation is outmatched by the dynamic updating, in particular as the size of the graph grows.The dynamic updating is consistently at least twice as fast as the ab initio computation, surpassing often that speedup factor.We evaluated also the memory requirements for dynamic updating.We only measured the memory consumption for the DELETE operation because the ADD operation is essentially reduced to an edge removal operation.Tables 2 to 5 show the memory usage comparison between the ab initio computation and the dynamic updating, averaged over 10 operations.In each table, the first column contains the % of the dataset being considered, the second column presents the memory usage for the ab initio computation, the third column presents the memory usage for the dynamic updating, and the fourth column presents the memory ratio between dynamic updating and ab initio computation.As an illustrative baseline, Figure 23 shows the memory usage for Yersinia.wgMLSTdataset as an increasing percentage of it added to the computation.Given these results, we can observe that both ab initio computation and dynamic updating require linear space on the size of the input.This is also consistent with the results for random graphs presented above.However the dynamic updating requires three times more memory on average than the ab initio computation, which is expected given that a more complex data structure needs to be managed.

Conclusions
We provided implementations of Edmonds algorithm and of Tarjan algorithm for determining optimal arborescences on directed and weighted graphs.Our implementation of Tarjan Algorithm incorporates the corrections by Camerini et al., and it runs in (m log n) time, where n is the number of vertices of the graph and m is the number o edges.We provide also an implementation for the dynamic updating of optimal arborescences based on the ideas by Pollatos, Telelis and Zissimopoulos, and that relies on Tarjan algorithm, running in O(n + ||ρ|| log |ρ|) per update operation and scaling linearly with respect to memory usage.We highlight the fact that our implementations are generic in the sense that a generic comparator is given as parameter and, hence, we are not restricted to weighted graphs; we can find the optimal arborescence on any graph equipped with a total order on the set of edges.To our knowledge, our implementation for optimal arborescence problem on dynamic graph is the first one to be publicly available.The code is available at https://gitlab.com/espadas/optimal-arborescences.
Experiment evaluation shows that our implementations comply with the expected theoretical bounds.Moreover, while multiple changes occur in G, the dynamic updating is at least twice as faster as the ab initio computation, requiring although more memory even if by a constant factor.Our experimental results corroborate also the results presented by Böther et al. and Pollatos et al.We found one shortcoming regarding the dynamic optimal arborescence, namely the high dependence between the time needed to recalculate the optimum arborescence and the affected level of the ATree.The lower the level, the larger the number of affected constituents will be.A prospect to achieve a more efficient dynamic algorithm could be relying on link-cut trees [36] which maintains a collection of node-disjoint forests of self-adjusting binary heaps (splay-trees [37]) under a sequence of LINK and CUT operations.Both operation take O(log n) time in worst-case.
With respect to the application in the phylogenetic inference context, we highlight the fact that the proposed implementation for dynamic updates allows to significantly improve the time required to update phylogenetic patterns as datasets grow in size.We note also that, due to the use of heuristics in the probable optimal tree inference, there are some algorithms that include a final step for further local optimizations [6].Although it may not be always the case, it seems that we can often incorporate such local optimization as total order over edges.Given that our implementations assume that such a total order is given as parameter, such optimizations can be incorporated.The challenge of combining these techniques to implement classes of local optimizations is also a path for future work.
Minimum weight incident edges in every vertex of graph G 0 , colored in red.
Contracted version of graph G 0 named G 1 .Minimum weight incident edges in every vertex of graph G 1 .
Cycle in G 2 colored in green.

Figure 3 .
Contraction of cycle Q 1 and identification of cycle Q 2 .
required, leading to G 4 with a single vertex, Figure4d.In this last iteration we have E 4 = ∅ and Q 4 = ∅, ending the contraction phase.

Figure 4 .
Figure 4. Contraction of Q 2 , identification of Q 3 , and final graph G 4 after Q 3 contraction.
Expansion subgraph H in iteration i = 3.
Expansion subgraph H in iteration i = 2.

Figure 5 .
Expansion subgraph H in iteration i = 1.Expansion subgraph H in iteration i = 0. Expansion phase and optimal arborescence.

Figure 8 .
Figure 8.First contraction of the input graph.

Figure 9 .
Figure 9. Adding directed edges from node (2, 1) to the nodes of cycle C in forest F.

Figure 10 .
Figure 10.Last cycle and final contracted graph.
(a) Path P 3 colored in red.
Running time for sparse graphs.
Memory for complete graphs.
Memory for sparse graphs.

Figure 15 .
Comparison between the Edmonds algorithm and the Tarjan algorithm (three different heap implementation) on complete and sparse graphs.

Figure 16 .
Running time for scale-free graphs.Running time for duplication model graphs.Memory for scale-free graphs.Memory for duplication model graphs.Comparison if three different heap implementations within the Tarjan algorithm sparse scale-free graphs.

Figure 17 .Figure 18 .
Optimal arborescence updating versus ab initio computation for DELETE and ADD operations on complete graphs.Running time averaged over 10 random operations.(a) DELETE operations.(b) ADD operations.Optimal arborescence updating versus ab initio computation for DELETE and ADD operations on clostridium.Griffiths dataset.Running time averaged over 10 random operations.(a) DELETE operations.(b) ADD operations.

Figure 19 .
Optimal arborescence updating versus ab initio computation for DELETE and ADD operations on Moraxella.Achtman7GeneMLST dataset.Running time averaged over 10 random operations.

Figure 20 .
Optimal arborescence updating versus ab initio computation for DELETE and ADD operations on Salmonella.Achtman7GeneMLST dataset.Running time averaged over 10 random operations.(a) DELETE operations.(b) ADD operations.

Figure 21 .Figure 22 .
Optimal arborescence updating versus ab initio computation for DELETE and ADD operations on Yersinia.cgMLSTv1dataset.Running time averaged over 10 random operations.(a) DELETE operations.(b) ADD operations.Optimal arborescence updating versus ab initio computation for DELETE and ADD operations on Yersinia.wgMLSTdataset.Running time averaged over 10 random operations.

Table 1 .
Phylogenetic datasets.The first three without missing data.The number of vertices n is the number of genetic profiles in each dataset.

Table 2 .
Memory usage comparison for the dynamic updating and ab initio computation of an optimal arborescence for the clostridium.Griffiths dataset.

Table 3 .
Memory usage comparison for the dynamic updating and ab initio computation of an optimal arborescence for the Moraxella.Achtman7GeneMLST dataset.

Table 4 .
Memory usage comparison for the dynamic updating and ab initio computation of an optimal arborescence for the Yersinia.cgMLSTv1dataset.

Table 5 .
Memory usage comparison for the dynamic updating and ab initio computation of an optimal arborescence for the Yersinia.wgMLSTdataset.

Table 6 .
Memory usage comparison for the dynamic updating and ab initio computation of an optimal arborescence for the Salmonella.Achtman7GeneMLST dataset.