View source for System of polynomial equations

{{Short description|Roots of multiple multivariate polynomials}}
A '''system of polynomial equations''' (sometimes simply a '''polynomial system''') is a set of simultaneous equations {{math|1=''f''<sub>1</sub> = 0, ..., ''f''<sub>''h''</sub> = 0}} where the {{math|''f''<sub>''i''</sub>}} are [[Polynomial|polynomial]]s in several variables, say {{math|''x''<sub>1</sub>, ..., ''x''<sub>''n''</sub>}}, over some [[Field (mathematics)|field]] {{mvar|''k''}}.

A ''solution'' of a polynomial system is a set of values for the {{math|''x''<sub>''i''</sub>}}s which belong to some algebraically closed field extension {{math|''K''}} of {{math|''k''}}, and make all equations true. When {{math|''k''}} is the field of [[Rational number|rational number]]s, {{math|''K''}} is generally assumed to be the field of [[Complex number|complex number]]s, because each solution belongs to a field extension of {{math|''k''}}, which is isomorphic to a subfield of the complex numbers.

This article is about the methods for solving, that is, finding all solutions or describing them. As these methods are designed for being implemented in a computer, emphasis is given on fields {{math|''k''}} in which computation (including equality testing) is easy and efficient, that is the field of [[Rational number|rational number]]s and [[Finite field|finite field]]s.

Searching for solutions that belong to a specific set is a problem which is generally much more difficult, and is outside the scope of this article, except for the case of the solutions in a given finite field. For the case of solutions of which all components are integers or rational numbers, see [[Diophantine equation]].

==Definition==
[[File:BarthSextic.png|thumb|The numerous singular points of the Barth sextic are the solutions of a polynomial system]]
A simple example of a system of polynomial equations is
:<math> \begin{align} x^2 + y^2 - 5&= 0 \\ xy - 2 &= 0. \end{align}</math>
Its solutions are the four pairs {{math|1=(''x'', ''y'') = (1, 2), (2, 1), (-1, -2), (-2, -1)}}. These solutions can easily be checked by substitution, but more work is needed for proving that there are no other solutions.

The subject of this article is the study of generalizations of such an examples, and the description of the methods that are used for computing the solutions.

A ''system of polynomial equations,''  or ''polynomial system'' is a collection of equations
:<math> \begin{align}
f_1\left(x_1, \ldots, x_m \right) &= 0 \\
& \;\;\vdots \\
f_n\left(x_1, \ldots, x_m \right) &= 0,
\end{align} </math>
where each {{math|''f<sub>h</sub>''}} is a [[Polynomial|polynomial]] in the [[Indeterminate (variable)|indeterminates]] {{math|''x''<sub>1</sub>, ..., ''x''<sub>''m''</sub>}}, with integer coefficients, or coefficients in some fixed [[Field (mathematics)|field]], often the field of [[Rational number|rational number]]s or a [[Finite field|finite field]].<ref name=Bates_p4>{{harvnb|Bates|Sommese|Hauenstein|Wampler|2013|p=4}}</ref> Other fields of coefficients, such as the [[Real number|real number]]s, are less often used, as their elements cannot be represented in a computer (only approximations of real numbers can be used in computations, and these approximations are always rational numbers).

A ''solution'' of a polynomial system is a [[Tuple|tuple]] of values of {{math|(''x''<sub>1</sub>, ..., ''x''<sub>''m''</sub>)}} that satisfies all equations of the polynomial system. The solutions are sought in the [[Complex number|complex number]]s, or more generally in an [[Algebraically closed field|algebraically closed field]] containing the coefficients. In particular, in characteristic zero, all [[Complex number|complex]] solutions are sought. Searching for the [[Real number|real]] or [[Rational number|rational]] solutions are much more difficult problems that are not considered in this article.

The set of solutions is not always finite; for example, the solutions of the system
:<math> \begin{align} x(x-1) &= 0 \\ x(y-1) &= 0 \end{align}</math>
are a point {{math|1=(''x'',''y'') = (1,1)}} and a line {{math|1=''x'' = 0}}.<ref>{{harvnb|Bates|Sommese|Hauenstein|Wampler|2013|p=8}}</ref> Even when the solution set is finite, there is, in general, no [[Closed-form expression|closed-form expression]] of the solutions (in the case of a single equation, this is [[Abel–Ruffini theorem]]).

The [[Barth surface]], shown in the figure is the geometric representation of the solutions of a polynomial system reduced to a single equation of degree 6 in 3 variables. Some of its numerous [[Singular point of an algebraic variety|singular points]] are visible on the image. They are the solutions of a system of 4 equations of degree 5 in 3 variables. Such an [[Overdetermined system|overdetermined system]] has no solution in general (that is if the coefficients are not specific). If it has a finite number of solutions, this number is at most {{math|1=5{{sup|3}} = 125}}, by [[Bézout's theorem]]. However, it has been shown that, for the case of the singular points of a surface of degree 6, the maximum number of solutions is 65, and is reached by the Barth surface.

==Basic properties and definitions==
A system is [[Overdetermined system|overdetermined]] if the number of equations is higher than the number of variables. A system is inconsistent if it has no [[Complex number|complex]] solution (or, if the coefficients are not complex numbers, no solution in an [[Algebraically closed field|algebraically closed field]] containing the coefficients). By [[Hilbert's Nullstellensatz]] this means that 1 is a [[Linear combination|linear combination]] (with polynomials as coefficients) of the first members of the equations. Most but not all overdetermined systems, when constructed with random coefficients, are inconsistent. For example, the system {{math|1=''x''<sup>3</sup> – 1 = 0, ''x''<sup>2</sup> – 1 = 0}} is overdetermined (having two equations but only one unknown), but it is not inconsistent since it has the solution {{math|1=''x'' = 1}}.

A system is [[Underdetermined system|underdetermined]] if the number of equations is lower than the number of the variables. An underdetermined system is either inconsistent or has infinitely many complex solutions (or solutions in an [[Algebraically closed field|algebraically closed field]] that contains the coefficients of the equations). This is a non-trivial result of [[Commutative algebra|commutative algebra]] that involves, in particular, [[Hilbert's Nullstellensatz]] and [[Krull's principal ideal theorem]].

A system is [[Zero-dimensional space|zero-dimensional]] if it has a finite number of complex solutions (or solutions in an algebraically closed field). This terminology comes from the fact that the [[Algebraic variety|algebraic variety]] of the solutions has [[Dimension of an algebraic variety|dimension]] zero. A system with infinitely many solutions is said to be ''positive-dimensional''.

A zero-dimensional system with as many equations as variables is sometimes said to be ''well-behaved''.<ref>Songxin Liang, J. Gerhard, D.J. Jeffrey, G. Moroz, ''[https://www.researchgate.net/profile/David_Stoutemyer/publication/235941679_Unit_normalization_of_multinomials_over_Gaussian_integers/links/02e7e532ddec858b2a000000.pdf#page=2 A Package for Solving Parametric Polynomial Systems]''. Communications in Computer Algebra (2009)</ref>
[[Bézout's theorem]] asserts that a well-behaved system whose equations have degrees {{math|''d''<sub>1</sub>, ..., ''d''<sub>''n''</sub>}} has at most {{math|''d''<sub>1</sub>⋅⋅⋅''d''<sub>''n''</sub>}} solutions. This bound is sharp. If all the degrees are equal to {{math|''d''}}, this bound becomes {{math|''d''<sup>''n''</sup>}} and is exponential in the number of variables. (The [[Fundamental theorem of algebra|fundamental theorem of algebra]] is the special case {{math|1= ''n'' = 1}}.)

This exponential behavior makes solving polynomial systems difficult and explains why there are few solvers that are able to automatically solve systems with Bézout's bound higher than, say, 25 (three equations of degree 3 or five equations of degree 2 are beyond this bound).{{Citation needed|date=July 2018}}

== What is solving? ==
The first thing to do for solving a polynomial system is to decide whether it is inconsistent, zero-dimensional or positive dimensional. This may be done by the computation of a [[Gröbner basis]] of the left-hand sides of the equations. The system is ''inconsistent'' if this Gröbner basis is reduced to 1. The system is ''zero-dimensional'' if, for every variable there is a [[Gröbner basis|leading monomial]] of some element of the Gröbner basis which is a pure power of this variable. For this test, the best [[Monomial order|monomial order]] (that is the one which leads generally to the fastest computation) is usually the [[Monomial order#Graded reverse lexicographic order|graded reverse lexicographic]] one (grevlex).

If the system is ''positive-dimensional'', it has infinitely many solutions. It is thus not possible to enumerate them. It follows that, in this case, solving may only mean "finding a description of the solutions from which the relevant properties of the solutions are easy to extract". There is no commonly accepted such description. In fact there are many different "relevant properties", which involve almost every subfield of [[Algebraic geometry|algebraic geometry]].

A natural example of such a question concerning positive-dimensional systems is the following: ''decide if a polynomial system over the [[Rational number|rational number]]s has a finite number of real solutions and compute them''. A generalization of this question is ''find at least one solution in each connected component of the set of real solutions of a polynomial system''.  The classical algorithm for solving these question is [[Cylindrical algebraic decomposition|cylindrical algebraic decomposition]], which has a [[Double exponential function|doubly exponential]] [[Computational complexity|computational complexity]] and therefore cannot be used in practice, except for very small examples.

For zero-dimensional systems, solving consists of computing all the solutions. There are two different ways of outputting the solutions. The most common way is possible only for real or complex solutions, and consists of outputting numeric approximations of the solutions. Such a solution is called ''numeric''. A solution is ''certified'' if it is provided with a bound on the error of the approximations, and if this bound separates the different solutions.

The other way of representing the solutions is said to be ''algebraic''. It uses the fact that, for a zero-dimensional system, the solutions belong to the algebraic closure of the field ''k'' of the coefficients of the system. There are several ways to represent the solution in an algebraic closure, which are discussed below. All of them allow one to compute a numerical approximation of the solutions by solving one or several univariate equations. For this computation, it is preferable to use a representation that involves solving only one univariate polynomial per solution, because computing the roots of a polynomial which has approximate coefficients [[Wilkinson's polynomial|is a highly unstable problem]].

==Extensions==
===Trigonometric equations===
A trigonometric equation is an equation {{math|1=''g'' = 0}} where {{math|''g''}} is a [[Trigonometric polynomial|trigonometric polynomial]]. Such an equation may be converted into a polynomial system by expanding the sines and cosines in it (using [[Trigonometric functions#Sum and difference formulas|sum and difference formulas]]), replacing {{math|sin(''x'')}} and {{math|cos(''x'')}} by two new variables {{math|''s''}} and {{math|''c''}} and adding the new equation {{math|1=''s''<sup>2</sup> + ''c''<sup>2</sup> – 1 = 0}}.

For example, because of the identity
:<math>\cos(3x)=4\cos^3(x)-3\cos(x),</math>
solving the equation
:<math> \sin^3(x)+\cos(3x)=0 </math>
is equivalent to solving the polynomial system
:<math> \begin{cases}
s^3+4c^3-3c&=0\\
s^2+c^2-1&=0.
\end{cases}
</math>

For each solution {{math|(''c''{{sub|0}}, ''s''{{sub|0}})}} of this system, there is a unique solution {{mvar|x}} of the equation such that {{math|0 ≤ ''x'' < 2{{pi}}}}.

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.

=== Solutions in a finite field ===
When solving a system over a finite field {{math|''k''}} with {{math|''q''}} elements, one is primarily interested in the solutions in {{math|''k''}}. As the elements of {{math|''k''}} are exactly the solutions of the equation {{math|1=''x''<sup>''q''</sup> – ''x'' = 0}}, it suffices, for restricting the solutions to {{math|''k''}}, to add the equation {{math|1=''x''<sub>''i''</sub><sup>''q''</sup> – ''x''<sub>''i''</sub> = 0}} for each variable&nbsp;{{math|''x''<sub>''i''</sub>}}.

=== Coefficients in a number field or in a finite field with non-prime order ===

The elements of an [[Algebraic number field|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 <math>\sqrt{2}</math>, a system over the rational numbers is obtained by adding the equation {{math|1=''r''<sub>2</sub><sup>2</sup> – 2 = 0}} and replacing <math>\sqrt{2}</math> by {{math|''r''<sub>2</sub>}} in the other equations.

In the case of a finite field, the same transformation allows always supposing that the field {{math|''k''}} has a prime order.

== Algebraic representation of the solutions==
=== Regular chains ===
{{main|Regular chain}}
The usual way of representing the solutions is through zero-dimensional regular chains. Such a chain consists of a sequence of polynomials {{math|''f''<sub>1</sub>(''x''<sub>1</sub>)}}, {{math|''f''<sub>2</sub>(''x''<sub>1</sub>, ''x''<sub>2</sub>)}}, ..., {{math|''f''<sub>''n''</sub>(''x''<sub>1</sub>, ..., ''x''<sub>''n''</sub>)}} such that, for every {{math|''i''}} such that {{math|1 ≤ ''i'' ≤ ''n''}}
* {{math|''f''<sub>''i''</sub>}} is a polynomial in {{math|''x''<sub>1</sub>, ..., ''x''<sub>''i''</sub>}} only, which has a degree {{math|''d''<sub>''i''</sub> > 0}} in {{math|''x''<sub>''i''</sub>}};
* the coefficient of {{math|''x''<sub>''i''</sub><sup>''d''<sub>''i''</sub></sup>}} in {{math|''f''<sub>''i''</sub>}} is a polynomial in {{math|''x''<sub>1</sub>, ..., ''x''<sub>''i'' −1 </sub>}} which does not have any common zero with {{math|''f''<sub>1</sub>}}, ..., {{math|''f''<sub>''i'' − 1</sub>}}.

To such a [[Regular chain|regular chain]] is associated a ''triangular system of equations''
:<math>
\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}
</math>

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 {{math|''f''<sub>''i''</sub>}} has degree {{math|''d''<sub>''i''</sub>}} and thus that the system has {{math|''d''<sub>1</sub> ... ''d''<sub>''n''</sub>}} solutions, provided that there is no multiple root in this resolution process ([[Fundamental theorem of algebra|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.
:<math>
\begin{cases}
x^2-1=0\\
(x-1)(y-1)=0\\
y^2-1=0.
\end{cases}
</math>

There are several algorithms for computing a [[Triangular decomposition|triangular decomposition]] of an arbitrary polynomial system (not necessarily zero-dimensional)<ref>{{cite journal | last1 = Aubry | first1 = P. | last2 = Maza | first2 = M. Moreno | title = Triangular Sets for Solving Polynomial Systems: a Comparative Implementation of Four Methods | journal = J. Symb. Comput. | volume = 28 | issue = 1–2| pages = 125–154 | year = 1999 | doi=10.1006/jsco.1999.0270| doi-access = free }}</ref> into [[Regular chain|regular chain]]s (or [[Regular semi-algebraic system|regular semi-algebraic system]]s).

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 [[Monomial order|graded reverse lexicographic order (grevlex)]], then deducing the lexicographical Gröbner basis by FGLM algorithm<ref>{{cite journal | last1 = Faugère | first1 = J.C. | last2 = Gianni | first2 = P. | last3 = Lazard | first3 = D. | last4 = Mora | first4 = T. | title = Efficient Computation of Zero-Dimensional Gröbner Basis by Change of Ordering | journal = Journal of Symbolic Computation | volume = 16 | issue = 4| year = 1993 | doi=10.1006/jsco.1993.1051 | pages=329–344| doi-access = free }}</ref> and finally applying the Lextriangular algorithm.<ref>{{cite journal | last1 = Lazard | first1 = D. | title = Solving zero-dimensional algebraic systems | journal = Journal of Symbolic Computation | volume = 13 | issue = 2| year = 1992 | doi=10.1016/S0747-7171(08)80086-7 | pages=117–131| doi-access =  }}</ref>

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:<ref>Xavier Dahan and Eric Schost. ''[http://www.csd.uwo.ca/~eschost/publications/DaSc.ps Sharp Estimates for Triangular Sets]''. Moreover, recent algorithms for decomposing polynomial systems into triangular decompositions produce regular chains with coefficients matching the results of Dahan and Schost.  In proc. ISSAC'04, pages 103--110, ACM Press, 2004</ref><ref>
{{cite book
| title= Proceedings of ISAAC 2005
| year = 2005
| first1=Xavier | last1=Dahan
| first2=Marc  | last2=Moreno Maza
| first3=Eric | last3=Schost
| first4=Wenyuan  | last4=Wu
| first5=Yuzhen  | last5=Xie
| chapter=Lifting techniques for triangular decompositions
| pages=108–105
| publisher=ACM Press
| chapter-url=http://www.cecm.sfu.ca/~monaganm/research/MITACS/papers/maza.pdf
}}</ref> 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 arithmetic|modular methods]] for computing efficiently the equiprojectable decomposition.<ref>Changbo Chen and Marc Moreno-Maza. ''[http://www.csd.uwo.ca/~moreno/Publications/Algorithms_for_computing_triangular_decomposition_of_polynomial_systems.pdf Algorithms for Computing Triangular Decomposition of Polynomial Systems]''.In proc. ISSAC'2011, pages 83-90, ACM Press, 2011 and Journal of Symbolic Computation (to appear)</ref>

The second issue is generally solved by outputting regular chains of a special form, sometimes called ''shape lemma'', for which all {{math|''d''<sub>''i''</sub>}} but the first one are equal to {{math|1}}. For getting such regular chains, one may have to add a further variable, called ''separating variable'', which is given the index {{math|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.

===Rational univariate representation===
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.<ref>
{{cite journal
| first=Fabrice
| last=Rouillier
| title=Solving Zero-Dimensional Systems Through the Rational Univariate Representation
|journal=Appl. Algebra Eng. Commun. Comput.
| volume=9
| pages=433–461
| number=9
| year=1999
| doi=10.1007/s002000050114
| s2cid=25579305
}}</ref>

A RUR of a zero-dimensional system consists in a linear combination {{math|''x''<sub>0</sub>}} of the variables, called ''separating variable'', and a system of equations<ref>{{cite book
 | author = Saugata Basu |author2=Richard Pollack |author3=Marie-Françoise Roy
 | year = 2006
 | title = Algorithms in real algebraic geometry, chapter 12.4
 | publisher = [[Company:Springer Science+Business Media|Springer-Verlag]]
 | url = http://perso.univ-rennes1.fr/marie-francoise.roy/bpr-ed2-posted1.html
}}</ref>
:<math>
\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}
</math>
where {{math|''h''}} is a univariate polynomial in {{math|''x''<sub>0</sub>}} of degree {{math|''D''}} and {{math|''g''<sub>0</sub>, ..., ''g''<sub>n</sub>}} are univariate polynomials in {{math|''x''<sub>0</sub>}} of degree less than {{math|''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 {{math|''h''}} and the {{math|''g''<sub>''i''</sub>}} are defined independently of any algorithm to compute them.
* The solutions of the system are in one-to-one correspondence with the roots of {{math|''h''}} and the [[Multiplicity (mathematics)|multiplicity]] of each root of {{math|''h''}} equals the multiplicity of the corresponding solution.
* The solutions of the system are obtained by substituting the roots of {{math|''h''}} in the other equations.
* If {{math|''h''}} does not have any multiple root then {{math|''g''<sub>0</sub>}} is the [[Formal derivative|derivative]] of {{math|''h''}}.

For example, for the system in the previous section, every linear combination of the variable, except the multiples of {{math|''x''}}, {{math|''y''}} and {{math|''x'' + ''y''}}, is a separating variable. If one chooses {{math|1=''t'' = {{sfrac|''x'' – ''y''|2}}}} as a separating variable, then the RUR is
:<math>
\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}
</math>

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 {{math|''h''(''x''<sub>0</sub>)}} 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|primary decomposition]] of the [[Radical of an ideal|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.

==Solving numerically{{anchor|Algorithms for numerically solving}}==
===General solving algorithms===
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.<ref>{{cite journal | last1 = Lazard | first1 = Daniel | year = 2009| title = Thirty years of Polynomial System Solving, and now? | journal = J. Symb. Comput. | volume = 44 | issue = 3| page = 2009 | doi=10.1016/j.jsc.2008.03.004| doi-access = free }}</ref> With the other known methods, this remains beyond the possibilities of modern technology, {{as of| 2022|lc=y}}. 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.

=== Homotopy continuation method ===
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.<ref name="Vers99">{{cite journal
| first=Jan | last=Verschelde
| title=Algorithm 795: PHCpack: A general-purpose solver for polynomial systems by homotopy continuation
| journal=ACM Transactions on Mathematical Software
| volume=25
| issue=2
| pages=251–276
| year=1999
| url=http://eaton.math.rpi.edu/csums/papers/homotopy/verschelde.pdf
| doi=10.1145/317275.317286
| s2cid=15485257
}}</ref>

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 <math>N</math>, is kept.

In the second step, a system <math>g_1=0,\, \ldots,\, g_n=0</math> of polynomial equations is generated which has exactly <math>N</math> solutions that are easy to compute. This new system has the same number <math>n</math> of variables and the same number <math>n</math> of equations and the same general structure as the system to solve, <math>f_1=0,\, \ldots,\, f_n=0</math>.

Then a [[Homotopy|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

:<math>(1-t)g_1+t f_1=0,\, \ldots,\, (1-t)g_n+t f_n=0</math>.

The homotopy continuation consists in deforming the parameter <math>t</math> from 0 to 1 and ''following'' the <math>N</math> solutions during this deformation. This gives the desired solutions for <math>t = 1</math>. ''Following'' means that, if <math>t_1<t_2</math>, the solutions for <math>t=t_2</math> are deduced from the solutions for <math>t=t_1</math> by Newton's method. The difficulty here is to well choose the value of <math>t_2-t_1:</math> 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.

===Numerically solving from the rational univariate representation===

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 [[Software:MPSolve|MPSolve]] computes all the complex roots to any precision.
* Uspensky's algorithm of Collins and Akritas,<ref>George E. Collins and Alkiviadis G. Akritas, ''[http://www.inf.uth.gr/~akritas/articles/1.pdf Polynomial Real Root Isolation Using Descartes' Rule of Signs]''. Proceedings of the 1976 ACM Symposium on Symbolic and Algebraic Computation</ref> improved by Rouillier and Zimmermann <ref>{{cite journal | last1 = Rouillier | first1 = F. | last2 = Zimmerman | first2 = P. | title = Efficient isolation of polynomial's real roots | journal = Journal of Computational and Applied Mathematics | volume = 162 | issue = 1| year = 2004 | doi=10.1016/j.cam.2003.08.015 | pages=33–50| bibcode = 2004JCoAM.162...33R | doi-access = free }}</ref> and based on [[Descartes' rule of signs]]. This algorithms computes the real roots, isolated in intervals of arbitrary small width. It is implemented in [[Software:Maple|Maple]] (functions ''fsolve'' and ''RootFinding[Isolate]'').

==Software packages==
There are at least four software packages which can solve zero-dimensional systems automatically (by automatically, one means that no human intervention is needed between input and output, and thus that no knowledge of the method by the user is needed). There are also several other software packages which may be useful for solving zero-dimensional systems. Some of them are listed after the automatic solvers.

The [[Software:Maple|Maple]] function ''RootFinding[Isolate]'' takes as input any polynomial system over the rational numbers (if some coefficients are floating point numbers, they are converted to rational numbers) and outputs the real solutions represented either (optionally) as intervals of rational numbers or as floating point approximations of arbitrary precision. If the system is not zero dimensional, this is signaled as an error.

Internally, this solver, designed by F. Rouillier computes first a Gröbner basis and then a Rational Univariate Representation from which the required approximation of the solutions are deduced. It works routinely for systems having up to a few hundred complex solutions.

The rational univariate representation may be computed with [[Software:Maple|Maple]] function ''Groebner[RationalUnivariateRepresentation]''.

To extract all the complex solutions from a rational univariate representation, one may use [[Software:MPSolve|MPSolve]], which computes the complex roots of univariate polynomials to any precision. It is recommended to run [[Software:MPSolve|MPSolve]] several times, doubling the precision each time, until solutions remain stable, as the substitution of the roots in the equations of the input variables can be highly unstable.

The second solver is PHCpack,<ref name="Vers99" /><ref>[http://www.math.uic.edu/~jan/download.html Release 2.3.86 of PHCpack]</ref> written under the direction of J. Verschelde. PHCpack implements the homotopy continuation method. This solver computes the isolated complex solutions of polynomial systems having as many equations as variables.

The third solver is Bertini,<ref>{{harvnb|Bates|Sommese|Hauenstein|Wampler|2013}}</ref><ref>[http://bertini.nd.edu/ Bertini: Software for Numerical Algebraic Geometry]</ref> written by D. J. Bates, J. D. Hauenstein, A. J. Sommese, and C. W. Wampler. Bertini uses numerical homotopy continuation with adaptive precision.  In addition to computing zero-dimensional solution sets, both PHCpack and Bertini are capable of working with positive dimensional solution sets.

The fourth solver is the [[Software:Maple|Maple]] library ''RegularChains'', written by Marc Moreno-Maza and collaborators. It contains various functions for solving polynomial systems by means of [[Regular chain|regular chain]]s.

==See also==
*[[Elimination theory]]
*Systems of polynomial inequalities
*[[Triangular decomposition]]
*[[Wu's method of characteristic set]]

== References ==
{{Reflist}}
{{Refbegin}}
*{{cite book |last1=Bates |first1=Daniel J. |last2=Sommese |first2=Andrew J. |last3=Hauenstein |first3=Jonathan D. |last4=Wampler |first4=Charles W. |title=Numerically solving polynomial systems with Bertini |date=2013 |publisher=Society for Industrial and Applied Mathematics |location=Philadelphia |isbn=978-1-61197-269-6}}
*{{cite book |last1=Cox |first1=David |first2=John |last2=Little |first3=Donal |last3=O'Shea |title=Ideals, varieties, and algorithms : an introduction to computational algebraic geometry and commutative algebra |year=1997 |publisher=Springer |location=New York |isbn=978-0387946801 |edition=2nd}}
*{{cite book |last1=Morgan |first1=Alexander |title=Solving polynomial systems using continuation for engineering and scientific problems |date=1987 |publisher=Society for Industrial and Applied Mathematics (SIAM, 3600 Market Street, Floor 6, Philadelphia, PA 19104) |isbn=9780898719031 |edition=SIAM}}
*{{cite book|last=Sturmfels|first=Bernd|title=Solving systems of polynomial equations|year=2002|publisher=American Mathematical Soc.|location=Providence, RI|isbn=0821832514}}
{{Refend}}

{{DEFAULTSORT:Systems Of Polynomial Equations}}
[[Category:Equations]]
[[Category:Algebra]]
[[Category:Computer algebra]]
[[Category:Polynomials]]
[[Category:Algebraic geometry]]

{{Sourceattribution|System of polynomial equations}}