numerically solving a system of multivariate polynomials I have a system of polynomial equations which I need to solve numerically. The one Im currently interested in has around 20 variables and a similar number of equations. All coefficients are rational numbers. In the end I am only interested in real solutions with all variables in the range $[-1,+1]$.
The most elegant approach to my knowledge involves finding a Groebner basis (ideally in diagonal form), which then can be solved by (numberically) solving a series of univariate polynomials, which is easy enough. This has worked for me for some smaller problems (around 10 variables, using the 'Singular' software). But the current problem seems to big for this approach (computation runs out of memory after many hours on quite powerful workstation). Any other ideas I could persue?

*

*I know that the runtime of Gröbner depends greatly on the choosen term order. By default, Singular uses 'degrevlex', which is usually good, but might possibly be improved by introducing some weights to the variables. Is there some heuristic to find good weights?


*Im only interested in finding some solution as numerical approximations. I dont need to find all of them or prove that I found all of them or anything like that (which a Gröbner basis would allow).


*I did try using Newtons Method to find a root. But if I start from a randomly chosen point, this is extremely unstable. Is there any more stable method for this? (most literature I could find deals with multi-dimension minimization, not with root-finding).
Any ideas would be welcome.
 A: In my former research group, we have been confronted with the same problem and then I fully agree with all your statements.
By the end (to make the story short), we concluded that the best way was to use optimization to minimize
$$\Phi=\sum_{n=1}^p \big[\text{equation}_n\big]^2$$
(no need to change any equation and no need to use Gröbner basis which are overkilling).
Concerning the problem of bounds, most otpimizers allow bound constraints (these are the simplest to handle). If your does not, for $a \leq x \leq b$, use the transformation
$$x=a +\frac {b-a}{1+e^{-X}}$$
The last question is the starting point : in our case, we knew that there was only one acceptable solution. So, we used to make multiple runs with randomly selected guesses first and then polish the solution.
The advantage of this approach is that, with polynomial equations, you can very easily write the analytical Jacobian and Hessian.
A: Maybe you can provide an example of such a system so that one can introduce you an approach that works the best for your case.
Anyways, one way to solve a system of polynomial equations numerically and not very time consuming is to use numerical homotopy method. It does not involve Grobner basis.
There is a nice book "Numerically solving polynomial systems with Bertini" by D. J. Bates, J. D. Hauenstein, A. J. Sommese, and C. W. Wampler, but you should be able to find about it in other basic numerical algebraic geometry books. If you are a Maple user, then you can use Homotopy command from RootFinding package (see https://www.maplesoft.com/support/help/maple/view.aspx?path=RootFinding%2fHomotopy). Another choice is the software Bertini which is free and has a very small size (less than 4 MB). Its homepage is https://bertini.nd.edu/ and if you are a Windows user, then here is a short video on how to install it on Windows and two simple examples https://youtu.be/9puCD-8XZxE . One thing that you may need to be careful is that if the coordinates of your solutions (for the 20 variables) are of different scales, say one is of scale $10^{-1}$ while the other one of scale $10^{-30}$, then you may need to play a bit with your input system such as rescaling. Otherwise, if they are in similar scale range, it should be usually ok. Also be careful of singularity issues ^_^
