# 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?

ps. If this is more of an issue than a question, please let me know and I will open an actual issue.

edit retag close merge delete

Just a few comments: First, .variety() is defined only for zero-dimensional ideals and so computing the dimension is somewhat unavoidable here. Second, computing both dimension and variety is based on computing Groebner basis, which is quite computationally expensive and is likely the cause of slowdown here. Third, sometimes you may improve computational performance by using a different term order - in particular, degrevlex often works faster than lex.