Speeding up iterative applications of the BUILD supertree algorithm

Context for Build in the current OToL pipeline

The current pipeline divides the full supertree problem into subproblems as described in Redelings & Holder (2017). These subproblems cover regions of the taxonomy tree in such a way that adjacent sub-problems contain common nodes (higher taxa) but do not contain common edges. Each subproblem also contains the relevant region of any input tree that coincides with the region of the taxonomy tree. Repeated calls to Build produce a supertree solution to each subproblem. Solving the subproblems is a time-consuming step for the entire pipeline. Increasing the speed of the subproblem-solver would speed up the pipeline, but would also allow it to handle larger subproblems. Larger subproblems occur when OToL curators add more input phylogenies to the pipeline and this large set of inputs conflict with more groupings found in the taxonomy. Larger subproblems also occur when the final location of an incertae sedis taxon is far from its initial location in the taxonomy (see Redelings & Holder, 2019 for a discussion of the handling of incertae sedis taxa). The increased speed might eliminate the need for decomposing the supertree problem into subproblems.

Overview of the Build algorithm

The Build algorithm is a recursive algorithm that determines if a set of rooted triplets of the form x, y|•z are jointly compatible (Aho et al., 1981). Here, we use the symbol • to stand for the root of the tree to emphasize the rooting of the triple. If they are compatible it constructs a tree that displays all the triplets.

The algorithm works by creating a graph for each recursive level. If the graph forms a single connected component, then the Build algorithm has detected an incompatibility between the set of triplets. The graph contains a node corresponding to each leaf that is relevant to the current recursion level. At the root-most level of the recursion every leaf is relevant. A rooted triple is relevant to a recursion level only if each of the three leaves is relevant at that level. For each relevant triple x, y|•z an edge is added between x and y. Since the leaves in the more closely related pair {x, y} are sometimes called the “cluster” of the triple, the resulting graph is called the “cluster graph” for the given set of triplets.

If no triplets are incompatible on a level, then the procedure applies Build to each connected component of the cluster graph. This is the recursive aspect of the algorithm. The nodes that form a connected component at one level will constitute the relevant nodes for the next recursive call. Thus the relevant leaf set is partitioned at each level. However, the edges between nodes are not passed to the next recursive level. Instead they are constructed from the set of relevant triplets at each level. A rooted triple that is relevant at one level may not be relevant at a subsequent recursion level. Thus, not all rooted triplets have a corresponding edge at the next level. A set of triplets is jointly compatible only if the Build algorithm finds no incompatibility at any recursive level.

Rooted splits

The Build algorithm was originally designed to work on sets of rooted triplets. However, when working with evolutionary trees, it is more convenient to work with sets of rooted splits because each branch corresponds to a rooted split. A rooted split σ on taxon set $\mathcal{T}$ is written σ = σ₁|•σ₂, where σ₁ and σ₂ are non-overlapping subsets of $\mathcal{T}$. We refer to σ₁ as the include group or “cluster”, and σ₂ as the exclude group. The root of the tree is on the side of the exclude group. Note that σ₁∪σ₂ may be smaller than $\mathcal{T}$.

It is straightforward to extend Build to operate on a set of rooted splits by treating each rooted split as shorthand for all the rooted triplets that it implies. So Build is commonly modified to take a set of rooted splits Σ instead of a set of rooted triplets.

The relationship of the current optimizations to previous work

Previous approaches to speed up Build have tried to improve the order of computation for the analysis at each recursive level. Deng & Fernández-Baca (2018) decreased the order to $\mathcal{O}\left(M{log}^{2}M\right)$ for a set of phylogenies with M = (number of edges) + (number of leaves). This is accomplished by decreasing the time for the analysis on each recursive level from $\mathcal{O}\left(M^2\right)$ to $\mathcal{O}\left({log}^{2}M\right)$. It also involves sharing work between different levels of the recursive algorithm by retaining a graph between successive levels instead of creating a new graph from scratch on each level.

Our approach differs from previous work because we attempt to share work between successive calls to Build. When calling Build with one set of splits followed by calling Build again with one additional split, we seek to reuse work from the first call while performing the second. This requires saving work from the first call in an object that represents the solution to Build, and passing that solution object as input to the second call. We do not attempt to decrease the order of the computation as a function of the number of leaves or edges. Instead we base our incremental approach on a naive algorithm that has total cost $\mathcal{O}\left(M^3\right)$. It may be possible to create an algorithm that is both incremental and has a more favorable order of computation for large M, but we do not attempt that here.

Data structures used to explain implementations of Build

In order to explain our incremental Build algorithm in an understandable manner, we seek to introduce the full complexity of the algorithm in stages.¹ We thus begin by describing a version of the traditional Build algorithm that saves its work in a solution data structure. This will allow us to focus on changes to the algorithm instead of changes to the data structures when we introduce the first incremental version of the algorithm at a later stage. Our version of the Build algorithm constructs the connected components of the cluster graph without explicitly constructing the cluster graph. That is, our algorithm does not directly represent the edges of the cluster graph in memory.

We begin by introducing the Solution and Component record types that we will use to store temporary work during the algorithm. Using the terminology common to object-oriented programming, we will refer to an instance of these record types residing in a computer’s memory as an “object.” The following sections will describe the distinct steps in Build using some common programming notational conventions: (1) parentheses after the name of an algorithm denote the arguments supplied to that algorithm and (2) the period (or “dot”) notation after the name of a data object is used to refer to a field within that record.

An initial Solution object can be created before applying the Build algorithm. Upon termination, the result of the algorithm will be emitted by storing it in that Solution object. As mentioned above in the overview of the algorithm, each level of Build’s recursion starts with a set of taxa and a set of relevant splits. Connected components of taxa are created and merged as each input split is considered. In our object-oriented description of the algorithm, this corresponds to creation of a Solution object for the current level of recursion. That Solution object will hold a set of Component objects –each of which will store information about the connected components created during the algorithm. Each Component represents a sub-problem that needs to be solved. If no incompatibility is detected at the current level, then the algorithm recursively calls Build to create a Solution object for each Component.

Thus the Solution object forms a tree by indirect recursion to mirror the recursive structure of Build: each Solution object can contain some Component objects, and each Component object contains a Solution object (see Section B.2 of the Supplementary Information File). We refer to the tree associated with a Solution object as a “solution data structure”.

The solution and component objects for a level of Build

The definition of the Solution and Component record types is given in Fig. 1. Creating a new Solution object is done with a procedure referred to as CreateBlankSolution(T), where T is a collection of taxon identifiers. A new solution object, S, has only its S.T field initialized; That field holds a copy of the taxon identifiers; thus initialization has computational complexity $\mathcal{O}\left(\left|T\right|\right)$.

Defining the Build algorithm

To start the algorithm, an initial Solution object can be created and initialized by filling its S.T field with the complete taxon set, $\mathcal{T}$. All other fields of the solution object are initially empty. We can think of the entire Build algorithm as taking a taxon set T and set of rooted splits, Σ. The full Build algorithm consists of the following steps:

1. initialize a Solution object, S to contain the taxon set T. We will refer to this set of operation as a function: CreateBlankSolution(T);

2. call a helper algorithm BuildA (S, Σ); and

3. return the result stored in object S.

The BuildA algorithm describes the set of operations to be performed at a single level of recursion as well as how to call the next levels of recursion.

Dissecting the operations required at each level of BuildA

The operations performed in BuildA (S, Σ) can be separated in several steps:

1. RemoveIrrelevantSplits(S, Σ) to remove splits from Σ that are not relevant at the current level. Splits that are removed here as irrelevant are also implied by the split S.T| •(T-S.I) and are recorded on the solution in field S.I.

2. MergeComponents(S, Σ) to find the connected components of the cluster graph by merging any two components that overlap the include group (“cluster”) of a split in Σ.

3. Fail(S). Return Failure if there is only one connected component.

4. AssignSplitsToComponents(S, Σ) to the single component on the next recursive level where they may be relevant.

5. Iterate. For each non-trivial component C ∈S.$\mathcal{C}$:

   1. use CreateBlankSolution(C.T) to create a new solution object for each Component object and store this in C.S.

   2. Recurse. Call BuildA (C.S, C.ΔΣ) on the sub-problem for each non-trivial component C, Return Failure if that call fails.

The RemoveIrrelevantSplits, MergeComponents, and AssignSplitsToComponents steps all have the form of a for-loop that iterates over splits in Σ and modifies S. This fact will be used later, when we seek to incrementalize these steps. It suggests a naive incrementalization strategy of simply iterating over any newly added splits ΔΣ to add their effects to S.

Simple optimization: trivial and non-trivial components

In a solution object S, we store an array S.M of pointers to components. We also implement S.T as an array, so that if t = T[i] then taxon t belongs to component S.M [i].

Trivial components are defined as components that contain only one taxon. An optimization to improve memory usage and speed is to simply use a Null reference in the S.M array to indicate that the corresponding taxon is in a trivial component, rather than creating a component object for every trivial component. This requires only minor tweaks to the MergeComponents step to create a new Component object on-the-fly whenever a set of trivial components are merged into each other with no non-trivial component involved in the merger.

Example #1: creating a solution

Consider calling Build $\left(\mathcal{T},\Sigma \right)$ with taxa $\mathcal{T}=\left\{A_1,A_2,B\right\}$ and splits Σ = {A₁A₂|•B}. This creates a solution data structure shown in Fig. 2. Build begins by creating a Solution object S₁ where S₁.T =[A₁, A₂, B] and then calling BuildA (S₁, Σ).

The BuildA creates a single non-trivial component C₁ containing {A₁, A₂} and a single trivial component containing B. BuildA does not fail, because there is more than one component. The single split A₁A₂|•B is assigned to C₁. BuildA then creates a new solution C₁.S where C₁.S.T = [A₁, A₂] and then calls BuildA (C₁.S, C₁.ΔΣ).

The second-level BuildA removes the split A₁A₂|•B from C₁.ΔΣ and adds it to C₁.S.I. It then calls MergeComponents, Fail, and AssignSplitsToComponents with no effects. No non-trivial components are created, and the call to BuildA succeeds because it contains the trivial components {A₁} and {A₂}. Since there are no non-trivial components, another round of recursion is not performed.

Stepwise construction of BuildInc

The goal of the incrementalized algorithm is to reuse previous work from computing Build $\left(\mathcal{T},\Sigma \right)$ when computing Build $\left(\mathcal{T},\Sigma +\Delta \Sigma \right)$. The incrementalized algorithm has the signature BuildInc (S, ΔΣ). Although the taxon set $\mathcal{T}$ and the previously added splits Σ are not explicitly present in the signature of BuildInc, they are contained within S. S may also contain intermediate results from a previous BuildInc call, so that BuildInc can reuse previous work. If BuildInc (S, ΔΣ) succeeds, then S is modified in-place to contain the new splits ΔΣ in addition to Σ. However, if BuildInc fails, then S is unmodified.

In order to simplify the description of the incrementalized algorithm BuildInc, we first introduce two partially-incrementalized algorithms BuildInc′ and BuildInc″. This allows us to simplify the explanation of BuildInc by introducing concepts in a step-wise fashion. The BuildInc′ and BuildInc″ algorithms both mutate the solution object even when the algorithm returns False. The fully incrementalized algorithm BuildInc adds the ability to track changes that are made to the solution object, and then roll them back if the algorithm ultimately returns False. We now describe key difference in the steps of the first incrementalized algorithm BuildInc′. We will use the ′and ″symbols to decorate the names of algorithms based on whether they are used in BuildInc′ or BuildInc″.

The incrementalized algorithm allows adding input splits either one-at-a-time, or in larger batches. We assume that the input splits have been partitioned into a series of non-overlapping subsets ΔΣ_i according to some strategy that is chosen in advance. Each subset constitutes a batch. The simplest strategy is to set ΔΣ_i = {σ_i} for some split σ_i, so that splits are added one-at-a-time. However, it is possible to incorporate input splits into the solution object in larger batches, and we explore this option in the Results section below.

Given a partition of the input splits, the general structure for all of the incrementalized algorithms is:

1. initialize a Solution object, S, to contain the taxon set T. We will refer to this set of operation as a function: CreateBlankSolution(T);

2. BuildIncA (S, 0̸)

3. for each subset ΔΣ_i

   1. call BuildIncA (S, ΔΣ_i)

4. return the result stored in object S.

In the descriptions that follow, Σ will be used to denote the set of splits that have been added in previous increments of BuildInc′ without failing. Thus, Σ + ΔΣ would be the set of splits included if the current round of BuildInc′ succeeds, and ΔΣ is the set of splits that are new to the current increment of BuildInc′.

Example #2: adding a split to an unmodified component

When incrementally adding a split that resolves a node below the top level, BuildIncA′ passes the split down the tree until it reaches the level of the resolved node. At levels higher than the resolved node, the split is assigned to one of the components, but does not modify it. At the resolved node, the split either modifies components or causes the creation of a new non-trivial component.

To see this, let us consider incrementally adding the split σ₂ = a₁a₂|•a₃ to a Solution object that contains σ₁ = a₁a₂a₃|•b (Fig. 3). Here the split σ₂ resolves a node in the original solution data structure.

Running BuildIncA′ (S₁, {σ₂} leaves the only non-trivial Component C₁ at the top level unchanged. σ₂ is Assign ed to C₁ and ends up in C₁.ΔΣ.

This leads to a call to BuildIncA′ (S₂, {σ₂}). σ₂ is not removed, so S₂.I is unchanged. However, σ₂ leads to the creation of a new non-trivial component C₂ in Merge. σ₂ is Assign ed to C₂ and ends up in C₂.Σ.

Finally, we get a call to BuildIncA′ (S₃, {σ₃}). Here σ₂ is finally removed, and placed into S₃.I. No non-trivial components are created, so we are done.

Example #3: merging two components

When incrementally adding a split that groups two children of a node in the solution tree, the components that correspond to the grouped children tend to re-emerge at a lower level. However, BuildIncA′ is not able to recognize this, and so must solve each of these re-emergent components from scratch.

To see this, let us consider incrementally adding the split a₁b₁|•c to a Solution object that contains a₁a₂|•b₁ and b₁b₂|•c (Fig. 4). The new split merges the original components C₁ (containing and {A₁, A₂}) and C₂ (containing {B₁, B₂}). The new component {A₁, A₂, B₁, B₂} is then stored in the modified Component object for C₁.

The BuildIncA handles merged components by extracting the splits from subsumed solutions. Therefore, C₁.ΔΣ contains all three splits, and not just the incrementally added split. We can see that the original component {A₁, A₂} re-emerges at a lower level as C₃. However, BuildIncA′ solves it from scratch. The component {B₁, B₂} does not re-emerge at a lower level. This is because the split b₁b₂|•c can be satisfied on the same level, and is not passed down.

BuildIncA″: MergeComponents ″(S, ΔΣ, O) step

We seek to construct the connected components of the cluster graph for the splits Σ(S)∪ΔΣ∪Σ(O) from the connected components for Σ(S). In addition to iterating over ΔΣ and merging components, BuildIncA″ must additionally handle splits in Σ(O). One possible approach would be to iterate over splits in Σ(O) and merge components that overlap each split. However, it turns out that there is a more efficient way.

We note that each original solution S′ ∈ O is a connected component of the original cluster graph for some higher recursive level. Additionally, the taxa in S′.T are connected solely by edges corresponding to splits in Σ(S′) (see Section ??). This is because any split whose include group overlaps the component is assigned to (and only relevant to) that component. Therefore, the effect of the splits Σ(S′) is only to connect the taxa in S′.T. We may therefore iterate over original solutions S′ ∈ O and merge any two components that both overlap S′.T. This makes it unnecessary to extract the splits σ ∈ S′.

BuildIncA″: AssignSplitsToComponents ″(S, ΔΣ, O) step

The BuildIncA″ needs to assign splits in Σ(S), ΔΣ, and Σ(O) to their corresponding components. Each split in ΔΣ is assigned to a component C and stored in C.ΔΣ, just as as in BuildA and BuildIncA′. However, there are two major differences from BuildIncA′.

First, BuildIncA″ does not need to do anything for splits in Σ(S), because they are already assigned to their correct component. BuildIncA′ needed to extract splits from C.O and add them to ΔΣ if C was a merged component. However, BuildIncA′ passes C.O down to the next recursive level directly, so this is unnecessary. Note that since BuildIncA″ no longer assigns splits from Σ(S) to C.ΔΣ, C.ΔΣ will consist only of splits from ΔΣ.

Second, BuildIncA′ must assign splits in Σ(O) that were recieved from the previous recursive level to their corresponding component. However, all splits for a subsolution S′ ∈O have to end up in the same component. This is because S′.T must be entirely contained in one component. The splits of S′ must then be in the same component. Therefore, we may simply assign original solutions S′ to components in their entirety. We can determine which component a solution S′ goes to by looking up the component for any element in S′.T.

BuildIncA″: recursive call to BuildIncA″

As before, BuildIncA″ (S, ΔΣ, O) must iterate over the non-trivial components C $\in S.\mathcal{C}$ and construct a solution C.S for each of them. In all cases this is done by calling BuildIncA″ (C.S, C.ΔΣ, C.O). Passing original solutions from modified components down the tree in C.O allows us to reuse previous work, as described in the section on RemoveIrrelevantSplits ′ below. The only difference from BuildIncA′ is the third argument, the set of original solutions, O.

BuildIncA″: RemoveIrrelevantSplits ″(S, ΔΣ, O) step

There are two differences between BuildIncA′ (S, ΔΣ) and BuildIncA″ (S, ΔΣ, O). First, O might contain a solution S′ to the problem we are trying to solve. In that case, we would like to reuse that solution. Second, when checking for implied splits, we need to check inside each original solution S′ ∈ O as well as inside ΔΣ. We now discuss these differences in greater detail.

Reusing previous work in BuildIncA″ (S, ΔΣ, O)

Suppose that O contains a solution S′ with the same taxon set as S. In that case, we would like to re-use the previous work in S′ instead of solving the problem from scratch.

We first note that Σ(S) must be empty, since O was non-empty, so it is safe to discard S and use S′ instead. Second, O is empty after removing S′. This is because S′ contains all the taxa in S.T and there are no remaining taxa in S.T for any other solutions to contain.

We therefore replace the call to BuildIncA″ (S, ΔΣ, O) with BuildIncA″ (S′, ΔΣ, 0̸). This preserves the invariant that S and O cannot both be non-empty.

Removing implied splits from original solutions

Recall that a split σ should be placed into S.I if and only if σ₂ does not intersect S.T. For an original solution S′ ∈ O, we claim that only splits in S′.I can satisfy this condition. If a split σ from S′ is not in S′.I, then σ₂ must intersect S′.T. But S′.T is a subset of S.T, so σ₂ intersects S.T as well and cannot go into S.I. Therefore, for each original solution S′, we only need to check the splits in S′.I.

If any split in S′.I meets the criterion, then we remove it from S′.I (conceptually) and move it to S.I. However, if one of the splits is removed from a solution, then the solution no longer has the property that its splits are all in the same connected component. It is no longer a solution. In such a case, we say that the solution is “punctured”.

In order to retain our invariants, we remove punctured solutions from O. However, we must retain the splits that they contain. For each punctured original solution S′, we copy the splits in S′.I that were not moved to S.I into the set ΔΣ. We move the child solutions S′. ${\mathcal{C}}_i$.S into O. In this way, all the splits of S′ are retained—some in S.I, some in ΔΣ, and some in O.

Original solutions in O that are not punctured may be retained unmodified.

Note that replacing an original solution S′ ∈O with its child solutions retains the invariants that solutions in O (i) have non-overlapping taxon sets and (ii) are contained within S.T. The taxon sets of child solutions to a solution S′ are contained within S′.T and are non-overlapping. Since S′.T is contained within S.T and does not overlap with any other solutions in O, its child solutions must also be contained within S.T and not overlap any other solutions in O.

Example #4: merging two components and re-using original solutions

When incrementally adding a split that groups two children of a node, the components that correspond to the grouped children tend to re-emerge at a lower level. While BuildIncA′ solves the solutions from scratch, BuildIncA″ is able to re-use the original solutions to these components.

To see this, we show how BuildIncA″ solves the same problem that BuildIncA′ solved in Example #3 (Fig. 4). Unlike BuildIncA′, BuildIncA″ leads to C₁.ΔΣ containing only the single incrementally-added split A₁B₁|•C. Instead of extracting the splits A₁A₂|•B₁ and B₁B₂|•C from the original solutions S₂ and S₃, BuildIncA′ passes down the original solutions themselves.

The original component {A₁, A₂} re-emerges at a lowever level, and is solved by re-using the original solution S₂. The other original solution S ₃ is punctured, and its split B₁B₂|•A is added to S₄.I.

An overview of the rollback approach

All modifications to each Solution object S are complete before any modifications are made to its children. We can therefore represent modifications to the solution data structure as a sequence R of modifications to individual Solution objects. If BuildIncA returns failure at the top level, we can then walk this sequence in reverse order, and roll back the changes that it describes. We create the record type RollbackInfo to record modifications made to an individual Solution object (Fig. 5).

Component mergers are the most interesting type of modifications that BuildIncA′ makes to Solution objects. We can represent component mergers for each Solution object as a sequence of individual mergers of two components. For each component merger, we record enough information about the two original components to reverse the merger. It is then possible to reverse the mergers by walking the sequence in reverse order, and undoing each merger in turn. We create the record type MergeRollbackInfo to record modifications made to an individual Solution object (Fig. 5).

Mergers are of two types: mergers of two non-trivial components, or mergers of a non-trivial component with a trivial component. We handle mergers of two trivial components, by creating an empty “non-trivial component” data structure, and merging it with each trivial component in turn.

Details of rollback info

Optimizations

Recording and replaying rollback info allows us to avoid discarding saved work. However, both recording and replaying rollback info also have a cost. In order to achieve the optimum speedup from rollback info, we must avoid paying this cost when we do not need to.

Modifications to BuildIncA

In order to record rollback info, we must make a few modifications to BuildIncA. In RemoveIrrelevantSplits, we record the original number of implied splits. In MergeComponents we record (i) the original number of components, (ii) a merge-record for each component merger, and (iii) a copy of $\mathcal{C}$ after new components are added, but before empty components are removed.

Rollback (S,r)

After running BuildIncA, BuildInc must run RollbackAll (S, R) in case of failure. This consists of running RollbackOne (S, r) on individual RollbackInfo objects r.

1. Truncate S.I to its previous length r.n_orig_implied_splits.

2. If r.n_orig_components is equal to 0, then clear S. $\mathcal{C}$ and set S.M [t] = Null for each taxon t.

3. If r.n_orig_components is more than 0, walk the list of merge records in reverse order, and undo each one.

4. If we recorded the original components vector S₀.$\mathcal{C}$ then

   1. swap(S.$\mathcal{C}$, S₀.$\mathcal{C}$).

   2. Truncate S.$\mathcal{C}$ to r.n_orig_components.

This process seems simple enough, but one aspect of it that is tricky. Some components are created during merging that (i) are not original components, and also (ii) end up being empty. So they are not final components. We need these components to survive (i.e not be deallocated) so that we can temporarily add elements to them during rollback. We will then move these elements out of them into original components.

Batching

Recall that our supertree algorithm works by considering an ordered list of trees. We seek to construct the set of splits from these trees that are jointly compatible.

The batching approach works by batch-adding the splits for the each tree in order. To batch-add a group of splits, we try and add the whole group of splits. If that succeeds, then we keep the whole group. If it fails, and the group has one split then we are done. If it fails, and the group has more than one split, then we batch-add the first half of the group, followed by the second half of the group. This ends up being both simpler and more efficient than trying fixed-size batches because we do not need to figure out the ideal batch size and it has exponential back-off.

Batching improves efficiency when Σ is large because the cost of determining the compatibility of Σ + ΔΣ does not increase much as the size of ΔΣ increases. If all the splits in ΔΣ will be accepted, it is thus substantially more efficient to add them in one batch.

Oracle

The oracle first runs conflict analysis on each input phylogeny T to identify branches of T that conflict with the tree of currently accepted splits. We then batch-add splits corresponding to the non-conflicting branches of T. This makes batching more efficient by making it more likely that large batches do not contain any conflicting splits.

Unfortunately, the oracle cannot filter splits for the taxonomy tree if there are any incertae sedis taxa. This is because taxonomy branches may correspond to partial splits in the presence of incertae sedis taxa. Our current conflict analysis does not handle partial splits.

Simulation experiments

The simulations script (gen_subproblem.py) can be found in the otcetera repository on GitHub (https://github.com/OpenTreeOfLife/otcetera and via DOI: https://zenodo.org/doi/10.5281/zenodo.10041275). The user of the script specifies: (a) a number of leaves in the full tree, (b) a number of phylogenetic input trees to simulate, (c) a tip inclusion probability for each phylogenetic input, (d) a number of edge-contract-refine (ECR) moves to conduct on each input tree, and (e) an edge contraction probability for the taxonomy. For each replicate, the script uses DendroPy (Sukumaran & Holder, 2010) to generate a pure birth (Yule) tree with the specified number of leaves as the true(model) full tree for that replicate. The specified number of phylogenetic inputs are created by sub-sampling the full tree (using the tip-inclusion probability to assess whether a tip remains in the sampled input); if all tips are deleted a new phylogenetic input is drawn, rather than emitting an empty tree. Then the specified number of ECR moves are applied to the tree to mimic phylogenetic estimation noise. The last tree emitted for each replicate is designed to mimic the taxonomic input in the problems used by the Open Tree of Life project. The taxonomic input is complete (lacks any sub-sampling based on taxon inclusion probabilities). In addition to having errors introduced by ECR moves, the taxonomic input undergoes branch collapsing (using the user-supplied edge-contraction probability on each internal edge independently) to mimic the unresolved character of most taxonomies.

We simulated a collection of supertree problems containing 50-1000 taxa, all with 20 phylo-inputs.

Each OTU was included in the phylogeny inputs with probability 0.5. Each edge in the taxonomy was collapsed with probability 0.75. We introduced two ECR errors per tree. For each simulation condition, we determined the run time by averaging across 15 simulated data sets.

Simulated data

We first examined the effect of the batching and oracle optimizations on simulated data sets as the number of taxa increased (Fig. 6). Batching yields a speedup that increases from 1.6-fold at 50 taxa to 17-fold at 1,000 taxa. Using the oracle to eliminate inconsistent splits shows no speedup when not paired with batching.

However, when combined with batching, the oracle yields an additional speedup that increases from 1.3-fold at 50 taxa to 2.4-fold at 1,000 taxa. This indicates that the oracle allows larger batch sizes to succeed.

Given that our simulation protocol performs two ECR edits to each phylogeny, we expect about four splits to be inconsistent with the underlying tree per phylo-input. Therefore larger trees have a smaller fraction of inconsistent splits. That may explain why larger trees recieve a bigger speedup from batching.

Incrementalizing Build yields a larger speedup than the oracle+batch optimization (Fig. 7). The speedup increases from 20-fold at 50 taxa to 398-fold at 1000 taxa. This represents an additional 9.9-fold speedup over batching + oracle at 1000 taxa. Much of this speedup relies on the abilty to save work via rollback: the incremental algorithm achieves only a 116-fold speedup at 1,000 taxa. Rollback ability thus provides a speedup of 3.4-fold for BuildInc over BuildInc″.

It is possible to combine the batching and oracle optimizations with BuildInc. This is because BuildInc allows ΔΣ to include a batch of splits instead of just adding one split at a time. Adding batching to BuildInc yields a slight speedup of 27% over BuildInc alone at 1,000 taxa. However adding batching + oracle yields a 2.1-fold slowdown.

OpenTree data

We also examined a data set taken from the OToL synthesis release 13.4. As mentioned above, the OToL project normally divides the full supertree problem into subproblems after taxonomy-only taxa are removed (Redelings & Holder, 2017). Here we consider the effect of optimizations on running the full supertree problem without dividing it. This is a much larger scale than the simulated data sets, which are designed to be similar to a single subproblem.

The data set includes 1,223 non-empty trees in addition to the taxonomy. The taxonomy tree has 94,028 leaves. To give an idea of the size of the input trees, the three largest input trees contain 11,217, 8,369, and 7,160 leaves respectively. All but 80 trees contain fewer than 331 leaves.

The total run time without optimization is 3,323 min = 55 h 23 min (Table 1). Batching + oracle decreases the runtime to 345 min, which is a 9.6-fold decrease. BuildInc decreases the runtime further to 5 min 44s. This is a 60-fold speedup over batching + oracle, and a 579-fold speedup overall. Surprisingly, if we disable rollback, the runtime is 527 min, which is even slower than batching + oracle. This indicates that the OToL supertree problem contains more conflicting splits than the simulated data set above. Additionally, BuildInc + batching + oracle achieves an additional 1.95-fold speedup over BuildInc. This is different than in the simulations above, where BuildInc + batching + oracle was slower than BuildInc.

When using the naive Build algorithm, our supertree algorithm considers a total of 183,850 splits and calls Build 183,850 times. Of these splits, 160,682 are accepted. So, 87.4% of splits are accepted and 12.6% are rejected.

When the oracle is enabled, 162,517 splits are considered. This lowers the fraction of rejected splits to 1.1%. When using the oracle + batch optimizations, Build is called only 14,481 times.

Figure 8 shows that some phylo-inputs take a lot more time than others. Processing splits from the taxonomy tree takes a large fraction of the total time. However, when graphed against the number of accepted splits, the relationship is more linear, indicating that the effect of large trees is at least partly driven by the fact that large trees have a large number of splits.

Incrementalizing faster version of Build

Our current approach uses a naive approach to Build that achieves $\mathcal{O}\left(M^3\right)$ time. Future research should consider whether it is possible to modify our incremental algorithm to incorporate recent improvements to Build that decrease the order, such as Deng & Fernández-Baca (2018). In fact, their algorithm relies on another incremental algorithm –incremental graph connectivity–which identifies new connected components that appear as edges are deleted or added (Holm, de Lichtenberg & Thorup, 2001). One approach might therefore be to replace our map of taxa to Component objects at each level with a graph that supports incremental connectivity queries. This would involve incrementally adding any required edges at each recursive level and identifying connected components that are merged.

One difficulty with this strategy is that we assume the data structures used to find connected components on each recursive level are separate. In constrast, Deng & Fernández-Baca (2018) construct a single graph at the beginning of the algorithm, and then remove additional edges from the single, shared graph at each recursive level. It might be possible to adapt this approach by saving a copy of the relevant portion of the graph at each recursive level to preserve a snapshot of previous work that is unmodified by work at deeper levels. However, even if this is feasible, it is not yet clear how much this would change the running time of the algorithm.