# Solving a System of Quadratic Equations

I have many polynomial equations in many variables which I want to jointly minimize (in a mean square sense, but you could pick a different reasonable measure which favors anything where all quantities go to zero). For example, I am looking to make the 6 equations below as "small" as possible (a-j are unknown real numbers).

This example probably actually has a solution where all equations are zero, but I also have cases which have no zero solution, so I'd rather not do the "repeatedly eliminate variables and solve for the quadratic root" approach (also, this approach takes too long; is there even any machine which could find a full zero for these equations within 10 minutes?). I'm thinking there might be some software tool that considers the "terrain" smartly and is locally minimizing on many global fronts...or maybe that is impractical. So, is there a free math tool (like Sage) which can minimize things for me (and be certain that no other point is better within some tolerance)? I'm open to theoretical advice, but feel like the options will all look like brute force.

Should I give up if I need to minimize a similar set of equations within 10 minutes on one machine?

a^2 + b^2 + c^2 - 4.52
0.136*a^2 - 0.15*a + d^2 + e^2 + f^2 - 3.84
1.12*a^2 + 0.415*a + b^2 - 2.0*b*d + c^2 - 2.0*c*e + d^2 + e^2 + f^2 - 0.593
0.602*a^2 - 0.0411*a*b - 0.851*a*d - 0.94*a + 0.634*b^2 - 0.489*b*d + 0.219*b + 0.407*d^2 + 0.588*d + h^2 + i^2 + j^2 - 0.612
0.0676*a^2 - 0.0495*a*b + 0.258*a*d + 0.155*a + 0.095*b^2 + 0.14*b*d + 0.157*b + c^2 - 2.0*c*h + 0.407*d^2 + 0.588*d + h^2 + i^2 + j^2 - 2.0
0.813*a^2 - 0.192*a*b - 0.843*a*d - 1.01*a + 0.634*b^2 - 2.03*b*d + 0.302*b + 2.04*d^2 + 0.212*d + e^2 - 2.0*e*h + f^2 - 2.0*f*i + h^2 + i^2 + j^2 - 2.76

• do you talk about real solutions only? Aug 21, 2013 at 7:22
• I am able to compute a Groebner basis for your quadratic equations, but it is too long to display here. Aug 21, 2013 at 8:37
• Real solutions only. @Dietrich Great, but are you able to find one solution from that basis in a reasonable time? Aug 21, 2013 at 14:03
• Have you tried free nonlinear programming solvers? There are many. Also, I suspect you would get more useful answers by asking on scicomp.stackexchange.com (since you are interested in approximate solutions). Oct 7, 2013 at 9:46
• You could also apply gloptipoly to the sum of squares of your quadratics. But you will need MATLAB. Oct 7, 2013 at 9:55

There is one fancy way specific for the quadratics. Suppose you want to minimize $\sum_n Q_n(x)^2$ where $Q_n$ are quadratic forms. Set the problem as $$\sum y_n^2\to \min, y_n-Q_n(x)=0$$ and consider both the objective and the conditions as quadratic forms of $z=(x,y)$. The Lagrange multiplier theorem tells that if you have a local minimum of $F(z)$ under the conditions $G_n(z)=0$, then we can find $q_n\in\mathbb R$ such that $$F_q(z)=F(z)+\sum_n q_n G_n(z)=(A(q)z,z)-2(B(q),z)+c$$ attains a global minimum at the same point. (note that $c$ does not depend on $q=(q_n)$!) Moreover, the global minimum in the original problem corresponds to $q$ such that $(A(q)^{-1}B(q),B(q))$ is minimal and both the domain of admissible $q$ ($A(q)\ge 0$) and the functional $q:\mapsto (A(q)^{-1}B(q),B(q))$ are convex, so the methods of convex optimization work. The catch is that the Lagrange problem degenerates when you have multiple solutions in the original problem, so you'll easily determine the value of the minimum of $F(z)$ but not immediately the value of $z$. That would require tracing the degenerate plane. However, this is quick enough and easy to program. Whether it has been done already in a standard package, I cannot tell. I suggested this scheme to some guy on MSE for a very particular case and he seemed to be happy with it: https://math.stackexchange.com/questions/468576/least-squares-fit-with-a-trick/

• Why doesn't this approach solve the general (NP-hard) quadratically constrained quadratic program (QCQP)? Aug 21, 2013 at 17:05
• Because "tracing the degenerate plane" is a nightmare for too many variables. This approach will detect that $\min\sum_n(x_n(1-x_n))^2=0$, but, of course, it won't list all vertices of the unit cube in subexponential time. However, with just 6-10 variables, it should be feasible. I'll have to do some actual programming and go over every little detail to make everything sure but the general idea may work in this particular setting (unless, of course, I miss some crucial point, which can also happen :)). Aug 21, 2013 at 17:28
• @Noah Stein The actual reason seems to be the special parabolic type of the constraints: for general quadratic constraints the Lagrange theorem gives much less. I think I should elaborate on this a bit not to make it sound too general. As to other details (degenerate plane tracing), just give me some time if you want either something working, or my confession that I was overly optimistic and my shameful retreat. Right now I see no principal obstacles but it'll take some effort to put a few things right :). Aug 21, 2013 at 23:20
• Interesting idea, but I would guess any improvement from the convexity is more than offset by the increase in functional complexity (from the inverse matrix). Still, I hope you try my example (or a similar one with higher constants so that there is no zero and the problem is truly a minimization problem) and prove me wrong...please keep us posted. Aug 22, 2013 at 2:29
• Part of what I said was gibberish (alas, I overlooked one essential difference between one constraint an many constraints) but I'm still fighting :) How many variables/equations do you actually have? I am looking at 6/7 case like the one you posted but if it is 200/300, that's a completely different story. Aug 25, 2013 at 22:57

In answer to the comment, it is possible to compute a Groebner basis for the given example in a very short time, and from here one finds several solutions, e.g.,

a= 0,

b= 0,

c= (2260*h)/783,

e= (863*h)/348,

i= 0,

f=( - sqrt(6016000*sqrt(852564769045) - 9064098136321))/(253200*sqrt(113)),

j=( - sqrt(1050472*sqrt(852564769045) - 741915647113))/(63300*sqrt(226)),

h=( - 783)/(100*sqrt(113)),

d=( - sqrt(852564769045) + 212440)/715290.

• I wasn't aware of the Groebner basis method. But, when I try this approach in Mathematica (version 6), the software is still thinking after hours, so it's been of no use to me. Are you using your own home-built software? I'll wait to close this question until others have had a chance to answer about minimization approaches, but I do believe your method is best in cases where a zero exists. Aug 22, 2013 at 2:37
• It should work with any CAS (I use reduce -reduce-algebra.com). The idea is, not to compute all solutions, which takes too long, but only some. For example, just set $a=0$. Aug 22, 2013 at 8:17
• Oh, I get it now and Mathematica then computes it quickly. Now that I think about things, I would prefer a minimal solution with j=0 and f=0 (or j=0 and i=0), but Mathematica can't compute this in reasonable time...do you see one? Sorry, I guess I underspecified by problem originally. Aug 22, 2013 at 19:50

What you have is an instance of a quadratically constrained quadratic program (QCQP). These problems are NP-hard in general (though it's possible your particular type of instance is not hard as fedja's answer seems to be suggesting). If you google that you'll find lots of work on approximation schemes, heuristics, etc.

Groebner basis methods have already been mentioned as an approach to exactly solving this kind of system of equations. They aren't a way to find the best least squares solution in a case where there isn't an exact solution.

These polynomial optimization problems are NP-Hard, so you shouldn't expect to be able to solve large instances in any reasonable amount of time. However, there are some approaches that have been developed that can be effective in solving (to some high accuracy) small to medium sized problems of this sort. I'd encourage you to look at the Gloptipoly package:

http://homepages.laas.fr/henrion/software/gloptipoly/