Software tools for medium-scale systems of polynomial equations

Question

I am attempting to find all real solutions of a system of 12 polynomial equations in 12 unknowns. The equations each have total degree 6 and contain up to 1700 terms. I am only interested in real solutions. The equations were derived as the gradients of a sum-of-squares cost function, which I am attempting to find all global optima of. I believe there are a finite number of real solutions but I have not confirmed this yet. I have floating point coefficients and I'm looking for numerical solutions (as opposed to symbolic solutions).

Which software packages (and which functions specifically) are generally most promising to solve such a problem?

I am aware of various functions in Maple, Matlab, and Mathematica that can solve systems of polynomial equations but there are a large number of options in each software package and I'm interested in advice on where I should be looking first for problems of this scale.

A numerical dump of the cost function is here:

https://raw.githubusercontent.com/alexflint/polygamy/master/out/epipolar_accel_bezier3/cost.txt

Answer 1

In Maple you can just do

with(Optimization):
g := (your function):
Minimize(g,iterationlimit = 200);

On my machine this takes only about 1.5 seconds to return the following:

[2.35579022955789696*10^(-9), 
[x0 = .696531801759957, x1 = .286105658731833, x10 = .342973444356395, 
 x11 = .728732510532874, x2 = .226824733028582, x3 = .551288843437034, 
 x4 = .719479494442298, x5 = .423120942717389, x6 = .980635895595386, 
 x7 = .684727337329935, x8 = .480773372241607, x9 = .391860913480735]]

If you ask for 30 digits of accuracy (so that Maple cannot just use hardware floats) then it gets a lot slower, but still only 40 seconds.

Incidentally, I tried this first as just Minimize(g), which gave me an answer with a warning that the maximum iteration limit had been reached. The default limit is 50 iterations, so I tried 200 and the warning went away.

Also, there is another function called LSSolve which is specifically for the case where your objective function is a sum of squares. I did not use that because I do not know how to express your objective in the relevant form.

Answer 2

To solve a polynomial system, I would try Bertini which is a homotopy-continuation numerical solver that parallelizes extremely well. You can also try to attack the optimization problem directly with semi-definite programming as explained by Dima Pasechnik.

Answer 3

As far as I can see, you'd like to find a global minimum of $F_0(x)=\sum_{k=1}^{M_0} f_k(x)^2$, where $x=(x_1,\dots,x_{12})$. Equivalently, the problem is to find maximal $\mu$ so that $F(\mu,x):=F_0(x)-\mu\geq 0$ for all $x_1,\dots,x_{12}$.

Now, a sufficient condition for $F(\mu,x)\geq 0$ is that $F(\mu,x)=\sum_{k=1}^{M} g_k(\mu,x)^2$. The latter condition can be rewritten in the form $F(\mu,x)=\tilde{x}^\top A\tilde{x}\geq 0$ for $\tilde{x}$ being the vector of monomials in $\mu$ and $x$ of degree at most 3, and $A$ being a symmetric positive semidefinite matrix. The condition $F(\mu,x)=\tilde{x}^\top A\tilde{x}$ is just a bunch of linear equations derived from the coefficients of the polynomials on the RHS and LHS being equal.

This can be converted into a "semidefinite programming problem" (SDP) to find such a maximal $\mu$, as proposed in Lasserre's paper; software to do this, and more, is readily available, too; say, YALMIP. The answer is, however, only a lower bound on true value of $\mu$. But often enough it is the exact value what one will get this way.

Answer 4

I believe INTLAB, an interval analysis package in MATLAB might be worth considering. It is explained here. Specifically, section 7.4 All Solutions of a Nonlinear System (with an implementation in Appendix A.3) is worth taking a look at.