Real Polynomial Systems

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 _x or _x 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 Q_i is or , and (x₁, ..., x_n;y₁, ..., y_m) 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 _{i, j} 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 all of the real polynomial system solving tutorials, we will always assume the system has been transformed to this form.

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

where B (x_k;x₁, ..., x_k-1) is one of

and r₁ and r₂ are algebraic functions (expressed using Root objects or radicals) such that for all x₁, ..., x_k-1 satisfying B (x₁;)B (x₂;x₁)...B (x_k-1;x₁, ..., x_k-2), r₁and r₂ are well defined (that is, denominators and leading terms of Root objects are nonzero), real valued, continuous, and satisfy inequality r₁<r₂.

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.

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.

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.

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.

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:

where C_k is a cell in ^k, r is a continuous algebraic function, r₁ and r₂ are continuous algebraic functions or - or , and r₁<r₂ on C_k.

By an algebraic function we mean a function r:C_k 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.

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 A_n+m of factors of the polynomials present in the quantifier-free part (x₁, ..., x_n;y₁, ..., y_m) of the system (2) and eliminate variables one by one using a projection operator P such that

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

This way the roots of polynomials of A₁, ..., A_n 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.

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 A₁ 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 t₁). Next we lift the cells to cells in ⁿ, one dimension at a time. Suppose we have lifted the cells to ^k. To lift a cell C^k to ^k+1, we find the real roots of A_k+1 with t₁, ..., t_k replaced with the coordinates of the sample point c in C. Since the polynomials of A_k+1 are delineable on C, each root r is a value of a continuous algebraic function at c, and the function can be represented as a p^th root of a polynomial fA_k+1 such that r is the p^th root of f (c, t_k+1). Now the lifting of the cell C to ^k+1 will consist of graphs of these algebraic functions and of the slices of C× between the subsequent graphs. The sample points in each of the new cells will be obtained by adding the k+1^st coordinate to c, 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 k≥n, t_k+1=y_l 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 Q_lylQ_l+1yl+1... Q_mym on C. If Q_l=, we know that the value is True as soon as the truth value of Q_l+1yl+1... Q_mym on one of the lifted cells is True. If Q_l=, we know that the value is False as soon as the truth value of Q_l+1yl+1... Q_mym 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, 24].

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

where x₁, ..., x_n 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 (x₁, ..., x_n) has solutions.

All algorithms used by Mathematica to solve real polynomial decision problems are capable of producing a point satisfying (x₁, ..., x_n) 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.

The primary method that allows Mathematica to solve arbitrary real polynomial decision problems is the Cylindrical Algebraic Decomposition (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]. If there is a variable that appears in the system only quadratically, the quadratic case of Weispfenning's quantifier elimination by virtual substitution algorithm [22, 23] could be used to eliminate the variable. For some examples this gives a substantial speedup; however, quite often it results in a significant slowdown. By default, the algorithm is not used as a preprocessor. Setting the system option QVSPreprocessor in the InequalitySolvingOptions group to True makes Mathematica use it.

There are 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.

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

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.

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 a or b is a nonzero constant, this eliminates the variable y.

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

If all polynomials in the system are at most quadratic in a variable from the innermost quantifier, the variable is eliminated using the quadratic case of Weispfenning's quantifier elimination by virtual substitution algorithm [22, 23]. With the default setting of the system option QuadraticQE, the algorithm is used for Resolve with no variables specified and with at least two parameters present, and for Reduce and Resolve with at least three variables as long as elimination of one variable at most doubles the LeafCount of the system.

The CAD algorithm is used when the previous three 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.

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.

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 the 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.

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

SetSystemOptions["ReduceOptions"->{"option name"->value}].

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 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

SetSystemOptions["InequalitySolvingOptions"->{"option name"->value}].

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 t=0, 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 acc=, we revert to RootReduce if eacc>$MaxPrecision. If t=1, RootReduce is always used. If t=2, 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.

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.

The EquationalConstraintsCAD option specifies whether the projection phase of the CAD algorithm should use equational constraints. With the default setting Automatic, Mathematica uses the projection operator proven correct in [11]. With EquationalConstraintsCAD->True, the smaller but unproven projection operator suggested in [4] is used.

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 x_i we find a univariate polynomial in x_i 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 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.

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.