# 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

edit retag close merge delete

The variety method can be turned into an iterator, but it is not clear whether it will be faster (it depends on where is the time-consuming part of the code). Could you please provide the code for the genuine larger example so that we can test ?

2

Can we find out what Maple's algorithm is? I'm also interested in this.

Sort by » oldest newest most voted

It is not currently possible to compute just a small subset of solutions, as far as I can tell. It would not be hard to improve the implementation of I.variety() to return the results lazily, though (for example using sage.misc.lazy_list). This could indeed be useful whenever the number of solutions is large.

However, the expensive part of the computation is the computation of a Gröbner basis which cannot be done lazily. I cannot reproduce the error you mentioned, but this step seems to be the one that fails. The implementation first computes a triangular decomposition which requires a Gröbner basis with respect to the lexicographical term ordering. It could help if you are able to compute a lex-ordered Gröbner basis beforehand, possibly using a different algorithm more suitable for your concrete example.

Also, I encourage you to take a look at the https://github.com/sagemath/sage/blob/8.8/src/sage/rings/polynomial/multi_polynomial_ideal.py#L2516-L2570 (source code) of this function to better understand what it does.

more

Thanks for your response! I understand the point that you are making. I was just wondering that, since I already have a Groebner basis computed, when I call the "solve" method on it, it shouldn't try to spend time recomputing the Groebner basis itself. Maybe we can consider having a lazy option for the variety computation which before everything queries the ideal if it is a Groebner basis first. In any case, for practical purposes, finding solutions to systems of polynomial equations more efficiently is always something we should aim for (and Maple is doing a great job in that front).

The results of Groebner basis computations are usually cached in Sage, so that the same expensive computation is not performed twice. In fact, the implementation of triangular_decomposition does also check whether the ideal generators form a Groebner basis, but if the term ordering is different from lex, the ordering needs to be changed to lex via transformed_basis, which potentially is an expensive operation as well.