NUMERICAL EVIDENCE FOR A CONJECTURE IN REAL ALGEBRAIC GEOMETRY

JAN VERSCHELDE

Abstract. Homotopies for polynomial systems provide computational evidence for a challenging instance of a conjecture about whether all solutions are real. The implementation of SAGBI homotopies involves polyhedral continuation, flat deformation and cheater’s homotopy. The numerical difficulties are overcome if we work in the true synthetic spirit of the Schubert calculus by selecting the numerically most favorable equations to represent the geometric problem. Since a well-conditioned polynomial system allows perturbations on the input data without destroying the reality of the solutions we obtain not just one instance, but a whole manifold of systems that satisfy the conjecture. Also an instance that involves totally positive matrices has been verified. The optimality of the solving procedure is a promising first step towards the development of numerically stable algorithms for the pole placement problem in linear systems theory. AMS Subject Classification. 14M15, 65H10, 68Q40, 93B55. Keywords. polynomial system, numerical Schubert calculus, SAGBI homotopy, polyhedral continuation, flat deformation, cheater’s homotopy, pole placement problem.

1. Introduction

Solving a polynomial system numerically means that approximations are computed to all isolated solutions of the system. Having an approximate zero (as in [6]) implies that Newton’s method doubles its accuracy in each step. Homotopy continuation methods provide paths to all isolated approximate zeros. The references [26], [22] and [9] treat polynomial homotopies respectively from within the fields of engineering, numerical analysis and computational algebraic geometry. Path-following methods are described in [1, 2].

Optimal homotopies for solving polynomial systems arising in the Schubert calculus of enumerative geometry were proposed by Birk Huber, Frank Sottile and Bernd Sturmfels in [18]. These homotopies are optimal in the sense that every path leads to a solution when applied to a generic problem instance, whereas the standard homotopies force to trace many diverging solution paths.

The type of polynomial system that needs to be solved is presented in the next section, followed by a survey on standard root-counting methods. Thereafter come implementational aspects for the homotopies and a derivation of the equations in the pole placement problem. A report on the main numerical difficulties and solutions is given in the sixth section. The last part of the paper contains a short description of the freely available software package PHC developed by the author. Execution times are listed illustrating the performance of the methods.

Research at MSRI is supported by a post-doctoral fellowship. Supported in part by NSF under Grant DMS - 9804846 at Michigan State University.

2. Testing a conjecture of Boris and Michael Shapiro

Our problem instance is encoded as (m, p) = (2, 8), see [14], [30], [33] and [34] for other cases and related work. The formulation used here is as a problem in enumerative geometry.

Given 16 2-planes, osculating a rational normal curve, we look for 8-planes in ℂ^{10} that meet these 16 given 2-planes nontrivially. The intersection is nontrivial if and only if

f_{i}(x) = det [ \begin{matrix}1 & 0 & x_{11} & x_{12} & ⋯ & x_{18} \\ s_{i} & 1 & x_{21} & x_{22} & ⋯ & x_{28} \\ s_{i}^{2} & 2s_{i} & 1 & 0 & ⋯ & 0 \\ s_{i}^{3} & 3s_{i}^{2} & 0 & 1 & ⋯ & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ... & ⋮ \\ s_{i}^{9} & 9s_{i}^{8} & 0 & 0 & ⋯ & 1 \end{matrix} ] = 0, i = 1, 2, …, 16.  (1)

The given 2-planes are spanned by the first two columns of (1). The second column is the derivative in the first one. This makes it tangent to a space curve, whence we call it ‘osculating’.

The last eight columns of (1) span the unknown 8-planes. Local coordinates are chosen fixing the lower-right block of the matrix as the identity matrix. The 16 unknowns x_{ij} are determined by 16 equations that come from choosing 16 distinct s_{i}-values and expanding the determinant (1). This distinct choice implies to have a generic problem instance with a finite number of solutions.

Schubert computed in [32] that there are generically 1,430 solutions to this problem. The Shapiro-Shapiro conjecture is that all solutions to the polynomial system (1) are real. Another version of this conjecture involving totally positive matrices is described at the end of section six.

3. Why we need new homotopies

The standard homotopies perform poorly when applied directly to our problem because the bounds on the number of solutions are too high. Next we review various standard rootcounting methods. These give a priori bounds which determine the number of paths in the homotopies.

Definition 3.1. The total degree D of F (x) = 0 is D = ∏_{i=1}^{n} deg(f_{i}).

Theorem 3.1. (Bézout) F (x) = 0 has no more than D isolated complex solutions, counted with multiplicities.

In [26], the theorem is proven by differential geometry. A proof based on projective elimination theory is in [8]. In [6] we find another proof with a treatment of the complexity.

Expanding (1) we obtain 16 quadratic polynomials. Applying Bézout’s theorem gives D = 2^{16} = 65, 536. A homotopy based on this bound requires thus the tracing of 65,536 solution paths.

The application of multi-homogeneous homotopies was introduced in [27] and has been applied with great success to various problems in mechanism design, see e.g. [43, 44].

Definition 3.2. Denote {x_{1}, x_{2}, …, x_{n}} by X. Consider S ⊆ X. The degree deg(f, S) of f in S is deg(f(x|_{x_{i}=1,i∉S})).

Definition 3.3. Let Z = {Z_{1}, Z_{2}, …, Z_{k}} be a partition of X. For F, the degree matrix M ∈ ℕ^{n×k} is defined as M_{ij} := deg(f_{i}, Z_{j}). The k-homogeneous Bézout number B is the permanent of M.

The permanent of the degree matrix M counts the number of solutions of a random linearproduct start system. Every equation is the product of #Z factors. The jth factor for the ith equation is the product of M_{ij} linear equations in the unknowns of Z_{j}, with random coefficients.

Theorem 3.2. For any partition Z of X, k = #Z. The corresponding k-homogeneous Bézout number B bounds the number of isolated complex solutions of F (x) = 0, counted with multiplicities.

A multi-homogeneous homotopy [27] reaches all isolated solutions of a polynomial system. Different partitions give different Bézout numbers B. A partition with the minimal B is

hard to find in general. See [42] for exhaustive search algorithms. In applications the optimal partition follows from the meaning of the variables. Here we naturally group according to rows or columns. For the (2,8)-case, the partition {{x_{11}, x_{12}, …, x_{18}}, {x_{21}, x_{22}, …, x_{28}}}, gives B = 12, 870.

Polyhedral methods have led to a computational breakthrough in solving polynomial systems. See [36] for developments and open problems on counting roots with polytopes.

Definition 3.4. Consider f(x) = ∑_{a∈A} c_{a}x^{a}, with c_{a} ∈ ℂ and x^{a} = x_{1}^{a_{1}}x_{2}^{a_{2}} ⋯ x_{n}^{a_{n}}. The set A = { a ∈ ℤ^{n} | c_{a} ≠ 0 } is the support of f. The convex hull of A is the Newton polytope of f.

Definition 3.5. Let P = (P_{1}, P_{2}, …, P_{n}) be a tuple of polytopes. The mixed volume V_{n}(P) of P is

V_{n}(P) = ∑_{I⊆{1,2,... ,n}}(−1)^{n−#I}vol_{n} (∑_{i∈I} P_{i}) ,  (2)

where vol_{n} is the usual volume.

Theorem 3.3. (Bernshteǐn) A system F (x) = 0 with Newton polytopes P has no more than V_{n}(P) isolated complex solutions in (ℂ^{∗})^{n}, ℂ^{∗} = ℂ \ {0}, counted with multiplicities.

The proof in [4] is constructive, and was implemented in [41]. See also [9].

We see from (1) that all supports are equal, so P = P_{i} for i = 1, 2, …, n, and then V_{n}(P) = n!vol_{n}(P), which is the case of Kushnirenko’s theorem [21]. The direct application of Bernshteǐn’s theorem leads to a homotopy with vol_{n}(P) = 12, 862 solution paths.

For this (2,8)-instance, only 1,430 solution paths converge using homotopies based on the above root counts, as predicted by the formula Schubert derived in [32]. Table 1 summarizes the performance of the root counts. More tables of this kind are in [33].

If p ≤ m, separating the unknowns that belong to different rows gives the minimal B. Otherwise, we better group unknowns that belong to the same columns. For m = 2, B is remarkably close to V. For m = 3, the gap widens. The calculation of B requires a permanent computation which is a #P -hard problem. Despite this, the algorithms in [42] outperform the volume computation.

m = 2  m = 3
p 2 3 4 5 6 7 8 9 2 3 4 5
D 16 64 256 1,024 4,096 16,384 65,536 262,144 64 19,683 531,441 14,348,907
B 6 20 70 252 924 3,432 12,870 48,620 20 1,680 34,650  756,756
V 4 17 66 247 918 3,425 12,862 48,611 17  642 22,148  ?
N 2 5 14  42 132  429 1,430 4,862 5  42  462  6,006

Table 1. Standard root counts and the actual number of roots N for various (m, p)-systems.

4. SAGBI Homotopies

The SAGBI homotopies are one of the three homotopies proposed in [18]. Here we summarize the method emphasizing the algorithmic aspects. The general principles are illustrated on a running example. At each stage characteristics in solving the (2,8)-case are listed.

To implement SAGBI homotopies we concatenate three different homotopies, see Figure 1. Table 2 defines the polynomial systems at each stage. Complex arithmetic in this concatenation is necessary because of singularities that otherwise may occur. In solving several real problem instances, we recycle the solutions of one generic complex problem using cheater’s homotopy.

Binomial Systems polyhedral homotopy

Generic Complex System

flat deformation

Generic Complex Problem cheater’s homotopy

Specific Real Problem

Figure 1. Concatenation of three homotopies. The central part in the chain is the SAGBI homotopy.

Binomial System  : exactly two monomials in every equation
Generic Complex System : coefficients are randomly chosen complex numbers
Generic Complex Problem : problem with randomly chosen complex planes
Specific Real Problem  : real osculating planes as in the conjecture

Table 2. The end points in the homotopies are four different types of polynomial systems.

Example 4.1. (Running example.) The system for (m, p) = (2, 2) consists of four equations:

f_{i}(x) = det [ \begin{matrix}1 & 0 & x_{11} & x_{12} \\ s_{i} & 1 & x_{21} & x_{22} \\ s_{i}^{2} & 2s_{i} & 1 & 0 \\ s_{i}^{3} & 3s_{i}^{2} & 0 & 1 \end{matrix} ] = 0, i = 1, 2, 3, 4.  (3)

To create a specific real problem we choose four different values for s_{i} and apply Laplace expansion in terms of 2-by-2 minors. The polynomials in the SAGBI homotopy are as follows:

f̂_{i}(x, t) = C_{34}^{(i)}(x_{11}x_{22} − t x_{12}x_{21}) − C_{24}^{(i)}(−x_{12}) + C_{23}^{(i)}x_{11} + C_{14}^{(i)}(−x_{22}) − C_{13}^{(i)}x_{21} + C_{12}^{(i)} = 0, (4)

for i = 1, 2, 3, 4. The lower indices in the coefficients C_{kl}^{(i)} refer to the choice of rows of the matrix in the elaboration of (3). The SAGBI homotopy starts at t = 0 and ends at t = 1. At t = 1 the generic complex problem has coefficients C_{kl}^{(i)} which are 2-by-2 minors of a 4-by-2 matrix of random complex coefficients.

Our running example is a classical problem in enumerative geometry. Note that 2-planes in ℂ^{4} correspond to lines in ℙ^{3}. The question is then to find all lines that intersect four given lines nontrivially. In general, two lines satisfy this condition, see Figure 2.

[figure: .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .] Figure 2: m = 2 = p.

The positive real orthant of projective 3-space corresponds to the interior of the tetrahedron pictured in Figure 2. In Figure 2, the four given lines are drawn in solid fine lines. The two thick lines meet those four given lines in a point. Note that this problem is fully real: both the input and output lines are real lines. The special position of two of the four input lines is adapted from the Pieri homotopy algorithm [18]. Another choice of local coordinates, unlike those in (3), for the output planes enables to solve this instance.

4.1. Newton Polytopes and Polyhedral Continuation. To compute the mixed volume we need a subdivision of the polytopes. A regular subdivision is obtained by lifting the polytopes into a higher-dimensional space and projecting the facets of the lower hull of the sum of the lifted polytopes down to the original space. The key idea of [19] is that this lifting induces the so-called polyhedral homotopy, starting at binomial systems whose Newton polytopes match the cells in the subdivision. We refer to [22] for numerical aspects.

For our problem we compute the volume of one Newton polytope, for which the dynamic lifting algorithm [40] is well-suited. This algorithm incrementally updates the triangulation each time selecting the lowest possible lifting value for the added point. The induced homotopies have a minimal power in the continuation parameter. For the (2,8)-case, the highest lifting value is 133.

Example 4.2. (Example 4.1 continued.) The polyhedral homotopy that is induced by the dynamic lifting algorithm consists of the equations

ĝ_{i}(x, t) = c_{1}^{(i)}x_{11}x_{22} + c_{2}^{(i)}x_{11}t + c_{3}^{(i)}x_{12} + c_{4}^{(i)}x_{21} + c_{5}^{(i)}x_{22} + c_{6}^{(i)} = 0, i = 1, 2, 3, 4. (5)

The coefficients c_{j}^{(i)} are random complex constants. The triangulation has two cells and correspondingly there are two homotopies. The first homotopy equals (5). At t = 1 we have the system we want to solve. At t = 0 we have a system of four equations, each with exactly five terms. By Gaussian elimination we obtain a binomial system that can be solved efficiently [41]. So at t = 0 we know the solutions and we can start following the paths defined by the homotopy (5).

Ordering the variables as (x_{11}, x_{12}, x_{21}, x_{22}, t), the support of (5) is the lifted point configuration Â = {(1, 0, 0, 1, 0), (0, 1, 0, 0, 0), (1, 0, 0, 0, 1), (0, 0, 0, 1, 0), (0, 0, 1, 0, 0), (0, 0, 0, 0, 0)}. The other cell in the triangulation consists of the points â ∈ Â for which the inner product ⟨â, v⟩ with v = (−1, 0, 0, 1, 1) is minimal. To set up the second homotopy we replace x^{a} by x^{a}t^{⟨â,v⟩}:

ĥ_{i}(x, t) = c_{1}^{(i)}x_{11}x_{22} + c_{2}^{(i)}x_{11} + c_{3}^{(i)}x_{12} + c_{4}^{(i)}x_{21} + c_{5}^{(i)}x_{22}t + c_{6}^{(i)} = 0, i = 1, 2, 3, 4. (6)

At t = 1 we have the system we want to solve. At t = 0 we can again reduce the system to a binomial system whose solutions are the starting points of the paths defined by the homotopy (6).

We have found all roots for the system we wanted to solve since (0, 0, 0, 0, 1) and (−1, 0, 0, 1, 1) are the only two vectors (unique upon scaling) for which the selection of points with minimal inner product leads to a cell with nonzero volume and a homotopy that has solutions in (ℂ^{∗})^{4} at t = 0.

The triangulation for the (2, 8)-case consists of 1,430 simplices because the Newton polytope is the order polytope [38] which is unimodular.

4.2. Flat Deformations arising from Gröbner and SAGBI Bases. The first homotopy of the three presented in [18] is the Gröbner homotopy. Although its efficiency is inferior compared to the SAGBI homotopy, it helps understanding the latter.

Gröbner basics [38] teach us that for any term order ≺ and ideal I, there exists a weight vector w ∈ ℕ^{n} defining a monomial order ≺_{w} (as: x^{a} ≺_{w} x^{b} ⇔ ⟨a, w⟩ > ⟨b, w⟩ or ⟨a, w⟩ = ⟨b, w⟩ and x^{a} ≺_{lex} x^{b}) that induces the same initial ideal: in_{≺}(I) = in_{≺_{w}}(I). In [13], given a Gröbner basis, a homotopy is defined to deform the zero set of the ideal I into the zero set of in_{≺_{w}}(I). In particular, the parameter t in the flat deformation is introduced substituting each x_{i} by x_{i}t^{w(x_{i})}, where w(x_{i}) is the weight of x_{i}. Hereby x^{a} is replaced by x^{a}t^{⟨a,w⟩}. All systems in that homotopy have the same initial ideal, whence the same Hilbert function, whence the same structure of the zero set ∀t.

The unknowns in the Gröbner homotopy are p-by-p minors, denoted by brackets. A bracket is a vector of indices to the rows that have been selected in constructing the minor. The monomials are products of brackets. We see monomials as tableaux that have brackets in their rows. A tableau is standard if its columns are sorted, otherwise it is nonstandard. The straightening algorithm provides relations to rewrite a nonstandard tableau as a linear combination of standard ones. The combinatorics of tableaux and the relation with the Schubert calculus are treated in [15].

Example 4.3. (Example of a straightening syzygy.) A 2-by-2 minor that selects the ith and jth row of a matrix, i < j, is represented by the bracket [i j]. Expanding the determinant in (7) along the first row, we see that we always obtain two equal columns. With Laplace expansion into 2-by-2 minors we obtain in bracket notation a straightening syzygy:

det [ \begin{matrix}x_{11} & x_{12} & 0 & 0 \\ x_{21} & x_{22} & x_{21} & x_{22} \\ x_{31} & x_{32} & x_{31} & x_{32} \\ x_{41} & x_{42} & x_{41} & x_{42} \end{matrix} ] ≡ 0. ⇔ [1 2][3 4] − [1 3][2 4] + [1 4][2 3] = 0. (7)

We can rewrite the nonstandard monomial [1 4][2 3] as a linear combination of standard ones.

The defining ideal of the Grassmann manifold is generated by all straightening syzygies needed to rewrite all possible nonstandard tableaux. In [37] it is proven that these rewriting rules constitute a Gröbner basis for the term order that selects the monomials corresponding to the nonstandard tableaux as the leading ones. In [18] we find a recipe to assign weights and to set up the flat deformation. The initial ideal is squarefree so that the start solutions are regular. Except for coefficients that belong to an algebraic set, all solutions are regular.

Example 4.4. (Example 4.1 continued.) The Gröbner homotopy consists of five equations, homogeneous in the six brackets:

\{ \begin{matrix}[1 4][2 3] − [1 3][2 4]t + [1 2][3 4]t^{2} = 0 \\ C_{34}^{(i)}[1 2] − C_{24}^{(i)}[1 3] + C_{23}^{(i)}[1 4] + C_{14}^{(i)}[2 3] − C_{13}^{(i)}[2 4] + C_{12}^{(i)}[3 4] = 0, i = 1, 2, 3, 4. \end{matrix} (8)

The first equation comes from (7). The coefficients C_{kl}^{(i)} are minors selecting the kth and lth row of a 4-by-2 matrix with random complex coefficients. At t = 0, there are two start solutions.

Note that the parameter t only appears in the equation that defines the Grassmann manifold. The linear equations of the intersection condition remain invariant under the flat deformation.

In the (2, 8)-case, the coefficients are 2-by-2 minors and the unknowns are brackets representing 8-by-8 minors, ( \begin{matrix}10 \\ 8 \end{matrix} ) = 45 unknowns in total. Expanding (1) leads to a 16-by-16 linear system in the brackets constrained by 210 quadratic equations that define Grassmann manifold. The Gröbner homotopy solves the (2,8)-problem as an overdetermined system of 226 equations in 45 unknowns.

When the minors are expanded in terms of the x_{ij}’s the equations that define the Grassmann manifold are automatically satisfied. Therefore we are left with fewer equations than in the Gröbner homotopy, which explains the superior efficiency of the SAGBI homotopy. In particular, for the (2,8)-case we are left with 16 equations in 16 unknowns.

After the expansion in the x_{ij}’s we no longer have an ideal but polynomials that form a subalgebra in the polynomial ring ℂ[x_{ij}]. SAGBI is the acronym for Subalgebra Analogue to Gröbner Bases for Ideals. The term order we use selects as leading monomial the product of elements on the diagonal of the matrices. In [38] flat deformations are extended to SAGBI bases. The general recipe for the flat deformation in [18] is given by substituting x_{ij} by x_{ij}t^{(i−1)(p−j)} and dividing out the lowest power of t, minor per minor. The leading monomials with this term order have minimal powers of t and are products of diagonal elements in the minors.

Example 4.5. (Example 4.1 continued.) The SAGBI homotopy (4) in determinantal form is

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

where the coefficients c_{kl}^{(i)} are random complex constants. The brackets are expanded as follows:

\begin{aligned}[1 2] & ← & x_{11}x_{22} − x_{21}x_{12} t \\ [1 3] & ← & −x_{12} \\ [1 4] & ← & x_{11} \end{aligned}  \begin{aligned}[2 3] & ← & −x_{22} \\ [2 4] & ← & x_{21} \\ [3 4] & ← & 1 \end{aligned}  (10)

Substituting the brackets into the linear equations of (8) gives the SAGBI homotopy (4).

The SAGBI homotopies are nonlinear in t. For the (2, 8)-case, the highest power of t equals 7.

4.3. The Cheater’s Homotopy. In [23] the so-called cheater’s homotopy was presented to solve repeatedly a polynomial system with coefficients as functions of parameters for several instances of these parameters. The procedure assumes that one has solved the polynomial system once – this is the cheating part – for a generic complex choice of the parameters. See [28] for a similar idea.

The parameters in our problem are the minors from the m-planes. To avoid repeated evaluation of minors we apply the result of [24]. The start system is a problem instance for a generic choice of the parameters c^{0} ∈ ℂ^{k}. For any c ∈ ℂ^{k}, the following homotopy guarantees to reach all solutions of F (c, x) = 0:

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

This homotopy avoids the evaluation of the parameters during path following.

Example 4.6. (Example 4.1 continued.) Take the homotopy H(x, t) = 0 in (11). At t = 0, the parameters c^{0} consist of all 2-by-2 minors C_{kl}^{(i)} of a 4-by-2 random complex matrix. This system has been solved by the SAGBI homotopy. At t = 1, the parameters c consist of the maximal minors of the first two columns of the matrix in (3). The two solution paths converge to two distinct real solutions that span 2-planes intersecting the two given 2-planes nontrivially.

With cheater’s homotopy, the conjecture can be tested systematically, deforming so that every system in the family is a real problem instance. If singular solutions occur where two real solution paths join into a complex conjugate pair of solutions, then the conjecture is false.

5. The Pole Placement Problem

In this section we describe how to translate the pole placement problem of linear systems theory (see [7] and [29]) into the geometric formulation (1).

Consider a system with input u ∈ ℝ^{m}, output y ∈ ℝ^{p}, and internal states x ∈ ℝ^{n}, whose evolution in time t is governed by the first-order linear differential equation

\{ \begin{matrix}ẋ(t) = Ax(t) + Bu(t), & A ∈ ℝ^{n×n}, B ∈ ℝ^{n×m}, \\ y(t) = Cx(t), & C ∈ ℝ^{p×n}, \\ u(t) = F y(t), & F ∈ ℝ^{m×p}, \end{matrix}  (12)

where the last part expresses the control of the input by constant output feedback F.

Substitution into the first equation of (12) yields ẋ(t) = (A + BF C)x(t), whose characteristic polynomial is

φ(s) = det(sI_{n} − A − BF C).  (13)

The roots s_{i}: φ(s_{i}) = 0, for i = 1, 2, …, n, are the natural frequencies of the controlled system.

The pole placement problem is an inverse problem: given A, B, C, and φ (determined by the frequencies s_{i}), compute the feedback laws F that satisfy φ(s_{i}) = 0, for i = 1, 2, …, n.

We rewrite the expression (13) for φ(s) as follows:

det(sI_{n} − A − BF C) = 0 ⇔ det [ \begin{matrix}sI_{n} − A & 0 & −B \\ −C & I_{p} & 0 \\ 0 & −F & I_{m} \end{matrix} ] = 0  (14)

⇔ det [( 0 −F I_{m} ) ( \begin{matrix}X(s) \\ Y (s) \\ U(s) \end{matrix} )] = 0  (15)

⇔ det [ \begin{matrix}U(s) & F \\ Y (s) & I_{m} \end{matrix} ] = 0  (16)

where I_{p} and I_{m} are the respective p-by-p and m-by-m identity matrices. The first step in (14) uses a determinantal identity. Secondly, we apply the fact that

det [ \begin{matrix}P (s) \\ M(s) \end{matrix} ] = c det [M(s)Q(s)] , c ∈ ℝ,  (17)

for

P (s) = [ \begin{matrix}sI_{n} − A & 0 & −B \\ −C & I_{p} & 0 \end{matrix} ] , Q(s) = [ \begin{matrix}X(s) \\ Y (s) \\ U(s) \end{matrix} ] and P (s)Q(s) = 0.  (18)

The third step is a simple elaboration to transform the equations of the pole placement problem into the familiar geometric form (16) we have used to solve it.

6. Symbolic-Numeric Calculations

Singularities do not occur for generic complex problems. To solve the real problem, we are choose orthogonal representations for the planes and obtain well-conditioned polynomial systems.

6.1. Localizing before Expanding, Expanding before Cheating. The equations (1) are already written in local coordinates. But any 8-by-8 submatrix using the last 8 columns can serve as unit matrix. Setting up the equations we better first fix this choice before expanding the determinants. Otherwise, see what happens in the (2,8)-case: every 8-by-8 minor has in its expansion 8! = 40, 320 terms and there are 45 such minors.

One can implement the cheater’s homotopy in two ways. Either one can introduce the continuation parameter t inside the matrix (1) or one can do this after the expansion of the minors. The first way [23] leads to homotopies that require the elaboration of determinants each time evaluation is needed. The second way [24] uses (11) which does not require the evaluation of the parameters which leads to a much more efficient cheating procedure.

To illustrate this last point, the determinantal cheater’s homotopy for the case (4,2) requires 6 min 24 sec CPU time versus 30 sec for the homotopy which avoids the evaluation of determinants in every step.

6.2. Use Chebychev Polynomials as Basis. Generating values for s in equation (1) quickly leads to huge numbers when s > 1 or very tiny numbers when s < 1. The calculation of minors in the expansion of (1) amplifies this effect. We cannot take all values for s too close to 1 because then the input planes lie too close to each other.

To approximate a function by a polynomial with minimal error the basis of Chebychev polynomials (see e.g. [16]) is commonly used instead of the basis of standard monomials {1, x, x^{2}, …, x^{n}}. The nth degree Chebychev polynomial is defined by T_{n}(x) = cos(n arccos(x)), whence T_{n}(x) ∈ [−1, +1], for x ∈ [−1, +1]. Therefore we obtain nice numbers for the coefficients in the polynomials generated by (1). We show the equivalence of problem formulations by example.

Example 6.1. (Example 4.1 reformulated.) Multiplying the 4-by-2 matrix by a nonsingular matrix gives an equivalent representation for the 2-plane. As result of the following multiplication

[ \begin{matrix}1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ −1 & 0 & 2 & 0 \\ 0 & −3 & 0 & 4 \end{matrix} ] [ \begin{matrix}1 & 0 \\ s & 1 \\ s^{2} & 2s \\ s^{3} & 3s^{2} \end{matrix} ] = [ \begin{matrix}1 & 0 \\ s & 1 \\ 2s^{2} − 1 & 4s \\ 4s^{3} − 3s & 12s^{2} − 3 \end{matrix} ]  (19)

we recognize the first four Chebychev polynomials and their derivatives.

We stress that this modification of the system (1) does not change the problem, but only its representation as polynomial system. Actually, this reformulation gives rise to a wider class of polynomial systems whose solutions are all expected to be real.

Polynomials in the Chebychev basis have more real roots on average than polynomials in the standard basis. In [5] we read that a polynomial of degree d with independent, normally distributed random coefficients written in the Chebychev basis is expected to have \frac{d}{\sqrt{3}} real roots, as d → ∞, whereas this asymptotic bound is only \frac{2}{π} log(d) for a random polynomial denoted in the standard basis. See [10] for more on the implication of the number of average real roots on complexity issues in computer algebra. In our case, Chebychev polynomials lead to a numerically better formulation of a problem that is expected to have all its solutions real.

6.3. Orthonormal Representations of the Planes. Using Chebychev polynomials suffices in case m = 2, because the function values are bounded. However, it is hard to get bounded derivatives. In the (8, 2)-case we must differentiate 7 times and compute 8-by-8 minors with the results of this differentiation. Without reformulation, only the cases where m = 2 are numerically tractable.

The column vectors in orthogonal matrices form an orthogonal basis. Orthonormal matrices are the best conditioned representations of the input planes. The QR factorization [17] of the matrix A ∈ ℂ^{m×n} rewrites A as the product

A = QR, Q ∈ ℂ^{m×m} and R ∈ ℂ^{m×n},  (20)

where Q is orthonormal and R upper triangular. If A has full column rank, then the first n columns of Q form an orthonormal basis for the space spanned by the columns of A. This factorization is computed by the QR algorithm [17]. Adopting this representation, there is no difference between the cases (3, 4) and (4, 3).

Carrying this much further would be to develop a Newton’s method for the Grassmann manifold, which would – instead of fixing the localization in advance – exploit the freedom of additional variables to work with orthogonal bases for the p-planes throughout the path following. As mentioned in [12], Newton’s method on the Grassmann manifold has never been studied before explicitly. In our case we would have to approximate the Jacobian matrix to avoid dealing with the huge expanded polynomials before the localization. This is an interesting application of the secant method for which an α-theory has recently been developed [46].

6.4. Total Positivity. Another conjecture of Boris and Michael Shapiro involves totally positive matrices, and is also described and tested in [33]. A real matrix is totally positive [3] if all its minors are positive. Upper triangular matrices are called totally positive when all minors are positive, except for those minors which vanish on all upper triangular matrices. For all upper triangular totally positive matrices T considered here we assume T_{ii} = 1.

The sequence of m-planes proposed in this Shapiro-Shapiro conjecture for which all solutions are expected to be real is defined as follows. The first m-plane can be spanned by any (m + p)-by-m matrix. The next m-plane is generated by the last m columns of the product of the previous m-plane with any (m+p)-by-(m+p) random totally positive matrix. Repeat this generation until mp m-planes are obtained.

Since totally positive matrices have all their entries positive, a multiplication with a totally positive matrix increases the size of the numbers. Geometrically, the limiting position of the planes in the sequence defined above corresponds to moving from the plane spanned by the first m standard basis vectors to the so-called opposite flag, spanned by the last m standard basis vectors. For example, in Figure 2, we see the first and last line in the sequence as those lines spanned by the corners of the tetrahedron. With a localization as in (3) the intersection condition for the opposite flag is

det [ \begin{matrix}1 & 0 & x_{11} & x_{12} \\ 0 & 1 & x_{21} & x_{22} \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{matrix} ] = 0 ⟺ 1 = 0.  (21)

This geometric interpretation exhibits the extremal cases of numerical troubles with this localization. So we use a localization that has in every row and column at least one free variable to represent any point of intersection. The expansion leads to polynomials of degree p.

Whitney’s reduction theorem [45], in the form of Loewner [25], is used to generate the upper triangular totally positive matrices. For every n-by-n matrix we generate \frac{n(n−1)}{2} positive random numbers. Tuning the size of those random numbers allows us to control the speed by which the sequence moves to the opposite flag. There is the following tradeoff. A selection of small random numbers leads to input planes that lie too close to each other, whereas with large numbers we end up with a cluster of planes too close to the opposite flag. The actual input to the cheater’s homotopy consists in the orthonormal part of the QR-factorization applied to the generated sequence of m-planes.

The case (3,4) verified by the SAGBI homotopies has total degree D = 4^{12} = 16, 777, 216 and 4-homogeneous Bézout number B = 369, 600 with the partition {{x_{11}, x_{51}, x_{71}}, {x_{22}, x_{62}, x_{72}}, {x_{33}, x_{63}, x_{73}}, {x_{44}, x_{64}, x_{74}}}, whereas we know there are only 462 solutions. All 462 solutions are different from each other, with condition numbers less than 10^{8}, and as conjectured, are all real. The components of the solution p-planes all have the same sign pattern. The matrices that represent those solutions are not totally positive, although all corresponding minors have the same sign.

7. Available Software and Execution Times

This section describes programming and computational experiences with the author’s software package PHC (Polynomial Homotopy Continuation), whose first version was publicly released in [39].

7.1. A new driver added to PHC. Homotopy continuation methods consist of two parts:

1. Bounds for the number of roots are obtained by applying the theorems of Bézout or Bernshteǐn. For sparse systems, the mixed volume yields a much sharper bound and polyhedral homotopies are invoked to solve a start system with random coefficients. For general dense systems or when the degrees are used to count the roots, the construction of a start system does not require continuation.

2. The path tracker exploits the fact that solution paths of polynomial homotopies do not turn back, so that the continuation parameter is fixed during the correction stage. Path crossing is avoided enforcing Newton’s method to converge quadratically when correcting the solutions. At the end of the solutions paths, condition numbers are computed and mathematical certificates of divergence [20] can be obtained.

The implemented homotopy methods are powerful and reliable enough to operate as a general-purpose solver in black-box mode for small and medium-sized problems. This is illustrated in [39] on the test suite of about eighty polynomial systems from a wide variety of sources that are in the distribution of the package. Of course, given the poor performance of the standard root counts, the black-box solver is not recommended for the (2,8)-problem.

The SAGBI homotopies are implemented in a separate module of PHC, which generates the equations, calls the dynamic lifting algorithm and the path-tracking routines. Besides the driver, the new implementation consists of Laplace expansion, the assignment of the powers of the continuation parameter in the flat deformation, and the generation of the osculating planes. The new public version of PHC offers this driver when calling phc -e. Another new feature in PHC is the multi-precision arithmetic to evaluate and refine the solutions of a polynomial system.

7.2. Computational experiences. In Table 3 the timing results for the main computation are summarized. Calculations are done on a 166 MHz Pentium II processor with 64 Mb internal memory running Linux.

stages of computation  (4, 3)  (2, 8)  (3, 4)TP
regular triangulation  0h 0m 7s 250ms 0h 0m 56s 860ms  0h 0m 6s 910ms
polyhedral continuation 1h 21m 21s 630ms 8h 1m 37s 330ms 1h 32m 28s 790ms
flat deformation 1h 2m 6s 0ms 4h 24m 49s 590ms 3h 51m 8s 50ms
cheater’s homotopy  2h 42m 15s 340ms 13h 23m 7s 330ms 13h 33m 29s 410ms
total time for solving  5h 6m 21s 110ms 25h 52m 55s 150ms 18h 57m 13s 160ms

Table 3. User CPU times for solving the cases (4, 3), (2, 8), and (3,4) with total positivity.

In Table 3 we see that the most expensive stage is the cheater’s homotopy to a specific real problem instance. The (4,3)-case is given to illustrate how the computations scale up when going from a 12-dimensional to a 16-dimensional system, from 462 to 1,430 roots. Dividing total time by 1,430 paths, the cost for each path is about 1 minute user CPU time. The (3,4)-case with totally positive matrices is harder because of the fourth-degree polynomials.

The systems with their solutions are publicly available via the author’s web page. The condition numbers of the Jacobian matrices evaluated at the end of the paths range between 10^{4} and 10^{8}, which is in the worst-case scenario still sufficient to guarantee the correctness to eight decimal places. The distance between two roots is at least 10^{−8}.

Moreover, these condition numbers are complex in the sense that they reflect the sensitivity of perturbations with complex numbers. The real condition numbers must be much smaller, because the largest imaginary parts of the solutions are of order 10^{−20}. After application of Newton’s method (using QR decomposition followed by least squares approximation to solve the linear systems) ten times to each solution, the imaginary parts shrunk below 10^{−100}.

8. Conclusions

For the (2,8)-case we computed a polynomial system with all 1,430 solutions real, wellconditioned and distinct from each other. Homotopy continuation extends this real problem instance to a manifold of polynomial systems that have all their solutions real. The results of the calculations described in this paper are relevant for the following reasons:

1. The class of systems constitutes an interesting benchmark for polynomial system solvers. We can exploit the synthetic nature of enumerative geometry questions to choose the most favorable numerical representation of the problem.

2. Besides providing computational evidences for conjecture in real algebraic geometry, these tools and experiences may become useful to the engineering community. Providing stable numerical algorithms for the pole placement problem is one of open problems listed in [31]. This paper shows that homotopy continuation is a first promising step.

3. Almost all computations in algebraic geometry are exact as implemented in computer algebra systems. The calculations reported in this paper show that meaningful results (in the spirit of [35]) are obtained from numerical approximations.

4. The practical efficiency of homotopies for solving polynomial systems that arise from expanding determinants is demonstrated. Jointly with Birk Huber, the author is working on an implementation of the Pieri homotopies presented in [18]. How to apply the theory of [11] to these overdetermined homotopies is an interesting future research project.

Acknowledgements. I thank Birk Huber for explaining flat deformations to me in his apartment in Princeton. Frank Sottile told me about the ‘osculating normal curves’ what triggered the idea of using Chebychev polynomials, also his explanation of ‘opposite flags’ was most inspiring. Frank gave several helpful comments on the first draft of this paper. I am grateful to both Frank and Bernd Sturmfels for encouraging the implementation and calculations. I thank Marie-Françoise Roy for pointing to the property of random real polynomials in the Chebychev basis. The discussions with T.Y. Li and Xiaoshen Wang on cheater’s homotopy were most fruitful. I thank Joachim Rosenthal for writing down for me the connections (14) between the equations in the pole placement problem. Last and certainly not least, the special programs at MSRI in the fall of 1998 brought us together which largely contributed to the writing of this paper.

References

[1] 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.

[2] 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.

[3] T. Ando. Totally positive matrices. Linear Algebra. Appl., 90:165–219, 1987.

[4] 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.

[5] A.T. Bharucha-Reid and M. Sambandham. Random Polynomials. Academic Press, 1986.

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

[7] 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.

[8] 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.

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

[10] F. Cucker and M.-F. Roy. A theorem on random polynomials and some consequences in average complexity. Journal of Symbolic Computation, 10:405–409, 1990.

[11] J.P. Dedieu and M. Shub. Newton and predictor-corrector methods for over-determined systems of equations. IBM Technical Report, IBM Research Division, 1998.

[12] A. Edelman, T. Arias, and S.T. Smith. The geometry of algorithms with orthogonality constraints. SIAM J. Matrix Anal. Appl., 20(2):303–353, 1998.

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

[14] J.-C. Faugère, F. Rouillier, and P. Zimmerman. Verifying the conjecture of Shapiro and Shapiro for m, p = 3, 4. In preparation, 1998.

[15] W. Fulton. Young Tableaux, volume 35 of London Mathematical Society Student Texts. Cambridge University Press, 1997.

[16] C.F. Gerald. Applied Numerical Analysis. Addison–Wesley, Reading, Mass., 2nd edition, 1978.

[17] G. Golub and C. Van Loan. Matrix Computations. The Johns Hopkins University Press, Baltimore and London, third edition, 1996.

[18] B. Huber, F. Sottile, and B. Sturmfels. Numerical Schubert calculus. To appear in Journal of Symbolic Computation.

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

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

[21] 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.

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

[23] 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.

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

[25] C. Loewner. On totally positive matrices. Math. Zeitschr., 63:338–340, 1963.

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

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

[28] 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).

[29] J. Rosenthal and J.M. Schumacher. Realization by inspection. IEEE Transactions on Automatic Control, 42(9):1257–1263, 1997.

[30] 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.

[31] 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.

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

[33] F. Sottile. Real Schubert calculus: polynomial systems and a conjecture of Shapiro and Shapiro. Manuscript, 1998.

[34] F. Sottile. The special Schubert calculus is real. Manuscript, 1998.

[35] H.J. Stetter. Numerical polynomial algebra. Notes for the tutorial Interaction between Numerical Analysis and Computer Algebra, held at ISSAC’98, Rostock, 1998.

[36] B. Sturmfels. Polynomial equations and convex polytopes. To Appear in American Mathematical Monthly, December 1998.

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

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

[39] J. Verschelde. PHCPACK: A general-purpose solver for polynomial systems by homotopy continuation.

    Paper and software available at http://www.math.msu.edu/~jan.

[40] 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.

[41] 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.

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

[43] 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.

[44] 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.

[45] A.M. Whitney. A reduction theorem for totally positive matrices. J. d’Analyse Math., pages 88–92, 1952.

[46] J.C. Yakoubsohn. Finding zeros of analytic functions: α-theory for secant type methods. To appear in J. Complexity.

Mathematical Sciences Research Institute, 1000 Centennial Drive, Berkeley, CA 94720-5070, USA.

Department of Mathematics, Michigan State University, East Lansing, MI 48824-1027, USA.

  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