Projections of Tropical Fermat-Weber Points

In the tropical projective torus, it is not guaranteed that the projection of a Fermat-Weber point of a given data set is a Fermat-Weber point of the projection of the data set. In this paper, we focus on the projection on the tropical triangle (the three-point tropical convex hull), and we develop one algorithm (Algorithm 1) and its improved version (Algorithm 4), such that for a given data set in the tropical projective torus, these algorithms output a tropical triangle, on which the projection of a Fermat-Weber point of the data set is a Fermat-Weber point of the projection of the data set. We implement these algorithms in R and test how it works with random data sets. The experimental results show that, these algorithms can succeed with a much higher probability than choosing the tropical triangle randomly, the succeed rate of these two algorithms is stable while data sets are changing randomly, and Algorithm 4 can output the results much faster than Algorithm 1 averagely.


INTRODUCTION
In this paper, we study the question: for a given data set X in the tropical projective torus, how to find a tropical polytope C, such that the projection of a Fermat-Weber point of X on C is a Fermat-Weber point of the projection of X on C.
This problem is motivated by the tropical principal component analysis (tropical PCA) proposed in [13,20], which is of great use in the analysis of phylogenetic trees in Phylogenetics. Phylogenetics is a subject that is very powerful for explaining genome evolution, processes of speciation and relationships among species. It offers a great challenge of analysing data sets that consist of phylogenetic trees.
Analysing data sets of phylogenetic trees with a fixed number of leaves is difficult because the space of phylogenetic trees is high dimensional and not Euclidean; it is a union of lower dimensional polyhedra cones in R ( n 2 ) , where n is the number of leaves [13]. Many multivariate statistical procedures have been applied to such data sets [17,4,5,3,7,19]. People also have done a lot of work to apply PCA on data sets that consist of phylogenetic trees. For instance, Nye showed an algorithm [11] to compute the first order principal component over the space of phylogenetic trees. Nye [11] used a two-point convex hull under the CAT(0)-metric as the first order principal component over the Billera-Holmes-Vogtman (BHV) tree space introduced in [2]. However, Lin et al. [8] showed that the three-point convex hull in the BHV tree space can have arbitrarily high dimension, which means that the idea in [11] cannot be generalized to higher order principal components (e.g., see [13]). In addition, Nye et al. [12] used the locus of the weighted Fréchet mean when the weights vary over the k-simplex as the k-th principal component in the BHV tree space, and this approach performed well in simulation studies.
On the other hand, the tropical metric in tree spaces is well-studied [10, Chapter 5] and well-behaved [8]. In 2019, Yoshida et al. [20] defined the tropical PCA under the tropical metric in two ways: the Stiefel tropical linear space of fixed dimension, and the tropical polytope with a fixed number of vertices. Page et al. [13] used tropical polytopes for tropical PCA to visualize data sets of phylogenetic trees, and used Markov Chain Monte Carlo (MCMC) approach to optimally estimate the tropical PCA. Their experimental results [13] showed that, this MCMC method of computing tropical PCA performed well on both simulated data sets and empirical data sets. This paper is motivated by a difference between classical PCA (in Euclidean spaces) and tropical PCA as follows. In classical PCA, the projection of the mean point of a data set X in the Euclidean space is the mean point of the projection of X (e.g., see [21,Page 188]). However, in tropical PCA defined by tropical polytopes, the projection of a tropical mean point (in this paper we call it a Fermat-Weber point) of a data set X is not necessarily a Fermat-Weber point of the projection of X (see Example 2.12). More specifically, it is known that, for a data set X in the Euclidean space, the mean point of X is unique. However, for a data set X in the tropical projective torus (denoted by R n /R1), the Fermat-Weber point of X is not necessarily unique [18,Proposition 20]. For a data set X ⊂ R n /R1 and a tropical convex hull C, the tropical projection [6,Formula 3.3] of the set of Fermat-Weber points of X on C are not exactly equal to the set of Fermat-Weber points of the projection of X on C. In addition, it is also known that, in R n /R1, if a set is the union of X and a Fermat-Weber point of X, then the union has exactly one Fermat-Weber point [9,Lemma 8]. So a natural question is, if a set is the union of X and a Fermat-Weber point of X, can the projection of the Fermat-Weber point of the union be a Fermat-Weber point of the projection of the union? By experiments we know that this is still not guaranteed, and it depends on the choice of the tropical convex hull C (see Example 2.13).
In this paper, we focus on tropical triangles (three-point tropical polytopes). We develop one algorithm (Algorithm 1) and its improved version (Algorithm 4), such that for a given data set X ⊂ R n /R1, these algorithms output a tropical triangle C, on which the projection of a Fermat-Weber point of X is a Fermat-Weber point of the projection of X. By sufficient experiments with random data sets, we show that Algorithm 1 and Algorithm 4 can both succeed with a much higher probability than choosing a tropical triangle C randomly (see Table 1 and Table 2). We also show that the succeed rate of these two algorithms is stable while data sets are changing randomly (see Table 3). Algorithm 4 can output the result much faster than Algorithm 1 does averagely (see Table 4), because in most cases, Algorithm 4 correctly terminates with less steps than Algorithm 1 does (see Figure 5). This paper is organized as follows. In Section 2, we remind readers of the basic definitions in tropical geometry. In Section 3, we prove Theorem 3.6 and Theorem 3.7 for the correctness of the algorithms developed in this paper. In Section 4, we present Algorithm 1 and Algorithm 4. We also explain how the algorithms work by two examples. In Section 5, we apply the algorithms developed in Section 4 on random data sets, and illustrate the experimental results.
For any point u ∈ R n , we define the equivalence class [u] := {u+c·1|c ∈ R}, where 1 = (1, . . . , 1). For instance, the vector (3, 3, 3) is equivalent to (0,0,0). In the rest of this paper, we consider the tropical projective torus R n /R1 := {[u]|u ∈ R n }. For convenience, we simply denote by u its equivalence class instead of [u], and we assume the first coordinate of every point in R n /R1 is 0. Because for any u = (u 1 , . . . , u n ) ∈ R n /R1, it is equivalent to Definition 2.4 (Tropical Distance). For any two points we define the tropical distance d tr (u, v) as: Note that the tropical distance is a metric in R n /R1 [8, Page 2030].  If |X| = 3, then the tropical convex hull of X is called a tropical triangle.
The tropical convex hull tconv(X) is shown in Figure 1. Note that R 3 /R1 is isomorphic to R 2 [16], so the points in Figure 1 are drawn on a plane.
We define the set of tropical Fermat-Weber points of X as The Fermat-Weber point of X is denoted by F X . . , x (t) } ⊂ R n /R1, the set of tropical Fermat-Weber points of X in R n /R1 is a convex polytope in R n−1 . It consists of all optimal solutions y = (y 1 , . . . , y n ) to the linear programming problem: for all 1 ≤ k < ≤ n, and for all i ∈ {1, 2, . . . , t}. 1 , . . . , u (t) n )} ⊂ R n /R1. Also let C = tconv(U ). For any point x = (x 1 , . . . , x n ) ∈ R n /R1, we define the projection of x on C as: where . , x (m) } has exactly one Fermat-Weber point, which is F X .

Examples.
Example 2.12. This example shows that, for a given data set X ⊂ R 3 /R1 and a given two-point tropical polytope C, the projection of a Fermat-Weber point of X on C is not necessarily a Fermat-Weber point of the projection of X on C.
Note that (0, 2, 2) is the unique Fermat-Weber point of P , while the projection of a Fermat-Weber point (0, 3, 3) of X is (0, 3, 2). So we can see that the projection of a Fermat-Weber point of X on C is not a Fermat-Weber point of the projection (see Figure 2). Example 2.13. This example shows that, in R n /R1, if a set X is the union of X and a Fermat-Weber point F X of X, then it is not guaranteed that the projection of the Fermat-Weber point F X of X is a Fermat-Weber point of the projection of X. Besides, whether the projection of the Fermat-Weber point F X of X is a Fermat-Weber point of the projection of X depends on the choice of the tropical convex hull C.
Points in X are red. The projection points of X are blue.

THEOREMS
In this section, we introduce Theorem 3.6 and Theorem 3.7 for proving the correctness of the algorithms developed in the next section.
Lemma 3.1. Suppose we have a data set . Let t be a number which is no more than For any two fixed integers d 1 and d 2 (2 ≤ d 1 < d 2 ≤ n), we define three points u (1) , u (2) , u (3) ∈ R n /R1 as follows.
d 2 are respectively located at the d 1 -th and d 2 -th coordinates of δ C (x (k) ). Proof. Recall that we assume the first coordinate of every point in R n /R1 is 0. For any n ) ∈ X, by Definition 2.10, we have that λ i in (2.4) should be: Then the conclusion follows from (2.4).
Suppose X is the data set stated in Lemma 3.1. For u (1) , u (2) and u (3) in Lemma 3.1, let C = tconv({u (1) , u (2) , u (3) }), we have the following remarks: the equalities (3.2)-(3.4) make sure that the tropical triangle C is big enough; the equalities (3.1) and (3.5) make sure that C parallels with a coordinate plane; the equality (3.5) makes sure that C is located under all points in X. Lemma 3.1 shows that we can project X vertically onto C (see Example 3.2 and Figure 4). Let t = −1. Fix d 1 = 2, and d 2 = 3. By (3.1)-(3.5), we can define three points u (1) , u (2) , and u (3) as: Figure 4). Then, by Lemma 3.1, the projection points of X on C are shown in Figure 4 (see the blue points). Definition 3.3 (Data Matrix). We define any matrix X with n columns as a data matrix, where each row of X is regarded as a point in R n /R1.
Below we denote by X m×n the data matrix X with size m × n.  . For a given data matrix X m×n , and for any two fixed integers d 1 and d 2 (2 ≤ d 1 < d 2 ≤ n), we define the projection matrix of X (denoted by P d 1 ,d 2 (X)) as a matrix with size m × n, such that for all k ∈ {1, . . . , m}, the k-th row of P d 1 ,d 2 (X) is d 2 are respectively the (k, d 1 )-entry and the (k, d 2 )-entry of X, and are respectively located at the (k, d 1 )-entry and the (k, d 2 )-entry of of P d 1 ,d 2 (X).
• t is a fixed number, such that t = min Note that the projection matrix P d 1 ,d 2 (X) is still a data matrix. Recall the Proposition 2.11 tells that, if F X is a Fermat-Weber point of X = {x (1) , . . . , x (m) } ⊂ R n /R1, then F X is the unique Fermat-Weber point of {F X , x (1) , . . . , x (m) }.
Theorem 3.6. Suppose we have a data matrix X (m+1)×n , where the last row of X is a Fermat-Weber point of the matrix made by the first m rows of X. We fix two integers d 1 and d 2 (2 ≤ d 1 < d 2 ≤ n). Let r be the last row of P d 1 ,d 2 (X).
If r is a Fermat-Weber point of P d 1 ,d 2 (X), and u (1) , u (2) and u (3) are defined by (3.1)-(3.5), then the projection of the Fermat-Weber point of X on tconv({u (1) Proof. By Lemma 3.1 and Definition 3.5 we know that, the projection of X on is P d 1 ,d 2 (X). Note that the last row of X is the unique Fermat-Weber point of X. Also note that r is the projection of the last row of X. Then by the assumption that r is a Fermat-Weber point of P d 1 ,d 2 (X) we know that, the projection of the Fermat-Weber point of X on C is a Fermat-Weber point of the projection of X on C.
Theorem 3.7. Suppose we have a data matrix X m×n . We fix two integers d 1 and d 2 (2 ≤ d 1 < d 2 ≤ n). Let r be a point (0, t, . . . , t, r d 1 , t, . . . , t, r d 2 , t, . . . , t), where r d 1 and r d 2 are undetermined numbers, and t is the smallest entry of X. Let f be a Fermat-Weber point of P d 1 ,d 2 (X). If r d 1 = f d 1 , and r d 2 = f d 2 , then r is a Fermat-Weber point of P d 1 ,d 2 (X).
Proof. By Definition 3.5, the i-th row of P d 1 ,d 2 (X) has the form It is easy to see that

ALGORITHMS
In this section, we develop Algorithm 1 and Algorithm 4, such that for a given data set X ⊂ R n /R1, these two algorithms output a tropical triangle C, on which the projection of a Fermat-Weber point of X is a Fermat-Weber point of the projection of X.
The input of Algorithm 1 and Algorithm 4 is a data set There are two main steps in each algorithm as follows.
Step 1. We define a data matrix X, such that for all i ∈ {1, . . . , m}, the i-th row of X is x (i) . We obtain a Fermat-Weber point F X by solving the linear programming (2.3). We define a matrix X with size (m + 1) × n, such that the last row of X is F X , and the first m rows of X come from X.
Step 2. We traverse all pairs (d 1 , d 2 ) such that 2 ≤ d 1 < d 2 ≤ n, and we calculate the projection matrix P d 1 ,d 2 ( X) by Definition 3.5. Check if the last row of P d 1 ,d 2 ( X) is a Fermat-Weber point of P d 1 ,d 2 ( X). If so, we calculate the three points u (1) , u (2) and u (3) by (3.1)-(3.5) in Lemma 3.1, return the output, and terminate. By Theorem 3.6 we know that, the projection of a Fermat-Weber point of X on C = tconv({u (1) , u (2) , u (3) }) is a Fermat-Weber point of the projection of X on C. If for all (d 1 , d 2 ), the last row of P d 1 ,d 2 ( X) is not a Fermat-Weber point of P d 1 ,d 2 ( X), then return FAIL.
Remark 4.1. It is not guaranteed that Algorithm 1 and Algorithm 4 will always succeed (return the tropical triangle). If the algorithms succeed, then by Theorem 3.6, Algorithm 1 is correct, and by Theorem 3.6 and Theorem 3.7, Algorithm 4 is correct.
Algorithm 1 and Algorithm 4 always succeed or fail simultaneously. But our experimental results in the next section show that, Algorithm 1 or Algorithm 4 succeeds with a much higher probability than choosing tropical triangles randomly (see Table 1 and Table 2). Our experimental results also show that, if Algorithm 1 and Algorithm 4 succeed, then with the probability more than 50%, Algorithm 4 would terminate in less traversal steps than Algorithm 1 does (see Figure 5).
Remark that, the difference between Algorithm 1 and Algorithm 4 is the traversal strategy, i.e., the Step 2. is different. Below we give more details about Step 2. Let 1. In Algorithm 1: Step 2.,we traverse all pairs (d 1 , d 2 ) (2 ≤ d 1 < d 2 ≤ n) in L one by one, i.e., we traverse the pairs in the lexicographical order. 2. In Algorithm 4: Step 2., we consider the same L defined in (4.1). Note that |L| = (n−1)(n−2) 2 .
Let W and S be two empty sets. In the future, we will record in W some indices that will be traversed in priority, and record in S the pairs that have been traversed. Let u (1) , u (2) and u (3) be null vectors. Now we start a loop (see lines 7-27 in Algorithm 4). In this loop, we traverse all pairs in L while |S| < (n−1)(n−2) 2 , and u (1) , u (2) and u (3) are null. For each pair (d 1 , d 2 ) ∈ L, if (d 1 , d 2 ) ∈ S, then we skip the pair. If (d 1 , d 2 ) / ∈ S, then we add the pair into S, and calculate the projection matrix P d 1 ,d 2 ( X) by Definition 3.5. Let r be the last row of P d 1 ,d 2 ( X). If r is a Fermat-Weber point of P d 1 ,d 2 ( X), then we calculate u (1) , u (2) and u (3) by formulas (3.1)- (3.5) in Lemma 3.1, return the output and terminate. By Theorem 3.6 we know that, the projection of a Fermat-Weber point of X on C := tconv({u (1) , u (2) , u (3) }) is a Fermat-Weber point of the projection of X on C. If r is not a Fermat-Weber point of P d 1 ,d 2 ( X), then by Theorem 3.7, at most one of the following two equalities holds: where f is a Fermat-Weber point of P d 1 ,d 2 ( X). So we have 3 cases.  2) is similar). Note that W is nonempty at this time, and u (1) , u (2) and u (3) are null. For each element ω ∈ W , we define We start traversing all pairs in L ω . For each pair (ω 1 , ω 2 ) ∈ L ω , if (ω 1 , ω 2 ) ∈ S, then we skip the pair. If (ω 1 , ω 2 ) / ∈ S, then we add the pair into S, and calculate the projection matrix P ω 1 ,ω 2 ( X) by Definition 3.5. Let r be the last row of P ω 1 ,ω 2 ( X). If r is a Fermat-Weber point of P ω 1 ,ω 2 ( X), then calculate u (1) , u (2) and u (3) by formulas (3.1)-(3.5) in Lemma 3.1, output u (1) , u (2) and u (3) , and terminate. If r is not a Fermat-Weber point of P ω 1 ,ω 2 ( X), then by Theorem 3.7, at most one of the following two equalities holds: where f is a Fermat-Weber point of P ω 1 ,ω 2 ( X). So we have 2 cases.
(Case 1.2) If only (4.6) holds, then add ω 2 into W . We move on to the next pair in L ω . If for any pair (ω 1 , ω 2 ) ∈ L, we have (ω 1 , ω 2 ) ∈ S, and the last row of P ω 1 ,ω 2 ( X) is not a Fermat-Weber point of P ω 1 ,ω 2 ( X), then we remove this ω from W . If W becomes empty again, then we continue the traversal of L we paused in (Case 1) (Page 10). If W is still nonempty after one element in W has been removed, then for the next ω ∈ W, we traverse L ω . Now we give two examples to better explain how Algorithm 1 and Algorithm 4 work.
By solving the linear programming (2.3), we obtain a Fermat-Weber point of X, which is F X = (0, −40, 4, 89, 15). Define a matrix X with size (m + 1) × n, such that the last row of X is F X , and the first m rows of X come from X. We have Let L be a list that contains all pairs (d 1 , d 2 )(2 ≤ d 1 < d 2 ≤ 5) in the lexicographical order, that is L = {(2, 3), (2, 4), (2,5), (3,4), (3,5), (4, 5)}. Also let W and S be two empty sets. We will record in W some indices that will be traversed in priority, and record in S the pairs that have been traversed. Now we start the traversal. 1. We first start traversing pairs in L. The first pair in L is (2, 3). Add (2, 3) into S. Note that 2. The next pair in L is (2,4). Add (2, 4) into S. Note that  (2, 3), which is in S already, so we skip it. Similarly we skip (2,4). The third pair in L 2 is (2, 5), which is not in S, so we do the following steps.
Add (2,5) into S. Note that Since for every pair (ω 1 , ω 2 ) ∈ L 2 , (ω 1 , ω 2 ) is in S, and the last row of P ω 1 ,ω 2 ( X) is not a Fermat-Weber point of P ω 1 ,ω 2 ( X), we remove 2 from W . 4. Note that, now W = {5} is nonempty. By (4.4), we have L 5 = {(2, 5), (3,5), (4,5)}. The first pair in L 5 is (2, 5), which is in S already, so we skip it. The second pair in L 5 is (3, 5), which is not in S, so we do the following steps. Add (3,5) into S. Note that Below we give the pseudo code of Algorithm 1 and Algorithm 4. Note that, Algorithm 2 and Algorithm 3 are sub-algorithms of Algorithm 1 and Algorithm 4. For a given data matrix X, Algorithm 2 calculates the summation of tropical distance between the last row of X and each row of X, and also calculates the summation of tropical distance between a Fermat-Weber point of X and each row of X. We will use Algorithm 2 to check if the last row of X is a Fermat-Weber point of X. Algorithm 3 calculates three points u (1) , u (2) and u (3) by (3.1)-(3.5).

IMPLEMENTATION AND EXPERIMENT
We implement Algorithm 1 and Algorithm 4, and test how Algorithm 1 and Algorithm 4 perform. Data matrices, R code and computational results are available online via: https://github.com/ DDDVE/the-Projection-of-Fermat-Weber-Points.git.
Software We implement Algorithm 1 and Algorithm 4 in R (version 4.0.4) [15], where we use the command lp() in the package lpSolve [1] to implement Line 1 in Algorithm 1 and Line 1 in Algorithm 4 for computing a Fermat-Weber point of a data matrix.
In our experiments, we use the command rmvnorm() in the package Rfast [14] to generate data matrices that obey multivariate normal distribution. Hardware and System We use a 3.6 GHz Intel Core i9-9900K processor (64 GB of RAM) under Windows 10. Now we present four tables and one figure to illustrate how Algorithm 1 and Algorithm 4 perform.
1. For a fixed data matrix X m×n , Table 1 shows the proportion of random tropical triangles, on which the projection of a Fermat-Weber point of X is a Fermat-Weber point of the projection of X. From Table 1 we can see that, for a fixed data matrix X, and for random tropical triangles, the "succeed rate" is low. Here, by "succeed rate", we mean the proportion of random tropical triangles, on which the projection of a Fermat-Weber point of X is a Fermat-Weber point of the projection of X. For instance, the highest proportion is 16%, and the lowest proportion is even only 1%. Besides, the succeed rate is extremely low when m and n are both big. (b) We record the proportion by "succeed rate". More specifically, for each pair (m, n), we generate one data matrix X m×n ∼ N (0, diag (10)) and 100 random tropical trian- Here, for all i = 1, 2, 3, we make the first coordinate of u (i) as 0, and all other coordinates of u (i) obey the uniform distribution on [−9999, 9999]. For each triangle C, we test if the projection of a Fermat-Weber point of X m×n on C is a Fermat-Weber point of the projection of X m×n on C.
2. Table 2 shows the succeed rate of Algorithm 1 or Algorithm 4 (recall Remark 4.1 tells that, Algorithm 1 and Algorithm 4 always succeed or fail simultaneously). From Table 1 and Table 2 we can see that, the succeed rates recorded in Table 2 are much higher than those in Table 1. For instance, the lowest rate in Table 2 is 34%, which is still higher than the highest rate in Table 1, and the highest rate in Table 2 is 94%, which is close to 100%. (b) We record the proportion as "succeed rate". More specifically, for each pair (m, n), we generate 100 data matrices X m×n ∼ N (0, diag (10)), run Algorithm 1 or Algorithm 4, and calculate the proportion of that Algorithm 1 or Algorithm 4 succeeds.

3.
We fix m = 120, and we fix n = 20. Table 3 shows how high the succeed rate of Algorithm 1 or Algorithm 4 would be when we change the data matrix X 120×20 . In order to change X, we change v, such that X ∼ N (0, diag(v)). We can see from Table 3 that, when v is changing from 1 to 800, the succeed rate of Algorithm 1 or Algorithm 4 is still around 70%. Note that v is the variance of each coordinate of data points, which means that, when the coordinate of data points fluctuates violently, the succeed rate of Algorithm 1 or Algorithm 4 is still stable. v 1 5 10 50 800 succeed rate 67% 65% 73% 66% 67% TABLE 3. (a) v is a real number such that X 120×20 ∼ N (0, diag(v)). (b) We record the proportion as "succeed rate". More specifically, for each v, we generate 100 random data matrices X 120×20 ∼ N (0, diag(v)), run Algorithm 1 or Algorithm 4, and calculate the proportion of that Algorithm 1 or Algorithm 4 succeeds. Table 4 shows the average computational time for Algorithm 1 and that for Algorithm 4. From Table 4 we can see that, Algorithm 1 and Algorithm 4 are both efficient. For instance, when there are 120 data points, and the dimension of each point is 20, the computational timings of Algorithm 1 and Algorithm 4 are still no more than 7 minutes (373.5734s and 291.9031s). In addition, in most cases, Algorithm 4 takes less time than Algorithm 1 does. For instance, when m is 120, and n is 20, Algorithm 4 takes around one and a half minutes less than Algorithm 1 does.  TABLE 4. (a) m represents the number of data points; n represents the dimension of data points.

4.
(b) We record the average computational time (in seconds) as "time". More specifically, for each pair (m, n), we run Algorithm 1 and Algorithm 4 for 100 random data matrices X m×n ∼ N (0, diag(10)), and record the average computational time for Algorithm 1 and that for Algorithm 4. (c) "A1" means the average computational time of Algorithm 1, and "A4" means the average computational time of Algorithm 4. Figure 5 compares the numbers of traversal steps of Algorithm 1 and Algorithm 4. From Figure 5 we can see that, with the proportion more than 50%, Algorithm 1 always takes more traversal steps than Algorithm 4 does. (e) For each pair (m, n), we run Algorithm 1 and Algorithm 4 with 100 random data matrices X m×n ∼ N (0, diag (10)). If Algorithm 1 and Algorithm 4 correctly terminate, then record the number of traversal steps that Algorithm 1 and Algorithm 4 respectively take.