VeryFastTree: speeding up the estimation of phylogenies for large alignments through parallelization and vectorization strategies

VeryFastTree: speeding up the estimation of phylogenies for large alignments through parallelization and vectorization strategies C´esar Pi˜neiro, Jos´e M. Abu´ın, and Juan C. Pichel Supplementary Material 1 Additional Experimental Results 1 2 4 8 10 12 Threads 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 Time (hours) Small (single precision) FT2 (SSE) VFT (SSE) VFT (AVX) 1 2 4 8 10 12 Threads 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 Time (hours) Medium (single precision) FT2 (SSE) VFT (SSE) VFT (AVX) 1 2 4 8 10 12 Threads 4 6 8 10 12 14 16 18 Time (hours) Large (single precision) FT2 (SSE) VFT (SSE) VFT (AVX) 1 2 4 8 10 12 Threads 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Time (hours) Small (double precision) FT2 (SSE not supported) VFT (SSE) VFT (AVX) 1 2 4 8 10 12 Threads 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 Time (hours) Medium (double precision) FT2 (SSE not supported) VFT (SSE) VFT (AVX) 1 2 4 8 10 12 Threads 4 5 10 15 20 25 30 35 Time (hours) Large (double precision) FT2 (SSE not supported) VFT (SSE) VFT (AVX) Figure S1: Running times of FastTree-2 and VeryFastTree using different number of threads, vector extensions and arithmetic precision (graphs in log-log scale). 1.1 Running times Experiments were conducted using a server with a 12-core Intel Xeon E5-2680v3 processor and 128 GB of memory. We have considered three datasets in our tests labeled as small, medium and large, according to their sizes. The small1dataset contains the alignment of 1,000 sequences from the ARB database including Eucarya, Bacteria and Archaea. The medium1dataset contains the alignment of 25,057 protobacteria, and the large2dataset includes 331,550 sequences of 16S ribosomal RNAs. Measurements were carried out using the default options (JTT model together with the CAT approximation). In addition, we used the 1https://cme.h-its.org/exelixis/web/software/raxml/index.html 2http://www.microbesonline.org/fasttree/ S1 Small (single precision) 1 2 4 8 10 12 Threads 0 0.5 1 1.5 2 2.5 3 Speedup FT2 (SSE) VFT (SSE) VFT (AVX) Medium (single precision) 1 2 4 8 10 12 Threads 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Speedup FT2 (SSE) VFT (SSE) VFT (AVX) Large (single precision) 1 2 4 8 10 12 Threads 0 0.5 1 1.5 2 2.5 3 3.5 4 Speedup FT2 (SSE) VFT (SSE) VFT (AVX) Small (double precision) 1 2 4 8 10 12 Threads 0 1 2 3 4 5 6 Speedup FT2 (SSE not supported) VFT (SSE) VFT (AVX) Medium (double precision) 1 2 4 8 10 12 Threads 0 1 2 3 4 5 6 Speedup FT2 (SSE not supported) VFT (SSE) VFT (AVX) Large (double precision) 1 2 4 8 10 12 Threads 0 1 2 3 4 5 6 7 8 Speedup FT2 (SSE not supported) VFT (SSE) VFT (AVX) Figure S2: Speedups obtained with different number of threads using as reference the sequential execution (one thread) of FastTree-2. 4 10 1000 20 100 Number of floats (array size) 0.95 1 1.05 1.1 1.15 1.2 Speedup vector_multiply vector_add_mult AVX vs. SSE (single precision) Figure S3: Speedup obtained by AVX in comparison to SSE using single precision for arrays of different size. ”-fastest” option for the large dataset, which is recommended for over 50,000 sequences. Figure S1 shows the running times for the three datasets using different number of threads with single and double precision arithmetic. We must take into account that FastTree-2 only uses SSE instructions with single precision calculations, while VeryFastTree supports SSE and AVX/AVX2 in all cases. Results for AVX512 are not displayed because they are similar to the ones obtained by AVX instructions. An analysis of this behavior can be found in Section 1.2. The scalability of FastTree-2 is limited to a maximum of 4 threads so that using more threads does not improve performance. For this reason graphs only display the FastTree-2 running times for 1, 2 and 4 threads. Figure S2 shows the previous results in terms of the speedup, calculated using as reference the sequential execution of FastTree-2. Performance results in Figures S1 and S2 show that under single precision there are small differences S2 4 10 20 100 1000 5000 Number of doubles (array size) 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 Speedup vector_multiply vector_add_mult AVX512 vs. AVX (double precision) Figure S4: Speedup obtained by AVX512 in comparison to AVX using double precision for arrays of different size. Experiments conducted on an Intel Xeon Gold 6248 processor (20 cores). between AVX and SSE, which means that computations are memory bandwidth bound. To demonstrate this, we implemented a benchmark that inserts scalar operations (non-useful work) between consecutive calls to a considered vector kernel to simulate the behavior of VeryFastTree when using SSE and AVX instructions. We measured the execution time after 10,000 vector operations, excluding the scalar time. Figure S3 shows the speedup obtained by AVX in comparison to SSE when considering arrays of different size and two VeryFastTree vector kernels (vector multiply and vector add mult). Arrays for vector operations when using protein alignments in VeryFastTree and FastTree-2 contain twenty floats (green line in the graph). Results demonstrate that AVX outperforms SSE except for small arrays, which is precisely the VeryFastTree case. Focusing on the green line, AVX is still the best option for the vector add mult kernel. The reason is that vector add mult contains more vector operations than vector multiply (codes of both kernels are available in the VeryFastTree repository). Therefore, improvements of the AVX version are related to the small performance differences detected for the most computing-intensive vector kernels. 1.2 AVX512 VeryFastTree supports AVX512 vector instructions. However, we do not recommend using them since the performance is always lower than the one obtained by AVX. The reason is that, even for double precision arithmetic, the likelihood calculations in VeryFastTree are entirely memory bandwidth bound. To demonstrate this we use the same benchmark detailed above but considering AVX and AVX512 instructions with double precision operands instead. Results varying the array size in terms of the AVX512 speedup are displayed in Figure S4. It can be observed that using AVX512 outpeforms AVX in most cases, improving the speedup as the array size increases. However, vector kernels in VeryFastTree and FastTree-2 operate on small arrays (green line in the figure), which results in a degradation in the AVX512 performance in comparison to AVX. In addition, there is another important factor to take into account. Modern Intel CPUs reduce their frequency when executing wide vector operations (AVX2 and AVX512 instructions), as these instructions increase power consumption [1]. The frequency is only increased again two milliseconds after the last code section containing such instructions has been executed in order to prevent excessive numbers of frequency changes. Due to this delay, intermittent use of wide vector operations can slow down the rest of the system significantly. That is exactly the pattern followed by VeryFastTree and FastTree-2, which constantly alternates vector and scalar operations. S3 Threads COG438 COG583 COG596 COG642 COG1028 COG1309 COG2814 Avg. 1 81.99% 86.25% 86.73% 80.06% 85.03% 79.45% 88.16% 83.95% 2 82.25% 86.25% 87.05% 80.00% 85.15% 79.35% 88.22% 84.04% 4 82.63% 86.25% 86.81% 80.08% 85.51% 79.75% 88.44% 84.21% 8 82.99% 86.45% 86.89% 80.16% 85.23% 79.17% 88.30% 84.17% 10 82.65% 86.09% 86.65% 80.10% 85.19% 79.19% 88.14% 84.00% 12 82.37% 88.17% 86.45% 79.96% 85.29% 79.45% 88.20% 84.27% Table S1: Topological accuracies obtained by VeryFastTree. Threads COG438 COG583 COG596 COG642 COG1028 COG1309 COG2814 Avg. 1 82.13% 86.43% 86.37% 80.05% 85.01% 79.45% 88.16% 83.94% 2 (min.) 82.29% 86.49% 86.75% 79.76% 84.97% 78.93% 88.48% 83.95% 2 (max.) 82.83% 86.87% 87.13% 80.52% 85.43% 79.55% 88.82% 84.45% 4 (min.) 82.19% 86.51% 86.67% 80.04% 85.05% 79.49% 88.40% 84.05% 4 (max.) 82.33% 86.73% 87.19% 80.80% 85.47% 79.85% 88.62% 84.43% Table S2: Topological accuracies obtained by FastTree-2. 1.3 Memory requirements Given that VeryFastTree was implemented in C++ instead of C, there is a small memory penalization caused by the allocation of objects. For instance, considering double precision, the memory requirements of both tools are: •Small dataset (sequential): FastTree-2 (1.75 GB) — VeryFastTree (1.90 GB) •Small dataset (8 threads): FastTree-2 (1.80 GB) — VeryFastTree (1.99 GB) •Medium dataset (sequential): FastTree-2 (5.64 GB) — VeryFastTree (6.09 GB) •Medium dataset (8 threads): FastTree-2 (5.70 GB) — VeryFastTree (6.22 GB) •Large dataset (sequential): FastTree-2 (60.13 GB) — VeryFastTree (61.99 GB) •Large dataset (8 threads): FastTree-2 (60.49 GB) — VeryFastTree (62.62 GB) We must also take into account that the VeryFastTree values were obtained using AVX instructions. Therefore, there is some extra memory needed for vector operations. In any case, differences are small. 1.4 Topological accuracy We have measured the topological accuracy of the trees generated using 5,000 sequences protein alignments2. Note that this dataset was used for the same purpose in the FastTree-2 original paper [2]. Accuracies were calculated using the Robinson-Foulds distance (ETE toolkit3). The script used in the experimentation, treecmp.py, is available in our repository. The VeryFastTree topological accuracies were calculated, for each number of threads, as the average of seven values (one for each alignment in the dataset). Table S1 shows these results. On the other hand, parallel executions of FastTree-2 are not deterministic. As a consequence, we provide a range of accuracies for each number of threads. This range 3http://etetoolkit.org S4 Total (hours) Tree estimation (hours) PastaSpark (FastTree-2) 26.1 18.0 PastaSpark (VeryFastTree) 23.4 15.3 Table S3: Running times of PASTASpark. was calculated as the average of the minimum and maximum values obtained after 5 executions for each alignment and thread count (see Table S2). According to the results, we conclude that VeryFastTree is deterministic and, at the same time, it produces trees with the same accuracy as FastTree-2. 1.5 Integrating VeryFastTree with other tools To facilitate the adoption from the Bioinformatics community, VeryFastTree keeps exactly the same command line arguments than FastTree-2. This way, it is only necessary to replace the call to FastTree-2 by a call to VeryFastTree using the same options to increase the overall performance. We can find many tools that use FastTree-2 at their core. A good example is PASTASpark [3], a Big Data-based multiple sequence alignment (MSA) tool which produces highly accurate alignments, improving the accuracy and scalability of other state-of-the-art methods. The PASTASpark iterative workflow consists of four steps. In the last one, FastTree-2 is executed to estimate a ML tree on the MSA produced on the previous steps. We will consider this application to illustrate the benefits of using VeryFastTree instead of FastTree-2. The following experiments were obtained using as input the largest dataset (200k RNASim, 200k sequences, 3.4 GB) evaluated in the PASTASpark paper, running on a 16-node Linux cluster with CentOS 7 and Spark 2.2.0. Each node consists of two 10-core Intel Xeon E5-2630v4 processors, 384 GB of memory and a 10 GbE network. The setup for Spark jobs consists of 64 executors (4 threads per executor) and 8 threads for the driver. Therefore, VeryFastTree is parallelized using 8 threads. PASTASpark limits to 4 threads the parallel executions when considering FastTree-2 due to its scalability issues. Performance results using double precision are displayed in Table S3. VeryFastTree decreases the total time in 2.7 hours, i.e., a reduction of 10.3%. 2 Changes and Optimizations VeryFastTree is a highly-tuned implementation of the FastTree-2 tool that takes advantage of parallelization and vectorization strategies to speed up the inference of phylogenies for huge alignments. It is important to highlight that VeryFastTree keeps unchanged the phases, methods and heuristics used by FastTree-2 to estimate the phylogenetic tree (the interested reader can find a detailed description in [2]). Next we list the most important changes and optimizations carried out in VeryFastTree with respect to FastTree-2: •Compilation: FastTree-2 was implemented in a single C file. It generates different binaries depending on the selected features when compiling (SSE2, OpenMP and/or single/double precision). SSE2 and OpenMP are not supported for Windows systems. VeryFastTree has been implemented in C++ and keeps total compatibility with Windows, Linux and MacOs operating systems. In addition, VeryFastTree generates a single binary where the different features can be enabled/disabled using program arguments. This way, it is only necessary to compile the source code one time. VeryFastTree uses the cmake tool to automate the compilation process and remove compilerdependent parameters. In addition, it allows to detect automatically the characteristics of the hardS5 ware platform on Linux and Mac systems, enabling in the executable only the supported features. For Windows machines, it is necessary to pass as argument a list of features before the compilation. •Although FastTree-2 has a parallel OpenMP implementation, the most important and time-consuming steps in FastTree-2 such as the ML-based NNI rearrangements are not parallelized at all or only in a small part. As a consequence, the scalability of FastTree-2 is very limited in such a way that using more than 3 or 4 threads do not provide any benefit. For this reason, VeryFastTree introduces a new tree partitioning method that allows to perform in parallel the rounds of ML NNIs. This method is guided by an objective function whose main goal is finding a good load balance among threads. A detailed explanation of the partitioning algorithm can be found in Section 3. This method also allows to parallelize the recomputation of the posterior distributions for each internal node in the tree before starting the rounds of ML NNIs. On the other hand, note that VeryFastTree keeps unchanged the remainder steps which were parallelized in FastTree-2 (see [4] for more information about the parallelization in FastTree-2). •The OpenMP version of FastTree-2 is not deterministic, which is unlikely to affect the quality of the final tree but there may be slight differences with respect to the sequential execution and even among parallel executions using the same number of threads. The cause is that there are several locations in the FastTree-2 parallel code where race conditions could arise. VeryFastTree solved this issue in such a way that parallel executions using the same number of threads produce always a tree with the same topology. •Over the years, x86 processors have added more and more vector capabilities, layered one on the other, starting with MMX through several versions of SSE to AVX512. Basically, vector instructions operate on multiple values contained in one large register at the same time. FastTree-2 supports SSE2 operations to accelerate the single precision (32 bits) computations, which allows to compute on 4 floats at the same time using 128-bit registers. However, SSE2 intrinsics are not considered with double precision arithmetic. On the other hand, SSE2 was introduced in 2001, while current processors support more modern SSE versions (SSE3 and SSE4) and the more powerful AVX/AVX2 and AVX512 vector operations. For this reason, VeryFastTree boosts the performance of the arithmetic calculations supporting SSE3, AVX/AVX2 and AVX512 vector operations for both single and double precision executions. •The exponential calculation (y=ex) is one of the most common and costly arithmetic operations when building the tree. To go a step further in terms of performance, VeryFastTree implements a very efficient and fast algorithm to compute an accurate approximation of ex. It is based on the well-known Cephes 4math library, but our implementation supports SSE, AVX and AVX2 vector operations in order to accelerate even more the computations. By default this option is not enabled. To use it, users should pass ”-fastexp” as argument in the command line (see Section 4 for a complete list of options in VeryFastTree). •A very efficient parallel implementation of the samplesort algorithm is used instead of the sequential quicksort approach considered by FastTree-2. While quicksort takes on average O(nlog n)comparisons to sort nitems and O(nlog n)and O(n2)comparisons in the best and worst cases, the sorting algorithm in VeryFastTree has the following complexities O(nlog n),O(n)and O(nlog n) for average, best and worst cases, respectively. 4https://github.com/nearform/node-cephes/tree/master/cephes S6 •To speed up hashing, VeryFastTree uses the robin-map 5library, a C++ implementation of a fast hash map and hash set using open-addressing and linear robin hood hashing with backward shift deletion to resolve collisions. •To facilitate the adoption from the Bioinformatics community, VeryFastTree keeps exactly the same command line arguments than FastTree-2. This way, it is only necessary to replace the call to FastTree-2 by a call to VeryFastTree using the same options. In any case, VeryFastTree changes how these arguments are parsed, replacing the non-flexible sequence of ”if-else” commands of FastTree-2 by a powerful parser. In particular, VeryFastTree uses CLI 6. •Other changes/optimizations: –All structures in the FastTree-2 code have been replaced by classes in VeryFastTree. Methods that initialize the structures are now the constructors of the classes. The destruction methods have been eliminated leaving the management to the class itself, and the objects are cleaned for reuse instead of creating new ones whenever possible. –Most of the pointers in the code have been replaced by references. 3 Tree partitioning in the Maximum-Likelihood phase The most time consuming phase to estimate the phylogenetic tree corresponds to the maximum-likelihood phase. In particular, the rounds of ML NNIs (Nearest-Neighbour-Interchanges) take most of the time of this phase. The FastTree-2 implementation optimizes likelihoods for 3 alternate topologies in parallel during ML NNIs in such a way that only 3 threads can be used, which limits significantly its performance and scalability. VeryFastTree uses a different approach to speed up ML NNIs. Specifically, it splits the tree into several subtrees in such a way that independent subtree computations are performed in parallel. Our methodology has two steps: 1. The partitioning method itself guided by an objective function. 2. The assignment of subtrees to threads. Before introducing the considered objective function, it is necessary to understand how the tree is traversed and the data dependencies during each round of NNIs. We will use the tree of Figure S5 as example. Each node in the tree is visited before its parents, which means that a depth-first post-order traversal is used [2]. Let’s assume that the tree is split into two subtrees assigned to threads t1and t2. Since the root node of the tree is not assigned to a thread, it will be labeled as not-assigned. At each node N, the likelihood of the trees AB|CD,AC|BD, and AD|BC is compared, where Aand Bare its children, Cis its sibling, and Dis the rest of the tree. FastTree-2 (and, therefore, VeryFastTree) uses the posterior distributions for A,B, and Cand an up-posterior D, which represents the rest of the tree. The posterior distribution for an internal node describes the state of the corresponding ancestor, given the branch lengths and the sequences beneath it. The up-posterior Dis the posterior distribution of the node’s parent N, given all of the nodes that are not children of N. Due to these data dependencies, the root node of each subtree and its children can not be computed in parallel since some required data should be obtained from other subtrees. As a consequence, the root node of each subtree and its children are also labeled as not-assigned. Each round of ML NNIs finishes processing all the not-assigned nodes sequentially. 5https://github.com/Tessil/robin-map 6https://github.com/CLIUtils/CLI11 S7 Figure S5: Example of the tree traversal and data dependencies when computing ML NNIs. The graph represents only internal nodes (no leaves). Blue and green nodes are processed in parallel by threads t1 and t2, respectively. Red nodes (labeled as not-assigned) are processed sequentially after processing all subtrees. Therefore, in the example of Figure S5, blue and green nodes are processed in parallel using two threads, while red nodes are processed afterwards using one thread. The total execution time will be the sum of the time spent by the slowest thread and the time required to process the not-assigned nodes. Obviously, the slowest thread will be the one that owns more internal nodes (thread t1in the example). On the other hand, we want to highlight that our partitioning method also allows to parallelize the recomputation of the posterior distributions for each internal node in the tree before starting the rounds of ML NNIs. In this case, data dependencies are only between each node and its children. Therefore, unlike ML NNIs, all the nodes in a subtree can be computed in parallel. Partitioning method Our partitioning method is guided by an objective function, whose goal is two-fold. First, it should balance the workload assigned to each thread. That is, the number of internal nodes to be processed by each thread. To model the workload we define the weight of a node as the number of internal nodes beneath it plus itself. This way, the workload associated to process a subtree is the weight of its root node. The second goal of the objective function should be to reduce to the minimum the amount of sequential work, which corresponds to processing the not-assigned nodes. According to the Amdahl’s law [5], a small percentage of sequential work harms noticeably the performance and scalability of a parallel application. Let Tbe the input tree to be partitioned, which only contains internal nodes. A partition P(T)is a set of disjoint binary subtrees siincluded in T. We must highlight that not all nodes in the tree should be included in a partition. Let |P(T)|denote the number of subtrees in the partition and nthe number of threads. Each subtree si∈P(T)is assigned to a particular thread in such a way that Pj(T)contains all subtrees in the partition P(T)assigned to thread j. This way, we can calculate the weight of the subtrees assigned to a particular thread jas weight(Pj(T)) = {Pweight(si)|si∈Pj(T)}, which represents its workload. We evaluate a partition P(T)using the following objective function: S8 max{weight(Pj(T)) |1≤j≤n}+nodes not assigned −3× |P(T)|(1) Note that nodes not assigned include those nodes labeled as not-assigned following the procedure detailed above, and those nodes not included in the partition P(T). Therefore, nodes not assigned models the amount of sequential work. Our objective function favors those partitions with a high number of subtrees because it facilitates the load balancing among threads. The goal of our partitioning method will be to find a partition P(T)that minimizes the value of the objective function, whose pseudocode is detailed in Algorithm 1. Our procedure to find the best partitioning is based on two operations on subtrees: splitting and removal. The first operation removes a subtree from the current partition P(T), and creates a new one P0(T)inserting the children of the removed subtree instead. In other words, it splits the considered subtree. The second operation creates a new partition P0(T)removing one subtree from P(T). If we apply all the possible combinations of both operations to the same partition with the aim of finding the solution which reduces the most the objective function, the search space would be huge. For this reason, our method uses a greedy approach. In particular, a new partition is obtained in each iteration of the algorithm after splitting only the largest subtree (that is, the most weighted), and removing only the smallest subtrees. At the end of each iteration, the new resulting partition is evaluated using Equation 1. If the quality of the new partition is better than the current best one, then the new partition can be considered the best until that moment. The method iterates until the splitting and removal operations cannot be applied to the considered partition or the most weighted subtree in the partition is below a size threshold. S9