# solving system of polynomial equations over reals using newton method

I have a set of polynomial equations and I want to find one of its real solutions close to some point, and I need only one solution. Here is an example:

This is the list of equations and variables:

Equations = [x_0*x_1*x_2*x_3 - x_0*x_1 - x_0*x_2 - x_0*x_3 - x_1*x_2 - x_1*x_3 + 2*x_0 + 2*x_1 - 448, -x_0*x_1*x_2 - x_0*x_1*x_3 - x_0*x_2*x_3 - x_1*x_2*x_3 + 3*x_0 +
3*x_1 + 2*x_2 + 2*x_3 + 452, x_0*x_1 + x_0*x_2 + x_0*x_3 + x_1*x_2 + x_1*x_3 + x_2*x_3 - 159, -x_0 - x_1 - x_2 - x_3 + 21]

Variables = [x_0, x_1, x_2, x_3]


If I ask Sage to solve this

S = solve(Equations,Variables)


it returns a bunch of solutions. But in some cases it doesn't give me any real solutions. I can prove that the above set of equations has a real solution close to [2,4,7,8]. Is there any way that I can perform an algorithm like the Newton's method with the start point [2,4,7,8], and find that real solution?

edit retag close merge delete

Sort by » oldest newest most voted

Here are two approaches using scipy.optimize routines.

def f(x):
q=[x[0]*x[1]*x[2]*x[3] - x[0]*x[1] - x[0]*x[2] - x[0]*x[3] - x[1]*x[2] - x[1]*x[3] + 2*x[0] + 2*x[1] - 448, -x[0]*x[1]*x[2] - x[0]*x[1]*x[3] - x[0]*x[2]*x[3] - x[1]*x[2]*x[3] + 3*x[0] + 3*x[1] + 2*x[2] + 2*x[3] + 452, x[0]*x[1] + x[0]*x[2] + x[0]*x[3] + x[1]*x[2] + x[1]*x[3] + x[2]*x[3] - 159, -x[0] - x[1] - x[2] - x[3] + 21]
def g(x):
q=[x[0]*x[1]*x[2]*x[3] - x[0]*x[1] - x[0]*x[2] - x[0]*x[3] - x[1]*x[2] - x[1]*x[3] + 2*x[0] + 2*x[1] - 448, -x[0]*x[1]*x[2] - x[0]*x[1]*x[3] - x[0]*x[2]*x[3] - x[1]*x[2]*x[3] + 3*x[0] + 3*x[1] + 2*x[2] + 2*x[3] + 452, x[0]*x[1] + x[0]*x[2] + x[0]*x[3] + x[1]*x[2] + x[1]*x[3] + x[2]*x[3] - 159, -x[0] - x[1] - x[2] - x[3] + 21]
return(q)


First, we try conjugate gradient minimization of the sum of squares.

from scipy.optimize import fmin_cg
#conjugate gradient minimization of the sum of squares
ans=fmin_cg(f,[2,4,7,8])
print ans
#value of the sum of squares
print f(ans)
#values of each element of the list
print g(ans)


Second, try truncated Newton conjugate gradient minimization of the sum of squares.

#truncated Newton Conjugate Gradient minimization
from scipy.optimize import fmin_tnc
print ans
#value of the sum of square
print f(ans[0])
#value of each element of the list
g(ans[0])


I've seeded these with your suggested initial values. I'm not sure the results are really giving a solution to your equations, given the values for each equations at the suggested point of minimization. However, you might be able to modify the code above to get what you need.

more

I worked out a Newton's method base technique and got some solutions that works for me, but the code is really ugly, so I am waiting to see if there are any better solutions out here.

( 2014-04-24 09:12:58 -0600 )edit

One more possible approach is to use a Groebner basis. Unfortunately, after finding such a basis, I'm still getting non-real solutions. But, again, maybe you can modify this approach to find what you want.

R = PolynomialRing(QQ, 4, 'abcd', order='lp')
a,b,c,d=R.gens()
I=(a*b*c*d - a*b - a*c - b*c - a*d - b*d + 2*a + 2*b - 448, -a*b*c - a*b*d - a*c*d - b*c*d + 3*a + 3*b + 2*c + 2*d + 452, a*b + a*c + b*c + a*d + b*d + c*d - 159, -a - b - c - d + 21)*R
B = I.groebner_basis()
print B

more

A Gorebner basis is not going to work for my problem, since my actual problem is much bigger, and even on 10 equations and 10 variables the Groebner basis gets really ugly. Thanks for the suggestions though.

( 2014-04-24 09:11:54 -0600 )edit