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 continuous-time 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 non-symmetric 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 .

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 hill-climbing 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 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 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 taking 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 , 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 associated to every edge e that reflect the change in the probabilities when moving from one vertex to another. The transition matrices have the formwhere stands for matrix exponentiation, represents time and 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 is incorporated in the rate matrix .

To define a group-based phylogenetic model, we first fix an abelian group , a finite set of labels and a labeling function . Let . A rate matrix is a rate matrix in the group-based model if it satisfies for a vector with whenever . Hence transition matrices of the group-based model form a subset of all the transition matrices that satisfy for a probability vector with whenever . This is because the matrix exponentiation is defined as and if a matrix M has the structure given by and L, then one can check that also has the structure given by and L for all . The phylogenetic models we consider are symmetric, which means . In the case of group-based models, this is equivalent to for all .

We will assume that the root distribution of a group-based model is uniform or the root distribution is such that the matrix defined by 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 whenever and it has to satisfy inequalities that guarantee that the transition matrix P defined by is a matrix exponential of a rate matrix. The latter problem is called the embedding problem and is studied for 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 at the n leaves can be written as polynomials in the root probabilities and in the entries of the transition matrices. Denote by the vector of joint probabilities . As it is common in phylogenetic algebraic geometry, we will use the discrete Fourier transform for the groups and 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 whose elements are the group homomorphisms from to the multiplicative group of complex numbers of magnitude one. Given a function , its discrete Fourier transform is the function defined byIt is an invertible linear transformation given by the matrix K, where . The group-based model being symmetric is equivalent to the vectors and being real, see (Matsen 2009, Section 2). If we regard the vector of joint probabilities as a function of , i.e., as an element of , then the image of under the Fourier transform of is denoted .

The map from the entries of the rate matrices to the joint probabilities at leaves can be seen as a composition of four maps:

The map from to is given by the discrete Fourier transform of . It is an invertible linear transformation given by the matrix K.

The map from to is given by by (Matsen 2009, Lemma 2.2). It is an isomorphism between and .

In the case when root distribution puts all the mass at the identity, the map from to is given by by (Székely et al. 1993, Theorem 3), where and 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 . Otherwise . This follows from (Sturmfels and Sullivant 2005, Lemma 4 and formula (12)).

On the domain , this map is injective: (Matsen 2009, Proposition 3.3 and Proposition 3.4) give a map from to . Taking nonnegative square roots results in a left inverse to the map (3).

The map from to is given by the inverse of the discrete Fourier transform of . It is an invertible linear transformation given by the matrix , 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 (Cavender 1978; Farris 1973; Neyman 1971) on the rooted claw tree . This example has been previously studied in (Sturmfels and Sullivant 2005, Example 3) and (Hosten et al. 2005, Example 14). The CFN model is a group-based model with and the labeling function L defined by and . Denote the root distribution by and the transition matrices at edges bySince are probabilities, they are real numbers in [0, 1], and . Moreover, the restriction on the root distribution that it is uniform or defines a valid transition matrix in the CFN model gives and ; however, in Example 2 we will show that for the CFN model we can consider arbitrary root distributions. The determinant of is positive, because is the matrix exponential of a rate matrix . Conversely, for every satisfying these constraints, there exists a rate matrix such that by (Kingman 1962, Proposition 2).

The joint probabilities at the leaves have the parametrizationIn Sect. 3, we characterize this model in joint probabilities and without parameters . This is called the implicit description of a model. It consists of polynomial equations and inequalities in 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 , these binomials areas 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 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 be a polynomial ring. A subset 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 , consists of all the polynomials such that some power 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 or . The affine variety defined by S isLet be the ideal generated by , i.e., the smallest ideal containing . ThenA point of the variety is a smooth point if the Jacobian of has maximal possible rank. Otherwise a point of the variety is called singular. Let T be a subset of . The Zariski closure 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 . 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 .

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 assign state to leaf x for and zero to all other leaf vertices. Then

Proposition 2

(Matsen (2009), Proposition 3.4) Given some internal edge e, let the two vertices incident to e be v and . Let i, j (respectively, ) be leaf vertices such that the path from i to j (respectively, the path from to ) contains v but not (respectively, but not v). Let assign state to leaf x for and zero to all other leaf vertices. Then

The next proposition will summarize the procedure in Matsen (2009) to construct inequalities that describe joint probabilities. We will denote by the row of the matrix labeled by g and by the Laurent monomial .

Proposition 3

Assume that the labeling function L satisfies for all . Consider the set of that satisfies , whenever and for all nonzero . The images of this set under the maps in (1) are:

The constraints for are obtained by substituting by in the constraints for . In particular, this gives , whenever and for all nonzero .

The constraints for are , whenever , for all nonzero and for all . This equation and inequalities are equivalent to , whenever , for all nonzero and for all . Here we have squared the inequalities .

The constraints for are given by phylogenetic invariants, equation , inequalities and inequalities that are obtained by substituting expressions for in Propositions 1 and 2 to inequalities in the previous item.

The constraints for are obtained by substituting by in the constraints for .

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.

Example 2

We will derive the implicit description of the CFN model on the rooted claw tree . We start with the case when and . In addition to phylogenetic invariants in Example 1, applying Proposition 3 gives the following inequalities in Fourier coordinates:The inequality is for the hidden leaf corresponding to the root. Since , we can multiply all the denominators by without changing the inequalities (4). Clearing denominators gives the following polynomial inequalitiesBy applying the discrete Fourier transformation, we get the implicit description in the original coordinatesIf and , we can switch 0 and 1 and apply the previously considered case. We obtain the implicit description by switching 0 and 1 in the subindices of the equations and inequalities (5)–(19). This operation leaves all the equations and the inequalities (11), (13), (14), (16)–(19) the same. It changes the inequalities (9), (10), (12) and (15). Explicitly, the implicit description isFinally, we consider the case when the root distribution is uniform. By (Sturmfels and Sullivant 2005, Proposition 31), one gets additional phylogenetic invariantsMoreover, we have the following equation and inequalities:In original coordinates, we get the following implicit description:The implicit description of the CFN model on the tree 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 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 for some .

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 satisfies one of the non-strict inequalities with equality if and only if it has a preimage that satisfies one of the inequalities with equality.

Example 3

We consider the CFN model. A joint probability vector satisfying the assumptions of Corollary 1 has in its parametrization the rate matrix for some . The transition matrix corresponding to the same edge is .

Maximum Likelihood Estimation via Numerical Algebraic Geometry

In this section, we use the terminology and notation introduced in Sect. 2. In particular, are the joint probability distributions at the n-leaves. Let be a vector of observations at leaves. The log-likelihood function of a phylogenetic model isMaximum 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 .

Example 4

In (Hosten et al. 2005, Example 14), maximum likelihood estimation on the Zariski closure of the CFN model on 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 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 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 problemIf is a local optimum and the optimization problem satisfies first-order constraint qualifications, then there exist , where , and , where , such that 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 problemIf is a local optimum of the optimization problem (55), then there exist , where , such that 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 be the ideal generated by the polynomials on the left-hand sides of the equations (50), (52) and (54) in the KKT conditions. For , let be the ideal generated by the polynomials in the Lagrange conditions for the optimization problemSpecifically, let be generated by the polynomialsWe denote by the ideal generated by the constraints in the above optimization problem, i.e., .

Theorem 1

Let L and 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).

First take an element of V(L). Let S be such that for all . Then , where is the projection of to the coordinates in S. Conversely, let . Let be such that for and otherwise. Then .

We have shown that , where is the projection of or on x. By the Closure Theorem (Cox et al. 1992, Theorem 3.2.3), is the smallest algebraic variety containing and is the smallest algebraic variety containing . The inclusion holds, because the right-hand side is a variety and contains and hence . On the other hand, since for every S, also for every S. Hence .

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).

Corollary 2

If V(L) is finite and the global maxima of the optimization problem (49) satisfy CRCQ, then Algorithm 1 outputs the global maxima.

Theorem 1 implies that . The variety V(L) being finite implies that and hence all are finite. Hence after Step 2, the list contains all solutions of Eqs. (50), (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 log-likelihood function and the constraints are polynomial equations and inequalities that describe a statistical model (written as for and for , respectively). Although Eq. (50) is not polynomial for , it can be made polynomial by multiplying the equationwith the variable .

One of the reasons why the variety 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.

Lemma 1

Assume that the ideal decomposes as , where and . If satisfies the Lagrange conditions for the optimization problem max F(x) subject to for , then satisfies the Lagrange conditions for the optimization problem max F(x) subject to for or for the optimization problem max F(x) subject to for .

Since , we have . Hence the optimization problem max F(x) subject to for is equivalent to max F(x) subject to for . The Lagrange conditions for the latter optimization problem areIf there exists k such that , then we must have for . Hence satisfiesDefining , we see that satisfies Lagrange conditions for the optimization problem max F(x) subject to for . Otherwise for and satisfies Lagrange conditions for the optimization problem max F(x) subject to for .

Lemma 2

Let and . If satisfies the Lagrange conditions for the optimization problem max F(x) subject to the generators of , then satisfies the Lagrange conditions for one of the optimization problems max F(x) subject to the generators of , where .

Assume , , and . Then and . The Lagrange conditions for the generators of areIf there exists such that and such that , then we must have for and for . Hence satisfiesDefining and , we see that satisfies Lagrange conditions for the optimization problem max F(x) subject to the generators of . If for all k and/or for all k, then we get other combinations , or .

Lemma 1 suggests that if S is a singleton in Step 2 of Algorithm 1, then we can replace the ideal of Lagrange conditions for in Step 2 of Algorithm 1 by the ideals of Lagrange conditions for minimal primes of . If , then . Hence by Lemmas 1 and 2, we can replace the ideal by the ideals of Lagrange conditions for the sum of minimal primes of , where .

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 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 for every .

Remark 4

In practice, it is crucial to know the degrees of the ideals 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 , 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 . 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 and the data vector . This data vector is obtained by generating 100 samples from the distribution inside the CFN model with rate parametersThe 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.

Fig. 1

The tree in Example 5 has two edges with short branch lengths and , one edge with a long branch length M 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 are given by the left-hand sides of the inequalities (9)–(19) and the polynomials are given by the left-hand sides of Eqs. (5)–(8); in the second case, the polynomials are given by the left-hand sides of the inequalities (24)–(34) and the polynomials are given by the left-hand sides of Eqs. (20)–(23); in the third case, the polynomials are given by the left-hand sides of the inequalities (41), (43), (44), (46)–(48) and the polynomials are given by the left-hand sides of Eqs. (35)–(40), (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 and are the same for the first two optimization problems; in the third optimization problem there is one polynomial less and some polynomials are among polynomials , 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 and 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 . 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.

Table 1

Table summarizing different boundary components

Dim I	Degree L	# of ideals
5	92	1
4	9	4
4	1	8
3	1	24
2	1	6
1	1	1
Total	167	44

Table 2

Critical points with the highest values of the log-likelihood function

p	l_u	MLE
(.183, .051, .256, .055, .147, .053, .204, .052)	- 188.451	No
(.183, .049, .243, .065, .156, .042, .207, .055)	- 188.722	No
(.191, .053, .243, .042, .156, .065, .199, .051)	- 188.803	No
(.165, .05, .23, .055, .165, .05, .23, .055)	- 188.927	No
(.17, .045, .225, .06, .17, .045, .225, .06)	- 189.042	No
(.174, .059, .221, .046, .174, .059, .221, .046)	- 189.303	No
(.22, .05, .22, .05, .175, .055, .175, .055)	- 189.488	Yes

The first six critical points in Table 2 satisfyandHence these critical points are not in the CFN model on 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 and the data vector 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 , which is equal to . Without having the implicit description of the CFN model on 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 are almost zero. We recall that the semialgebraic description of the CFN model includes strict inequalities . The global maximum of the likelihood function on the closure of the CFN model on satisfies . 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.