System of polynomial equations

A trigonometric equation is an equation g = 0 where g is a trigonometric polynomial. Such an equation may be converted into a polynomial system by expanding the sines and cosines in it (using sum and difference formulas), replacing sin(x) and cos(x) by two new variables s and c and adding the new equation s² + c² – 1 = 0.

For example, because of the identity

${\displaystyle \cos(3x)=4\cos ^{3}(x)-3\cos(x),}$

solving the equation

${\displaystyle \sin ^{3}(x)+\cos(3x)=0}$

is equivalent to solving the polynomial system

${\displaystyle {\begin{cases}s^{3}+4c^{3}-3c&=0\\s^{2}+c^{2}-1&=0.\end{cases}}}$

For each solution (c₀, s₀) of this system, there is a unique solution x of the equation such that 0 ≤ x < 2π.

In the case of this simple example, it may be unclear whether the system is, or not, easier to solve than the equation. On more complicated examples, one lacks systematic methods for solving directly the equation, while software are available for automatically solving the corresponding system.

When solving a system over a finite field k with q elements, one is primarily interested in the solutions in k. As the elements of k are exactly the solutions of the equation x^q – x = 0, it suffices, for restricting the solutions to k, to add the equation x_i^q – x_i = 0 for each variable x_i.

The elements of an algebraic number field are usually represented as polynomials in a generator of the field which satisfies some univariate polynomial equation. To work with a polynomial system whose coefficients belong to a number field, it suffices to consider this generator as a new variable and to add the equation of the generator to the equations of the system. Thus solving a polynomial system over a number field is reduced to solving another system over the rational numbers.

For example, if a system contains ${\displaystyle {\sqrt {2}}}$, a system over the rational numbers is obtained by adding the equation r₂² – 2 = 0 and replacing ${\displaystyle {\sqrt {2}}}$ by r₂ in the other equations.

In the case of a finite field, the same transformation allows always supposing that the field k has a prime order.

The usual way of representing the solutions is through zero-dimensional regular chains. Such a chain consists of a sequence of polynomials f₁(x₁), f₂(x₁, x₂), ..., f_n(x₁, ..., x_n) such that, for every i such that 1 ≤ i ≤ n

• f_i is a polynomial in x₁, ..., x_i only, which has a degree d_i > 0 in x_i;
• the coefficient of x_i^d_i in f_i is a polynomial in x₁, ..., x_{i −1} which does not have any common zero with f₁, ..., f_{i − 1}.

To such a regular chain is associated a triangular system of equations

${\displaystyle {\begin{cases}f_{1}(x_{1})=0\\f_{2}(x_{1},x_{2})=0\\\quad \vdots \\f_{n}(x_{1},x_{2},\ldots ,x_{n})=0.\end{cases}}}$

The solutions of this system are obtained by solving the first univariate equation, substituting the solutions in the other equations, then solving the second equation which is now univariate, and so on. The definition of regular chains implies that the univariate equation obtained from f_i has degree d_i and thus that the system has d₁ ... d_n solutions, provided that there is no multiple root in this resolution process (fundamental theorem of algebra).

Every zero-dimensional system of polynomial equations is equivalent (i.e. has the same solutions) to a finite number of regular chains. Several regular chains may be needed, as it is the case for the following system which has three solutions.

${\displaystyle {\begin{cases}x^{2}-1=0\\(x-1)(y-1)=0\\y^{2}-1=0.\end{cases}}}$

There are several algorithms for computing a triangular decomposition of an arbitrary polynomial system (not necessarily zero-dimensional)[4] into regular chains (or regular semi-algebraic systems).

There is also an algorithm which is specific to the zero-dimensional case and is competitive, in this case, with the direct algorithms. It consists in computing first the Gröbner basis for the graded reverse lexicographic order (grevlex), then deducing the lexicographical Gröbner basis by FGLM algorithm[5] and finally applying the Lextriangular algorithm.[6]

This representation of the solutions are fully convenient for coefficients in a finite field. However, for rational coefficients, two aspects have to be taken care of:

• The output may involve huge integers which may make the computation and the use of the result problematic.
• To deduce the numeric values of the solutions from the output, one has to solve univariate polynomials with approximate coefficients, which is a highly unstable problem.

The first issue has been solved by Dahan and Schost:[7][8] Among the sets of regular chains that represent a given set of solutions, there is a set for which the coefficients are explicitly bounded in terms of the size of the input system, with a nearly optimal bound. This set, called equiprojectable decomposition, depends only on the choice of the coordinates. This allows the use of modular methods for computing efficiently the equiprojectable decomposition.[9]

The second issue is generally solved by outputting regular chains of a special form, sometimes called shape lemma, for which all d_i but the first one are equal to 1. For getting such regular chains, one may have to add a further variable, called separating variable, which is given the index 0. The rational univariate representation, described below, allows computing such a special regular chain, satisfying Dahan–Schost bound, by starting from either a regular chain or a Gröbner basis.

The rational univariate representation or RUR is a representation of the solutions of a zero-dimensional polynomial system over the rational numbers which has been introduced by F. Rouillier.[10]

A RUR of a zero-dimensional system consists in a linear combination x₀ of the variables, called separating variable, and a system of equations[11]

${\displaystyle {\begin{cases}h(x_{0})=0\\x_{1}=g_{1}(x_{0})/g_{0}(x_{0})\\\quad \vdots \\x_{n}=g_{n}(x_{0})/g_{0}(x_{0}),\end{cases}}}$

where h is a univariate polynomial in x₀ of degree D and g₀, ..., g_n are univariate polynomials in x₀ of degree less than D.

Given a zero-dimensional polynomial system over the rational numbers, the RUR has the following properties.

• All but a finite number linear combinations of the variables are separating variables.
• When the separating variable is chosen, the RUR exists and is unique. In particular h and the g_i are defined independently of any algorithm to compute them.
• The solutions of the system are in one-to-one correspondence with the roots of h and the multiplicity of each root of h equals the multiplicity of the corresponding solution.
• The solutions of the system are obtained by substituting the roots of h in the other equations.
• If h does not have any multiple root then g₀ is the derivative of h.

For example, for the system in the previous section, every linear combination of the variable, except the multiples of x, y and x + y, is a separating variable. If one chooses t = x – y/2 as a separating variable, then the RUR is

${\displaystyle {\begin{cases}t^{3}-t=0\\x={\frac {t^{2}+2t-1}{3t^{2}-1}}\\y={\frac {t^{2}-2t-1}{3t^{2}-1}}.\\\end{cases}}}$

The RUR is uniquely defined for a given separating variable, independently of any algorithm, and it preserves the multiplicities of the roots. This is a notable difference with triangular decompositions (even the equiprojectable decomposition), which, in general, do not preserve multiplicities. The RUR shares with equiprojectable decomposition the property of producing an output with coefficients of relatively small size.

For zero-dimensional systems, the RUR allows retrieval of the numeric values of the solutions by solving a single univariate polynomial and substituting them in rational functions. This allows production of certified approximations of the solutions to any given precision.

Moreover, the univariate polynomial h(x₀) of the RUR may be factorized, and this gives a RUR for every irreducible factor. This provides the prime decomposition of the given ideal (that is the primary decomposition of the radical of the ideal). In practice, this provides an output with much smaller coefficients, especially in the case of systems with high multiplicities.

Contrarily to triangular decompositions and equiprojectable decompositions, the RUR is not defined in positive dimension.

The general numerical algorithms which are designed for any system of nonlinear equations work also for polynomial systems. However the specific methods will generally be preferred, as the general methods generally do not allow one to find all solutions. In particular, when a general method does not find any solution, this is usually not an indication that there is no solution.

Nevertheless, two methods deserve to be mentioned here.

• Newton's method may be used if the number of equations is equal to the number of variables. It does not allow one to find all the solutions nor to prove that there is no solution. But it is very fast when starting from a point which is close to a solution. Therefore, it is a basic tool for the homotopy continuation method described below.
• Optimization is rarely used for solving polynomial systems, but it succeeded, circa 1970, in showing that a system of 81 quadratic equations in 56 variables is not inconsistent.[12] With the other known methods, this remains beyond the possibilities of modern technology, as of 2022. This method consists simply in minimizing the sum of the squares of the equations. If zero is found as a local minimum, then it is attained at a solution. This method works for overdetermined systems, but outputs an empty information if all local minimums which are found are positive.

This is a semi-numeric method which supposes that the number of equations is equal to the number of variables. This method is relatively old but it has been dramatically improved in the last decades.[13]

This method divides into three steps. First an upper bound on the number of solutions is computed. This bound has to be as sharp as possible. Therefore, it is computed by, at least, four different methods and the best value, say ${\displaystyle N}$, is kept.

In the second step, a system ${\displaystyle g_{1}=0,\,\ldots ,\,g_{n}=0}$ of polynomial equations is generated which has exactly ${\displaystyle N}$ solutions that are easy to compute. This new system has the same number ${\displaystyle n}$ of variables and the same number ${\displaystyle n}$ of equations and the same general structure as the system to solve, ${\displaystyle f_{1}=0,\,\ldots ,\,f_{n}=0}$.

Then a homotopy between the two systems is considered. It consists, for example, of the straight line between the two systems, but other paths may be considered, in particular to avoid some singularities, in the system

${\displaystyle (1-t)g_{1}+tf_{1}=0,\,\ldots ,\,(1-t)g_{n}+tf_{n}=0}$.

The homotopy continuation consists in deforming the parameter ${\displaystyle t}$ from 0 to 1 and following the ${\displaystyle N}$ solutions during this deformation. This gives the desired solutions for ${\displaystyle t=1}$. Following means that, if ${\displaystyle t_{1}<t_{2}}$, the solutions for ${\displaystyle t=t_{2}}$ are deduced from the solutions for ${\displaystyle t=t_{1}}$ by Newton's method. The difficulty here is to well choose the value of ${\displaystyle t_{2}-t_{1}:}$ Too large, Newton's convergence may be slow and may even jump from a solution path to another one. Too small, and the number of steps slows down the method.

To deduce the numeric values of the solutions from a RUR seems easy: it suffices to compute the roots of the univariate polynomial and to substitute them in the other equations. This is not so easy because the evaluation of a polynomial at the roots of another polynomial is highly unstable.

The roots of the univariate polynomial have thus to be computed at a high precision which may not be defined once for all. There are two algorithms which fulfill this requirement.

• Aberth method, implemented in MPSolve computes all the complex roots to any precision.
• Uspensky's algorithm of Collins and Akritas,[14] improved by Rouillier and Zimmermann [15] and based on Descartes' rule of signs. This algorithms computes the real roots, isolated in intervals of arbitrary small width. It is implemented in Maple (functions fsolve and RootFinding[Isolate]).