Maximum Likelihood Estimation of Symmetric Group-Based Models via Numerical Algebraic Geometry

Phylogenetic models admit polynomial parametrization maps in terms of the root distribution and transition probabilities along the edges of the phylogenetic tree. For symmetric continuous-time group-based models, Matsen studied the polynomial inequalities that characterize the joint probabilities in the image of these parametrizations (Matsen in IEEE/ACM Trans Comput Biol Bioinform 6:89–95, 2009). We employ this description for maximum likelihood estimation via numerical algebraic geometry. In particular, we explore an example where the maximum likelihood estimate does not exist, which would be difficult to discover without using algebraic methods.


Introduction
A phylogenetic tree is a rooted tree that depicts evolutionary relationships between species. A phylogenetic model is a statistical model describing the evolution of species on a phylogenetic tree. There is a discrete random variable associated with every vertex of the tree. The random variables associated with interior vertices are hidden and correspond to extinct species; the random variables associated with leaves are observed and correspond to extant species. The model parameters are the root distribution and the rate or transition matrices at the edges of the phylogenetic tree. There are different constraints on the model parameters depending on the phylogenetic model. The joint probabilities of random variables associated with leaves (leaf probabilities) are polynomials in the model parameters. Cavender and Felsenstein (1987), and, separately, Lake (1987), introduced an algebraic approach to study phylogenetic models focusing on the search for phylogenetic invariants. A phylogenetic invariant of the model is a polynomial in the leaf probabilities which vanishes for every choice of model parameters. However, phylogenetic invariants alone do not describe the image of the parametrization map. One needs to include inequalities in order to obtain a complete description of the set of leaf probabilities corresponding to phylogenetic tree models.
This paper focuses on the study of continuous-time group-based models. In the rest of the paper, a phylogenetic model is always continuous-time unless written otherwise. Transition matrices of continuous-time phylogenetic models come from continuoustime Markov processes and they are matrix exponentials of rate matrices. Rate matrices of group-based models have a special structure that is determined by an abelian group. A symmetric group-based model assumes that the rate matrices along every edge are symmetric. In particular, a symmetric group-based model can be a submodel of a nonsymmetric group-based model with extra symmetricity conditions on rate matrices. The precise definitions are given in Sect. 2.
Generating sets for phylogenetic invariants for group-based models are described in Sturmfels and Sullivant (2005), Casanellas et al. (2015). These papers consider discrete-time group-based models that require transition matrices to have a special structure determined by an abelian group, but they do not require transition matrices to be matrix exponentials of rate matrices. Generating sets derived in these papers are also valid under the continuous-time approach. However, inequalities defining both models differ, because the set of transition matrices is smaller under the continuous-time approach. A method for deriving the inequalities under the continuous-time approach is given in Matsen (Matsen 2009, Proposition 3.5). We will explicitly derive the semialgebraic description of the leaf probabilities of the CFN model on the tripod tree K 1,3 .
Identifying the equation and inequality characterization of the leaf probabilities is only one part of the problem. The maximum likelihood estimation aims to find parameters that maximize the likelihood of observing the data for the given phylogenetic tree and phylogenetic model. Estimating the tree topology is another part of phylogenetic inference not considered here, see for example Dhar and Minin (2016) for a general overview on phylogenetic inference. Standard methods for the maximum likelihood estimation of the model parameters are the Newton-Raphson method (Schadt et al. 1998;Kenney and Gu 2012), quasi-Newton methods Olsen et al. (1994) and the EM algorithm (Felsenstein 1981;Friedman et al. 2002;Holmes and Rubin 2002;Hobolth and Jensen 2005). It is shown in Steel (1994), Chor et al. (2000) that likelihood functions on phylogenetic trees can have multiple local and global maxima, and thus none of the above methods can guarantee finding the global MLE as these methods are hillclimbing methods. It is stated in Dhar and Minin (2016) that currently no optimization method can guarantee to solve the optimization of the likelihood function over model parameters.
We suggest an alternative method that theoretically gives the solution to the maximum likelihood estimation problem with probability one. This method is based on numerical algebraic geometry (Sommese and Wampler 2005;Bates et al. 2013). The main idea behind this method is to use a numerical algebraic geometry package to compute all critical points of a likelihood function and then choose the critical point with the highest likelihood value. A similar method has been previously applied in optimal control (Rostalski et al. 2011) and in the life sciences (Gross et al. 2016).
Since phylogenetic models are not necessarily compact, the MLE might not even exist. We will use the proposed method to study an example for which the MLE does not exist for the CFN model on the tripod K 1,3 and a particular data vector. In this example, the global maximum is achieved when one of the model parameters goes to infinity. The nonexistence of the MLE would be very difficult to discover without the algebraic methods that we use in this paper, because standard numerical solvers output a solution close to the boundary of the model as we will demonstrate by solving the same MLE problem in Mathematica. One should see the example for the CFN model on the tripod K 1,3 as an illustration of a concept. It will be the subject of future work to develop a package that automatizes the computation in the phylogenetics setting, so that it can be easily used for studying further examples.
In Sect. 2, we introduce the preliminaries of phylogenetic models and present tools from Matsen (2009). Based on Matsen (2009), we state in Sect. 3 Proposition 3 that gives an algorithm for deriving the semialgebraic description of the leaf probabilities of a symmetric group-based model. A proof of Proposition 3 is given in "Appendix A". Algorithm 1 in Sect. 4 outlines how to use numerical algebraic geometry to theoretically give the MLE with probability one. This algorithm is applied on the CFN model on the tripod in Example 5.

Preliminaries of Group-Based Models
The exposition in this section largely follows Matsen (2009). A phylogenetic tree T is a rooted tree with n labeled leaves and it represents the evolutionary relationship between different species. Its leaves correspond to current species and the internal nodes correspond to common ancestors. There is a discrete random variable X v taking k ∈ N possible values associated to each vertex v of the tree T . Typical values for k are two, four or twenty, corresponding to a binary feature, the number of nucleotides and the number of amino acids. For example, if k = 4, the random variable at a leaf represents the probability of observing A, C, G or T in the DNA of the species corresponding to the leaf.
A phylogenetic model assumes a collection of random variables under a Markov process (see Norris (1998) for a detailed introduction on Markov chains). The Markov process on the tree is determined entirely by the probability distribution at the root and the transition matrices P (e) associated to every edge e that reflect the change in the probabilities when moving from one vertex to another. The transition matrices have the form where exp stands for matrix exponentiation, t e ≥ 0 represents time and Q (e) is a rate matrix. The non-diagonal entries of a rate matrix are nonnegative and each row sums to zero. In the rest of the paper, we assume that t e is incorporated in the rate matrix Q (e) .
To define a group-based phylogenetic model, we first fix an abelian group G, a finite set of labels L and a labeling function L : Hence transition matrices of the group-based model form a subset of all the transition matrices that satisfy P . This is because the matrix exponentiation is defined as e M = ∞ i=0 1 i! M i and if a matrix M has the structure given by G, L and L, then one can check that also M i has the structure given by G, L and L for all i ∈ N. The phylogenetic models we consider are symmetric, which means Q h,g . In the case of group-based models, this is equivalent to L(g) = L(−g) for all g ∈ G.
We will assume that the root distribution π of a group-based model is uniform or the root distribution π is such that the matrix P ∈ R G×G defined by P g,h := π(h − g) is a transition matrix in the group-based model (i.e., it is exponential of a rate matrix in the group-based model). In the latter case, we add a new edge starting from the root and re-root the tree at the additional leaf. Instead of the previous root distribution, we use a new root distribution that puts all the mass at the identity and a new transition matrix which is the transition matrix P defined above. We will consider the new leaf as a hidden vertex while other leaves are considered as observed vertices. The same rerooting procedure is used in Sturmfels and Sullivant (2005), Matsen (2009). This approach does not allow completely arbitrary root distributions. In particular, a root distribution has to satisfy π(g 1 ) = π(g 2 ) whenever L(g 1 ) = L(g 2 ) and it has to satisfy inequalities that guarantee that the transition matrix P defined by P g,h := π(h − g) is a matrix exponential of a rate matrix. The latter problem is called the embedding problem and is studied for 2 × 2 matrices in Kingman (1962) and for the Kimura 3-parameter model in Roca-Lacostena and Fernández-Sánchez (2017). In (Sturmfels and Sullivant (2005), Section 6), a workaround is described for deriving phylogenetic invariants for arbitrary root distributions for discrete-time group-based models. We will describe a workaround for deriving inequalities describing the CFN model for arbitrary root distributions; however, we do not know how to generalize this approach to other models.
The joint probability distributions p i 1 ,...,i n = Pr(X 1 = i 1 , . . . , X n = i n ) at the n leaves can be written as polynomials in the root probabilities and in the entries of the transition matrices. Denote by p the vector of joint probabilities p i 1 ,...,i n . As it is common in phylogenetic algebraic geometry, we will use the discrete Fourier transform for the groups G and G n to study the set of transition matrices and the set of joint probabilities at the leaves for a given phylogenetic tree and a group-based model. The reason for this is that phylogenetic invariants are considerably simpler in the Fourier coordinates (see Sturmfels and Sullivant 2005).
Denote byĜ the dual group of G whose elements are the group homomorphisms from G to the multiplicative group of complex numbers of magnitude one. Given a function a : G → C, its discrete Fourier transform is the functionǎ :Ĝ → C defined byǎ It is an invertible linear transformation given by the matrix K , where K g,h =ĝ(h).
The group-based model being symmetric is equivalent to the vectorsψ (e) andf (e) being real, see (Matsen 2009, Section 2). If we regard the vector p of joint probabilities as a function of G n , i.e., as an element of Hom (G n , C), then the image of p under the Fourier transform of G n is denoted q.
The map from the entries of the rate matrices to the joint probabilities at leaves can be seen as a composition of four maps: (1) • The map from {ψ (e) } e∈E to {ψ (e) } e∈E is given by the discrete Fourier transform of G. It is an invertible linear transformation given by the matrix K . • The map from {ψ (e) } e∈E to {f (e) } e∈E is given by by (Matsen 2009, Lemma 2.2). It is an isomorphism between R E×G and R E×G >0 . • In the case when root distribution puts all the mass at the identity, the map from {f (e) } e∈E to q is given by by (Székely et al. 1993, Theorem 3), where * g e = i∈ (e) g i and (e) is the set of observed leaves below e. See also (Sturmfels and Sullivant 2005, Sections 2 and 3) for a nice exposition of this result. In the case of the uniform root distribution, the identity (3) holds whenever g 1 + · · · + g n = 0. Otherwise q g = 0. This follows from (Sturmfels and Sullivant 2005, Lemma 4 and formula (12)). On the domain R E×G >0 , this map is injective: (Matsen 2009, Proposition 3.3 and Proposition 3.4) give a map from q to {[f (e) ] 2 } e∈E . Taking nonnegative square roots results in a left inverse to the map (3).
• The map from q to p is given by the inverse of the discrete Fourier transform of G n . It is an invertible linear transformation given by the matrix H −1 , where H is the n-fold Kronecker product of the matrix K .
Since π i , α e i , β e i are probabilities, they are real numbers in [0, 1], π 0 + π 1 = 1 and α e i + β e i = 1. Moreover, the restriction on the root distribution that it is uniform or defines a valid transition matrix in the CFN model gives 1 ≥ π 0 ≥ 1 2 and 1 2 ≥ π 1 ≥ 0; however, in Example 2 we will show that for the CFN model we can consider arbitrary root distributions. The determinant of P (e i ) is positive, because P (e i ) is the matrix exponential of a rate matrix Q (e i ) . Conversely, for every P (e i ) satisfying these constraints, there exists a rate matrix Q (e i ) such that P (e i ) = exp(t e i Q (e i ) ) by (Kingman 1962, Proposition 2).
The joint probabilities at the leaves have the parametrization In Sect. 3, we characterize this model in joint probabilities p i jk and without parameters π i , α e i , β e i . This is called the implicit description of a model. It consists of polynomial equations and inequalities in p i jk that describe the joint probabilities that come from a parametrization by rate matrices. In the Fourier coordinates, these equations can always be chosen to be binomials for any group-based model and tree (Evans and Speed 1993;Székely et al. 1993). These binomials are characterized in (Sturmfels and Sullivant 2005, Theorem 1). In the case of the CFN model on K 1,3 , these binomials are as was shown in (Sturmfels and Sullivant 2005, Example 3). The equations defining the model in the original coordinates can be obtained by applying the Fourier transformation of (Z 2 ) 3 on these binomials: Finally, we introduce basic notions from commutative algebra and algebraic geometry. A good introduction is given in Cox et al. (1992). Let R = R[x 1 , . . . , x n ] be a polynomial ring. A subset I ⊆ R is an ideal if it is an additive subgroup of R and is closed under multiplication by elements of the ring. The radical of an ideal I , denoted by √ I , consists of all the polynomials f ∈ R such that some power f m of f is in I . Let S be a set of polynomials in R and let k be a field. In this article, k is always R or C. The affine variety defined by S is Let f 1 , . . . , f s be the ideal generated by f 1 , . . . , f s , i.e., the smallest ideal containing f 1 , . . . , f s . Then . . , f s has maximal possible rank. Otherwise a point of the variety is called singular. Let T be a subset of k n . The Zariski closure T of T is the smallest affine variety containing T .

Implicit Descriptions of Symmetric Group-Based Models
Phylogenetic invariants are polynomials that vanish at joint probabilities at leaves for a given model and tree. They were introduced in Cavender and Felsenstein (1987) and Lake (1987) and have been characterized for group-based phylogenetic models in (Sturmfels and Sullivant 2005, Theorem 1). Phylogenetic varieties are algebraic varieties derived from phylogenetic models and were first introduced in Allman and Rhodes (2003,2004). In this paper, an algebraic variety is not necessarily irreducible. Phylogenetic invariants are elements of the ideal of a phylogenetic variety. Specifying a system of generators of the ideal of a phylogenetic variety is an important problem in phylogenetic algebraic geometry. However, the set of probability distributions forms only a (real, semialgebraic) subset of the phylogenetic variety, therefore providing a complete system of generators might have no biological interest. In Casanellas et al. (2015), a minimal set of phylogenetic invariants is constructed that defines the intersection of a phylogenetic variety with a Zariski open set. In the case of the Kimura 3-parameter model, all the leaf probabilities that are images of real parameters in the phylogenetic model (not in the complexification of the model) lie in this Zariski open set. The number of polynomials in this set is equal to the codimension of the phylogenetic variety and each polynomial has degree at most |G|. This reduces drastically the number of phylogenetic invariants used: For the Kimura 3-parameter model on a quartet tree, it drops from 8002 generators of the ideal to the 48 polynomials described in (Casanellas and Fernández-Sánchez 2008, Example 4.9).
Besides phylogenetic invariants, polynomial inequalities are needed to give an exact characterization of joint probabilities at leaves for a given model and a tree. For general symmetric group-based models, polynomial inequalities that describe joint probabilities at leaves are studied in Matsen (2009). We recall (Matsen 2009, Propositions 3.3 and 3.4) that give the left inverse to the map (3) on the domain R E×G >0 .
Proposition 1 (Matsen (2009), Proposition 3.3) Given some leaf edge e, let i denote the leaf vertex incident to e and let v be the internal vertex incident to e. Let j, k be leaf vertices different from i such that the path from j to k contains v. Let w(g i , g j , g k ) ∈ G n assign state g x to leaf x for x ∈ {i, j, k} and zero to all other leaf vertices. Then Proposition 2 (Matsen (2009) The next proposition will summarize the procedure in Matsen (2009) to construct inequalities that describe joint probabilities. We will denote by (K −1 ) g,: the row of the matrix K −1 labeled by g and by (f (e) ) (K −1 ) g,: the Laurent monomial Proposition 3 Assume that the labeling function L satisfies L(g) = L(−g) for all g ∈ G. Consider the set of {ψ (e) } e∈E that satisfies g∈G ψ (e) (g) = 0, ψ (e) (g 1 ) = ψ (e) (g 2 ) whenever L(g 1 ) = L(g 2 ) and ψ (e) (g) ≥ 0 for all nonzero g ∈ G. The images of this set under the maps in (1) are: This equation and inequalities are equivalent tof (e) Here we have squared the inequalities (f (e) ) (K −1 ) g,: ≥ 1. (iii) The constraints for q are given by phylogenetic invariants, equation q 00...0 = 1, inequalities q > 0 and inequalities that are obtained by substituting expressions for [f (e) ] 2 in Propositions 1 and 2 to inequalities (f (e) ) 2(K −1 ) g,: ≥ 1 in the previous item. (iv) The constraints for p are obtained by substituting q by H p in the constraints for q.
For the sake of completeness, a proof of Proposition 3 is given in "Appendix A".

Remark 1
In Proposition 3 item (iii), one applies Propositions 1 and 2 to obtain inequalities in the Fourier coordinates. However, in Propositions 1 and 2 one has a choice in choosing the leaf vertices. Since the Fourier coordinates are strictly positive, then any choice of leaf vertices in Propositions 1 and 2 gives equivalent inequalities in Proposition 3 item (iii) and it does not matter which choice we make.
The implicit description of the CFN model on the tree K 1,3 for an arbitrary root distribution is given as the union of three sets: the set defined by equations and inequalities (5)

-(19), the set defined by equations and inequalities (20)-(34) and the set defined by equations and inequalities (35)-(48).
Remark 2 Identifiability of parameters of a phylogenetic model means that if for a fixed tree two sets of parameters map to the same joint probabilities at leaves, then these sets of parameters must be equal. Generic identifiability means that this statement is true with probability one. The identifiability of the CFN model was shown in (Hendy 1991, Theorem 1), of the Kimura 3-parameter model in (Steel et al. 1998, Theorem 7) and the generic identifiability of the general Markov model in Chang (1996). The identifiability of any group-based model follows also from the proof of Proposition 3, since each of the maps in (1) is an isomorphism in the region we are interested in.

Corollary 1 Consider a symmetric group-based model. Any p satisfying the equations and inequalities described in Proposition 3 that satisfies one of the inequalities with equality comes from a parametrization with an off-diagonal zero in the rate matrix Q (e) for some e ∈ E.
Proof There are two different kinds of inequalities in item (4) of Proposition 3. The strict inequalities can never be satisfied with equality. The non-strict inequalities in each step are obtained by substituting the inverse map to the inequalities in the previous step. Hence p satisfies one of the non-strict inequalities with equality if and only if it has a preimage {ψ (e) } e∈E that satisfies one of the inequalities ψ (e) (g) ≥ 0 with equality.

Example 3
We consider the CFN model. A joint probability vector p satisfying the assumptions of Corollary 1 has in its parametrization the rate matrix Q (e) = 0 0 0 0 for some e ∈ E. The transition matrix corresponding to the same edge is P (e) = 1 0 0 1 .

Maximum Likelihood Estimation via Numerical Algebraic Geometry
In this section, we use the terminology and notation introduced in Sect. 2. In particular, p i 1 ,...,i n are the joint probability distributions at the n-leaves. Let u = (u i 1 ,...,i n ) (i 1 ,...,i n )∈G n be a vector of observations at leaves. The log-likelihood function of a phylogenetic model is Maximum likelihood estimation aims to find a vector of joint probability distributions at leaves or model parameters (if the joint probabilities are considered as polynomials in model parameters) that lies in the model and maximizes the log-likelihood function for a given observation u.

Example 4
In (Hosten et al. 2005, Example 14), maximum likelihood estimation on the Zariski closure of the CFN model on K 1,3 is considered. This is the model that is defined by the equations in Example 2. For generic data, the number of complex critical points of the likelihood function on the Zariski closure of a model is called the ML degree. It is shown in (Hosten et al. 2005, Example 14) that the ML degree of the CFN model on K 1,3 is 92. Using tools from numerical algebraic geometry, one can compute the 92 critical points and among the real critical points choose the one that gives the maximal value of the log-likelihood function. However, the MLE can lie on the boundary of a statistical model or even not exist. Neither of this can be detected by considering only the Zariski closure of the model. We will see the latter happening for the CFN model on K 1,3 in Example 5.
In practice, the MLE is solved using numerical methods such as the Newton-Raphson method (Schadt et al. 1998;Kenney and Gu 2012), quasi-Newton methods Olsen et al. (1994) and the EM algorithm (Felsenstein 1981;Friedman et al. 2002;Holmes and Rubin 2002;Hobolth and Jensen 2005). However, since these methods are hill-climbing methods and the likelihood function on phylogenetic trees can have multiple local maxima (Steel 1994;Chor et al. 2000), they are only guaranteed to give a local maximum or a saddle point of the log-likelihood function and not necessarily the global maximum. Usually one uses a heuristic to find a good initialization for these methods or runs them for different starting points and chooses the output that maximizes the log-likelihood function.
We suggest a global method based on numerical algebraic geometry that theoretically gives the solution to the maximum likelihood estimation problem on phylogenetic trees with probability one. The main idea behind numerical algebraic geometry is homotopy continuation. Homotopy continuation finds isolated complex solutions of a system of polynomial equations starting from the known solutions of another system of polynomial equations. Numerical algebraic geometry methods give theoretically correct results with probability one, meaning that bad phenomena can happen when certain parameters are chosen from a measure zero set. An introduction to numerical algebraic geometry can be found in Sommese and Wampler (2005), Bates et al. (2013). In our context, the system of polynomial equations that we wish to solve comes from the Karush-Kuhn-Tucker (KKT) conditions (Karush 1939;Kuhn and Tucker 1951) for the optimization problem that maximizes the likelihood function on a phylogenetic model. The set of solutions of this polynomial system contains all the critical points of the likelihood function. The global maximum of the likelihood function is the solution of the polynomial system that maximizes the likelihood function among all the solutions that lie in the model. This global approach for solving a nonconvex optimization problem on a set that is described by polynomial equations and inequalities has been previously employed in optimal control Rostalski et al. (2011) and in the life sciences (Gross et al. 2016). Our setup and algorithm are similar to those in Rostalski et al. (2011), although we provide further lemmas that allow us to decompose the system of polynomial equations that we want to solve to simpler systems of polynomial equations. The article Gross et al. (2016) uses Fritz John conditions instead of KKT conditions and focuses mostly on optimization problems on sets that are described by polynomial equations only. Sets that are described by polynomial equations and inequalities are considered in Section 3 of the supplementary material of Gross et al. (2016). In particular, the ideas for Theorem 1 and Remark 3 appear there.
More specifically, consider the optimization problem If x * is a local optimum and the optimization problem satisfies first-order constraint qualifications, then there exist μ i , where i = 1, . . . , m, and λ j , where j = 1, . . . , l, such that x * satisfies the KKT conditions: One first-order constraint qualification is the constant rank constraint qualification (CRCQ) defined in Janin (1984). A point satisfies the CRCQ if there is a neighborhood of the point where gradients of the equations and gradients of the active inequalities, i.e., inequalities that the point satisfies with equality, have constant rank. We also consider the optimization problem If x * is a local optimum of the optimization problem (55), then there exist λ j , where j = 1, . . . , l, such that x * satisfies the Lagrange conditions: In the rest of the section, we assume that the KKT conditions (50)-(54) and the Lagrange conditions (56)-(57) are polynomial. In this case, a point satisfies the CRCQ if it is a smooth point of the variety defined by the equations and active inequalities.
Let L ⊆ C[μ, λ, x] be the ideal generated by the polynomials on the left-hand sides of the equations (50), (52) and (54) in the KKT conditions. For S ⊆ [m], let L S ⊆ C[μ S , λ, x] be the ideal generated by the polynomials in the Lagrange conditions for the optimization problem Specifically, let L S ⊆ C[μ S , λ, x] be generated by the polynomials We denote by I S ⊆ C[x] the ideal generated by the constraints in the above optimization problem, i.e., I S = G i , H j : i ∈ S, j = 1, . . . , l .

Theorem 1 Let L and L S be as defined above. Then
The idea behind Theorem 1 is that instead of optimizing a function over a semialgebraic set, one can optimize the function over the Zariski closure of the semialgebraic set and the Zariski closures of each of the boundaries of the semialgebraic set. This concept is discussed in Section 3 of the supplementary material of Gross et al. (2016).
We have shown that π x (V (L)) = ∪π x (V (L S )), where π x is the projection of (μ, λ, x) or (μ S , λ, x) on x. By the Closure Theorem (Cox et al. 1992, Theorem 3.2 ) holds, because the right-hand side is a variety and contains ∪π x (V (L S )) and hence π x (V (L)). On the other hand, since π Theorem 1 suggests Algorithm 1 for solving the equations in the KKT conditions. Algorithm 1 is related to (Rostalski et al. 2011, Algorithm 1) and (Gross et al. 2016, Algorithm 3).

Algorithm 1 Global maximum of a polynomial optimization problem
Step 1: Let C = {}.
Step 2: For every S ⊆ [m], if dim(V (L S )) = 0, then add all elements of V (L S ) to C.
Step 3: Remove the elements of C that are not real or do not satisfy G i (x) ≥ 0 or μ i ≥ 0 for i = 1, . . . , m.
Step 4: Find the element (μ * S , λ * , x * ) of C that maximizes F. Output: The element x * from Step 4.  (52) and (54) in the KKT conditions. Since the global maxima satisfy the CRCQ, they must be solutions of these equations. By choosing among the real solutions that satisfy inequalities (51) and (53) in the KKT conditions the ones that maximize the value of the cost function F, we get the global maxima.
We are interested in the optimization problem, when the cost function is the loglikelihood function l u (p) = u i 1 ,...,i n log p i 1 ,...,i n and the constraints are polynomial equations and inequalities that describe a statistical model (written as H j (p) = 0 for j = 1, . . . , l and G i (p) ≥ 0 for i = 1, . . . , m, respectively). Although Eq. (50) is not polynomial for F = l u , it can be made polynomial by multiplying the equation One of the reasons why the variety V (L S ) in Step 2 of Algorithm 1 might not be finite is that the Lagrange conditions for MLE might be satisfied by higher-dimensional components where some variable is identically zero. For MLE, Gross and Rodriguez have defined a modification of the Lagrange conditions, known as Lagrange likelihood equations (Gross and Rodriguez 2014, Definition 2), whose solution set does not contain solutions with some variable equal to zero if the original data does not contain zeros (Gross and Rodriguez 2014, Proposition 1). However, the Lagrange likelihood equations can be applied only to homogeneous prime ideals. This motivates us to study Lagrange conditions for decompositions of ideals.
k (x) = 0 for j = 1, . . . , m 1 , k = 1, . . . , m 2 . The Lagrange conditions for the latter optimization problem are k , we see that x * satisfies Lagrange conditions for the optimization problem max F(x) subject to G The Lagrange conditions for the generators of J + K are If there exists k 1 such that G k , we see that x * satisfies Lagrange conditions for the optimization problem max F(x) subject to the generators of J 1 + K 1 . If G (2) k (x * ) = 0 for all k and/or H (2) k (x * ) = 0 for all k, then we get other combinations J 1 + K 2 , J 2 + K 1 or J 2 + K 2 .
Lemma 1 suggests that if S is a singleton in Step 2 of Algorithm 1, then we can replace the ideal L S of Lagrange conditions for I S in Step 2 of Algorithm 1 by the ideals of Lagrange conditions for minimal primes of I S . If S = {i 1 , . . . , i |S| }, then I S = I {i 1 } + . . . + I {i |S| } . Hence by Lemmas 1 and 2, we can replace the ideal L S by the ideals of Lagrange conditions for the sum of minimal primes of I {i j } , where 1 ≤ j ≤ |S|.
Remark 3 As discussed in Section 3.2 of the supplementary material to Gross et al. (2016), one can ignore all the components where one of the constraints is x k = 0 or the sum of some variables is zero. If one of the variables is zero, then the value of the log-likelihood function is −∞. If the sum of some variables is zero, then all of them have to be zero, because none of them can be negative.
We summarize the results of Lemmas 1, 2 and Remark 3 in Algorithm 2. The output of Algorithm 2 is a list of ideals. For each of the ideals consider the optimization problem where equation constraints are given by the generators of the ideal. The ideals generated by the Lagrange conditions for the optimization problems can be used in Step 2 of Algorithm 1 instead of the ideals L S for every S ⊆ [m].

Algorithm 2 A list of ideals for
Step 2 of Algorithm 1 Input: An optimization problem Step 1: Let P be the power set of [m].
Step 2: For each S ∈ P associate a list of ideals: • If S = {}, then the list of ideals associated to S consists of these minimal primes of H j (x) : j = 1, . . . , l that do not contain any sums of the variables. • If S = {s} for some 1 ≤ s ≤ m, then the list of ideals associated to S consists of these minimal primes of G s (x), H j (x) : j = 1, . . . , l that do not contain any sums of the variables. • If |S| > 1, write S 2 = {max(S)} and S 1 = S\S 2 . The list of ideals associated to S consists of these minimal primes of the pairwise sums of the ideals in the list associated to S 1 and in the list associated to S 2 that do not contain any sums of the variables.
Step 3: Take the union of all lists in Step 2 and remove repeated ideals. Output: The list of ideals from Step 3.

Remark 4
In practice, it is crucial to know the degrees of the ideals L S of Lagrange conditions. We recall that these degrees are also known as ML degrees. Although in theory, polynomial homotopy continuation finds all solutions of a system of polynomial equations with probability one, in practice, this can depend on the settings of the program. Without knowing the ML degree, there is no guarantee that any numerical method finds all critical points. For the CFN model on K 1,3 , we experimented with Bertini Bates et al. (2006), NumericalAlgebraicGeometry package in Macaulay2 Leykin (2011) and PHCpack Verschelde (1999). We ran these three programs with default settings to find the critical points of the log-likelihood function on the Zariski closure of the CFN model on K 1,3 . For our example, only PHCpack found all 92 critical points discussed in Example 4.

Example 5
We aim to compute the MLE for the CFN model on K 1,3 and the data vector u = (17,5,27,5,16,5,19,6). This data vector is obtained by generating 100 samples from the distribution inside the CFN model with rate parameters The corresponding tree is depicted in Figure 1. It has two short edges, one long edge and the root distribution is very close to the uniform distribution.
To find the MLE, we have to consider three different optimization problems corresponding to the three different cases in Example 2. In each of the cases, we relax the implicit characterization given in Example 2 by replacing strict inequalities with non-strict inequalities. Specifically, in the first case, the polynomials G i are given by the left-hand sides of the inequalities (9)-(19) and the polynomials H j are given by the left-hand sides of Eqs. (5)-(8); in the second case, the polynomials G i are given by Fig. 1 The tree in Example 5 has two edges with short branch lengths 1 and 3 , one edge with a long branch length M and the root distribution is very close to the uniform distribution the left-hand sides of the inequalities (24) (42) and (45). We apply the modified version of Algorithm 1 that uses the output of Algorithm 2 in Step 2. It is enough to run Algorithm 2 and Step 2 of Algorithm 1 for the first optimization problem only as the polynomials G i and H j are the same for the first two optimization problems; in the third optimization problem there is one polynomial less and some polynomials G i are among polynomials H j , but all ideals considered in Algorithm 2 and Step 2 of Algorithm 1 for the third optimization problem are among the ideals for the first optimization problem. In Step 3 we have to check whether elements satisfy G i (x) ≥ 0 and H j (x) = 0 for any of the three optimization problems. The code for this example can be found at the link: https://github.com/kaiekubjas/phylogenetics As a result, we obtain 44 ideals summarized in Table 1. The first row of this table corresponds to the Zariski closure of the CFN model on K 1,3 . It has degree 92 which agrees with the ML degree 92 computed in (Hosten et al. 2005, Example 14). However, to find the MLE one has to consider critical points of the likelihood function in the interior and on all the boundary components, in total 167 of them. We compute all the 167 complex critical points using numerical algebraic geometry software PHCpack. Out of the 167 complex critical points 99 are real and 51 are positive. We list the seven points among them that have the highest value of the log-likelihood function in Table 2. The first six critical points in Table 2 satisfy p 000 − p 001 + p 010 − p 011 + p 100 − p 101 + p 110 − p 111 > 0 and p 000 + p 001 − p 010 − p 011 + p 100 + p 101 − p 110 − p 111 < 0.
Hence these critical points are not in the CFN model on K 1,3 as in all three cases in Example 2, the two linear inequalities are satisfied with the same sign. The critical point with the seventh highest log-likelihood value is in the image of the following parameters: This implies that the MLE for the CFN model on K 1,3 and the data vector u = (17,5,27,5,16,5,19,6) does not exist-the global maximum of the log-likelihood function is achieved when we allow one of the parameters to go to infinity. Strictly speaking this statement is true for the set of points in the model that satisfy CRCQ. We believe that for random data the global maximum will satisfy CRCQ with probability one. When we run the same optimization problem in Mathematica, then we get a solution with similar value for the log-likelihood function and all parameters besides ψ (e 2 ) , which is equal to ψ (e 2 ) = (−8.120, 8.120). Without having the implicit description of the CFN model on K 1,3 and using numerical algebraic geometry to study the MLE, it would be very difficult to say that the MLE does not exist.

Remark 5
In Example 5, we chose the rate parameters of the true data generating distribution such that the joint leaf probabilities of this distribution would be close to the boundary of the model. In particular, the Fourier leaf probabilities q 010 , q 011 , q 110 , q 111 are almost zero. We recall that the semialgebraic description of the CFN model includes strict inequalities q > 0. The global maximum of the likelihood function on the closure of the CFN model on K 1,3 satisfies q 010 = q 011 = q 110 = q 111 = 0. Since this global maximum is not in the model, the MLE does not exist. We expect the similar phenomenon that if our true data generating distribution is close to the boundary, then the MLE does not exist to happen with nonzero probability. In particular, if the normalized data vector lies on the part of the boundary that is not in the model, then we know that the MLE does not exist. the same inequalities in item (ii) and it clearly maps to the same q. Indeed, since the inequalities are of the form (f (e) ) 2(K −1 ) g,: ≥ 1, it means that in the product (f (e) ) 2(K −1 ) g,: minus signs cancel out and hence the absolute values give the same product.
The map (3)  to { q ag q bg } g∈G n for some vectors a g , b g ∈ R G n . Since q g = q ag q bg , or equivalently q 2 g q b g = q a g , for all q in the image, the same equation must be satisfied for all elements in the Zariski closure of the image. Moreover, q ag q bg is well-defined on the positive part of the Zariski closure, hence we have the isomorphism. It follows that on the positive part of the Zariski closure we get the inequalities for q by substituting expressions for [f (e) ] 2 to inequalities (f (e) ) 2(K −1 ) g,: ≥ 1.
This completes the proof that the equations and inequalities in items (i)-(iv) are correct.