ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Fri, 29 Apr 2022 14:54:42 +0200Calculating variety of ideal does not finishhttps://ask.sagemath.org/question/62210/calculating-variety-of-ideal-does-not-finish/Hello everyone,
I'm currently working with systems of multivariate polynomials and trying to find their roots. The way I'm doing this is quite straightforward. First I make a polynomial ring, then I form the ideal from my polynomials and subsequently call the ``variety`` function on this ideal. The problem is that when crossing a certain threshold of the number of variables in the polynomial, it does not seem to finish even when running for multiple hours. Consider the code below (please note I'm using sage directly in python so the syntax might be slightly different):
from sage.all import *
p0_0, p0_1, p1_0, p1_1, p2_0, p2_1, p2_2, p3_0, p3_1, p3_2 = var('p0_0, p0_1, p1_0, p1_1, p2_0, p2_1, p2_2, p3_0,
p3_1, p3_2')
polys = [-4*p1_0*p2_0*p3_0 - p1_1*p2_0*p3_0 + 2*p1_0*p2_1*p3_0 + 4*p1_1*p2_1*p3_0 + p1_0*p2_2*p3_0 -
2*p1_1*p2_2*p3_0 + 3*p1_0*p2_0*p3_1 - 2*p1_1*p2_0*p3_1 + 2*p1_0*p2_1*p3_1 - 5*p1_1*p2_1*p3_1 + 5*p1_0*p2_2*p3_1
+ 6*p1_1*p2_2*p3_1 - 3*p1_0*p2_0*p3_2 + 8*p1_1*p2_0*p3_2 + 7*p1_0*p2_1*p3_2 - 5*p1_1*p2_1*p3_2 +
4*p1_0*p2_2*p3_2 - 4*p1_1*p2_2*p3_2, p0_0 + p0_1 - 1, -2*p0_0*p2_0*p3_0 + 2*p0_1*p2_0*p3_0 - p0_0*p2_1*p3_0 +
8*p0_1*p2_1*p3_0 - 5*p0_0*p2_2*p3_0 - 2*p0_1*p2_2*p3_0 - 5*p0_0*p2_0*p3_1 + 7*p0_1*p2_0*p3_1 - 3*p0_0*p2_1*p3_1
+ 5*p0_1*p2_1*p3_1 - 7*p0_0*p2_2*p3_1 - 2*p0_1*p2_2*p3_1 - p0_0*p2_0*p3_2 + 2*p0_1*p2_0*p3_2 - 9*p0_0*p2_1*p3_2
- p0_1*p2_1*p3_2 - 7*p0_0*p2_2*p3_2 - p0_1*p2_2*p3_2, p1_0 + p1_1 - 1, 3*p0_0*p1_0*p3_0 + 4*p0_1*p1_0*p3_0 -
4*p0_0*p1_1*p3_0 - p0_1*p1_1*p3_0 - 2*p0_0*p1_0*p3_1 + 5*p0_1*p1_0*p3_1 - 2*p0_0*p1_1*p3_1 + p0_1*p1_1*p3_1 -
3*p0_0*p1_0*p3_2 + p0_1*p1_0*p3_2 + 9*p0_0*p1_1*p3_2 - 5*p0_1*p1_1*p3_2, 4*p0_0*p1_0*p3_0 + 5*p0_1*p1_0*p3_0 -
p0_0*p1_1*p3_0 + p0_1*p1_1*p3_0 - 3*p0_0*p1_0*p3_1 + 7*p0_1*p1_0*p3_1 - p0_0*p1_1*p3_1 + 5*p0_1*p1_1*p3_1 -
6*p0_0*p1_0*p3_2 + 4*p0_1*p1_0*p3_2 + 3*p0_0*p1_1*p3_2 - 7*p0_1*p1_1*p3_2, p2_0 + p2_1 + p2_2 - 1,
4*p0_0*p1_0*p2_0 - 4*p0_1*p1_0*p2_0 + p0_0*p1_1*p2_0 + p0_1*p1_1*p2_0 - 3*p0_0*p1_0*p2_1 + p0_1*p1_0*p2_1 -
3*p0_0*p1_1*p2_1 + 4*p0_1*p1_1*p2_1 + 9*p0_0*p1_0*p2_2 - 3*p0_1*p1_0*p2_2 - 3*p0_0*p1_1*p2_2 -
3*p0_1*p1_1*p2_2, -4*p0_0*p1_0*p2_0 - 3*p0_1*p1_0*p2_0 - 6*p0_0*p1_1*p2_0 - 6*p0_1*p1_1*p2_0 + 6*p0_0*p1_0*p2_1
- 3*p0_1*p1_0*p2_1 + p0_0*p1_1*p2_1 - 3*p0_1*p1_1*p2_1 + 3*p0_0*p1_0*p2_2 + 4*p0_1*p1_0*p2_2 + 5*p0_0*p1_1*p2_2
- p0_1*p1_1*p2_2, p3_0 + p3_1 + p3_2 - 1]
poly_ring = PolynomialRing(QQ, [p0_0, p0_1, p1_0, p1_1, p2_0, p2_1, p2_2, p3_0, p3_1, p3_2])
poly_ideal = poly_ring.ideal(*polys)
sols = poly_ideal.variety(RR)
print(sols)
I'm very much aware that the list of polynomials is very long, but for the purpose of this post it is not important what these polynomials actually represent.
Now, when investigating this issue, I have found that where it gets stuck is line 1061 in `` sage/rings/polynomial/multi_polynomial_ideal.py``. What it does is the following:
if I.dimension() != 0:
raise TypeError("dimension must be zero")
When commenting these two lines out, the code runs in a matter of seconds. My question now has multiple parts:
1. Is there something wrong with my approach?
2. If not, is there a way to avoid taking such a long time to calculate the dimension of the ideal?
3. Why does it calculate the dimension of the transformed ideal, and not the original ideal? Can these dimensions be different? If these dimensions can not be different, it would make way more sense to try and get the dimension from the beginning, as it might already be cached.
4. What are the dangers of commenting out the two lines above. If all it does is catch the case where the dimension is not zero, can't I simply wrap my call to ``variety`` in a try-catch statement and avoid the costly dimension check?
Thanks in advance!
ps. If this is more of an issue than a question, please let me know and I will open an actual issue.WillemRFri, 29 Apr 2022 14:54:42 +0200https://ask.sagemath.org/question/62210/Find solution to Polynomial Sequences without going through varietyhttps://ask.sagemath.org/question/47357/find-solution-to-polynomial-sequences-without-going-through-variety/Hi everyone,
I'm computing the Groebner basis of an ideal defined over the QQ ring. Once I have this Groebner basis, I would like to obtain a set of values that satisfy the equations in the Groebner basis. I know that the full set is going to be the variety of the ideal, but since this object is huge, I might not be interested in finding those values.
In Maple, after computing the Groebner basis, I'm able to call the solve() method on it and even set a maximum number of solutions I want to obtain. As a small example:
P.<x,y,z,t>=PolynomialRing(QQ,4)
I = P.ideal(x*(x-1), y*(y-1), z*(z-1), t - x*y*z)
gb = I.groebner_basis()
Here I could have called `I.variety()` or `gb.variety()` and obtained the same set of solutions:
sage: gbI.variety()
[{y: 0, z: 0, t: 0, x: 0},
{y: 0, z: 0, t: 0, x: 1},
{y: 1, z: 0, t: 0, x: 0},
{y: 1, z: 0, t: 0, x: 1},
{y: 0, z: 1, t: 0, x: 0},
{y: 0, z: 1, t: 0, x: 1},
{y: 1, z: 1, t: 0, x: 0},
{y: 1, z: 1, t: 1, x: 1}]
But I would like to know if it is possible to call a method like (I can do this in Maple):
solve(gb,[max_sol=2])
Such that I can obtain a subset of the variety instead of the whole set. My motivation is that the size of the initial system of polynomials that I have is considerably larger than this example, and finding the feasible solutions on the reduced Groebner basis is more manageable. I might also not be interested in all the elements in the variety. Finally, if I transform the Groebner basis in an ideal itself and try to compute the variety on that object
gbI = ideal(gb)
gbI.variety()
I find the following error
RuntimeError: error in Singular function call 'groebner':
int overflow in hilb 3
error occurred in or before standard.lib::stdhilb line 299: ` intvec hi = hilb( Id(1),1,W );`
expected intvec-expression. type 'help intvec;'
leaving standard.lib::stdhilb
leaving standard.lib::groebner
bernaldeSat, 03 Aug 2019 03:27:39 +0200https://ask.sagemath.org/question/47357/