Real Polynomial Systems

Introduction

A real polynomial system is an expression constructed with polynomial equations and inequalities

combined using logical connectives and quantifiers

An occurrence of a variable x inside is called a bound occurrence; any other occurrence of x is called a free occurrence. A variable x is called a free variable of a real polynomial system if the system contains a free occurrence of x. A real polynomial system is quantifier free if it contains no quantifiers.

An example of a real polynomial system with free variables x, y, and z is the following

Any real polynomial system can be transformed to the prenex normal form

where each is or , and is a quantifier-free formula called the quantifier-free part of the system.

Any quantifier-free real polynomial system can be transformed to the disjunctive normal form

where each is a polynomial equation or inequality.

Reduce, Resolve, and FindInstance always put real polynomial systems in the prenex normal form, with quantifier-free parts in the disjunctive normal form, and subtract sides of equations and inequalities to put them in the form

In the following sections we will always assume the system has been transformed to the previous form.

Reduce can solve arbitrary real polynomial systems. For a system with free variables , the solution (possibly after expanding with respect to ) is a disjunction of terms of the form

where is one of

and and are algebraic functions (expressed using Root objects or radicals) such that for all satisfying , and are well defined (that is, denominators and leading terms of Root objects are nonzero), real valued, continuous, and satisfy inequality .

The subset of described by formula (4) is called a cell. The cells described by different terms of solution of a real polynomial system are disjoint.

This solves the system (1). The cells are represented in a nested form.

In[1]:=

This defines a function expanding with respect to .

In[2]:=

Here is the solution of the system (1) written explicitly as a union of disjoint cells.

In[5]:=

Resolve can eliminate quantifiers from arbitrary real polynomial systems. If no variables are specified in the input and all input polynomials are at most linear in the bound variables, Resolve may be able to eliminate the quantifiers without solving the resulting system. Otherwise, Resolve uses the same algorithm and gives the same answer as Reduce.

This eliminates quantifiers from the system (1).

In[6]:=

FindInstance can handle arbitrary real polynomial systems, giving instances of real solutions or an empty list for systems that have no solutions. If the number of instances requested is more than one, the instances are randomly generated from the full solution of the system and therefore may depend on the value of the RandomSeed option. If one instance is requested and the system does not contain general () quantifiers, a faster algorithm producing one instance is used and the instance returned is always the same.

This finds a solution of the system (1).

In[7]:=

The main general tool used in solving real polynomial systems is the Cylindrical Algebraic Decomposition (CAD) algorithm (see, for example, [1]). CAD for quantifier-free systems is available in Mathematica directly as CylindricalDecomposition. There are also several other algorithms used to solve special case problems.

Cylindrical Algebraic Decomposition

Semi-algebraic Sets and Cell Decomposition

A subset of is semi-algebraic if it is a solution set of a quantifier-free real polynomial system. According to Tarski's theorem [2], solution sets of arbitrary (quantified) real polynomial systems are semi-algebraic.

Every semi-algebraic set can be represented as a finite union of disjoint cells [3] defined recursively as follows.

A cell in is a point or an open interval.

A cell in has one of the two forms

where is a cell in , r is a continuous algebraic function, and are continuous algebraic functions, -, or , and on .

By an algebraic function we mean a function for which there is a polynomial

such that

In Mathematica algebraic functions can be represented as Root objects or radicals.

The CAD algorithm, introduced by Collins [4], computes a cell decomposition of solution sets of arbitrary real polynomial systems. The objective of the original Collins algorithm was to eliminate quantifiers from a quantified real polynomial system and to produce an equivalent quantifier-free polynomial system. After finding a cell decomposition, the algorithm performed an additional step of finding an implicit representation of the semi-algebraic set in terms of polynomial equations and inequalities in the free variables. The objective of Reduce is somewhat different. Given a semi-algebraic set presented by a real polynomial system, quantified or not, Reduce finds a cell decomposition of the set, explicitly written in terms of algebraic functions.

While Reduce may use other methods to solve the system, CylindricalDecomposition gives a direct access to the CAD algorithm. For a quantifier-free real polynomial system, CylindricalDecomposition gives a nested formula representing disjunction of cells in the solved form (4). As in the output of Reduce, the cells are disjoint and additionally are always ordered lexicographically with respect to ranges of the subsequent variables.

This finds a cell decomposition of an annulus.

In[8]:=

The Projection Phase of the CAD Algorithm

Finding a cell decomposition of a semi-algebraic set using the CAD algorithm consists of two phases, projection and lifting. In the projection phase, we start with the set of factors of the polynomials present in the quantifier-free part of the system (2) and eliminate variables one by one using a projection operator such that

Generally speaking, if all polynomials of have constant signs on a cell , then all polynomials of are delineable over , that is, each has a fixed number of real roots on as a polynomial in , the roots are continuous functions on , they have constant multiplicities, and two roots of two of the polynomials are equal either everywhere or nowhere in . Variables are ordered so that

This way the roots of polynomials of are the algebraic functions needed in the construction of the cell decomposition of the semi-algebraic set.

Several improvements have reduced the size of the original Collins projection. The currently best projection operator applicable in all cases is due to Hong [5]; however, in most situations we can use a smaller projection operator given by McCallum [6, 7], with an improvement by Brown [8]. There are even smaller projection operators that can be applied in some special cases. When equational constraints are present, we can use the projection operator suggested by Collins [9], and developed and proven by McCallum [10, 11]. When there are no equations and only strict inequalities, and there are no free variables or we are interested only in the full-dimensional part of the semi-algebraic set, we can use an even smaller projection operator described in [12, 13]. For systems containing equational constraints that generate a zero-dimensional ideal, Gröbner bases are used to find projection polynomials.

Mathematica uses the smallest of the previously mentioned projections that is appropriate for the given example. Whenever applicable, we use the equational constraints; otherwise, we attempt to use McCallum's projection with Brown's improvement. When the system does not turn out to be well oriented, we compute Hong's projection.

The Lifting Phase of the CAD Algorithm

In the lifting phase, we find a cell decomposition of the semi-algebraic set. Generally speaking, although the actual details depend on the projection operator used, we start with cells in consisting of all distinct roots of and the open intervals between the roots. We find a sample point in each of the cells and remove the cells whose sample points do not satisfy the system describing the semi-algebraic set (the system may contain conditions involving only ). Next we lift the cells to cells in , one dimension at a time. Suppose we have lifted the cells to . To lift a cell to , we find the real roots of with replaced with the coordinates of the sample point in . Since the polynomials of are delineable on , each root is a value of a continuous algebraic function at , and the function can be represented as a th root of a polynomial such that is the th root of . Now the lifting of the cell to will consist of graphs of these algebraic functions and of the slices of between the subsequent graphs. The sample points in each of the new cells will be obtained by adding the st coordinate to , equal to one of the roots, or to a number between two subsequent roots. As in the first step, we remove those lifted cells whose sample points do not satisfy the system describing the semi-algebraic set.

If , is a quantifier variable and we may not need to construct all the lifted cells. All we need is to find the (necessarily constant) truth value of on . If , we know that the value is True as soon as the truth value of on one of the lifted cells is True. If , we know that the value is False as soon as the truth value of on one of the lifted cells is False.

The coefficients of sample points computed this way are in general algebraic numbers. To save costly algebraic number computations, Mathematica uses arbitrary-precision floating point number (Mathematica "bignum") approximations of the coefficients, whenever the results can be validated. Note that using approximate arithmetic may be enough to prove that two roots of a polynomial or a pair of polynomials are distinct, and to find a nonzero sign of a polynomial at a sample point. What we cannot prove with approximate arithmetic is that two roots of a polynomial or a pair of polynomials are equal, or that a polynomial is zero at a sample point. However, we can often use information about the origins of the cell to resolve these problems. For instance, if we know that the resultant of two polynomials vanishes on the cell, and these two polynomials have exactly one pair of complex roots that can be equal within the precision bounds, we can conclude that these roots are equal. Similarly, if the last coordinate of a sample point was a root of a factor of the given polynomial, we know that this polynomial is zero at the sample point. If we cannot resolve all the uncertainties using the collected information about the cell, we compute the exact algebraic number values of the coordinates. For more details, see [14].

Decision Problems, FindInstance, and Assumptions

A decision problem is a system with all variables existentially quantified, that is, a system of the form

where are all variables in . Solving a decision problem means deciding whether it is equivalent to True or to False, that is, deciding whether the quantifier-free system of polynomial equations and inequalities has solutions.

All algorithms used by Mathematica to solve real polynomial decision problems are capable of producing a point satisfying if the system has solutions. Therefore the algorithms discussed in this section are used not only in Reduce and Resolve for decision problems, but also in FindInstance, whenever a single instance is requested and the system is quantifier free or contains only existential quantifiers. The algorithms discussed here are also used for inference testing by Mathematica functions using assumptions such as Simplify, Refine, Integrate, and so forth.

Solving this decision problem proves that the set contains the disk of radius 4/5 centered at the origin.

In[9]:=

This shows that does not contain the unit disk and provides a counterexample: a point in the unit disk that does not belong to .

In[10]:=

The primary method that allows Mathematica to solve arbitrary real polynomial decision problems is the CAD algorithm. There are, however, several other special case algorithms that provide much better performance in cases in which they are applicable.

When all polynomials are linear with rational number or floating point number coefficients, Mathematica uses a method based on the Simplex linear programming method. For other linear systems, Mathematica uses a variant of the Loos-Weispfenning linear quantifier elimination algorithm [15]. When the system contains no equations and only strict inequalities, a faster "generic" version of CAD is used [12, 13]. For systems containing equational constraints that generate a zero-dimensional ideal, Mathematica uses Gröbner bases to find a solution. For nonlinear systems with floating point number coefficients, an inexact coefficient version of CAD [16] is used.

There are also some special case methods that can be used as preprocessors to other decision methods. When the system contains an equational constraint linear with a constant coefficient in one of the variables, the constraint is used to eliminate the linear variable. If there is a variable that appears in the system only linearly with constant coefficients, the variable is eliminated using the Loos-Weispfenning linear quantifier elimination algorithm [15].

There are a two other special cases of real decision algorithms available in Mathematica. An algorithm by Aubry, Rouillier, and Safey El Din [17] applies to systems containing only equations. There are examples for which the algorithm performs much better than CAD; however, for randomly chosen systems of equations, it seems to perform significantly worse; therefore, it is not used by default. Setting the system option ARSDecision in the InequalitySolvingOptions group to True causes Mathematica to use the algorithm. Another algorithm by G.X. Zeng and X.N. Zeng [18] applies to systems that consist of a single strict inequality. Again, the algorithm is faster than CAD for some examples, but slower in general; therefore, it is not used by default. Setting the system option ZengDecision in the InequalitySolvingOptions group to True causes Mathematica to use the algorithm.

Arbitrary Real Polynomial Systems

Solving Real Polynomial Systems

According to Tarski's theorem [2], the solution set of an arbitrary (quantified) real polynomial system is a semi-algebraic set. Reduce gives a description of this set in the solved form (4).

This shows for what the set contains the disk of radius centered at the origin.

In[11]:=

This gives the projection of on the plane along the -axis.

In[12]:=

This finds the projection of Whitney's umbrella on the plane along the -axis.

In[13]:=

Here we find the interior of the previous projection set by directly using the definition.

In[14]:=

Quantifier Elimination

The objective of Resolve with no variables specified is to eliminate quantifiers and produce an equivalent quantifier-free formula. The formula may or may not be in a solved form, depending on the algorithm used.

Producing a fully solved quantifier-free formula here is difficult because of the complexity of polynomials in , , and appearing in the input. However, since appears in the input polynomials only linearly, the quantifier can be quickly eliminated using the Loos-Weispfenning linear quantifier elimination algorithm, which depends very little on the complexity of coefficients.

In[15]:=

Algorithms

The primary method used by Mathematica for solving real polynomial systems and real quantifier elimination is the CAD algorithm. There are, however, simpler methods applicable in special cases.

If the system contains an equational constraint in a variable from the innermost quantifier, the constraint is used to simplify the system using the identity

Note that if or is a nonzero constant, this eliminates the variable .

If all polynomials in the system are linear in a variable from the innermost quantifier, the variable is eliminated using the Loos-Weispfenning linear quantifier elimination algorithm [15].

The CAD algorithm is used when the previous two special case methods are no longer applicable, but there are still quantifiers left to eliminate or a solution is required.

For systems containing equational constraints that generate a zero-dimensional ideal, Mathematica uses Gröbner bases to find the solution set.

Options

The Mathematica functions for solving real polynomial systems have a number of options that control the way that they operate. This section gives a summary of these options.

Cubics and Quartics

By default, Reduce does not use the Cardano formulas for solving cubics or quartics over the reals.

In[16]:=

Setting options Cubics and Quartics to True makes Reduce use the Cardano formulas to represent numeric solutions of cubics and quartics.

In[17]:=

Solutions of cubics and quartics involving parameters will still be represented using Root objects.

In[18]:=

This is because the Cardano formulas do not separate real solutions from nonreal ones. For instance, in this case, for the third radical solution is real, but for the first radical solution is real.

In[19]:=

In[20]:=

WorkingPrecision

The setting of WorkingPrecision affects the lifting phase of the CAD algorithm. With a finite working precision prec, sample points in the first variable lifted are represented as arbitrary-precision floating point numbers with prec digits of precision. When we compute sample points for subsequent variables, we find roots of polynomials whose coefficients depend on already computed sample point coordinates and therefore may be inexact. Hence coordinates of sample points will have precision prec or lower. Determining sign of polynomials at sample points is simply done by evaluating Sign of the floating point number obtained after the substitution. Using a finite WorkingPrecision may allow getting the answer faster; however, the answer may be incorrect or the computation may fail due to loss of precision.

This problem is too hard for Reduce working in infinite WorkingPrecision, due to high degrees of algebraic numbers involved. Using 30 digits of precision sample points solves it in under two seconds.

In[21]:=

ReduceOptions Group of System Options

Here are the system options from the ReduceOptions group that may affect the behavior of Reduce, Resolve, and FindInstance for real polynomial systems. The options can be set with

FactorInequalities

Using transformations

at the input preprocessing stage may speed up the computations in some cases. In general, however, it does not make the problem easier to solve, and, in some cases, it may make the problem significantly harder. By default, these transformations are not used.

Here Reduce does not use transformations (7).

In[22]:=

Using transformations (7) speeds up the first example; however, it makes the other two examples significantly slower. The second example suffers from exponential growth of the number of inequalities. By replacing with in the third example, we get a degree 21 system in instead of a degree 3 system in .

In[26]:=

In[31]:=

ReorderVariables

By default, Reduce is not allowed to reorder the specified variables. Variables appearing earlier in the variable list may be used to express solutions for variables appearing later in the variable list, but not vice versa.

In[32]:=

Setting the system option ReorderVariables->True allows Reduce to pick a variable order that makes the system easier to solve.

In[33]:=

In[35]:=

InequalitySolvingOptions Group of System Options

Here are the system options from the InequalitySolvingOptions group that may affect the behavior of Reduce, Resolve, and FindInstance for real polynomial systems. The options can be set with

ARSDecision

The option ARSDecision specifies whether Mathematica should use the algorithm by Aubry, Rouillier, and Safey El Din [17]. The algorithm applies to decision problems containing only equations. There are examples for which the algorithm performs much better than the CAD algorithm; however, for randomly chosen systems of equations it seems to perform significantly worse. Therefore it is not used by default. Here is a decision problem (referred to as butcher8 in the literature), which is not done by CAD in 1000 seconds, but which can be done quite fast by the algorithm given in [17].

In[36]:=

In[37]:=

BrownProjection

By default, the Mathematica implementation of the CAD algorithm uses Brown's improved projection operator [8]. The improvement usually speeds up computations substantially. There are some cases where using Brown's projection operator results in a slight slowdown. The option BrownProjection specifies whether Brown's improvement should be used. In the first example [21], using Brown's improved projection operator results in a speedup by a factor of 3; in the second, it results in a 20% slowdown.

In[38]:=

In[43]:=

In[44]:=

In[47]:=

CAD

The option CAD specifies whether Mathematica is allowed to use the CAD algorithm. With CAD set to False, computations that require CAD will fail immediately instead of attempting the high complexity CAD computation. With CAD enabled, this computation is not done in 1000 seconds.

In[48]:=

In[50]:=

CADDefaultPrecision

By default, Mathematica uses validated numeric computations in the lifting phase of the CAD algorithm, reverting to exact algebraic number computations only if the numeric computations cannot be validated [14]. The option CADDefaultPrecision specifies the initial precision with which the sample point coordinates are computed. Choosing the value of CADDefaultPrecision is a tradeoff between speed of numeric computations and the number of points where the algorithm reverts to exact computations due to precision loss. With the default value of 100 bits, the cases where the algorithm needs to revert to exact computations due to precision loss seem quite rare. Setting CADDefaultPrecision to Infinity causes Mathematica to use exact algebraic number computations in the lifting phase of CAD. Here is an example that runs fastest with the lowest CADDefaultPrecision setting. (Specifying values lower than 16.2556 (54 bits) results in CADDefaultPrecision being set to 16.2556.) With CADDefaultPrecision->Infinity, the example did not finish in 1000 seconds.

In[51]:=

In[52]:=

In[54]:=

CADSortVariables

The performance of the CAD algorithm often depends quite strongly on the order of variables used. Some aspects of the variable ordering are fixed by the problem we are solving: quantifier variables need to be projected before free variables, and variables from innermost quantifiers need to be projected first. Variables specified in Reduce and Resolve cannot be reordered unless ReorderVariables is set to True. This, however, still leaves some freedom in ordering of variables: variables from the same quantifier can be reordered, and so can be variables given to FindInstance. By default, Mathematica uses a variable ordering heuristic to determine the order of these variables. In most cases the heuristic improves the performance of CAD; in some examples, however, the heuristic does not pick the best ordering. Setting CADSortVariables to False disables the heuristic and the order of variables used is as given in the quantifier variable list or in the variable list argument to FindInstance. Here is an example [21] that without reordering of quantified variables does not finish in 1000 seconds.

In[55]:=

This shows the optimal variable ordering for the example.

In[56]:=

In[57]:=

In[58]:=

CADZeroTest

One of the most time-consuming operations in the lifting phase of the CAD algorithm is determining the sign of a polynomial evaluated at a sample point with algebraic number coordinates. We try to avoid the problem by using sample points with arbitrary-precision floating point number coordinates and keeping track of the "genealogy" of projection polynomials and sample points in order to validate the results. However, if some of the results cannot be validated, we have to revert to computations with exact algebraic number coordinates. To determine the sign of a polynomial evaluated at a sample point with algebraic number coordinates, we first evaluate the polynomial at numeric approximations of the algebraic numbers. If the result is nonzero (that is zero is not within the error bounds of the resulting bignum), we know the sign. Otherwise, we need to test whether a polynomial expression in algebraic numbers is zero. The value of the CADZeroTest option specifies what zero testing method should be used at this moment. The value should be a pair {t, acc}. With the default value Mathematica computes an accuracy eacc such that if the expression is zero up to this accuracy, it must be zero. If eacc≤acc, the value of the expression is computed up to accuracy eacc and its sign is checked. Otherwise, the expression is represented as a single Root object using RootReduce and the sign of the Root object is found. With the default value we revert to RootReduce if eacc>$MaxPrecision. If , RootReduce is always used. If , expressions that are zero up to accuracy acc are considered zero. This is the fastest method, but, unlike the other two, it may give incorrect results because expressions that are nonzero but close to zero may be treated as zero.

This example runs faster with the CAD algorithm using the 30 digits of accuracy numeric zero test. The result in this example is correct; however, this setting of CADZeroTest may lead to incorrect results.

In[59]:=

In[63]:=

ContinuedFractionRootIsolation

To isolate real roots of polynomials, Mathematica uses methods based on Descartes' rule of sign. There are two interval subdivision strategies implemented, one based on interval bisection and another based on continued fractions (see [19] for details). The variant based on continued fractions is generally faster and is used by default. Setting ContinuedFractionRootIsolation to False causes Mathematica to use the interval bisection variant.

Here is an example where the speed difference between the two root isolation methods affects Reduce timing. We need to clear the Root cache between the Reduce calls; otherwise, the second call would save time on factoring the 400th degree polynomial when Root objects are created.

In[64]:=

In[71]:=

FGLMBasisConversion

For systems with equational constraints generating a zero-dimensional ideal I, Mathematica uses a variant of the CAD algorithm that finds projection polynomials using Gröbner basis methods. If the lexicographic order Gröbner basis of I does not contain linear polynomials with constant coefficients in every variable but the last one, then for every variable we find a univariate polynomial in that belongs to I. Mathematica can do this in two ways. By default, it uses a method based on GroebnerWalk computations. Setting FGLMBasisConversion to True causes Mathematica to use a method based on [20].

The method based on [20] seems to be slightly slower in general.

In[72]:=

In[76]:=

FGLMElimination

The FGLMElimination option specifies whether Mathematica should use a special case heuristic applicable to systems with equational constraints generating a zero-dimensional ideal I. The heuristic uses a method based on [20] to find in I polynomials that are linear (with a constant coefficient) in one of the quantified variables and uses such polynomials for elimination. The method can be used both in the decision algorithm and in quantifier elimination. With the default Automatic setting, it is used only in Resolve with no "solve" variables specified and for systems with at least two free variables.

This by default uses the elimination method based on [20], and returns a quantifier-free system in an unsolved form.

In[77]:=

With FGLMElimination set to False, the example takes longer to compute and the answer is in a solved form. (We show N of the answer for better readability.)

In[78]:=

If there is only one free variable, Resolve by default does not use the elimination method based on [20]. (We show N of the answer for better readability.)

In[80]:=

With FGLMElimination set to True, the example takes longer to compute and the answer is given in an unsolved form.

In[82]:=

In[84]:=

GenericCAD

Mathematica uses a simplified version of the CAD algorithm described in [13] to solve decision problems or find solutions of real polynomial systems that do not contain equations. The method finds a solution or proves that there are no solutions if all inequalities in the system are strict (< or >). The method is also used for systems containing weak ( or ) inequalities. In this case, if it finds a solution of the strict inequality version of the system, it is also a solution of the original system. However, if it proves that the strict inequality version of the system has no solutions, the full version of the CAD algorithm is needed to decide whether the original system has solutions. The system option GenericCAD specifies whether Mathematica should use the method.

Here the GenericCAD method finds a solution of the strict inequality version of the system.

In[85]:=

Without GenericCAD, finding a solution of the system takes much longer.

In[86]:=