\UseRawInputEncoding

Computing linear extensions for polynomial posets subject to algebraic constraints

Shane Kepley [email protected] Lun Zhang [email protected] and Konstantin Mischaikow [email protected]

Abstract

In this paper we consider the classical problem of computing linear extensions of a given poset which is well known to be a difficult problem. However, in our setting the elements of the poset are multivariate polynomials, and only a small “admissible” subset of these linear extensions, determined implicitly by the evaluation map, are of interest. This seemingly novel problem arises in the study of global dynamics of gene regulatory networks in which case the poset is a Boolean lattice. We provide an algorithm for solving this problem using linear programming for arbitrary partial orders of linear polynomials. This algorithm exploits this additional algebraic structure inherited from the polynomials to efficiently compute the admissible linear extensions. The biologically relevant problem involves multilinear polynomials and we provide a construction for embedding it into an instance of the linear problem.

1 Introduction

Consider a set of real polynomials $\mathcal{P}$ , defined on a domain $\Xi\subset\mathbb{R}^{d}$ , equipped with a partial order $\prec$ . We are interested in identifying linear extensions (total orders compatible with $\prec$ ) that are satisfied by $\mathcal{P}$ under evaluation at a point in $\Xi$ . To be more precise consider a semi-algebraic set $\Xi\subset\mathbb{R}^{d}$ , called the evaluation domain, and a collection of polynomials $\mathcal{P}:=\left\{p_{0},\ldots,p_{K}\right\}\subset\mathbb{R}[x_{1},\ldots,x_{d}]$ . Let $\prec$ denote a partial order on $\mathcal{P}$ such that if $p\prec q$ , then

p(\xi)<q(\xi)\quad\text{for all}\ \xi\in\Xi.

(1)

Let $S_{K+1}$ denote the set of permutations on $K+1$ symbols. We identify linear extensions of $\mathcal{P}$ with a subset of $S_{K+1}$ as follows. Given $\sigma\in S_{K+1}$ , let $\prec_{\sigma}$ denote the linear order

p_{\sigma(0)}\prec_{\sigma}p_{\sigma(1)}\prec_{\sigma}\cdots\prec_{\sigma}p_{\sigma(K)}.

We define the realizable set associated to $\sigma$ by

\Xi_{\sigma}:=\left\{\xi\in\Xi:p_{\sigma(k)}(\xi)<p_{\sigma(k+1)}(\xi)\ \text{for all}\ 0\leq k\leq K-1\right\}.

(2)

Observe that if $\Xi_{\sigma}\neq\emptyset$ , then $\prec_{\sigma}$ is a linear extension of $\prec$ . The algebraically constrained linear extension problem (AC-LEP) defined by $(\mathcal{P},\prec,\Xi)$ is to determine

\displaystyle\mathcal{T}(\mathcal{P},\prec,\Xi):=\left\{\sigma\in S_{K+1}:\Xi_{\sigma}\ \text{is nonempty}\right\}.

Notice that in the formulation of the AC-LEP we have identified the partial order $\prec_{\sigma}$ with the associated element in $S_{K+1}$ . We will use this identification throughout the remainder of the paper. We say that each $\sigma\in\mathcal{T}(\mathcal{P},\prec,\Xi)$ is an admissible linear extension.

As is discussed in Section 2 our immediate motivation for studying the AC-LEP comes from modeling the dynamics of regulatory networks in biology and in particular characterizing relevant subsets of the parameter space. For the moment we attempt to put this problem into a broader mathematical context as the problem itself, as well as our solutions for some special cases, have elements of both classical real algebraic geometry and order theory.

Quantifier elimination and real algebraic geometry

Observe that if $\Xi=\mathbb{R}^{d}$ , and $\mathcal{P}$ is an arbitrary collection of polynomials, then $\sigma\in S_{K+1}$ is admissible if and only if there exists $\xi\in\mathbb{R}^{d}$ such that $p_{\sigma(k)}(\xi)-p_{\sigma(k+1)}(\xi)<0$ for all $0\leq k\leq K$ . These inequalities define a semi-algebraic set. Therefore, if $\prec$ is the trivial partial order (i.e. $\mathcal{P}$ is a single anti-chain), then this instance of AC-LEP is equivalent to the classical problem of decidability for real semi-algebraic sets.

The previous example illustrates a major challenge in solving the AC-LEP. The first general algorithm for solving the quantifier elimination/decidability problem for polynomials in $\mathbb{R}^{d}$ with feasible running time was the cylindrical algebraic decomposition (CAD) given by Collins [13] in 1975. The CAD algorithm works by subdividing $\Xi$ into subsets on which the polynomials are sign invariant. Given such a decomposition, decidability is reduced to simply evaluating each polynomial at a sample point located in each subset and checking if it satisfies the necessary inequalities. Unfortunately, the computational complexity of the algorithm grows like

\mathcal{O}\left(\left(2D\right)^{2^{d}-1}(K+1)^{2^{d}-1}2^{2^{d-1}-1}\right)\quad\text{where}\quad D=\max\left\{\deg p:p\in\mathcal{P}\right\}.

(3)

This worst case running time is known to be sharp even for classes of “nice” polynomials e.g. linear [10], and moreover, the worst case is also typical [4]. As a result, the question of whether or not $\sigma\in\mathcal{T}\left(\mathcal{P},\prec,\Xi\right)$ , even for a single $\sigma\in S_{K+1}$ , is often intractable for problems of practical interest.

In addition, the CAD 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, we note that if additional algebraic constraints are added e.g. we assume $\Xi_{0}\subset\Xi\subset\mathbb{R}^{d}$ is a strict semi-algebraic subset, then the CAD algorithm can handle this by simply appending the polynomial constraints which define $\Xi_{0}$ to the set of polynomials. However, this dramatically increases the complexity of the CAD algorithm, despite the fact that the number of admissible linear extensions can only decrease.

Some improved algorithms have been proposed which aim to reduce the complexity of specific aspects of the problem or for special classes of polynomials (e.g. [35, 8, 9]). These improvements often provide dramatic algorithmic speedups for checking whether $\sigma\in\mathcal{T}\left(\mathcal{P},\prec,\Xi\right)$ for a single linear extension. However, these algorithms have the same worse case running time as the general CAD algorithm and understanding which classes of polynomials benefit is still a very active area of research. Therefore, these improved algorithms alone are not sufficient to handle instances of AC-LEP since we are interested in determining which of the $(K+1)!$ possible semi-algebraic sets are nonempty. An efficient algorithm that does not produce a decomposition of $\Xi$ into sign invariant subsets would still need to be called $(K+1)!$ times. However, we will make use of these improved algorithms 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 $\mathcal{P}$ is a collection of polynomials. Hence, we focus only on the poset structure $(\mathcal{P},\prec)$ , and consider the problem of computing all linear extensions. The related problem of counting all linear extensions of a partial order is a well studied problem in order theory. Its importance is due in large part to its connection with the complexity of sorting elements in a list. If one considers a list of $(K+1)$ distinct values which have been partially sorted by making pairwise comparisons on a subset of its elements, then these comparisons induce a partial order. Therefore, the linearly ordered values of the fully sorted list are given by one of the possible linear extensions for the partial order. As a result, the complexity of completely sorting a list is intimately connected to counting linear extensions for posets.

Observe that computing the set of all linear extensions of $(\mathcal{P},\prec)$ is not easier than counting them which is known to be $\#P$ -complete [6]. In particular, a polynomial time algorithm for computing all possible linear extensions for arbitrary posets would imply that $P=NP$ by Toda’s theorem [51]. Moreover, we are interested not only in counting linear extensions, but explicitly computing them. Therefore, we are also concerned with how fast the number of admissible linear extensions grows.

For reasons we discuss in Section 4, we are specifically interested in the case that $(\mathcal{P},\prec)$ is a Boolean lattice. Specifically, for fixed $n\in\mathbb{N}$ , define $S_{n}:=\left\{1,\dots,n\right\}$ and let $2^{S_{n}}$ denote its power set. The standard $n$ -dimensional Boolean lattice is the poset, $(2^{S_{n}},\prec_{B})$ , where $\prec_{B}$ is the partial order defined by inclusion. We say a poset, $(\mathcal{P},\prec)$ , is an $n$ -dimensional Boolean lattice if $(\mathcal{P},\prec)$ is order isomorphic to the standard $n$ -dimensional Boolean lattice and we write $\prec_{B}$ in place of $\prec$ .

Estimating the number of linear extensions for Boolean lattices was first considered in [45] which established a nontrivial upper bound. Later, Brightwell and Tetali [7] proved the following result that essentially settles the question for all practical considerations. If $Q(n)$ denotes the number of linear extensions of an $n$ -dimensional Boolean lattice, then

\frac{\log Q(n)}{2^{n}}=\log\binom{n}{\lfloor\nicefrac{{n}}{{2}}\rfloor}-\frac{3}{2}\log e+\mathcal{O}\left(\frac{\ln n}{n}\right).

(4)

The estimate in Equation (4) illustrates a major challenge in solving the instances of AC-LEP of interest in this paper. Consider an instance of AC-LEP given by $\left(\mathcal{P},\prec_{B},\Xi\right)$ where $(\mathcal{P},\prec_{B})$ is an $n$ -dimensional Boolean lattice and $\Xi$ is any evaluation domain. Suppose we had a black box for efficiently computing all linear extensions of a Boolean lattice denoted by $L\subset S_{K+1}$ . Furthermore, assume we also had a “CAD”-like algorithm which could efficiently check if $\Xi_{\sigma}\neq\emptyset$ . Then, one would need to call this algorithm only $\#L$ -many times as opposed to $(K+1)!$ as we argued above. However, the growth estimate in Equation (4) implies that 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. The first is the linearly constrained linear extension problem (LC-LEP), in which $\mathcal{P}$ is a set of linear polynomials, $\Xi$ is the interior of a polyhedral cone (i.e. $\Xi$ is defined by a finite number of linear inequalities), and $\prec$ is an arbitrary partial order. We present an efficient algorithm for 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 precise definition), is motivated by an application from systems biology described in Section 2.

Definition 1.1.

For $n\in\mathbb{N}$ , an interaction function of order $n$ is a polynomial in $n$ variables, $z=\left(z_{1},\dotsc,z_{n}\right)$ , of the form

f(z)=\prod_{j=1}^{q}f_{j}(z)

(5)

where each factor has the form

f_{j}(z)=\sum_{i\in I_{j}}z_{i}

and the indexing sets $\left\{I_{j}:1\leq j\leq q\right\}$ form a partition for $\left\{1,\dots,n\right\}$ . We define the interaction type of $f$ to be $\mathbf{n}:=\left(n_{1},\dotsc,n_{q}\right)$ 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 not matter for any of the analysis in this paper. Therefore, for convenience in reporting results (see Section 5) we will always assume that

n_{1}\geq n_{2}\geq\cdots\geq n_{q}.

To define an instance of the PSD problem, fix an interaction function $f$ of order $n$ and let $\mathcal{P}$ be the collection of polynomials in the $2n$ positive real variables, $\left\{\ell_{i},\delta_{i}:1\leq i\leq n\right\}$ , obtained by evaluating $f(z)$ with each $z_{i}\in\left\{\ell_{i},\ell_{i}+\delta_{i}\right\}$ . Taking all possible combinations of $z_{i}$ for $1\leq i\leq n$ produces the polynomials for the PSD problem,

\mathcal{P}=\left\{p_{0},\ldots,p_{2^{n}-1}\right\}\subset\mathbb{R}[\ell_{1},\ldots,\ell_{n},\delta_{1},\ldots,\delta_{n}].

(6)

In Section 4, we will define an indexing map between the $2^{n}$ elements of $\mathcal{P}$ and the standard $n$ -dimensional Boolean lattice. Let $\prec_{B}$ denote the Boolean lattice partial order with respect to this index map, and set $\Xi=(0,\infty)^{2n}$ . The PSD problem is the instance of the AC-LEP defined by $\left(\mathcal{P},\prec_{B},\Xi\right)$ . In Section 4, we prove that $(\mathcal{P},\prec_{B},\Xi)$ satisfies Equation (1). However, we present some examples before continuing.

Example 1.3.

The simplest nonlinear PSD problem arises from the interaction function

f(z)=(z_{1}+z_{2})z_{3}

which has interaction type, $\mathbf{n}=(2,1)$ . The PSD polynomials for this interaction function are given by

$\displaystyle p_{0}$	$\displaystyle=(\ell_{1}+\ell_{2})\ell_{3}$	$\displaystyle p_{4}$	$\displaystyle=(\ell_{1}+\ell_{2}+\delta_{1})\ell_{3}$
$\displaystyle p_{1}$	$\displaystyle=(\ell_{1}+\ell_{2})(\ell_{3}+\delta_{3})$	$\displaystyle p_{5}$	$\displaystyle=(\ell_{1}+\ell_{2}+\delta_{1})(\ell_{3}+\delta_{3})$
$\displaystyle p_{2}$	$\displaystyle=(\ell_{1}+\ell_{2}+\delta_{2})\ell_{3}$	$\displaystyle p_{6}$	$\displaystyle=(\ell_{1}+\ell_{2}+\delta_{1}+\delta_{2})\ell_{3}$
$\displaystyle p_{3}$	$\displaystyle=(\ell_{1}+\ell_{2}+\delta_{2})(\ell_{3}+\delta_{3})$	$\displaystyle p_{7}$	$\displaystyle=(\ell_{1}+\ell_{2}+\delta_{1}+\delta_{2})(\ell_{3}+\delta_{3}).$

The PSD evaluation domain is $\Xi=(0,\infty)^{6}$ and the partial order, $\prec_{B}$ , is imposed on $\mathcal{P}$ by identifying $p_{i}$ with the vertex of a unit cube whose coordinates are $(i)_{2}\in\mathbb{F}_{2}^{3}$ where $(i)_{2}$ is the binary expansion of $i$ . The solution to this PSD problem is the set of admissible linear extensions of $(\mathcal{P},\prec_{B})$ , such that $\sigma\in\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{6}\right)$ if and only if $\Xi_{\sigma}\neq\emptyset$ . We note that there are $8!=40,320$ linear orders on $\mathcal{P}$ . However, only $48$ of these are linear extensions of $\prec_{B}$ , and of these, only the following $20$ linear extensions are admissible.

\begin{aligned} (0,1,2,3,4,5,6,7)\\ (0,1,2,3,4,6,5,7)\\ (0,1,2,4,3,5,6,7)\\ (0,1,2,4,3,6,5,7)\\ (0,4,2,6,1,5,3,7)\end{aligned}\quad\begin{aligned} (0,1,2,4,6,3,5,7)\\ (0,1,4,2,5,3,6,7)\\ (0,1,4,2,5,6,3,7)\\ (0,1,4,2,6,5,3,7)\\ (0,4,1,5,2,6,3,7)\end{aligned}\quad\begin{aligned} (0,1,4,5,2,3,6,7)\\ (0,1,4,5,2,6,3,7)\\ (0,2,1,3,4,6,5,7)\\ (0,2,1,4,3,6,5,7)\\ (0,4,2,1,6,5,3,7)\end{aligned}\quad\begin{aligned} (0,2,1,4,6,3,5,7)\\ (0,2,4,1,6,3,5,7)\\ (0,2,4,6,1,3,5,7)\\ (0,4,1,2,5,6,3,7)\\ (0,4,1,2,6,5,3,7)\end{aligned}

The $28$ “missing” linear extensions are those which do not satisfy certain algebraic constraints which are imposed by the polynomial structure. For example, observe that for fixed $\xi\in(0,\infty)^{6}$ , if $p_{3}(\xi)<p_{6}(\xi)$ , then $p_{1}(\xi)<p_{4}(\xi)$ must also hold.

Unlike the partial order which constrains all possible linear extensions, this order relation is conditional. Indeed, there exist choices of $\xi$ such that $p_{3}(\xi)>p_{6}(\xi)$ in which case there is no requirement imposed on the order of $p_{1}(\xi),p_{4}(\xi)$ , and in fact, there are admissible linear extensions which satisfy both choices, e.g. the first two orders in column four. As another example, observe that $p_{5}(\xi)<p_{6}(\xi)$ , if and only if $p_{1}(\xi)<p_{2}(\xi)$ . This algebraic constraint is bi-conditional. However, it too can not be represented in the partial order since both choices occur in at least one admissible order.

To emphasize the role of the interaction function in determining the algebraic constraints, we consider a similar PSD problem that is also an instance of LC-LEP.

Example 1.4.

Consider the interaction type, $\mathbf{n}=(3)$ with corresponding interaction function

f=z_{1}+z_{2}+z_{3}.

As in Example 1.3, we obtain $8$ PSD polynomials given explicitly by

$\displaystyle p_{0}$	$\displaystyle=\ell_{1}+\ell_{2}+\ell_{3}$	$\displaystyle p_{4}$	$\displaystyle=\ell_{1}+\ell_{2}+\ell_{3}+\delta_{1}$
$\displaystyle p_{1}$	$\displaystyle=\ell_{1}+\ell_{2}+\ell_{3}+\delta_{3}$	$\displaystyle p_{5}$	$\displaystyle=\ell_{1}+\ell_{2}+\ell_{3}+\delta_{1}+\delta_{3}$
$\displaystyle p_{2}$	$\displaystyle=\ell_{1}+\ell_{2}+\ell_{3}+\delta_{2}$	$\displaystyle p_{6}$	$\displaystyle=\ell_{1}+\ell_{2}+\ell_{3}+\delta_{1}+\delta_{2}$
$\displaystyle p_{3}$	$\displaystyle=\ell_{1}+\ell_{2}+\ell_{3}+\delta_{2}+\delta_{3}$	$\displaystyle p_{7}$	$\displaystyle=\ell_{1}+\ell_{2}+\ell_{3}+\delta_{1}+\delta_{2}+\delta_{3}$

The evaluation domain and partial order are identical to the PSD problem in Example 1.3. Nevertheless, only the following $12$ linear extensions are admissible

\begin{aligned} (0,1,2,3,4,5,6,7)\\ (0,1,2,4,3,5,6,7)\\ (0,1,4,2,5,3,6,7)\end{aligned}\quad\begin{aligned} (0,1,4,5,2,3,6,7)\\ (0,2,1,3,4,6,5,7)\\ (0,2,1,4,3,6,5,7)\end{aligned}\quad\begin{aligned} (0,2,4,1,6,3,5,7)\\ (0,2,4,6,1,3,5,7)\\ (0,4,1,2,5,6,3,7)\end{aligned}\quad\begin{aligned} (0,4,1,5,2,6,3,7)\\ (0,4,2,1,6,5,3,7)\\ (0,4,2,6,1,5,3,7)\end{aligned}

Similarly, the missing $36$ linear extensions in this example fail to satisfy some algebraic constraints. In both cases, the set of admissible linear extensions is a fraction of the set of all linear extensions of the Boolean lattice. In other words, the algebraic structure implies that the admissible linear extensions are a sparse subset of all linear extensions. The algorithm in this paper exploits the algebraic and order theoretic aspects of the PSD problem to overcome the computational complexity limitations which plague both problems in general. Furthermore, we prove that this algorithm finds all possible linear extensions. For both examples we obtained the ( $12$ and $20$ respectively) admissible solutions without first computing the linear extensions of $(\mathcal{P},\prec_{B})$ and then checking which are admissible.

Related work

As is indicated above our original motivation for this paper arises from problems in systems biology for which explicit complete solutions to the PSD problem are required. As such the majority of this introduction has focused on the question of efficacy of computation. However, there is another direction in which the work of this paper overlaps with other efforts. In particular, observe that the case where the interaction function is linear, i.e. has interaction type $\mathbf{n}=(n)$ , solving the AC-LEP is equivalent to identifying all the cells of a hyperplane arrangement. This latter problem has been the subject of considerable study (see [47] for an introduction) and in particular, Maclagan [33] provides the number of solutions for the linear PSD problem for $n=1,\ldots,7$ . Our computations (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,1,1,1)$ , $(1,1,1,1,1)$ , and $(1,1,1,1,1,1)$ , reported in column two of Table 1, we obtain

\frac{336}{4!}=14=:a_{4},\quad\frac{61920}{5!}=516=:a_{5},\quad\frac{89414640}{6!}=124,187=:a_{6}

which align with sequence A009997 in the OEIS [44]. From [22], we know this sequence represents the number of comparative probability orderings on all subsets of $n$ elements that can arise by assigning a probability distribution to the individual elements. The equivalence of comparative probability orderings and solutions to the linear PSD problem follows directly from the definition of comparative probability.

Organization of paper

The remainder of this paper is organized as follows. In Section 2 we briefly describe how the PSD problem arises naturally in the study of global dynamics for gene regulatory networks. In Section 3, we present an efficient algorithm for solving instances of the LC-LEP. In Section 4, we show that the LC-LEP is related to the PSD problem in the following way. If $\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ is an instance of the PSD problem, then we construct an associated instance of LC-LEP, denoted by $\left(\mathcal{P}^{\prime},\prec_{B},\Xi^{\prime}\right)$ , which satisfies the inclusion

\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)\subseteq\mathcal{T}\left(\mathcal{P}^{\prime},\prec_{B},\Xi^{\prime}\right).

(7)

We refer to this instance of LC-LEP as the linearized PSD problem associated to $\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ . We exploit this construction and the algorithm for solving the LC-LEP presented in Section 3, to provide a means of efficiently computing a collection of candidates that contains the solution to the PSD problem. We prove that in some cases the inclusion in Equation (7) is actually an equality. More generally, this inclusion is strict, but the candidate set is a sparse subset of the collection of all linear extensions of $(\mathcal{P},\prec_{B})$ . In this case we describe algorithms for removing the non-admissible solutions.

Finally, in Section 5 we present the results for all PSD solutions with order up to four. Additionally, we have some results for PSDs of order five and six. For the remaining cases and PSDs of higher order the computations become too large, i.e. we do not have immediate use for them in applications and they require significant computing resources.

2 Dynamic Signatures Generated by Regulatory Networks

This section provides a brief description of how the AC-LEP arises in the context of mathematical modeling of problems from systems biology with a few comments at the end indicating its potential application in more general settings of nonlinear dynamics. Biologists often describe regulatory networks in terms of annotated directed graphs, such as that shown in Figure 1(a) where the labeling of the edges, $n\to m$ or $n\dashv m$ indicates whether node $n$ activates or represses node $m$ . Our goal is to describe the type of dynamics that can be expressed by the regulatory network. This requires two things: (i) a mathematical framework in which we can rigorously discuss dynamics when an analytic expression of the nonlinearity is unknown, and more specifically, (ii) imposing a mathematical interpretation on the regulatory network that is compatible with its use as a biological model.

We discuss both of these topics in the following sections, but hasten to add that there are other complementary approaches to this problem. Perhaps the strongest dichotomy is between Boolean models and ordinary differential equation (ODE) models (see [1] for a brief overview and references). Recent efforts in the ODE direction seek to go from networks to explicit polynomial vector fields (see [27] and references therein). The approach we describe lies somewhere between these two. Our approach is partially combinatorial, and thus we obtain efficacy similar to that of the Boolean models, but we maintain continuous parameterization (hence the necessity of the results of this paper), and therefore our computations can be interpreted in the context of ODEs. However, unlike methods based on explicit vector fields, we derive the structure of the dynamics via order theory and algebraic topology. Understanding the relative strengths and weaknesses of these various methods is still an active area of research.

Figure 1: (a) Example of a regulatory network. (b) Model for edge

n\to m

where

\ell_{m,n}

indicates low level of growth rate of

x_{m}

induced by

x_{n}

and

\ell_{m,n}+\delta_{m,n}

indicates high level of growth rate of

x_{m}

induced by

x_{n}

\theta_{m,n}

provides information about the value of

x_{n}

that lies between low values inducing low and high expression levels.

2.1 Mathematical framework

Classical nonlinear dynamics focuses on invariant sets and bifurcations, both of which are extremely sensitive to parameters, e.g. nonlinearities. Unfortunately, this is precisely the information that is unknown in multiscale systems such as those associated with biology. C. Conley proposed that in these situations one should focus on isolated invariant sets [14] as there is a sheaf structure associated to them [43] and thus are not subject to the above mentioned sensitivity. Furthermore, he developed an algebraic topological invariant for isolated invariant sets, called the Conley index, from which one can recover considerable information about the existence and structure of the invariant set [42]. In particular, using the homology Conley index one can identify the existence of invariant sets [14], equilibria [46, 36], periodic orbits [38], heteroclinic orbits [15], chaotic dynamics [40, 49], and semi-conjugacies onto nontrivial dynamics [39, 37].

One indication of the strength of Conley’s proposition is that it is possible to study a differential equation (ordinary or partial) or a continuous map (finite or infinite dimensional) numerically and with appropriate bounds on the errors, to draw rigorous conclusions about the dynamics from homological Conley index computations [41, 18, 28, 19]. The weakness of the above mentioned results is that they are perturbative in nature, i.e. the results are computed at given parameter values and the results are guaranteed to be true for a sufficiently small neighborhood of that parameter value. In practice one can extract numerical bounds for the size of the neighborhood, but they tend to be small.

For multiscaled systems such as those arising in systems and synthetic biology we need to be able to investigate dynamics over large ranges of parameter values. This was done in the context of low dimensional maps with a low dimensional parameter space [3, 12, 11]. The strategy is to subdivide parameter space uniformly, for each region of parameter space compute dynamics with error bounds valid over the entire region, and once again interpret the structure of the dynamics using the Conley index. While the approach is extremely effective for low dimensional problems it is difficult to extend to higher dimensions for two reasons:

R1: the cost of computing effective numerical error bounds becomes prohibitive, and
R2: the subdivision of parameter space is not chosen according to the underlying dynamics and refinement of the subdivision leads to exponential growth with dimension in the number of subdivisions.

With R1 in mind W. Kalies, R. Vandervorst, and the third author have developed an approach to dynamics that is based on combinatorics, order theory (posets and lattices), and algebraic topology. The combinatorics arises from the choice of a decomposition of phase space into closed regular sets and the dynamics is characterized by a relation on this collection of closed regular sets. It is proven that Conley’s fundamental characterization of dynamics in terms of attractor-repeller pairs or equvalently Morse decompositions can be completely recovered via this combinatorial approach [29, 30, 31]. The algebra of the order theory, in particular Birkhoff’s theorem relating finite posets with finite distributive lattices, provides an algorithmic approach for the identification of index pairs from which the Conley index can be computed [32]. This in turn allows one to obtain a chain complex that codifies gradient-like structure of the dynamics and the Conley indices of all the associated isolated invariant sets [25]. Observe that, at least conceptually, we have replaced the numerical computations by combinatorial computations. In practice, at least for the biological regulatory networks discussed below, this combinatorial approach is extremely efficient, both with respect to time and memory [17].

Dealing with R2 is the raison d’etre for this paper. As is made clear below for the regulatory networks of interest our closed regular sets take the form of cubes determined by parameter values. The relation on these cubes that represents the dynamics is obtained by determining the directions of inequalities. All possible sets of directions are in turn determined by understanding the orders of the polynomials $\mathcal{P}$ defined by (6). Thus, solving the AC-LEP provides us with an a priori optimal potentially nonlinear decomposition of parameter space.

2.2 Mathematical interpretation of regulatory networks

We now turn to our mathematical interpretation of a regulatory network that is compatible with its use as a biological model. With this in mind, we assign to node $m$ a state variable, $x_{m}>0$ , that corresponds to the quantity of a protein, mRNA, or a signaling molecule. Precise nonlinear expressions for the interactions of the variables are not assumed to be known, but we do assume that the sign of the rate of change of $x_{m}$ is determined by the expression

-\gamma_{m}x_{m}+\Lambda_{m}(x)

(8)

where $\gamma_{m}$ indicates the decay rate and $\Lambda_{m}$ is a parameter dependent function that characterizes the rate of growth of $x_{m}$ . Note that $\Lambda_{m}$ is a function of $x_{n}$ if and only if there exists an edge from $n$ to $m$ in the regulatory network.

Since the biological model provides minimal information about the effect of $x_{n}$ on $x_{m}$ we want to 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<\ell_{m,n}<\ell_{m,n}+\delta_{m,n}$ . Figure 1(b) corresponds to an edge $n\to m$ and $\theta_{m,n}$ indicates that the rate of increase $\ell_{m,n}$ must occur at some lesser value of $x_{n}$ and the rate of increase $\ell_{m,n}+\delta_{m,n}$ must occur at some greater value of $x_{n}$ . An arrow of the form $n\dashv m$ leads to the opposite relation.

This introduces three positive parameters, $\ell_{m,n}$ , $\delta_{m,n}$ , and $\theta_{m,n}$ , for each edge in the regulatory network. Note that this is the minimal number of parameters that allows one to quantify the assumption that $x_{n}$ activates $x_{m}$ (or equivalently that $x_{n}$ represses $x_{m}$ ). We encode this information with the following functions

\lambda^{+}_{m,n}(x_{n}):=\begin{cases}\ell_{m,n}&\text{if $x_{n}<\theta_{m,n}$}\\ \ell_{m,n}+\delta_{m,n}&\text{if $x_{n}>\theta_{m,n}$}\end{cases}\quad\text{and}\quad\lambda^{-}_{m,n}(x_{n}):=\begin{cases}\ell_{m,n}+\delta_{m,n}&\text{if $x_{n}<\theta_{m,n}$}\\ \ell_{m,n}&\text{if $x_{n}>\theta_{m,n}$.}\end{cases}

(9)

We do not assume that the values of $\ell_{m,n},\delta_{m,n}$ , or $\theta_{m,n}$ are known (see Figure 1(b)). This is intentional as many of these parameters do not have an easy biological interpretation and/or correspond to physical constants which are difficult or impossible to precisely measure. Thus, the goal is not to determine the dynamics at any particular choice of parameters, but to determine the range and robustness of the qualitative dynamics exhibited by a network.

A regulatory network such as that of Figure 1(a) does not indicate how multiple inputs to a particular node should be processed. An approach that is used is to assume a simple algebraic relationship involving sums and products of the $\lambda^{\pm}$ . As an example, a reasonable choice for $x_{1}$ of Figure 1(a) is

\Lambda_{1}(x_{2},x_{3},x_{4})=\left(\lambda^{+}(x_{2})+\lambda^{+}(x_{3})\right)\lambda^{-}(x_{4}).

(10)

Observe that each $\lambda^{\pm}$ takes only two values and therefore, generically Equation (10) takes $8$ distinct values which are precisely the values of the elements of $\mathcal{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 $\lambda^{\pm}$ , or $\Lambda$ as literal expressions of the nonlinear interactions, but rather, that the associated parameter values are landmarks of whatever of the “true” nonlinear function is. This has several consequences.

1.

We cannot expect Equation (8) 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 $\gamma$ , $\theta$ , $\ell$ and $\delta$ .
2.

The only values of $x_{m}$ at which the dynamics of $x_{n}$ change are of the form $\theta_{m,*}$ . The associated hyperplanes $x_{j}=\theta_{k,j}$ decompose phase space $[0,\infty)^{N}$ , where $N$ is the number of nodes in the network, into $N$ -dimensional rectangular regions called domains.
3.

Since we have determined $\mathcal{T}\left(\left\{p_{0},\ldots,p_{7}\right\},\prec,(0,\infty)^{6}\right)$ for (10) we can determine all possible signs of (8) associated with (10) by cataloguing the relative values of $\gamma_{1}\theta_{j,1}$ , $j=4,5$ with respect to $\left\{p_{0},\ldots,p_{7}\right\}$ .

The Dynamic Signatures Generated by Regulatory Networks (DSGRN) library contains software that, given the information of the form provided by consequence 3, is capable of efficiently building a database of the global dynamics of a regulatory network over all of parameter space [17, 34]. In [17] the authors considered networks whose nodes have at most three in-edges, and at most three out-edges. This constraint was due to the difficulty in solving the PSD problem and the results of this paper provide a means to vastly expand the capabilities of DSGRN [34]. 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.

We remark that DSGRN is proving to be a versatile tool in the realm of systems and synthetic biology with applications in the realm of regulatory networks relevant to oncology [17, 24, 54], developmental biology [20], yeast cell cycle [16], and synthetic biology [23].

2.3 Extensions

The focus of this paper is on the PSD problem because of the immediate need, arising from applications, to expand the capabilities of DSGRN to handle regulatory networks with nodes that have more than three in and out edges. However, as should be clear from Section 2.1 the mathematical foundations for our approach is quite general. In particular, the dynamics associated with functions described in Section 2.2 or more generally with the PSD problem (See Definition 4.6) represent a rather special subclass. Two immediate generalizations that are being pursued in other work, but rely on application of the techniques of this work are as follows.

The first generalization is to replace the constant $\gamma_{m}$ in (8) by a function $\gamma_{m}(x)$ . This allows the DSGRN strategy to be applied in the systems biology setting to regulatory networks that involve post-transcriptional regulation. The second generalization is to replace the functions $\lambda^{\pm}_{m,n}$ given in (9) by a step function involving multiple steps without necessarily assuming monotonicity. In this setting we no longer have the Boolean lattice as the minimal partial order, but the essential arguments of this paper still apply.

3 Solving the LC-LEP

In this Section, we provide an efficient algorithm to solve the LC-LEP defined in Section 1. Note that if $q\in\mathbb{R}[x_{1},\dotsc,x_{d}]$ is a linear polynomial, then evaluation of $q$ defines a linear functional on $\mathbb{R}^{d}$ . Thus, there exists a unique vector $\mathbf{u}_{q}\in\mathbb{R}^{d}$ , that we call the representation vector for $q$ , satisfying

q(\xi)=\mathbf{u}_{q}\cdot\xi\qquad\text{for all}\quad\xi\in\mathbb{R}^{d}.

Recall that the evaluation domain for the LC-LEP is the interior of a polyhedral cone. Thus, there exists a collection of linear polynomials $\mathcal{Q}_{\Xi}$ , such that

\Xi=\left\{\xi\in\mathbb{R}^{d}:\xi\cdot\mathbf{u}_{q}>0\ \text{for all}\quad q\in\mathcal{Q}_{\Xi}\right\}.

(11)

We assume that (11) is satisfied for the remainder of this section.

To foreshadow our approach recall that by definition $\sigma\in\mathcal{T}(\mathcal{P},\prec,\Xi)$ if and only if $\Xi_{\sigma}\neq\emptyset$ . Our approach to determining the latter is to recast it in the language of linear algebra on cones in $\mathbb{R}^{d}$ . From this perspective, the problem is equivalent to rigorously solving a linear programming problem 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.

3.1 Cones

Definition 3.1.

A subset $C\subset\mathbb{R}^{d}$ is a cone if $v\in C$ and $\theta\in[0,\infty)$ implies that $\theta v\in C$ . The cone $C$ is pointed if it is closed, convex, and satisfies

C\cap-C=C\cap\left\{-v:v\in C\right\}　=0.

(12)

Observe that (12) implies that a pointed cone does not contain any lines. A vector $v\in\mathbb{R}^{d}$ is a conic combination of $\{v_{1},\dots,v_{k}\}\subset\mathbb{R}^{d}$ if $v=\theta_{1}v_{1}+\dots+\theta_{k}v_{k}$ where $\theta_{1},\dotsc,\theta_{k}\geq 0$ . Suppose $V=\left\{v_{1},\dotsc,v_{k}\right\}\subset\mathbb{R}^{d}$ . The conic hull of $V$ is given by

\operatorname{cone}(V):=\left\{\sum_{i=1}^{k}\theta_{i}v_{i}:0\leq\theta_{i},\ i=1,\dotsc,k\right\}.

The following result is left to the reader to check.

Proposition 3.2.

Given $V=\left\{v_{1},\dotsc,v_{k}\right\}\subset\mathbb{R}^{d}$ , $\operatorname{cone}(V)$ is the smallest closed convex cone that contains $V$ .

We make use of the following propositions.

Proposition 3.3.

Suppose $V=\{v_{0},\dots,v_{m}\}\subset\mathbb{R}^{d}$ is a collection of nonzero vectors such that $\operatorname{cone}(V)$ is a pointed cone. Then, there exists some $v^{\prime}\in\mathbb{R}^{d}$ such that $v^{\prime}\cdot v_{i}>0$ for all $0\leq i\leq m$ .

Proof.

Observe that $-v_{m}\notin\operatorname{cone}(\left\{v_{0},\dots,v_{m-1}\right\})\subset\operatorname{cone}(V)$ since $\operatorname{cone}(V)$ is pointed. Hence $\left\{-v_{m}\right\}$ and $\operatorname{cone}(\left\{v_{0},\dots,v_{m-1}\right\})$ are disjoint, convex, closed subsets of $\mathbb{R}^{d}$ .

Therefore, by the hyperplane separation theorem [5], there exists $v^{\prime}\in\mathbb{R}^{d}$ such that $v^{\prime}\cdot-v_{m}<0$ and $v^{\prime}\cdot v>0$ for any $v\in\operatorname{cone}(\left\{v_{0},\dotsc,v_{m-1}\right\})$ . ∎

Proposition 3.4.

Suppose $V=\left\{v_{0},\dots,v_{m}\right\}\subset\mathbb{R}^{d}$ is a collection of nonzero vectors such that $\operatorname{cone}(V)$ is a pointed cone. If $-v\notin\operatorname{cone}(V)$ , then $\operatorname{cone}(V\cup\{v\})$ is pointed.

Proof.

Suppose that $\operatorname{cone}(V\cup\{v\})$ is not pointed. Then, there exists $w\neq 0$ such that $w,-w\in\operatorname{cone}(V\cup\{v\})$ or equivalently

w=\sum_{i=0}^{m}\alpha_{i}v_{i}+\alpha v\quad\text{and}\quad-w=\sum_{i=0}^{m}\beta_{i}v_{i}+\beta v

where $\alpha_{i},\beta_{i},\alpha,\beta$ are all nonnegative. Note that if $\alpha=\beta=0$ , then $\pm w\in\operatorname{cone}(V)$ , which contradicts the assumption that $V$ is pointed. The sum of the two equations above is

-(\alpha+\beta)v=\sum_{i=0}^{m}(\alpha_{i}+\beta_{i})v_{i}.

This implies that $-v\in\operatorname{cone}(V)$ , contradicting the assumption that $\operatorname{cone}(V)$ is pointed. ∎

The previous propositions illustrate the importance of solving the cone inclusion problem: given a vector $v\in\mathbb{R}^{d}$ , and finite set of vectors $V\subset\mathbb{R}^{d}$ , determine whether or not $v\in\operatorname{cone}(V)$ . Algorithm 1, stated below, solves this problem. Observe that checking if $v\in\operatorname{cone}(V)$ is equivalent to solving the following linear programming feasibility problem.

$\displaystyle\text{Does there exist}\quad\alpha$
$\displaystyle\text{such that}\quad\textbf{V}\alpha$	$\displaystyle=v$	(13)
$\displaystyle\text{and}\quad\alpha$	$\displaystyle\geq 0?$

where $\mathbf{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 [52]. The results can be made rigorous by performing computations using interval arithmetic [48] or rational linear programming [2] in the case that $\mathbf{V}$ is rational. Observe that the PSD problem defined in Section 1 satisfies this constraint. As is made clear in Section 5, we use different solvers depending on the machine employed to do the computations. In any case, we take for granted the existence of a rigorous solver for the feasibility problem in Equation 3.1 as a “black box” which we call LPSolver which is employed in the following algorithm.

Input:

v,V=\{v_{1},\dots,v_{m}\},{\tt LPSolver}

Output:

\textbf{True},\textbf{False}

Result : Return True if

v\in\operatorname{cone}(V)

otherwise False

1 Function InCone( $v,V,\textnormal{{LPSolver}}$ ):

2 Return LPSolver( $v,V$ )

End Function

Algorithm 1 Cone inclusion

The next algorithm uses Proposition 3.4 (see line 4) and Algorithm 1 to determine if a set of vectors defines a pointed cone.

Input:

V=\{v_{1},\dots,v_{m}\}

Output:

\textbf{True},\textbf{False}

Result : Return True if

\operatorname{cone}(V)

is pointed otherwise False

1 Function CheckCone( $V$ ):

V^{\prime}=\{v_{1}\}

3 for i = 2 … m do

4 if InCone( $-v_{i},V^{\prime}$ ) then

5 Return False

6 else

V^{\prime}=V^{\prime}\cup\{v_{i}\}

8 end if

10 end for

11 Return True

End Function

Algorithm 2 Cone pointedness

We now show that the LC-LEP can be reformulated as a problem of identifying whether some specific subsets of vectors generate pointed cones.

Definition 3.5.

Given an instance of the LC-LEP, $(\mathcal{P},\prec,\Xi)$ , we define the base cone as $\operatorname{cone}(V_{0}):=\operatorname{cone}(V_{\Xi}\cup V_{\prec})$ where $V_{\Xi}$ and $V_{\prec}$ are defined as follows. Set

V_{\Xi}:=\left\{\mathbf{u}_{q}:q\in\mathcal{Q}_{\Xi}\right\},

where $\mathcal{Q}_{\Xi}$ are the representation vectors as defined in Equation (11). Applying Algorithm 2 to $V_{\Xi}$ (and the fact that we assume $\Xi\neq\emptyset$ ) shows that $\operatorname{cone}(V_{\Xi})$ is pointed. Define

V_{\prec}:=\left\{\mathbf{u}:\text{$\mathbf{u}$ is the representing vector of $p-q$ where $q\prec p$ and $p,q\in\mathcal{P}$}\right\}.

Observe that if $(\mathcal{P},\prec)$ satisfies Equation (1), then the representation vector for $p-q$ is an element of $V_{\prec}$ by definition. Therefore, by Proposition 3.3, $V_{\prec}$ is pointed.

The motivation behind our definition of $V_{0}$ is that it characterizes the algebraic constraints in the LC-LEP in terms of linear algebra that can be efficiently checked. The next proposition shows that the same idea works for linear extensions.

Given $\sigma\in S_{K+1}$ , we define

V_{\sigma}:=V_{0}\cup\left\{\mathbf{u}_{p_{\sigma(i+1)}}-\mathbf{u}_{p_{\sigma(i)}}:p_{\sigma(i)}\in\mathcal{P},\ i=0,\dotsc,K-1\right\}.

(14)

Proposition 3.6.

For any $\sigma\in S_{K+1}$ , $\Xi_{\sigma}\neq\emptyset$ if and only if $\operatorname{cone}(V_{\sigma})$ is a pointed cone.

This Proposition is a trivial application of a well known result in convex optimization. If $\mkern 1.5mu\overline{\mkern-1.5mu\Xi\mkern-1.5mu}\mkern 1.5mu$ denotes the closure of $\Xi$ , then $\mkern 1.5mu\overline{\mkern-1.5mu\Xi\mkern-1.5mu}\mkern 1.5mu$ is a convex cone. Its dual cone is defined to be the set

\left\{y\in\mathbb{R}^{d}:y\cdot\xi\geq 0\ \text{for all}\ \xi\in\mkern 1.5mu\overline{\mkern-1.5mu\Xi\mkern-1.5mu}\mkern 1.5mu\right\}

and we observe that this set is nothing more than $\operatorname{cone}(V_{\Xi})$ . Similarly, for each $\sigma\in S_{K+1}$ , $\mkern 1.5mu\overline{\mkern-1.5mu\Xi_{\sigma}\mkern-1.5mu}\mkern 1.5mu$ is again a convex cone as its simply a restriction of $\mkern 1.5mu\overline{\mkern-1.5mu\Xi\mkern-1.5mu}\mkern 1.5mu$ obtained by imposing finitely many additional linear constraints and $\operatorname{cone}(V_{{}_{\sigma}})$ is its associated dual cone. Consequently, Proposition 3.6 follows from the fact that a convex cone is pointed if and only if its dual cone is nonempty. A proof of this result can be found in [21]. We emphasize that the importance of Proposition 3.6 is the implied equivalence

\mathcal{T}(\mathcal{P},\prec,\Xi)=\left\{\sigma:\Xi_{\sigma}\neq\emptyset\right\}=\left\{\sigma:\operatorname{cone}(V_{\sigma})\ \text{is pointed}\right\}.

3.2 An algorithm for identifying $\mathcal{T}(\mathcal{P},\prec,\Xi)$

In this section we present an algorithm for solving an arbitrary instance of the LC-LEP.

Input:

\sigma_{\operatorname{part}}=[\text{ }],\mathcal{P},V=V_{0},{\tt Ret}=\{\}

Output:

\mathcal{T}(\mathcal{P},\prec,\Xi)

Result : Ret: collection of all linearly realizable total order under restriction of

V

3 Function OrderingGenerator( $\sigma_{\operatorname{part}},\mathcal{P},V,{\tt Ret}$ ):

5 if $\sigma_{\operatorname{part}}==[\text{ }]$ and CheckCone(V) is not True then

6 Return

7 end if

l+1

= length of

\sigma_{\operatorname{part}}

10 if l == K then

11 add

\sigma_{\operatorname{part}}

{\tt Ret}

12 Return

13 end if

15 for i = 0 .. K do

16 if $i\not\in\sigma_{\operatorname{part}}$ then

\mathbf{u}^{\prime}=\mathbf{u}_{p_{i}}-\mathbf{u}_{p_{\sigma_{\operatorname{part}}(l)}}

18 if not InCone( $-\mathbf{u}^{\prime},V$ ) then

19 OrderingGenerator( $\sigma_{\operatorname{part}}+[i],\mathcal{P},V\cup\{v^{\prime}\},{\tt Ret}$ )

20 end if

22 end if

24 end for

End Function

Algorithm 3 LC-LEP solver

In the Algorithm 3, for convenience, we take $\mathbf{u}_{p_{\sigma_{\operatorname{part}}(-1)}}=0$ . To prove the correctness of the algorithm it is useful to denote the return of Algorithm 3 given input $(\mathcal{P},\prec,\Xi)$ as $\mathcal{T}_{alg}(\mathcal{P},\prec,\Xi)$ .

Definition 3.7.

For fixed $(\mathcal{P},\prec,\Xi)$ , $\sigma\in S_{K+1}$ and for $k=1,\cdots,K$ , define

V_{\sigma,k}=\{\mathbf{u}_{p_{\sigma(i)}}-\mathbf{u}_{p_{\sigma(i-1)}}\}_{i=1,\dots,k}\cup V_{0},

where $V_{0}$ is the base cone for $(\mathcal{P},\prec,\Xi)$ as in Definition 3.5, and $\mathbf{u}_{p_{\sigma(j)}},j=0,\ldots,K$ is the representation vector of $p_{j}\in\mathcal{P}$ . For convenience, we define $V_{\sigma,0}=V_{0}$ and we observe that $V_{\sigma,K}=V_{\sigma}$ from Equation (14).

Theorem 3.8.

Algorithm 3 solves the LC-LEP.

Proof.

Given $(\mathcal{P},\prec,\Xi)$ , we need to show that $\mathcal{T}(\mathcal{P},\prec,\Xi)=\mathcal{T}_{alg}(\mathcal{P},\prec,\Xi)$ . We may assume that $\operatorname{cone}(V_{0})$ is pointed since if not, then both $\mathcal{T}_{alg}(\mathcal{P},\prec,\Xi)$ and $\mathcal{T}(\mathcal{P},\prec,\Xi)$ are empty.

We first show that $\mathcal{T}_{alg}(\mathcal{P},\prec,\Xi)\subset\mathcal{T}(\mathcal{P},\prec,\Xi)$ , i.e. for any $\sigma\in\mathcal{T}_{alg}(\mathcal{P},\prec,\Xi)$ we show that the set $\Xi_{\sigma}\neq\emptyset$ . As indicated above we assume $\operatorname{cone}(V_{0})=\operatorname{cone}(V_{\sigma,0})$ is pointed. For Algorithm 3, lines 2-4 returns the empty set if $\operatorname{cone}(V_{0})$ is not pointed. Otherwise, it passes to lines 5-9 which checks if $\sigma_{\operatorname{part}}$ is a total order over $\{0,\dots,K\}$ . If so, it is added to the return variable, Ret. If $\sigma_{\operatorname{part}}$ is not a total order, then lines 10-17 extend it to a total order by recursively constructing $V_{\sigma,i}$ from $V_{\sigma,i-1}$ for $1\leq i\leq K$ .

Therefore, it suffices to show that $V_{\sigma,k}$ are all pointed for $k=1,\dots,K$ .

Fix $k\in\left\{1,\dotsc,K\right\}$ . In lines 11-12, we find a candidate $i\in\{0,\dots,K\}$ which is not in the image of $\sigma_{\operatorname{part}}$ , and define $V_{\sigma,k+1}=V_{\sigma,k}\cup\{\mathbf{u}^{\prime}\}$ where $\mathbf{u}^{\prime}=\mathbf{u}_{p_{i}}-\mathbf{u}_{p_{\sigma(k)}}$ . In line 13, we verify that $-\mathbf{u}^{\prime}\notin V=V_{\sigma,k-1}$ and it follows from Proposition 3.4, that $\operatorname{cone}(V_{\sigma,k})$ is pointed. For each $\sigma$ appended to Ret, we have $\operatorname{cone}(V_{\sigma,k})$ is pointed for $k=0,\dots,K$ . In particular, $\operatorname{cone}(V_{\sigma,K})=\operatorname{cone}(V_{\sigma})$ is pointed, and from Proposition 3.6, we have $\Xi_{\sigma}\neq\emptyset$ .

We now prove that $\mathcal{T}(\mathcal{P},\prec,\Xi)\subset\mathcal{T}_{alg}(\mathcal{P},\prec,\Xi)$ . Assume that $\sigma\in\mathcal{T}(\mathcal{P},\prec,\Xi)$ . By definition this means that $\Xi_{\sigma}\neq\emptyset$ and from Proposition 3.6, $V_{\sigma}$ is pointed. For each $k=1,\dots,K$ , we have $V_{\sigma,k}\subset V_{\sigma}$ , and thus $V_{\sigma,k}$ is pointed. As $V_{\sigma,k}$ is pointed, we know $-(\mathbf{u}_{p_{\sigma(k)}}-\mathbf{u}_{p_{\sigma(k-1)}})\notin V_{\sigma,k-1}$ for $k=1,\dots,K$ . Therefore, line 13 in Algorithm 3 will not fail to extend $\sigma$ at each step in the recursion and after $K$ recursive extensions, $\sigma$ will be appended to Ret and thus, $\sigma\in\mathcal{T}_{alg}(\mathcal{P},\prec,\Xi)$ . ∎

3.3 Complexity analysis

As discussed in Section 1, computing linear extensions of a poset is at least as hard as counting them, which has been established to be a difficult problem. While the algebraic constraints induced by the evaluation domain reduces the number of admissible linear extensions, it is not clear whether or not the constrained problem is easier (i.e. whether or not the LC-LEP is also $\#P$ -complete). The answer may depend on the partial order, the evaluation domain, or both. In any case, this appears to be a difficult and open problem which is outside the scope of this work. Consequently, we do not attempt to provide any rigorous asymptotic estimates of Algorithm 3.

However, for the instances of LC-LEP solved in this work, our algorithm is dramatically more efficient than a brute force algorithm as evidenced by the results in Section 5. By brute force we mean an algorithm which checks every linear order on $\mathcal{P}$ to determine whether or not it is admissible. Thus we conjecture that Algorithm 3 is more efficient in the typical case which explains our results.

A heuristic, but not rigorous, justification is as follows. We start by assuming a fixed implementation for the InCone function in Algorithm 1. This is the core algorithm in the sense that it dominates the computational cost of Algorithm 3. Let $(\mathcal{P},\prec,\Xi)$ be a given instance of the LC-LEP where $\mathcal{P}=\left\{p_{0},\dotsc,p_{K}\right\}\subset\mathbb{R}[x_{1},\dots,x_{d}]$ and let $\mathcal{T}(\mathcal{P},\prec,\Xi)\subseteq S_{K+1}$ denote the solution of the LC-LEP.

We consider the number of InCone function calls required to solve this instance (i.e. we suppose that each call to InCone has unit computational cost). It is trivial to count the number of calls to InCone necessary for a brute force search. By Proposition 3.6, $\Xi_{\sigma}$ is nonempty for any $\sigma\in S_{K+1}$ if and only $\operatorname{cone}(V_{\sigma})$ is pointed which is determined by a single InCone call. Therefore, a brute force search recovers $\mathcal{T}(\mathcal{P},\prec,\Xi)$ in exactly $(K+1)!$ calls to InCone.

To contrast this with Algorithm 3, consider a fixed $\sigma\in S_{K+1}$ which is not admissible. Hence, there exists a least index, $i\in\left\{0,\dotsc,K\right\}$ , such that

\sigma^{\prime}\big{|}_{\left\{0,\dotsc,i-1\right\}}=\sigma\big{|}_{\left\{0,\dotsc,i-1\right\}}\qquad\text{for at least one }\sigma^{\prime}\in\mathcal{T}(\mathcal{P},\prec,\Xi)

and

\sigma^{\prime}\big{|}_{\left\{0,\dotsc,i\right\}}\neq\sigma\big{|}_{\left\{0,\dotsc,i\right\}}\qquad\text{for any }\sigma^{\prime}\in\mathcal{T}(\mathcal{P},\prec,\Xi).

Consequently, in the $i^{\rm th}$ recursive call, line $13$ of Algorithm 3 returns False and $\sigma$ is “pruned” (i.e. removed from consideration as a candidate solution). In other words, $\sigma$ has already been ruled incompatible with the algebraic constraints imposed by $\prec$ and $\Xi$ using only its partial image obtained by restriction to the subset $\left\{0,\dotsc,i\right\}$ . The key observation is that this applies to every such incompatible linear extension. In other words, every $\sigma^{\prime}\in S_{K+1}$ satisfying

\sigma^{\prime}\big{|}_{\left\{0,\dotsc,i\right\}}=\sigma\big{|}_{\left\{0,\dotsc,i\right\}}

is also pruned in this step. Evidently, there are $(K-i)!$ such extensions which are simultaneously removed which results in saving $(K-i)!$ calls to InCone compared to the brute force search.

Pruning incompatible partial images early leads to an exponential reduction in the number of InCone calls which explains the results shown in Table 1. However, a more rigorous analysis for the complexity of Algorithm 3 amounts to estimating how early a given incompatible partial image is pruned. It is likely that this is heavily dependent on either $\prec$ or $\Xi$ , or both, and may also depend on the order in which the partial images are extended.

4 Solving the general PSD problem

In this section we present a solution for the PSD problem described in Section 1. The solution is based on the observation that the PSD problem naturally has a Boolean lattice structure. Therefore, the linear PSD problem is an instance of the LC-LEP and and reduces to a simple application of the algorithms presented in Section 3. For nonlinear PSD problems, we construct a map that “embeds” it into an instance of LC-LEP (of higher dimension) in the sense that the inclusion in Equation (7) holds for the Boolean partial order. We prove a sufficient condition for which this inclusion is equality and describe a method for disqualifying spurious solutions when it is strict.

4.1 The PSD as an instance of AC-LEP

Throughout this section $\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ denotes a PSD problem for a fixed interaction function $f$ of order type $\mathbf{n}\in\mathbb{N}^{q}$ as defined in Equation (6) where $\prec_{B}$ is the Boolean lattice partial order. Our first goal is to show that $\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ satisfies Equation (1), and in particular, that every $\sigma\in\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ is a linear extension of a Boolean lattice. We start by defining appropriate indices for the elements of $\mathcal{P}$ .

Definition 4.1.

Suppose $\mathbf{n}\in\mathbb{N}^{q}$ is the interaction type for $f\in\mathbb{R}\left[z_{1},\dotsc,z_{n}\right]$ . As in Definition 1.1 let $\left\{I_{1},\dotsc,I_{q}\right\}$ denote the indexing sets for each summand of $f$ which form an integer partition of $I:=\left\{1,\dotsc,n\right\}$ . For each $1\leq j\leq q$ , we consider $I_{j}$ as an ordered set and and for $1\leq k\leq n_{j}$ , we let $I_{j}(k)$ denote the $k^{\rm th}$ largest element of $I_{j}$ . Let $E:=\left\{\alpha:\left\{1,\dotsc,n\right\}\to\left\{0,1\right\}\right\}$ be the set of all Boolean functions defined on $I$ . The Boolean indexing map, denoted by $B:E\to\left\{0,\dotsc,2^{n-1}\right\}$ , is defined by the formula

B(\alpha):=\sum_{j=1}^{q}\sum_{k=1}^{n_{j}}\alpha(I_{j}(k))2^{n-I_{j}(k)}.

In other words, $E$ and $B$ are defined so that for each $d\in I$ , $\alpha(d)$ is the $d^{\rm th}$ binary digit of $B(\alpha)$ in Little-endian order.

We will also consider, $\alpha\in E$ , as a vector of Boolean functions defined as follows. Let $E_{j}$ denote the set of Boolean functions defined on $I_{j}$ . Then, elements of $E$ can be represented as vectors of the form

\alpha=\left(\alpha_{1},\dotsc,\alpha_{q}\right)\quad\text{where}\ \alpha_{j}:=\alpha\Big{|}_{I_{j}}\in E_{j}\quad\text{for}\ 1\leq j\leq q.

Note that under this identification, $E$ has the equivalent representation as $E=E_{1}\times\dots\times E_{q}$ .

Definition 4.2.

Suppose $\mathbf{n}\in\mathbb{N}^{q}$ is the interaction type for an interaction function, $f\in\mathbb{R}\left[z_{1},\dotsc,z_{n}\right]$ as in Definition 1.1, and $E$ denotes the corresponding Boolean indices. For $\alpha\in E$ , define $p_{\alpha}\in\mathcal{P}\subset\mathbb{R}\left[\ell_{1},\dotsc,\ell_{n},\delta_{1},\dotsc,\delta_{n}\right]$ , by the formula

p_{\alpha}:=\prod_{j=1}^{q}\left(\sum_{k\in I_{j}}\ell_{k}+\alpha(k)\delta_{k}\right).

(15)

When convenient, we use a linear indexing scheme for elements of $\mathcal{P}$ which we define via the Boolean indexing map by identifying $p_{i}:=p_{\alpha}$ where $\alpha=B^{-1}(i)$ . To avoid confusion, we exclusively use Greek subscripts when referring to elements of $\mathcal{P}$ by their Boolean indices, and Latin subscripts when referring to elements of $\mathcal{P}$ by their linear indices. We leave it to the reader to check that the linearly indexed polynomials in Examples 1.3 and 1.4 are in agreement with that of Definition 4.2 via this identification.

Definition 4.3.

Let $\alpha,\beta\in E$ be a pair of Boolean indices corresponding to $\mathbf{n}\in\mathbb{N}^{q}$ . An ordered pair $(\alpha,\beta)$ satisfies the one bit condition if $\alpha(I_{j}(k))\leq\beta(I_{j}(k))$ , for all $1\leq j\leq q$ and $0\leq k\leq n_{j}-1$ , with equality for all but exactly one $(j,k)$ pair.

Remark 4.4.

Observe that if $(\alpha,\beta)$ satisfy the one bit condition and $(j_{0},k_{0})$ is the unique pair for which $\alpha$ and $\beta$ take different values, then $\alpha(I_{j_{0}}(k))=0$ and $\beta(I_{j_{0}}(k))=1$ .

Remark 4.5.

The one bit condition induces a poset structure on $E$ by setting $\alpha\prec\beta$ for each $(\alpha,\beta)$ satisfying the one bit condition, and extending the relation transitively. The one bit condition is a covering relation for the partial order.

Definition 4.6.

Suppose $\mathbf{n}\in\mathbb{N}^{q}$ is the interaction type for an interaction function, $f\in\mathbb{R}\left[z_{1},\dotsc,z_{n}\right]$ as in Definition 1.1, and $E$ denotes the corresponding Boolean indices. Let $\mathcal{P}$ be the set of polynomials indexed as in Definition 4.2. The PSD problem is defined by the triple, $(\mathcal{P},\prec,(0,\infty)^{2n})$ where $\prec$ is given by Remark 4.5.

The next proposition proves that $(\mathcal{P},\prec,(0,\infty)^{2n})$ satisfies Equation (1), and furthermore that $(\mathcal{P},\prec)$ is a Boolean partial order which justifies expressing the PSD problem as $\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ .

Proposition 4.7.

Consider a PSD problem $(\mathcal{P},\prec,(0,\infty)^{2n})$ . Then,

1.

$(\mathcal{P},\prec)$ is a Boolean lattice.
2.

For any $\alpha,\beta\in E$ , if $\alpha\prec\beta$ , then

$p_{\alpha}(\xi)<p_{\beta}(\xi)\quad\text{for all}\quad\xi\in(0,\infty)^{2n}.$

Proof.

To prove the first claim, let $S_{n}=\left\{1,\dotsc,n\right\}$ and let $(2^{S_{n}},\prec_{B})$ denote the standard Boolean lattice. Define a map, $\varphi:E\to 2^{S_{n}}$ , by the formula

\varphi(\alpha)=\left\{j\in S_{n}:\alpha(j)=1\right\},

and we note that $\varphi$ is a bijection since $E$ is defined to be the collection of all Boolean maps defined on $S_{n}$ . Furthermore, for any $\alpha,\beta\in E$ , we have by Definition 4.3 that $\alpha\prec\beta$ if and only if

\left\{j\in S_{n}:\alpha(j)=1\right\}\subset\left\{j\in S_{n}:\beta(j)=1\right\}

implying that $\varphi$ is an order isomorphism.

To establish the second claim, we must show that if $\alpha\prec\beta$ , then $p_{\alpha}(\xi)<p_{\beta}(\xi)$ holds for all $\xi\in(0,\infty)^{2n}$ . Note that by transitivity, it suffices to prove this holds for $(\alpha,\beta)$ satisfying the one bit condition. In this case we have

\beta(I_{j}(k))-\alpha(I_{j}(k))=\begin{cases}1&\text{if}\ j=j_{0}\text{ and }k=k_{0}\\ 0&\text{otherwise}.\end{cases}

for some $j_{0}\in\left\{1,\dotsc,q\right\}$ , $k_{0}\in\left\{0,\dotsc,n_{j_{0}-1}\right\}$ . If $\xi=\left(\ell_{1},\dotsc,\ell_{n},\delta_{1},\dotsc,\delta_{n}\right)\in\Xi$ , then from Equation (15) we have

	$\displaystyle p_{\beta}(\xi)$	$\displaystyle=\left(\left(\sum_{k\in I_{j_{0}}}\ell_{k}+\alpha(I_{j}(k))\delta_{k}\right)+\delta_{k_{0}}\right)\prod_{j\neq j_{0}}\sum_{k\in I_{j}}\ell_{k}+\alpha(I_{j}(k))\delta_{k}$
		$\displaystyle=p_{\alpha}(\xi)+\delta_{k_{0}}\prod_{j\neq j_{0}}\sum_{k\in I_{j}}\ell_{k}+\alpha(I_{j}(k))\delta_{k}$
		$\displaystyle>p_{\alpha}(\xi)$

as required. ∎

With Proposition 4.7 in mind, we return to writing $\prec_{B}$ in place of $\prec$ for the PSD problem where $\prec_{B}$ is the partial order of a Boolean lattice inherited by $\mathcal{P}$ from the one bit condition.

4.2 The linear PSD problem

We consider two cases of the PSD problem: the interaction type $\mathbf{n}\in\mathbb{N}^{q}$ for the function $f\in\mathbb{R}\left[z_{1},\dotsc,z_{n}\right]$ has the form $\mathbf{n}=(n)$ or $\mathbf{n}=(1,1,\ldots,1)$ . 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 a simple change of variables, we obtain an instance of LC-LEP with equivalent solutions since $\log$ is monotone, hence order preserving. We focus on the first case, leaving it to the reader to check that the second case is same modulo the evaluation domain ( $\mathbb{R}^{2n}$ versus $(0,\infty)^{2n})$ . Following the algorithm described in Section 3, we encode the partial order defined by $\prec_{B}$ as a set of linear constraints defined by a base cone which we must show is pointed. We begin by denoting the set of representation vectors for $\mathcal{P}$ as

\mathcal{V}:=\left\{\mathbf{u}_{p_{\alpha}}\in\left\{0,1\right\}^{2n}:\mathbf{u}_{p_{\alpha}}\ \text{is the representation vector of}\ p_{\alpha},\ \alpha\in E\right\}.

We define the set,

V_{\prec_{B}}:=\left\{\mathbf{u}_{p_{\beta}}-\mathbf{u}_{p_{\alpha}}:\mathbf{u}_{p_{\alpha}},\mathbf{u}_{p_{\beta}}\in\mathcal{V},\ (\alpha,\beta)\ \text{satisfies the one bit condition}\right\}.

(16)

which encodes the $\prec_{B}$ partial order into the representation vectors. These vectors will be the generators of the base cone for the algorithm in Section 3. Thus, we must show that $\operatorname{cone}(V_{\prec_{B}})$ generates a pointed cone.

Lemma 4.8.

Let $C_{0}:=\operatorname{cone}(V_{\prec_{B}})$ denote the cone generated by $V_{\prec_{B}}$ , then $C_{0}$ is pointed.

Proof.

By Proposition 3.2, $C_{0}$ is closed and convex so it suffices to prove that if $v\in C_{0}$ and $-v\in C_{0}$ , then $v=0$ . Fix $\xi\in(0,\infty)^{2n}$ and suppose $(\alpha,\beta)$ satisfies the one bit condition. By the formula in Equation (15) it follows that

p_{\beta}-p_{\alpha}=\delta_{i}

for some $i\in\left\{1,\dotsc,n\right\}$ . Since $\delta_{i}=\xi_{n+i}>0$ for all $\xi\in\Xi$ , it follows that

p_{\beta}(\xi)-p_{\alpha}(\xi)>0,

for every $(\alpha,\beta)$ satisfying the one bit condition. Passing to the representation vectors, it follows that for every $v\in V_{\prec_{B}}$ , we have $v\cdot\xi>0$ . Taking the conic hull, we have that if $v\in C_{0}\setminus\left\{0\right\}$ , then $v\cdot\xi>0$ . It follows that if $v,-v\in C_{0}$ simultaneously, then $v\cdot\xi\geq 0$ and $-v\cdot\xi\geq 0$ implying $v=0$ . ∎

4.3 The general PSD problem

Given a general PSD problem $\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ we present the construction of a LC-LEP denoted by $(\mathcal{P}^{\prime},\prec_{B},\mathbb{R}^{m})$ with the property that $\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)\subseteq\mathcal{T}(\mathcal{P}^{\prime},\prec_{B},\mathbb{R}^{m})$ . The importance of this is that $(\mathcal{P}^{\prime},\prec_{B},\mathbb{R}^{m})$ can be solved using Algorithm 3 and hence we obtain a rigorous upper bound on $\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ .

Definition 4.9.

Given an interaction type $\mathbf{n}\in\mathbb{N}^{q}$ , let $E=E_{1}\times\dots\times E_{q}$ denote the corresponding Boolean indices. Set $m:=\sum_{j=1}^{q}2^{n_{j}}$ and define the linearized evaluation domain to be

\mathbb{R}^{2^{n_{1}}}\times\dots\times\mathbb{R}^{2^{n_{q}}}\cong\mathbb{R}^{m}.

(17)

Define a polynomial ring in $m$ indeterminates with Boolean indexing by

\mathcal{R}:=\mathbb{R}\left[\left\{x_{j,\alpha_{j}}:\alpha_{j}\in E_{j},1\leq j\leq q\right\}\right],

(18)

and define a collection of linear polynomials by

\mathcal{P}^{\prime}:=\left\{p^{\prime}_{\alpha}:\alpha\in E\right\}\subset\mathcal{R}\qquad\text{where}\quad p^{\prime}_{\alpha}:=\sum_{j=1}^{q}x_{j,\alpha_{j}}.

The linearized PSD problem determined by $\mathbf{n}$ is to compute $\mathcal{T}\left(\mathcal{P}^{\prime},\prec_{B},\mathbb{R}^{m}\right)$ .

Theorem 4.10.

Fix an interaction type, $\mathbf{n}\in\mathbb{N}^{q}$ and let $\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ and $\mathcal{T}\left(\mathcal{P}^{\prime},\prec_{B},\mathbb{R}^{m}\right)$ denote the corresponding PSD and linearized PSD problems, respectively. The following are true.

1.

Let $\alpha,\beta\in E$ and $\xi\in(0,\infty)^{2n}$ . If $p_{\alpha}(\xi)<p_{\beta}(\xi)$ and

$\xi^{\prime}_{j,\alpha_{j}}=\log{\left(\sum_{k\in I_{j}}\xi_{k}+\alpha_{j}(k)\xi_{n+k}\right)}\in\mathbb{R}^{m},$

then $p^{\prime}_{\alpha}(\xi^{\prime})<p^{\prime}_{\beta}(\xi^{\prime})$ .
2.

$\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)\subseteq\mathcal{T}\left(\mathcal{P}^{\prime},\prec_{B},\mathbb{R}^{m}\right)$ .

Proof.

To prove the first claim, we define a map, $T:(0,\infty)^{2n}\to\mathbb{R}^{m}$ , by $\xi\mapsto\xi^{\prime}:=T(\xi)$ where the coordinates of $\xi^{\prime}$ are given by the formula

\xi^{\prime}_{j,\alpha_{j}}=\log{\left(\sum_{k\in I_{j}}\xi_{k}+\alpha_{j}(k)\xi_{n+k}\right)}.

(19)

Observe that $T$ is defined to satisfy the functional equation

\log\circ p_{\alpha}(\xi)=p^{\prime}_{\alpha}\circ T(\xi)\qquad\text{for all}\quad\alpha\in E,\ \xi\in(0,\infty)^{2n}.

(20)

Therefore, if $\alpha,\beta\in E$ and $\xi\in(0,\infty)^{2n}$ satisfies $p_{\alpha}(\xi)<p_{\beta}(\xi)$ , then $\log\left(p_{\alpha}(\xi)\right)<\log\left(p_{\beta}(\xi)\right)$ and it follows from Equation (20) that $p^{\prime}_{\alpha}(\xi^{\prime})<p^{\prime}_{\beta}(\xi^{\prime})$ where $\xi^{\prime}=T(\xi)$ as required.

To prove the second claim, consider $\mathcal{P}$ and $\mathcal{P}^{\prime}$ equipped with the linear indices as in Definition 4.2, and suppose $\sigma\in\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ . Then, by definition there exists $\xi\in(0,\infty)^{2n}$ satisfying

p_{\sigma(0)}(\xi)<p_{\sigma(1)}(\xi)<\dots<p_{\sigma(2^{n}-1)}(\xi).

Let $\xi^{\prime}=T(\xi)$ and apply the first result to successive pairs in the ordering which implies that for all $0\leq k\leq 2^{n}-2$ , we have

p^{\prime}_{\sigma(k)}(\xi^{\prime})=\log\left(p_{\sigma(k)}(\xi)\right)<\log\left(p_{\sigma(k+1)}(\xi)\right)=p^{\prime}_{\sigma(k+1)}(\xi^{\prime}).

Thus, we $\xi^{\prime}\in\mathbb{R}^{m}$ satisfies

p^{\prime}_{\sigma(0)}(\xi^{\prime})<p^{\prime}_{\sigma(1)}(\xi^{\prime})<\dots<p^{\prime}_{\sigma(2^{n}-1)}(\xi^{\prime}),

and it follows that $\sigma\in\mathcal{T}\left(\mathcal{P}^{\prime},\prec_{B},\mathbb{R}^{m}\right)$ which completes the proof. ∎

Example 4.11.

Recall the PSD in Example 1.3 with interaction function $f(z)=(z_{1}+z_{2})z_{3}$ corresponding to interaction type $\mathbf{n}=(2,1)$ . The corresponding polynomials for the linearized PSD problem are the following polynomials in $m=6$ variables

$\displaystyle p^{\prime}_{0}$	$\displaystyle=x_{1,(0,0)}+x_{2,0}$	$\displaystyle p^{\prime}_{4}$	$\displaystyle=x_{1,(0,1)}+x_{2,0}$
$\displaystyle p^{\prime}_{1}$	$\displaystyle=x_{1,(0,0)}+x_{2,1}$	$\displaystyle p^{\prime}_{5}$	$\displaystyle=x_{1,(0,1)}+x_{2,1}$
$\displaystyle p^{\prime}_{2}$	$\displaystyle=x_{1,(1,0)}+x_{2,0}$	$\displaystyle p^{\prime}_{6}$	$\displaystyle=x_{1,(1,1)}+x_{2,0}$
$\displaystyle p^{\prime}_{3}$	$\displaystyle=x_{1,(1,0)}+x_{2,1}$	$\displaystyle p^{\prime}_{7}$	$\displaystyle=x_{1,(1,1)}+x_{2,1}$

where we have used linear indexing for elements of $\mathcal{P}$ to match the polynomials in Example 1.3.

As defined in Equation (18), each variable of the form, $x_{1,(x,y)}$ , denotes an indeterminate in $\mathcal{R}$ which is identified with the element, $\alpha_{1}\in E_{1}$ , which satisfies $\alpha_{1}(1)=x$ and $\alpha_{1}(2)=y$ . In other words, the binary vector subscript, $(x,y)$ , denotes the two values which $\alpha_{1}\in E_{1}$ takes for the inputs in $I_{1}=\left\{1,2\right\}$ . Similarly, $x_{2,x}$ is identified with $\alpha_{2}\in E_{2}$ satisfying $\alpha_{2}(3)=x$ .

4.4 Solving the PSD problem for interaction type $\mathbf{n}=(2,1,\ldots,1)$ .

In this section we prove the following theorem.

Theorem 4.12.

Let $f$ be an interaction function with interaction type, $\mathbf{n}=(2,1,\ldots,1)$ . Let $\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ denote the corresponding PSD problem and $\left(\mathcal{P}^{\prime},\prec_{B},\mathbb{R}^{m}\right)$ the associated linearized PSD problem. Then $\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)=\mathcal{T}(\mathcal{P}^{\prime},\prec_{B},\Xi^{\prime})$ where $\Xi^{\prime}=\mathbb{R}^{m}\cap\{-\xi^{\prime}_{1,0}+\xi^{\prime}_{1,1}+\xi^{\prime}_{1,2}-\xi^{\prime}_{1,3}>0\}$ .

The proof of the theorem is based on the following lemma

Lemma 4.13.

Fix parameters, $x_{0},x_{1},x_{2},x_{3}\in\mathbb{R}$ , and define the function, $g:\mathbb{R}\to\mathbb{R}$ by the formula

g(t)=\exp(tx_{0})-\exp(tx_{1})-\exp(tx_{2})+\exp(tx_{3}).

If $x_{0}<x_{1}\leq x_{2}<x_{3}$ , then $g$ has a positive root if and only if $g^{\prime}(0)<0$ .

Proof.

Suppose first that $t_{0}$ is a root of $g$ . Expanding $\exp(t_{0}x_{1})$ and $\exp(t_{0}x_{2})$ to first order about $x_{0}$ and $x_{3}$ , respectively, and applying the mean value theorem yields the formula

g(t_{0})=-t_{0}\exp(t_{0}c_{1})(x_{1}-x_{0})-t_{0}\exp(t_{0}c_{2})(x_{2}-x_{3})=0

(21)

for some $c_{1}\in(x_{0},x_{1})$ and $c_{2}\in(x_{2},x_{3})$ . We define $k=c_{2}-c_{1}$ and multiply Equation (21) by $t_{0}e^{-kt_{0}}$ to obtain

e^{kt_{0}}(x_{3}-x_{2})-(x_{1}-x_{0})=0.

Noting that $c_{1}<x_{2}<c_{2}$ , it follows that $k>0$ . Therefore if $t_{0}>0$ , then $x_{3}-x_{2}<x_{1}-x_{0}$ or equivalently, $g^{\prime}(0)=x_{0}-x_{1}-x_{2}+x_{3}<0$ .

Conversely, if $g^{\prime}(0)<0$ then $g$ has at least one positive root since clearly $g(0)=0$ and $\lim\limits_{t\to\infty}g(t)=\infty$ . ∎

Proof of Theorem 4.12.

Suppose $\sigma\in\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ , $(0,\infty)^{2n}_{\sigma}$ is the associated realizable set as defined in Equation (2), and $\xi\in(0,\infty)^{2n}_{\sigma}$ , then by Theorem 4.10 we have $\xi^{\prime}=T(\xi)\in\mathbb{R}^{m}_{\sigma}$ . Note that by definition the first four coordinates of $\xi^{\prime}$ are given by the formulas

	$\displaystyle\xi^{\prime}_{1,0}$	$\displaystyle=\log(\xi_{1}+\xi_{2})$
	$\displaystyle\xi^{\prime}_{1,1}$	$\displaystyle=\log(\xi_{1}+\xi_{2}+\xi_{n+2})$
	$\displaystyle\xi^{\prime}_{1,2}$	$\displaystyle=\log(\xi_{1}+\xi_{2}+\xi_{n+1})$
	$\displaystyle\xi^{\prime}_{1,3}$	$\displaystyle=\log(\xi_{1}+\xi_{2}+\xi_{n+1}+\xi_{n+2}).$

Since $\xi_{i}>0$ for $i\in\left\{1,2,n+1,n+2\right\}$ , it follows that

-\xi^{\prime}_{1,0}+\xi^{\prime}_{1,1}+\xi^{\prime}_{1,2}-\xi^{\prime}_{1,3}>0,

so we have $\sigma\in\mathcal{T}(\mathcal{P}^{\prime},\prec_{B},\Xi^{\prime})$ .

Conversely, suppose $\sigma\in\mathcal{T}(\mathcal{P}^{\prime},\prec_{B},\Xi^{\prime})$ and $\xi^{\prime}\in\Xi^{\prime}_{\sigma}$ . From the Boolean lattice $\prec_{B}$ we have $\xi^{\prime}_{1,0}<\xi^{\prime}_{1,1}\leq\xi^{\prime}_{1,2}<\xi^{\prime}_{1,3}$ or $\xi^{\prime}_{1,0}<\xi^{\prime}_{1,2}\leq\xi^{\prime}_{1,1}<\xi^{\prime}_{1,3}$ . Moreover, $\xi^{\prime}$ also satisfies, $-\xi^{\prime}_{1,0}+\xi^{\prime}_{1,1}+\xi^{\prime}_{1,2}-\xi^{\prime}_{1,3}>0$ . Hence, Lemma 4.13 implies that there exists $t^{\prime}>0$ such that $\hat{\xi}^{\prime}:=t^{\prime}\xi^{\prime}$ satisfies

\exp(\hat{\xi}^{\prime}_{1,0})-\exp(\hat{\xi}^{\prime}_{1,1})-\exp(\hat{\xi}^{\prime}_{1,2})+\exp(\hat{\xi}^{\prime}_{1,3})=0.

Next, we define $\hat{\xi}\in(0,\infty)^{2n}$ by

\ \hat{\xi}_{j}=\begin{cases}\exp(\hat{\xi^{\prime}}_{j,0})&2<j\leq n\\ \exp(\hat{\xi^{\prime}}_{j,1})-\exp(\hat{\xi^{\prime}}_{j,0})&n+2<j<2n\\ \frac{1}{2}\exp(\hat{\xi^{\prime}}_{1,0})&j=1,2\\ \exp(\hat{\xi^{\prime}}_{1,2})-\exp(\hat{\xi^{\prime}}_{1,0})&j=n+1\\ \exp(\hat{\xi^{\prime}}_{1,1})-\exp(\hat{\xi^{\prime}}_{1,0})&j=n+2\\ \end{cases}

One easily verifies that $\hat{\xi}_{j}>0$ for all $1\leq j\leq 2n$ , and that $T(\hat{\xi})=\hat{\xi^{\prime}}$ . From Theorem 4.10, we have $\hat{\xi}\in(0,\infty)^{2n}_{\sigma}=\{\xi\in(0,\infty)^{2n}:p_{\sigma(0)}(\xi)<p_{\sigma(1)}(\xi)<\dots<p_{\sigma(2^{n}-1)}(\xi)\}$ which implies that $\sigma\in\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ . ∎

4.5 Solving the general PSD problem

In the general case, we have $\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)\subsetneq\mathcal{T}\left(\mathcal{P}^{\prime},\prec_{B},\mathbb{R}^{m}\right)$ , and thus, computing $\mathcal{T}\left(\mathcal{P}^{\prime},\prec_{B},\mathbb{R}^{m}\right)$ provides only a set of candidates for $\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ . This candidate set contains spurious linear 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, $\Xi^{\prime}\subsetneq\mathbb{R}^{m}$ , such that we still have the inclusion

\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)\subseteq\mathcal{T}(\mathcal{P}^{\prime},\prec_{B},\Xi^{\prime}).

(22)

Restricting to a smaller evaluation domain amounts to imposing more of the algebraic constraints a-priori which results in improved efficiency. In order for the candidate set on the right hand side to be efficiently computable using the algorithm in Section 3, it must be an instance of the LC-LEP i.e. $\Xi^{\prime}$ should be the interior of a polyhedral cone. For example, for the PSD with interaction type $\textbf{n}=(2,1,\dotsc,1)$ , analyzed in Section 4.4, we computed on the restricted domain

\Xi^{\prime}=\mathbb{R}^{m}\cap\left\{\xi^{\prime}\in\mathbb{R}^{m}:-\xi^{\prime}_{1,0}+\xi^{\prime}_{1,1}+\xi^{\prime}_{1,2}-\xi^{\prime}_{1,3}>0\right\}.

In terms of the algorithm in Section 3, this domain restriction amounts to taking our base cone in Algorithm 3 to be $\operatorname{cone}(V_{0})$ where

V_{0}=V_{\prec_{B}}\cup\left\{\mathbf{u}\right\}

and $\mathbf{u}$ is the representation vector for the linear functional defined by the formula

x\mapsto-x_{1,0}+x_{1,1}+x_{1,2}-x_{1,3}.

The requirement that this linear functional must be strictly positive is a special case of the following Lemma whose proof is a trivial computation.

Lemma 4.14.

Suppose $\alpha,\alpha^{\prime},\beta,\beta^{\prime}$ are Boolean indices such that for any $\xi\in(0,\infty)^{2n}$ , the following equations are satisfied.

	$\displaystyle p_{\alpha}(\xi)<p_{\beta}(\xi)<p_{\beta^{\prime}}(\xi)<p_{\alpha^{\prime}}(\xi)$
	$\displaystyle p_{\alpha}(\xi)+p_{\alpha^{\prime}}(\xi)=p_{\beta}(\xi)+p_{\beta^{\prime}}(\xi).$

Then,

\log(p_{\alpha}(\xi))+\log(p_{\alpha^{\prime}}(\xi))-\log(p_{\beta}(\xi))-\log(p_{\beta^{\prime}}(\xi))>0.

Lemma 4.14 provides a means to restrict the evaluation domain for the general linearized PSD problem as follows. Fix $j\in\left\{1,\dotsc,q\right\}$ and suppose $\left\{\alpha,\alpha^{\prime},\beta,\beta^{\prime}\right\}\subset E$ differ only in the $j^{\rm th}$ coordinate with $\alpha\prec_{B}\beta\prec_{B}\beta^{\prime}\prec_{B}\alpha^{\prime}$ , and also assume that $B(\alpha)+B(\alpha^{\prime})=B(\beta)+B(\beta)^{\prime}$ where $B$ is the Boolean indexing map. Then, it follows that for any $\xi\in\Xi$ , the values, $\left\{p_{\alpha}(\xi),p_{\alpha^{\prime}}(\xi),p_{\beta}(\xi),p_{\beta^{\prime}}(\xi)\right\}$ , satisfy both equations in Lemma 4.14. Therefore, if $\mathbf{u}(\left\{\alpha,\alpha^{\prime},\beta,\beta^{\prime}\right\})$ is the representation vector for the linear functional defined by

x\mapsto x_{j,B(\beta)}+x_{j,B(\beta^{\prime})}-x_{j,B(\alpha)}-x_{j,B(\alpha^{\prime})},

then $v(\left\{\alpha,\alpha^{\prime},\beta,\beta^{\prime}\right\})$ lies in $V_{\sigma}$ for any $\sigma\in\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ . Equivalently, we may impose the required linear constraint, $x_{j,B(\beta)}+x_{j,B(\beta^{\prime})}-x_{j,B(\alpha)}-x_{j,B(\alpha^{\prime})}>0$ on the evaluation domain of the linearized problem. Hence, for each $1\leq j\leq q$ , we define

V_{j}:=\left\{\mathbf{u}(\left\{\alpha,\alpha^{\prime},\beta,\beta^{\prime}\right\}):B(\alpha)+B(\alpha^{\prime})=B(\beta)+B(\beta)^{\prime},\ \alpha\prec_{B}\beta\prec_{B}\beta^{\prime}\prec_{B}\alpha^{\prime}\right\}

and for an arbitrary PSD problem, we may take our base cone to be

V_{0}=V_{\prec_{B}}\cup V_{\Xi}\qquad\text{where}\quad V_{\Xi}=\bigcup_{j=1}^{q}V_{j}.

Applying Algorithm 3 with the base cone generated by $V_{0}$ is equivalent to solving the instance of LC-LEP defined by $(\mathcal{P}^{\prime},\prec_{B},\Xi^{\prime})$ where $\Xi^{\prime}$ is the restriction of $\mathbb{R}^{m}$ to the subset for which the linear functionals defined by each $v\in V_{j}$ are strictly positive for each $1\leq j\leq q$ .

In addition to restricting the computation to the polyhedral cones discussed above, we can reuse solutions of smaller PSD problems in some larger computations. As an example, suppose $\mathcal{P}^{\prime}=\left\{p^{\prime}_{0},\dotsc,p^{\prime}_{7}\right\}$ is the set of interaction polynomials for the PSD with interaction type $\mathbf{n}^{\prime}=(2,1)$ and $\mathcal{P}:=\left\{p_{0},\dotsc,p_{15}\right\}$ the polynomials for the PSD problem with interaction type $\mathbf{n}=(2,1,1)$ . Observe that each admissible linear order on $\mathcal{P}^{\prime}$ induces an imposed linear order on the even indexed polynomials, $\mathcal{P}_{\operatorname{even}}:=\left\{p_{0},p_{2},\dots,p_{14}\right\}\subset\mathcal{P}$ . A similar linear order is induced on the odd indexed polynomials, $\mathcal{P}_{\operatorname{odd}}:=\left\{p_{1},p_{3},\dots,p_{15}\right\}\subset\mathcal{P}$ . Hence, a necessary condition to have an admissible linear extension for $\mathcal{P}$ is that the order of $\mathcal{P}_{\operatorname{even}}$ and $\mathcal{P}_{\operatorname{odd}}$ must both be consistent with one of the PSD solutions in $\mathcal{T}\left(\mathcal{P}^{\prime},\prec_{B},(0,\infty)^{6}\right)$ . This implies the inclusion

\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{8}\right)\subseteq\bigcup_{\sigma^{\prime}\in\mathcal{T}\left(\mathcal{P}^{\prime},\prec_{B},(0,\infty)^{6}\right)}\mathcal{T}(\mathcal{P},\prec_{B}\cup\prec_{\sigma^{\prime}},(0,\infty)^{8})

(23)

where $\prec_{B}\cup\prec_{\sigma^{\prime}}$ represents the refinement of the Boolean lattice partial order, and the partial order induced by $\sigma^{\prime}$ on the even/odd subsets.

To exploit this in general, we say that the PSD problem of type $\mathbf{{n}^{\prime}}$ is a sub-problem for the PSD problem of type $\mathbf{n}$ whenever the polynomials for $\mathbf{n}$ must obey an implied partial order determined by the solutions of $\mathbf{n^{\prime}}$ . Notice that the preceding discussion as well as Equation (23) applies also to an arbitrary polyhedral cone. Therefore, if there are a total of $k$ admissible linear extensions for all sub-problems of the PSD problem of type $\mathbf{n}$ which we have previously computed, then we bootstrap those results when computing $\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ via the inclusion

\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)\subseteq\bigcup_{i=1}^{k}\mathcal{T}(\mathcal{P},\prec_{B}\cup\prec_{\sigma_{i}^{\prime}},\Xi^{\prime})\subseteq\mathcal{T}(\mathcal{P},\prec_{B},\Xi^{\prime})

where $\prec_{B}\cup\prec_{\sigma^{\prime}_{i}}$ represents the refinement of the Boolean lattice partial order, and the partial order induced by $\sigma^{\prime}_{i}$ on the corresponding subsets obtained from any sub-problem. This technique has been used in the computation for all the cases of order $\geq 4$ . Observe that the computation of $\mathcal{T}(\mathcal{P},\prec_{B}\cup\prec_{\sigma_{i}^{\prime}},\Xi^{\prime})$ can be done distributively for $i=1,\dots,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 (22) is actually equality when $\Xi^{\prime}$ is constructed as we have described. However, in the typical case, these additional algebraic constraints are not sufficient to remove all spurious linear extensions except in the case $\textbf{n}=(2,1,\dots,1)$ . It remains an open problem to determine a smaller set $\Xi^{\prime}$ such that $\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)=\mathcal{T}\left(\mathcal{P}^{\prime},\prec_{B},\Xi^{\prime}\right)$ for other interaction types. However, in the remainder of this section we consider the problem of extracting $\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ from $\mathcal{T}\left(\mathcal{P}^{\prime},\prec_{B},\Xi^{\prime}\right)$ when they are not equal.

Observe that we may obtain large subsets of $\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ simply by sampling. The particular strategy that we adopted is as follows. We uniformly sampled between $10^{8}$ and $10^{9}$ points

\xi=(l_{1},\dots,l_{n},\delta_{1},\dots,\delta_{n})\in\mathbb{Z}_{+}^{2n}\cap B^{2n}_{\infty}(r),

where $B^{2n}_{\infty}(r)=\left\{\|\xi\|_{\infty}\leq r\right\}$ . We chose $r=1000$ . Mathematically the particular choice of $r$ is not important since the PSD polynomials are homogeneous, though in practice it does have an effect on sampling precision and speed. For each such $\xi$ we evaluated $\left\{p_{\alpha}(\xi):p\in\mathcal{P}\right\}$ . If $\sigma\in S_{2^{n}}$ denotes the linear order of these values, then $\xi$ serves as a “witness” for the claim that $\Xi_{\sigma}\neq\emptyset$ . This produces

\mathcal{S}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right):=\left\{\sigma\in\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right):\sigma\ \text{is witnessed by at least one sample}\right\}.

Obviously,

\mathcal{S}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)\subseteq\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)\subseteq\mathcal{T}(\mathcal{P},\prec_{B},\Xi^{\prime}).

In general, sampling is relatively efficient and in cases where $\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ is not too large (see Table 1 for details), we recover the entire solution.

Once we have constructed the set $\mathcal{S}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right),\mathcal{T}(\mathcal{P},\prec_{B},\Xi^{\prime})$ from sampling and algorithms in Section 3 respectively, we apply a CAD algorithm to check whether or not the semi-algebraic set,

\Xi_{\sigma}=\{\xi\in\Xi:p_{\sigma(0)}(\xi)<p_{\sigma(1)}(\xi)<\dots<p_{\sigma(2^{n}-1)}(\xi)\},

is empty for each $\sigma\in\mathcal{T}(\mathcal{P},\prec_{B},\Xi^{\prime})\setminus\mathcal{S}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ and then $\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ is recovered. The CAD algorithm implementation we are using is CylindricalAlgebraicDecomposition in Mathematica 11 [26].

4.6 Computational complexity

Analyzing the efficiency of our solver for the PSD problem can be broken into three cases. The first case is the linear PSD problem. This is an instance of the LC-LEP in $2n$ variables, and consequently, the complexity analysis in Section 3.3 applies.

The second case is the PSD problem with interaction type $\mathbf{n}=(2,1,\dotsc,1)$ . By Theorem 4.12, the solution set is identical to the solutions of the associated linearized PSD problem. The latter is an instance of LC-LEP in $m=4+2(q-2)=2(q+1)$ variables and since $q=n-1$ , we have $m=2n$ . In other words, the computational cost of solving this PSD problem is again equivalent to solving an instance of the LC-LEP in the same number of variables and the analysis in Section 3.3 applies.

The final case is the general PSD problem for interaction type $\mathbf{n}\in\mathbb{N}^{q}$ , which we assume does not satisfy either of the two previous cases. Here we cannot guarantee that the complexity is any better than that of a CAD-based approach since the final step in our algorithm is to determine the admissibility of every linear extension

\sigma\in\mathcal{S}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)\setminus\mathcal{T}(\mathcal{P},\prec_{B},\Xi^{\prime}).

That is, any total order which is admissible for the linearized PSD problem, but is not witnessed when sampling the nonlinear PSD problem, must be rigorously checked for admissibility. This is done using a CAD-based algorithm to certify whether or not the semi-algebraic set, $\Xi_{\sigma}$ , is empty for each such $\sigma$ . As discussed in Section 1, solving this problem requires, again in the worst case, computing a full CAD for $\mathcal{P}$ subject to the additional algebraic constraints $\ell_{i}>0,\delta_{i}>0$ for $1\leq i\leq n$ . Consequently, our algorithm for solving the PSD inherits the double exponential complexity of the CAD algorithm in Equation (3). Based on the actual computations we have carried out, we conjecture that this worst case bound is not typical and may not even be sharp. There are several heuristics which make this conjecture plausible.

The first is the fact that sampling alone was sufficient to recover the full solution of the PSD problem in every case examined in this work (compare the first and last columns of Table 1). Consequently, every candidate ordering which was not ruled out by sampling was in fact non-admissible and the CAD-based computation only needed to certify this face. This is advantageous due to the fact that many speedups for the general CAD algorithm apply specifically for proving that a semi-algebraic set is in fact empty. Once again, each of these improvements maintains the same worst case running time, but they often produce much faster typical running times.

The second heuristic is due to the structure of the PSD problem itself. Not only are the polynomials in $\mathcal{P}$ linear with respect to each variable, but $\mathcal{P}$ itself contains a “complete” set of polynomials arising from the interaction function which forms a sort of template. This imposes stricter algebraic constraints on the linearized PSD problem as a consequence of Lemma 4.14. Additionally, this structure causes the PSD problems to nest into one another neatly in the sense of Equation (23). These properties combine to produce sets of candidate solutions of the linearized PSD problem which are remarkably smaller than the candidates produced by solving the linearized PSD problem naively (compare the second and third columns in Table 1).

These heuristics combined with the results for the cases we have computed make a compelling case that this algorithm is already reasonably efficient, but this is by no means the last word on the subject. Further improvements to this approach via imposing additional algebraic constraints, improving sampling techniques, or applying the linearization technique to other classes of polynomials is the subject of ongoing research.

5 Results for some PSD problems

In this section we provide (see Table 1) the results of our computations for interaction functions of orders 4, 5, and 6. A slightly different approach was taken to compute orders 5 and 6, from that used for 4. This had to do with the machines being used, but highlights the flexibility of our method.

For interaction functions of order 4, we applied Algorithm 1 using a rational linear programming algorithm. In particular, we used the implementation MixedIntegerLinearProgram from SageMath 8 [50]. This implies that the output of Algorithm 3 is correct. Observe that interaction type $(4)$ is linear and type $(1,1,1,1)$ is log linear, and therefore Algorithm 3 produces $\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ . The fact that our output agrees with that of [33] suggests that our code is functioning as desired. To compute the interaction type $(2,1,1)$ we apply Algorithm 3 to obtain $\mathcal{T}\left(\mathcal{P}^{\prime},\prec_{B},\mathbb{R}^{m}\right)$ . By Theorem 4.12 this determines $\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ .

To solve the PSD problem from interaction types $(2,2)$ and $(3,1)$ requires that we make use of the strategy discussed in Section 4.5. Again, we use Algorithm 3 to obtain $\mathcal{T}\left(\mathcal{P}^{\prime},\prec_{B},\mathbb{R}^{m}\right)$ . By Theorem 4.10, $\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)\subset\mathcal{T}\left(\mathcal{P}^{\prime},\prec_{B},\mathbb{R}^{m}\right)$ . As indicated in Column 7 of Table 1, we chose $10^{8}$ samples from $(0,\infty)^{8}$ and identified $5344$ and $3084$ linear orders, respectively. We ran CylindricalAlgebraicDecomposition in Mathematica 11 [26] on each element of $\mathcal{T}\left(\mathcal{P}^{\prime},\prec_{B},\mathbb{R}^{m}\right)\setminus\mathcal{S}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ . As can be seen by comparing Columns 6 and 3, none of these elements were admissible.

$\mathbf{n}$	$\#\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$	$\#\mathcal{T}(\mathcal{P}^{\prime},\prec_{B},\Xi^{\prime})$	$\#\mathcal{T}\left(\mathcal{P}^{\prime},\prec_{B},\mathbb{R}^{m}\right)$	$\#\mathcal{S}(\mathcal{P}^{\prime},\prec_{B},(0,\infty)^{2n})$
(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)	-	11,035,808	*	6,570,952
(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,640

Table 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 the number of elements in an associated LC-LEP. This is not relevant where the AC-LEP problem of interest is a LC-LEP problem and is indicated by -. The * indicates that the computation was too large to complete. Column 4 provides the number of elements in the linearized PSD problem without 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 number of cells identified via sampling. We used

10^{8}

samples for all

n=4

cases and

10^{9}

samples for the

n=5,6

cases. The symbol ^† indicates that our sampling was not sufficient.

We now turn to the computations of interaction functions of order 5 and 6. As these problems are 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 package scipy [53] with the default numerical error $10^{-13}$ , in Algorithm 1. The interaction type type $(5)$ and $(6)$ are linear and $(1,1,1,1,1)$ and $(1,1,1,1,1,1)$ are log linear, and therefore via Algorithm 3 we obtain $\mathcal{T}_{alg}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ . We use the sampling technique (see Columns 7 and 8 of Table 1) to verify each of the elements of $\mathcal{T}_{alg}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ , thereby obtaining $\mathcal{T}\left(\mathcal{P},\prec_{B},(0,\infty)^{2n}\right)$ .

The computation for each order 4 case was done on a Mac Pro laptop (2.7 Hz Intel i5 and memory 8GB) with computation time under 4 hours. The computation of the remaining cases were done using a computing server with CentOs, intel 17.1, memory 32 GB, and less than 30 nodes. The computation time for both $(1,1,1,1,1)$ and $(5)$ was less than 4 hours, while the computation time for $(2,1,1,1)$ was on the order of 7 days. The codes which produced all of the computations in Table 1 are available on GitHub.

Acknowledgments

The authors would like to thank Shaun Harker, Sandra Di Rocco, Tomas Gedeon, and Mike Saks for helpful 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-0033 and FA8750-17-C- 0054, and NIH grant R01 GM126555-01.

References

[1] Reka Albert, James J. Collins, and Leon Glass. Introduction to Focus Issue: Quantitative approaches to genetic networks. Chaos, 23(2):025001, JUN 2013.
[2] David L Applegate, William Cook, Sanjeeb Dash, and Daniel G Espinoza. Exact solutions to linear programming problems. Operations Research Letters, 35(6):693–699, 2007.
[3] Zin Arai, William Kalies, Hiroshi Kokubu, Konstantin Mischaikow, Hiroe Oka, and Paweł Pilarczyk. A database schema for the analysis of global dynamics of multiparameter systems. SIAM J. Appl. Dyn. Syst., 8(3):757–789, 2009.
[4] Saugata Basu, Richard Pollack, and Marie-Françoise Roy. Algorithms in real algebraic geometry, volume 10 of Algorithms and Computation in Mathematics. Springer-Verlag, Berlin, 2003.
[5] Stephen Boyd and Lieven Vandenberghe. Convex Optimization, 2004.
[6] Graham Brightwell and Peter Winkler. Counting Linear Extensions is #P-complete. In Proceedings of the Twenty-third Annual ACM Symposium on Theory of Computing, STOC ’91, pages 175–181, New York, NY, USA, 1991. ACM.
[7] Graham R Brightwell. The Number of Linear Extensions of Ranked Posets. Technical report, 2003.
[8] Christopher W. Brown. Improved projection for cylindrical algebraic decomposition. Journal of Symbolic Computation, 32(5):447–465, 2001.
[9] 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.
[10] Christopher W. Brown and James H. Davenport. The complexity of quantifier elimination and cylindrical algebraic decomposition. Proceedings of the International Symposium on Symbolic and Algebraic Computation, ISSAC, pages 54–60, 2007.
[11] Justin Bush, Wes Cowan, Shaun Harker, and Konstantin Mischaikow. Conley-Morse databases for the angular dynamics of Newton’s method on the plane. SIAM J. Appl. Dyn. Syst., 15(2):736–766, 2016.
[12] Justin Bush, Marcio Gameiro, Shaun Harker, Hiroshi Kokubu, Konstantin Mischaikow, Ippei Obayashi, and PawełPilarczyk. Combinatorial-topological framework for the analysis of global dynamics. Chaos, 22(4):047508, 16, 2012.
[13] George E. Collins. Quantifier elimination for real closed fields by cylindrical algebraic decompostion. In H. Brakhage, editor, Automata Theory and Formal Languages, pages 134–183, Berlin, Heidelberg, 1975. Springer Berlin Heidelberg.
[14] Charles Conley. Isolated invariant sets and the Morse index, volume 38 of CBMS Regional Conference Series in Mathematics. American Mathematical Society, Providence, R.I., 1978.
[15] Charles C. Conley and Joel A. Smoller. The existence of heteroclinic orbits, and applications. pages 511–524. Lecture Notes in Phys., Vol. 38, 1975.
[16] Bree Cummins, Tomas Gedeon, Shaun Harker, and Konstantin Mischaikow. Model Rejection and Parameter Reduction via Time Series. SIAM J. Appl. Dyn. Syst., 17(2):1589–1616, 2018.
[17] Bree Cummins, Tomas Gedeon, Shaun Harker, Konstantin Mischaikow, and Kafung Mok. Combinatorial Representation of Parameter Space for Switching Networks. SIAM J. Appl. Dyn. Syst., 15(4):2176–2212, 2016.
[18] Sarah Day, Yasuaki Hiraoka, Konstantin Mischaikow, and Toshi Ogawa. Rigorous numerics for global dynamics: a study of the Swift-Hohenberg equation. SIAM J. Appl. Dyn. Syst., 4(1):1–31 (electronic), 2005.
[19] Sarah Day, Oliver Junge, and Konstantin Mischaikow. A rigorous numerical method for the global analysis of infinite-dimensional discrete dynamical systems. SIAM J. Appl. Dyn. Syst., 3(2):117–160 (electronic), 2004.
[20] Rocky Diegmiller, Lun Zhang, Marcio Gameiro, Justinn Barr6, Jasmin Imran Alsous, Paul Schedl, Stanislav Y. Shvartsman, and Konstantin Mischaikow. Mapping parameter spaces of biological switches. 2020.
[21] Guenter Ewald. Combinatorial convexity and algebraic geometry. GTM, 168, 1996.
[22] Terrence Fine and John Gill. The enumeration of comparative probability relations. Ann. Prob., 4:667–673, 1976.
[23] Marcio Gameiro, Tomas Gedeon, Shane Kepley, and Konstantin Mischaikow. Rational design of complex phenotype via network models. arXiv e-prints, page arXiv:2010.03803, October 2020.
[24] Tomas Gedeon, Bree Cummins, Shaun Harker, and Konstantin Mischaikow. Identifying robust hysteresis in networks. PLOS Computational Biology, 14(4):1–23, 04 2018.
[25] Shaun Harker, Konstantin Mischaikow, and Kelly Spendlove. A Computational Framework for the Connection Matrix Theory. arXiv e-prints, page arXiv:1810.04552, October 2018.
[26] Wolfram Research, Inc. Mathematica, Version 11.
[27] Abdul Salam Jarrah, Reinhard Laubenbacher, Brandilyn Stigler, and Michael Stillman. Reverse-engineering of polynomial dynamical systems. Advances in Applied Mathematics, 39(4):477 – 489, 2007.
[28] Tomasz Kaczynski, Konstantin Mischaikow, and Marian Mrozek. Computational homology, volume 157 of Applied Mathematical Sciences. Springer-Verlag, New York, 2004.
[29] William D. Kalies, Konstantin Mischaikow, and Robert C. A. M. van der Vorst. An algorithmic approach to chain recurrence. Found. Comput. Math., 5(4):409–449, 2005.
[30] William D. Kalies, Konstantin Mischaikow, and Robert C. A. M. van der Vorst. Lattice structures for attractors I. J. Comput. Dyn., 1(2):307–338, 2014.
[31] William D. Kalies, Konstantin Mischaikow, and Robert C. A. M. van der Vorst. Lattice structures for attractors II. Found. Comput. Math., 16(5):1151–1191, 2016.
[32] William D. Kalies, Konstantin Mischaikow, and Robert C. A. M. van der Vorst. Lattice Structures for Attractors III. arXiv e-prints, page arXiv:1911.09382, November 2019.
[33] Diane Maclagan. Boolean Term Orders and the Root System $B_{n}$ . Order, 15:279–295, 1999.
[34] Bree Cummins Marcio Gameiro, Shaun Harker. Dsgrn: Dynamic signatures generated by regulatory networks. https://github.com/marciogameiro/DSGRN, 2020.
[35] Scott McCallum. An improved projection operation for cylindrical algebraic decomposition. Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), 204 LNCS:277–278, 1985.
[36] Christopher McCord. Mappings and homological properties in the Conley index theory. Ergodic Theory Dynam. Systems, 8^∗(Charles Conley Memorial Issue):175–198, 1988.
[37] Christopher McCord and Konstantin Mischaikow. On the global dynamics of attractors for scalar delay equations. J. Amer. Math. Soc., 9(4):1095–1133, 1996.
[38] Christopher McCord, Konstantin Mischaikow, and Marian Mrozek. Zeta functions, periodic trajectories, and the Conley index. J. Differential Equations, 121(2):258–292, 1995.
[39] Konstantin Mischaikow. Global asymptotic dynamics of gradient-like bistable equations. SIAM J. Math. Anal., 26(5):1199–1224, 1995.
[40] Konstantin Mischaikow and Marian Mrozek. Chaos in the Lorenz equations: a computer-assisted proof. Bull. Amer. Math. Soc. (N.S.), 32(1):66–72, 1995.
[41] Konstantin Mischaikow and Marian Mrozek. Chaos in the {L}orenz equations: a computer-assisted proof. Bull. Amer. Math. Soc. (N.S.), 32(1):66–72, 1995.
[42] Konstantin Mischaikow and Marian Mrozek. Conley index. In Handbook of dynamical systems, Vol. 2, pages 393–460. North-Holland, Amsterdam, 2002.
[43] John Montgomery. Cohomology of isolated invariant sets under perturbation. J. Differential Equations, 13:257–299, 1973.
[44] OEIS Foundation Inc. The On-Line Encyclopedia of Integer Sequences, 2020. https://www.sagemath.org.
[45] Jichang Sha and D. J. Kleitman. The number of linear extensions of subset ordering. Discrete Mathematics, 63(2-3):271–278, 1987.
[46] Roman Srzednicki. On rest points of dynamical systems. Fund. Math., 126(1):69–81, 1985.
[47] 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.
[48] Herry Suprajitno and I Bin Mohd. Linear programming with interval arithmetic. Int. J. Contemp. Math. Sciences, 5(7):323–332, 2010.
[49] Andrzej Szymczak. The Conley index and symbolic dynamics. Topology, 35(2):287–299, 1996.
[50] The Sage Developers. SageMath, the Sage Mathematics Software System (Version 8). https://www.sagemath.org.
[51] Seinosuke Toda. PP is as hard as the polynomial-time hierarchy. SIAM Journal on Computing, 20(5):865–877, 1991.
[52] Robert J Vanderbei. Linear programming: foundations and extensions. Springer, 2020.
[53] Pauli Virtanen et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17:261–272, 2020.
[54] Ying Xin, Bree Cummins, and Tomas Gedeon. Multistability in the epithelial-mesenchymal transition network. BMC BIOINFORMATICS, 21(1), FEB 24 2020.

Computing linear extensions for polynomial posets subject to algebraic constraints

Abstract

1 Introduction

Quantifier elimination and real algebraic geometry

Computing linear extensions of Boolean lattices

This work

Definition 1.1.

Remark 1.2.

Example 1.3.

Example 1.4.

Related work

Organization of paper

2 Dynamic Signatures Generated by Regulatory Networks

2.1 Mathematical framework

2.2 Mathematical interpretation of regulatory networks

2.3 Extensions

3 Solving the LC-LEP

3.1 Cones

Definition 3.1.

Proposition 3.2.

Proposition 3.3.

Proof.

Proposition 3.4.

Proof.

Definition 3.5.

Proposition 3.6.

3.2 An algorithm for identifying 𝒯​(𝒫,≺,Ξ)\mathcal{T}(\mathcal{P},\prec,\Xi)

Definition 3.7.

Theorem 3.8.

Proof.

3.3 Complexity analysis

4 Solving the general PSD problem

4.1 The PSD as an instance of AC-LEP

Definition 4.1.

Definition 4.2.

Definition 4.3.

Remark 4.4.

Remark 4.5.

Definition 4.6.

Proposition 4.7.

Proof.

4.2 The linear PSD problem

Lemma 4.8.

Proof.

4.3 The general PSD problem

Definition 4.9.

Theorem 4.10.

Proof.

Example 4.11.

4.4 Solving the PSD problem for interaction type 𝐧=(2,1,…,1)\mathbf{n}=(2,1,\ldots,1).

Theorem 4.12.

Lemma 4.13.

Proof.

Proof of Theorem 4.12.

4.5 Solving the general PSD problem

Lemma 4.14.

4.6 Computational complexity

5 Results for some PSD problems

Acknowledgments

References

3.2 An algorithm for identifying $\mathcal{T}(\mathcal{P},\prec,\Xi)$

4.4 Solving the PSD problem for interaction type $\mathbf{n}=(2,1,\ldots,1)$ .