POLYNOMIAL HOMOTOPIES FOR DENSE, SPARSE AND DETERMINANTAL SYSTEMS

JAN VERSCHELDE

Abstract. Numerical homotopy continuation methods for three classes of polynomial systems are presented. For a generic instance of the class, every path leads to a solution and the homotopy is optimal. The counting of the roots mirrors the resolution of a generic system that is used to start up the deformations. Software and applications are discussed.

AMS Subject Classification. 14N10, 14M15, 52A39, 52B20, 52B55, 65H10, 68Q40.

Keywords. polynomial system, numerical algebraic geometry, homotopy, continuation, deformation, path following, dense, sparse, determinantal, Bézout bound, Newton polytope, mixed volume, root count, enumerative geometry, numerical Schubert calculus.

Contents

1. Introduction 1 2. Three Classes of Polynomial Systems 2 3. The Principles of Polynomial Homotopy Continuation Methods 5 4. The Geometry of the Deformations 8 5. Root Counts and Start Systems 10 5.1. Dense Polynomials modeled by Highest Degrees 10 5.2. Mixed Subdivisions of Newton Polytopes to compute Mixed Volumes 13 5.3. Sparse Polynomial Systems solved by Polyhedral Homotopies 16 5.4. Determinantal Polynomials arising in Enumerative Geometry 18 6. Numerical Software for Solving Polynomial Systems 21 7. The Database of Applications 22 8. Closing Remarks and Open Problems 24 References 25

1. Introduction

Solving polynomial systems numerically means computing approximations to all isolated solutions. Homotopy continuation methods provide paths to approximate solutions. The idea is to break up the original system into simpler problems. To solve the original system, the solutions of the simpler systems are deformed into the solutions of the original problem.

Date: July 11, 1999.

Research at MSRI is supported in part by NSF grant DMS-9701755, benefited from a post-doctoral fellowship at MSRI and is also supported in part by NSF grant DMS-9804846 at MSU.

This paper presents optimal homotopies for three different classes of polynomial systems. Optimal means that for generic instances of the classes there are no diverging solution paths, whence the amount of computational work is linear in the number of solutions. In the next section we list the principal key words, definitions and main theorems for dense, sparse and determinantal polynomial systems. The proofs of these theorems follow from the correctness of the homotopies.

Path-following methods are standard numerical techniques ([3, 4, 5], [71], [123, 125]) to achieve global convergence when solving nonlinear systems. For polynomial systems we can reach all isolated solutions. In the third section we describe the paradigm of Cheater’s homotopy ([57], [59]) or coefficient-parameter polynomial continuation ([75], [76]). This paradigm allows to construct homotopies for which singularities only occur at the end of the paths. To deal with components of solutions we use an embedding method that leads to generic points on each component. This method is essential to numerical algebraic geometry [93].

From [48] we cite: “Algebraic geometry studies the delicate balance between the geometrically plausible and the algebraically possible”. By a choice of coordinates we set up an algebraic formulation for a geometric problem that is then solved by automatic computations. While this approach is extremely powerful, we might get trapped into tedious wasted computations after loosing the original geometric meaning of the problem. In section four we stress the geometric intuition of homotopy methods. Compactifications and homogeneous coordinates provide us the tools to generate the numerically most favorable representations for the solutions to our problem. In section five we arrive at the heart of modern homotopy methods where we outline specific algorithms to implement the root counts^{1}. The counting of the roots mirrors the resolution of a system in generic position that is used as starting point in the deformations.

Polyhedral methods occupy the central part of current research, as they are responsible for a computational breakthrough in numerical general-purpose solvers for polynomial systems. Section six is devoted to numerical software with an emphasis on the structure of the package PHC, developed by the author during the past decade. Another novel and exciting research development concerns the numerical Schubert calculus, which is one of the major new features in the second public release of PHC. The author has gathered more than one hundred polynomial systems that arose in various application fields. This collection serves as a test suite for software and a gallery to demonstrate the importance of polynomial systems to mathematical modelling. In section seven we sample some interesting cases.

The reference list contains a compilation of the most relevant technical contributions to polynomial homotopy continuation. Besides those we want to point at some other works in the literature that are of special interest. Some user-friendly introductions to algebraic geometry appeared in recent years: see [1], [27], [37], with computational aspects in [15] and [16]. As Newton polytopes have become extremely important, we recommend [130] and the handbook chapters [35]. See also [104] for the interplay between the combinatorics of polytopes and the (real) roots of polynomials. A recent survey that also covers polyhedral homotopies along with other polynomial continuation methods appeared in [64].

1The term “root count” was coined by Canny and Rojas [13] while introducing mixed volumes to computational algebraic geometry.

2. Three Classes of Polynomial Systems

The classification in Table 1 is inspired by [44]. The dense class is closest to the common algebraic description, whereas the determinantal systems arise in enumerative geometry.

system  model  theory  space
dense  highest degrees  Bézout  ℙ^{n}  projective
sparse  Newton polytopes Bernshteǐn (ℂ^{∗})^{n}  toric
determinantal localization posets Schubert G_{mr} Grassmannian

Table 1. Key words of the three classes of polynomial systems.

For the vector of unknowns x = (x_{1}, x_{2}, …, x_{n}) and exponents a = (a_{1}, a_{2}, …, a_{n}) ∈ ℕ^{n}, denote x^{a} = x_{1}^{a_{1}}x_{2}^{a_{2}} ⋯ x_{n}^{a_{n}}. A polynomial system P (x) = 0 is given by P = (p_{1}, p_{2}, …, p_{n}), a tuple of polynomials p_{i} ∈ ℂ[x], i = 1, 2, …, n.

The complexity of a dense polynomial p is measured by its degree d:

p(x) = ∑_{0≤a_{1}+a_{2}+···+a_{n}≤d} c_{a}x^{a},  d = deg(p),  (1)

where at least one monomial of degree d should have a nonzero coefficient. The total degree D of a dense system P is D = ∏_{i=1}^{n} deg(p_{i}).

Theorem 2.1. (Bézout [15]) The system P (x) = 0 has no more than D isolated solutions, counted with multiplicities.

Consider for example

P (x_{1}, x_{2}) = \{ \begin{matrix}x_{1}^{4} + x_{1}x_{2} + 1 = 0 \\ x_{1}^{3}x_{2} + x_{1}x_{2}^{2} + 1 = 0 \end{matrix} with total degree D = 4 × 4 = 16.  (2)

Although D = 16, this system has only eight solutions because of its sparse structure.

The support A of a sparse polynomial p collects all exponents of those monomials whose coefficients are nonzero. Since we allow negative exponents (a ∈ ℤ^{n}), we restrict x ∈ (ℂ^{∗})^{n}, ℂ^{∗} = ℂ \ {0}.

p(x) = ∑_{a∈A} c_{a}x^{a},  ∀a ∈ A : c_{a} ≠ 0, A ⊂ ℤ^{n}, #A < ∞.  (3)

The Newton polytope Q of p is the convex hull of the support A of p. We model the structure of a sparse system P by a tuple of Newton polytopes Q = (Q_{1}, Q_{2}, …, Q_{n}), spanned by A = (A_{1}, A_{2}, …, A_{n}), the so-called supports of P.

The volume of a positive linear combination of polytopes is a homogeneous polynomial in the multiplication factors. The coefficients are mixed volumes. For instance, for (Q_{1}, Q_{2}), we write:

2!vol_{2}(λ_{1}Q_{1} + λ_{2}Q_{2}) = V_{2}(Q_{1}, Q_{1})λ_{1}^{2} + 2 · V_{2}(Q_{1}, Q_{2})λ_{1}λ_{2} + V_{2}(Q_{2}, Q_{2})λ_{2}^{2},  (4)

normalizing V_{2}(Q, Q) = 2!vol_{2}(Q). For the Newton polytopes of the system (2): 2!vol_{2}(λ_{1}Q_{1}+ λ_{2}Q_{2}) = 4λ_{1}^{2} + 2 · 8λ_{1}λ_{2} + 5λ_{2}^{2}. To interpret this we look at Figure 1 and see that multiplying P_{1} and P_{2} respectively by λ_{1} and λ_{2} changes their areas respectively with λ_{1}^{2} and λ_{2}^{2}. The cells in the subdivision of Q_{1} + Q_{2} whose area is scaled by λ_{1}λ_{2} contribute to the mixed volume. So, for the example in (2), the root count is eight.

Figure 1. Newton polytopes Q_{1}, Q_{2}, a mixed subdivision of Q_{1} + Q_{2} with volumes.

Theorem 2.2. (Bernshteǐn [7]) A system P (x) = 0 with Newton polytopes Q has no more than V_{n}(Q) isolated solutions in (ℂ^{∗})^{n}, counted with multiplicities.

The mixed volume was nicknamed [13] as the BKK bound to honor Bernshteǐn [7], Kushnirenko [51], and Khovanskiǐ [49].

For the third class of polynomial systems we consider a matrix [C|X] where C ∈ ℂ^{(m+r)×m} and X ∈ ℂ^{(m+r)×r} respectively collect the coefficients and indeterminates. Laplace expansion of the maximal minors of [C|X] in m-by-m and r-by-r minors yields a determinantal polynomial

p(x) = ∑ ; \begin{matrix}I ∪ J = U \\ I ∩ J = ∅ \end{matrix} sign(I, J)C[I]X[J], U = \{1, 2, …, m + r\},  (5)

where the summation runs over all distinct choices I of m elements of U. The partition {I, J} of U defines the permutation U ↦ (I, J) with sign(I, J) its sign. The symbols C[I] and X[I] respectively represent coefficient minors and minors of indeterminates. Note that for more general intersection conditions, the matrices [C|X] are not necessarily square.

The vanishing of a polynomial as in (5) expresses the condition that the r-plane X meets a given m-plane nontrivially. The counting and finding of all figures that satisfy certain geometric conditions is the central theme of enumerative geometry. For example, consider the following.

Theorem 2.3. (Schubert [91]) Let m, r ≥ 2. In ℂ^{m+r} there are

d_{m,r} = \frac{1! 2! 3! ⋯ (r−2)! (r−1)! · (mr)!}{m! (m+1)! (m+2)! ⋯ (m+r−1)!}  (6)

r-planes that nontrivially meet mr given m-planes in general position.

This root count d_{m,r} is sharp compared to other root counts, see [98] and [114] for examples.

We can picture the simplest case, using the fact that 2-planes in ℂ^{4} represent lines in ℙ^{3}. In Figure 2 the positive real projective 3-space corresponds to the interior of the tetrahedron.

Figure 2: m = 2 = r.

In Figure 2, we see two thick lines meeting four given skew fine lines in a point. When not all input planes have the same dimension, but when the number of solutions is still finite, Pieri’s formula [82] provides a root count [96, 97]. In [44] the problem is solved in chains of nested subspaces, using a cellation of the Grassmannian G_{mr} of r-planes in ℂ^{m+r}. A localization poset models [46] the specialization of the solution r-plane when the input is specialized.

Algorithmic proofs for the above theorems consist in two steps. First we show how to construct a generic start system that has exactly as many regular solutions as the root count. Then we set up a homotopy for which all isolated solutions of any particular target system lie at the end of some solution path originating at some solution of the constructed start system.

3. The Principles of Polynomial Homotopy Continuation Methods

Homotopy continuation methods operate in two stages. Firstly, homotopy methods exploit the structure of P to find a root count and to construct a start system P^{(0)}(x) = 0 that has exactly as many regular solutions as the root count. This start system is embedded in the homotopy

H(x, t) = γ(1 − t)P^{(0)}(x) + tP (x) = 0, t ∈ [0, 1],  (7)

with γ ∈ ℂ a random number. In the second stage, as t moves from 0 to 1, numerical continuation methods trace the paths that originate at the solutions of the start system towards the solutions of the target system. The good properties we expect from a homotopy H(x, t) = 0 are (borrowed from [64]):

1. (triviality) The solutions for t = 0 are trivial to find.

2. (smoothness) No singularities along the solution paths occur.

3. (accessibility) All isolated solutions can be reached.

Continuation or path-following methods are standard numerical techniques ([3, 4, 5], [71], [123, 125]) to trace the solution paths defined by the homotopy using predictor-corrector methods. The smoothness property of complex polynomial homotopies implies that paths never turn back, so that during correction the parameter t stays fixed, which simplifies the set up of path trackers. A pseudo-code description of a path tracker is in Algorithm 3.1.

The predictor delivers at each step of the method a new value for the continuation parameter and predicts an approximate solution of the corresponding new system in the homotopy. Figure 3 shows two common predictor schemes. The predicted approximate solution is adjusted by applying Newton’s method as corrector. The third ingredient in path-following methods is the adaptive step size control. The step length is determined to enforce quadratic convergence in the corrector to avoid path crossing.

Algorithm 3.1. Following one solution path by an increment-and-fix predictor-corrector method with an adaptive step size control strategy.

Input: H(x, t), x^{∗} ∈ ℂ^{n}: H(x^{∗}, 0) = 0,  homotopy and start solution
ϵ > 0, max it, max steps.  accuracy and upper bounds
Output: x^{∗}, success if ||H(x^{∗}, 1)|| ≤ ϵ.  approximate solution if success
t := 0; k := 0;  initialization
h := max step size;  step length
old t := t; old x^{∗} := x^{∗}  back up values for t and x^{∗}
previous x^{∗} := x^{∗};  previous approximate solution
stop := false;  combines stopping criteria

while t < 1 and not stop loop

t := min(1, t + h);  secant predictor for t
x^{∗} := x^{∗} + h(x^{∗} − previous x^{∗});  secant predictor for x^{∗}
Newton(H(x, t), x^{∗}, ϵ, max it,success);  correct with Newton’s method
if success  step size control
then h := min(Expand(h), max step size);  enlarge step length
previous x^{∗} := old x^{∗};  go further along path
old t := t; old x^{∗} := x^{∗};  new back up values
else h := Shrink(h);  reduce step length
t := old t; x^{∗} := old x^{∗};  step back and try again

end if;

k := k + 1;  augment counter
stop := (h < min step size) or (k > max steps);  stopping criteria

end loop;

success := (||H(x^{∗}, 1)|| ≤ ϵ).  report success or failure

Following all paths can be done sequentially, one path at a time, or in parallel, with for each solution path the same sequence of values of the continuation parameter. The sequential path-following method has the advantage that the low overhead of communication [6] makes it very suitable to run on multi-processor environments. Note that the memory requirements are optimal.

Figure 3. The secant and tangent predictor with step length h.

To solve repeatedly a polynomial system with the same coefficient structure P (c, x) = 0, the homotopy (7) is applied with P^{(0)} = P (c^{0}, x) = 0 a system with random coefficients c^{0}. Solving P (c^{0}, x) = 0 is no longer trivial, so the name cheater’s homotopy [57] is appropriate. A similar idea appeared in [75, 76]. For coefficients given as functions of parameters, a refined version of cheater’s homotopy in [59] avoids repeated evaluation of those functions during path following:

H(x, t) = P ((1 − [t − t(1 − t)γ])c^{0} + (t − t(1 − t)γ)c, x) = 0, t ∈ [0, 1], γ ∈ ℂ. (8)

In [59] it is proven that with (8) all isolated solutions of P (c, x) = 0 can be reached and that singularities can only occur at the end of the paths.

Typically, when using a cheater’s homotopy, the computational effort spent towards the end of the paths often accounts for most of the work. The main numerical problem is then to distinguish irrelevant solutions at infinity from ill-conditioned but possibly meaningful solutions. End games [43], [77, 78, 79], [95] provide several procedures to approximate the winding number of a path. Recently, Zeuthen’s rule was applied in [50] to determine numerically the multiplicity of an isolated solution. Multi-precision facilities are useful for evaluation of residuals and root refinement for badly scaled solutions.

In most applications, the polynomial systems have real coefficients and invite the use of real homotopies. In [11] it was conjectured and proven in [60] that generically, real homotopies contain no singular points other than a finite number of quadratic turning points. At those bifurcation points pairs of real solution paths become imaginary or conversely, complex conjugated solution paths join to yield two real solution paths. We refer to [2], [38], [64] and [60, 61] for a discussion of numerical techniques to deal with quadratic turning points. A remarkable application of real homotopies in the real world consists in the finding of the relevant parameters of a polynomial system to maximize the number of real roots, see [18] for the 40 real solutions for the Stewart-Gough platform in mechanics.

In [93] the use of homotopy continuation to deal with overdetermined and components of solutions is discussed. Geometrically one slices the components of solutions with as many random hyperplanes as the dimension of the components. The solutions to the original polynomial system augmented with these random linear equations for the hyperplanes are generic points of the components, constituting the main numerical data to study those components. In particular, the number of generic points one obtains by this slicing procedure equals the sum of the degrees over all top-dimensional components of solutions.

To make the algorithms of [93] more efficient, in [94], the following embedding of the polynomial system P (x) = 0 is proposed:

\{ \begin{matrix}p_{i}(x) + λ_{i}z = 0, & i = 1, 2, …, n \\ ∑_{j=1}^{n} c_{j}x_{j} + z = 0 \end{matrix}  (9)

where the λ_{i}’s and c_{j}’s are random complex numbers. This embedding has the advantage over the algorithms in [93] that fewer solution paths diverge. Solutions to the system (9) with z = 0 lie on a component of solutions. By Bertini’s theorem, all solutions with z ≠ 0 are regular. In [94], it is proven that those solutions can be used as start solutions to reach all isolated solutions of the original polynomial system P (x) = 0.

The embedding (9) is performed repeatedly in the routine ‘Embed’ in the algorithm (copied from [94]) below.

Algorithm 3.2. Cascade of homotopies between embedded systems.

Input: P, n.  system with solutions in ℂ^{n}
Output: (E_{i}, X_{i}, Z_{i})_{i=0}^{n}.  embeddings with solutions
E_{0} := P;  initialize embedding sequence
for i from 1 up to n do  slice and embed
E_{i} := Embed(E_{i−1}, z_{i});  z_{i} = new added variable
end for;  homotopy sequence starts
Z_{n} := Solve(E_{n});  all roots are isolated, nonsingular, with z_{n} ≠ 0
for i from n − 1 down to 0 do  countdown of dimensions
H_{i+1} := tE_{i+1} + (1 − t) ( \begin{matrix}E_{i} \\ z_{i+1} \end{matrix} );

homotopy continuation t : 1 → 0 to remove z_{i+1}

X_{i} := limits of solutions of H_{i+1}
as t → 0 with z_{i} = 0;  on component
Z_{i} := H_{i+1}(x, z_{i} ≠ 0, t = 0);  not on component: these solutions

are isolated and nonsingular

end for.

This embedding allows the efficient treatment of overdetermined systems and other nonproper intersections. By perturbing the added hyperplanes and extending the generic points by continuation, interpolation methods can lead to equations for the components.

4. The Geometry of the Deformations

Homotopy methods have an intuitive geometric interpretation. In this section we illustrate the geometry of the three types of moving into special position: product, toric, and Pieri deformations. These can be regarded as three applications of the principle of continuity or conservation of number in enumerative geometry.

Product homotopies deform polynomial equations into products of linear equations. In Figure 4 we see the line configuration at the start and the ellipse-parabola intersection in the end. Note that complex space is the natural space for deformations. The other two complex conjugated intersection points could not be displayed in Figure 4.

The sparser a system, the easier it can be solved. In Figure 5 we illustrate the idea of making a system sparser by setting up a so-called polyhedral homotopy that reduces this particular system at t = 0 to a linear system. The lower hull of the Newton polytope of this homotopy induces a triangulation, which is used to count the roots. In particular, every cell in the triangulation gives rise to a homotopy with as many paths to follow as the volume of the cell. The other root for the example in Figure 5 can be computed with a homotopy obtained from P̂ by the substitution of variables x_{1} ← x̃_{1}t^{−1} and x_{2} ← x̃_{2}t^{−1}. This transformation pushes the constant monomial up, so that at t = 0 we have the nonconstant monomials in the start system to compute the other root.

Figure 6 displays a special and a general configuration of four lines. The basis has been chosen such that two of the four input lines are spanned by standard basis vectors. To Figure 4. Intersection of quadrics: a degenerate and a target configuration.

Figure 5. Triangulation of the Newton polytope of P with polyhedral homotopy P̂.

compute all lines that meet four given lines, one of the four given lines is moved into special position so that it intersects two other given lines, see the left of Figure 6. The solution lines must then originate at those two intersection points and reach to the other opposite line while meeting the line left in general position.

The constructions above are in a sense [1] “heuristic proofs”. With the general position assumption we cheat a bit, avoiding the hard problem of assigning multiplicities. Making this so-called [127] “method of degeneration” rigorous was an important development in algebraic geometry.

To deal with solution paths diverging to ill-conditioned roots or to infinity we need to compactify our space. Instead of polynomials in n variables we consider homogeneous forms with coordinates subject to equivalence relations. While mathematically all coordinate choices are equivalent, we select the numerically most favorable representations of the solutions.

The usual projective transformation consists in the change of variables x_{i} := \frac{z_{i}}{z_{0}} , for i = 1, 2, …, n, which leads to the homogeneous system P (z) = 0. To have as many equations as unknowns, we add to this system a random hyperplane. Except for an algebraic set of Figure 6. In ℙ^{3} two thick lines meet four given lines L_{1}, L_{2}, L_{3}, and L_{4} in a point. At the left we see a special configuration and the general configuration is at the right.

the coefficients of this added hyperplane, all solution paths are guaranteed to stay inside ℂ^{n+1} when homotopy continuation is applied. We refer to [64] for numerical techniques that dynamically restrict the computations to n dimensions.

For sparse polynomial systems, we introduce as in [116] a new variable for every facet of the Newton polytopes. The advantage of this more involved compactification is based on the observation that when paths diverge to infinity certain coefficients of the polynomial system become dominant. With toric homogenization the added variables that become zero identify the faces of the Newton polytopes for the parts of the system that become dominant. This compactification works in conjunction with polyhedral end games [43] which are summarized in Section 5.3.

The natural way to compactify G_{mr} is to consider a multi-projective homogenization according to rows or columns of the matrix representations for the planes. In addition, we have that the planes are equivalent upon a linear change of basis. Choosing orthonormal matrices to represent the input planes leads to drastic improvements in the conditioning of the solution paths, see [46] and [114] for experimental data.

5. Root Counts and Start Systems

The main principle is that counting roots corresponds to solving start systems. Algorithms to illustrate this principle will be shown for little examples for the three classes of polynomial systems.

For dense polynomial systems, the computation of generalized permanents model the resolution of linear-product start systems. The algorithms to compute mixed volumes lead to polyhedral homotopies to solve sparse polynomial systems. The localization posets describe the structure of the cellation of the Grassmannian used to set up the Pieri deformations.

5.1. Dense Polynomials modeled by Highest Degrees. A polynomial in one variable has as many complex solutions as its degree. A linear system has either infinitely many solutions or exactly one isolated solution in projective space. By this analogy [53] we see that Bézout’s theorem generalizes these last two statements: a polynomial system has either infinitely many solutions or exactly as many isolated solutions in complex projective space as the total degree.

As the above presentation of Bézout’s theorem suggests, the simplest cases are univariate and linear systems, which are used as start systems. For the example system (2), a start system P^{(0)}(x) = 0 based on the total degree D is given by two univariate quartic equations x_{1}^{4} − c_{1} = 0 = x_{2}^{4} − c_{2}, where c_{1} and c_{2} are randomly chosen complex numbers. Note that the computation of D = 4 × 4 models the structure of the solutions of P^{(0)} as four solutions for x_{1} crossed with four solutions for x_{2}.

The earliest approaches of this homotopy appear in [14], [19], [31], [32], and were further developed in [52], [70], [129]. The book [71] contains a very good introduction to the practice of solving polynomial system by homotopy continuation. Regularity results can be found in [54] and [131]. While this homotopy algorithm has a sound theoretical basis, the total degree is a too crude upper bound for the number of affine roots to be useful in most applications.

Multi-homogeneous homotopies were introduced in [72, 73] and applied to various problems in mechanism design, see e.g. [118, 119]. Similar are the random product homotopies [55, 56], applying intersection theory in [58], but less suitable for automatic procedures. For our running example (2), we follow the approach of multi-homogenization and we group the unknowns according to the partition Z = {{x_{1}}, {x_{2}}}. The corresponding degree matrix M_{Z} has in its (i, j)-th entry the degree of the i-th polynomial in the variables of the j-th set of Z. The 2-homogeneous Bézout bound B_{Z} is the permanent of M_{Z}.

\begin{matrix}P (x_{1}, x_{2}) = \\ \{ & \begin{matrix}x_{1}^{4} + x_{1}x_{2} + 1 = 0 \\ x_{1}^{3}x_{2} + x_{1}x_{2}^{2} + 1 = 0 \end{matrix} \end{matrix}  \begin{aligned}M_{Z} & = & M_{\{\{x_{1}\},\{x_{2}\}\}} \\ & = & [ \begin{matrix}4 & 1 \\ 3 & 2 \end{matrix} ] \end{aligned}  \begin{aligned}B_{Z} & = & per(M_{Z}) \\ & = & 4 × 2 + 3 × 1 \\ & = & 11 \end{aligned}  (10)

The computation of the permanent follows the expansion for the determinant, except for the permanence of signs, as it corresponds to adding up the roots when solving the corresponding linear-product start system:

P^{(0)}(x) = \{ \begin{matrix}∏_{i=1}^{4}(x_{1} − α_{1i}) ∏_{i=1}^{1}(x_{2} − β_{1i}) = 0 \\ ∏_{i=1}^{3}(x_{1} − α_{2i}) ∏_{i=1}^{2}(x_{2} − β_{2i}) = 0 \end{matrix}  (11)

In most applications the grouping of variables follows from their meaning, e.g.: for eigenvalue problems Ax = λx, Z = {{λ}, {x_{1}, x_{2}, …, x_{n}}}. Efficient permanent evaluations in conjunction with exhaustive searching algorithms for finding an optimal grouping were developed in [120]. In case the number of independent roots equals the Bézout bound, interpolation methods [105] are useful.

Partitioned linear-product start systems were developed in [109] elaborating the idea that several different partitions can be used for the polynomials in the system. Motivated by symmetric applications [107], general linear-product start system were proposed in [106]. These start systems are based on a supporting set structure S which provides a more refined model of the degree structure of a polynomial system.

S = \begin{matrix}\{x_{1}\}\{x_{1}, x_{2}\}\{x_{1}\}\{x_{1}\} \\ \{x_{1}\}\{x_{1}, x_{2}\}\{x_{1}\}\{x_{2}\} \end{matrix}  (12)

To compute the bound formally, one collects all admissible n-tuples of sets, picking one set out of every row in the set structure.

\begin{aligned}B_{S} = & #\{(\{x_{1}\}, \{x_{1}, x_{2}\}), (\{x_{1}\}, \{x_{2}\}), (\{x_{1}, x_{2}\}, \{x_{1}\}), \\ & (\{x_{1}, x_{2}\}, \{x_{1}, x_{2}\}), (\{x_{1}, x_{2}\}, \{x_{1}\}), (\{x_{1}, x_{2}\}, \{x_{2}\}), \\ & (\{x_{1}\}, \{x_{2}\}), (\{x_{1}, x_{2}\}, \{x_{2}\})(\{x_{1}\}, \{x_{2}\})(\{x_{1}\}, \{x_{2}\})\} \end{aligned}  (13)

Each admissible pair corresponds to a linear system that leads to a solution of a generic start system:

P^{(0)}(x) = \{ \begin{matrix}(x_{1} + c_{11})(x_{1} + c_{12}x_{2} + c_{13})(x_{1} + c_{14})(x_{1} + c_{15}) = 0 \\ (x_{1} + c_{21})(x_{1} + c_{22}x_{2} + c_{23})(x_{1} + c_{24})(x_{2} + c_{25}) = 0 \end{matrix}  (14)

This start system has B_{S} = 10 solutions. In [106], the following theorems were proven.

Theorem 5.1. Except for a choice of coefficients belonging to an algebraic set, there are exactly B_{S} regular solutions to a random linear-product system based on the set structure S.

The proof of the theorem consists in collecting the determinants that express the degeneracy conditions. These determinants are polynomials in the coefficients and vanish at an algebraic set.

Theorem 5.2. All isolated solutions to P (x) = 0 lie at the end of some solution path defined by a convex-linear homotopy originating at a solution of a random linear-product start system, based on a supporting set structure for P.

The idea of the proof is to embed the homotopy into an appropriate projective space and to consider the projection of the discriminant variety as an algebraic set for the continuation parameter. See [63] for an alternative proof.

A general approach to exploit product structures was developed in [80]. For systems whose polynomials are sums of products one may arrive at a much tighter bound replacing the products by one simple product. An efficient homotopy to solve the nine-point problem in mechanical design was obtained in this way.

The complexity of this homotopy based on the total degree is addressed in [10] where α-theory is applied to give bounds on the number of steps that is needed to trace the solution paths. A major result is that one can decide in polynomial time whether an average polynomial system has a solution. A similar analysis of Newton’s method in multi-projective space was recently done in [17].

While the above complexity results apply to random systems, the problem of automatically extracting and exploiting the degree structure of a polynomial system is a much harder problem. Finding an optimal multi-homogeneous grouping essentially requires the enumeration of all partitions [120]. With supporting set structures one may obtain a high success rate, see [63] for a efficient heuristic algorithm. Recent software extensions for finding optimal partitioned linear-product start systems are in [128].

5.2. Mixed Subdivisions of Newton Polytopes to compute Mixed Volumes. For (2), we collect the exponent vectors of the system P in the supports A:

\begin{matrix}P (x_{1}, x_{2}) = \\ \{ & \begin{matrix}x_{1}^{4} + x_{1}x_{2} + 1 = 0 \\ x_{1}^{3}x_{2} + x_{1}x_{2}^{2} + 1 = 0 \end{matrix} \end{matrix}  \begin{matrix}A = (A_{1}, A_{2}) \\ A_{1} = \{(0, 0), (1, 1), (4, 0)\} \\ A_{2} = \{(0, 0), (1, 2), (3, 1)\} \end{matrix}  (15)

The supports A_{1} and A_{2} span the respective Newton polytopes Q_{1} and Q_{2}.

The Cayley trick [33, Proposition 1.7, page 274] is a method to rewrite a certain resultant as a discriminant of one single polynomial with additional variables. The polyhedral version of this trick as in [102, Lemma 5.2] is due to Bernd Sturmfels. It provides a one-to-one correspondence between the cells in a mixed subdivision and a triangulation of the so-called Cayley polytope spanned by the points of A_{i} embedded in a (2n − 1)-dimensional space. See [45] for another application besides mixed-volume computation. As in [45], Figure 7 gives a “one-picture proof” of this trick, displaying the Cayley polytope for the supports A in (15). Note that this construction provides a definition for mixed subdivisions.

The Cayley polytope is spanned by the points in A_{1}, where each point of A_{1} is extended with n − 1 zero coordinates, and the points in A_{i} where each point in A_{i} is extended with the respective i-th standard basis vector, for i = 1, 2, …, n − 1. Omitting the added coordinates of this Cayley embedding, every cell in a triangulation of the Cayley polytope is identified with a cell in a mixed subdivision of the original tuple of polytopes. We can see this identification geometrically when slicing the Cayley polytope with a hyperplane that separates the embedded polytopes. As in Figure 7, the slice contains λ_{1}Q_{1} + λ_{2}Q_{2} and the cells of a mixed subdivision are cut out by the cells in a triangulation of the Cayley polytope.

Figure 7 : Cayley polytope of Q_{1} and Q_{2}.

The Cayley trick was implemented in [112] as an application of the dynamic lifting algorithm to construct regular triangulations. This method calculates the volume polynomial (4) completely. When one is only interested in the mixed volume, the method is only efficient when the supports do not differ much from each other.

To compute only the mixed volume, the lift-and-prune approach was presented in [24], using a primal model to prune in the tree of edge-edge combinations. This approach operates in two stages. First the polytopes are lifted by adding one coordinate to every point in the supports. In the second stage, one computes the facets of the lower hull of the Minkowski sum that are spanned by sums of edges. These facets constitute the mixed cells in a mixed subdivision. On the supports A in (15), we consider the lifted supports

Â = (Â_{1}, Â_{2}) \begin{matrix}Â_{1} = \{(0, 0, 1), (1, 1, 0), (4, 0, 0)\} \\ Â_{2} = \{(0, 0, 0), (1, 2, 0), (3, 1, 1)\} \end{matrix}  (16)

The lower hulls of the lifted polytopes are displayed in Figure 8.

Figure 8. Lifted polytopes Q̂_{1}, Q̂_{2}, and a regular mixed subdivision of Q̂_{1} + Q̂_{2}.

The two cells that contribute to the mixed volume are identified by inner normals α and β that satisfy systems of linear equations and inequalities:

\begin{matrix}α = (0, 0, 1) \\ \{ & \begin{aligned}4α_{1} = α_{1} + α_{2} & < & 1 \\ α_{1} + 2α_{2} = 0 & < & 3α_{1} + α_{2} + 1 \end{aligned} \end{matrix}  \begin{matrix}β = (2, −1, 1) \\ \{ & \begin{aligned}β_{1} + β_{2} = 1 & < & 4β_{1} \\ β_{1} + 2β_{2} = 0 & < & 3β_{1} + β_{2} + 1 \end{aligned} \end{matrix} (17)

These systems express that the cells correspond to facets spanned by the sum of two edges on the lower hulls of Q̂_{1} and Q̂_{2} respectively. The lift-and-prune method with a dual version of the linear inequality constrains as in (17) was elaborated in [112], exploiting the fact that several polynomials can share the same Newton polytope (see [40]) and with dimension reductions.

The geometric dual construction to Figure 8 is displayed in Figure 9.

Figure 9. On the left we see the projection of a regular mixed subdivision of Q̂_{1} + Q̂_{2}. On the right, we have the dual construction with complexes N_{∨}^{1}(Q_{1}) and N_{∨}^{1}(Q_{2}) collecting the cones of all vectors normal to the edges on the lower hulls of Q̂_{1} and Q̂_{2} respectively. The intersection of the cones contain the normals to the mixed cells.

As in [40], we assume that there are r different Newton polytopes. Given a tuple of lifted point sets Â = (Â_{1}, Â_{2}, …, Â_{r}), any lifted cell Ĉ_{v} of a regular subdivision can be characterized by its inner normal as

Ĉ_{v} = (∂_{v}Â_{1}, ∂_{v}Â_{2}, …, ∂_{v}Â_{r}).  (18)

Since conv(Ĉ_{v}) = conv(∑_{i=1}^{r} ∂_{v}Â_{i}) is a facet of the lower hull, the inner product ⟨., v⟩ attains its minimum over Â_{i} at ∂_{v}Â_{i}, i.e.,

∀â, b̂ ∈ ∂_{v}Â_{i} : ⟨â, v⟩ = ⟨b̂, v⟩, i = 1, 2, …, r,  (19)
∀â ∈ Â_{i} \ ∂_{v}Â_{i}, ∀b̂ ∈ ∂_{v}Â_{i} : ⟨â, v⟩ > ⟨b̂, v⟩, i = 1, 2, …, r.  (20)

Algorithm 5.1 (presented in [112]) gives a way to compute all mixed cells by searching for feasible solutions to the constraints (19) and (20). The algorithm generates a tree of all possible combinations of k_{i}-faces, with feasibility tests to prune branches that do not lead to mixed cells. The order of enumeration is organized so that mixed cells which share some faces also share a part of the factorization work to be done to solve the system defined by (19).

Algorithm 5.1. Pruning algorithm with shared factorizations subject to inequality constraints:

Input: (Â_{1}, Â_{2}, …, Â_{r}),  lifted point sets
k = (k_{1}, k_{2}, …, k_{r}), n = ∑_{i=1}^{r} k_{i},  A_{i} appears k_{i} times
(F̂_{1}, F̂_{2}, …, F̂_{r}).  k_{i}-faces of lower hull of conv(Â_{i})
Output: 𝔖̂_{ω} = { Ĉ ∈ Ŝ_{ω} | V_{n}(C, k) > 0 }.  collection of lifted mixed cells

At level i, 1 ≤ i < r:

DATA and INVARIANT CONDITIONS:

(M_{1}, κ): M_{1}v = 0 ⇏ v_{n+1} = 0, κ = ∑_{j=1}^{i−1} k_{j}

equalities (19)

upper triangular up to row κ

\begin{aligned}(M_{2}, κ): & M_{2}v ≥ 0 ⇏ −v_{n+1} ≥ 0 \\ & dim(M_{2}) = n − κ \end{aligned} inequalities (20)

still feasible and reduced

ALGORITHM:

for each Ĉ_{i} ∈ F̂_{i} loop  enumerate over all k_{i}-faces
Triangulate(M_{1}, κ, Ĉ_{i});  ensure invariant conditions
if M_{1}v = 0 ⇏ v_{n+1} = 0  test for feasibility w.r.t. (19)
then Eliminate(M_{1}, M_{2}, κ, Ĉ_{i}, Â_{i});  eliminate unknowns
if M_{2}v ≥ 0 ⇏ −v_{n+1} ≥ 0  test for feasibility w.r.t. (20)

then proceed to next level i + 1; end if;

end if; end for.

At level i = r:

Compute v: M_{1}v = 0;

Merge the new cell C_{v} with the list 𝔖̂_{ω}.

Note that (20) has to be weakened to ≥ type inequalities in order to be able to compute also subdivisions that are not fine. This also explains the merge operation at the end. The feasibility tests in the algorithm allow an efficient computation of the mixed cells. The conditions (19) and (20) are verified incrementally. After choosing a k_{i}-face Ĉ_{i} = {ĉ_{0i}, ĉ_{1i}, …, ĉ_{k_{i}i}} of Â_{i}, linear programming is used to check whether (Ĉ_{1}, …, Ĉ_{i}) can lead to a mixed cell in the induced subdivision.

We end this section with complexity results. The complexity of computing mixed volumes is proven [21] to be #P -hard. This complexity class is typical for all enumerative problems, since, unlike the class NP, there exists no algorithm that runs in polynomial time for arbitrary dimensions to verify that a guessed answer is correct. Although the current algorithmical practice suggests that computing mixed volumes is harder than computing volumes of polytopes (which is also known as a #P -hard problem [20]), this is not the case from a complexity point of view as shown in [21]. In [25] it is shown that the mixed volume V_{n}(Q) is bounded from below by n!vol_{n}(Q_{μ}), Q_{μ} being the polytope of minimum volume in Q.

5.3. Sparse Polynomial Systems solved by Polyhedral Homotopies. The simplest system in the polytope model that still has isolated solutions in (ℂ^{∗})^{n} has exactly two terms in every equation. Polyhedral homotopies [40] solve systems with random complex coefficients starting from sparser subsystems. For (2), the homotopy with supports Â as in (16) is

P̂(x_{1}, x_{2}, t) = \{ \begin{matrix}c_{1}x_{1}^{4}t^{0} + c_{2}x_{1}x_{2}t^{0} + c_{3}t^{1} = 0 \\ c′_{1}x_{1}^{3}x_{2}t^{1} + c′_{2}x_{1}x_{2}^{2}t^{0} + c′_{3}t^{0} = 0 \end{matrix} with c_{i}, c′_{i} ∈ ℂ^{∗}.  (21)

The exponents of t are the values of the lifting ω applied to the supports.

To find the start systems, we look at Figure 8, at the subdivision that is induced by this lifting process. The start systems have Newton polytopes spanned by one edge of the first and one edge of the second polytope. Since the two cells that contribute to the mixed volume are characterized by their inner normals α and β satisfying (17) we denote the start systems respectively by P^{α} and P^{β}. To compute start solutions, unimodular transformations make the system triangular as follows. After dividing the equations so that the constant term is present, we apply the substitution x_{1} = y_{2}, x_{2} = y_{1}^{−1}y_{2}^{3} on P^{α} as follows:

P^{α}(x) = \{ \begin{matrix}x_{1}^{3}x_{2}^{−1} + c′′_{1} = 0 \\ x_{1}x_{2}^{2} + c′′_{2} = 0 \end{matrix} P^{α}(x = y^{U}) = \{ \begin{matrix}y_{1} + c′′_{1} = 0 \\ y_{1}^{−2}y_{2}^{7} + c′′_{2} = 0 \end{matrix}  (22)

The substitution in (22) is apparent in the notation (as used in [64]) x^{V} = (y^{U})^{V} = y^{V U} = y^{L} elaborated as

\begin{aligned}( \begin{matrix}x_{1}^{3} · x_{2}^{−1} \\ x_{1}^{1} · x_{2}^{2} \end{matrix} ) & = & ( \begin{matrix}(y_{1}^{0}y_{2}^{1})^{3} · (y_{1}^{−1}y_{2}^{3})^{−1} \\ (y_{1}^{0}y_{2}^{1})^{1} · (y_{1}^{−1}y_{2}^{3})^{2} \end{matrix} ) \\ & = & ( \begin{matrix}y_{1}^{3·0−1·(−1)} · y_{2}^{3·1−1·3} \\ y_{1}^{1·0+2·(−1)} · y_{2}^{1·1+2·3} \end{matrix} ) = ( \begin{matrix}y_{1}^{1} · y_{2}^{0} \\ y_{1}^{−2} · y_{2}^{7} \end{matrix} ) . \end{aligned}  (23)

The exponents are calculated by the factorization V U = L:

[ \begin{matrix}3 & −1 \\ 1 & 2 \end{matrix} ] [ \begin{matrix}0 & 1 \\ −1 & 3 \end{matrix} ] = [ \begin{matrix}1 & 0 \\ −2 & 7 \end{matrix} ] .  (24)

Since det(U) = 1, the matrix U is called unimodular.

The polyhedral homotopy (21) directly extends the solutions of P^{α} to the target system. To obtain a homotopy starting at the solutions of P^{β}, we substitute in (21) x_{1} ← x̃_{1}t^{β_{1}}, x_{2} ← x̃_{2}t^{β_{2}} and clear out the lowest powers of t. This construction appeared in [40] and provides an algorithmic proof of the following theorem.

Theorem 5.3. (Bernshteǐn [7, Theorem A]) For a general choice of coefficients for P, the system P (x) = 0 has exactly as many regular solutions as its mixed volume V_{n}(Q).

The original algorithm Bernshteǐn used in his proof was implemented in [110].

For the numerical stability of polyhedral continuation, it is important to have subdivisions induced by low lifting values, since those influence the power of the continuation parameter. In [112] explicit lower bounds on integer lifting values were derived, but unfortunately the dynamic lifting algorithm does not generalize that well [65] if one is only interested in the mixed cells of a mixed subdivision. A balancing method was proposed in [30] to improve the stability of homotopies induced by random floating-point lifting values.

Once all solutions to a polynomial system with randomly generated coefficients are computed, we use cheater’s homotopy to solve any specific system with the same Newton polytopes. One could say that polyhedral homotopies have removed the cheating part. The main advantage of polyhedral methods is that the mixed volume is a much sharper root count in most applications, leading to fewer paths to trace. They also allow more flexibility to exploit symmetry as demonstrated in [108].

In case the system has fewer isolated solutions than the mixed volume, we consider the face systems. Define the face of a polynomial p with support A as follows:

p(x) = ∑_{a∈A} c_{a}x^{a} has faces ∂_{v}p(x) = ∑_{a∈∂_{v}A} c_{a}x^{a} for v ≠ 0.  (25)

For v ≠ 0, the corresponding face system of P is ∂_{v}P = (∂_{v}p_{1}, ∂_{v}p_{2}, …, ∂_{v}p_{n}).

Theorem 5.4. (Bernshteǐn [7, Theorem B]) Suppose V_{n}(Q) > 0. Then, P (x) = 0 has fewer than V_{n}(Q) isolated solutions if and only if ∂_{v}P (x) = 0 has a solution in (ℂ^{∗})^{n}, for v ≠ 0.

As is the case for our running example (2), the Newton polytopes may be in generic position such that for any nonzero choice of the coefficients, the system has exactly as many isolated solutions as the mixed volume. In practical applications however, how can we decide whether paths are really going towards infinity? Relying on the actual computed values is arbitrarily, because 10^{4} is as far from infinity as 10^{8}, so we need algebraic structural data to certify the divergence.

In the polyhedral end game [43] solution paths are represented by power series expansions:

\{ \begin{matrix}x_{i}(s) & = & a_{i}s^{v_{i}}(1 + O(s)) \\ t(s) & = & 1 − s^{m} \end{matrix}  t ≈ 1, s ≈ 0.  (26)

The winding number m is lower than or equal to the multiplicity of the solution. For a solution diverging to infinity or to a zero-component solution we observe that v_{i} ≠ 0. According to Theorem 5.4, this solution vanishes at a face system ∂_{v}P (same v with components v_{i} as in (26)), certifying the divergence.

To check whether a solution path really diverges is equivalent to the test on the value for v_{i}. A first-order approximation of v_{i} can be computed by

\frac{log |x_{i}(s_{1})| − log |x_{i}(s_{0})|}{log(s_{1}) − log(s_{0})} = v_{i} + O(s_{0}),  (27)

with 0 < s_{1} < s_{0}. The above formula assumes the correct value for m. To compute m, solution paths are sampled geometrically with ratio h as s_{k} = h^{k/m}s_{0}. The errors on the estimates for v_{i} are

e_{i}^{(k)} = (log |x_{i}(s_{k})| − log |x_{i}(s_{k+1})|) − (log |x_{i}(s_{k+1})| − log |x_{i}(s_{k+2})|)  (28)
= c_{1}h^{k/m}s_{0}(1 + O(h^{k/m})).  (29)

An estimate for m is derived from two consecutive errors e_{i}^{(k)}. Extrapolation improves this estimate. So, by an inexpensive side calculation at the end of the paths, we obtain important structural algebraic information about the system.

Recall that V_{n}(Q) count the roots in (ℂ^{∗})^{n}. Using Newton polytopes to count affine roots (i.e.: in ℂ^{n} instead of (ℂ^{∗})^{n}) was proposed in [85] with the notion of shadowed polytopes obtained by the substitution x_{i} ← x_{i} + c_{i} for arbitrary constants c_{i}. To arrive at sharper bounds, it suffices (see [62] and [86]) to add a random constant to every equation. Stable mixed volumes [42] provide a generically sharp affine root count. The constructions in [26] and [29] avoid the use of recursive liftings to compute stable mixed volumes. Further developments and generalizations can be found in [87].

5.4. Determinantal Polynomials arising in Enumerative Geometry. Homotopies for solving problems in enumerative geometry appeared in [44]. The algorithms in the numerical Schubert calculus originated from questions in real enumerative geometry [96, 97] and have their main application to the pole placement problem [12], [83, 84], [88], [90] in control theory.

The enumerative problems are formalized in some “finiteness” theorems, avoiding the explicit but involved (as in (6)) formulas for the root counts.

Theorem 5.5. The number of r-planes meeting mr general m-planes in ℂ^{m+r} is a finite constant.

The first homotopy presented in [44] uses a Gröbner basis for the ideal that defines G_{mr}, as is derived in [101]. By Gröbner bases questions concerning any polynomial system are solved by relation to monomial equations. Every Gröbner basis defines a flat deformation, which preserves the structure of the solution set [22]. Geometrically, this type of deformation is used to collapse the solution set in projective space to the coordinate hyperplanes, or in the opposite direction, to extend the solutions of the monomial equations to those of the original system. The flat deformations that are obtained in this way are similar to toric deformations in the sense that one moves from the solutions of a subsystem to the solutions of the whole system.

The Gröbner homotopies of [44] work in the synthetic Plücker embedding, and need to take the large set of defining equations of G_{mr} into account. When expanding the minors into local coordinates, these equations are automatically satisfied, which leads to a much smaller polynomial system. Consequently, the second type of homotopies of [44], the socalled SAGBI homotopies are more efficient. Instead of an ideal, we now have a subalgebra and work with a SAGBI basis, i.e.: the Subalgebra Analogue to a Gröbner Basis for an Ideal. The term order selects the monomials on the diagonal as the dominant ones. This implies that in the flat deformation (see [103] for a general description) only the diagonal monomials remain at t = 0.

For m = 2 = r, the equations of the SAGBI homotopy in determinantal form are

p_{i}(x) = det [ \begin{matrix}c_{11}^{(i)} & c_{12}^{(i)} \\ c_{21}^{(i)} & c_{22}^{(i)} \\ c_{31}^{(i)} & c_{32}^{(i)} \\ c_{41}^{(i)} & c_{42}^{(i)} \end{matrix} | \begin{matrix}x_{11} & x_{12} \\ x_{21}t & x_{22} \\ 1 & 0 \\ 0 & 1 \end{matrix} ] = 0, i = 1, 2, 3, 4,  (30)

where the coefficients c_{kl}^{(i)} are random complex constants. In expanding the minors of (30), the lowest power of t is divided out, minor per minor. The system at t = 0 is solved by polyhedral continuation. The system at t = 1 serves as start system in the cheater’s homotopy to solve any system with particular real values for the coefficients c_{kl}^{(i)}. Figure 10 outlines the structure of the general solver.

Binomial Systems polyhedral homotopy

Generic Complex System flat deform.

Generic Complex Problem cheater’s homotopy

Specific Real Problem

Figure 10. The SAGBI homotopy is at the center of the concatenation.

SAGBI homotopies have been implemented [114] to verify some large instances of input planes for which it was conjectured that all solution planes would be real. We refer to [89] and [98] for related work on these conjectures. In [99] an asymptotic choice of inputs is generated for which all solutions are proven to be real.

The third type of homotopies presented in [44] are the so-called Pieri homotopies. Since they are closer to the intrinsic geometric nature, they are applicable to a broader class of problems. In particular, we obtain an effective proof for the following.

Theorem 5.6. The number of r-planes meeting n general (m+1−k_{i})-planes in ℂ^{m+r}, with k_{1} + k_{2} + ⋯ + k_{n} = mr, is a finite constant.

Note that when all k_{i} = 1, we arrive at Theorem 5.5. For general k_{i} ≠ 1, we are not aware of any explicit formulas for the number of roots.

Figure 11 shows a part of a cellular decomposition of G_{22} with the determinantal equations. We specialize the pattern X that represents a solution line by setting some of its coordinates to zero. This specialization determines a specialization of the input lines as follows: take those basis vectors not indexed by rows of X where zeroes have been introduced. The special line S_{X} for this example is as in (31) spanned by the first and third standard basis vector.

Figure 11 pictures patterns of the moving 2-planes in the Pieri homotopy algorithm for the case (m, r) = (2, 2), see Figure 6. The bottom matrix is the general representation of a solution that intersects already the two input lines spanned by the standard basis vectors. At the leaves of the tree by linear algebra operations we can intersect with a third input line.

[figure: det [S_{X} | \begin{matrix}x_{11} & 0 \\ 0 & 0 \\ 0 & x_{32} \\ 0 & x_{42} \end{matrix} ] = 0; det [S_{X} | \begin{matrix}x_{11} & 0 \\ x_{21} & 0 \\ 0 & x_{32} \\ 0 & 0 \end{matrix} ] = 0; det[L_{3}|X] = det [L_{3} | \begin{matrix}x_{11} & 0 \\ x_{21} & 0 \\ 0 & x_{32} \\ 0 & x_{42} \end{matrix} ] = 0; [2 4]; [1 4]; [2 3]]

Figure 11. Part of a cellular decomposition of the Grassmannian of all 2planes. At the right we have the short-hand notation with brackets. The bracket [2 4] contains the row indices to the lowest nonzero entries in X.

Moving down the poset, we deform form the left configuration in Figure 6 to the general problem.

Denote by L_{1} and L_{2} the lines already met by X. At the leaves of the tree in Figure 11 we intersect with the fourth line L_{4}. The special position for the third line L_{3} is represented by the matrix S_{X}, which intersects any X with coordinates as at the leaves of the tree. In the homotopy H(X, t) = 0 we deform the line spanned by the columns of S_{X} to line L_{3}, for t = 0 to t = 1.

S_{X} = [ \begin{matrix}1 & 0 \\ 0 & 0 \\ 0 & 1 \\ 0 & 0 \end{matrix} ]  H(X, t) = \{ \begin{matrix}det(L_{4}|X) = 0 \\ det((1 − t)S_{X} + tL_{3}|X) = 0 \end{matrix} t ∈ [0, 1].  (31)

Every solution X(t) of H(X(t), t) = 0 intersects already three lines: L_{1}, L_{2} and L_{3}. At the end, for t = 1, X also meets the line L_{3} in a point.

The homotopy (31) deforms two solution lines, starting at patterns which have their row indices for the lowest nonzero entries respectively as in [1 4] and in [2 3]. The correctness of this homotopy (proven in [44] and [46]) justifies the formal root count using the localization poset. This combinatorial root count proceeds in two stages. First we build up the poset, from the bottom up, diminishing the entries in the brackets under the restriction that the same entry never occurs twice or more. Secondly, we descend from the top of the poset, collecting and adding up the counts at the nodes in the poset. More examples and variations are in [46].

To solve the general intersection problem of Theorem 5.6, the special (m + 1 − k_{i})-planes lie in the intersection of special m-planes. In the construction of the poset one has to follow additional rules as to ensure a solution that meets the intersection of special m-planes. We refer to [46] for details.

The third enumerative problem we can solve is formalized as follows.

Theorem 5.7. The number of all maps of degree q that produce r-planes in ℂ^{m+r} meeting mr + q(m + r) general m-planes at specified interpolation points is a finite constant.

In [83, 84] and [122] explicit formulas are given for this finite constant along with other combinatorial identities. Following a hint of Frank Sottile (see also [100]) and reverse engineering on the root counts in [83], Pieri homotopies were developed in [46] whose correctness yields a proof for Theorem 5.7.

The analogue to Figure 11 for maps of degree one into G_{22} is displayed in Figure 12.

[figure: det [S_{X} | \begin{matrix}x_{11}^{0} & x_{12}^{1}s \\ x_{21}^{0} & x_{22}^{0}t \\ 0 & x_{32}^{0}t \\ 0 & x_{42}^{0}t \end{matrix} ] = 0; det [S_{X} | \begin{matrix}x_{11}^{0} & 0 \\ x_{21}^{0} & x_{22}^{0} \\ x_{31}^{0} & x_{32}^{0} \\ 0 & x_{42}^{0} \end{matrix} ] = 0; det[L_{n}|X(s, t)] = det [L_{n} | \begin{matrix}x_{11}^{0} & x_{12}^{1}s \\ x_{21}^{0} & x_{22}^{0}t \\ x_{31}^{0} & x_{32}^{0}t \\ 0 & x_{42}^{0}t \end{matrix} ] = 0; [3 5]; [2 5]; [3 4]]

Figure 12. Part of a cellular decomposition of the Grassmannian of maps of degree 1 that produce 2-planes in projective 3-space. The bracket notation at the right corresponds to a matrix representation of the coefficients of the map X(s, t).

To solve the problem in Theorem 5.7 we need a special position for the interpolation points. By moving those to infinity, the dominant monomials in the maps allow to re-use the same special m-planes, whose entries should be considered modulo m + r. The homotopy to satisfy the n-th intersection condition is:

H(X(s, t), s, t) = \{ \begin{matrix}det(L_{i}|X(s_{i}, t_{i})) = 0 & i = 1, 2, …, n − 1 \\ det((1 − t)S_{X} + tL_{n}|X(s, t)) = 0 \\ (s − 1)(1 − t) + (s − s_{n})t = 0 & t ∈ [0, 1] \end{matrix}  (32)

Note that the continuation parameter t moves the interpolation point from infinity, at (s, t) = (1, 0), to the specific value (s, t) = (s_{n}, 1).

See [100] for information on the selection of the input planes so that all maps are real.

As an example of another problem in enumerative geometry we mention the 27 lines on a cubic surface in 3-space. According to [81], this is one of the gems hidden in the rag-bag of projective geometry. In [92], the 27 lines are determined by breaking up the cubic surface into three planes in a continuous way such that each intermediate position is nonsingular. It is shown that this continuous variation is also valid in the real field.

6. Numerical Software for Solving Polynomial Systems

In computer algebra one wants to compute exactly as long as possible and to defer the approximate calculations to the very last end. Exactly the opposite way is taken in homotopy methods: here we use floating-point arithmetic and only increase the precision when needed.

Next we mention programs with special features for polynomial systems. See [5, Chapter VIII] for a list of available software for path following. HOMPACK [74, 124] is a general continuation package with a polynomial driver. It has been parallelized [6, 36], extended with an end game [95], and upgraded [126] to Fortran 90. POLSYS PLP [128] provides linear-product root counts to be used in conjunction with HOMPACK90. The Fortran code for CONSOL is contained in [71, Appendix 6]. The C-program pss [66] applies homotopy continuation with verification by α-theory. Pelican [39, 41] implements in C the polyhedral methods of [40]. Efficient Fortran software for polyhedral continuation with facilities to compute all affine roots is used in [29]. The computation of mixed volumes with the C program mvlp [23, 24] is a crucial step for sparse resultants [117]. A distributed version has been created in [34].

PHC is written in Ada and originated during the doctoral research of the author [113]. Executable versions were first released at the PoSSo open workshop on software [111]. The public release of the sources is described in [115]. The package is organized as a tool box, organized along four stages of the solver. Figure 13 presents the flow of the solver. The package is menu-driven and file-oriented. A general-purpose black-box solver is available.

1. Preconditioning

⋄ Coefficient Scaling ⋄ Reduction of degrees 4. Validation

⋄ Refining of the roots ⋄ Analysis of condition #s

2. Root Counting

⋄ Bézout : degrees

⋄ Bernshtein : polytopes ⋄ Schubert : SAGBI/Pieri 3. Homotopy Continuation

⋄ Fix continuation parameters ⋄ Choose Predictor-Corrector

[figure]

Figure 13. The four stages in the flow of the PHC solver.

The new second release of PHC uses Ada 95 concepts in the construction of the mathematical library. It is developed with the freely available gnu-ada compiler (currently at version 3.11p) on various platforms. To run the software no compilation is needed, as binaries are available for Unix Workstations: running SUN Solaris and SGI Irix, and Pentium PCs: running Linux and Solaris. The portability of PHC is ensured by the gnu-ada compiler.

Another main feature of the second release are the homotopy methods for the Schubert calculus. Implementing those homotopies was a matter of plugging in the equations and calling the path trackers. The third release of the package should offer a more comprehensive environment to construct homotopies, providing an easier access to the two main computational engines: mixed-volume computation and polynomial continuation.

7. The Database of Applications

The polynomial systems in scientific and engineering models are a continuing source of open problems. Systems that come from academic questions are often conjectures providing computational evidence in a developing theory. In various engineering disciplines polynomial systems represent a modeling problem, e.g.: a mechanical device. The origin of a polynomial system matters when the original problem formulation does not admit well-conditioned solutions. As a general method to deal with badly scaled systems to compute equilibria of chemical reaction systems, coefficient and equation scaling was developed in [69], see also [71, Chapter 5] and [124].

The collection of test systems is organized as a database and available via the author’s web pages, A good test example reveals properties of the solution method and has a meaningful application. Besides the algebraic formulation it contains the fields: title (meaningful description), references (problem source), root counts (Bézout bounds and mixed volume), and solution list.

Instead of producing a huge list with an overview, we pick some important case studies.

katsura-n: (magnetism problem [47]) The number of solutions equals the total degree D = 2^{n}, so the homotopy based on D is optimal to solve this problem. Because the constant term is missing in all except one equation, the system is an interesting test problem for affine polyhedral methods.

camera1s: (computer vision [28]) The system models the displacement of a camera between two positions in a static environment [23]. The multi-homogeneous homotopy is optimal for this problem, requiring 20 solution paths to trace instead of D = 64.

gamentwo: (totally mixed Nash equilibria for n players with two strategies [67, 68]) This is another instance where multi-homogeneous homotopies are optimal. The number of solutions grows like n!e^{−1} as n → ∞. The largest system that is currently completely solvable is for n = 8 requiring 14,833 paths to trace. Situations exist for which all solutions are meaningful.

cassou: (real algebraic geometry) This system illustrates the success story of polyhedral homotopies: the total degree equals 1,344, best known Bézout bound is 312 (see [63]), whereas the mixed volume gives 24. Still eight paths are diverging to infinity and polyhedral end games [43] are needed to separate those diverging paths from the from the other finite ill-conditioned roots.

cyclic-n: (Fourier transforms [8, 9]) For n = 7, polyhedral homotopies are optimal, with all 924 paths leading to finite solutions. For n ≥ 8, the mixed volume overestimates the number of roots and there are components of solutions. In [94] the degrees of the components were computed for n = 8, 9. There are 34940 cyclic 10-roots, generated by 1747 solutions.

pole28sys: (pole placement problem [12]) This system illustrates the efficiency of SAGBI homotopies for verifying a conjecture in real algebraic geometry [98]. With the input planes chosen to osculate a rational normal curve, an instance with all 1,430 solutions real and isolated was solved in [114]. The problem is relevant to control theory [90].

stewgou40: (mechanism design [18]) Whether the Stewart-Gough parallel platform in robotics could have all its 40 solutions real was a notorious open problem, until recently, as it was solved by numerical continuation methods [18]. The problem formulation in [18] is highly deficient: the mixed volume equals 1,536 whereas only 40 solution paths will converge.

We emphasize that we have optimal homotopies for three classes of polynomial systems, but not for all possible structures. Although one can solve a modelling problem by a blackbox polynomial-system solver, knowing the origin of the problem leads in most cases to more favorable algebraic formulations that help the resolution of a polynomial system. To produce really meaningful solutions one often has to be close to the source of the problems and be able to interact with the people who formulate the polynomial systems.

In closing this section we list some notable usages of PHC. Charles Wampler [121] used a preliminary version of PHC to count the roots of various systems in mechanical design. Root counts for linear subspace intersections in the Schubert calculus were computed by Frank Sottile, see [98] for various tables. A third example comes from computer graphics. To show that the 12 lines tangent to four given spheres can all be real, Thorsten Theobald used PHC, choosing appropriate parameters in the algebraic formulation set up by Cassiano Durand.

8. Closing Remarks and Open Problems

The three classes presented in this paper are by no means exhaustive, but give an idea of what can be done with homotopies to solve polynomial systems. The root counts constitute the theoretical backbone for general-purpose black-box solving. Yet, the homotopy methods are flexible enough to exploit a particular geometrical situation, with guaranteed optimal complexity when applied to generic instances.

From algebraic geometry formal procedures based on intersection theory count the number of solutions to classes of polynomial systems. Examples are the theorems of Bézout, Bernshteǐn and Schubert. For these situations we construct a start system and have a homotopy to deform the solutions to this start system to the solutions to any specific problem. There are many other cases for which one knows how to count but not how to deform and solve efficiently. Research in homotopy methods is aimed at turning the formal root counts into effective numerical methods. As open problem we can ask for a meta-homotopy method to connect formal root counting methods to solving generic systems and deformation procedures.

In most applications, only the real solutions are important. Once we know an optimal homotopy to solve the problem in the complex case, we would like to know whether all solutions can be real and how the real solutions are distributed. The reality question appears for instance in the theory of totally mixed Nash equilibria and in the pole placement problem. Finding well-conditioned instances of fully real problems can be done by homotopy methods. The finding of 40 real solutions to the Stewart-Gough platform [18] is perhaps the most striking example. The question is to find an efficient procedure to deform from the complex case to the fully real case.

Acknowledgments. The interest in homotopy methods by the FRISCO project has stimulated the author’s research. The author is deeply indebted to all his co-authors. The interactions with Birk Huber and Frank Sottile at MSRI were influential in this current treatment of homotopies.

References

[1] S.S. Abhyankar. Algebraic Geometry for Scientists and Engineers, volume 35 of Mathematical Surveys and Monographs. American Mathematical Society, Providence, Rhode Island, 1990.

[2] E.L. Allgower. Bifurcation arising in the calculation of critical points via homotopy methods. In T. Kupper, H.D. Mittleman, and N. Weber, editors, Numerical Methods for Bifurcation Problems, pages 15–28. Birkhäuser, Boston, 1984.

[3] E.L. Allgower and K. Georg. Numerical Continuation Methods, an Introduction, volume 13 of Springer Ser. in Comput. Math. Springer–Verlag, Berlin Heidelberg New York, 1990.

[4] E.L. Allgower and K. Georg. Continuation and path following. Acta Numerica, pages 1–64, 1993.

[5] E.L. Allgower and K. Georg. Numerical Path Following. In P.G. Ciarlet and J.L. Lions, editors, Techniques of Scientific Computing (Part 2), volume 5 of Handbook of Numerical Analysis, pages 3–203. North-Holland, 1997.

[6] D.C.S. Allison, A. Chakraborty, and L.T. Watson. Granularity issues for solving polynomial systems via globally convergent algorithms on a hypercube. J. of Supercomputing, 3:5–20, 1989.

[7] D.N. Bernshteǐn. The number of roots of a system of equations. Functional Anal. Appl., 9(3):183–185, 1975. Translated from Funktsional. Anal. i Prilozhen., 9(3):1–4,1975.

[8] G. Björck. Functions of modulus one on Z_{n} whose Fourier transforms have constant modulus, and ”cyclic n-roots”. In J.S. Byrnes and J.F. Byrnes, editors, Recent Advances in Fourier Analysis and its Applications, volume 315 of NATO Adv. Sci. Inst. Ser. C: Math. Phys. Sci., pages 131–140. Kluwer, 1989.

[9] G. Björck and R. Fröberg. A faster way to count the solutions of inhomogeneous systems of algebraic equations, with applications to cyclic n-roots. J. Symbolic Computation, 12(3):329–336, 1991.

[10] L. Blum, F. Cucker, M. Shub, and S. Smale. Complexity and Real Computation. Springer–Verlag, New York, 1997.

[11] P. Brunovský and P. Meravý. Solving systems of polynomial equations by bounded and real homotopy. Numer. Math., 43(3):397–418, 1984.

[12] C.I. Byrnes. Pole assignment by output feedback. In H. Nijmacher and J.M. Schumacher, editors, Three Decades of Mathematical Systems Theory, volume 135 of Lecture Notes in Control and Inform. Sci., pages 13–78. Springer–Verlag, Berlin, 1989.

[13] J. Canny and J.M. Rojas. An optimal condition for determining the exact number of roots of a polynomial system. In Proceedings of the 1991 International Symposium on Symbolic and Algebraic Computation, pages 96–101. ACM, 1991.

[14] S.N. Chow, J. Mallet-Paret, and J.A. Yorke. Homotopy method for locating all zeros of a system of polynomials. In H.O. Peitgen and H.O. Walther, editors, Functional differential equations and approximation of fixed points, volume 730 of Lecture Notes in Mathematics, pages 77–88, Berlin Heidelberg, 1979. Springer–Verlag.

[15] D. Cox, J. Little, and D. O’Shea. Ideals, Varieties and Algorithms. An Introduction to Computational Algebraic Geometry and Commutative Algebra. Undergraduate Texts in Mathematics. Springer–Verlag, New York, 2nd edition, 1997.

[16] D. Cox, J. Little, and D. O’Shea. Using Algebraic Geometry, volume 185 of Graduate Texts in Mathematics. Springer–Verlag, New York, 1998.

[17] J.P. Dedieu and M. Shub. Multihomogeneous Newton methods. To appear in Math. Comp.

[18] P. Dietmaier. The Stewart-Gough platform of general geometry can have 40 real postures. In J. Lenarcic and M.L. Husty, editors, Advances in Robot Kinematics: Analysis and Control, pages 1–10. Kluwer Academic Publishers, Dordrecht, 1998.

[19] F.J. Drexler. Eine Methode zur Berechnung sämtlicher Lösungen von Polynomgleichungssystemen. Numer. Math., 29(1):45–58, 1977.

[20] M.E. Dyer and A.M. Frieze. On the complexity of computing the volume of a polyhedron. SIAM J. Comput., 17(5):967–974, 1988.

[21] M. Dyer, P. Gritzmann, and A. Hufnagel. On the complexity of computing mixed volumes. SIAM J. Comput., 27(2):356–400, 1998.

[22] D. Eisenbud. Commutative Algebra with a View Toward Algebraic Geometry, volume 150 of Graduate Texts in Mathematics. Springer–Verlag, New York, 1995.

[23] I.Z. Emiris. Sparse Elimination and Applications in Kinematics. PhD thesis, Computer Science Division, Dept. of Electrical Engineering and Computer Science, University of California, Berkeley, 1994.

[24] I.Z. Emiris and J.F. Canny. Efficient incremental algorithms for the sparse resultant and the mixed volume. J. Symbolic Computation, 20(2):117–149, 1995. Software available at

http://www.inria.fr/saga/emiris.

[25] I.Z. Emiris. On the complexity of sparse elimination. J. Complexity, 12(2):134–166, 1996.

[26] I.Z. Emiris and J. Verschelde. How to count efficiently all affine roots of a polynomial system. Discrete Applied Mathematics, 93(1):21–32, 1999.

[27] P. Falb. Methods of Algebraic Geometry in Control Theory: Part I, Scalar Linear Systems and Affine Algebraic Geometry, volume 4 of Systems & Control: Foundations & Applications. Birkhäuser, Boston, 1990.

[28] O.D. Faugeras and S. Maybank. Motion from point matches: multiplicity of solutions. Intern. J. Comp. Vision, 4:225–246, 1990.

[29] T. Gao, T.Y. Li, and X. Wang. Finding isolated zeros of polynomial systems in C^{n} with stable mixed volumes. To appear in J. of Symbolic Computation.

[30] T. Gao, T.Y. Li, J. Verschelde, and M. Wu. Balancing the lifting values to improve the numerical stability of polyhedral homotopy continuation methods. To appear in Applied Math. Comput.

[31] C.B. Garcia and W.I. Zangwill. Finding all solutions to polynomial systems and other systems of equations. Math. Programming, 16(2):159–176, 1979.

[32] C.B. Garcia and T.Y. Li. On the number of solutions to polynomial systems of equations. SIAM J. Numer. Anal., 17(4):540–546, 1980.

[33] I.M. Gel’fand, M.M. Kapranov, and A.V. Zelevinsky. Discriminants, Resultants and Multidimensional Determinants. Birkhäuser, Boston, 1994.

[34] T. Giordano. Implémention distribuée du calcul du volume mixte. Master’s thesis, University of Nice, Sophia-Antipolis, 1996.

[35] J.E. Goodman and J. O’Rourke (editors). Handbook of Discrete and Computational Geometry. CRC Press, New York, 1997.

[36] S. Harimoto and L.T. Watson. The granularity of homotopy algorithms for polynomial systems of equations. In G. Rodrigue, editor, Parallel processing for scientific computing, pages 115–120. SIAM, 1989.

[37] J. Harris. Algebraic Geometry, A first Course, volume 133 of Graduate Texts in Mathematics. Springer– Verlag, New York, 1992. Corrected third printing, 1995.

[38] M.E. Henderson and H.B. Keller. Complex bifurcation from real paths. SIAM J. Appl. Math., 50:460– 482, 1990.

[39] B. Huber. Pelican manual. Available at http://www.msri.org/people/members/birk.

[40] B. Huber and B. Sturmfels. A polyhedral method for solving sparse polynomial systems. Math. Comp., 64(212):1541–1555, 1995.

[41] B.T. Huber. Solving Sparse Polynomial Systems. PhD thesis, Cornell University, 1996. Available at http://www.msri.org/people/members/birk.

[42] B. Huber and B. Sturmfels. Bernstein’s theorem in affine space. Discrete Comput. Geom., 17(2):137– 141, 1997.

[43] B. Huber and J. Verschelde. Polyhedral end games for polynomial continuation. Numerical Algorithms, 18(1):91–108, 1998.

[44] B. Huber, F. Sottile, and B. Sturmfels. Numerical Schubert calculus. J. of Symbolic Computation, 26(6):767–788, 1998.

[45] B Huber, J. Rambau, and F. Santos. The Cayley trick, lifting subdivisions and the Bohne-Dress theorem on zonotopal tilings. Submitted to the Journal of the European Mathematical Society., 1999.

[46] B. Huber and J. Verschelde. Pieri homotopies for problems in enumerative geometry applied to pole placement in linear systems control. MSRI Preprint #1998-064, 1999.

[47] S. Katsura. Users posing problems to posso. PoSSo Newsletter, no. 2, July 1994, edited by L. Gonzalez-Vega and T. Recio. Available at http://posso.dm.unipi.it/.

[48] G.R. Kempf. Algebraic Varieties, volume 172 of London Mathematical Society Lecture Note Series. Cambridge University Press, 1993.

[49] A.G. Khovanskiǐ. Newton polyhedra and the genus of complete intersections. Functional Anal. Appl., 12(1):38–46, 1978. Translated from Funktsional. Anal. i Prilozhen., 12(1),51–61,1978.

[50] H. Kobayashi, H. Suzuki, and Y. Sakai. Numerical calculation of the multiplicity of a solution to algebraic equations. Math. Comp., 67(221):257–270, 1998.

[51] A.G. Kushnirenko. Newton Polytopes and the Bézout Theorem. Functional Anal. Appl., 10(3):233–235, 1976. Translated from Funktsional. Anal. i Prilozhen., 10(3),82–83,1976.

[52] T.Y. Li. On Chow, Mallet-Paret and Yorke homotopy for solving systems of polynomials. Bulletin of the Institute of Mathematics. Acad. Sin., 11:433–437, 1983.

[53] T.Y. Li. Solving polynomial systems. The Mathematical Intelligencer, 9(3):33–39, 1987.

[54] T.Y. Li and T. Sauer. Regularity results for solving systems of polynomials by homotopy method. Numer. Math., 50(3):283–289, 1987.

[55] T.Y. Li, T. Sauer, and J.A. Yorke. Numerical solution of a class of deficient polynomial systems. SIAM J. Numer. Anal., 24(2):435–451, 1987.

[56] T.Y. Li, T. Sauer, and J.A. Yorke. The random product homotopy and deficient polynomial systems. Numer. Math., 51(5):481–500, 1987.

[57] T.Y. Li, T. Sauer, and J.A. Yorke. The cheater’s homotopy: an efficient procedure for solving systems of polynomial equations. SIAM J. Numer. Anal., 26(5):1241–1251, 1989.

[58] T.Y. Li and X. Wang. Solving deficient polynomial systems with homotopies which keep the subschemes at infinity invariant. Math. Comp., 56(194):693–710, 1991.

[59] T.Y. Li and X. Wang. Nonlinear homotopies for solving deficient polynomial systems with parameters. SIAM J. Numer. Anal., 29(4):1104–1118, 1992.

[60] T.Y. Li and X. Wang. Solving real polynomial systems with real homotopies. Math. Comp., 60:669–680, 1993.

[61] T.Y. Li and X. Wang. Higher order turning points. Appl. Math. Comput., 64:155–166, 1994.

[62] T.Y. Li and X. Wang. The BKK root count in C^{n}. Math. Comp., 65(216):1477–1484, 1996.

[63] T.Y. Li, T. Wang, and X. Wang. Random product homotopy with minimal BKK bound. In J. Renegar, M. Shub, and S. Smale, editors, The Mathematics of Numerical Analysis, volume 32 of Lectures in Applied Mathematics, pages 503–512, 1996. Proceedings of the AMS-SIAM Summer Seminar in Applied Mathematics, Park City, Utah, July 17-August 11, 1995, Park City, Utah.

[64] T.Y. Li. Numerical solution of multivariate polynomial systems by homotopy continuation methods. Acta Numerica, 6:399–436, 1997.

[65] T.Y. Li and X. Wang. Counterexamples to the connectivity conjecture of the mixed cells. Discrete Comput. Geom., 20(4):515–521, 1998.

[66] G. Malajovich. pss 2 – Polynomial System Solver. Available at http://www.labma.ufrj.br:80/˜ gregorio.

[67] R.D. McKelvey and A. McLennan. The maximal number of regular totally mixed Nash equilibria. Journal of Economic Theory, 72:411–425, 1997.

[68] A. McLennan. The maximal generic number of pure Nash equilibria. Journal of Economic Theory, 72:408–410, 1997.

[69] K. Meintjes and A.P. Morgan. A methodology for solving chemical equilibrium systems. Appl. Math. Comput., 22(4):333–361, 1987.

[70] A.P. Morgan. A method for computing all solutions to systems of polynomial equations. ACM Trans. Math. Softw., 9(1):1–17, 1983.

[71] A. Morgan. Solving polynomial systems using continuation for engineering and scientific problems. Prentice-Hall, Englewood Cliffs, N.J., 1987.

[72] A. Morgan and A. Sommese. A homotopy for solving general polynomial systems that respects mhomogeneous structures. Appl. Math. Comput., 24(2):101–113, 1987.

[73] A. Morgan and A. Sommese. Computing all solutions to polynomial systems using homotopy continuation. Appl. Math. Comput., 24(2):115–138, 1987. Errata: Appl. Math. Comput. 51 (1992), p. 209.

[74] A.P. Morgan, A.J. Sommese, and L.T. Watson. Finding all isolated solutions to polynomial systems using HOMPACK. ACM Trans. Math. Softw., 15(2):93–122, 1989.

[75] A.P. Morgan and A.J. Sommese. Coefficient-parameter polynomial continuation. Appl. Math. Comput., 29(2):123–160, 1989. Errata: Appl. Math. Comput. 51:207(1992).

[76] A.P. Morgan and A.J. Sommese. Generically nonsingular polynomial continuation. In E.L. Allgower and K. Georg, editors, Computational Solution of Nonlinear Systems of Equations, pages 467–493, Providence, R.I., 1990. AMS.

[77] A.P. Morgan, A.J. Sommese, and C.W. Wampler. Computing singular solutions to nonlinear analytic systems. Numer. Math., 58(7):669–684, 1991.

[78] A.P. Morgan, A.J. Sommese, and C.W. Wampler. Computing singular solutions to polynomial systems. Adv. Appl. Math., 13(3):305–327, 1992.

[79] A.P. Morgan, A.J. Sommese, and C.W. Wampler. A power series method for computing singular solutions to nonlinear analytic systems. Numer. Math., 63:391–409, 1992.

[80] A.P. Morgan, A.J. Sommese, and C.W. Wampler. A product-decomposition theorem for bounding Bézout numbers. SIAM J. Numer. Anal., 32(4):1308–1325, 1995.

[81] D. Mumford. Algebraic Geometry I; Complex Projective Varieties, volume 221 of Grundlehren der mathematischen Wissenschaften. Springer–Verlag, Berlin Heidelberg New York, 1976.

[82] M. Pieri. Formule di coincidenza per le serie algebriche ∞^{n} di coppie di punti dello spazio a n dimensioni. Rend. Circ. Mat. Palermo, 5:252–268, 1891.

[83] M.S. Ravi, J. Rosenthal, and X. Wang. Dynamic pole placement assignment and Schubert calculus. SIAM J. Control and Optimization, 34(3):813–832, 1996.

[84] M.S. Ravi, J. Rosenthal, and X. Wang. Degree of the generalized Plücker embedding of a quot scheme and quatum cohomology. Math. Ann., 311:11–26, 1998.

[85] J.M. Rojas. A convex geometric approach to counting the roots of a polynomial system. Theoret. Comput. Sci., 133(1):105–140, 1994. Extensions and corrections available at http://www.cityu.edu.hk/ma/staff/rojas/.

[86] J.M. Rojas and X. Wang. Counting affine roots of polynomial systems via pointed Newton polytopes. J. Complexity, 12:116–133, 1996.

[87] J.M. Rojas. Toric intersection theory for affine root counting. Journal of Pure and Applied Algebra, 136(1):67–100, 1999.

[88] J. Rosenthal. On dynamic feedback compensation and compactifications of systems. SIAM J. Control and Optimization, 32(1):279–296, 1994.

[89] J. Rosenthal and F. Sottile. Some remarks on real and complex output feedback. Systems and Control Lett., 33(2):73–80, 1998. See http://www.nd.edu/~rosen/pole for a description of computational aspects of the paper.

[90] J. Rosenthal and J.C. Willems. Open problems in the area of pole placement. In V.D. Blondel, E.D. Sontag, M. Vidyasagar, and J.C. Willems, editors, Open Problems in Mathematical Systems and Control Theory, Communication and Control Engineering Series. Springer–Verlag, London, 1999.

[91] H. Schubert. Beziehungen zwischen den linearen Räumen auferlegbaren charakteristischen Bedingungen. Math. Ann., 38:588–602, 1891.

[92] B. Segre. The Non-Singular Cubic Surfaces. A New Method of Investigation with Special Reference to Questions of Reality. Oxford University Press, 1942.

[93] A.J. Sommese and C.W. Wampler. Numerical algebraic geometry. In J. Renegar, M. Shub, and S. Smale, editors, The Mathematics of Numerical Analysis, volume 32 of Lectures in Applied Mathematics, pages 749–763, 1996. Proceedings of the AMS-SIAM Summer Seminar in Applied Mathematics, Park City, Utah, July 17-August 11, 1995, Park City, Utah.

[94] A.J. Sommese and J. Verschelde. Numerical homotopies to compute generic points on positive dimensional algebraic sets. MSRI Preprint #1999-038, 1999.

[95] M. Sosonkina, L.T. Watson, and D.E. Stewart. Note on the end game in homotopy zero curve tracking. ACM Trans. Math. Softw., 22(3):281–287, 1996.

[96] F. Sottile. Enumerative geometry for real varieties. In J. Kollár, R. Lazarsfeld, and D. R. Morrison, editors, Algebraic Geometry - Santa Cruz 1995 (University of California, Santa Cruz, July 1995), volume 62, Part I of Proceedings of Symposia in Pure Mathematics, pages 435–447. AMS, 1997.

[97] F. Sottile. Pieri’s formula via explicit rational equivalence. Can. J. Math., 49(6):1281–1298, 1997.

[98] F. Sottile. Real Schubert calculus: polynomial systems and a conjecture of Shapiro and Shapiro. MSRI Preprint #1998-066, 1998. Accepted for publication in Experimental Mathematics.

[99] F. Sottile. The special Schubert calculus is real. Electronic Research Announcements of the AMS, 5:35–39, 1999. Available at http://www.ams.org/era.

[100] F. Sottile. Real rational curves in Grassmannians. MSRI Preprint #1999-025, 1999.

[101] B. Sturmfels. Algorithms in Invariant Theory. Texts and Monographs in Symbolic Computation. Springer–Verlag, New York, 1993.

[102] B. Sturmfels. On the Newton polytope of the resultant. Journal of Algebraic Combinatorics, 3:207–236, 1994.

[103] B. Sturmfels. Gröbner Bases and Convex Polytopes, volume 8 of University Lecture Series. AMS, Providence, R.I., 1996.

[104] B. Sturmfels. Polynomial equations and convex polytopes. Amer. Math. Monthly, 105(10):907–922, 1998.

[105] J. Verschelde, M. Beckers, and A. Haegemans. A new start system for solving deficient polynomial systems using continuation. Appl. Math. Comput., 44(3):225–239, 1991.

[106] J. Verschelde and R. Cools. Symbolic homotopy construction. Applicable Algebra in Engineering, Communication and Computing, 4(3):169–183, 1993.

[107] J. Verschelde and R. Cools. Symmetric homotopy construction. J. Comput. Appl. Math., 50:575–592, 1994.

[108] J. Verschelde and K. Gatermann. Symmetric Newton polytopes for solving sparse polynomial systems. Adv. Appl. Math., 16(1):95–127, 1995.

[109] J. Verschelde and A. Haegemans. The GBQ-Algorithm for constructing start systems of homotopies for polynomial systems. SIAM J. Numer. Anal., 30(2):583–594, 1993.

[110] J. Verschelde, P. Verlinden, and R. Cools. Homotopies exploiting Newton polytopes for solving sparse polynomial systems. SIAM J. Numer. Anal., 31(3):915–930, 1994.

[111] J. Verschelde. PHC and MVC: two programs for solving polynomial systems by homotopy continuation. In J.C. Faugère, J. Marchand, and R. Rioboo, editors, Proceedings of the PoSSo Workshop on Software. Paris, March 1–4, 1995, pages 165–175, 1995.

[112] J. Verschelde, K. Gatermann, and R. Cools. Mixed-volume computation by dynamic lifting applied to polynomial system solving. Discrete Comput. Geom., 16(1):69–112, 1996.

[113] J. Verschelde. Homotopy Continuation Methods for Solving Polynomial Systems. PhD thesis, K.U.Leuven, Dept. of Computer Science, 1996. Postscript version available from the author’s web pages.

[114] J. Verschelde. Numerical evidence for a conjecture in real algebraic geometry. MSRI Preprint #1998-064, 1998. Accepted for publication in Experimental Mathematics. Paper and software available at the author’s web-pages.

[115] J. Verschelde. PHCpack: A general-purpose solver for polynomial systems by homotopy continuation. Accepted for publication in ACM Trans. Math. Softw.. Paper and software available at the author’s web-pages., 1999.

[116] J. Verschelde. Toric Newton method for polynomial homotopies. MSRI Preprint #1999-005, 1999. Revised for special issue of J. Symbolic Computation.

[117] A. Wallack, I.Z. Emiris, and D. Manocha. MARS: A Maple/Matlab/C resultant-based solver. In O. Gloor, editor, Proceedings of ISSAC-98, pages 244–251, Rostock, Germany, 1998. ACM.

[118] C.W. Wampler, A.P. Morgan, and A.J. Sommese. Numerical continuation methods for solving polynomial systems arising in kinematics. ASME J. of Mechanical Design, 112(1):59–68, 1990.

[119] C.W. Wampler, A.P. Morgan, and A.J. Sommese. Complete solution of the nine-point path synthesis problem for four-bar linkages. ASME J. of Mechanical Design, 114(1):153–159, 1992.

[120] C.W. Wampler. Bezout number calculations for multi-homogeneous polynomial systems. Appl. Math. Comput., 51(2–3):143–157, 1992.

[121] C.W. Wampler. Isotropic coordinates, circularity and Bezout numbers: planar kinematics from a new perspective. Proceedings of the 1996 ASME Design Engineering Technical Conference. Irvine, CA, Aug 18–22, 1996. Available on CD-ROM.

[122] X. Wang, M.S. Ravi, and J. Rosenthal. Algebraic and combinatorial aspects of the dynamic pole assignment problem. In U. Helmke, R. Mennicken, and J. Saurer, editors, Systems and Networks: Mathematical Theory and Applications, volume 79 of Mathematical Research, pages 547–550. Akademie Verlag, Berlin, 1994. Proc. of the International Symposium MTNS 93 held in Regensburg Germany. Vol II: Invited and contributed papers.

[123] L.T. Watson. Numerical linear algebra aspects of globally convergent homotopy methods. SIAM Rev., 28(4):529–545, 1986.

[124] L.T. Watson, S.C. Billups, and A.P. Morgan. Algorithm 652: HOMPACK: a suite of codes for globally convergent homotopy algorithms. ACM Trans. Math. Softw., 13(3):281–310, 1987.

[125] L.T. Watson. Globally convergent homotopy methods: a tutorial. Appl. Math. Comput., 31(Spec. Issue):369–396, 1989.

[126] L.T. Watson, M. Sosonkina, R.C. Melville, A.P. Morgan, and H.F. Walker. HOMPACK90: A suite of Fortran 90 codes for globally convergent homotopy algorithms. ACM Trans. Math. Softw., 23(4):514– 549, 1997. Available at http://www.cs.vt.edu/˜ ltw/.

[127] A. Weil. Foundations of Algebraic Geometry, volume XXIX of AMS Colloquium Publications. AMS, Providence, Rhode Island, 1962. Revised and Enlarged Edition.

[128] S.M. Wise, A.J. Sommese, and L.T. Watson. POLSYS PLP: A partitioned linear product homotopy code for solving polynomial systems of equations. Available at http://www.cs.vt.edu/˜ ltw/, 1998.

[129] A.H. Wright. Finding all solutions to a system of polynomial equations. Math. Comp., 44(169):125–133, 1985.

[130] G.M. Ziegler. Lectures on Polytopes, volume 152 of Graduate Texts in Mathematics. Springer–Verlag, New York, 1995.

[131] W. Zulehner. A simple homotopy method for determining all isolated solutions to polynomial systems. Math. Comp., 50(181):167–177, 1988.

Mathematical Sciences Research Institute, 1000 Centennial Drive, Berkeley, CA 94720-5070, U.S.A.

Department of Mathematics, Michigan State University, East Lansing, MI 48824-1027, U.S.A.

  E-mail address: jan@msri.org, jan@math.msu.edu or jan.verschelde@na-net.ornl.gov, and urls:
http://www.msri.org/people/members/jan http://www.mth.msu.edu/∼jan