Lectures in Applied Mathematics Volume 00, 19xx

Numerical Algebraic Geometry

Andrew J. Sommese and Charles W. Wampler

Abstract. We discuss the well developed theory for finding isolated complex solutions of polynomial systems, and show how it leads to probabilistic algorithms for analyzing the full possibly positive dimensional solution sets of polynomial systems.

Contents

Introduction 1 1. A brief review of algebraic sets 2 2. Finding isolated solutions by continuation 5 3. Probabilistic algorithms about algebraic sets 7 4. Closing remarks and open problems 13 References 15

Introduction

There is a well developed theory for finding all smooth or at worst isolated solutions of a system of N polynomial equations on ℂ^{N} which is “small” in an appropriate sense. Almost all papers on finding the numerical solutions of a system of polynomials assume that we have the same number of equations as unknowns. This restriction is not necessary and in fact obscures what we actually know how to do numerically. This wouldn’t matter so much except that in practice it very often happens that the number of equations and the number of unknowns differs. Moreover there are natural engineering problems which require us to “find” all irreducible components of an algebraic set.

In this article we don’t make the assumption that the number of equations equals the number of unknowns, nor do we restrict ourselves to smooth or isolated solutions of systems. This numerical analysis of algebraic sets we call numerical algebraic geometry. Numerical algebraic geometry is to algebraic geometry what numerical linear algebra is to linear algebra.

\begin{matrix}c \\ ◯ \end{matrix}0000 American Mathematical Society 0075-8485/00 $1.00 + $.25 per page

1991 Mathematics Subject Classification. Primary 65H10; Secondary 14Q99.

The first author was supported in part by the University of Notre Dame, General Motors Research and Development Center, and by NSF Grant DMS 93-02121.

Our algorithms are probabilistic—as are most of the fast numerical methods for solving systems of n polynomials in n variables. They hold except for a probability zero “set” of choices. This is acceptable from the point of view of engineering problems where we are trying to find all solutions and are willing to give up complete mathematical certainty that we have found all of them.

In §1 we discuss algebraic sets and their decomposition into irreducible algebraic sets. We discuss also general or generic points. This concept is modeled very precisely by the concept of random points. It is through this modeling that probability enters into the numerical algorithms.

In §2 we discuss what we can compute using current numerical methods for finding the smooth or isolated solutions of a system of polynomials. We start by considering what we call square systems, i.e., systems with the same number of equations as variables. These are the systems that most of the software and papers on continuation treat. We then turn to systems with more equations than variables. We discuss a probabilistic way of reducing the system to a square system, which lets us compute much of the information contained in the original system.

The emphasis on the number of variables and equations is usual, but somewhat unnatural. Indeed a square system might have a number of equations that are functions of the preceding equations. In this case we will have no isolated solutions. Therefore we take a more intrinsic look in §3 at the algebraic set X defined by the equations. From this point of view what we want to find are a finite set of zeroes of an arbitrary system that includes at least one “general” point of each irreducible component of X. Since every isolated solution of the original system is an irreducible component of X we see that the solution set we are searching for includes the classical solution set found for square systems by continuation.

In §3.1 we give a probabilistic algorithm to do this based on the standard algorithms for solving N equations in N unknowns with probability 1. We use the classical technique of slicing by linear spaces. The slicing technique is a common tool in algebraic geometry, e.g., see [SS, BS].

In §3.2 we use the results of §2 and §3.1 to give a probabilistic algorithm for finding the dimension of an arbitrary algebraic set in ℂ^{N}.

In §3.3 we use the results of §2 and §3.1 to give a probabilistic algorithm for finding the maximum dimension of an irreducible component of a given algebraic set X ⊂ ℂ^{N} at a point p ∈ X.

In §3.4 we give an algorithm for deciding inclusion of reduced algebraic sets. In §3.5 we discuss briefly some of the problems from the theory of mechanisms

that motivate our numerical study of positive dimensional algebraic sets.

In §4 we discuss the methods and some of the questions that should be looked at next.

We would like to thank General Motors Research and Development Center for making this work possible. The first author would also like to thank the the University of Notre Dame.

1. A brief review of algebraic sets

First let us review the concept of irreducible components of an algebraic set. For simplicity we work only with affine algebraic sets, i.e., closed subsets of ℂ^{N}. This is sufficient for the purposes of this paper, and for numerical analysis problems there is a simple device we discuss in §4 to reduce general algebraic sets to affine algebraic sets. For a more detailed introduction to this material with references we refer to the appendix of [MS].

A reduced algebraic set is the subset of common zeroes in ℂ^{N} of a set of polynomials:

( f_{1}(x_{1}, …, x_{N}) ; ⋮ ; f_{n}(x_{1}, …, x_{N}) ) = 0. (1)

An algebraic set in ℂ^{N} is a reduced algebraic set with addition of a possibly nonreduced structure, a sheaf of germs of algebraic functions, that is derived from the defining equations f_{i}. This extra structure is related to multiplicities.

By a Zariski open subset U of a reduced algebraic set X we mean a set of the form U = X − A where A is a closed subset of X which is a reduced algebraic set also. The set of smooth (i.e., manifold) points of a reduced algebraic set X is denoted X_{reg}. X_{reg} is Zariski open in X and the closure of X_{reg} in the usual topology inherited from Euclidean space is X. We define the dimension of a reduced algebraic set X to be dim X_{reg}. Here we follow the convention that the dimension of a manifold is the maximum of the dimensions of the connected components of the manifold. We define the dimension of an arbitrary algebraic set to be the dimension of the underlying reduced algebraic set.

The next result is the analogue of the decomposition of a manifold into its connected components.

1. The set, X_{reg}, of smooth points of X is a union of finitely many connected components, X′_{i}, i = 1, …, r;

2. each X′_{i} is Zariski open in its closure, X_{i}, which is a reduced algebraic set.

3. X = ∪X_{i}. The X_{i} are called the irreducible components of X, and if X consists of one component then X is called irreducible.

In actual practice, it is important to work not just with reduced algebraic sets, but with algebraic sets and keep track of the extra structure. In this case, we still have the decomposition of the underlying reduced algebraic set, but in addition there will be multiplicities for the components.

To understand this extra structure let us give a few typical examples.

Let the system of polynomials be

z = 0

on ℂ. In this case, the algebraic set is the point 0 and the structure is the set of relevant functions on the point, i.e., the complex numbers ℂ. Usually we think of the algebraic set as the point 0, and forget about the functions on the point. We don’t lose anything by doing this.

As a second more complicated example, let the system of polynomials be

z^{2} = 0

on ℂ. In this case the algebraic set is the point 0 and the structure is the set of relevant “functions on the point,” i.e., ℂ[z]/(z^{2}), the Taylor series, c_{0} +c_{1}z, of order 1 at 0. We say the point 0 has multiplicity 2, and we think of the set as a fuzzy point. In this case we lose nothing by forgetting about the functions and thinking of the algebraic set as a point counted twice.

In higher dimensions things get a bit more complicated, and a single number like the multiplicity, doesn’t capture everything. A good example is the system of polynomials

( \begin{matrix}z^{2} & = & 0 \\ w & = & 0. \end{matrix} )

on ℂ^{2}. Here the algebraic set is the point (0, 0) with the set of “functions,” ℂ[z, w]/(z^{2}, w), i.e., c_{0} + c_{1}z. Again we have a point of multiplicity 2, but that single number doesn’t uniquely specify the algebraic set. Indeed consider the system of polynomials

( \begin{matrix}w^{2} & = & 0 \\ z & = & 0. \end{matrix} )

on ℂ^{2}. Here the algebraic set is still the point (0, 0) but this time with the set of “functions,” ℂ[z, w]/(w^{2}, z), i.e., c_{0} + c_{1}w. The different multiplicity 2 algebraic set structures on (0, 0) are parameterized by a ℙ^{1} (for much more on this type of example see [G]).

The irreducible decomposition with the multiplicity structure can be looked at as a generalization of the factorization of a polynomial p(z) on ℂ. Indeed if x_{1}, …, x_{r} are distinct points of ℂ with p(z) = cΠ_{i=1}^{r}(x − x_{i})^{n_{i}} for some c ≠ 0 and with n_{i} > 0 for all i, then the irreducible decomposition of the algebraic set associated with p(z) = 0 is the set x_{1}, …, x_{r} with multiplicities n_{1}, …, n_{r}. For a single polynomial on ℂ^{N} the corresponding picture is the similar. A polynomial p(z_{1}, …, z_{N}) cannot be factored into lower degree polynomials if and only if the algebraic set associated with p(z_{1}, …, z_{N}) = 0 is irreducible with multiplicity 1. In general every polynomial p(z_{1}, …, z_{N}) has a factorization into polynomials p(z_{1}, …, z_{N}) = Π_{i=1}^{r}p_{i}(z_{1}, …, z_{N})^{n_{i}} with n_{i} > 0 for all i, where for each i the polynomial p_{i} cannot be factored into polynomials of lower degrees (the factorization is unique up to multiplying the factors by nonzero constants). The irreducible components of the algebraic set associated to p(z_{1}, …, z_{N}) = 0 are the sets p_{i}(z_{1}, …, z_{N}) = 0. The extra structure is given by the multiplicities of the irreducible components which are just the numbers n_{i} in the factorization.

The irreducible components of the set of zeroes of a single nonconstant polynomial f(z_{1}, …, z_{n}) each have dimension N − 1. In general, if we have a system of the form (1) giving rise to an algebraic set X, we can only say that if X is not empty, then dim X ≥ N − n.

When we talk about points on an algebraic set X ⊂ ℂ^{N} being near and don’t specify otherwise, we are using the complex topology, i.e., the topology induced on the closed set X from the usual topology on ℂ^{N} induced by the Euclidean metric.

The notion of a general or a generic point is basic but also a bit slippery. The desire is for the concept to capture the idea of “a point of an algebraic set which has no special properties not possessed by the whole algebraic set.” For example, we would like to say that if a function vanishes at such a point, it vanishes everywhere on the algebraic set. As stated, this is nonsensical, but there are a number of ways to make the concept precise.

First it is a concept about irreducible algebraic sets. For an arbitrary algebraic set we talk only about generic points of the irreducible components.

In modern scheme theory the problem is finessed by defining the generic point to be a nonclosed point corresponding to the whole irreducible algebraic set. Such a point has no coordinates and doesn’t directly translate into anything numerical about the variety.

Another approach which is closer to our needs is to say that a generic (respectively general) point on an irreducible algebraic set X is any point that is contained in the complement U in X of a proper algebraic subset of X (respectively a countable union of proper algebraic subsets of X). This makes good sense. Also X − U is very thin, i.e., real codimension 2 and of measure zero. Often this captures what we want, but it still doesn’t give us a good numerical interpretation of the point.

The third and most classical meaning of generic is very close to numerical analysis. Here every irreducible algebraic set X in ℂ^{N} is defined over a subfield of the complex numbers. Indeed there are finitely many equations defining the set and there are finitely many coefficients for the equations and adding those to the rationals generates a field much smaller than ℂ. We can further add all coefficients of the equations defining each of the finite number of algebraic objects entering the discussion at hand to get a subfield F ⊂ ℂ relevant to the given discussion. Any point on the set X ⊂ ℂ^{N} with at least one coordinate a point in ℂ not algebraic over F is a generic point. This is the classical approach taken to algebraic geometry by A. Weil in his very useful introductory text [W]. There is also a brief but helpful discussion of this in [Mu]. From this point of view, a generic or general point of X is modeled very nicely by a point in X with random coordinates. Moreover in numerical problems we often have equations with real coefficients and work with irreducible algebraic sets (such as ℂ^{N}) in which the real points are Zariski dense. In this case we can model a generic point by a point in X with random real coordinates. Since the points on a computer are only finite we have a possible problem, but it is a problem common to much of scientific computing. We discuss this further in §4.

We use the term generic point throughout this paper. An objection that might be raised to using generic points in numerical analysis is that we have no way of knowing that a point is a generic point. Using generic points puts us in the same “state of sin” as users of random numbers [Kn, quote of Von Neumann on page 1]. What we can do is design “random number generators” which output streams of numbers satisfying properties random numbers would satisfy. This is of course nontrivial but random numbers are too useful not to use [Kn]. Based on algebraic geometry our algorithms output “generic points” using random number generators. We considered using the term random points instead of generic points, but decided against it because generic point is a useful classical concept in algebraic geometry with a clear connotation to an algebraic geometer. Our concept of generic points is about as close to the algebraic geometry notion as you can get using a system of finite arithmetic.

2. Finding isolated solutions by continuation

In this section we discuss the system (1) of n polynomials on ℂ^{N}. We first deal with square systems, i.e., systems where n = N. These are the systems which the current software solves. Next we turn to overdetermined systems, i.e., n > N, and show how to reduce them to square systems—this is easy when working with probabilistic methods. In §3 we turn to the study of algebraic sets, which includes the study of underdetermined systems, i.e., systems with n < N.

2.1. Finding isolated and smooth solutions of square systems. In this section we only deal with “square systems,” i.e., systems of N polynomials in N variables:

f(x) = ( f_{1}(x_{1}, …, x_{N}) ; ⋮ ; f_{N}(x_{1}, …, x_{N}) ) = 0.

Assuming that the system is not too “large” the current state of the art lets us find a set solutions on ℂ^{N} which include all isolated and in particular smooth solutions on ℂ^{N}. We refer to [WMS1] for a survey of the continuation methods we use. Our basic reference for continuation methods is [MS]. One important point to note is that when the continuation method leads to a solution that is isolated but singular there are methods (“endgames”) based on Cauchy’s integral and fractional power series for accurately computing the solution [MSW1]. When we are talking about large, we are talking about a priori bounds on the size of the Bezout number. What is too large is a function of time, dependent on the given technology. Currently a Bezout number bound by total degree or multihomogeneous degree of a few hundred thousand is on the boundary of what can be done. See [WMS2] for the solution of a large problem by this technology. The references [AG, MSW2, VVC] give more up to date information on what can be done. Since the continuation algorithms we use parallelize naturally, it is reasonable to expect in the near future an increase in the maximum size of a system that can be handled by a few orders of magnitude.

2.2. Overdetermined systems. In this section we assume that the system (1) on ℂ^{N} satisfies n > N. There is a natural procedure for obtaining a square system from the above system, f(x) = 0. Given a matrix of complex numbers

Λ = ( λ_{1,1} ⋯ λ_{1,n} ; ⋮ ; λ_{N,1} ⋯ λ_{N,n} )

we can form a square system

Λ · f = ( λ_{1,1}f_{1} + ⋯ + λ_{1,n}f_{n} ; ⋮ ; λ_{N,1}f_{1} + ⋯ + λ_{N,n}f_{n} ) = 0.

This system is equivalent to the system

g · Λ · f = 0

where g is an invertible N × N matrix. Thus the Grassmannian Gr(N, n) of N dimensional vector subspaces of ℂ^{n} parameterizes a natural family of square systems derived from our original system. Systems of the form

( f_{1} + λ_{1,N+1}f_{N+1} + ⋯ + λ_{1,n}f_{n} ; ⋮ ; f_{N} + λ_{N,N+1}f_{N+1} + ⋯ + λ_{N,n}f_{n} ) = 0. (2)

with

( λ_{1,N+1} ⋯ λ_{1,n} ; ⋮ ; λ_{N,N+1} ⋯ λ_{N,n} )

a matrix of complex (respectively real) numbers form a Zariski open dense (respectively Zariski dense) set in the systems parameterized by Gr(N, n).

What is the relation between the solutions of the original system (1) and the square systems (2) derived from it?

The next result is a straightforward consequence of the results in [MS].

Proposition 2.2.1. There are nonempty Zariski open dense sets of parameters λ_{i,j} ∈ ℂ^{N(n−N)} or λ_{i,j} ∈ ℝ^{N(n−N)} such that every isolated (respectively nonsingular) solution of

f = ( f_{1}(x_{1}, …, x_{N}) ; ⋮ ; f_{n}(x_{1}, …, x_{N}) ) = 0

is an isolated (respectively nonsingular) solution of

g = ( f_{1}(x_{1}, …, x_{N}) + ∑_{j=N+1}^{n} λ_{1,j}f_{j}(x_{1}, …, x_{N}) ; ⋮ ; f_{N}(x_{1}, …, x_{N}) + ∑_{j=N+1}^{n} λ_{N,j}f_{j}(x_{1}, …, x_{N}) ) = 0.

Moreover for any k > 0, a variety V is an irreducible k dimensional component of g = 0, if and only if it is a k dimensional component of f = 0. If V occurs with multiplicity 1 as an irreducible k dimensional component of f = 0, then it occurs with multiplicity 1 as an irreducible k dimensional component of g = 0.

We can choose which of our equations we want to play the role of f_{1}, …, f_{N}. This choice is dictated by the form of the equations so as to achieve the lowest Bezout number.

The randomization introduced by the parameters λ_{i,j} is necessary. For example consider the system:

\begin{aligned}xy & = & 0 \\ x(x + y) & = & 0 \\ y(x + y) & = & 0. \end{aligned}

Any two of the 3 equations have a 1 dimensional solution set, but all three together have the origin (with multiplicity 3) as the solution set.

It is important to note that although an irreducible component of f = 0 is an irreducible component of the randomized square system, g = 0, it’s multiplicity as an irreducible component of g (if not 1) might be larger than as an irreducible component of f = 0. The following system, which is equivalent to the above system, illustrates this:

\begin{aligned}xy & = & 0 \\ x^{2} & = & 0 \\ y^{2} & = & 0. \end{aligned}

The origin is an isolated solution of multiplicity 3. The randomized square system is:

\begin{aligned}x(y + μ_{1}x) & = & 0 \\ y(x + μ_{2}y) & = & 0. \end{aligned}

It has the origin as an isolated solution of multiplicity 4.

3. Probabilistic algorithms about algebraic sets

3.1. An algorithm for finding generic points. Let X ⊂ ℂ^{N} be the algebraic subset defined by the system of polynomials (1), and let X = ∪_{i∈I}X_{i} be the decomposition of X into irreducible algebraic sets. We would like to give a prescription to choose at least one general point from each of the X_{i}. Noether’s normalization theorem combined with Bertini’s theorem shows that if Y is an irreducible k-dimensional closed algebraic subset of ℂ^{N}, then a general affine linear subspace ℂ^{m} of ℂ^{N} meets Y only if m + k ≥ N and in this case ℂ^{m} meets Y in a set of dimension m + k − N. Moreover, ℂ^{m} is transverse to Y in the sense that, if letting Sing(Z) denote the singular set of an algebraic set Z, we stratify Y using Y ⊇ Sing(Y) ⊇ Sing(Sing(Y)) …, then ℂ^{m} is transverse to Y − Sing(Y), to Sing(Y) − Sing(Sing(Y)),…

Assume we are looking for generic points of the A-dimensional irreducible components of X. Assume we have a generic affine linear ℂ^{N−A}. This will meet X in generic points of the A-dimensional irreducible components of X. To find one of these points it suffices by the above to find the isolated solutions of the system (1) restricted to a generic affine linear ℂ^{N−A}. We refer to this as “solving the system restricted to ℂ^{N−A}.” But what does this mean in numerical terms?

The answer is straightforward. A generic ℂ^{N−A} is defined by

λ_{1,1}x_{1} + ⋯ + λ_{1,N}x_{N} = λ_{1}
⋮
λ_{A,1}x_{1} + ⋯ + λ_{A,N}x_{N} = λ_{A}

As earlier we can assume that the defining equations are of the form we would have using Gaussian elimination

x_{1} + λ_{1,A+1}x_{A+1} + ⋯ + λ_{1,N}x_{N} = λ_{1} ; ⋮ ; x_{A} + λ_{A,A+1}x_{A+1} + ⋯ + λ_{A,N}x_{N} = λ_{A} (3)

with all the λ’s random numbers. By the above we must find the isolated solutions of:

( f_{1}(x_{1}, …, x_{N}) = 0 ; ⋮ ; f_{n}(x_{1}, …, x_{N}) = 0 ; x_{1} + λ_{1,A+1}x_{A+1} + ⋯ + λ_{1,N}x_{N} = λ_{1} ; ⋮ ; x_{A} + λ_{A,A+1}x_{A+1} + ⋯ + λ_{A,N}x_{N} = λ_{A} ) .

We can first use the linear equations to reduce the system to the the system of the equations f_{i} = 0 restricted to the ℂ^{N−A} defined by (3). This is a system of n equations in N − A unknowns. We can use Proposition 2.2.1 to pass to a square system.

We give pseudocode below for this algorithm. We work from irreducible components of X of highest dimension down to lowest dimension.

Algorithm to find a finite set of points on an algebraic set including at least one generic point of each irreducible component of the algebraic set

Given polynomials f_{1}, …, f_{n} on ℂ^{N} defining an algebraic set X.

1. If 0 < n < N, then use N −n random linear equations to reduce to a system of n polynomials in n unknowns. Set N := n.

2. If n > N > 0, then using a generic (μ_{j,i}) ∈ ℂ^{N(n−N)} replace the system with a new system (after possibly reordering to lower the Bezout number)

( f_{1}(x_{1}, …, x_{N}) + ∑_{i=N+1}^{n} μ_{1,i}f_{i}(x_{1}, …, x_{N}) ; ⋮ ; f_{N}(x_{1}, …, x_{N}) + ∑_{i=N+1}^{n} μ_{N,i}f_{i}(x_{1}, …, x_{N}) ) = 0.

of N equations in N unknowns. Rename these equations f_{1}, …, f_{N}.

3. If N = n > 0, then set A = N and do the following procedure. While A ≥ 0:

a) “Solve” ( f_{1}(x_{1}, …, x_{N}) ; ⋮ ; f_{N}(x_{1}, …, x_{N}) ) = 0 restricted to a general affine linear

ℂ^{N−A}.

b) Store those points from a) that give a solution to the original system.

c) Set A := A − 1.

3.2. An algorithm for the dimension of an algebraic set. Let X ⊂ ℂ^{N} be an d-dimensional algebraic subset of ℂ^{N}.

How do we find X’s dimension numerically? To do this we simply need a stripped down version of the algorithm of §3.1 to find generic points.

Geometrically we know that a generic ℂ^{N−A} meets X if and only if d ≥ A. Thus we can successively intersect X with ℂ^{N−A}’s with A = N, N − 1, …, 0 until a solution is found or A = 0 and no solution is found. Note that this means that

3.2.1. for the largest A with ℂ^{N−A} ∩ X ≠ ∅ for a generic ℂ^{N−A} we have a finite intersection.

This is important because our solution procedure finds all isolated solutions of a square system.

Numerically assume we have a system of the form (1) consisting of n polynomial equations in N unknowns.

If n < N and there are any solutions, then they will form a set of dimension at least N − n. Thus, cutting with a generic ℂ^{N−n} we will have solutions if and only if the original system had a solution. Thus we can assume without loss of generality that n ≥ N. Starting with A = N and decrementing down to A = 0 it follows from 3.2.1 that for the first A where solutions appear, the solutions are isolated. If n = N, we have solution procedures guaranteed to find at least the isolated solutions. Thus we can run these procedures to check if there is a solution.

Algorithm to find the dimension of an algebraic set

Given polynomials f_{1}, …, f_{n} on ℂ^{N} defining an algebraic set X.

1. Set d := 0

2. If 0 < n < N, then use N −n random linear equations to reduce to a system of n polynomial equations in n unknowns. Set d := N − n and then set N := n.

3. If n > N > 0, then using a generic (μ_{j,i}) ∈ ℂ^{N(n−N)} replace the system with a new system (after possibly reordering to lower the Bezout number)

( f_{1}(x_{1}, …, x_{N}) + ∑_{i=N+1}^{n} μ_{1,i}f_{i}(x_{1}, …, x_{N}) ; ⋮ ; f_{N}(x_{1}, …, x_{N}) + ∑_{i=N+1}^{n} μ_{N,i}f_{i}(x_{1}, …, x_{N}) ) = 0.

of N equations in N unknowns. Rename these equations f_{1}, …, f_{N}.

4. Set A = N and do the following procedure. While A ≥ 0:

a) Replacing A of the variables using random linear equations, “solve”

( f_{1}(x_{1}, …, x_{N}) ; ⋮ ; f_{N}(x_{1}, …, x_{N}) ) = 0.

b) If there is at least one of these points from a) that gives a solution to the original system, then stop: the algebraic set X has dimension d + A.

c) Set A := A − 1.

5. X is empty.

3.3. An algorithm for the dimension of an algebraic set at a point. Let X be an algebraic subset of ℂ^{N} defined by a system of polynomial equations f = 0. Let p ∈ X, i.e., p ∈ ℂ^{N} and f(p) = 0. In this section we give an algorithm to compute the dimension of X at p, i.e., if X = ∪_{i=1}^{r}X_{i} is the decomposition of X into irreducible components, then the dimension of X at p ∈ X is max_{{i|p∈X_{i}}} dim X_{i}. In particular:

1. if p is a generic point of an irreducible component X_{i} of X, then this algorithm computes dim X_{i};

2. this algorithm lets us decide whether a solution p of a system f = 0 is

isolated. To check if p is a smooth solution, a simple check of whether the Jacobian matrix of f at p is of rank N suffices, but to our knowledge the algorithm below gives the first numerical algorithm in the literature to check whether a possibly singular solution is isolated.

The algorithm proceeds as follows. If X_{i} is an irreducible component of X containing p, then affine ℂ^{N−dim X_{i}}’s near a generic affine ℂ^{N−dim X_{i}} containing p meets X_{i} in at least one point near p. Moreover if dim X_{i} is the maximum dimension of any irreducible component of X containing p, then for A > dim X_{i} affine ℂ^{N−A}’s near a generic affine ℂ^{N−A} containing p don’t meet X in any points near p. An affine ℂ^{A} containing p := (p_{1}, …, p_{N}) is specified by

x_{1} + ∑_{j=A+1}^{N} λ_{1,j}x_{j} = p_{1} + ∑_{j=A+1}^{N} λ_{1,j}p_{j} ; ⋮ ; x_{A} + ∑_{j=A+1}^{N} λ_{A,j}x_{j} = p_{A} + ∑_{j=A+1}^{N} λ_{A,j}p_{j}. (4)

To get a generic affine ℂ^{N−A} we can choose a matrix

( λ_{1,A+1} ⋯ λ_{1,N} ; ⋮ ; λ_{A,A+1} ⋯ λ_{A,N} )

with random complex entries (or real entries).

At the moment checking nearness and choosing nearby affine ℂ^{N−A}’s requires judgement based on the given problem. This is not satisfactory: see Problem 4.1.

Algorithm to find the dimension of an algebraic set at a point

Given polynomials f_{1}, …, f_{n} on ℂ^{N} defining an algebraic set X and a point p ∈ X. 1. Set d := 0 2. If 0 < n < N, then use N − n random linear equations vanishing on p to reduce to a system of n polynomials in n unknowns. Set d := N − n and then set N := n.

3. If n > N > 0, then using a generic (μ_{j,i}) ∈ ℂ^{N(n−N)} replace the system with a new system (after possibly reordering to lower the Bezout number)

( f_{1}(x_{1}, …, x_{N}) + ∑_{i=N+1}^{n} μ_{1,i}f_{i}(x_{1}, …, x_{N}) ; ⋮ ; f_{N}(x_{1}, …, x_{N}) + ∑_{i=N+1}^{n} μ_{N,i}f_{i}(x_{1}, …, x_{N}) ) = 0.

of N equations in N unknowns. Rename these equations f_{1}, …, f_{N}.

4. Set A = N and do the following procedure. While A ≥ 1:

a) Using a system of random linear equations near a system of random linear equations vanishing at p to replace A of the variables, “solve”

( f_{1}(x_{1}, …, x_{N}) ; ⋮ ; f_{N}(x_{1}, …, x_{N}) ) = 0.

b) If there is at least one of these points from a) that gives a solution to the original system and is “near” to p then stop: the dimension of the algebraic set X at p has dimension d + A.

c) Set A := A − 1.

5. the dimension of the algebraic set X at p has dimension d.

3.4. An algorithm for deciding inclusion and equality of reduced algebraic sets. Let X be an algebraic subset of ℂ^{N} defined by a system of polynomial equations f = 0. Let Y be a second algebraic subset of ℂ^{N} defined by a system of polynomial equations g = 0. Using the algorithm of §3.1

Algorithm for deciding if every solution of f = 0 is a solution of g = 0

1. Find a finite set F of points containing generic points of each irreducible component of the algebraic set defined by f = 0.

2. Every solution of f = 0 is a solution of g = 0 if and only if g(x) = 0 for each x ∈ F.

To check if f = 0 and g = 0 have the same set of solutions (not counting multiplicities) it suffices to use the above algorithm twice, first to decide if every solution of f = 0 is a solution of g = 0, and if this is true then to decide if every solution of g = 0 is a solution of f = 0.

We have not dealt with multiplicities in this algorithm. Thus this algorithm gives a way of deciding if the reduced algebraic set defined by f = 0 is an algebraic subset of the reduced algebraic set defined by g = 0.

3.5. Some illustrative applications. As motivation for the algorithms discussed above, we briefly outline some questions in kinematics where the methods could prove useful. First, consider the singularity problem for a serial-link robot arm consisting of a chain of n moving links connected by rotational joints. The first link in the chain is additionally connected to an immobile base link by a rotational joint. By rotating the joints, the final “hand” link is moved around in space. Let θ ∈ T^{n} (the n-dimensional torus) be the joint angles and let x ∈ SE(3) = R^{3} × SO(3) be the position and orientation of the hand. The forward kinematic function f : T^{n} → SE(3) can be readily written down in terms of trigonometric functions using standard kinematic formulations. It is of considerable importance in the control of such a robot to identify points in the joint space where the Jacobian matrix f_{θ} becomes rank deficient, because at such points the robot cannot produce endeffector velocities in arbitrary directions in the tangent space. In the common case that n = 6, f_{θ} ∈ R^{6×6} is square and the singularity condition is det f_{θ} = 0. This condition is trigonometric, but it can be converted to an algebraic condition with the change of variables cos θ_{i} = (1 − t_{i}^{2})/(1 + t_{i}^{2}), sin θ_{i} = 2t_{i}/(1 + t_{i}^{2}) (i = 1, …, 6). This is an example where the algebraic set under study is given by one equation in six unknowns.

To be more specific, we consider the special case of a two-link planar robot arm having link lengths r_{1}, r_{2} and joint angles θ_{1}, θ_{2}. The first joint is at the origin and the forward kinematic function gives the endpoint of the arm as

f(θ) = (r_{1} cos θ_{1} + r_{2} cos θ_{2}, r_{1} sin θ_{1} + r_{2} sin θ_{2}).

The Jacobian matrix is

f_{θ} = ( \begin{matrix}−r_{1} sin θ_{1} & −r_{2} sin θ_{2} \\ r_{1} cos θ_{1} & r_{2} cos θ_{2} \end{matrix} )

and the singularity condition is

(5) det f_{θ} = −r_{1}r_{2}(sin θ_{1} cos θ_{2} − cos θ_{1} sin θ_{2}) = 0.

This problem is easily treated by hand since by a trigonometric identity Eq.(5) becomes

sin(θ_{1} − θ_{2}) = 0,

which implies θ_{1} = θ_{2} or θ_{1} = θ_{2}+π. These describe the two irreducible components of the singularity set.

To illustrate our algorithms, we treat the example numerically as follows. First, convert Eq.(5) to a polynomial expression using the change of variables given above to get

−2r_{1}r_{2} \frac{t_{1}(1 − t_{2}^{2}) − (1 − t_{1}^{2})t_{2}}{(1 + t_{1}^{2})(1 + t_{2}^{2})} = 0. (6)

Now, we use algorithms 3.1 and 3.2 to find generic points on the algebraic set defined by this equation and to determine its dimension. We begin with n = 1 and N = 2. Step 1 is to add one random linear equation to make the system square:

.2338t_{1} + .8374t_{2} + .4721 = 0.

In our accounting we now set d = 1 and N = 1. The linear equation is used to eliminate one variable to obtain a cubic in the remaining variable. Next, we set A = 1 and test a random point in ℂ^{0}. Since the cubic is not trivially zero, there are no solutions and the singularity set is not 2-dimensional. Hence, we decrement A to 0 and solve the cubic equation to get three solutions (rounded to 12 decimal places):

(t_{1}, t_{2}) = ( \begin{matrix}−3.154625658213, & .316994822316 \\ −.440720672006, & −.440720672006 \\ 1.135378430734, & −.880763605271 \end{matrix} ) .

From this result, we conclude that the singularity set is of dimension d + A = 1 and that the three solution points are close approximations to at least one point on each irreducible component of the singularity set. An astute analyst might easily guess from the second solution that one component is t_{1} = t_{2} and upon converting the first and third solutions back to angles might also guess the second component as θ_{1} = θ_{2} + π. (The latter case corresponds to the component t_{1}t_{2} = −1.) Algorithm 3.4 can be used to test this hypothesis. Repeating the solution for a different random linear equation (to ensure independence from our hypothesis) we will find that the solutions do indeed fall on the hypothesized sets. This shows that the singularity set is included in the union of t_{1} = t_{2} and t_{1}t_{2} = −1. To test the opposite inclusion, we check random points from these two sets against Eq.(6). It will be zero to high precision, thus indicating equality of the singularity set to the two component sets. As always, these numerical conclusions do not carry the same level of certitude as the analytical proof, but at a minimum, one may say that they provide strong insight to the truth.

Another question in kinematics concerns the study of overconstrained mechanisms. Consider again a serial chain of links, but this time additionally connect the final link rigidly to the base link. The closure equation for the closed-loop linkage is f(θ, p) = id, where id is the identity element in SE(3), and where p represents the geometric parameters of the links, such as link lengths, twists, and offsets. All single-loop mechanisms with n ≥ 7 joints will have internal motion, i.e., a positive dimensional solution set in θ with p held fixed. For n ≤ 6, the closure equation gives at least as many algebraic equations as there are joint variables, so the device may or may not have internal motion depending on whether the link parameters p are especially fortuitous. If the device has a 1-dimensional joint motion, it is called an overconstrained linkage. It turns out that for n ≤ 3 there are no overconstrained linkages of any interest. Delassus [D] showed that the three known families of overconstrained linkages for n = 4 are exhaustive, and the classification has been generalized to include other types of joints by Waldron [W2]. For n = 5, 6, many examples have been found, but it is not known yet if the classification is complete [W1, MR1, MR2]. In such studies, it is useful to be able to verify a candidate mechanism given by a set of parameters, say p^{∗}. If the dimension of f(θ, p^{∗}) = id as computed by Algorithm 3.2 is 1, then the device is a valid overconstrained linkage. If the dimension is greater than 1, the device might still have a 1-dimensional irreducible component. This can be checked by applying Algorithm 3.3 to each of the generic points found by Algorithm 3.1.

4. Closing remarks and open problems

One question is how to deal with general algebraic sets (in this paper we dealt with affine sets, i.e., algebraic sets whose underlying set of points embeds as a closed subset of ℂ^{N}). This is straightforward. For example, here is a simple approach that covers most questions that come up in practice. Let X be an arbitrary projective algebraic set, i.e., a not necessarily reduced or irreducible algebraic subset X of complex projective space ℙ^{N}. If we choose any linear hyperplane H of ℙ^{N}, then ℙ^{N} − H, the complement of H in ℙ^{N}, is identified with ℂ^{N}. If H is chosen to be general, then all isolated points of X belong to ℙ^{N} − H, and in fact no irreducible component of X lies in H. Thus every irreducible component of H ∩ X has real codimension 2 in some irreducible component of X. Thus the affine algebraic set X − H ∩ X is—for numerical analysis questions about isolated points or irreducible components—almost identical to X. This is in fact a standard and useful device, e.g., see Morgan’s projective transformation [Mo]. A variant of this procedure works for a quasi-projective algebraic set X, i.e., X is a projective algebraic set Y minus a closed algebraic subset Z. There is some degree d such that the ideal sheaf of Z is generated by degree d homogeneous polynomials. Choose a general polynomial p of that degree which vanishes on Z and note that Y − p^{−1}(0) is an algebraic set which can be embedded as a closed subset of a Euclidean space, Y − p^{−1}(0) ⊂ X, and Y − p^{−1}(0) contains all isolated points of X. Moreover no irreducible component of X lies in p^{−1}(0). Thus replacing X by Y − p^{−1}(0) reduces us for most numerical analysis purposes to the types of sets we have been considering.

Problem 4.1. How do we decide in our algorithms that some quantity is “zero,” or that some point is “near” another point.

This problem, which is analogous to deciding when a matrix is “singular,” has no easy answer. Large problems often require individual inspection of unclear cases based on the special properties of a given problem. This is not an acceptable solution of the problem.

Problem 4.2. How secure are we using probabilistic algorithms of the sort we describe?

Our experience with solving systems of polynomials using probabilistic algorithms has been very good. Nevertheless a generic point might be close enough to a point that we want to avoid that numerical analysis difficulties ensue. This points out the need for some more quantitative measures of the size of numerically bad sets. Though such bounds will likely not be so useful in practice, they will be useful theoretically.

Problem 4.3. Enumerate the irreducible components of an algebraic set defined by the system (1) of n polynomials on ℂ^{N}.

We know by Algorithm 3.1 how to find generic points of irreducible components, and by Algorithm 3.3 how to find the dimension of the irreducible component that a generic point belongs to. What we don’t know how to do is to decide when two generic points are on the same irreducible component.

Problem 4.4. Let p ∈ X be a generic point of an irreducible component Z of the algebraic set X defined by the polynomial system (1) on ℂ^{N}. Give a numerical algorithm to compute the multiplicity of Z in X.

It suffices to do this for isolated points. To see this note that Algorithm 3.3 finds the dimension A of Z. Restricting f = 0 to a generic ℂ^{N−A} containing p we obtain a system with p an isolated solution. The multiplicity of p for this new system is equal to the multiplicity of Z for the original system. For square systems the computation of the multiplicity is a byproduct of finding all isolated solutions. The problem 4.4 thus reduces to the following problem.

Problem 4.5. Let p ∈ ℂ^{N} be an isolated solution of the system (1) of n polynomials on ℂ^{N} with n > N. Give a probabilistic numerical algorithm to compute the multiplicity of p.

Using continuation and an associated randomized square system computes an upper bound for the multiplicity, but as the example system at the end of §2.2 shows this upper bound can be strictly larger.

References

[AG] E. Allgower and K. Georg, Numerical Continuation Methods: An Introduction, Computational Mathematics, 13 (1990), Springer-Verlag, Berlin.

[BS] M. Beltrametti and A. Sommese, The Adjunction Theory of Complex Projective Varieties, Expositions in Mathematics, 16 (1995), Walter De Gruyter, Berlin.

[D] E. Delassus, Les chaines articuléres fermées et déformables ã quatre membres, Bull. Sci. Math. et Astron. 46 (1922), 283–304.

[G] L. Göttsche, Hilbert Schemes of Zero-dimensional Subschemes of Smooth Varieties, Lecture Notes in Mathematics, 1572 (1994), Springer-Verlag, New York.

[Kn] D.E. Knuth, The Art of Computer Programming. Volume 2: Seminumerical Algorithms, second edition, (1981) Addison-Wesley, Reading, Mass.

[Mo] A. Morgan, Solving Polynomial Systems Using Continuation for Scientific and Engineering Problems, (1987) Prentice-Hall, Englewood Cliffs, N.J.

[MS] A. Morgan and A. Sommese, Coefficient-parameter polynomial homotopies, Appl. Math. Comput. 29 (1989), 123–160; Erratum, 51 (1992), p. 207.

[MSW1] A. Morgan, A. Sommese, and C. Wampler, A power series method for computing singular solutions to nonlinear analytic systems, Numerische Math. 63 (1992), 391–409.

[MSW2] , A generic product-decomposition bound for Bezout numbers, SIAM J. Numer. Anal. 32 (1995), 1308–1325.

[Mu] D. Mumford, Algebraic Geometry I, Grundlehren Math. Wiss. 221 (1976), Springer-Verlag, New York.

[MR1] C. Mavroidis and B. Roth, Analysis of overconstrained mechanisms, ASME J. Mechanical Design 117 (1995), 69–74.

[MR2] , New and revised overconstrained mechanisms, ASME J. Mechanical Design 117 (1995), 75–82.

[SS] B. Shiffman and A. Sommese, Vanishing Theorems on Complex Manifolds, Progress in Mathematics, 56 (1985), Birkhäuser, Boston.

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

[W1] K. Waldron, A study of overconstrained linkage geometry by solution of the closure equations—Part I. Method of study, Mechanism and Machine Theory, 8 (1973), 95–104.

[W2] , A study of overconstrained linkage geometry by solution of the closure equations—Part II. Four-bar linkages with lower pair joints other than screw joints, Mechanism and Machine Theory 8 (1973), 233–247.

[WMS1] C. Wampler, A. Morgan, and A. Sommese, Numerical continuation methods for solving polynomial systems arising in kinematics, ASME Journal of Mechanical Design 112 (1990), 59–68.

[WMS2] , Complete solution of the nine-point path synthesis problem for four-bar linkages, ASME Journal of Mechanical Design 114 (1992), 153–159.

[W] A. Weil, Foundations of Algebraic Geometry, revised and enlarged edition (1962), Colloquium Publications, 29, American Mathematical Society, Providence, Rhode Island.

Department of Mathematics, University of Notre Dame, Notre Dame, Indiana 46556, U.S.A.

    E-mail address: sommese.1@nd.edu

Mathematics Department, General Motors Research and Development Center, Bldg. 1-6, 30500 Mound Rd., Warren, Michigan 48090-9055, U.S.A.

    E-mail address: cwampler@math.ma.gmr.com