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:
- Is there something wrong with my approach?
- If not, is there a way to avoid taking such a long time to calculate the dimension of the ideal?
- 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.
- 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
varietyin 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.
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,degrevlexoften works faster thanlex.Thanks for your comment. I always check the dimensionality before calling the variety function. The point is that if you calculate the dimension once, it is normally saved in the object. However, the way calling variety internally works now is that first a transformation of the ideal is performed and then the dimension is calculated again. This second dimension calculation is taking extremely long, and I think it can simply be sidestepped by doing it before the transformation as I believe the dimension can't have changed.
As per the degrevlex vs lex comment, do you know a good heuristic or test to decide which to use? Or alternatively if Sage already has such heuristics in the background?
You may like to play with a custom implementation of variety: https://ask.sagemath.org/question/555...
As for the term order, degrevlex seems to be used by default and that's not without a reason. The choice of term order is fixed (to default or given one) at the time of polynomial ring creation and it's not a subject to any heuristic.