ASKSAGE: Sage Q&A Forum - Individual question feedhttp://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Thu, 24 Apr 2014 09:12:58 -0500solving system of polynomial equations over reals using newton methodhttp://ask.sagemath.org/question/11359/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?Sat, 19 Apr 2014 12:12:09 -0500http://ask.sagemath.org/question/11359/solving-system-of-polynomial-equations-over-reals-using-newton-method/Answer by calc314 for <p>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:</p>
<p>This is the list of equations and variables:</p>
<pre><code>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]
</code></pre>
<p>If I ask Sage to solve this </p>
<pre><code>S = solve(Equations,Variables)
</code></pre>
<p>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 <code>[2,4,7,8]</code>. Is there any way that I can perform an algorithm like the Newton's method with the start point <code>[2,4,7,8]</code>, and find that real solution?</p>
http://ask.sagemath.org/question/11359/solving-system-of-polynomial-equations-over-reals-using-newton-method/?answer=16117#post-id-16117One 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
Sun, 20 Apr 2014 16:34:09 -0500http://ask.sagemath.org/question/11359/solving-system-of-polynomial-equations-over-reals-using-newton-method/?answer=16117#post-id-16117Comment by k1 for <p>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.</p>
<pre><code>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
</code></pre>
http://ask.sagemath.org/question/11359/solving-system-of-polynomial-equations-over-reals-using-newton-method/?comment=16194#post-id-16194A 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.Thu, 24 Apr 2014 09:11:54 -0500http://ask.sagemath.org/question/11359/solving-system-of-polynomial-equations-over-reals-using-newton-method/?comment=16194#post-id-16194Answer by calc314 for <p>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:</p>
<p>This is the list of equations and variables:</p>
<pre><code>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]
</code></pre>
<p>If I ask Sage to solve this </p>
<pre><code>S = solve(Equations,Variables)
</code></pre>
<p>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 <code>[2,4,7,8]</code>. Is there any way that I can perform an algorithm like the Newton's method with the start point <code>[2,4,7,8]</code>, and find that real solution?</p>
http://ask.sagemath.org/question/11359/solving-system-of-polynomial-equations-over-reals-using-newton-method/?answer=16116#post-id-16116Here 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]
return(add([q[i]^2 for i in range(0,4)]))
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
ans=fmin_tnc(f,[2,4,7,8],bounds=[(-10,10),(-10,10),(-10,10),(-10,10)],approx_grad=True)
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.
Sun, 20 Apr 2014 15:34:25 -0500http://ask.sagemath.org/question/11359/solving-system-of-polynomial-equations-over-reals-using-newton-method/?answer=16116#post-id-16116Comment by k1 for <p>Here are two approaches using <code>scipy.optimize</code> routines.</p>
<pre><code>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]
return(add([q[i]^2 for i in range(0,4)]))
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)
</code></pre>
<p>First, we try conjugate gradient minimization of the sum of squares.</p>
<pre><code>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)
</code></pre>
<p>Second, try truncated Newton conjugate gradient minimization of the sum of squares.</p>
<pre><code>#truncated Newton Conjugate Gradient minimization
from scipy.optimize import fmin_tnc
ans=fmin_tnc(f,[2,4,7,8],bounds=[(-10,10),(-10,10),(-10,10),(-10,10)],approx_grad=True)
print ans
#value of the sum of square
print f(ans[0])
#value of each element of the list
g(ans[0])
</code></pre>
<p>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.</p>
http://ask.sagemath.org/question/11359/solving-system-of-polynomial-equations-over-reals-using-newton-method/?comment=16193#post-id-16193I 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.Thu, 24 Apr 2014 09:12:58 -0500http://ask.sagemath.org/question/11359/solving-system-of-polynomial-equations-over-reals-using-newton-method/?comment=16193#post-id-16193