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