{'fraction_non_alphanumeric': 0.05000142130248159, 'fraction_numerical': 0.022968248102561185, 'mean_word_length': 3.973141080011309, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Geometry of symmetric group-based models

Phylogenetic models have polynomial parametrization maps. For symmetric group-based models, Matsen studied the polynomial inequalities that characterize the joint probabilities in the image of these parametrizations. 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. We also study the embedding problem for symmetric group-based models, i.e. we identify which mutation matrices are matrix exponentials of rate matrices that are invariant under a group action.

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 probability and the rate or mutation 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.
In 1987, Cavender and Felsenstein [3] and, separately, Lake [17], 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 cut out all the biologically meaningful points. 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 group-based models which require the mutation matrices to be invariant under the action of an abelian group G. A symmetric group-based model assumes that the mutation matrices are symmetric. In particular, a symmetric group-based model can be a submodel of a non-symmetric group-based model with extra symmetricity conditions on mutation matrices. Phylogenetic invariants for group-based models are described in [23]. A smaller set of phylogenetic invariants that cut out biologically meaningful points is given in [5]. Our first aim is to understand polynomial inequalities that characterize leaf probabilities for symmetric group-based models. A method for deriving the inequalities in the Fourier coordinates is given in Matsen [19,Proposition 3.5]. For the sake of completeness, we present a proof of [19,Proposition 3.5] here, and we add certain positivity constraints that appear for the Kimura 3-parameter model in [4]. We explicitly derive the semialgebraic description of the leaf probabilities of the CFN model on the tripod tree K 1,3 .
Identifying the equality 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. In practice, numerical methods are used to find the MLE. However, as MLE is a non-convex optimization problem, there is no guarantee of outputting the global optimum. Since phylogenetic models are not necessarily compact, the MLE might even not exist. We will use methods from computational and numerical algebraic geometry similar to [11] 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 have been very difficult to discover without the algebraic methods that we use in this paper. One should see the example for the CFN model on the tripod K 1,3 as an illustration how to use numerical algebraic geometry for MLE in phylogenetics. 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.
Finally, we consider the embedding problem for symmetric group-based models. Phylogenetic models are originally defined using rate matrices, e.g. matrices whose off-diagonal entries are nonnegative and the entries in each row sum to one. Mutation matrices are matrix exponentials of rate matrices. We furthermore assume that both the rate and mutation matrices satisfy the assumptions of a symmetric group-based model, i.e. they are symmetric and invariant under the action of an abelian group G. We call all mutation matrices that are invariant under the action of G and that can be written as matrix exponentials of rate matrices that are invariant under the action of G, G-embeddable. We will characterize G-embeddable mutation matrices for symmetric group-based models. For the special cases of the CFN and the Kimura 3-parameter model, such characterizations are given in [16,21]. Our characterization is for any number of states k (for example k = 20 corresponds to amino acids).
In Section 2, we introduce the preliminaries of phylogenetic models and present tools from [19]. Proposition 3 in Section 3 gives an algorithm for deriving the semialgebraic description of the leaf probabilities of a symmetric group-based model. This is mostly application of [19], however, it also considers additional positivity constraints. The main result in Section 4 is Theorem 1 that gives a characterization of G-embeddable mutation matrices for symmetric group-based models. Finally, Algorithm 1 in Section 5 outlines how to use numerical algebraic geometry to compute the MLE with probability one. This algorithm is applied on the CFN model on the tripod in Example 8.

Preliminaries
The exposition in this section largely follows [19]. A phylogenetic tree T is a rooted tree with n labeled leaves and it represents the genetic 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 that the random variables at vertices evolve according to a Markov process, i.e. there is a transition (mutation) matrix P (e) associated to every edge e that reflects 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 non-negative 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) .
Without loss of generality, we will assume that the distribution at the root is uniform. If the distribution at the root is arbitrary, then the following procedure can be applied to reduce to the uniform case: One adds a new edge starting from the root and reroots the tree at the additional leaf. The previous root distribution is obtained by using a suitable transition matrix for the new edge. The only difference is that the new leaf is hidden while other leaves are observed.
In this paper, we investigate group-based phylogenetic models. This means that we fix an abelian group G, set k = |G| and assume that rate matrices are invariant under the action of G, i.e. there exists a vector ψ (e) ∈ R G such that Q for a probability vector f (e) ∈ R G . The phylogenetic models we consider are symmetric, which means Q h,g . In the case of group-based models, this is equivalent to ψ (e) (g) = ψ (e) (−g) and implies f (e) (g) = f (e) (−g).
The joint probabilities p i1,...,in = 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 mutation matrices. Denote by p the vector of joint probabilities p i1,...,in . As it is common in phylogenetic algebraic geometry, we will use the discrete Fourier transform to study the set of mutation 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.
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 functioň a :Ĝ → C defined byǎ (ĝ) = h∈Gĝ (h)a(h).
It is an invertible linear transformation given by the matrix K, where K g,h = g(h). The image of the vector p of joint probabilities 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.
(2) The map from {ψ (e) } e∈E to {f (e) } e∈E is given by by [19,Lemma 2.2]. It is an isomorphism between R E×G and R E×G >0 . (3) The map from {f (e) } e∈E to q is given by by [24,Theorem 3], where * g e = i∈Λ(e) g i and Λ(e) is the set of observed leaves below e. On the domain R E×G >0 , this map is injective: [19, 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). (4) 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.
Example 1 We will consider in detail the Cavender-Farris-Neyman (CFN) model [2,10,20] on the rooted claw tree T = K 1,3 . This example has been previously studied in [23,Example 3] and [14,Example 14]. The CFN model is a group-based model with G = Z 2 and k = 2. Denote the root distribution by π = (π 0 , π 1 ) and the mutation matrices at edges e 1 , e 2 , e 3 by Since π i , α ei , β ei are probabilities, they are real numbers in [0, 1] and α ei + β ei = 1. Moreover, the determinant of P (ei) is positive, because P (ei) comes from a rate matrix Q (ei) . Conversely, for every P (ei) satisfying these constraints, there exists a rate matrix Q (ei) such that P (ei) = exp(t ei Q (ei) ) by [16,Proposition 2]. In Section 4, we study constraints on mutation matrices for general symmetric group-based models.
In Section 3, we characterize this model in joint probabilities p ijk and without parameters π i , α ei , β ei . This is called the implicit description of a model. It consists of polynomial equations and inequalities in p ijk that cut out 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 [9,24]. These binomials are characterized in [23,Theorem 1]. In the case of the CFN model on K 1,3 , these binomials are {q 001 q 110 − q 000 q 111 , q 010 q 101 − q 000 q 111 , q 100 q 011 − q 000 q 111 }.
The equations defining the model in the original coordinates can be obtained by applying the Fourier transformation of (Z 2 ) 3 on these binomials: q 000 = p 000 + p 001 + p 010 + p 011 + p 100 + p 101 + p 110 + p 111 , q 001 = p 000 − p 001 + p 010 − p 011 + p 100 − p 101 + p 110 − p 111 , Finally, we introduce basic notions from commutative algebra and algebraic geometry. A good introduction is given in [7]. Let R = R[x 1 , . . . , x n ] be a polynomial ring. A subset I ⊆ R is an ideal, if it is closed under addition and multiplication by scalars. 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 V (S) = {(a 1 , . . . , a n ) ∈ k n : f (a 1 , . . . , a n ) = 0 for all f ∈ S}.
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 description
Phylogenetic invariants are polynomials that vanish at joint probabilities at leaves for a given model and tree. They were introduced in [3,17] and have been characterized for group-based phylogenetic models in [23,Theorem 1]. A subset of them containing the biologically meaningful points is given in [5]. This subset forms a local complete intersection and each polynomial in this subset has degree at most |G|. It 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 [4,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 the Kimura 3-parameter model, a set of inequalities is given in [4]. For general symmetric group-based models, polynomial inequalities that cut out joint probabilities at leaves are investigated in [19]. This section is essentially an exposition of results in [19]. There are three reasons why we present it here: For the sake of completeness, we will give a proof of [19, Proposition 3.5], we will include the condition that Fourier coordinates need to be positive, and we will apply part of this exposition to the study of the embedding problem in Section 4.
We recall [19, Propositions 3.3 and 3.4] that give the left inverse to the map 3 on the domain R E×G >0 .

Proposition 1 ([19], 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 ([19], Proposition 3.4)
Given some internal edge e, let the two vertices incident to e be v and v ′ . Let i, j (respectively i ′ , j ′ ) be leaf vertices such that the path from i to j (respectively the path from Proposition 3 Consider the set of {ψ (e) } e∈E that satisfies g∈G ψ (e) (g) = 0 and ψ (e) (g) ≥ 0 for all non-zero g ∈ G. The images of this set under the maps in (1) are: (1) The constraints for {ψ (e) } e∈E are obtained by substituting ψ (e) by K −1ψ(e) in the constraints for {ψ (e) }. In particular, this givesψ (e) (0) = 0 and (K −1ψ(e) )(g) ≥ 0 for all non-zero g ∈ G.
(3) The constraints for q are given by phylogenetic invariants, equality 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 multiplicative inequalities in the previous item. (4) The constraints for p are obtained by substituting q by Hp in the constraints for q.
Here (f (e) ) K −1 g,: denotes the Laurent monomial h∈G (f Proof The constraints in items (1) and (4) are obtained, because the corresponding maps are invertible linear transformations. We will prove that the constraints in items (2) and (3) are correct.
Proof The positivity constraints come from the exponentiation map. Moreover, a T x ≥ 0 is equivalent to e a T x ≥ 1, and since e a T x = (e x ) a T , it is also equivalent to (e x ) a T ≥ 1. Hence the equalities and inequalities for {f (e) } e∈E are the multiplicative versions of the equalities and inequalities for {ψ (e) } e∈E together withf (e) (g) > 0 for all g ∈ G. (3) is cut out by the constraints in item (3).
Lemma 2 is very similar to [19,Proposition 3.5], however, for the sake of completeness, we will give a proof here. We also include the positivity constraints that do not appear in [19,Proposition 3.5].
Proof The inequalities q > 0 are clearly valid inequalities. We will show that we do not have to additionally consider the inequalitiesf (e) > 0 to construct inequalities for q. Assume there is {f (e) } e∈E with not all entries positive that satisfies all other inequalities in item (2) and maps to q > 0. We claim that {|f (e) |} e∈E also satisfies the same inequalities in item (2) and it clearly maps to the same q. Indeed, since the inequalities are of the form (f (e) ) 2K −1 g,: ≥ 1, it means that in the product (f (e) ) 2K −1 g,: minus signs cancel out and hence the absolute values give the same product.
The map (3) (3), map q 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 bg = q ag , 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 multiplicative inequalities forf (e) .
This completes the proof that inequalities in items (1)-(4) are correct.
Example 2 We will derive the implicit description of the CFN model on the rooted claw tree T = K 1,3 . In addition to phylogenetic invariants in Example 1, applying Proposition 3 gives the following inequalities in Fourier coordinates: The inequality q100q010 q110 ≤ 1 is for the hidden leaf corresponding to the root. Since q 000 = 1, we can multiply all the denominators by q 000 without changing the inequalities (4). Clearing denominators gives the following polynomial inequalities q 000 = 1, q > 0, q 000 q 110 − q 100 q 010 ≥ 0, q 000 q 011 − q 110 q 101 ≥ 0, q 000 q 001 − q 011 q 010 ≥ 0, q 000 q 100 − q 001 q 101 ≥ 0.
By applying the discrete Fourier transformation, we get the implicit description of the CFN model on K 1,3 in the original coordinates p 000 − p 001 − p 010 + p 011 + p 100 − p 101 − p 110 + p 111 > 0, p 000 + p 001 + p 010 p 000 + p 001 − p 010 − p 011 − p 100 − p 101 + p 110 + p 111 > 0, p 000 − p 001 − p 010 + p 011 − p 100 + p 101 + p 110 − p 111 > 0, Remark 1 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 [13, Theorem 1], of the Kimura 3-parameter model in [22,Theorem 7] and the generic identifiability of the general Markov model in [6]. The identifiability of any group-based model follows also from the proof of Proposition 3, since each of the maps in (1) is isomorphism in the region we are interested in.
Corollary 1 Consider a symmetric group-based model. Any p satisfying the implicit constraints of the model obtained by 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 kind 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 .

Embedding problem
Deciding whether a stochastic matrix P is the matrix exponential of a rate matrix Q, is called the embedding problem. A stochastic matrix that can be represented as the matrix exponential of a rate matrix, is said to be embeddable. An overview on embeddable stochastic matrices is given in [8]. We call a stochastic matrix that is invariant under the action of G G-embeddable if it can be written as a matrix exponential of a rate matrix that are invariant under the action of G. The aim of this section is to give an exact characterization of G-embeddable stochastic matrices for any symmetric group-based model.
Example 4 Without loss of generality, one may assume that a rate matrix in the CFN model takes the form Applying P (e) = exp(t e Q (e) ) for t e ≥ 0 gives Studying this matrix, one obtains that P (e) is the matrix exponential of a rate matrix in the CFN model if and only if f (e) (0) + f (e) (1) = 1, 1 ≥ f (e) (0) > 1 2 and 1 2 > f (e) (1) ≥ 0. These conditions are equivalent to P (e) being a stochastic matrix invariant under the Z 2 -action and satisfying det(P (e) ) > 0, or equivalently tr(P (e) ) > 1. This result is stated for 2 × 2 stochastic matrices (not necessarily invariant under the Z 2 -action) in [16, Proposition 2].
For n × n stochastic matrices, a necessary condition for being embeddable is given in [16,Proposition 3]: The set set of embeddable matrices is relatively closed as a subset of the space of all stochastic n × n matrices with positive determinant. An exact characterization for stochastic matrices embeddable in the Kimura 3-parameter model (the group-based model with G = Z 2 × Z 2 ), is presented in [21,Theorem 3.2]. We present an exact characterization for stochastic matrices embeddable in any symmetric group-based model (for any number of states k).
Theorem 1 Fix a group G and let K be the matrix of the Fourier transform of G. Let P (e) g,h = f (e) (h − g) for some stochastic vector f (e) ∈ R G satisfying f (e) (g) = f (e) (−g) for every g ∈ G. Then P (e) = exp(Q (e) ) where Q (e) g,h = ψ (e) (h − g) is a rate matrix for ψ (e) satisfying ψ (e) (g) = ψ (e) (−g) if and only if (1) Re(K g,· ), f (e) > 0 for every non-zero g ∈ G; (2) for every non-zero g ∈ G.
Here Re(·) denotes the real part of a number or a vector. Proof The conditions on f (e) follow from the conditions onf (e) by substitutinǧ f (e) by Kf (e) . By Proposition 3 item (2), the conditions forf (e) aref (e) (0) = 1, (f (e) ) K −1 g,: ≥ 1 for all non-zero g ∈ G andf (e) (g) > 0 for all g ∈ G. The equalityf (e) (0) = 1 gives g∈G f (e) (g) = 1. The inequalities (1) follow fromf (e) (g) > 0 for all non-zero g ∈ G. The equation K g,−h = K g,h always holds for a discrete Fourier transform and f (e) (h) = f (e) (−h) holds because of symmetry. Hence The inequalities (2) follow from (f (e) ) K −1 g,: ≥ 1 for all non-zero g ∈ G. We can consider the real part for exponents, since
All factors in this product are positive by inequalities (1) in Theorem 1 and thus det(P (e) ) > 0. More precisely, the set of G-embeddable matrices for a symmetric group-based model is a relatively closed subset of one connected component of the complement of det(P (e) ) = 0 which is given by inequalities (1).

Example 5
The discrete Fourier transformation for the CFN model is given by the matrix By Theorem 1, the conditions on mutation matrices P (e) in the CFN model are An easy check verifies that these conditions are equivalent to the conditions in Example 4.
By Theorem 1, the conditions on the mutation matrices to be Z 2 × Z 2embeddable are This characterization for the Kimura 3-parameter model is first given in [21,Theorem 3.2]. Moreover, it is shown in [21] that matrices satisfying (16) constitute 3 32 of all the stochastic matrices invariant under the Z 2 × Z 2 action and 1 2 of all such stochastic matrices with positive eigenvalues.
Remark 3 By [16, Corollary on page 18], the map from rate matrices to mutation matrices is locally homeomorphic except possibly when the rate matrix has a pair of eigenvalues differing by a non-zero multiple of 2πi. Since for symmetric group-based models rate matrices are real symmetric, then all their eigenvalues are real and hence the map from rate matrices to mutation matrices is a homeomorphism. This can be also seen by analyzing maps in Maximum likelihood estimation aims to find a vector of joint probabilities at leaves or model parameters (if the joint probabilities are considered as polynomials in model parameters) that lie in the model and maximize the log-likelihood function for a given observation u.

Example 7
In [14,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. It is shown that for generic data, the likelihood function has 92 complex critical points on the model. This is called the ML degree of the model. Using tools from numerical algebraic geometry as in [11], 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 8.
In practice, the MLE is solved using numerical methods, however, these methods are only guaranteed to give a local maxima of the log-likelihood function and not necessarily the global maximum. Usually one runs these methods for different starting points and chooses the output that maximizes the log-likelihood function.
An alternative is to solve the Karush-Kuhn-Tucker (KKT) conditions to find all the critical points of the log-likelihood function and then choose the real solution that maximizes the log-likelihood function. Consider the optimization problem max F (x) subject to (17) G i (x) ≥ 0 for i = 1, . . . , m, H j (x) = 0 for j = 1, . . . , l.
The KKT conditions are H j (x) = 0 for j = 1, . . . , l, 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 above. One first-order constraint qualification is the constant rank constraint qualification (CRCQ) defined in [15]. A point satisfies the CRCQ if there is a neighborhood of the point where gradients of the equality constraints and gradients of the active inequality constraints have constant rank.
In the rest of the section, we assume that conditions (18)- (22) are polynomial. In this case, a point satisfies the CRCQ if it is a smooth point of the variety defined by the equality and active inequality constraints. Finding all solutions of the equations (18), (20) and (22) in the KKT conditions is computationally heavier than using numerical methods for finding the MLE, however it is still desirable since it provides the guarantee of finding the global maximum. Symbolic solving of the equations (18), (20) and (22) using Gröeber basis methods is possible only for small instances. An alternative is to use numerical algebraic geometry and homotopy continuation methods that find approximations of isolated complex solutions of a system of polynomial equations with probability one. This approach is taken in [11] for optimization problems in the life sciences. Furthermore, we suggest a "decomposition" of the KKT conditions into parts that makes solving them easier.
Let L be the ideal generated by the equations (18), (20) and (22) in the KKT conditions. For S ⊆ [m], let L S be the ideal generated by the Lagrange conditions for the optimization problem max F (x) subject to G i (x) = 0 for i ∈ S, H j (x) = 0 for j = 1, . . . , l.
Specifically, let L S be generated by the equations We denote by I S 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 2 Let L and L S be as defined above. Then where √ · denotes the radical of an ideal.
. Then it will follow that First take an element (µ, λ, x) of V (L). Let S be such that g i (x) = 0 for all i ∈ S. Then (µ S , λ, x) ∈ V (L S ), where µ S is the projection of µ to the coordinates in S. Conversely, let (µ S , λ, x) ∈ V (L S ). Let µ ∈ C m be such that µ i = (µ S ) i for i ∈ S and µ i = 0 otherwise. Then (µ, λ, x) ∈ V (L).
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 [7, Theorem 3.2.3], V (L ∩ C[x]) is the smallest algebraic variety containing π x (V (L)) and ) holds, because the right hand side is a variety and contains ∪π x (V (L S )) and hence π x (V (L)). On the other hand, Step 1: Let C = {}.
Step 2: For every S ⊆ [m], if dim(L S ) = 0 find V (L S ) and add all its elements 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.  (18), (20) and (22) in the KKT conditions. Since the global maximum satisfies the CRCQ, it must be a solution of these equations. By choosing among the real solutions that satisfy inequalities (19) and (21) in the KKT conditions the one that maximizes the value of the cost function F , we get the global optimum.
We are interested in the optimization problem, when the cost function is the log-likelihood function l u and the constraints are polynomials that define a statistical model. Although the equations (18) are not automatically polynomial for F = l u , they 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 equal to zero. For MLE, Gross and Rodriguez have defined a modification of the Lagrange conditions, known as Lagrange likelihood equations [12,Definition 2], whose solution set does not contain solutions with some variable equal to zero, if the original data does not contain zeros [12,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.

Remark 4
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.

Remark 5
In practice, it is crucial to know the ML degree, i.e. the degree of the ideal of KKT or Lagrange conditions. If the ideal of KKT or Lagrange conditions has a finite number of solutions, then the number of solutions is equal to the ML degree. 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 [1], NumericalAlgebraicGeometry package in Macaulay2 [18] and PHCpack [25]. Only PHCpack found all 92 solutions with initial settings and also running times of different programs differed by several hours.

Example 8
We aim to compute the MLE for the CFN model on K 1,3 and the data vector (100, 11,85,55,56,7,75,8). To do so, we relax the implicit characterization of the CFN model on K 1,3 given in Example 2 replacing strict inequalities by non-strict inequalities. We apply the modified version of Algorithm 1 described after Lemma 4. We automatically remove all ideals that contain a variable p ijk , a sum of two or four such variables. The code for this example can be found at the link:   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 [14,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 97 are real and 49 are positive. We list the seven points among them that have the highest value of the log-likelihood function in Table 2.
The first, second, fifth and sixth do not satisfy the relaxed inequalities and hence are not in the model. The third and fourth satisfy the inequalities of the relaxed model and hence the third one is the MLE of the model we defined at the beginning of this example. However, neither of them come from a parametrization by rate matrices. Furthermore, neither of them come from a parametrization by rate matrices where time is allowed to go to infinity. This would imply that somef (e) (g) and q ijk are zero. For this reason we have to consider a larger set of inequalities than in Proposition 1. Instead of picking any j and k such that v is on the path between them, we have to consider all such j and k. Although the third and fourth solution both satisfy in the Fourier coordinates the inequality q 000 q 110 − q 100 q 010 ≥ 0, they both fail the inequality q 000 q 101 − q 100 q 001 ≥ 0, which is obtained by choosing a different k in Proposition 1. If none of q ijk would be zero, then if one of the two inequalities is satisfied, also the other one is satisfied. This does not have to be true anymore when some q ijk = 0.
The seventh critical point 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 (100, 11,85,55,56,7,75,8) does not exist -the global maximum of the loglikelihood 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 ψ (e2) , which is equal to ψ (e2) = (−8.983, 8.983). 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.