Make .variety() work for ideals of any dimension over finite fields

[maxale]

Problem Description

Clearly .variety() cannot produce a result for an ideal of positive dimension over an infinite field since such a variety is infinite. However, in the case of finite fields, the variety is always finite (as a subset of the cartesian power of the underlying field) irrespectively of the ideal dimension.

This RFE asks for an extention of the .variety() functionality to the ideals of positive dimension over finite fields.

Proposed Solution

When we have an ideal of positive dimension d generated by polynomials from a set S over a finite field F, we can fix values of some variables (with a smart choice of variables, their number would be d) in all possible ways to reduce the ideal dimension to 0 and enable application of the current functionality of .variety().

Fixing a value x0 for the variable x can be done by simply adding the polynomial x - x0 to the set S. An example is given at https://ask.sagemath.org/question/81379?answer=81382#post-id-81382

Overall, this calls for a recursive implementation of .variety(), which in the case of dimension > 0, would find a suitable variable x , iterate its possible values x0 over the elements of F, and call itself on the ideal generated by S and x - x0, combining all the varieties from those calls into the result.

For the choice of x, I think at very least it should (i) be present in some polynomial in S; and (ii) not be the sole variable of any polynomial in S (or, even better, S here be replaced by the Groebner basis of the ideal if it's computed).

Alternatives Considered

Alternatively to adding polynomial x - x0, it can make substitution x = x0 in all polynomials from S, which will effectively eliminate x, and thus reduce the number of variables and the ideal dimension (with a suitable choice of x). This may lead to a better performance than working with all the original variables, and may help in selecting right variables for elimination later on.

Additional Information

No response

Is there an existing issue for this?

• I have searched the existing issues for a bug report that matches the one I want to file, without success.

[vneiger]

If the ideal has dimension 1 or more, in many cases (at least when the field is not small), the variety will have very large cardinality (at least, larger than the field itself). So, this would mean providing a method which, for seemingly simple instances, will consume an unreasonable amount of time and memory. Could you please clarify what would be a use case that would justify to add such a functionality despite this "danger"? And in these use cases are there no other algebraic representations of the set of solutions that may be much more concise to store in memory and that would somehow give access to the sought solutions (lexicographic bases, ...) ?

Also, another approach could be to add the so-called "field equations" to your input polynomials. Is this for some reason not an acceptable solution to the problem (that would not require modifying the existing method)? To make this more clear consider the code below based on the example at the given link. When I run this, I get the variety about 10 times faster (for GF(2^4)) or 100 times faster (for GF(13^3)) than when I use the suggested solution; however note that for other instances this variant may be slower (this depends a lot on the input polynomials and on the field cardinality).

# Define a finite field and a bivariate polynomial ring on this field:
q = 13**3
Fq.<z> = GF(q)
R.<x, y> = Fq[]

# Define a bunch of bivariate polynomials:
polynomials = [
        x^15 + (z^2 + z)*x^12*y + x^9*y^2 + (z^2 + z + 1)*x^3*y^4 + y^5,
        x^28*y + x^25*y^2 + x^19*y^4 + (z^2 + z + 1)*x^16*y^5 + x^7*y^8 + (z^2 + z + 1)*x*y^10,
        x^48*y^5 + x^36*y^9 + x^33*y^10 + x^12*y^17 + x^9*y^18 + x^3*y^20,
        x^64*y^21 + x^16*y^37 + x^4*y^41 + x*y^42
]

import time

tt = time.time()
variety1 = sum( (R.ideal(polynomials+[x-x0]).variety() for x0 in Fq), [] )
tt = time.time() - tt
print(tt)

tt = time.time()
variety2 = R.ideal(polynomials + [x**q-x, y**q-y]).variety()
tt = time.time() - tt
print(tt)

print(sorted(variety1, key=str) == sorted(variety2, key=str))

[maxale]

@vneiger: I agree that there may be a "danger" of a solution being prohibitively large and/or requiring enormous resources to compute it. However, similar "dangers" exist here and there, starting with a potential impossibility of computing Groebner bases in reasonable time. Still, we want to solve such problems in the cases that are within our reach. Same thing applies to the current RFE - we won't be able to compute the variety in all the cases, but sometimes we can and it'd be nice to have a suitable machinery that will either get us a solution, or will provide us with an evidence that getting such a solution may be out of our reach.
Talking specifically about the memory requirement, we may want to alleviate it by providing a generator for the variety elements that will produce them one at a time but not store them altogether in memory. Also, I suspect that in some applications people may not really want to get all the solutions but just one / a few particular ones.

As for a particular approach to the problem, I do not have a strong preference for either one, but just described a couple of them. I'm fine if any other approach(es) or their combination is implemented.

[vneiger]

@vneiger: I agree that there may be a "danger" of a solution being prohibitively large and/or requiring enormous resources to compute it. However, similar "dangers" exist here and there, starting with a potential impossibility of computing Groebner bases in reasonable time.

Pushing this comparison further we could say that Gaussian elimination becomes infeasible for extremely large matrices. But for Gaussian elimination and for Gröbner bases, there are known complexity bounds (especially for zero-dimensional systems), which do not depend on the size of the base field. And these objects provide algebraic representations (basis of vector space or of ideal) of solutions which precisely allow one to represent all solutions without having to enumerate them.

Computing a Gröbner basis will surely be feasible on small examples (few polynomials, few variables, small degree). This is not the case (unless the field is small) of listing all solutions when there are more solutions than the cardinality of the field.

Still, we want to solve such problems in the cases that are within our reach. Same thing applies to the current RFE - we won't be able to compute the variety in all the cases, but sometimes we can and it'd be nice to have a suitable machinery that will either get us a solution, or will provide us with an evidence that getting such a solution may be out of our reach. Talking specifically about the memory requirement, we may want to alleviate it by providing a generator for the variety elements that will produce them one at a time but not store them altogether in memory.

Yes, such a generator would be a nice way around some of the issues. Yet this generator itself has to be reasonably efficient and memory-friendly, and computed in a reasonable time. By efficient I mean that a too naive version could sometimes give the next solution quite fast and some other times after a long computation time, which does not seem ideal. Did you have any approach in mind for the computation of this generator?

Also, I suspect that in some applications people may not really want to get all the solutions but just one / a few particular ones.

A few (or at least one) specific examples would help to understand the aim of this functionality. This is why I asked in my first message "Could you please clarify what would be a use case that would justify to add such a functionality". Without such understanding of relevance, as far as I'm concerned I will probably not spend time incorporating such a method into SageMath (maybe others, including yourself, will, and then don't hesitate to put me as a reviewer).

As for a particular approach to the problem, I do not have a strong preference for either one, but just described a couple of them. I'm fine if any other approach(es) or their combination is implemented.

Another related question, for the "specialize one variable" technique: how do you make this an algorithm? i.e. how do you choose how many variables to specialize, and how do you select which ones to specialize?

[DannyFeldman]

Hello,
I am new to SageMath and hope this is the right place for this message.

A related task, that might also be relevant to infinite varieties in real/complex fields, is to return any point in the variey when there are infinitely many such points. Such as FindInstance() in Mathematica
https://reference.wolfram.com/language/ref/FindInstance.html.en

Students in my group implemented this in SageMath and we interested in contributing the code to SageMath. Not sure yet how to do it yet.

[maxale]

@vneiger: I do not have specific examples, besides a couple of questions I've seen at ask.sagemath.org. As for the algorithm, I did not put much thought into it, but maybe any variable present in the Groebner basis (computed in the course of determining the dimension) would be ok, or maybe a random choice of a variable would do a job, etc.

@DannyFeldman: That's a great suggestion, but I think it should go into its own issue - could you please create one? As for contributing the code SageMath, please refer to SageMath's developer guide.

[DannyFeldman]

@vneiger, Done! Thanks for the encouragement.
#39682

We are about to integrate the solution to SageMath, but currently have compilations errors when compiling SageMath. Any hint for which forum takes care of these technical issues? I could not find an explicit one in the developer's guide.

[maxale]

@DannyFeldman: Thanks for creating the issue. Technical questions can be asked at https://ask.sagemath.org/ or sage-support google group https://groups.google.com/g/sage-support

[vneiger]

@vneiger: I do not have specific examples, besides a couple of questions I've seen at ask.sagemath.org. As for the algorithm, I did not put much thought into it, but maybe any variable present in the Groebner basis (computed in the course of determining the dimension) would be ok, or maybe a random choice of a variable would do a job, etc.

@maxale Ok, then let's wait and see if anyone is interested in actually crafting an algorithm (and SageMath implementation) for this.

@DannyFeldman Interesting variant of this issue, thanks! (most of my concerns expressed above about the output size and computing time should not apply to this variant)

[DannyFeldman]

@vneiger, @maxale, thank you for the valuable feedback. I just want to emphasize few points:

• Our code in Sagemath is ready. My group already wrote the pseudo code and then implemented. We now wish to integrate to the github of SageMath, but never done it before and new in this business. I hope it will not be too hard.
• The function is available in commercial software, including Maple and Mathematica, with many example usage on the web. So I guess there is a proof of demand.
• The algorithm we implemented (from the classic papers that I mentioned) can also be found in the tutorials of Mathematica here:
  https://reference.wolfram.com/language/tutorial/ComplexPolynomialSystems.html