Tropical medians by transportation

Fermat-Weber points with respect to an asymmetric tropical distance function are studied. It turns out that they correspond to the optimal solutions of a transportation problem. The results are applied to obtain a new method for computing consensus trees in phylogenetics. This method has several desirable properties; e.g., it is Pareto and co-Pareto on rooted triplets.


Introduction
The following optimization problem can be studied in any metric space.Given a finite number of points, sometimes called sites, find a point which minimizes the sum of the distances to the sites.Such a point is called a Fermat-Weber point, and this is some version of a geometric median of the sites, which is known to be robust in a certain sense [7, §21].Computing Fermat-Weber points is a rich topic with a remarkable history; see [7,Chapter II].Here we consider a specific distance function, d △ , which occurs in tropical geometry [2,32].This function is asymmetric, i.e., d △ (a, b) may differ from d △ (b, a).So we call d △ the asymmetric tropical distance, to differentiate from the symmetric tropical distance, which is more common [20, §5.3].The symmetric tropical Fermat-Weber problem was investigated by Lin, Sturmfels, Tang and Yoshida [28] and Lin and Yoshida [29].
As our main result we prove that one (asymmetric tropical) Fermat-Weber point can be computed by solving a transportation problem.The transportation problem is an optimization classic, with numerous applications, both in theory and practice.For an overview we refer to Schrijver's monograph [39, §21.6]; see also the survey by De Loera and Kim for the polyhedral geometry point of view [14].Efficient methods for solving transportation problems include algorithms by Tokuyama and Nakano [42], Kleinschmidt and Schannath [22], and Brenner [8].In general, Fermat-Weber points are not unique, so one part of the present work is devoted to understanding the entire (asymmetric tropical) Fermat-Weber set for a given set of sites.This is tightly related to the study of tropical hyperplane arrangements and tropical convex hulls, which are at the core of tropical combinatorics [20].The latter subject is concerned with the rich interplay between tropical geometry and optimization.
One motivation for studying the Fermat-Weber problem in the setting of tropical geometry comes from phylogenetics [28,29].In that field, a part of computational biology, the goal is to associate trees to input data.A typical example are trees encoding ancestral relations among species, where the data originates from strands of DNA.In tropical geometry spaces of metric trees with n labeled leaves occur naturally as the tropical Grassmannians TGr (2, n); see [31, §4.3] and [20, §10.6].In phylogenetics many different methods are known to construct trees from a fixed data set.Since those methods usually do not come up with the same tree, there is a need to find the common ground.This gives rise to some consensus tree, which describes where the several methods agree.Finding consensus trees is a topic of its own [9], and the authors of [28] argue that "tropical convexity and tropical linear algebra . . .behave better" than other methods.Here we show that passing from the symmetric tropical distance function to its asymmetric sibling leads to a new method for computing metric consensus trees which is even better behaved.This is because the asymmetric tropical Fermat-Weber sets are nicer geometrically.In particular, we show that our approach leads to a consensus method which is regular in the sense of [10].Such a procedure is not known for symmetric tropical distances.For the purpose of finding a consensus method in tree space, there is no immediate disadvantage of employing an asymmetric distance rather than a symmetric one.In fact, asymmetric distances are common in location theory [34].
Our paper is organized as follows.We start out with a brief summary of facts from tropical combinatorics which are relevant for the Fermat-Weber problem.Then we prove that the Fermat-Weber set arises as a cell in the covector decomposition of the tropical torus induced by the sites.That covector cell is then characterized in several ways.The first approach employs regular subdivisions of products of simplices; see [15, §6.2].That leads to a linear programming formulation of the Fermat-Weber problem, and the dual linear program is a transportation problem.The latter then provides efficient algorithms.The final chapter is devoted to computing Fermat-Weber sets in spaces of equidistant trees.In addition to theoretical results we report on computational experiments with polymake [19] and mcf [30].Related work.As an early contribution of tropical geometry to data science Gärtner and Jaggi [18] developed a concept for "tropical support vector machines", with applications to classification in mind.A different train of thought was developed by Pachter and Sturmfels [35, §2.4] who related tropical geometry to phylogenetic trees.Later, Lin et al. [28] connected these ideas to the geometry of tree spaces studied by Billera, Holmes and Vogtman [5].Yoshida, Zhang and Zhang [43] proposed a method to analyze data, which they call "tropical principal component analysis".In a way, the latter may be viewed as a synthesis of the above.A key contribution here is work of Ardila and Klivans, who saw that spaces of equidistant trees form the Bergman fans of the graphic matroids of complete graphs [3].The term "tropical convexity" was coined by Develin and Sturmfels [16] to connect tropical geometry with the older topic of (max, +)-linear algebra [11].

Tropical convexity
The purpose of this section is to set our notation and to collect a few facts which are key to our methods; for the details we refer to [20].We consider the tropical semiring T max = (R∪{−∞}, ⊕, ⊙) with ⊕ = max as the tropical addition and ⊙ = + as the tropical multiplication.The additive neutral element is −∞, and 0 is neutral with respect to the multiplication.Usually, we abbreviate T = T max .The set T n inherits the structure of a semimodule by componentwise tropical addition and tropical scalar multiplication.
A tropical cone in T n is a nonempty subset C which contains each tropical linear combination λ⊙x ⊕ µ⊙y for λ, µ ∈ T and x, y ∈ C. Each tropical cone contains the point −∞1 and the entire set R1, where 1 is the all-ones vector.Therefore it is convenient to pass to the quotient TP n−1 := (T n \ {−∞1})/R1, which is called the tropical projective space.A subset of TP n−1 is tropically convex if it arises as the image of a tropical cone under the canonical projection.A tropical polytope is a finitely generated tropically Figure 1.Five sites in R 3 /R1, their max-tropical convex hull, the induced min-tropical hyperplane arrangement, and the unique Fermat-Weber point (marked white) convex set.The tropical projective torus R n /R1 is the subset of points in TP n−1 with finite coordinates.We say that a tropical polytope is bounded if it lies in R n /R1.
Tropical convexity is intimately related to ordinary convexity, polyhedral geometry and (linear) optimization.For instance, tropical polytopes arise as the images of ordinary convex polytopes over ordered fields of real Puiseux series under the valuation map; see [20,Observation 5.10].Yet, here the following less algebraic description is more relevant.
We consider an arbitrary m×n-matrix V = (v ij ) with real coefficients.Each row v i = (v i1 , . . ., v in ) is a point in R n (or R n /R1, if we ignore the tropical scaling).So V may be viewed as a configuration of m labeled points in R n /R1.Technically, it is convenient to assume that each ordinary row sum equals zero, i.e., each row lies in the set Observe that each point in R n /R1 has a unique representative in R n which lies in H.So we can identify the tropical projective torus R n /R1 with the ordinary linear hyperplane H in R n .This also works topologically, since the quotient vector space R n /R1 is homeomorphic to R n−1 (and thus with H considered as a subset of R n ).The tropical projective space TP n−1 is a compactification of R n /R1, where the boundary comprises those points which have at least one infinite coordinate.
Adding vectors tropically works coefficient-wise, and there is also tropical multiplication by scalars.With these notions, the max-tropical convex hull of (the rows of) V is which is a subset of R n /R1.The rows of the matrix V also define an arrangement of m tropical hyperplanes in R n /R1, with respect to min as the tropical addition, and we denote this T (V ).In this context each row arises as the apex of a min-tropical hyperplane.Everything that we explained above also works for the min-tropical semiring T min = (R∪{∞}, min, +), which is isomorphic to T max as a semiring via x → −x.Observe that this involution leaves the hyperplane H invariant.
The min-tropical hyperplane arrangement with the rows of V as their apices induces an ordinary polyhedral subdivision, CovDec(V ), of R n /R1 (or, equivalently, H), called the covector subdivision induced by (the rows of) V .Its cells, which are called covector cells, are convex in three different senses: they are max-tropically convex, min-tropically convex and convex in the ordinary sense (as subsets of H).Such polyhedra are known as polytropes, and they may be bounded or unbounded.As ordinary polyhedra, the polytropes are characterized by the property that their facet normal directions are e i − e j for i, j ∈ [n] distinct.The union of the bounded covector cells equals the tropical convex hull tconv max (V ).The rows of V (equivalently, the columns of V ⊤ ) encode five points in R 3 /R1; see Fig. 1.
The covector decomposition has 15 regions of maximal dimension 2, and six of them are bounded.The max-tropical convex hull tconv max (V ) consists of the union of bounded cells, which are shaded gray; it also contains the green line segment extending from (3, −3, 0) to the lower left.
We consider the envelope which is an unbounded ordinary polyhedron.By [20,Theorem 6.14] the cells of CovDec(V ) arise as the images of faces of E(V ) under the coordinate projection (t, x) → x.Moreover, CovDec(V ) is dual to the regular subdivision, Σ(V ), of the product of simplices } obtained from lifting (e i , e j ) to the height v ij .Here we take regular subdivisions induced by upper convex hulls since max is our tropical addition.For the same reason the inequality sign "≥" is reversed in comparison with the min-tropical version in [20, (6.1)].Via the Cayley trick the subdivision Σ(V ) of ∆ m−1 × ∆ n−1 corresponds to a mixed subdivision, S(V ), of the dilated simplex m • ∆ n−1 ; see [15, §9.2] and [20, §4.5].For instance, this is convenient for properly visualizing Σ(V ), which is a polyhedral complex of dimension (m − 1)(n − 1).Fig. 2 shows the mixed subdivision S(V ) of 5 • ∆ 2 for the matrix V from Example 1.

Fermat-Weber sets
We examine the Fermat-Weber problem through tropical combinatorics and polyhedral geometry.As our first key observation we show that asymmetric tropical Fermat-Weber sets arise as specific cells in the covector subdivisions induced by the sites.
The asymmetric tropical distance in R n is given by . Throughout this section we fix a finite set of sites V = {v 1 , v 2 , . . ., v m } in H, which we identify with R n /R1.Definition 2. An (asymmetric tropical) Fermat-Weber point with respect to V is a point in H which minimizes the sum of the asymmetric tropical distances from these sites.
In general such a point is not unique.Hence, we let FW(V ) denote the set of all asymmetric tropical Fermat-Weber points and call it the (asymmetric tropical) Fermat-Weber set with respect to V .This is the asymmetric analog of the symmetric tropical Fermat-Weber set studied in [29].
Fixing the site v i ∈ V , the distance function from v i , which reads is convex in the ordinary sense and piecewise linear.Its regions of linearity are precisely the n closed sectors of the min-tropical hyperplane with apex v i ; see [20, §5.5].Consequently, the common subdivision of the regions of linearity of the sites is exactly the covector decomposition CovDec(V ).Our first main theorem shows that the Fermat-Weber set FW(V ) is a bounded cell of that subdivision.
For the sake of a precise formulation of that result, we pass to the regular triangulation

is dual to the covector subdivision CovDec(V ). The relatively open cells of Σ(V ) partition the product of the ordinary polytope ∆
So there is a unique cell, C V , which contains that ( 1 m 1, 1 n 1) in its relatively interior.This will play an important role in our study of Fermat-Weber points.Definition 3. We call the unique cell of Σ(V ) which contains that ( 1 m 1, 1 n 1) in its relatively interior the central cell of Σ(V ), and denote it by C V .Its dual in CovDec(V ) will be called the central covector cell.
Theorem 4. The Fermat-Weber set FW(V ) agrees with the central covector cell in CovDec(V ).In particular, FW(V ) is a bounded polytrope in H, and it is contained in the tropical polytope tconv max (V ).
Proof.Consider the linear program with mn+1 constraints in the m+n variables t 1 , t 2 , . . ., t m , x 1 , x 2 , . . ., x n .The coefficients v ij are the coordinates of the sites.The constant factor n in the objective function is not relevant here, but it does make the dual linear program (5) studied below look more natural.If (t * , x * ) is an optimal solution of (4), then x * ∈ FW(V ), and Conversely, each Fermat-Weber point arises in this way.The constraints v ij ≤ t i + x j describe the max-tropical version of the envelope E(V ); see [20, §6.1].In that reference that inequality would look like "v ij ≤ −t i + x j ".Yet the linear substitution t i → −t i is natural here, because ( 4) is a minimization problem; see also [20,Remark 6.28].The cells of CovDec(V ) are precisely the projections of the faces of the unbounded ordinary polyhedron E(V ) ⊂ R m × R n onto the x-coordinates; see [20,Proposition 6.11].Let F be the optimal face of the linear program (4).The set FW(V ) is the projection of F .
The face F is gotten by minimizing A sublevel set of a function f : R n /R1 → R is a set of the form {x ∈ R n /R1 | f (x) ≤ α} for some α ∈ R. If all of its sublevel sets are bounded, then the set of minima is bounded.We use this property to show that FW(V ) is bounded.
The sublevel sets of d △ (v i , •) are simplices if non-empty; in particular, they are bounded.Consequently, the function x → i∈[m] d △ (v i , x) has bounded sublevel sets.The latter implies that FW(V ) is bounded, as it is the set of minimizers of the aforementioned function.The max-tropical convex hull of the rows of V equals the union of the bounded covector cells in CovDec(V ); see [20,Corollary 6.17].
Example 5.For the matrix V from Example 1 the unique optimal solution of the primal linear program (4) reads t * = (5, 4, 5, 7, 3) and x * = (9, −6, −3) , with optimal value 3 • 24 = 72.We have FW(V ) = {x * }; see Fig. 1.As gcd(5, 3) = 1, the uniqueness is implied by Theorem 7 below.The point x * , which is a pseudovertex of CovDec(V ), is dual to the central cell C V = conv{113, 122, 212, 221}; see Fig. 2. Remark 6. Theorem 4 reveals a similarity to the Euclidean case: by [7,Proposition 19.1] any Fermat-Weber point with respect to the Euclidean distance is contained in the convex hull of the sites.The analogous result to Theorem 4 for the symmetric tropical distance is [28,Proposition 26].As shown in [28,Example 27], in general, the symmetric tropical Fermat-Weber set is not contained in the tropical convex hull.
Via the Cayley trick, the central cell C V in Σ(V ) corresponds to the central covector cell in CovDec(V ), which is FW(V ).The dimension of the latter is the codimension of the former.
Theorem 7. The dimension of FW(V ) is at most gcd(m, n) − 1.In particular, if m and n are relatively prime, the Fermat-Weber point is unique.
Proof.We consider the regular subdivision, Σ(V ), of ∆ m−1 × ∆ n−1 induced by the matrix V .Let C V be the central cell.By [15, §6.2] the vertices of C V are in bijection with the edges in a subgraph G(C V ) of the complete bipartite graph K m,n with c := codim(C V )+1 connected components; see also [20, §4.7].We will show that c ≤ gcd(m, n).Here we may assume that V is generic, whence C V is a simplex; note that any refinement of Σ(V ) can only increase the codimension of the cell containing a specific point in its relative interior.If c = 1, there is nothing to prove.Hence, we assume that c ≥ 2.
We consider a maximal simplex, U, of Σ(V ) which contains C V , and let T be the subtree of K m,n corresponding to U. A facet, F , of U corresponds to the removal of some edge e of T , yielding disjoint unions [m] = I ′ ⊔ I ′′ and [n] = J ′ ⊔ J ′′ such that the connected components of T \ e are the restrictions on I ′ ⊔ J ′ and I ′′ ⊔ J ′′ and e is incident to a vertex in I ′′ and one in J ′ .The hyperplane i ′ ∈I ′ x i ′ + j ′′ ∈J ′′ y j ′′ = 1 contains the facet F .If I ′ were empty, then j ′ ∈J ′ y j ′ = 0 for every (x, y) ∈ C V .But C V contains the point ( 1 m 1, 1 n 1), so J ′ must be empty as well.This contradicts the fact that J ′ contains a vertex incident to e. Hence, I ′ = ∅ and, similarly, J ′′ = ∅.
Consequently, there exist two proper partitions a tree, and there exist edges e 2 , . . ., e c in K m,n between a vertex of I i and J i−1 such that adding these edges we obtain the tree T .
A supporting hyperplane for the facet corresponding to T \{e i } is given by the equation where 2 ≤ i ≤ c.Now the point ( 1 m 1, 1 n 1) is contained in these hyperplanes, yielding the c − 1 equalities Multiplying with the least common multiple of m and n we obtain

Putting those relations together with k∈
Restricting to the one-dimensional case (i.e., n = 2), we recover the known fact that an odd number of points have a unique median, while the median can be selected from an interval for an even number of points.The following example shows that the inequality in Theorem 7 is tight for all m and n; see also Example 1.

Example 8. Our example employs the matrix
The rows are points on the tropical moment curve, and their (max-)tropical convex hull is a tropical cyclic polytope; see [6, §4] and [20,Example 5.18].The dual triangulation Σ(V ) of ∆ m−1 × ∆ n−1 is the staircase triangulation; see [15, §6.2.3].We explain the construction.In this case, we represent the simplices as points in an m×n grid instead of subgraphs of K m,n .The staircase triangulation consists of all the paths in an m×n grid starting at (1, 1) and ending at (m, n) obtained by going right or down.Figure 3 portrays such a path.Due to its shape, it is called a staircase.We will show that in the staircase triangulation the simplex containing Each remaining facet of U corresponds to the removal of a grid point from some block I k ×J k .In particular, each facet induces a partition [m] = I ′ ⊔I ′′ such that Moreover, at least one of the intersections I k ∩ I ′ ∩ I ′′ and J k ∩ J ′ ∩ J ′′ is nonempty, which implies that at least one of |I ′ | and |J ′′ | is not a multiple of d.As in the proof of Theorem 7, the facet defining equation The staircase for U is obtained by using the Northwest Corner Rule with breaking ties by going East.If we break the ties randomly, then 2 d−1 staircases appear with nonzero probability.These staircases are in bijection with the ordinary vertices of FW(V ), which is a (d−1)-dimensional cube, seen as an ordinary polytope.
Remark 9.In the special case m = n−1 computing tropical Fermat-Weber points reduces to the tropical Cramer rule [20, §4.9].Adding one more site, p, the new Fermat-Weber set, F , can have higher dimension, but contains the tropical Cramer point, c.Perturbing c in the direction of p, we obtain a point in the relative interior of F .This also agrees with results of Gärtner and Jaggi [18, §4.1] in the context of "tropical support vector machines" on computing a "separating hyperplane for n points".
Remark 10.As a consequence of [20, Theorem 6.14] the Fermat-Weber set FW(V ⊤ ) in R m /R1 of the n columns of V is affinely isomorphic to the Fermat-Weber set FW(V ) in R n /R1 of the m rows.
Remark 11.Our analysis rests on the decision, in (3), to look at the distances from the sites to the Fermat-Weber points.This leads to min-tropical hyperplane arrangements and max-tropical convex hulls.Reversing the direction, i.e., considering distances from the Fermat-Weber points to the sites, amounts to exchanging min and max throughout.Conceptually, the results remain the same; cf.[20, §1.3].

Transportation
This section comprises the algorithmic core of this paper.The basic ingredient is the transportation problem which already occurred in Example 8.
Again let us fix a matrix V = (v ij ) ∈ R m×n .Whenever it will suit us, we may also view the rows of V is a m labeled points v 1 , . . ., v m in H or R n /R1.The following is the dual of the linear program (4) with variables λ and y ij for i ∈ [m] and j ∈ [n]: As the right hand sides are the integral constants m and n, the feasible region of ( 5) is a central transportation polytope, and we denote it T (m, n).The polytope T (m, n) is known to be integral [14,Lemma 2.13].The nonzero entries of any vertex, which is an m×n-matrix, defines a subgraph of K m,n by picking edges [14,Lemma 2.9].This is the support graph of that vertex, and this is a forest; see [39,Theorem 21.15].
Example 12.For the 5×3-matrix V from Example 1 the transpose of the unique optimal solution of ( 5) is The support graph a spanning tree; see Fig. 4. By Theorem 15 below that tree encodes the covector of the unique Fermat-Weber point x * .The degree sequence for the column nodes is (223).This is the coarse type of x * and also the componentwise maximum of the four vertices of the central cell in the mixed subdivision; cf.Example 5 and [20, §4.5].We borrow some of their ideas.For J ⊆ [n] and u ∈ H consider the set We have S ∅ (u) = ∅, and S [n] (u) = H, where we use the convention that the maximum of the empty set is −∞, the neutral element of the tropical addition.If J is a nonempty proper subset of [n], then S J (u) is a max-tropical halfspace with apex u; see [20, §7.1].In the special case where J = {j} is a singleton that tropical halfspace is a closed sector; in general, S J (u) = j∈J S {j} (u).Definition 13.We say that u ∈ H evenly splits V if for every subset J of [n] we have n Tokuyama and Nakano [42] call the point u a "1-splitter", and the sectors are the "regions split by u".Theorem 15 below may be seen as our interpretation of their results in the setting of tropical convexity.
Example 14.Consider the points V from Example 1.The Fermat-Weber point u = (9, −6, −3) evenly splits them.This is illustrated in Figure 5, where we have drawn the max-tropical hyperplane based at u with dotted lines.In particular, we see the subdivision of R 3 /R1 in three convex regions.
Theorem 15.A point u ∈ H evenly splits V if and only if u ∈ FW(V ).The support graph of an optimal dual solution y * is the covector of the Fermat-Weber point u(y * ).
Proof.Assume that u evenly splits V and consider Thus (u, t) is a feasible solution of the primal linear program (4).Moreover, Now [42, Theorem 2.2] says that there exists a solution (y ij ) of the transportation problem (5), which is dual to (4), such that (u, t) and (y ij ) satisfy the complementary slackness condition.Indeed, if t i + u j = v ij , then v i / ∈ S {j} (u), so the aforementioned result gives y ij = 0.So it follows from linear programming duality that (u, t) is an optimal solution of (4).In particular, u is a Fermat-Weber point for V .
For the converse, we denote by φ the convex function 1 n s∈V d △ (s, •).Also, for every subset J of [n], denote by σ J the cardinality of V ∩ S J (u).Abbreviating f J := (n − |J|) j∈J e j − |J| i∈[n]\J e i , which is a point in H, we obtain: • if s ∈ S J (u), then s ∈ S J (u − δf J ) for any δ ≥ 0; • if s / ∈ S J (u), then s / ∈ S J (u − δf J ) for any δ ≥ 0 sufficiently small.The condition s ∈ S J (u) implies d △ (s, u) = n(s j − u j ) for some j ∈ J. Therefore, for δ > 0 sufficiently small, we obtain . Summing up and dividing by n yields (7) φ(u for any δ > 0 sufficiently small.If u ∈ FW(V ), then u is a minimizer, so φ(u−δf J ) ≥ φ(u) for any δ > 0 and J ⊆ [n].Equation ( 7) implies nσ J ≥ m|J| for every J ⊆ [n], under this assumption.Consequently, u evenly splits V .
The argument in the proof of Theorem 15 leads to the following algorithm for obtaining a tropical Fermat-Weber point from an optimal solution y * of (5).By complementary slackness each edge (i, j) of the support graph, B(y * ), imposes the equality (8) for any optimal dual solution.If we assume x j = 0 for some column node j, then we can perform a depth-first search from j and recover all the the other values of (t, x) using the equalities in (8).There may be more connected components, in which case we start a depth-first search from every unvisited column node.In this way all the values are recovered, as B(y * ) is spanning-no row or column of y * can be zero as its elements must sum up to a positive number.Note that the point x obtained this way may not lie in H. Yet, by adding (x 1 + • • • + x n )/n to every entry of x and subtracting the same value from every entry of t, we obtain a feasible solution (t * , x * ) of ( 4) that satisfies the equations (8).In particular, that solution is optimal, whence x * ∈ FW(V ).So the method of [42] gives the following complexity result.
Corollary 16.Assuming m ≥ n, one tropical Fermat-Weber point can be computed in O(n 2 m log 2 m) time.
The complexity bounds of the algorithms by Kleinschmidt and Schannath [22] and Brenner [8] are similar, but slightly different.Naturally, they carry over as well.
It is also natural to ask for explicit representations of the entire Fermat-Weber set.A first idea could be to list the vertices of FW(V ), seen as an ordinary polytope.Example 8 shows that cubes occur as Fermat-Weber sets, and their number of vertices is exponential in the dimension.Yet, the polytropal structure allows for more efficient choices.Namely, a polytope in R n /R1 has at most n tropical vertices and at most n 2 − n ordinary facets.
Let us start with the latter.Since we know the possible directions of the (outer) facet normal vectors, e k − e ℓ for k = ℓ, we can find the ordinary facets by solving n 2 − n linear programs like: where p * is the optimal value of (4).The constraint matrix has only 0 and ±1 entries, and so a linear program of the form ( 9) can be solved in strongly polynomial time [40,Corollary 15.3a].
Each such linear program yields one tight inequality , where a k,ℓ is the optimal value.Then, the tropical vertices can be found as , where a k,k = 0. From Corollary 16 we thus infer the following result.

Tropical median consensus trees
Phylogenetics is a branch of (computational) biology that seeks to associate trees to mark ancestral relations among given taxa, e.g., species [41].The taxa correspond to the leaves.Here we show how asymmetric tropical Fermat-Weber sets give rise to a new algorithm for finding consensus trees.Since there is no particular shortage on consensus tree methods, we compare its features to other approaches, and we discuss the practical applicability.
holds for all i, j, k ∈ [n].Since the (zero) diagonal is implicit, and the matrix is required to be symmetric, we may view a dissimilarity map as an element of R n 2 .A rooted metric tree with n labeled leaves is equidistant if the distance from any leaf to the root is the same.It is known that the ultrametrics are precisely the distance functions among the leaves in an equidistant tree; see [41,Theorem 7.2.5].Note that the all-ones vector 1 of length n 2 is an ultrametric.The corresponding equidistant tree has interior edges of length zero, while each leaf edge has length 1  2 .Ardila and Klivans [3,Theorem 3] showed that a dissimilarity map is an ultrametric if and only if it corresponds to a point on the Bergman fan B(K n ) of the complete graph K n .The Bergman fan of a matroid is a special case of a tropical linear space, i.e., with constant coefficients.The ultrametric inequality ( 10) is a tropical convexity condition with respect to max.From a dissimilarity map D and any real constant c we can construct the dissimilarity map D + c1.Moreover, if D is an ultrametric, and D + c1 is nonnegative, then it is an ultrametric, too.In this way, we may view B(K n ) as a maxtropical linear space in the tropical projective torus R n 2 /R1; see [28,Proposition 16] and [43,Theorem 3].We abbreviate T n := B(K n )/R1 and call it the space of equidistant trees on n labeled leaves.Now we can apply our results from Sections 3 and 4 to points in T n .Our first observation says that Fermat-Weber points of equidistant trees are again equidistant trees.
Theorem 18.Let V ⊂ T n be a finite set of equidistant trees on n leaves.Then the tropical polytope FW(V ) is contained in T n .Moreover, any two trees in FW(V ) share the same tree topology.
Proof.According to Theorem 4 the set FW(V ) is contained in the max-tropical convex hull of V .The space of equidistant trees T n is a max-tropical linear space and thus maxtropically convex; see [20,Proposition 10.33].This is the first claim.Page, Yoshida and Zhang showed [36, Theorem 3.2] that the trees in any relatively open cell of the covector decomposition of V in T n share the same tree topology.With this observation the second claim follows also from Theorem 4.
As we know, Fermat-Weber points are not unique, in general.Here is a more precise statement.
Corollary 19.Let V ⊂ T n be a set of m equidistant trees on n leaves.Then dim FW(V ) ≤ min n − 1 , gcd(m, n 2 ) − 1 .Proof.This follows from Theorem 7 and the fact that the graphic matroid of K n has rank n − 1, so dim T n = n − 2.

Consensus trees.
Our goal now is to employ the results obtained so far to study the consensus problem for metric trees.Formally, a consensus method on equidistant trees, with n taxa, is a function c : T * n → T n where T * n = m≥1 T m n .For surveys on the subject see [9] and, for a more recent account, [10].We say that a consensus method c is tropically convex if c(D 1 , . . ., D m ) ∈ tconv max (D 1 , . . ., D m ) for every m ≥ 1 and every D 1 , . . ., D m ∈ T n .Theorem 18 says that selecting an arbitrary tree in FW(D 1 , . . ., D m ) yields a consensus method which is tropically convex.Recall that FW(D 1 , . . ., D m ) is a polytrope in R n 2 /R1, which thus has at most n 2 tropical vertices.Definition 20.We define the tropical median consensus tree of D 1 , . . ., D m ∈ T n as the ordinary average of the tropical vertices of FW(D 1 , . . ., D m ).
That ordinary average is the barycenter of the ordinary simplex spanned by the tropical vertices.As polytropes are convex in the ordinary sense, the tropical median consensus tree method is tropically convex.By Corollary 17, the tropical vertices of the Fermat-Weber set can be found in strongly polynomial time and hence also their ordinary average.The tropical median consensus tree is depicted in Fig. 6 as (d).This is the unique Fermat-Weber point of the three input trees in T 9 .This can be verified using version 4.9 of polymake [19].
The Newick format is a standard for encoding phylogenetic trees [17, p. 590]; it is supported by polymake.A tree is represented as a string, where leaves are given by their (text) labels, and internal nodes correspond to matching pairs of parentheses.Recursively, such a pair of parentheses encloses the Newick representation of the subtree rooted at that internal node.Further, each node is followed by a numerical value, after a colon, and this denotes the length of the edge between the node and its parent.For example, the Newick representation of (a) is (A:8,((B:2,(C:1,D:1):1):5,((E:1,F:1):3,(G:2,(H:1,I:1):1):2):3):1).
Bryant, Francis and Steel [10] impose three conditions for a consensus method to be regular : (U) The consensus of any number of copies of the same tree, T , is T ; Whatever convention we may fix for the representatives, the pointwise maximum exhibits a drawback.To exemplify it, consider the trees (a), (b), and one million copies of the tree (c) from Fig. 6.Then, the pointwise maximum consensus is still the one from Fig. 7, whereas the tropical median consensus tree looks like (c).So the pointwise maximum consensus is highly sensitive to outliers.In contrast, the tropical median consensus is robust.
Most known consensus tree methods deal with unweighted phylogenetic trees, so they may be seen as discrete analogues of our approach.For instance, unweighted "median consensus trees" were defined by Bathélémy and McMorris [4], and the asymmetric case was analyzed by Phillips and Warnow [38].Bryant [9] presents only one consensus method for rooted trees that takes the edge lengths into consideration: the "average consensus tree" of Lapointe and Cucumel [25].In [9, §2.4.2] two drawbacks of the average consensus tree method are explicitly mentioned: no efficient algorithm is known, and the (co-)Pareto properties are unclear.The average consensus method involves the Euclidean distance, and the unconstrained optima might lie outside the tree space T n .Therefore, ultrametric conditions must be imposed to the solution, making it difficult to obtain a regular consensus method; see [25] for details.Further, Lapointe and Cucumel [24] show that the procedure proposed is NP-hard, and the solution may not be unique.A similar complexity result exists for the median consensus method developed by Lavasseur and Lapointe [26]; see [13].
In the remainder of this section we compare the tropical median consensus method to algorithms proposed by Lin and Yoshida [29].In fact, our approach is very similar and was inspired by that article.Lin and co-authors [28] studied the Fermat-Weber problem for the symmetric tropical distance function, with a focus on the tree space.Crucially, the symmetric tropical Fermat-Weber points may lie outside the max-tropical convex hull; see [28,Example 27].The symmetric Fermat-Weber set of T 1 and T 2 , denoted FW sym (T 1 , T 2 ), contains the tropical segment between them, which exhibits seven distinct tree topologies, four of which are binary.In contrast, the asymmetric setting is trivial: the tropical median is (A : 10, D : 10, (B : 4, C : 4) : 6), and this is the unique asymmetric tropical Fermat-Weber point.
In view of [36,Lemma 3.5] the symmetric tropical distance function might lead to a robust tropical consensus method.However, examples like the above form a challenge to defining a method which is globally consistent.The next case shows more differences between the symmetric and the asymmetric distances.
Remark 26.The max-tropical convexity of the tropical median consensus method ultimately rests on the specific formulation of the Fermat-Weber problem in (3).Exchanging the arguments in the asymmetric tropical distance function gives the min-tropical analog.5.3.Computational experiments.We compare running times for experiments concerning Fermat-Weber sets in tree space with respect to the symmetric and the asymmetric tropical distance functions.As input data we take random trees which were produced using the function rmtree from the R library ape [37].This is similar to the experiment reported in [43,Example 8].Most of the trees generated are not equidistant, so we adjust the lengths of the leaf edges to make them equidistant.In this way we get any number of trees in T n , for various values of n, the number of taxa.Recall that, by [28,Proposition 26], the symmetric tropical Fermat-Weber set FW sym (T 1 , . . ., T m ) of trees T i ∈ T n is a convex polytope in R n 2 /R1.The asymmetric tropical Fermat-Weber set FW(T 1 , . . ., T m ) is a polytrope, and thus a polytope, too; cf.Theorem 4. All timings are obtained with polymake, version 4.9, running on a quad core Intel Core i5-4590 processor (6599.89While the trees are generated with edge lengths given by floating point numbers, we convert them to exact rationals.We are not aware of a published implementation for computing symmetric tropical Fermat-Weber sets, so we implemented it in polymake. The algorithm suggested by [28,Proposition 26] allows for an improvement by exploiting properties of polyhedral L-convex functions in the sense of Murota [33,§7.8].In this way, a facet description of FW sym (T 1 , . . ., T m ) can be obtained from solving n(n − 1) linear programs similar to those in (9).The timings for up to ten trees on up to ten leaves are given in Table 1.Lin and Yoshida report about computations of symmetric tropical Fermat-Weber sets in [29, §4]; however, they do not give timings.The parameters leading to [29, Table 1] are much smaller than the parameters in our Table 1.The combinatorial description of the asymmetric tropical Fermat-Weber set FW(T 1 , . . ., T m ) via optimal dual variables of a transportation problem allows us to exploit complementary slackness.So we can compute the description of FW(T 1 , . . ., T m ) in terms of its ordinary facets much faster than its symmetric counterpart; see Table 2. Tropical median consensus trees.Our most interesting experiment is concerned with computing tropical median consensus trees of m trees in T n , again for various values of m and n.This time we use mcf [30] (through our polymake interface), which is a standard implementation of the network simplex algorithm, using floating point computations.Table 3 has the timings for up to 25 trees on up to 300 leaves.
The last row of Table 3, for m = 25 trees, is particularly interesting.Observe that the timings in that row, with an increasing number of leaves, are not monotone.The most likely explanation comes from Corollary 19, which gives an upper bound for the dimension of FW(T 1 , . . ., T m ); the critical contribution is gcd(m, n 2 ).We have 25  2 = 300, and the running times are almost proportional to gcd(m, 300).
In order to see what is going on, we conducted a more refined experiment for n = 25 taxa, trying all values of m in the set {1, 2, . . ., 299}.The results are shown in Fig. 8.The timing for m = 300, being more than 380 seconds, has been omitted for better readability of the other results.There are only ten computations which last more than four seconds.These are the values m ∈ {50, 60, 75, 100, 120, 150, 200, 225, 240, 300}; all are integers with gcd(m, 300) ≥ 50.A more detailed analysis is beyond the scope of the present article.
To the best of our knowledge, no regular consensus method arising from the symmetric tropical distance function on tree space has been proposed.So there is no direct way to make a comparison.Via simulated annealing, the latter authors generate three trees whose tropical convex hull fits best the input data; then they project the input trees on this tropical triangle and display the result in [36,Fig. 6].So this line of research is more about dimension reduction techniques for analyzing data rather than obtaining location statistics; see also [43].In particular, this does not seem to lead to a regular consensus method in the sense of [10].Nonetheless, our method applies to their data.The trees discussed in [36, §6.2] have been reconstructed by Kuo, Wares and Kissinger [23] from 268 orthologous sequences with eight species of protozoa.Seven species among the taxa are Babesia bovis (Bb), Cryptosporidium parvum (Cp), Eimeria tenella (Et), Plasmodium falciparum (Pf), Plasmodium vivax (Pv), Theileria annulata (Ta) and Toxoplasma gondii (Tg).The eighth taxon is Tetrahymena thermophila (Tt), which forms Tt Cp Et Tg Pf Pv Bb Ta Figure 9.The tropical median consensus of the apicomplexa data the outgroup.For this data we obtain the tropical median consensus tree (Cp:0.570333,Et:0.570333,Tg:0.570333,Tt:0.570333,(Pf:0.43862,Pv:0.43862):0.131713,(Bb:0.57033,Ta:0.57033):0.000003), which is displayed in Fig. 9.Only two cherries are resolved; the outgroup was not detected.So our method is quite conservative, making an effort to avoid false positive results.This may be an advantage, in particular since the 268 input trees were generated by a diverse range of methods.

Conclusion
There is a considerable amount of work on the tropical metric geometry of the space of equidistant trees [28,29,43,36,27].This is quite natural, because the symmetric tropical distance between two equidistant trees can be interpreted naturally, as the cost of changing one tree into the other along a tropical line segment, where the cost is measured in the ℓ ∞ norm; see [27] for details.Here we study the asymmetric tropical distance, which has a similar interpretation, but for the ℓ 1 norm.What makes the asymmetric tropical distance attractive is the fact that it leads to a regular consensus method, while this is not known to exist in the symmetric case.Further, the tropical median consensus method is robust and fast in practice; it can be applied to hundreds of trees with dozens of taxa.The robustness is computationally valuable as it makes floating-point computations reliable.
The tropical median consensus method seems to be rather conservative, in the sense that our consensus trees sometimes show little resolution.This is not necessarily bad.For instance, the apicomplexa gene trees studied in Section 5 come from a diverse mix of tree building methods, which should probably make it difficult to find a fully resolved consensus; cf.Fig. 9.

Figure 2 .
Figure 2. Mixed subdivision S(V ) of 5•∆ 2 for V as in Example 1.The 21 lattice points of 5 • ∆ 2 are marked with their coordinates.The min-tropical hyperplane arrangement T (V ) admits a piecewise-linear embedding

( 1 m 1, 1 n 1 )m 1 and marginal row 1 n 1 . 1 , 1 n 1 ) 1 , 1 n 1 )
lies in a cell of codimension gcd(m, n) − 1.For improved readability, we abbreviate d = gcd(m, n), a = m/d, and b = n/d.If d = 1, then there is a unique maximal simplex in the staircase triangulation containing ( 1 m 1, 1 n 1) in its interior.To find the precise staircase, we use the Northwest Corner Rule [12, §8.3.1] in the standard transportation array with marginal column 1 The visited cells form the staircase, which we call the central staircase.Note that this method provides also barycentric coordinates for ( 1 m in this simplex.When d ≥ 2, consider the partitions in d subsets [m] = I 1 ⊔• • •⊔I d and [n] = J 1 ⊔• • •⊔J d where I i = {(i − 1)a + 1, (i − 1)a + 2, . . ., ia} and J i = {(i − 1)b + 1, (i − 1)b + 2, . . ., ib}.Consider on I i × J i the central staircase on an a × b grid and add the points (ia, ib + 1) for i = 1, . . ., d − 1: in Figure3the blocks on I i × J i correspond to the gray areas whereas the added points are those in the white area.This staircase corresponds to a maximal simplex U in the staircase triangulation, which contains ( 1 m 1, 1 n 1).The removal of the grid point (ia, ib + 1) yields the facet defined by k≤ia x k + ℓ>ib y ℓ = 1 .Using d = m/a = n/b, we see that ( 1 m is contained in this facet.In total, there are d − 1 facets of this form. we get λ = −m for every feasible point.By substituting that value in(5) we obtain the standard linear programming formulation of a transportation problem; see[39,  §21.6].The primal linear program (4) and the envelope E(V ) are related to transportation via Hitchcock's theorem[39, Theorem 21.13].

Figure 4 .( 3 ,Figure 5 .
Figure 4. Spanning subtree of K 5,3 encoding the covector of the unique Fermat-Weber point from Example 1. Column nodes at the bottom

Corollary 17 .
The ordinary facet description and the tropical vertices of FW(V ) can be found in strongly polynomial time.If the tropical vertices are known, then we can check if a given point lies in FW(V ) in O(n 2 ) time by checking the criterion[20, Proposition 5.37].

5. 1 .
Equidistant trees.We recall known facts about ultrametrics and equidistant trees to fix our notation; see [41, §7.2] and [20, §10.9].A dissimilarity map is a symmetric n×n-matrix D = (d ij ) with zero diagonal.It is called an ultrametric if D is nonnegative, and the ultrametric inequality(10)

Figure 7 .
Figure 7. Pointwise maximum of the three trees from Example 23

Figure 8 .
Figure 8. Number of trees vs. time in seconds for up to 299 trees on 25 leaves this induces a directed distance function in the (n−1)-dimensional quotient vector space R n /R1.We do not distinguish between the function d △ : R n × R n → R ≥0 and the induced function on R n /R1.The asymmetric tropical distance is a "polyhedral norm" with respect to the standard simplex ∆ n−1 := conv{e 1 , ..., e n } + R1 in R n /R1; see[7,  §20].

Table 2 .
Timings (in seconds) for computing asymmetric tropical Fermat-Weber sets of equidistant trees.Each entry is the average running time from 100 individual experiments, which show only very small variance Linux 6.1.0).For details see our data repository at https://github.com/micjoswig/TropicalDataAnalysis/tree/main/Tropical_medians_by_transpEntire Fermat-Weber sets.First, we compute exact facet descriptions of the polytopes FW sym (T 1 , . . ., T m ) and FW(T 1 , . . ., T m ), where T i ∈ T n , for various values of m and n.