On the optimality of the neighbor-joining algorithm

The popular neighbor-joining (NJ) algorithm used in phylogenetics is a greedy algorithm for finding the balanced minimum evolution (BME) tree associated to a dissimilarity map. From this point of view, NJ is "optimal" when the algorithm outputs the tree which minimizes the balanced minimum evolution criterion. We use the fact that the NJ tree topology and the BME tree topology are determined by polyhedral subdivisions of the spaces of dissimilarity maps R+(n2) to study the optimality of the neighbor-joining algorithm. In particular, we investigate and compare the polyhedral subdivisions for n ≤ 8. This requires the measurement of volumes of spherical polytopes in high dimension, which we obtain using a combination of Monte Carlo methods and polyhedral algorithms. Our results include a demonstration that highly unrelated trees can be co-optimal in BME reconstruction, and that NJ regions are not convex. We obtain the l2 radius for neighbor-joining for n = 5 and we conjecture that the ability of the neighbor-joining algorithm to recover the BME tree depends on the diameter of the BME tree.


Introduction
The popular neighbor-joining algorithm used for phylogenetic tree reconstruction [1] has recently been "revealed" to be a greedy algorithm for finding the balanced minimum evolution tree associated to a dissimilarity map [2]. This means the following: Let be a dissimilarity map (this is an n × n symmetric matrix with zeroes on the diagonals and nonnegative real entries). The balanced minimum evolution problem is to find the unrooted binary tree T with n leaves that minimizes Here where if i ≠ j and . In [3], Day shows that choosing a minimizing tree for (2) from among the (2n-5)!! unrooted binary trees is an NPhard problem. Yet it is desirable to find algorithms for minimizing (2) because of the following statistical interpretation:

Theorem 1.2 ([4])Let T be a binary tree with edge lengths given by l: E(T)
→ + and a dissimilarity map. If the variance of d ij is proportional to (i. e., var(d ij ) = for some constant c) then (2) is the minimum variance tree length estimator of T. Moreover, the weighted least squares tree length estimate is equal to (2).
This result provides a weighted least squares rationale for the minimization of (2), and highlights the importance of understanding the balanced minimum evolution polytope:

Definition 1.3 The balanced minimum evolution polytope is the convex hull of the vectors
Example. There are four trees with n = 4 leaves. They are the 3 binary trees and the star-shaped tree. In this case the balanced minimum evolution polytope is the convex hull of the vectors: The balanced minimum evolution polytope in this case is a triangle in 6 . Note that the star-shaped tree is in the interior of the triangle.
For any dissimilarity map, the trees which minimize (2) will be vertices of the balanced minimum evolution polytope; these are always the binary trees. In fact, for such trees ; this is Pauplin's formula [5].
The BME polytope lies in and has dimension -n. The normal fan [6] of the BME polytope gives rise to BME cones which form a polyhedral subdivision of the space of dissimilarity maps . They describe, for each tree T, those dissimilarity maps for which T minimizes (2). We provide an introduction to the necessary polyhedral combinatorics in Section 2, and discuss the polytope in more detail in Section 3.
The neighbor-joining algorithm is a greedy algorithm for finding an approximate solution to (2). We omit a detailed description of the algorithm here -readers can consult [2] -but we do mention the crucial fact that the selection criterion is linear in the dissimilarity map [7]. Thus, the NJ algorithm will pick pairs of leaves to merge in a particular order and output a particular tree T if and only if the pairwise distances satisfy a system of linear inequalities, whose solution set forms a polyhedral cone in . We call such a cone a neighbor-joining cone. or NJ cone. The NJ algorithm will output a particular tree T if and only if the distance data lies in a union of NJ cones. In Section 4 we show that the NJ cones partition , but do not form a fan. This has important implications for the behavior of the NJ algorithm.
Our main result is a comparison of the neighbor-joining cones with the normal fan of the balanced minimum evolution polytope. This means that we characterize those dissimilarity maps for which neighbor-joining, despite being a greedy algorithm, is able to identify the balanced minimum evolution tree. These results are discussed in Section 5.

Polyhedral preliminaries
In this section we will introduce some of the elementary polyhedral combinatorics necessary for this paper. For more details see [8].
where V ⊆ d is a subspace and a ∈ A. V is uniquely determined by A and the affine dimension of A is defined to be the dimension of V.
Given two distinct points x, y ∈ d , the set [x, y] = {αx + (1 -α)y : 0 ≤ α ≤ 1} of all convex combinations of x and y is called the interval with endpoints x and y.
Then the set is called a polyhedron. The convex hull of a finite set of points in d is called a polytope and the Weyl-Minkowski Theorem says that a polytope is a bounded polyhedron [9]. Polytopes are familiar objects in geometry. In the plane, polytopes are precisely the convex polygons. In 3 , examples of polytopes are shown in Figure 1. The dimension dim P of a polytope or polyhedron P is defined to be the dimension of the affine hull of P.
A (d -1) dimensional affine set in d is called a hyperplane and every hyperplane can be represented as {x ∈ d : n·x = b} for some n ≠ 0 ∈ d and b ∈ , where n·x is the dotproduct of n and x. We call n a normal vector of this hyperplane.
A polyhedron C is a cone if it can be written as The four types of facets of P Figure 1 The four types of facets of P.
for some y 1 , ..., y N ∈ d . This is equivalent to the existence Given a face F of a polytope P, the normal cone N(F) is the set of all vectors c for which c·v = max x∈P c·x for all v ∈ F.
The collection of relative interiors of normal cones of faces of P partition d , and for each face we have dim(F) + dim(N(F)) = d. The collection of normal cones of faces of P is called the normal fan of P.
Given a polyhedron P, the lineality space of P is the set of vectors v for which y + c·v ∈ P for all y ∈ P and c ∈ R. The largest such subspace is called lineality space of P. If a polyhedron P has lineality space V, we can let V' be the orthogonal complement V' (i.e. V ⊕ V' = d ) and consider the polyhedron P' := P ʝ V', which has lineality space {0}.

The balanced minimum evolution polytope
Throughout this paper we work with binary unrooted trees on n leaves labeled {1, ..., n}. Such trees are also known as phylogenetic X-trees. We refer the reader to [10] for more detail about such trees, and for related definitions. Recall there are 2n -3 edges in an unrooted tree with n leaves. For a fixed tree topology T, let B T be the × (2n -3) matrix with rows indexed by pairs of leaves and columns indexed by edges in T defined as follows: For example, for the tree in Figure 2, where its rows are indexed by pairs of leaves (1, 2), (1, 3), (2,3), (1,4), (2,4), (3,4), (1,5), (2,5), (3,5), (4,5) and its columns are indexed by edges (1, a), (2, a), (3, b), (4, c), with a is an internal node adjacent to leaves 1 and 2, c is an internal node adjacent to leaves 4, 5, and b is an internal node adjacent to nodes 3, a and c. Given edge lengths l : E(T) → + we let b be the vector with components l(e) as e ranges over E(T). Any dissimilarity map d (encoded as a row vector) can now be written as where e is a vector of "error" terms that are zero when d is a tree metric.
The weighted least squares solution for the edge lengths b assuming a variance matrix V with off-diagonal entries (as defined in the introduction) and dissimilarity map d is given by where · t denotes matrix transpose. The length of T with respect to the least squares edge lengths is then where 1 and 1 is the vector of all 1's. We call the vectors v T the balanced minimum evolu-  tion vectors (or BME vectors). In the case of Figure 2, the BME vector is The BME method is equivalent to minimizing the linear functional v T ·d over all BME vectors for all tree topologies T. The BME polytope is the convex hull of all BME vectors in . The following facts follow from the definition of the balanced minimum evolution tree: The vertices of the BME polytope are the BME vectors of binary trees. The BME vector of the star phylogeny lies in the interior of the BME polytope, and all other BME vectors lie on the boundary of the BME polytope.
The normal fan of a BME polytope partitions the space of dissimilarity maps into cones, one for each tree. We call these BME cones. They completely characterize the BME method: T is the BME tree topology if and only if the dissimilarity map D lies in the BME cone of T.
For a leaf node a in a binary unrooted tree, the shift vector s a is the dissimilarity map in which a is at distance 1 from all other leaves, and all other distances are 0 (see [11] for the description of shift vectors). According to [5], for a tree T, (v T ) ab gives the probability that a will immediately precede b in a random circular ordering of T. Thus the dotproduct of a BME vector with a shift vector must necessarily equal 1, and in fact the lineality space of BME cones is spanned by shift vectors. So when we describe a BME cone we will always describe just the pointed component, i.e. modulo the lineality space of shift vectors.
As part of our computational study, we computed the BME polytope and BME cones for trees with n = 4, 5, 6, 7, 8 leaves using the software polymake [12]. In Table 1 we display some of the components of f-vectors we were able to compute. This provides information about the polytopes: Recall that the ith component of the f-vector of a polytope is the number of faces of dimension i -1. For example, the first component in each vector in Table 1 is the number of 0-dimensional faces (vertices) of the corresponding BME polytope, i.e., the number of binary trees.
We found that the edge graph of the BME polytope is the complete graph for n = 4, 5, 6 which means that for every pair of trees T 1 and T 2 with the same number (≤ 6) of leaves, there is a dissimilarity map for which T 1 and T 2 are (the only) co-optimal BME trees. However, for n = 7, the BME polytope does in fact have one combinatorial type of non-edge. Namely, two bifurcating trees with seven leaves and three cherries (two leaves adjacent to the same node in the tree) will form a non-edge if and only if they are related by two leaf exchanges as depicted in Figure 3. This completely characterizes the non-edges for n = 7. It is an interesting open problem to characterize the non-edges of the BME polytope in general.

Neighbor-joining cones
The neighbor-joining algorithm takes as input a dissimilarity map and outputs a tree. The tree is constructed "one cherry at a time". In each step the algorithm chooses a pair of leaves a and b that minimize the Q-criterion, which is defined by the formula The non-edges on the BME polytope for n = 7 Figure 3 The non-edges on the BME polytope for n = 7. Two trees will form a non-edge if and only if they are trees that have three cherries, and differ by the pair of leaf exchanges shown in the figure. There are two ways to perform each leaf-exchange, so each binary tree with three cherries is not adjacent to 4 trees. cones is a subset -but not necessarily a face -of the boundary of each of the cones. Given an input from the interior of C σ , the NJ algorithm will pick the cherries in the order σ and output the corresponding tree. For inputs d on the boundary of one (and therefore at least two) of the cones, the order in which NJ picks cherries is undefined, because at some point there will be two cherries both of which have minimal Q-criterion. We call the cones C σ neighbor-joining cones, or NJ cones. See [11] for the hyperplane representation of NJ cones and descriptions how to construct each cone.
Example. There is only one unlabeled binary tree with 5 leaves and there are 15 distinct labeled trees. For each labeled tree, there are two ways in which a cherry might be picked by the NJ algorithm in the first step. For instance, neighbor-joining applied to any dissimilarity map in C 12,45 or C 45,12 will produce the tree in Figure 2. There are a total of 30 NJ cones for n = 5.
We note that all Q-criteria for shift vectors equal -2, so adding any linear combination of shift vectors to a dissimilarity map does not change the relative values of the Qcriteria. Also, after picking a cherry, the reduced distance matrix of a shift vector is again a shift vector. Thus, for any input vector d, the behavior of the NJ algorithm on d will be the same as on d + s if s is any linear combination of shift vectors. In fact it can be shown that the lineality space of NJ cones is spanned by shift vectors, just as for BME cones [11]. So from now on, when we refer to NJ cones, we will mean the pointed portion of the cone, i.e. modulo the lineality space.

Theorem 4.1
The cones in do not form a fan. In particular, they are not the normal fan of any polytope for n ≥ 5.
The theorem follows from that fact that the NJ cones have rays which are on the boundary of other cones but not rays of them. Thus there are pairs of cones whose intersection is not a face of both cones. We describe the case n = 5 in detail; it also suffices to prove the theorem.
We begin by noting that all of the NJ cones are equivalent under the action of the symmetric group on five elements (S 5 ), where an element of S 5 permutes the five taxa or, equivalently, the rows and columns of the input distance matrix. Each NJ cone is defined by inequalities that are implied by the Q-criteria as the NJ algorithm picks the two cherries. The cones are 5-dimensional, and their intersection with a suitable hyperplane leaves a four dimensional polytope P. The f-vector of P is (14,32,27,9).
The 30 cones share many of their rays, giving a total of 82 rays which decompose into three orbits under the action of S 5 . We refer to the types of rays as Type I, Type II and Type III. Each cone has 6 rays of type I, 4 rays of type II and 4 rays of type III. Each ray of type I is the common ray of 3 cones, and belongs to 2 other cones of which it is not a ray (i.e. it is in the interior of a face). Note that this implies that the cones cannot form a fan. The type II rays are contained in 10 cones each, and the type III rays in 12. Type II and III rays are rays of all cones which contain them. For the cone C 23,45 , this information is tabulated in Table 2.
We note that the rays of NJ cones are minimal intersections of NJ cones, and thus give dissimilarity maps for which the NJ algorithm is least stable.
Example. Consider two alignments of 5 sequences that are to be used to construct a tree. These may consist of two different genes and for each of them the homologs among 5 genomes. Suppose that distances are estimated using the Jukes-Cantor correction [6,13] where g ij is the fraction of different nucleotides between sequences i and j in the second set.
If the fractions f ij and g ij are given by then we obtain Notice that the vector representation of D 1 lies in the cone C 12,45 and the vector representation of D 2 lies in the cone C 45,12 . Thus NJ returns the same tree topology for both D 1 and D 2 .
If we concatenate the alignments and combine the data to build one tree, then we estimate the distances using the average of f and g:

. Each ray is determined by a vector shown in the second column. The third column shows, for each ray, which cones it belongs to. If a cone is starred then the ray is on the boundary of that cone, but not a ray of it.
Type rays Cones  Using this frequency matrix we obtain the distance matrix D 3 via the Jukes-Cantor correction: However, the vector representation of D 3 lies in the cone C 24, 15 , which means that neighbor-joining returns a different tree topology for D 3 . This example provides a distancebased recon-struction analog to the recent mixture model results of [14].
An analysis of the rays of suffices to prove Theorem 4.1. but the facet structure of each cone is also informative, and we were able to obtain complete information for n = 5. The types of facets constituting each cone are shown in Figure 1. We used our description of the NJ cones to examine the l 2 distance between tree metrics and the boundaries of NJ cones. Without loss of generality, by shifting the leaves in the cherries, we can assume the tree metric is of the form where α and β are the internal branch lengths, α ≥ and α + β = 1. It is easy to see that D T ∈ C 12,45 confirming the consistency of neighbor-joining. The cone C 12,45 contains 9 faces, but we may ignore one of them (namely the one shared with C 45,12 ) because it is shared with a cone result-ing in the same tree topology. The distance to the closest of the remaining eight faces is The l 2 radius is obtained by dividing (4) by min(α, β), so the minimum is attained at α = β =

Theorem 4.2
The l 2 radius of neighbor-joining for 5 taxa is ≈ 0.5773. This is slightly larger than the l ∞ radius of given by Atteson's theorem [15]. It is an interesting problem to compute the l 2 radius for neighbor-joining with more taxa.
The description of the NJ cones we have provided can also be used in practice to evaluate the robustness of the algorithm when used with a specific dataset. For n = 5, we examined data simulated from subtrees of the two tree models T 1 and T 2 in [16] with the Jukes-Cantor model and the Kimura 2-parameter models [6]. For each of 40, 000 simulations, we calculated the ᐍ 2 -distance between the NJ cone of the given tree and the maximum likelihood estimates for the pairwise distances (see supplementary material). These show that in many cases the maximum likelihood estimates lie very close to the boundary. In such cases, one must conclude that the NJ tree is possibly incorrect due to the variance in the distance estimates.

Optimality of the neighbor-joining algorithm
In order to study the optimality of the neighbor-joining algorithm, we compared the BME cones with the NJ cones. Such a comparison involves intersecting the cones with the ( -1)-sphere (in the first orthant) and then studying the volumes of their intersection by computing the standard Euclidean volume of the resulting surfaces. These surfaces are an intersection of closed hemispheres, i.e. spherical polytopes. Computing Euclidean volumes of (non-spherical) polytopes is a standard problem that is usually solved by triangulating and summing the volumes of the simplices. However there has been no publicly available software developed for computing or approximating volumes of spherical polytopes of dimension > 3 using this method. One possible reason for this is that in higher dimensions the volumes of spherical simplices are given by complicated analytical formulas [17] whose computational complexities are unknown.
We implemented two approaches in MATLAB (using polymake as a preprocessing step) for approximating the volume of a spherical polytope P. One approach is trivial: it simply samples uniformly from the sphere, and counts how many points are inside P. This approach is particularly suitable if P has large volume, or if many spherical polytopes are being simultaneously measured which partition the sphere, as is the case for NJ and BME cones. The second approach is suitable for spherical polytopes having small volume. We used this approach for computing the volumes of consistency cones [18] which we discuss briefly in the Discussion section. Our main results on the optimality of NJ for n = 5, 6, 7, 8 taxa are summarized in Table 3. Each row of the table describes one type of tree. Trees are classified by their topology. A k-cherry tree is a tree with k cherries. The NJ volume column shows the volume of that part of the positive orthant of dissimilarity maps for which the NJ tree is of the specified type. Similarly, the BME volume column shows the same statistic for BME trees. Finally, NJ accuracy shows the fraction of the BME cone that overlaps the NJ cone. In other words, NJ accuracy is a measure of how frequently NJ will find the BME tree for a dissimilarity map that is chosen at random.
We also classified and measured the intersections of NJ and BME cones in which the NJ tree differs from the BME tree. Many of these intersection cones are equivalent under the action of S n on the leaf labels, particularly as the stabilizer of the BME tree permutes the leaf labels in the NJ tree. In fact, for n = 5 taxa there are only three types of mistakes that the NJ algorithm can make when it fails to reproduce the BME tree. These are depicted in Figure 4 and Frequencies of the all three possible types of NJ trees that may picked instead of the BME tree for n = 5 leaves Figure 4 Frequencies of the all three possible types of NJ trees that may picked instead of the BME tree for n = 5 leaves. Neighbor-joining agrees with the BME tree 98.06% of the time. the normalized spherical volumes of corresponding NJ/ BME intersection cones are given. Figure 4 can be interpreted as follows: For a random dissimilarity map, if the NJ algorithm does not produce the BME tree, then with probability 0.67 it produces the tree on the right, and if not then it almost always produces the tree in the middle. This tree differs from the BME tree significantly. A surprising result is that the tree on the left is almost never the NJ tree. We believe that a deeper understanding of the "mistakes" NJ makes when it does not optimize the balanced minimum evolution criterion may be important in interpreting the results, especially for large trees.
We also computed analogous results for n = 6, 7, 8, 9, 10. They are available, together with the software for computing volumes at [19].

Discussion
Theoretical studies of the neighbor-joining algorithm have focused on statistical consistency and the robustness of the algorithm to small perturbations of tree metrics. The paper by [20] established the consistency of NJ, that is, if D T is a tree metric then NJ outputs the tree T. This result was then extended in [15] and more recently by [18] who show that if D is "close" to a tree metric D T for some T, then NJ outputs T on input D.
Our results provide a different perspective on the NJ algorithm. Namely, we address the question of the accuracy of the greedy approach for the underlying linear programming problem of BME optimization. This led us to the study of BME polytopes, and the combinatorics of these polytopes is interesting in its own right:  [21]for an example of a local search approach to finding minimum evolution trees.
Similarly, a better understanding of the combinatorics of the NJ cones will lead to a clearer view of the strengths and weaknesses of the neighbor-joining algorithm. A basic problem is the following:

Question 6.2
Find a combinatorial description of the NJ cones for general n. How many facets/rays are there?
Our computational results lend new insights into the performance of the NJ and BME algorithms for small trees. We have measured the relative sizes of cones for different shapes of trees, and measured the frequencies of all combinatorial types of discrepancies between BME and NJ trees. In particular, we have observed that the NJ algorithm is least likely to reproduce the BME tree when the BME tree is the caterpillar tree.

Conjecture 6.3
For n > 6, it is the caterpillar tree that yields the smallest ratio of spherical cone volumes vol(NJ ʝ BME)/vol(BME) where NJ is the spherical cone volume of a union of the NJ cones and BME is the spherical cone volume of the BME cone for a fixed tree. In other words, the caterpillar tree is the most difficult BME tree topology for the NJ algorithm to reproduce.
Another problem we believe is very important is to extend the results shown in Figure 4 to large trees. In other words, to understand how neighbor-joining can fail when it does not succeed in finding the balanced minimum evolution tree.

Question 6.4
What tree topologies is neighbor-joining likely to pick when it fails to construct the balanced minimum evolution tree?
There are many other interesting cones related to distancebased methods that can be considered in this context. For example, in [18], it is shown that the quartet consistency condition is sufficient for neighbor-joining to reconstruct a tree from a dissimilarity map for n ≤ 7 leaves. The quartet consistency conditions define polyhedral cones (consistency cones) in ; see [18] for details. For n = 4 taxa the consistency cones cover all of showing that quartet consistency explains the behavior of neighbor-joining for all dissimilarity maps. Using the second method outlined in Section 4 we succeeded in computing the volumes of the consistency cones intersected with the first orthant of the sphere for n = 5 taxa. There are 15  Such computations are pushing the boundary of computational polyhedral geometry. For n ≥ 6 taxa, triangulating a consistency cone is too unwieldy, although we are confident that spherical volumes could still be computed using polynomial time hit-and-run sampling methods for volume approximation [22]. Such methods are complicated and not yet implemented.
Finally, we comment on the example in Section 3 that shows how different alignments may lead to the same neighbor-joining tree, whereas the neighbor-joining tree constructed from a concatenation of the alignments is different. This result has significant implications for studies where species trees are constructed from multiple gene families by combining the data.