Computing linear extensions for polynomial posets subject to algebraic constraints
CComputing linear extensions for Boolean lattices withalgebraic constraints
Shane Kepley ∗ , Konstantin Mischaikow † , and Lun Zhang ‡ Abstract
In this paper we consider the classical problem of computing linear extensions of a givenposet which is well known to be a difficult problem. However, in our setting the elementsof the poset are multi-variate polynomials, and only a small “admissible” subset of these lin-ear extensions, determined implicitly by the evaluation map, are of interest. This seeminglynovel problem arises in the study of global dynamics of gene regulatory networks in which casethe poset is a Boolean lattice. We provide an algorithm for solving this problem using linearprogramming for arbitrary partial orders of linear polynomials. This algorithm exploits thisadditional algebraic structure inherited from the polynomials to efficiently compute the admis-sible linear extensions. The biologically relevant problem involves multi-linear polynomials andwe provide a construction for embedding it into an instance of the linear problem.
Consider a set of real polynomials P , defined on a domain Ξ ⊂ R d , equipped with a partial order ≺ .We are interested in identifying linear extensions (total orders compatible with ≺ ) that are satisfiedby P under evaluation at a point in Ξ.To be more precise consider a semi-algebraic set Ξ ⊂ R d , called the evaluation domain , and acollection of polynomials P := { p , . . . , p K } ⊂ R [ x , . . . , x d ]. Let ≺ denote a partial order on P such that if p ≺ q , then p ( ξ ) < q ( ξ ) for all ξ ∈ Ξ . (1)Let S K +1 denote the set of permutations on K + 1 symbols. We identify linear extensions of P with a subset of S K +1 as follows. Given σ ∈ S K +1 , let ≺ σ denote the linear order p σ (0) ≺ σ p σ (1) ≺ σ · · · ≺ σ p σ ( K ) . We define the realizable set associated to σ byΞ σ := (cid:8) ξ ∈ Ξ : p σ ( k ) ( ξ ) < p σ ( k +1) ( ξ ) for all 0 ≤ k ≤ K − (cid:9) . (2) ∗ [email protected] † [email protected] ‡ [email protected] a r X i v : . [ m a t h . C O ] J un bserve that if Ξ σ (cid:54) = ∅ , then ≺ σ is a linear extension of ≺ . The algebraically constrained linearextension problem (AC-LEP) defined by ( P , ≺ , Ξ) is to determine T ( P , ≺ , Ξ) := { σ ∈ S K +1 : Ξ σ is nonempty } . Notice that in the formulation of the AC-LEP we have identified the partial order ≺ σ with theassociated element in S K +1 . We will use this identification throughout the remainder of the paper.We say that each σ ∈ T ( P , ≺ , Ξ) is an admissible linear extension.As is discussed in Section 2 our motivation for study the AC-LEP comes from modeling thedynamics of regulatory networks in biology and in particular characterizing relevant subsets ofparameter space. For the moment we attempt to put this problem into a broader mathematicalcontext as the problem itself, as well as our solutions for some special cases, have elements of bothclassical real algebraic geometry and order theory.
Quantifier elimination and real algebraic geometry
Observe that if Ξ = R d , P is an arbitrary collection of polynomials, then σ ∈ S K +1 is admissibleif and only if there exists ξ ∈ R d such that p σ ( k ) ( ξ ) − p σ ( k +1) ( ξ ) < ≤ k ≤ K . Theseinequalities define a semi-algebraic set and therefore, taking ≺ to be the trivial partial order (i.e. P is a single anti-chain), this instance of AC-LEP is equivalent to the classical problem of decidabilityfor real semi-algebraic sets.The previous example illustrates a major challenge in solving the AC-LEP. The first generalalgorithm for solving the quantifier elimination/decidability problem for polynomials in R d withfeasible running time was the cylindrical algebraic decomposition (CAD) given by Collins [9] in1975. The CAD algorithm works by subdividing Ξ into subsets on which the polynomials are signinvariant. Given such a decomposition, decidability is reduced to simply evaluating each polynomialat a sample point located in each subset and checking if it satisfies the necessary inequalities.Unfortunately, the computational complexity of the algorithm grows like O (cid:16) (2 D ) d − ( K + 1) d − d − − (cid:17) where D = max { deg p : p ∈ P} . This worst case running time is known to be sharp even for classes of “nice” polynomials e.g. linear[8], and moreover, the worst case is also typical [2]. As a result, the question of whether or not σ ∈ T ( P , ≺ , Ξ), even for a single σ ∈ S K +1 , is often intractable for problems of practical interest.In addition, this algorithm does not provide partial information. It either runs to completion,in which case it is guaranteed to provide an answer, or it provides no information. Furthermore, wenote that if additional algebraic constraints are added e.g. we assume Ξ ⊂ Ξ ⊂ R d is a strict semi-algebraic subset, then the CAD algorithm can handle this by simply appending the polynomialconstraints which define Ξ to the set of polynomials. However, this dramatically increases thecomplexity of the CAD algorithm, despite the fact that the number of admissible linear extensionscan only decrease.Some improved algorithms have been proposed which aim to reduce the complexity of specificaspects of the problem or for special classes of polynomials (e.g. [15, 6, 7]). These improvementsoften provide dramatic algorithmic speedups for checking whether σ ∈ T ( P , ≺ , Ξ) for a singlelinear extension. However, these algorithms have the same worse case running time as the generalCAD algorithm and understanding which classes of polynomials benefit is still a very active areaof research. Therefore, these improved algorithms alone are not sufficient to handle instances of2C-LEP since we are interested in determining which of the ( K +1)! possible semi-algebraic sets arenonempty. An efficient algorithm which does not produce a decomposition of Ξ into sign invariantsubsets would still need to be called ( K + 1)! times. However, we will make use of these improvedalgorithms as a post-processing step which we discuss further in Section 4. Computing linear extensions of Boolean lattices
Let us momentarily ignore the algebraic structure in the AC-LEP by forgetting that P is a collectionof polynomials. Hence, we focus only on the poset structure ( P , ≺ ), and consider the problem ofcomputing all linear extensions. The related problem of counting all linear extensions of a partialorder is a well studied problem in order theory. Its importance is due in large part to its connectionwith the complexity of sorting elements in a list. If one considers a list of ( K + 1) distinct valueswhich have been partially sorted by making pairwise comparisons on a subset of its elements, thenthese comparisons induce a partial order. Therefore, the linearly ordered values of the fully sortedlist are given by one of the possible linear extensions for the partial order. As a result, the complexityof completely sorting a list is intimately connected to counting linear extensions for posets.Observe that computing the set of all linear extensions of ( P , ≺ ) is not easier than counting themwhich is known to be P -complete [4]. In particular, a polynomial time algorithm for computingall possible linear extensions for arbitrary posets would imply that P = N P by Toda’s theorem[21]. Moreover, we are interested not only in counting linear extensions, but explicitly computingthem. Therefore, we are also concerned with how fast the number of admissible linear extensionsgrows.For reasons we discuss in Section 4, we are specifically interested in the case that ( P , ≺ ) is aBoolean lattice. Specifically, for fixed n ∈ N , define S n := { , . . . , n } and let 2 S n denote its powerset. The standard n -dimensional Boolean lattice is the poset, (2 S n , ≺ B ), where ≺ B is the partialorder defined by inclusion. We say a poset, ( P , ≺ ), is an n -dimensional Boolean lattice if ( P , ≺ ) isorder isomorphic to the standard n -dimensional Boolean lattice and we write ≺ B in place of ≺ .Estimating the number of linear extensions for Boolean lattices was first considered in [17] whichestablished a nontrivial upper bound. Later, Brightwell and Tetali [5] proved the following resultthat essentially settles the question for all practical considerations. If Q ( n ) denotes the number oflinear extensions of an n -dimensional Boolean lattice, thenlog Q ( n )2 n = log (cid:18) n (cid:98) n / (cid:99) (cid:19) −
32 log e + O (cid:18) ln nn (cid:19) . (3)The estimate in Equation (3) illustrates a major challenge in solving the instances of AC-LEPof interest in this paper. Consider an instance of AC-LEP given by ( P , ≺ B , Ξ) where ( P , ≺ B ) is an n -dimensional Boolean lattice and Ξ is any evaluation domain.Suppose we had a black box for efficiently computing all linear extensions of a Boolean latticedenoted by L ⊂ S K +1 . Furthermore, assume we also had a “CAD”-like algorithm which couldefficiently check if Ξ σ (cid:54) = ∅ . Then, one would need to call this algorithm only L -many times asopposed to ( K + 1)! as we argued above. However, the growth estimate in Equation (3) impliesthat the number of calls to this algorithm would still grow exponentially. This work
In this work we present efficient algorithms for solving two specific instances of the AC-LEP. Thefirst is the linearly constrained linear extension problem (LC-LEP), in which P is a set of linear3olynomials, Ξ is a polytope, and ≺ is an arbitrary partial order. We present an efficient algorithmfor solving the LC-LEP in Section 3. The second instance of the AC-LEP, which we call the parameter space decomposition (PSD) problem and describe now (see Definition 4.6 for a precisedefinition), is motivated by an application from systems biology described in Section 2. Definition 1.1.
For n ∈ N , an interaction function of order n is a polynomial in n variables, z = ( z , . . . , z n ), of the form f ( z ) = q (cid:89) j =1 f j ( z ) (4)where each factor has the form f j ( z ) = (cid:88) i ∈ I j z i and the indexing sets { I j : 1 ≤ j ≤ q } form a partition for { , . . . , n } . We define the interactiontype of f to be n := ( n , . . . , n q ) where n j denotes the number of elements in I j . Remark 1.2.
We leave it to the reader to check that the order of the indexing sets I j does notmatter for any of the analysis in this paper. Therefore, for convenience in reporting results (seeSection 5) we will always assume that n ≥ n ≥ · · · ≥ n q . To define an instance of the PSD problem, fix an interaction function f of order n and let P be the collection of polynomials in the 2 n positive real variables, { (cid:96) i , δ i : 1 ≤ i ≤ n } , obtained byevaluating f ( z ) with each z i ∈ { (cid:96) i , (cid:96) i + δ i } . Taking all possible combinations of z i for 1 ≤ i ≤ n produces the polynomials for the PSD problem, P = { p , . . . , p n − } ⊂ R [ (cid:96) , . . . , (cid:96) n , δ , . . . , δ n ] . (5)In Section 4, we will define an indexing map between the 2 n elements of P and the standard n -dimensional Boolean lattice. Let ≺ B denote the Boolean lattice partial order with respect to thisindex map, and set Ξ = (0 , ∞ ) n . The PSD problem is the instance of the AC-LEP defined by( P , ≺ B , Ξ). In Section 4, we prove that ( P , ≺ B , Ξ) satisfies Equation (1). However, we presentsome examples before continuing.
Example 1.3.
The simplest nonlinear PSD problem arises from the interaction function f ( z ) = ( z + z ) z which has interaction type, n = (2 , p = ( (cid:96) + (cid:96) ) (cid:96) p = ( (cid:96) + (cid:96) + δ ) (cid:96) p = ( (cid:96) + (cid:96) )( (cid:96) + δ ) p = ( (cid:96) + (cid:96) + δ )( (cid:96) + δ ) p = ( (cid:96) + (cid:96) + δ ) (cid:96) p = ( (cid:96) + (cid:96) + δ + δ ) (cid:96) p = ( (cid:96) + (cid:96) + δ )( (cid:96) + δ ) p = ( (cid:96) + (cid:96) + δ + δ )( (cid:96) + δ ) . The PSD evaluation domain is Ξ = (0 , ∞ ) and the partial order, ≺ B , is imposed on P by identifying p i with the vertex of a unit cube whose coordinates are ( i ) ∈ F where ( i ) is the binary expansion4f i . The solution to this PSD problem is the set of admissible linear extensions of ( P , ≺ B ), suchthat σ ∈ T (cid:0) P , ≺ B , (0 , ∞ ) (cid:1) if and only if Ξ σ (cid:54) = ∅ . We note that there are 8! = 40 ,
320 linear orderson P . However, only 48 of these are linear extensions of ≺ B , and of these, only the following 20linear extensions are admissible.(0 , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,
7) (0 , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,
7) (0 , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,
7) (0 , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ξ ∈ (0 , ∞ ) , if p ( ξ ) < p ( ξ ), then p ( ξ ) < p ( ξ ) must also hold.Unlike the partial order which constrains all possible linear extensions, this order relation isconditional. Indeed, there exist choices of ξ such that p ( ξ ) > p ( ξ ) in which case there is norequirement imposed on the order of p ( ξ ) , p ( ξ ), and in fact, there are admissible linear extensionswhich satisfy both choices, e.g. the first two orders in column four. As another example, observethat p ( ξ ) < p ( ξ ), if and only if p ( ξ ) < p ( ξ ).This algebraic constraint is bi-conditional, however, it also can not be represented in the partialorder since both choices occur in at least one admissible order.To emphasize the role of the interaction function in determining the algebraic constraints, weconsider a similar PSD problem that is also an instance of LC-LEP. Example 1.4.
Consider the interaction type, n = (1 , ,
1) with corresponding interaction function f = z + z + z . As in Example 1.3, we obtain 8 PSD polynomials given explicitly by p = (cid:96) + (cid:96) + (cid:96) p = (cid:96) + (cid:96) + (cid:96) + δ p = (cid:96) + (cid:96) + (cid:96) + δ p = (cid:96) + (cid:96) + (cid:96) + δ + δ p = (cid:96) + (cid:96) + (cid:96) + δ p = (cid:96) + (cid:96) + (cid:96) + δ + δ p = (cid:96) + (cid:96) + (cid:96) + δ + δ p = (cid:96) + (cid:96) + (cid:96) + δ + δ + δ The evaluation domain and partial order are identical to the PSD problem in Example 1.3. Never-theless, only the following 12 linear extensions are admissible(0 , , , , , , , , , , , , , , , , , , , , ,
7) (0 , , , , , , , , , , , , , , , , , , , , ,
7) (0 , , , , , , , , , , , , , , , , , , , , ,
7) (0 , , , , , , , , , , , , , , , , , , , , , without first computing the linear extensions of ( P , ≺ B ) and thenchecking which are admissible. Related work
As is indicated above our original motivation for this paper arises from problems in systems biologyfor which explicit complete solutions to the PSD problem are required. As such the majority ofthis introduction has focused on the question of efficacy of computation. However, there is anotherdirection in which the work of this paper overlaps with other efforts. In particular, observe thatthe case where the interaction function is linear, i.e. has interaction type n = (1 , . . . , n = 1 , . . . ,
7. Ourcomputations (see Table 1) lead to the same numbers, as expected.After accounting for symmetry in the number of linear PSD solutions for interaction types(1 , , , , , , , , , , , , a , a , a which align with sequence A009997 in the OEIS [16]. From [11], we know this sequence representsthe number of comparative probability orderings on all subsets of n elements that can arise byassigning a probability distribution to the individual elements. The equivalence of comparativeprobability orderings and solutions to the linear PSD problem follows directly from the definitionof comparative probability. Organization of paper
The remainder of this paper is organized as follows. In Section 2 we briefly describe how the PSDproblem arises naturally in the study of global dynamics for gene regulatory networks. In Section3, we present an efficient algorithm for solving instances of the LC-LEP. In Section 4, we show thatthe LC-LEP is related to the PSD problem in the following way. If (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) is an instanceof the PSD problem, then we construct an associated instance of LC-LEP, denoted by ( P (cid:48) , ≺ B , Ξ (cid:48) ),which satisfies the inclusion T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) ⊆ T ( P (cid:48) , ≺ B , Ξ (cid:48) ) . (6)We refer to this instance of LC-LEP as the linearized PSD problem associated to (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) .We exploit this construction and the algorithm for solving the LC-LEP presented in Section 3, toprovide a means of efficiently computing a collection of candidates that contains the solution to thePSD problem. We prove that in some cases the inclusion in Equation (6) is actually an equality.More generally, this inclusion is strict, but the candidate set is a sparse subset of the collection of alllinear extensions of ( P , ≺ B ). In this case we describe algorithms for removing the non-admissiblesolutions. 6inally, in Section 5 we present the results for all PSD solutions with order up to four. Addi-tionally, we have some results for PSDs or order five and size. For the remaining cases and PSDsof higher order the computations become too large. This section provides a brief description of how the AC-LEP arises in the context of mathematicalmodeling of problems from systems biology. Biologists often describe regulatory networks in termsof annotated directed graphs, such as that shown in Figure 1(a) where the labeling of the edges, n → m or n (cid:97) m indicates whether node n activates or represses node m . Our goal is to describethe type of dynamics that can be expressed by the regulatory network. This requires imposing amathematical interpretation on the regulatory network that is compatible with its use as a biologicalmodel. With this in mind, we assign to node m a state variable, x m >
0, that corresponds withthe quantity of a protein, mRNA, or a signaling molecule. Precise nonlinear expressions for theinteractions of the variables are not assumed to be known, but we do assume that the sign of therate of change of x m is determined by the expression − γ m x m + Λ m ( x ) (7)where γ m indicates the decay rate and Λ m is a parameter dependent function that characterizesthe rate of growth of x m . Note that Λ m is a function of x n if and only if there exists an edge from n to m in the regulatory network.(a) (b) x n (cid:96) m,n (cid:96) m,n + δ m,n θ m,n Figure 1: (a) Example of a regulatory network. (b) Model for edge n → m where (cid:96) m,n indicateslow level of growth rate of x m induced by x n and (cid:96) m,n + δ m,n indicates high level of growth rateof x m induced by x n . θ m,n provides information about the value of x n that lies between low valuesinducing low and high expression levels.Since the biological model provides minimal information about the effect of x n on x m we wantto choose a mathematical expression with a minimal set of assumptions. The rates of growth of x m due to x n are labeled as 0 < (cid:96) m,n < (cid:96) m,n + δ m,n . Figure 1(b) corresponds to an edge n → m and θ m,n indicates that the rate of increase (cid:96) m,n must occur at some lesser value of x n and the rate of7ncrease (cid:96) m,n + δ m,n must occur at some greater value of x n . An arrow of the form n (cid:97) m leads tothe opposite relation.This introduces three positive parameters, (cid:96) m,n , δ m,n , and θ m,n , for each edge in the regulatorynetwork. Note that this is the minimal number of parameters that allows one to quantify theassumption that x n activates x m (or equivalently that x n represses x m ). We encode this informationwith the following functions λ + m,n ( x n ) = (cid:40) (cid:96) m,n if x n < θ m,n (cid:96) m,n + δ m,n if x n > θ m,n and λ − m,n ( x n ) = (cid:40) (cid:96) m,n + δ m,n if x n < θ m,n (cid:96) m,n if x n > θ m,n .We do not assume that the values of (cid:96) m,n , δ m,n , or θ m,n are known. This is intentional as many ofthese parameters do not have an easy biological interpretation and/or correspond to physical con-stants which are difficult or impossible to precisely measure. Thus, the goal is not to determine thedynamics at any choice of parameters, but to determine the range and robustness of the qualitativedynamics exhibited by a network.A regulatory network such as that of Figure 1(a) does not indicate how multiple inputs to aparticular node should be processed. An approach that is used is to assume a simple algebraicrelationship involving sums and products of the λ ± . As an example, a reasonable choice for x ofFigure 1(a) is Λ ( x , x , x ) = (cid:0) λ + ( x ) + λ + ( x ) (cid:1) λ − ( x ) . (8)Observe that each λ ± takes only two values and therefore, generically Equation (8) takes 8 distinctvalues which are precisely the values of the elements of P in the PSD of Example 1.3.As is suggested in the caption of Figure 1(b), we do not interpret the values of λ ± , or Λ asliteral expressions of the nonlinear interactions, but rather, that the associated parameter valuesare landmarks of whatever of the “true” nonlinear function is. This has several consequences.1. We cannot expect Equation (7) to provide precise information about the growth rate of x m .Therefore we restrict our attention to asking whether x m is increasing or decreasing. However,we wish to answer this question over all the possible parameter values γ , θ , (cid:96) and δ .2. The only values of x m at which the dynamics of x n change are of the form θ m, ∗ . The associatedhyperplanes x j = θ k,j decompose phase space [0 , ∞ ) N , where N is the number of nodes inthe network, into N -dimensional rectangular regions called domains .3. Since we have determined T (cid:0) { p , . . . , p } , ≺ , (0 , ∞ ) (cid:1) for (8) we can determine all possiblesigns of (7) associated with (8) by cataloguing the relative values of γ θ j, , j = 4 , { p , . . . , p } .The Dynamic Signatures Generated by Regulatory Networks (DSGRN) library is software that,given the information of the form provided by consequence 3, is capable of efficiently building adatabase of the global dynamics of a regulatory network over all of parameter space [10, 14]. In[10] the authors considered networks whose nodes have at most three in-edges, and at most threeout-edges. This constraint was due to the difficulty in solving the PSD problem and the resultsof this paper provide a means to vastly expand the capabilities of DSGRN [14]. In particular,DSGRN can now handle the algebraic combinations of 4 to 6 in-edges as is indicated in Section 5,and arbitrarily many out-edges. 8 Solving the LC-LEP
In this Section, we provide an efficient algorithm to solve the LC-LEP defined in Section 1. Notethat if q ∈ R [ x , . . . , x d ] is a linear polynomial, then evaluation of q defines a linear functionalon R d . Thus, there exists a unique vector u q ∈ R d , that we call the representation vector for q ,satisfying q ( ξ ) = u q · ξ for all ξ ∈ R d . Equivalently, u q is the vector of coefficients of q . Since the evaluation domain for the LC-LEP is apolytope, there exists a collection of linear polynomials Q Ξ , such thatΞ = (cid:8) ξ ∈ R d : ξ · u q > q ∈ Q Ξ (cid:9) . (9)We assume that (9) is satisfied for the remainder of this section.To foreshadow our approach recall that by definition σ ∈ T ( P , ≺ , Ξ) if and only if Ξ σ (cid:54) = ∅ .Our approach to determining the latter is to recast it in the language of linear algebra on cones in R d . From this perspective, the problem is equivalent to rigorously solving a linear programmingproblem and the efficacy of our algorithm is based on the fact that this can be done efficiently.With this goal in mind, we begin with a few remarks concerning cones and ordered vector spaces. Definition 3.1.
A subset C ⊂ R d is a cone if v ∈ C and θ ∈ [0 , ∞ ) implies that θv ∈ C . The cone C is pointed if it is closed, convex, and satisfies C ∩ − C = C ∩ {− v : v ∈ C } = 0 . (10)Observe that (10) implies that a pointed cone does not contain any lines. A vector v ∈ R d isa conic combination of { v , . . . , v k } ⊂ R d if v = θ v + · · · + θ k v k where θ , . . . , θ k ≥
0. Suppose V = { v , . . . , v k } ⊂ R d . The conic hull of V is given bycone( V ) := (cid:40) k (cid:88) i =1 θ i v i : 0 ≤ θ i , i = 1 , . . . , k (cid:41) . The following result is left to the reader to check.
Proposition 3.2.
Given V = { v , . . . , v k } ⊂ R d , cone( V ) is the smallest closed convex cone thatcontains V . We make use of the following propositions.
Proposition 3.3.
Suppose V = { v , . . . , v m } ⊂ R d is a collection of nonzero vectors such that cone( V ) is a pointed cone. Then, there exists some v (cid:48) ∈ R d such that v (cid:48) · v i > for all ≤ i ≤ m .Proof. Observe that − v m / ∈ cone( { v , . . . , v m − } ) ⊂ cone( V ) since cone( V ) is pointed. Hence {− v m } and cone( { v , . . . , v m − } ) are disjoint, convex, closed subsets of R d .Therefore, by the hyperplane separation theorem [3], there exists v (cid:48) ∈ R d such that v (cid:48) · − v m < v (cid:48) · v > v ∈ cone( { v , . . . , v m − } ). 9 roposition 3.4. Suppose V = { v , . . . , v m } ⊂ R d is a collection of nonzero vectors such that cone( V ) is a pointed cone. If − v / ∈ cone( V ) , then cone( V ∪ { v } ) is pointed.Proof. Suppose that cone( V ∪ { v } ) is not pointed. Then, there exists w (cid:54) = 0 such that w, − w ∈ cone( V ∪ { v } ) or equivalently w = m (cid:88) i =0 α i v i + αv and − w = m (cid:88) i =0 β i v i + βv where α i , β i , α, β are all nonnegative. Note that if α = β = 0, then ± w ∈ cone( V ), which contradictsthe assumption that V is pointed. The sum of the two equations above is − ( α + β ) v = m (cid:88) i =0 ( α i + β i ) v i . This implies that − v ∈ cone( V ), contradicting the assumption that cone( V ) is pointed.The previous propositions illustrate the importance of solving the cone inclusion problem: givena vector v ∈ R d , and finite set of vectors V ⊂ R d , determine whether or not v ∈ cone( V ). Algorithm1, stated below, solves this problem. Observe that checking if v ∈ cone( V ) is equivalent to solvingthe following linear programming feasibility problem.Does there exist α such that V α = v (11)and α ≥ V is the column matrix of V .Linear programming is a powerful tool that is widely used in convex optimization, and as a result,there are many available solvers/algorithms for solving the linear programming feasibility problem[22]. The results can be made rigorous by performing computations using interval arithmetic [19]or rational linear programming [1] in the case that V is rational. Observe that the PSD problemdefined in Section 1 satisfies this constraint. As is made clear in Section 5, we use different solversdepending on the machine employed to do the computations. In any case, we take for granted theexistence of a rigorous solver for the feasibility problem in Equation 11 as a “black box” which wecall LPSolver which is employed in the following algorithm.
Algorithm 1:
Cone inclusion
Input: v, V = { v , . . . , v m } , LPSolver
Output: True , FalseResult:
Return
True if v ∈ cone( V ) otherwise False Function
InCone( v, V,
LPSolver) : Return
LPSolver( v, V ) End Function
The next algorithm uses Proposition 3.4 (see line 4) and Algorithm 1 to determine if a set of10ectors defines a pointed cone.
Algorithm 2:
Cone pointedness
Input: V = { v , . . . , v m } Output: True , FalseResult:
Return
True if cone( V ) is pointed otherwise False Function
CheckCone( V ) : V (cid:48) = { v } for i = 2 ... m do if InCone( − v i , V (cid:48) ) then Return
False else V (cid:48) = V (cid:48) ∪ { v i } end end Return
True End Function
We now show that the LC-LEP can be reformulated as a problem of identifying whether somespecific subsets of vectors generate pointed cones.
Definition 3.5.
Given an instance of the LC-LEP, ( P , ≺ , Ξ), we define the base cone as cone( V ) :=cone( V Ξ ∪ V ≺ ) where V Ξ and V ≺ are defined as follows. Set V Ξ := { u q : q ∈ Q Ξ } where Q Ξ are the representation vectors as defined in Equation (9). Applying Algorithm 2 to V Ξ (and the fact that we assume Ξ (cid:54) = ∅ ) shows that cone( V Ξ ) is pointed. Define V ≺ := { u : u is the representing vector of p − q where q ≺ p and p, q ∈ P} . Observe that if ( P , ≺ ) satisfies Equation (1), then the representation vector for p − q is an elementof V ≺ by definition. Therefore, by Proposition 3.3, V ≺ is pointed.The motivation behind our definition of V is that it characterizes the algebraic constraints inthe LC-LEP in terms of linear algebra that can be efficiently checked. The next proposition showsthat the same idea works for linear extensions.Given σ ∈ S K +1 , we define V σ := V ∪ (cid:8) u p σ ( i +1) − u p σ ( i ) : p σ ( i ) ∈ P , i = 0 , . . . , K − (cid:9) . (12) Proposition 3.6.
For any σ ∈ S K +1 , Ξ σ (cid:54) = ∅ if and only if cone( V σ ) is a pointed cone.Proof. First we assume Ξ σ (cid:54) = ∅ and that ξ ∈ Ξ σ . On one hand, if cone( V σ ) is not pointed, thenthere are nonzero vectors − v, v ∈ cone( V σ ). However, the definition of Ξ σ implies that − v · ξ > v · ξ >
0, which is a contradiction.On the other hand, if cone( V σ ) is pointed, then by Proposition 3.3 there exists a ξ ∈ R d suchthat ξ · v >
0, for all v ∈ V σ . Thus, the definition of V σ implies that ξ ∈ Ξ σ and hence Ξ σ (cid:54) = ∅ .We emphasize that the importance of Proposition 3.6 is the implied equivalence T ( P , ≺ , Ξ) = { σ : Ξ σ (cid:54) = ∅} = { σ : cone( V σ ) is pointed } . .2 An algorithm for identifying T ( P , ≺ , Ξ) In this section we present an algorithm for solving an arbitrary instance of the LC-LEP.
Algorithm 3:
LC-LEP solver
Input: σ part = [ ] , P , V = V , Ret = {} Output: T ( P , ≺ , Ξ) Result:
Ret : collection of all linearly realizable total order under restriction of V Function
OrderingGenerator( σ part , P , V, Ret) : if σ part == [ ] and CheckCone( V ) is not True then Return end l + 1 = length of σ part if l == K then add σ part to Ret Return end for i = 0 .. K do if i (cid:54)∈ σ part then u (cid:48) = u p i − u p σ part( l ) if not InCone( − u (cid:48) , V ) then OrderingGenerator( σ part + [ i ] , P , V ∪ { v (cid:48) } , Ret) end end end End Function
To prove the correctness of the algorithm it is useful to denote the return of Algorithm 3 giveninput ( P , ≺ , Ξ) as T alg ( P , ≺ , Ξ).
Definition 3.7.
For fixed ( P , ≺ , Ξ), σ ∈ S K +1 and for k = 1 , · · · , K , define V σ,k = { u p σ ( i ) − u p σ ( i − } i =1 ,...,k ∪ V , where V is the base cone for ( P , ≺ , Ξ) as in Definition 3.5, and u p σ ( j ) , j = 0 , . . . , K is the represen-tation vector of p j ∈ P . For convenience, we define V σ, = V and we observe that V σ,K = V σ fromEquation (12). Theorem 3.8.
Algorithm 3 solves the LC-LEP.Proof.
Given ( P , ≺ , Ξ), we need to show that T ( P , ≺ , Ξ) = T alg ( P , ≺ , Ξ). We may assume thatcone( V ) is pointed since if not, then both T alg ( P , ≺ , Ξ) and T ( P , ≺ , Ξ) are empty.We first show that T alg ( P , ≺ , Ξ) ⊂ T ( P , ≺ , Ξ), i.e. for any σ ∈ T alg ( P , ≺ , Ξ) we show that theset Ξ σ (cid:54) = ∅ . As indicated above we assume cone( V ) = cone( V σ, ) is pointed. For Algorithm 3, lines2-4 returns the empty set if cone( V ) is not pointed. Otherwise, it passes to lines 5-9 which checksif σ part is a total order over { , . . . , K } . If so, it is added to the return variable, Ret . If σ part isnot a total order, then lines 10-17 extend it to a total order by recursively constructing V σ,i from V σ,i − for 1 ≤ i ≤ K .Therefore, it suffices to show that V σ,k are all pointed for k = 1 , . . . , K .12ix k ∈ { , . . . , K } . In lines 11-12, we find a candidate i ∈ { , . . . , K } which is not in the imageof σ part , and define V σ,k +1 = V σ,k ∪ { u (cid:48) } where u (cid:48) = u p i − u p σ ( k ) . In line 13, we verify that − u (cid:48) / ∈ V = V σ,k − and it follows from Proposition 3.4, that cone( V σ,k ) is pointed. For each σ appendedto Ret , we have cone( V σ,k ) is pointed for k = 0 , . . . , K . In particular, cone( V σ,K ) = cone( V σ ) ispointed, and from Proposition 3.6, we have Ξ σ (cid:54) = ∅ .We now prove that T ( P , ≺ , Ξ) ⊂ T alg ( P , ≺ , Ξ). Assume that σ ∈ T ( P , ≺ , Ξ). By definitionthis means that Ξ σ (cid:54) = ∅ and from Proposition 3.6, V σ is pointed. For each k = 1 , . . . , K , we have V σ,k ⊂ V σ , and thus V σ,k is pointed. As V σ,k is pointed, we know − ( u p σ ( k ) − u p σ ( k − ) / ∈ V σ,k − for k = 1 , . . . , K . Therefore, line 13 in Algorithm 3 will not fail to extend σ at each step in the recursionand after K recursive extensions, σ will be appended to Ret and thus, σ ∈ T alg ( P , ≺ , Ξ).
In this section we present a solution for the PSD problem described in Section 1. The solutionis based on the observation that the PSD problem has a natural Boolean lattice structure. Thus,for the linear PSD problem, the LC-LEP solver described in Section 3 provides a solution. Fornonlinear PSD problems, we construct a map that “embeds” it into an instance of LC-LEP (ofhigher dimension) in the sense that the inclusion in Equation (6) holds. We prove a sufficientcondition for which this inclusion is equality and describe a method for disqualifying spurioussolutions when it is strict.
Throughout this section (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) denotes a PSD problem for a fixed interaction function f of order type n ∈ N q as defined in Equation (5) where ≺ B is the Boolean lattice partial order.Our first goal is to show that (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) satisfies Equation (1), and in particular, that every σ ∈ T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) is a linear extension of a Boolean lattice. We start by defining appropriateindices for the elements of P . Definition 4.1.
Suppose n ∈ N q is the interaction type for f ∈ R [ z , . . . , z n ]. As in Defini-tion 1.1 let { I , . . . , I q } denote the indexing sets for each summand of f . Setting I := (cid:83) qj =1 I j we denote a typical element of I by I j ( k ) which we define as the k th largest element of I j . Let E := { α : { , . . . , n } → { , }} be the set of all Boolean functions defined on I . The Booleanindexing map , denoted by B : E → (cid:8) , . . . , n − (cid:9) , is defined by the formula B ( α ) := q (cid:88) j =1 n j − (cid:88) k =0 α ( I j ( k ))2 κ j,k κ j,k = k + j − (cid:88) j (cid:48) =1 n j (cid:48) . We will also consider, α ∈ E , as a vector of Boolean functions defined as follows. Let E j denotethe set of Boolean functions defined on I j . Then, elements of E can be represented as vectors ofthe form α = ( α , . . . , α q ) where α j := α (cid:12)(cid:12)(cid:12) I j ∈ E j for 1 ≤ j ≤ q. Note that under this identification, E has the equivalent representation as E = E × · · · × E q .13 efinition 4.2. Suppose n ∈ N q is the interaction type for an interaction function, f ∈ R [ z , . . . , z n ]as in Definition 1.1, and E denotes the corresponding Boolean indices. For α ∈ E , define p α ∈ P ⊂ R [ (cid:96) , . . . , (cid:96) n , δ , . . . , δ n ], by the formula p α := q (cid:89) j =1 (cid:88) k ∈ I j (cid:96) k + α ( k ) δ k . (13)When convenient, we use a linear indexing scheme for elements of P which we define via theBoolean indexing map by identifying p i := p α where α = B − ( i ). To avoid confusion, we exclusivelyuse Greek subscripts when referring to elements of P by their Boolean indices, and Latin subscriptswhen referring to elements of P by their linear indices. We leave it to the reader to check that thelinearly indexed polynomials in Examples 1.3 and 1.4 are in agreement with that of Definition 4.2via this identification. Definition 4.3.
Let α, β ∈ E be a pair of Boolean indices corresponding to n ∈ N q . An orderedpair ( α, β ) satisfies the one bit condition if α ( I j ( k )) ≤ β ( I j ( k )), for all 1 ≤ j ≤ q and 0 ≤ k ≤ n j − exactly one ( j, k ) pair. Remark 4.4.
Observe that if ( α, β ) satisfy the one bit condition and ( j , k ) is the unique pair forwhich α and β take different values, then α ( I j ( k )) = 0 and β ( I j ( k )) = 1. Remark 4.5.
The one bit condition induces a poset structure on E by setting α ≺ β for each( α, β ) satisfying the one bit condition, and extending the relation transitively. Definition 4.6.
Suppose n ∈ N q is the interaction type for an interaction function, f ∈ R [ z , . . . , z n ]as in Definition 1.1, and E denotes the corresponding Boolean indices. Let P be the set of poly-nomials indexed as in Definition 4.2. The PSD problem is defined by the triple, ( P , ≺ , (0 , ∞ ) n )where ≺ is given by Remark 4.5.The next proposition proves that ( P , ≺ , (0 , ∞ ) n ) satisfies Equation (1), and furthermore that( P , ≺ ) is a Boolean partial order which justifies expressing the PSD problem as (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) . Proposition 4.7.
Consider a PSD problem ( P , ≺ , (0 , ∞ ) n ) . Then,1. ( P , ≺ ) is a Boolean lattice.2. For any α, β ∈ E , if α ≺ β , then p α ( ξ ) < p β ( ξ ) for all ξ ∈ (0 , ∞ ) n . Proof.
To prove the first claim, let S n = { , . . . , n } and let (2 S n , ≺ B ) denote the standard Booleanlattice. Define a map, ϕ : E → S n , by the formula ϕ ( α ) = { j ∈ S n : α ( j ) = 1 } , and we note that ϕ is a bijection since E is defined to be the collection of all Boolean maps definedon S n . Furthermore, for any α, β ∈ E , we have by Definition 4.3 that α ≺ β if and only if { j ∈ S n : α ( j ) = 1 } ⊂ { j ∈ S n : β ( j ) = 1 } ϕ is an order isomorphism.To establish the second claim, we must show that if α ≺ β , then p α ( ξ ) < p β ( ξ ) holds for all ξ ∈ (0 , ∞ ) n . Note that by transitivity, it suffices to prove this holds for ( α, β ) satisfying the onebit condition. In this case we have β ( I j ( k )) − α ( I j ( k )) = (cid:40) j = j and k = k . for some j ∈ { , . . . , q } , k ∈ { , . . . , n j − } . If ξ = ( (cid:96) , . . . , (cid:96) n , δ , . . . , δ n ) ∈ Ξ, then from Equation(13) we have p β ( ξ ) = (cid:88) k ∈ I j (cid:96) k + α ( I j ( k )) δ k + δ k (cid:89) j (cid:54) = j (cid:88) k ∈ I j (cid:96) k + α ( I j ( k )) δ k = p α ( ξ ) + δ k (cid:89) j (cid:54) = j (cid:88) k ∈ I j (cid:96) k + α ( I j ( k )) δ k > p α ( ξ )as required.With Proposition 4.7 in mind, we return to writing ≺ B in place of ≺ for the PSD problem where ≺ B is the partial order of a Boolean lattice inherited by P from the one bit condition. We consider two cases of the PSD problem: the interaction type n ∈ N q for the function f ∈ R [ z , . . . , z n ] has the form n = (1 , , . . . ,
1) or n = ( n ). In the first case, f is linear (see Example 1.4)and the PSD problem is an instance of LC-LEP. In the second case, log f is linear so after asimple change of variables, we obtain an instance of LC-LEP with equivalent solutions since log ismonotone, hence order preserving(cid:32)l. We focus on the first case, leaving it to the reader to checkthat the second case is same modulo the evaluation domain ( R n versus (0 , ∞ ) n ). Following thealgorithm described in Section 3, we encode the partial order defined by ≺ B as a set of linearconstraints defined by a base cone which we must show is pointed. We begin by denoting the setof representation vectors for P as V := (cid:110) u p α ∈ { , } n : u p α is the representation vector of p α , α ∈ E (cid:111) . We define the set, V ≺ B := (cid:8) u p β − u p α : u p α , u p β ∈ V , ( α, β ) satisfies the one bit condition (cid:9) . (14)which encodes the ≺ B partial order into the representation vectors. These vectors will be thegenerators of the base cone for the algorithm in Section 3. Thus, we must show that cone( V ≺ B )generates a pointed cone. Lemma 4.8.
Let C := cone( V ≺ B ) denote the cone generated by V ≺ B , then C is pointed. roof. By Proposition 3.2, C is closed and convex so it suffices to prove that if v ∈ C and − v ∈ C ,then v = 0. Fix ξ ∈ (0 , ∞ ) n and suppose ( α, β ) satisfies the one bit condition. By the formula inEquation (13) it follows that p β − p α = δ i for some i ∈ { , . . . , n } . Since δ i = ξ n + i > ξ ∈ Ξ, it follows that p β ( ξ ) − p α ( ξ ) > , for every ( α, β ) satisfying the one bit condition. Passing to the representation vectors, it followsthat for every v ∈ V ≺ B , we have v · ξ >
0. Taking the conic hull, we have that if v ∈ C \ { } ,then v · ξ >
0. It follows that if v, − v ∈ C simultaneously, then v · ξ ≥ − v · ξ ≥ v = 0. Given a general PSD problem (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) we present the construction of a LC-LEP denotedby ( P (cid:48) , ≺ B , R m ) with the property that T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) ⊆ T ( P (cid:48) , ≺ B , R m ). The importance ofthis is that ( P (cid:48) , ≺ B , R m ) can be solved using Algorithm 3 and hence we obtain a rigorous upperbound on T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) . Definition 4.9.
Given an interaction type n ∈ N q , let E = E × · · · × E q denote the correspondingBoolean indices. Set m := (cid:80) qj =1 n j and define the linearized evaluation domain to be R n × · · · × R nq ∼ = R m . (15)Define a polynomial ring in m indeterminates with Boolean indexing by R := R (cid:2)(cid:8) x j,α j : α j ∈ E j , ≤ j ≤ q (cid:9)(cid:3) , (16)and define a collection of linear polynomials by P (cid:48) := { p (cid:48) α : α ∈ E } ⊂ R where p (cid:48) α := q (cid:88) j =1 x j,α j . The linearized PSD problem determined by n is to compute T ( P (cid:48) , ≺ B , R m ). Theorem 4.10.
Fix an interaction type, n ∈ N q and let T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) and T ( P (cid:48) , ≺ B , R m ) denote the corresponding PSD and linearized PSD problems, respectively. The following are true.1. Let α, β ∈ E and ξ ∈ (0 , ∞ ) n . If p α ( ξ ) < p β ( ξ ) and ξ (cid:48) j,α j = log (cid:88) k ∈ I j ξ k + α j ( k ) ξ n + k ∈ R m , then p (cid:48) α ( ξ (cid:48) ) < p (cid:48) β ( ξ (cid:48) ) .2. T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) ⊆ T ( P (cid:48) , ≺ B , R m ) . roof. To prove the first claim, we define a map, T : (0 , ∞ ) n → R m , by ξ (cid:55)→ ξ (cid:48) := T ( ξ ) where thecoordinates of ξ (cid:48) are given by the formula ξ (cid:48) j,α j = log (cid:88) k ∈ I j ξ k + α j ( k ) ξ n + k . (17)Observe that T is defined to satisfy the functional equationlog ◦ p α ( ξ ) = p (cid:48) α ◦ T ( ξ ) for all α ∈ E, ξ ∈ (0 , ∞ ) n . (18)Therefore, if α, β ∈ E and ξ ∈ (0 , ∞ ) n satisfies p α ( ξ ) < p β ( ξ ), then log ( p α ( ξ )) < log ( p β ( ξ )) andit follows from Equation (18) that p (cid:48) α ( ξ (cid:48) ) < p (cid:48) β ( ξ (cid:48) ) where ξ (cid:48) = T ( ξ ) as required.To prove the second claim, consider P and P (cid:48) equipped with the linear indices as in Definition4.2, and suppose σ ∈ T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) . Then, by definition there exists ξ ∈ (0 , ∞ ) n satisfying p σ (0) ( ξ ) < p σ (1) ( ξ ) < · · · < p σ (2 n − ( ξ ) . Let ξ (cid:48) = T ( ξ ) and apply the first result to successive pairs in the ordering which implies that forall 0 ≤ k ≤ n −
2, we have p (cid:48) σ ( k ) ( ξ (cid:48) ) = log (cid:0) p σ ( k ) ( ξ ) (cid:1) < log (cid:0) p σ ( k +1) ( ξ ) (cid:1) = p (cid:48) σ ( k +1) ( ξ (cid:48) ) . Thus, we ξ (cid:48) ∈ R m satisfies p (cid:48) σ (0) ( ξ (cid:48) ) < p (cid:48) σ (1) ( ξ (cid:48) ) < · · · < p (cid:48) σ (2 n − ( ξ (cid:48) ) , and it follows that σ ∈ T ( P (cid:48) , ≺ B , R m ) which completes the proof. Example 4.11.
Recall the PSD in Example 1.3 with interaction function f ( z ) = ( z + z ) z corresponding to interaction type n = (2 , p (cid:48) = x , + x , p (cid:48) = x , + x , p (cid:48) = x , + x , p (cid:48) = x , + x , p (cid:48) = x , + x , p (cid:48) = x , + x , p (cid:48) = x , + x , p (cid:48) = x , + x , where we have used linear indexing to match the polynomials in Example 1.3. = (2 , , . . . , . In this section we prove the following theorem.
Theorem 4.12.
Let f be an interaction function with interaction type, n = (2 , , . . . , . Let T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) denote the corresponding PSD problem and ( P (cid:48) , ≺ B , R m ) the associated lin-earized PSD problem. Then T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) = T ( P (cid:48) , ≺ B , Ξ (cid:48) ) where Ξ (cid:48) = R m ∩ {− ξ (cid:48) , + ξ (cid:48) , + ξ (cid:48) , − ξ (cid:48) , > } . The proof of the theorem is based on the following lemma17 emma 4.13.
Fix parameters, x , x , x , x ∈ R , and define the function, g : R → R by the formula g ( t ) = exp( tx ) − exp( tx ) − exp( tx ) + exp( tx ) . If x < x ≤ x < x , then g has a positive root if and only if g (cid:48) (0) < .Proof. Suppose first that t is a root of g . Expanding exp( t x ) and exp( t x ) to first order about x and x , respectively, and applying the mean value theorem yields the formula g ( t ) = − t exp( t c )( x − x ) − t exp( t c )( x − x ) = 0 (19)for some c ∈ ( x , x ) and c ∈ ( x , x ). We define k = c − c and multiply Equation (19) by t e − kt to obtain e kt ( x − x ) − ( x − x ) = 0 . Noting that c < x < c , it follows that k >
0. Therefore if t >
0, then x − x < x − x orequivalently, g (cid:48) (0) = x − x − x + x < g (cid:48) (0) < g has at least one positive root since clearly g (0) = 0 and lim t →∞ g ( t ) = ∞ . Proof of Theorem 4.12.
Suppose σ ∈ T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) and ξ ∈ (0 , ∞ ) nσ , then by Theorem 4.10we have ξ (cid:48) = T ( ξ ) ∈ R mσ . Note that by definition the first four coordinates of ξ (cid:48) are given by theformulas ξ (cid:48) , = log( ξ + ξ ) ξ (cid:48) , = log( ξ + ξ + ξ n +2 ) ξ (cid:48) , = log( ξ + ξ + ξ n +1 ) ξ (cid:48) , = log( ξ + ξ + ξ n +1 + ξ n +2 ) . Since ξ i > i ∈ { , , n + 1 , n + 2 } , it follows that − ξ (cid:48) , + ξ (cid:48) , + ξ (cid:48) , − ξ (cid:48) , > , so we have σ ∈ T ( P (cid:48) , ≺ B , Ξ (cid:48) ).Conversely, suppose σ ∈ T ( P (cid:48) , ≺ B , Ξ (cid:48) ) and ξ (cid:48) ∈ Ξ (cid:48) σ . From the Boolean lattice ≺ B we have ξ (cid:48) , <ξ (cid:48) , ≤ ξ (cid:48) , < ξ (cid:48) , or ξ (cid:48) , < ξ (cid:48) , ≤ ξ (cid:48) , < ξ (cid:48) , . Moreover, ξ (cid:48) also satisfies, − ξ (cid:48) , + ξ (cid:48) , + ξ (cid:48) , − ξ (cid:48) , > t (cid:48) > ξ (cid:48) := t (cid:48) ξ (cid:48) satisfiesexp( ˆ ξ (cid:48) , ) − exp( ˆ ξ (cid:48) , ) − exp( ˆ ξ (cid:48) , ) + exp( ˆ ξ (cid:48) , ) = 0 . Next, we define ˆ ξ ∈ (0 , ∞ ) n byˆ ξ j = exp( ˆ ξ (cid:48) j, ) 2 < j ≤ n exp( ˆ ξ (cid:48) j, ) − exp( ˆ ξ (cid:48) j, ) n + 2 < j < n exp( ˆ ξ (cid:48) , ) j = 1 , ξ (cid:48) , ) − exp( ˆ ξ (cid:48) , ) j = n + 1exp( ˆ ξ (cid:48) , ) − exp( ˆ ξ (cid:48) , ) j = n + 2One easily verifies that ˆ ξ j > ≤ j ≤ n , and that T ( ˆ ξ ) = ˆ ξ (cid:48) . From Theorem 4.10, wehave ˆ ξ ∈ (0 , ∞ ) nσ = { ξ ∈ (0 , ∞ ) n : p σ (0) ( ξ ) < p σ (1) ( ξ ) < · · · < p σ (2 n − ( ξ ) } which implies that σ ∈ T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) . 18 .5 Solving the general PSD problem In the general case, we have T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) (cid:40) T ( P (cid:48) , ≺ B , R m ), and thus, computing T ( P (cid:48) , ≺ B , R m )provides only a set of candidates for T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) . This candidate set contains spurious lin-ear extensions so we consider the problem of removing linear extensions which are non-admissible.We have two strategies for doing this efficiently.The first is to restrict the evaluation domain to a strict subset, Ξ (cid:48) (cid:40) R m , such that we still havethe inclusion T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) ⊆ T ( P (cid:48) , ≺ B , Ξ (cid:48) ) . (20)Restricting to a smaller evaluation domain amounts to imposing more of the algebraic constraintsa-priori which results in improved efficiency. In order for the candidate set on the right hand side tobe efficiently computable using the algorithm in Section 3, it must be an instance of LC-LEP i.e. Ξ (cid:48) should be a polytope. For example, for the PSD with interaction type n = (2 , , . . . , (cid:48) = R m ∩ (cid:8) ξ (cid:48) ∈ R m : − ξ (cid:48) , + ξ (cid:48) , + ξ (cid:48) , − ξ (cid:48) , > (cid:9) . In terms of the algorithm in Section 3, this domain restriction amounts to taking our base cone inAlgorithm 3 to be cone( V ) where V = V ≺ B ∪ { u } and u is the representation vector for the linear functional defined by the formula x (cid:55)→ − x , + x , + x , − x , . The requirement that this linear functional must be strictly positive is a special case of the followingLemma whose proof is a trivial computation.
Lemma 4.14.
Suppose α, α (cid:48) , β, β (cid:48) are Boolean indices such that for any ξ ∈ (0 , ∞ ) n , the followingequations are satisfied. p α ( ξ ) < p β ( ξ ) < p β (cid:48) ( ξ ) < p α (cid:48) ( ξ ) p α ( ξ ) + p α (cid:48) ( ξ ) = p β ( ξ ) + p β (cid:48) ( ξ ) . Then, log( p α ( ξ )) + log( p α (cid:48) ( ξ )) − log( p β ( ξ )) − log( p β (cid:48) ( ξ )) > . Lemma 4.14 provides a means to restrict the evaluation domain for the general linearized PSDproblem as follows. Fix j ∈ { , . . . , q } and suppose { α, α (cid:48) , β, β (cid:48) } ⊂ E differ only in the j th coordinatewith α ≺ B β ≺ B β (cid:48) ≺ B α (cid:48) , and also assume that B ( α ) + B ( α (cid:48) ) = B ( β ) + B ( β ) (cid:48) where B is theBoolean indexing map. Then, it follows that for any ξ ∈ Ξ, the values, { p α ( ξ ) , p α (cid:48) ( ξ ) , p β ( ξ ) , p β (cid:48) ( ξ ) } ,satisfy both equations in Lemma 4.14. Therefore, if u ( { α, α (cid:48) , β, β (cid:48) } ) is the representation vector forthe linear functional defined by x (cid:55)→ x j,B ( β ) + x j,B ( β (cid:48) ) − x j,B ( α ) − x j,B ( α (cid:48) ) , then v ( { α, α (cid:48) , β, β (cid:48) } ) lies in V σ for any σ ∈ T ( P , ≺ B , (0 , ∞ ) n . Equivalently, we may impose therequired linear constraint, x j,B ( β ) + x j,B ( β (cid:48) ) − x j,B ( α ) − x j,B ( α (cid:48) ) > ≤ j ≤ q , we define V j := { u ( { α, α (cid:48) , β, β (cid:48) } ) : B ( α ) + B ( α (cid:48) ) = B ( β ) + B ( β ) (cid:48) , α ≺ B β ≺ B β (cid:48) ≺ B α (cid:48) } V = V ≺ B ∪ V Ξ where V Ξ = q (cid:91) j =1 V j . Applying Algorithm 3 with the base cone generated by V is equivalent to solving the instance ofLC-LEP defined by ( P (cid:48) , ≺ B , Ξ (cid:48) ) where Ξ (cid:48) is the restriction of R m to the subset for which the linearfunctionals defined by each v ∈ V j are strictly positive for each 1 ≤ j ≤ q .In addition to restricting the computation to the polytopes discussed above, we can reuse so-lutions of smaller PSD problems in some larger computations. As an example, suppose P (cid:48) = { p (cid:48) , . . . , p (cid:48) } is the set of interaction polynomials for the PSD with interaction type n (cid:48) = (2 ,
1) and P := { p , . . . , p } the polynomials for the PSD problem with interaction type n = (2 , , P (cid:48) induces an imposed linear order on the even indexedpolynomials, P even := { p , p , . . . , p } ⊂ P . A similar linear order is induced on the odd indexedpolynomials, P odd := { p , p , . . . , p } ⊂ P . Hence, a necessary condition to have an admissiblelinear extension for P is that the order of P even and P odd must both be consistent with one of thePSD solutions in T (cid:0) P (cid:48) , ≺ B , (0 , ∞ ) (cid:1) . This implies the inclusion T (cid:0) P , ≺ B , (0 , ∞ ) (cid:1) ⊆ (cid:91) σ (cid:48) ∈T ( P (cid:48) , ≺ B , (0 , ∞ ) ) T ( P , ≺ B ∪ ≺ σ (cid:48) , (0 , ∞ ) ) (21)where ≺ B ∪ ≺ σ (cid:48) represents the refinement of the Boolean lattice partial order, and the partial orderinduced by σ (cid:48) on the even/odd subsets.To exploit this in general, we say that the PSD problem of type n (cid:48) is a sub-problem for the PSDproblem of type n whenever the polynomials for n must obey an implied partial order determinedby the solutions of n (cid:48) . Notice that the preceding discussion as well as Equation (21) applies alsoto an arbitrary polytope. Therefore, if there are a total of k admissible linear extensions for allsub-problems of the PSD problem of type n which we have previously computed, then we bootstrapthose results when computing T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) via the inclusion T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) ⊆ k (cid:91) i =1 T ( P , ≺ B ∪ ≺ σ (cid:48) i , Ξ (cid:48) ) ⊆ T ( P , ≺ B , Ξ (cid:48) )where ≺ B ∪ ≺ σ (cid:48) i represents the refinement of the Boolean lattice partial order, and the partialorder induced by σ (cid:48) i on the corresponding subsets obtained from any sub-problem. This techniquehas been used in the computation for all the cases of order ≥
4. Observe that the computation of T ( P , ≺ B ∪ ≺ σ (cid:48) i , Ξ (cid:48) ) can be done distributively for i = 1 , . . . , k on different computational nodes,which, as is indicated in Section 5, we employed for the PSD problems of orders 5 and 6.In the special case of Section 4.4, we proved that inclusion in Equation (20) is actually equalitywhen Ξ (cid:48) is constructed as we have described. However, in the typical case, these additional alge-braic constraints are not sufficient to remove all spurious linear extensions except in the case n =(2 , , . . . , (cid:48) such that T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) = T ( P (cid:48) , ≺ B , Ξ (cid:48) ) for other interaction types. However, in the remainder of this section we considerthe problem of extracting T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) from T ( P (cid:48) , ≺ B , Ξ (cid:48) ) when they are not equal.Observe that we may obtain large subsets of T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) simply by sampling. Theparticular strategy that we adopted is as follows. We uniformly sampled between 10 and 10 points ξ = ( l , . . . , l n , δ , . . . , δ n ) ∈ Z n + ∩ B n ∞ ( r ) , B n ∞ ( r ) = {(cid:107) ξ (cid:107) ∞ ≤ r } . We chose r = 1000. Mathematically the particular choice of r is notimportant since the PSD polynomials are homogeneous, though in practice it does have an effect onsampling precision and speed. For each such ξ we evaluated { p α ( ξ ) : p ∈ P} . If σ ∈ S n denotes thelinear order of these values, then ξ serves as a “witness” for the claim that Ξ σ (cid:54) = ∅ . This produces S (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) := (cid:8) σ ∈ T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) : σ is witnessed by at least one sample (cid:9) . Obviously, S (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) ⊆ T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) ⊆ T ( P , ≺ B , Ξ (cid:48) ) . In general, sampling is relatively efficient and in cases where T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) is not too large(see Table 1 for details), we recover the entire solution.Once we constructed the set S (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) , T ( P , ≺ B , Ξ (cid:48) ) from sampling and algorithmsin section 3 respectively, we can apply CAD algorithm to check the set { ξ ∈ Ξ : p σ (0) ( ξ ) < p σ (1) ( ξ ) < · · · < p σ (2 n − ( ξ ) } is empty or not for each σ ∈ T ( P , ≺ B , Ξ (cid:48) ) \ S (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) and then T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) isrecovered. The CAD algorithm implementation we are using is CylindricalAlgebraicDecompositionin Mathematica 11 [12]. In this section we provide (see Table 1) the results of our computations for interaction functionsof orders 4, 5, and 6. A slightly different approach was taken to compute orders 5 and 6, fromthat used for 4. This had to do with the machines being used, but highlights the flexibility of ourmethod.For interaction functions of order 4, we applied Algorithm 1 using a rational linear programmingalgorithm. In particular, we used the implementation MixedIntegerLinearProgram from SageMath 8[20]. This implies that the output of Algorithm 3 is correct. Observe that interaction type (1 , , , T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) . The factthat our output agrees with that of [13] suggests that our code is functioning as desired. To computethe interaction type (2 , ,
1) we apply Algorithm 3 to obtain T ( P (cid:48) , ≺ B , R m ). By Theorem 4.12 thisdetermines T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) .To solve the PSD problem from interaction types (2 ,
2) and (3 ,
1) requires that we make useof the strategy discussed in Section 4.5. Again, we use Algorithm 3 to obtain T ( P (cid:48) , ≺ B , R m ). ByTheorem 4.10, T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) ⊂ T ( P (cid:48) , ≺ B , R m ). As indicated in Column 7 of Table 1, wechose 10 samples from (0 , ∞ ) and identified 5344 and 3084 linear orders, respectively. We ranCylindricalAlgebraicDecomposition in Mathematica 11 [12] on each element of T ( P (cid:48) , ≺ B , R m ) \S (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) . As can be seen by comparing Columns 6 and 3, none of these elements wereadmissible.We now turn to the computations of interaction functions of order 5 and 6. As these problemsare too big to be done on a laptop we turned to a server for which SageMath was not installed.Thus, we made use of a numerical linear programing algorithm, linprog from Python 3.5 packagescipy [23] with the default numerical error 10 − , in Algorithm 1. The interaction type (1 , , , , , , , , ,
1) are linear and type (5) and (6) are log linear, and therefore via Algorithm 3 we21 T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) T ( P (cid:48) , ≺ B , Ξ (cid:48) ) T ( P (cid:48) , ≺ B , R m ) S ( P (cid:48) , ≺ B , (0 , ∞ ) n )(1,1,1,1) 336 - - -(4) 336 - - -(2,1,1) 1,344 1,344 2,352 -(2,2) 5,344 7,920 26,640 5,344(3,1) 3,084 5,112 68,641 3,084(1,1,1,1,1) 61,920 - - 61,920(5) 61,920 - - 61,920(2,1,1,1) 790,200 790,200 * 790,200(2,2,1) - 1,1035,808 * 6570952(3,2) - * * 71,959,088 † (4,1) - * * 11,213,616 † (1,1,1,1,1,1) 89,414,640 - - 89,414,640(6) 89,414,640 - - 89,414,640Table 1: Computational results for several PSD problems. Column 1 indicates the interaction type.Column 2 provides the number of elements in the AC-LEP of interest. Column 3 provides thenumber of elements in an associated LC-LEP. This is not relevant where the AC-LEP problem ofinterest is a LC-LEP problem and is indicated by -. The * indicates that the computation wastoo large to complete. Column 4 provides the number of elements in the linearized PSD problemwithout additional constraints. Again the irrelevance for linear problems is indicated by - and *indicates that the computation is large to be performed. The last column indicates the numberof cells identified via sampling. We used 10 samples for all n = 4 cases and 10 samples for the n = 5 , † indicates that our sampling was not sufficient.obtain T alg (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) . We use the sampling technique (see Columns 7 and 8 of Table 1) toverify each of the elements of T alg (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) , thereby obtaining T (cid:0) P , ≺ B , (0 , ∞ ) n (cid:1) .The computation for each order 4 case was done on a Mac Pro laptop (2.7 Hz Intel i5 andmemory 8GB) with computation time under 4 hours. The computation of the remaining cases weredone using a computing server with CentOs, intel 17.1, memory 32 GB, and less than 30 nodes.The computation time for both (1 , , , ,
1) and (5) was less than 4 hours, while the computationtime for (2 , , ,
1) was on the order of 7 days. The codes which produced all of the computationsin Table 1 are available on GitHub.
Acknowledgments
The authors would like to thank Shaun Harker, Sandra Di Rocco, Tomas Gedeon, and Mike Saks forhelpful conversations. The authors acknowledge support from NSF DMS-1521771, DMS-1622401,DMS-1839294, the NSF HDR TRIPODS award CCF-1934924, DARPA contracts HR0011-16-2-0033and FA8750-17-C- 0054, and NIH grant R01 GM126555-01.22 eferences [1] David L Applegate, William Cook, Sanjeeb Dash, and Daniel G Espinoza. Exact solutions tolinear programming problems.
Operations Research Letters , 35(6):693–699, 2007.[2] Saugata Basu, Richard Pollack, and Marie Roy. Algorithms in Real Algebraic Geometry. pages1–676.[3] Stephen Boyd and Lieven Vandenberghe.
Convex Optimization , 2004.[4] Graham Brightwell and Peter Winkler. Counting Linear Extensions is
Pro-ceedings of the Twenty-third Annual ACM Symposium on Theory of Computing , STOC ’91,pages 175–181, New York, NY, USA, 1991. ACM.[5] Graham R Brightwell. The Number of Linear Extensions of Ranked Posets. Technical report,2003.[6] Christopher W. Brown. Improved projection for cylindrical algebraic decomposition.
Journalof Symbolic Computation , 32(5):447–465, 2001.[7] Christopher W. Brown. Constructing a single open cell in a cylindrical algebraic decomposition.
Proceedings of the International Symposium on Symbolic and Algebraic Computation, ISSAC ,pages 133–140, 2013.[8] Christopher W. Brown and James H. Davenport. The complexity of quantifier elimination andcylindrical algebraic decomposition.
Proceedings of the International Symposium on Symbolicand Algebraic Computation, ISSAC , pages 54–60, 2007.[9] George E. Collins. Quantifier elimination for real closed fields by cylindrical algebraic decom-postion. In H. Brakhage, editor,
Automata Theory and Formal Languages , pages 134–183,Berlin, Heidelberg, 1975. Springer Berlin Heidelberg.[10] Bree Cummins, Tomas Gedeon, Shaun Harker, Konstantin Mischaikow, and Kafung Mok.Combinatorial Representation of Parameter Space for Switching Networks.
SIAM Journal onApplied Dynamical Systems , 15(4):2176–2212, 2016.[11] T. Fine and J. Gill. The enumeration of comparative probability relations.
Ann. Prob. , 4:667–673, 1976.[12] Wolfram Research, Inc. Mathematica, Version 11.[13] D. Maclagan. Boolean Term Orders and the Root System B n . Order , 15:279–295, 1999.[14] Bree Cummins Marcio Gameiro, Shaun Harker. Dsgrn: Dynamic signatures generated byregulatory networks. https://github.com/marciogameiro/DSGRN , 2020.[15] Scott McCallum. An improved projection operation for cylindrical algebraic decomposition.
Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligenceand Lecture Notes in Bioinformatics) , 204 LNCS:277–278, 1985.[16] OEIS Foundation Inc.
The On-Line Encyclopedia of Integer Sequences , 2020. . 2317] Jichang Sha and D. J. Kleitman. The number of linear extensions of subset ordering.
DiscreteMathematics , 63(2-3):271–278, 1987.[18] Richard P. Stanley. An introduction to hyperplane arrangements. In
Geometric combinatorics ,volume 13 of
IAS/Park City Math. Ser. , pages 389–496. Amer. Math. Soc., Providence, RI,2007.[19] Herry Suprajitno and I Bin Mohd. Linear programming with interval arithmetic.
Int. J.Contemp. Math. Sciences , 5(7):323–332, 2010.[20] The Sage Developers.
SageMath, the Sage Mathematics Software System (Version 8) . .[21] Seinosuke Toda. PP is as hard as the polynomial-time hierarchy. SIAM Journal on Computing ,20(5):865–877, 1991.[22] Robert J. Vanderbei.
Linear Programming Foundations and Extensions , Third Edition.[23] Pauli Virtanen et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.