Dynamic Enumeration of All Mixed Cells

Discrete & Computational Geometry, Mar 2007

Tomohiko Mizutani, Akiko Takeda, Masakazu Kojima

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1007%2Fs00454-006-1300-9.pdf

Dynamic Enumeration of All Mixed Cells

Discrete Comput Geom Geometry Discrete & Computational Tomohiko Mizutani 0 Akiko Takeda 0 Masakazu Kojima 0 0 Department of Mathematical and Computing Sciences, Tokyo Institute of Technology , 2-12-1 Oh-Okayama, Meguro-ku, Tokyo 152-8552 , Japan The polyhedral homotopy method, which has been known as a powerful numerical method for computing all isolated zeros of a polynomial system, requires all mixed cells of the support of the system to construct a family of homotopy functions. The mixed cells are reformulated in terms of a linear inequality system with an additional combinatorial condition. An enumeration tree is constructed among a family of linear inequality systems induced from it such that every mixed cell corresponds to a unique feasible leaf node, and the depth-first search is applied to the enumeration tree for finding all the feasible leaf nodes. How to construct such an enumeration tree is crucial in computational efficiency. This paper proposes a dynamic construction of an enumeration tree, which branches each parent node into its child nodes so that the number of feasible child nodes is expected to be small; hence we can prune many subtrees which do not contain any mixed cell. Numerical results exhibit that the proposed dynamic construction of an enumeration tree works very efficiently for large scale polynomial systems; for example, it generated all mixed cells of the cyclic-15 problem for the first time in less than 16 hours. 1. Introduction for some nonempty finite subset Ai of Zn and some nonzero ci (a) ∈ C, a ∈ Ai . Here xa = x1a1 x2a2 · · · xnan for a = (a1, a2, . . . , a+n) ∈ Zn . The set Ai consists of mi elements, + and is called the support of fi (x). Also, A = (A1, A2, . . . , An) is called the support of f (x). In this paper we focus on a fully mixed polynomial system where all supports A1, A2, . . . , An are all distinct. See [ 9 ] and [ 14 ] for the semi-mixed case where some of them are identical to each other. Enumeration of all mixed cells of the support A of a polynomial system f (x), which is the subject of this paper, plays an essential role in the polyhedral homotopy method. Using the mixed cells, we construct a family of polyhedral homotopy functions between start systems, which are auxiliary polynomial systems whose zeros can be computed easily, and the target system f (x). Starting from zeros of every start system, we then trace all curves of zeros, so-called homotopy paths, of every polyhedral homotopy function. Bernshtein’s theorem [ 1 ] guarantees that the mixed volume, which is equal to the total volume of all mixed cells and coincides with the total number of solutions of the start systems, provides an upper bound for the number of isolated solutions of the target system f (x) in (C\{0})n. Especially if the target system is in general position, i.e., the coefficients of the target system are generic, this bound is tight and all homotopy paths reach isolated solutions of the target system. Some formulations were proposed for finding all mixed cells. Among others, the formulation of finding all mixed cells as a system of linear inequalities with a certain additional combinatorial condition [ 8 ], [ 9 ], [ 17 ], [ 21 ] (or a family of systems of linear inequalities) is more efficient in computational time and memory requirement than the geometric formulation used in [ 23 ]. This paper is founded on the former formulation. In this formulation, an enumeration tree is constructed among a family of systems of linear inequalities, which are induced from the system of linear inequalities with a combinatorial condition describing all mixed cells. The enumeration tree satisfies the following properties: (i) A leaf node describes a mixed cell if and only if it is feasible (or more precisely the system of linear inequalities attached to the leaf node is feasible). (ii) Each mixed cell is corresponding to a unique feasible leaf node. (iii) Each node different from leaf nodes is a common subsystem of its child nodes, so that if it is infeasible then so are all of its descendant nodes. (iv) The root node is an empty system, which is always feasible. We then apply an enumeration method for finding all feasible leaf nodes; if a node is determined to be infeasible then so is its descendant nodes; hence the subtree having the node as a root can be pruned because it does not contain any mixed cell. There are two important issues in the efficient implementation of the enumeration of all mixed cells which correspond to feasible leaf nodes of an enumeration tree satisfying properties (i)–(iv). One is how we check the feasibility of each node. For this purpose, [ 8 ], [ 9 ], [ 17 ] and [ 21 ] utilize a linear programming (LP) problem having a linear system of inequalities attached to each node as its constraint and check the feasibility of the linear system. Papers [ 8 ], [ 9 ] and [ 17 ] apply the primal simplex method to the LP problem while [ 21 ] employs the dual simplex method. If we take account of the effective use of information obtained at a node for its child nodes, the dual simplex method has an advantage. Specifically, we can easily choose a feasible solution of the dual of the child LP from an optimal solution of the dual of the parent LP. At least, application of the dual simplex method is popular in the field of optimization to deal with such a situation effectively [ 18 ]. This paper also applies the dual simplex method to an LP problem attached to each node to check its feasibility, which is described as the application of the primal simplex method to the dual of the LP problem. The other important issue is how we construct enumeration trees. Enumeration trees need to satisfy properties (i)–(iv) as we mentioned above. Specifically, the root node is fixed to be an empty system of linear inequalities by property (iv), and properties (i) and (ii) determine the collection of leaf nodes. One has much freedom in choosing and allocating systems of linear inequalities, which are induced from the system of linear inequalities with a combinatorial condition describing the mixed cells, for intermediate level nodes. In the existing works [ 8 ], [ 9 ], [ 17 ], [ 21 ], the structure of the enumeration trees is determined and fixed before enumerating mixed cells. In such a static construction of an enumeration tree, any information obtained at a node during the execution of enumeration is never utilized at all to branch the node into its child nodes because a branching rule is fixed in advance. For numerical efficiency, however, it is ideal to branch the node into its child nodes so that a larger portion of its child nodes are infeasible and are pruned. To pursue this idea, this paper proposes dynamic enumeration where branching at a node is carried out with the effective use of information which is obtained from the dual simplex method applied to some “child LP problems” at the node; hence the dual simplex method plays an essential role in this situation too. We note that dynamic enumeration is often utilized in the branch-and-bound method for integer programs [ 18 ]. Numerical results exhibit that our dynamic enumeration method works very efficiently for finding all mixed cells in comparison with existing static enumeration methods [ 7 ], [ 9 ], [ 10 ], [ 17 ], [ 21 ], [ 23 ]. For instance, our dynamic enumeration method generated all mixed cells of the cyclic-14 problem and computed the mixed volume, which had been the largest one in cyclic-n problems solved by the existing methods, in 1 hour 36 minutes, while MixedVol [ 9 ], [ 10 ], which is known as the fastest software among the existing ones, solves the same problem in 7 hours 14 minutes. Furthermore, our method computed the mixed volume for the cyclic-15 problem through finding all mixed cells of the problem in 15 hours 45 minutes. It appears that this is the first report for the mixed volume of this problem. As for noon-n and chandra-n problems, it is shown that the speedup ratio between the computational times of our method and MixedVol increases as the size of these problems becomes larger. This paper is organized as follows. In Section 2 we describe a procedure for the construction of an enumeration tree satisfying properties (i)–(iv), and then outline our dynamic enumeration algorithm. Section 3 is devoted to technical details of the algorithm. We first show an LP formulation for checking the feasibility of each node in an enumeration tree, and then discuss the size of the primal–dual pair of LP problems. We then present how we branch each parent node into its child nodes so that the number of feasible child nodes is expected to be small. This is a core part of the dynamic enumeration algorithm. In Section 4 we show numerical results for some benchmark polynomial systems. 2. An Outline of the Dynamic Enumeration Algorithm for Finding All Mixed Cells 2.1. The Mixed Volume and Mixed Cells of the Support Set Let N := {1, 2, . . . , n}. For every i ∈ N , let Qi = conv(Ai ) denote the convex hull of the support set Ai of the polynomial fi (x). For all positive number λ1, λ2, . . . , λn, we consider the n-dimensional volume of the Minkowski sum λ1 Q1 + λQ2 + · · · + λn Qn = {λ1q1 + λ2q2 + · · · + λnqn: qi ∈ Qi , i ∈ N }. This volume is known to be a homogeneous polynomial of degree n in λ1, λ2, . . . , λn. The mixed volume of A = (A1, A2, . . . , An) is defined as the coefficient of λ1λ2 · · · λn in the polynomial. Huber and Sturmfels [ 14 ] introduced a fine mixed subdivision of the Minkowski sum Q = Q1 + Q2 +· · ·+ Qn satisfying the properties described below to compute the mixed volume of A. Suppose that the subdivision of Q consists of n-dimensional polytopes R1, R2, . . . , Rs , which satisfy (i) dim(Rj ) = n for every j ∈ {1, 2, . . . , s}, (ii) sj=1 Rj = Q, (iii) Rj1 ∩ Rj2 is a face of both Rj1 and Rj2 for j1 = j2. In addition, in order for the collection of R1, R2, . . . , Rs to be a fine mixed subdivision of Q, each piece polytope Rj needs to be represented as the Minkowski sum conv(C1j ) + conv(C2j ) + · · · + conv(Cnj ) for some subset C = (C1j , C2j , . . . , Cnj ) of A = (A1, A2, . . . , An), such that (iv) for each j ∈ {1, 2, . . . , s}, conv(Cij ) (i ∈ N ) is a simplex of dimension #Cij − 1 and in=1 dim(conv(Cij )) = n. Let Rj = conv(C1j ) + conv(C2j ) + · · · + conv(Cnj ) be a polytope in a fine mixed subdivision of Q. Then we call it a mixed cell if dim(conv(C1j )) = dim(conv(C2j )) = · · · = dim(conv(Cnj )) = 1. For a fully mixed polynomial system, the set of all mixed cells in a fine mixed subdivision of Q provides the mixed volume of A as the summation of the volume of these cells. See [ 9 ], [ 16 ] and [ 17 ] for a detail description of the computation of the mixed volume via mixed cells. Furthermore, [ 14 ] presents how to construct a fine mixed subdivision of Q. We choose a real-valued function ωi : Ai → R, and call ω = (ω1, ω2, . . . , ωn) a lifting on A. A lifting ω lifts Ai to Aˆi = a ωi (a) : a ∈ Ai . the set of lower facets of Qˆ forms the fine mixed subdivision of Q. We use the notation Aˆ = (Aˆ1, Aˆ2, . . . , Aˆn) and Qˆ i = conv(Aˆi ), Qˆ = Qˆ 1 + · · · + Qˆ n. We consider the lifted polytope Qˆ ⊆ Rn+1. A face of Qˆ is said to be a lower facet if its ianrnaenrdnoomrmnaulmhbaesrafoproωsiit(ivae) slaostthcaotoωrdi(ian)aties.gFeonreerivce.rTyhien∈ thNe apnrdojeevcteiroynaR∈n+A1i , choose → Rn of Li and Li [ 17 ] developed an efficient algorithm for the enumeration of all lower facets, which correspond to mixed cells of the support set A, of the lifted polytope Qˆ via LP problems. All mixed cells of A are in one-to-one correspondence to the set C = ({a1, a1}, {a2, a2}, . . . , {an, an}) ∈ in=1 Ai such that the linear inequality system (1) where aˆi , αˆ = aˆi , αˆ , aˆi , αˆ ≤ aˆ, αˆ , ∀aˆ ∈ Aˆi \{aˆi , aˆi } (i ∈ N ), aˆ = a ωi (a) ∈ Aˆi and αˆ = is feasible. The feasibility check of this system can be represented using an LP formulation. 2.2. Enumeration Tree We now describe mixed cells in terms of systems of linear inequalities based on the Li–Li algorithm [ 17 ], and the basic idea of our dynamic enumeration method. For every L ⊆ N , define (L) = C = (C1, C2, . . . , Cn): Ci ⊆ Ai , #Ci = 2 (i ∈ L) , Cj = ∅ ( j ∈ L) The set serves as candidates of the nodes of enumeration trees. Specifically, ∅n ∈ (∅) = {∅n} is the root node, and (N ) ⊂ the leaf nodes. In general, the th level nodes of an enumeration tree are chosen from (L). For every C ∈ , let L(C) = {i ∈ N : Ci = ∅}. For every C ∈ and L L⊆⊆NN,#,Lw=e denote the vector consisting of Ci (i ∈ L) by CL = (Ci : i ∈ L). For every i ∈ N and every a ∈ Ai , let ωi (a) denote a random number chosen from some bounded interval of R. For every C ∈ with Ci = {api , aqi } (i ∈ L(C)), we consider a linear inequality system in a variable vector α ∈ Rn: I(C): api − aqi , α api − a, α = ωi (aqi ) − ωi (api ), ≤ ωi (a) − ωi (api ), (a ∈ Ai \{api , aqi }, i ∈ L(C)). We say that C ∈ is feasible when I(C) is feasible. Let ∗(L) = {C ∈ (L): C is feasible}. Then we know from (1) that ∗(N ) forms the set of all mixed cells. Note that a leaf node C ∈ (N ) is a mixed cell if and only if C is feasible, so properties (i) and (ii), which are described in the previous section, are satisfied. Thus, for the enumeration of all mixed cells, we need to find all elements in ∗(N ). For every C ∈ with a proper subset L(C) of N and every t ∈ N \L(C), define a set of child nodes of C by W (C, t ) = {C¯ ∈ (L(C) ∪ {t }): C¯ L(C) = CL(C)}. We can build an enumeration tree if we successively choose t ∈ N \L(C) at each node C of the tree starting from the root node C = ∅n with L(C) = ∅ (see Algorithm 2.2 for more details). In the static enumeration method employed in [ 8 ], [ 9 ], [ 17 ] and [ 21 ], we first choose a permutation of N or a one-to-one mapping π : N → N , and restrict nodes of the enumeration trees to C ∈ ({π(1), π(2), . . . , π( )}) with some ∈ N . Note that we have the set (N ) of leaf nodes if we take = n, and the root node is a feasible node ∅n ∈ (∅). First, the static enumeration method constructs the set W (∅n, π(1)) as child nodes of the root node ∅n ∈ (∅). Next, if a node C ∈ ({π(1), π(2), . . . , π( )}) for some 1 ≤ < n has been found to be feasible, the static enumeration generates W (C, π( +1)) as the set of child nodes of C. Thus the structure of the static enumeration tree is completely determined by a permutation π : N → N . In the enumeration method that we propose in this paper, an enumeration tree is constructed dynamically as the enumeration of nodes proceeds. To explain a procedure for the construction of such a tree, we define some notation. Let T = (V , E ) be a rEoo=ted tnr=ee0 Esuc.hV0thcaotrtrheespvoenrdtesxtosethteVroaontdneoddgee∅sne,taEndawreewderifittneenEa0s aVs a=n emnp=t0yVsetafnodr consistency with discussions below. The procedure for the construction of a tree T with allocating nodes dynamically is written as follows: Algorithm 2.1. (Construction of a Tree T ) Input: A support A = (A1, A2, . . . , An). Output: A tree T = (V = n=0 V , E = n =0 E ). V ← ∅n ( = 0, . . . , n), E ← ∅ ( = 0, . . . , n) and ← 0. while < n do for all C ∈ V do Choose t from N \L(C). V +1 ← V +1 ∪ W (C, t ) and E +1 ← E +1 ∪ {(C, C¯ ) ∈ V ×V +1: C¯ ∈ W (C, t )}. end for ← end while + 1. If two nodes C ∈ V and C¯ ∈ V +1 of a tree are joined with an edge, we say that C¯ is a child node of the parent node C. Any node on all paths from C to reachable leaf nodes, which are elements in Vn, is a descendant node of C. Each feasible leaf node corresponds to mixed cells. Given an input data A, Algorithm 2.1 produces various types of trees T depending on the choice of an index t from N \L(C). For instance, a static enumeration tree proposed in the existing algorithms [ 8 ], [ 9 ], [ 17 ], [ 21 ] is constructed if we fix a permutation π of N in advance and if we choose π( + 1) for any C ∈ V as the index t in Algorithm 2.1. As in the static enumeration, ∅n ∈ (∅) serves as the root node and (N ) as the set of leaf nodes. Suppose that a node C with some proper subset L(C) of N has been found to be feasible. Then we try to choose a t ∈ N \L(C) so that only a small portion of its child nodes W (C, t ) are expected to be feasible. The important issue here is how inexpensively we estimate the number of feasible child nodes in W (C, s) for all s ∈ N \L(C). For this purpose, we propose a simple technique of feasibility check which applies a criterion of unboundedness detection in the simplex method. We also utilize the relation table given in [ 9 ] to find some infeasible child nodes in W (C, s). 2.3. The Dynamic Enumeration Algorithm By deleting worthless nodes which do not contain any mixed cell, we can efficiently enumerate all feasible leaf nodes of a tree. Indeed, if a node is infeasible, all of its descendant nodes are infeasible, and thus we need not check the feasibility of the descendant nodes. Furthermore, we employ a depth-first order when applying the feasibility check to all nodes of an enumeration tree to save memory requirement during the execution of enumeration. Taking account of these factors, a depth-first search algorithm for the enumeration of all mixed cells is constructed. It is convenient to introduce the word “list” which will be used in this algorithm. For a finite set A, we denote list( A) as an ordered sequence of the elements in A, where the actual order is not relevant in our succeeding discussions but it is fixed. For a pair of list( A) and list(B), where A and B are finite sets, list( A)+list(B) stands for the list which is generated by connecting list(B) with list( A) by “stacking” list(B) on list( A); for example, if list( A) = (a, b, c) and list(B) = (d, e), then list( A) + list(B) = (a, b, c, d, e). The dynamic enumeration algorithm is outlined below. Va is a list of nodes, initialized as an empty set, and ν∗ denotes the total number of nodes generated by the algorithm. The total amount of work to generate all mixed cells by the algorithm is measured by ν∗. Algorithm 2.2. (The Dynamic Enumeration Algorithm) Input: A lifted support Aˆ = (Aˆ1, Aˆ2, . . . , Aˆn). Output: All mixed cells C ∈ ∗(N ) and the total number ν∗ of nodes generated. Recall that for each node C, the system of linear inequalities I(C) is solved to see whether C is feasible or infeasible. Since the efficiency of the algorithm depends on the size of Va, we utilize the one point test [ 8 ], [ 9 ], [ 17 ], [ 21 ] to decrease the number of elements to be added to Va. Suppose that C ∈ Va and t ∈ N \L(C). For every a ∈ At , the one point test checks the feasibility of the system of linear inequalities in α ∈ Rn, I(C, t, a): I(C), a − b, α ≤ ωt (b) − ωt (a) (b ∈ At \{a}). (2) If I(C, t, a) is infeasible, we can remove C¯ with C¯ t = {a, a } for any a ∈ At \{a} from W (C, t ) because I(C¯ ) is also infeasible. Therefore, as feasible node candidates we only consider W1(C, t ) = {C¯ ∈ (L(C) ∪ {t }): C¯ L = CL and C¯ t ⊆ At (C)}, where At (C) = {a ∈ At : I(C, t, a) is feasible}. After mt (= #At ) feasibility checks of linear inequality systems for constructing At (C), we have W1(C, t ) satisfying {C¯ ∈ W (C, t ): C¯ is feasible} ⊆ W1(C, t ) ⊆ W (C, t ). Thus we can replace (C) by Va ← Va + list(W1(C, t )). This one point test is known to be very effective to increase the computational efficiency of enumeration [ 8 ], [ 9 ], [ 17 ], [ 21 ]. Ideally we would like to choose a t ∈ N \L(C) at (B) so that the size of W1(C, t ) is the smallest among the sizes of W1(C, s) (s ∈ N \L(C)). Finding such a t ∈ N \L(C) exactly, however, is expensive because all W1(C, s) (s ∈ N \L(C)) need to be constructed. Therefore, we propose to replace W1(C, s) by another set which can be obtained easily. Utilizing a feasible solution xinit of I(C, t, a) which is generated from a solution of I(C), our method computes Wˆ 1(C, s, xinit) (s ∈ N \L(C)) satisfying W1(C, s) ⊆ Wˆ 1(C, s, xinit) ⊆ W (C, s) (s ∈ N \L(C)), and chooses a t ∈ N \L(C) such that the size of Wˆ 1(C, t, xinit) attains the minimum among the sizes of Wˆ 1(C, s, xinit) (s ∈ N \L(C)). We call this method the dynamic enumeration method. In Section 3.2, we explain how to generate the set Wˆ 1(C, s, xinit) (s ∈ N \L(C)). Also the relation table proposed in [ 9 ] can be used to find some infeasible child nodes in W (C, s). In our numerical experiments, the relation table is applied to remove infeasible child nodes from W (C, s) before Wˆ 1(C, s, xinit) is constructed. 3. Technical Details of the Algorithm 3.1. Formulation for Checking the Feasibility of a System of Linear Inequalities The feasibility check of C ∈ , conducted at (A) of Algorithm 2.2, can be formulated via an LP problem. Namely, we test the feasibility of the following problem in the vector α ∈ Rn of decision variables: P(C): max. γ, α s.t. I(C), where γ ∈ Rn is some fixed vector. For every C ∈ the dual problem is written as with Ci = {api , aqi } (i ∈ L(C)), D(C): min. (x; C) s.t. (x; C) = γ, xa ≥ 0 (a ∈ Ai \{api , aqi }), −∞ < xaqi < +∞ (i ∈ L(C)). Here, a vector of decision variables is given by the column vector x = (xa: a ∈ Ai \{api }, i ∈ L(C)) ∈ RdD , where dD := (mi − 1), (3) i∈L(C) and the symbols (x; C) and (x; C) are linear functions in x such that (x; C) = (x; C) = i∈L(C) a∈Ai \{api } i∈L(C) a∈Ai \{api } (ωi (a) − ωi (api )) xa and (api − a) xa. Any real vector γ in P(C) can be taken for the cost vector. Accordingly we set γ so that D(C) becomes feasible. Since this primal–dual pair satisfies the duality theorem, P(C) is feasible if and only if D(C) is bounded below, and P(C) is infeasible if and only if D(C) is unbounded. Therefore, to determine the feasibility of C, we need to see whether D(C) is bounded or not. Now we consider a formulation of an LP as min. c, x s.t. Gx = h, xi ≥ 0 (i ∈ I ), −∞ < xj < +∞ ( j ∈ J ), where a coefficient matrix G ∈ Rk×d , cost vector c ∈ Rd and constant vector h ∈ Rk are given, and x ∈ Rd is a vector of decision variables. These index sets I and J of decision variables satisfy I ∩ J = ∅ and I ∪ J = {1, 2, . . . , d}. Here, xi (i ∈ I ) and xj ( j ∈ J ) are called the nonnegative variable and a free variable, respectively. The primal–dual pair P(C) and D(C) can be transformed into (4) by introducing slack variables to the inequalities of P(C) and replacing the cost vector γ of P(C) by −γ. In consequence of these transformations, P(C) has dP variables and kP equalities such that dP = n + (mi − 2) and kP = On the other hand, D(C) has dD variables in (3) and kD equalities such that kD = n. Here dD is not greater than dP for any L(C) ⊆ N . Also kD is constant whereas kP is (4) monotonically increasing with respect to the cardinality of L(C). From the definition of a mixed cell of the support set, Ai (i ∈ N ) has at least two elements, i.e., mi ≥ 2. Therefore there exists L ⊆ N such that kP ≥ kD. Consequently for any L such that L ⊆ L ⊆ N , the number of constraints and that of variables in D(C) are not greater than those of P(C). Therefore, it is reasonable to observe whether D(C) is bounded for checking the feasibility of C. We formulate the one point test, stated in the previous section, via an LP problem. For every C ∈ , t ∈ N \L(C) and a ∈ At , checking the feasibility of (2) can be written as the following LP problem in the vector α ∈ Rn of decision variables: P1(C, t, a): max. γ, α s.t. I(C, t, a), where γ ∈ Rn is some fixed vector. As the one point test, we check the feasibility of this problem for all a ∈ At . The dual problem of P1(C, t, a) is given by D1(C, t, a): min. (x; C) + s.t. (x; C) + Using x ∈ RdD of (3) and the column vector y = (yb: b ∈ At \{apt }) ∈ R(mt −1), the vector of decision variables in this problem is represented as x y ∈ Rd¯D , where d¯D := Since we can choose γ such that D1(C, t, a) is feasible, this primal–dual pair satisfies the duality theorem. Similar to P(C) and D(C), we can say that the size of D1(C, t, a) is not larger than that of P1(C, t, a) for any C ∈ , t ∈ N \L(C) and a ∈ At . Therefore, we deal with D1(C, t, a) as the one point test, and check whether this problem is bounded or not. From the results of the one point test, we can generate the set W1(C, t ) which satisfies W1(C, t ) ⊆ W (C, t ). We refer to how to fix a right-hand constant vector γ on D(C) and D1(C, t, a). For C ∈ (L) with Ci = {api , aqi } (i ∈ L(C)) and a proper subset L(C) of N , we consider the problem D(C). Using the arbitrary nonnegative vector xˆ ∈ RdD , we compute and set this γˆ as a right-hand vector γ of the problem D(C) and D1(C, t, a) for some a ∈ At and t ∈ N \L(C). As a result, D(C) is feasible. We suppose that the problem D(C) is bounded, and denote an optimal solution of this problem by x∗ ∈ RdD . Then the vector γˆ = (xˆ; C) xinit = is feasible solution in D1(C, t, a) (a ∈ At and t ∈ N \L(C)) because D1(C, t, a) with fixed y = (yb: b ∈ At \{apt }) = 0 is equivalent to D(C). Also, we consider the problem D(C¯ ) (C¯ ∈ (L(C) ∪ {t })) with C¯ L = CL , and use γˆ as a right-hand vector γ of this problem. We easily see that this problem is feasible. Furthermore, an optimal solution of D1(C, t, a) is a feasible solution of D(C¯ ). Since the simplex method is suitable for solving many LP problems with a similar structure, we employ the method to solve problems arising from checking the feasibility of linear inequality systems. 3.2. How to Choose an Index t from N \L(C) The choice of a t ∈ N \L(C) at (B) of Algorithm 2.2 has a major effect on the computational efficiency of this algorithm. As stated in Section 2.3, we want to choose a t ∈ N \L(C) such that the size of W1(C, s) is the smallest among s ∈ N \L(C). However, this task is expensive in general because we need to check the feasibility of I(C, t, a) for every a ∈ At and t ∈ N \L(C) at (B) additionally. We employ a less expensive technique to choose a t ∈ N \L(C) so that an evaluation measure ν∗ of the efficiency of our dynamic enumeration method becomes smaller. To check the feasibility of C ∈ (L) with a proper subset L of N , we have solved the dual problem D(C), and in (B) we have an optimal solution x∗ ∈ RdD of D(C). As stated in the previous subsection, if we set xinit ∈ Rd¯D by (5) using x∗, the vector xinit is a feasible solution of D1(C, t, a) for any a ∈ At and t ∈ N \L(C). Because the structures of D1(C, t, a) and D(C) are similar to each other, we usually require only a few iterations to solve D1(C, t, a) by the simplex method when using xinit as an initial feasible solution. Thus we expect that xinit is incident to an unbounded direction if D1(C, t, a) is unbounded. Accordingly, instead of applying the simplex method to check the feasibility of P1(C, t, a) (a ∈ At and t ∈ N \L(C)) , we propose to test whether the feasible solution xinit of D1(C, t, a) has unbounded directions or not. At (B) we consider Wˆ 1(C, s, xinit) = {C¯ ∈ (L(C) ∪ {s}): C¯ L = CL and C¯ s ⊆ Aˆs (C, xinit)}, where Aˆt (C, xinit) = {a ∈ At : xinit of D1(C, t, a) has no unbounded direction} and we choose an index tˆ ∈ N \L(C) which attains min s∈N\L(C) #Wˆ 1(C, s, xinit). In general, this index tˆ does not coincide with the index t which achieves the minimum number of elements in W1(C, s) (s ∈ N \L(C)). We will observe from the numerical results, however, that the evaluation measure ν∗ is much smaller for our dynamic enumeration method than for the static enumeration method, and the total computational time for finding all mixed cells is reduced dramatically. We next explain how to compute elements in Wˆ 1(C, t, xinit) (t ∈ N \L(C)) more precisely, using an LP form of (4) instead of D1(C, t, a) for simplicity of notation. For (4), we assume that the number of variables d is not less than that of constraints k and the vmearttreixxxG∈=R(dgo1,ngth2,e.f.e.a,sigbdl)e hreagsifounllwrohwichracnokn.sIifstthsiosfptrwooblceomm(p4o)nisenfetsa:stihbelev,ethcteorreoefxbisatssica variables xB = G−B1h ∈ Rk and nonbasic variables xN = 0 ∈ Rd−k where GB ∈ Rk×k is a basic matrix. Note that xB = (xb1 , xb2 , . . . , xbk )T and xN = (xn1 , xn2 , . . . , xnd−k )T . In particular, the set of basic variables and nonbasic variables are called a basis and nonbasis. An adjacent vertex x˜ of x is represented as x˜ = x + θ d by using a nonnegative scalar θ and a direction vector d. Let D = {1, 2, . . . , d}, and we denote the set of basic indices B = {b1, b2, . . . , bk } ⊂ D. Note that (d − k) extreme rays extend from a vertex x and one direction d is chosen by fixing an index j ∈ D\B. The direction d is composed of two component vectors dB and dN such that dB = −G−B1gj ∈ Rk and dN = di = 1, i = j, di = 0, i ∈ D\(B ∪ { j }), ∈ Rd−k , where di represents a component of d. When we move from x to x˜, the cost change per unit θ is c, d . Using a component cB of a cost vector c which corresponds to a basis B, this amount can be written as c, d = cj − cTB G−B1gj (6) and called a reduced cost for j ∈ D. To obtain an adjacent vertex x˜ of x so that the value of a cost function decreases from x, we search for a direction d with j ∈ D such that its reduced cost is negative, and determine the step size θ ≥ 0 such that a new vertex x˜ = x + θ d satisfies its constraints; if components di of a direction vector d are nonnegative for all i ∈ B ∩ I where I is the index set of nonnegative variables in (4), this problem is unbounded and we say that x has an unbounded direction. Otherwise, we compute the largest θ allowed by constraints for variables. Now we provide the criteria for detecting that the feasible solution xinit of D1(C, t, a) (maat∈rixAGtaBnd∈ tR∈n×Nn \ofL(DC()C))hiasseaqnuaulntbootuhnedebdasdicirmecatitorinx. oNnoxtiicneit tohfaDtt1h(eC,otp,tiam) aflorbaasniyc a ∈ At and t ∈ N \L(C) because the number of constraints is equal to each other, and a feasible solution xinit can be represented as (5) using an optimal solution x∗ of D(C). Also, note that the vector (G−B1)T cB in (6) is an optimal solution of its dual problem when the duality theorem holds for this primal–dual pair. Let α∗ ∈ Rn be an optimal solution of P(C). For some apt ∈ At and t ∈ N \L(C), reduced costs on xinit of D1(C, t, apt ) with Ci = {api , aqi } (i ∈ L(C)) are written as ωi (b) − ωi (api ) − api − b, α∗ , for i ∈ L(C) ∪ {t } and b ∈ Ai \{api }. Since α∗ is an optimal solution of P(C) which has I(C) as a constraint, these reduced costs are nonnegative for every b ∈ Ai \{api } and i ∈ L(C). Consequently, if there are b ∈ At \{apt } such that (i) ωt (b) − ωt (apt ) − apt − b, α∗ < 0, (ii) all components of a vector −G−B1(apt − b), which corresponds to the nonnegative variables in the basis, are nonnegative, then we see that the feasible solution xinit of D1(C, t, apt ) for apt ∈ At has an unbounded direction. Conversely, if the feasible solution xinit of D1(C, t, apt ) has no such b ∈ At \{apt }, then apt is added to Aˆt (C, xinit). The problem D(C) with #L(C) = has free variables, and the optimal basis contains all free variables if this problem is bounded. Therefore, since the basis on xinit of D1(C, t, a) (a ∈ At and t ∈ N \L(C)) has free variables, it is enough to check (n − ) components of a vector −G−B1(apt − b) in (ii). 4. Numerical Results The proposed algorithm has been implemented and coded in C++ language. All numerical experiments were executed on a 2.4 GHz Opteron 850 with 8 GB memory, running Linux. First, we observe the total number of nodes ν∗ generated by the static and dynamic enumeration, described in Section 2.3, for the cyclic-n [ 2 ] and noon-n [ 20 ] problems. In the cyclic-n problem, one polynomial has two monomials and others have n monomials such as In the noon-n problem, all polynomials have (n + 1) monomials such as x1 + x2 + · · · + xn−1 + xn, x1x2 + x2x3 + · · · + xn−1xn + xn x1, x1x2x3 + x2x3x4 + · · · + xn−1xn x1 + xn x1x2, . . . x1x2 · · · xn − 1. For each polynomial system, we denote the support set of the i th polynomial from the top as Ai , and for every i ∈ N and a ∈ Ai choose ωi (a) randomly from some bounded interval of R. We set Aˆ = (Aˆ1, Aˆ2, . . . , Aˆn) as the input data of Algorithm 2.2. When the static enumeration method is employed, we choose a one-to-one mapping π : N → N such that π(i ) = i (i ∈ N ). Table 1 shows the total number of nodes ν∗ generated by the static and dynamic enumeration (abbreviated by “Static enum.” and “Dynamic enum.”, respectively) for the cyclic-n and noon-n problem, and the next row “Ratio” indicates the ratio between ν∗ given by these two methods. For these systems, this table reveals the efficiency of the dynamic enumeration method in comparison with the static one, and we can see that the ratio increases as the size of each system becomes larger. Next, in terms of the computational time, we compare our dynamic enumeration algorithm with some existing ones. The following software packages have been developed for enumerating all mixed cells of a polynomial equation system: MVLP [ 7 ], MixedVol [ 9 ], [ 10 ], PHCpack [ 23 ], PHoM [ 21 ] and mvol [ 17 ]. In particular, [ 9 ] and [ 10 ] report the superiority of MixedVol, which is coded in C++ language, in computational time over the other existing software packages. Therefore, we compare our algorithm with MixedVol for the benchmark systems: the economic-n [ 19 ], chandra-n [ 4 ] and katsura-n [ 3 ] problems in addition to the cyclic-n and noon-n problems. Notice that the katsuran problem consists of (n + 1) polynomials with (n + 1) variables, and the others n polynomials with n variables. Here, we have omitted the description of the economic-n, chandra-n and katsura-n problems, which are found on the web site [ 22 ]. We summarize these computational times in Table 2. The column “Mixed volume” presents the mixed volume of the support A of the system, and “Mixed cells” is the total number of mixed cells generated by our algorithm. The column “Speed-up ratio” indicates the ratio between the cpu time of our algorithm and MixedVol, and “—” means that the software is not applied to the corresponding system. The last column “CPU time/cells” is the cpu time of our algorithm divided by the total number of mixed cells generated by the algorithm. For all tested systems, “CPU time/cells” increases as their sizes become larger. The increasing rate is more than linear and varies depending on each system. From the column “Speed-up ratio”, we can see that our dynamic enumeration method improves the cpu time necessary for finding all mixed cells considerably, and solves large-scale polynomial systems such as the cyclic-15, noon-19,20,21, economic-20, chandra-20,21,22 and katsura-15,16 problems. This table shows the first numerical results for enumerating mixed cells of these large-scale problems. For finding all mixed cells of a polynomial system efficiently, the following two issues are essential: (1) how we construct an enumeration tree among a family of linear inequalities induced from the polynomial system, and also (2) how we check the feasibility of a linear inequality system. In this paper we proposed a depth-first search algorithm for a dynamic construction of an enumeration tree. At each iteration, the algorithm checks the feasibility of a linear inequality system by solving an LP problem with the effective use of information obtained at the previous iteration. Our numerical results show that the proposed algorithm considerably decreases the number of the LP problems to be solved, compared with the existing algorithms which utilize a static construction of an enumeration tree. In consequence, finding all mixed cells of large-scale polynomial systems becomes possible. Indeed, the proposed algorithm generated all mixed cells of the cyclic-15 problem for the first time in less than 16 hours. Enumeration of all mixed cells of a polynomial system, which is the subject of this paper, plays an essential role in the polyhedral homotopy method, known as a powerful numerical method for computing all isolated zeros of a polynomial system. We expect that the polyhedral homotopy method utilizing the proposed dynamic enumeration technique successfully solves large-scale polynomial systems which have not been processed. Acknowledgments References The author is grateful to the referees for many helpful comments. 1. D. N. Bernshtein , The number of roots of a system of equations, Funct . Anal. Appl . 9 , 183 - 185 ( 1975 ). 2. G. Bjo¨rk and R. Fro ¨berg, A faster way to count the solutions of inhomogeneous systems of algebraic equations , J. Symbolic Comput . 12 ( 3 ), 329 - 336 ( 1991 ). 3. W. Boege , R. Gebauer and H. Kredel , Some examples for solving systems of algebraic equations by calculating Groebner bases , J. Symbolic Comput . 2 , 83 - 98 ( 1986 ). 4. S. Chandrasekhar , Radiative Transfer, Dover, New York, 1960 . 5. V. Chva´tal, Linear Programming , Freeman, San Francisco, CA, 1983 . 6. Y. Dai , S. Kim and M. Kojima , Computing all nonsingular solutions of a cyclic-n polynomial using polyhedral homotopy continuation methods , J. Comput. Appl . Math. 152 , 83 - 97 ( 2003 ). 7. I. Z. Emiris and J. F. Canny , Efficient incremental algorithms for the sparse resultant and the mixed volume , J. Symbolic Comput . 20 , 117 - 149 ( 1995 ). Software available at http://cgi.di.uoa.gr/∼emiris/indexeng.html. 8. T. Gao and T. Y. Li , Mixed volume computation via linear programming , Taiwanese J. Math. 4 , 599 - 619 ( 2000 ). 9. T. Gao and T. Y. Li , Mixed volume computation for semi-mixed systems , Discrete Comput. Geom . 29 ( 2 ), 257 - 277 ( 2003 ). 10. T. Gao , T. Y. Li and M. Wu , MixedVol: a software package for mixed volume computation , ACM Trans. Math. Software 31 ( 4 ), ( 2005 ). Software available at http://www.csulb.edu/∼tgao/. 11. T. Gunji , S. Kim , M. Kojima , A. Takeda , K. Fujisawa and T. Mizutani , PHoM-a polyhedral homotopy continuation method , Computing 73 ( 1 ), 57 - 77 ( 2004 ). 12. T. Gunji , S. Kim , K. Fujisawa and M. Kojima , PHoMpara-parallel implementation of the polyhedral homotopy continuation method , Computing 77 ( 4 ), 387 - 411 ( 2006 ). 13. L. Gurvits and A. Samorodnitsky , A deterministic algorithm for approximating the mixed discriminant and mixed volume, and a combinatorial corollary , Discrete Comput. Geom . 27 ( 4 ), 531 - 550 ( 2002 ). 14. B. Huber and B. Sturmfels , A polyhedral method for solving sparse polynomial systems , Math. Comp. 64 , 1541 - 1555 ( 1995 ). 15. S. Kim and M. Kojima , Numerical stability of path tracing in polyhedral homotopy continuation methods , Computing 73 , 329 - 348 ( 2004 ). 16. T. Y. Li, Solving polynomial systems by polyhedral homotopies , Taiwanese J. Math. 3 , 251 - 279 ( 1999 ). 17. T. Y. Li and X. Li , Finding mixed cells in the mixed volume computation , Found. Comput. Math. 1 , 161 - 181 ( 2001 ). Software available at http://www.math.msu.edu/∼li/. 18. J. T. Linderoth and M. W. P. Savelsbergh , A computational study of branch and bound search strategies for mixed integer programming , INFORMS J. Comput . 11 , 173 - 187 ( 1999 ) 19. A. Morgan, Solving Polynomial Systems Using Continuation for Engineering and Scientific Problems , Prentice-Hall, Englewood Cliffs, NJ, 1987 . 20. V. W. Noonburg , A neural network modeled by an adaptive Lotka-Volterra system , SIAM J. Appl. Math . 49 , 1779 - 1792 ( 1989 ). 21. A. Takeda , M. Kojima and K. Fujisawa , Enumeration of all solutions of a combinatorial linear inequality system arising from the polyhedral homotopy continuation method , J. Oper. Res. Soc. Japan 45 , 64 - 82 ( 2002 ). Software available at http://www.is.titech.ac.jp/∼kojima/index.html. 22. J. Verschelde , The database of polynomial systems is on his web site : http://www.math.uic.edu/∼jan/. 23. J. Verschelde , Algorithm 795: PHCPACK: a general-purpose solver for polynomial systems by homotopy continuation , ACM Trans. Math. Software 25 , 251 - 276 ( 1999 ). Software available at http://www.math.uic.edu/∼jan/. 24. J. Verschelde , P. Verlinden and R. Cools , Homotopies exploiting Newton polytopes for solving sparse polynomial systems , SIAM J. Numer. Anal . 31 , 915 - 930 ( 1994 ).


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs00454-006-1300-9.pdf

Tomohiko Mizutani, Akiko Takeda, Masakazu Kojima. Dynamic Enumeration of All Mixed Cells, Discrete & Computational Geometry, 2007, 351-367, DOI: 10.1007/s00454-006-1300-9