# Software tools for medium-scale systems of polynomial equations

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

• your link does not seem to work; it leads to the URL raw.githubusercontent.com/alexflint/polygamy/master/out/… which does not show up. Aug 3, 2014 at 22:35
• as you have a sum squares cost function to optimise, I suppose the most natural method will be one based on semidefinite programming (a.k.a. Lasserre hierarchies). Did you try these? Aug 3, 2014 at 22:39
• @DimaPasechnik the actual link is here: raw.githubusercontent.com/alexflint/polygamy/master/out/…
– Ryan
Aug 4, 2014 at 14:51
• @Ryan thanks for the catch - have updated the post Aug 4, 2014 at 14:53
• would you only be interested to find all solutions or would it also be ok to find one solution? Aug 5, 2014 at 17:09

In Maple you can just do

with(Optimization):
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.

• What is Maple doing internally? For an arbitrary function, random search for a minimum is as good as any other method, so is it using random search? Aug 5, 2014 at 7:18
• maplesoft.com/support/help/Maple/… Aug 5, 2014 at 7:25
• In this case it looks like it will use a general SQP approach, then. I somehow doubt that it will take only a few seconds to find a reasonable approximation to a global minimum of a long sum of squares in 12 variables, but I have not used Maple in a couple of years... Aug 5, 2014 at 7:50
• SQP most probably will terminate in a local minimum. I'd be surprised if any restarts from such a minimum are attempted. Aug 5, 2014 at 8:45
• @Thomas: the solution that I gave suggests that the correct global minimum with exact coefficients is quite likely to be zero. This in turn suggests that the problem may have hidden structure that makes it unusually easy to solve. Aug 5, 2014 at 8:49

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.

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.

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.