Macaulay matrix for Sequence Multivariate Polynomials

[MercedesHaiech]

We implemented the Macaulay matrix for Sequence of multivariate polynomials

📝 Checklist

• The title is concise and informative.
• The description explains in detail what this PR is about.
• I have linked a relevant issue or discussion.
• I have created tests covering the changes.
• I have updated the documentation and checked the documentation preview.

[github-actions]

Documentation preview for this PR (built with commit 50c83b7; changes) is ready! 🎉
This preview will update shortly after each push to this PR.

[vneiger]

Thanks for this nice addition! I will review soon.

A first (too long) note on a non-programming-related point, but following our discussion concerning what reference to show for Macaulay matrices. The Macaulay 1902 paper that you refer to in your method's documentation may look like a sufficient reference by itself. Another reference, similar but maybe more to the point concerning the construction of the matrices, could be the chapter I of the book with SageMath reference [Mac1916] (F.S. Macaulay, The algebraic theory of modular systems, 1916, chapter I available here). As we discussed, picking a more recent presentation of these matrices could have some advantage, for example one that would also cover non-homogeneous polynomials... but finding "the right" reference is tricky (among others, because putting a 21th century reference for something already described more than 100 years before seems somewhat odd). The reference you introduced (Bardet & Faugère & Salvy, 2015) in fact refers directly to Macaulay's 1902 paper you listed: It is classical that the computation of Gröbner bases can be performed by linear algebra on a large matrix that has been described precisely by Macaulay (1902). A suggestion could be Lazard's 1983 paper "Gröbner bases, Gaussian elimination and resolution of systems of algebraic equations", which defines the Macaulay matrix for non-homogeneous polynomials in section 2; another suggestion (my preference, I think) could be to put no reference other than Macaulay 1902/1916.

[vneiger]

Depending on the decision for the suggestion about references, this might have to be removed for the moment (to avoid a reference that is currently not being referred to in Sage).

[vneiger]

Mentioning this to @robinkouba who might have some comments since he recently dealt with Macaulay matrices in SageMath.

[robinkouba]

I think it would be great to add an option for ordering the columns and rows.

For the columns, this is straightforward: one can simply provide a monomial order as input, possibly along with a choice between increasing or decreasing order.

For the rows, the situation is more complex since they are indexed by elements of the form (m,i), where m is a monomial and i indexes the polynomials in the input. This indexing corresponds to monomials in a free K[x1​,…,xn​]-module, so perhaps an ordering could also be provided as input, but I think sagemath currently does not have monomial orderings for modules. At least, providing the choice between position over term and term over position could be easy to implement (the current ordering seems to correspond to position over term).

[MercedesHaiech]

I think it would be great to add an option for ordering the columns and rows.

For the columns, this is straightforward: one can simply provide a monomial order as input, possibly along with a choice between increasing or decreasing order.

For the rows, the situation is more complex since they are indexed by elements of the form (m,i), where m is a monomial and i indexes the polynomials in the input. This indexing corresponds to monomials in a free K[x1​,…,xn​]-module, so perhaps an ordering could also be provided as input, but I think sagemath currently does not have monomial orderings for modules. At least, providing the choice between position over term and term over position could be easy to implement (the current ordering seems to correspond to position over term).

Thanks for the feedback. I could indeed add an option for the order.
For the rows I'll see if there's already something coded in sagemath for the module, if not I'll think about adding at least one sub-function to handle the case of the Macaulay matrix.

I was also hesitating to add an option to make a Macaulay matrix with the reduced equations modulo the field equations. This is often useful in cryptography on F2, for example.

[vneiger]

There might be an input restriction to indicate here (this set must consist of variables from the parent of the input polynomials). More importantly: it is not very clear what this set is impacting. Is this about variables for the monomials that are constructed to multiply the input polynomials? or is this also about monomials that index the columns and appear in the input polynomials?

[MercedesHaiech]

Thank you very much for all the suggestions. I modified the word 'system' in 'sequence' and specified that set_variables was used to multiply the polynomials of the starting system by this set.
I hope I got the push right.
I'll have a look at your other suggestions.
Thanks again!

[vneiger]

Here are some more comments with suggested changes or improvements. (I have not yet tried to compile the documentation to see the html form, but I think there will be some things to fix there as well.)

[vneiger]

Two comments:

• What is the expected use of remove_zero? At first I thought it would remove possible zero columns from the output matrix. Instead, it hides all columns corresponding to monomials that do not appear in input polynomials. But these hidden monomials may appear in multiples constructed for the matrix, so it seems that we are hiding potentially important information: is there a (somehow important) use case that was thought of when designing this?

• The doc example shows a usage of the set_variables option which is set to ['x']: I would suggest to restrict the possibilities to actual variables, e.g. forbidding strings/characters (or maybe better: do some mechanism like what is in the method remove_var, which would also ensure that the monomial order is preserved if possible). Otherwise, this yields the kind of confusing things that you can observe below. If you agree on this I can make a relevant change suggestion in the code.

sage: ring = PolynomialRing(GF(7), 3, 'x')
sage: x2,x1,x0 = ring.gens()
sage: L = Sequence([ring.random_element() for i in range(2)])
sage: L.macaulay_matrix(1, variables=["x0"])
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 4 0 4]
[0 0 0 0 0 0 0 1 0 0 2 4 0 0 0 0 4 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 1 0 4 4 5 0 1 0]
[0 0 1 0 0 4 0 4 0 0 5 0 0 1 0 0 0 0 0 0]
sage: L.macaulay_matrix(1, variables=[x0])
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 4 0 4]
[0 0 0 0 0 0 0 0 0 1 0 0 0 2 4 0 0 0 4 0]
[0 0 0 0 0 0 0 0 0 0 0 0 1 0 4 4 5 0 1 0]
[0 0 0 0 0 0 1 0 4 4 0 0 0 5 0 1 0 0 0 0]
sage: L.macaulay_matrix(1, variables=[x0]) == L.macaulay_matrix(1, variables=["x2"])
True
sage: L.macaulay_matrix(1, variables=[x0]) == L.macaulay_matrix(1, variables=["x0"])
False

[vneiger]

There were two redundant computations that impacted performance:

• The monomials_of_degree sets were recomputed typically two times, but potentially many times (depending on homogeneous or not, and on the degrees of the input polynomials). I suggest a fix based on precomputing the required sets once for all, see MercedesHaiech@9a95aef

• The multiplications "monomials * polynomial" are done many times at the step where the matrix is constructed: with how it was written, it was done once for each column monomial. I suggest to modify the matrix construction so that this multiplication is performed only once for each product. See MercedesHaiech@19e4e2d

Combining both improvements does speed things up:

sage: ring = PolynomialRing(GF(2), 4, 'x')
sage: L = Sequence([ring.random_element() for i in range(4)])
# BEFORE:
sage: %time M = L.macaulay_matrix(5)
CPU times: user 109 ms, sys: 0 ns, total: 109 ms
Wall time: 108 ms
sage: %time M = L.macaulay_matrix(10)
CPU times: user 4.5 s, sys: 21.9 ms, total: 4.52 s
Wall time: 4.52 s
sage: %time M = L.macaulay_matrix(13)
CPU times: user 23.8 s, sys: 191 ms, total: 24 s
Wall time: 24 s
# AFTER:
sage: %time M = L.macaulay_matrix(5)
CPU times: user 37 ms, sys: 0 ns, total: 37 ms
Wall time: 36.7 ms
sage: %time M = L.macaulay_matrix(10)
CPU times: user 1.22 s, sys: 0 ns, total: 1.22 s
Wall time: 1.22 s
sage: %time M = L.macaulay_matrix(13)
CPU times: user 5.96 s, sys: 3.81 ms, total: 5.97 s
Wall time: 5.97 s

[vneiger]

• What is the expected use of remove_zero? At first I thought it would remove possible zero columns from the output matrix. Instead, it hides all columns corresponding to monomials that do not appear in input polynomials. But these hidden monomials may appear in multiples constructed for the matrix, so it seems that we are hiding potentially important information: is there a (somehow important) use case that was thought of when designing this?

Please forget this comment, it seems I had not understood well what this option does, judging from the documentation. I'll make a suggestion for improving the description.

[vneiger]

Concerning performance, it makes a big difference to store in some dictionary the indices of column monomials, and then use this to populate the matrix (see MercedesHaiech@8496a77 ). This would allow the build matrices for much larger instances in a reasonable time, e.g. the timings below improve from ~1sec to ~40msec, and from ~5sec to ~90msec.

sage: ring = PolynomialRing(GF(2), 4, 'x')
sage: L = Sequence([ring.random_element() for i in range(4)])
# BEFORE:
sage: %time M = L.macaulay_matrix(5)
CPU times: user 37 ms, sys: 0 ns, total: 37 ms
Wall time: 36.7 ms
sage: %time M = L.macaulay_matrix(10)
CPU times: user 1.22 s, sys: 0 ns, total: 1.22 s
Wall time: 1.22 s
sage: %time M = L.macaulay_matrix(13)
CPU times: user 5.96 s, sys: 3.81 ms, total: 5.97 s
Wall time: 5.97 s
# AFTER:
sage: %time M = L.macaulay_matrix(5)
CPU times: user 10.8 ms, sys: 0 ns, total: 10.8 ms
Wall time: 10.4 ms
sage: %time M = L.macaulay_matrix(10)
CPU times: user 40.5 ms, sys: 883 µs, total: 41.4 ms
Wall time: 40.9 ms
sage: %time M = L.macaulay_matrix(13)
CPU times: user 91.2 ms, sys: 1.13 ms, total: 92.3 ms
Wall time: 91.7 ms

[vneiger]

Short summary: I have collected many suggestions in MercedesHaiech#1, to finalize this review/comments/suggestions, my plan is to give a try to the suggestion about ordering the rows and columns, and also make the option of specifying variables more robust. The third remaining point (if I don't forget one?) is about the field equations, which I will leave aside for the moment (it could also be a task for later in another PR).