add a method to test whether a polynomial is symmetric

[videlec]

As mentioned in this sage-devel thread it is desirable to have a simple check to test whether a given polynomial is symmetric (with respect to a given permutation group).

This ticket aims to implement a generic is_symmetric(self, group=None) on multivariate polynomials (where the default group is the full symmetric group).

We also update the check of symmetry in SymmetricFunctions so that construction from polynomials get faster. The computation time for

sage: n = 10
sage: f = x / (1-exp(-x))
sage: Sym = SymmetricFunctions(QQ)
sage: P = PolynomialRing(QQ, 'x', n)
sage: seq = prod(f.subs({f.default_variable(): var}) for var in P.gens())
sage: sage: inp = [(var, 0) for var in P.gens()]
sage: seq_taylor = seq.taylor(*inp, n // 2)
sage: Sym.e().from_polynomial(P(seq_taylor))

went down from a minute to half a second.

CC: @nbruin @tscrim

Component: algebra

Author: Vincent Delecroix

Branch/Commit: c53ac9b

Reviewer: Travis Scrimshaw, Markus Wageringel

Issue created by migration from https://trac.sagemath.org/ticket/29553

[videlec]

Commit: 41faa3e

[videlec]

Branch: u/vdelecroix/29553

[videlec]

Author: Vincent Delecroix

[videlec]

New commits:

41faa3e

is_symmetric for multivariate polynomials

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 41faa3e to 98db2ca

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

98db2ca

documentation fix

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

7c79429	fix doctest
4cabb57	speed up check in SymmetricFunctions

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 98db2ca to 4cabb57

[videlec]

Description changed:

--- 
+++ 
@@ -1,3 +1,17 @@
 As mentioned in [this sage-devel thread](https://groups.google.com/forum/#!topic/sage-devel/Z9YS7Hh146I) it is desirable to have a simple check to test whether a given polynomial is symmetric (with respect to a given permutation group).

 This ticket aims to implement a generic `is_symmetric(self, group=None)` on multivariate polynomials (where the default group is the full symmetric group).
+
+We also update the check of symmetry in `SymmetricFunctions` so that construction from polynomials get faster. The computation time for
+
+```
+sage: sage: n = 10
+sage: sage: f = x / (1-exp(-x))
+sage: sage: Sym = SymmetricFunctions(QQ)
+sage: sage: P = PolynomialRing(QQ, 'x', n)
+sage: sage: seq = prod(f.subs({f.default_variable(): var}) for var in P.gens())
+sage: sage: inp = [(var, 0) for var in P.gens()]
+sage: seq_taylor = seq.taylor(*inp, n // 2)
+sage: Sym.e().from_polynomial(P(seq_taylor))
+```
+went down from a minute to half a second.

[videlec]

Description changed:

--- 
+++ 
@@ -5,11 +5,11 @@
 We also update the check of symmetry in `SymmetricFunctions` so that construction from polynomials get faster. The computation time for

-sage: sage: n = 10 -sage: sage: f = x / (1-exp(-x)) -sage: sage: Sym = SymmetricFunctions(QQ) -sage: sage: P = PolynomialRing(QQ, 'x', n) -sage: sage: seq = prod(f.subs({f.default_variable(): var}) for var in P.gens()) +sage: n = 10 +sage: f = x / (1-exp(-x)) +sage: Sym = SymmetricFunctions(QQ) +sage: P = PolynomialRing(QQ, 'x', n) +sage: seq = prod(f.subs({f.default_variable(): var}) for var in P.gens()) sage: sage: inp = [(var, 0) for var in P.gens()] sage: seq_taylor = seq.taylor(*inp, n // 2) sage: Sym.e().from_polynomial(P(seq_taylor))

[mwageringel]

comment:6

It might be good to convert the polynomial to a dictionary first in order to speed up the looking up of coefficients, e.g.:

-        for e in self.exponents():
-            coeff = self[e]
-            for g in gens:
-                if self[e.permuted(g)] != coeff:
-                    return False
-        return True
+        pd = self.dict()
+        return all(pd[e.permuted(g)] == c for e, c in pd.items() for g in gens)

For the example at hand, this is noticeably faster:

sage: sage: n = 12
....: sage: f = x / (1-exp(-x))
....: sage: Sym = SymmetricFunctions(QQ)
....: sage: P = PolynomialRing(QQ, 'x', n)
....: sage: seq = prod(f.subs({f.default_variable(): var}) for var in P.gens())
....: sage: sage: inp = [(var, 0) for var in P.gens()]
....: sage: seq_taylor = seq.taylor(*inp, n // 2)
....: sage: g = (P(seq_taylor))
....: %time g.is_symmetric()
....:
CPU times: user 2.48 s, sys: 5.07 ms, total: 2.48 s
Wall time: 2.49 s                                       # before
CPU times: user 121 ms, sys: 2.11 ms, total: 123 ms
Wall time: 123 ms                                       # after
True

and even for small polynomials this seems to be at least as fast:

sage: R.<x,y,z> = QQ[]
....: f = (1 + x^2 + y^2 + z^2)^2
....: %timeit f.is_symmetric()
....:
10000 loops, best of 5: 177 µs per loop    # before
10000 loops, best of 5: 64.3 µs per loop   # after

Secondly, in permgroup_element.pyx, there are methods _act_on_list_on_position and _act_on_array_on_position. Do you think it would be more consistent to add a PermutationGroupElement._act_on_etuple_on_position, rather than ETuple.permuted?

Also, you will want to merge #29540.

[videlec]

comment:8

Replying to @mwageringel:

It might be good to convert the polynomial to a dictionary first in order to speed up the looking up of coefficients, e.g.:

It is sad this is the case. The polynomial datastructure is supposed to be fast in accessing coefficients... I believe it strongly depends on the base ring. But given the time difference I agree that it makes sense.

<SNIP

Secondly, in permgroup_element.pyx, there are methods _act_on_list_on_position and _act_on_array_on_position. Do you think it would be more consistent to add a PermutationGroupElement._act_on_etuple_on_position, rather than ETuple.permuted?

Sure. That will also be convenient for action of permutations on polynomials.

Also, you will want to merge #29540.

Indeed.

[tscrim]

comment:9

Replying to @videlec:

Replying to @mwageringel:

It might be good to convert the polynomial to a dictionary first in order to speed up the looking up of coefficients, e.g.:

It is sad this is the case. The polynomial datastructure is supposed to be fast in accessing coefficients... I believe it strongly depends on the base ring. But given the time difference I agree that it makes sense.

The MPolynomial_libsignular.__getitem__ code looks like this:

        m = p_ISet(1,r)
        i = 1
        for e in x:
            overflow_check(e, r)
            p_SetExp(m, i, int(e), r)
            i += 1
        p_Setm(m, r)

        while(p):
            if p_ExpVectorEqual(p, m, r) == 1:
                p_Delete(&m,r)
                return si2sa(p_GetCoeff(p, r), r, self._parent._base)
            p = pNext(p)

So it looks more like it is going through a list rather than a dict. I don't know how singular does this, but it looks like a very different data structure than the naïve implementation. IIRC, their data structure is solely to be efficient at computing Gröbner bases.

TL;DR Converting to a dict is definitely the best option with the current implementation.

As a more broader question, it might be worthwhile to consider reimplementing generic multivariate polynomials in Cython and only convert to (lib)singular when wanting a Gröbner basis.

[dimpase]

comment:10

How about generalising it to computing actions of linear group elements and of linear groups on polynomials?

It seems to be a bit restrictive to only have methods for fixed point computation, whereas it's only slightly less general, and you compute the action anyway!

[videlec]

comment:11

Replying to @dimpase:

How about generalising it to computing actions of linear group elements and of linear groups on polynomials?

It seems to be a bit restrictive to only have methods for fixed point computation, whereas it's only slightly less general, and you compute the action anyway!

Feel free to open a ticket for that. The code here is permuting exponents and do not touch the coefficients. The linear action might better be done via a proper action.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 4cabb57 to 4e3cffe

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

ccdd618	is_symmetric for multivariate polynomials
4e3cffe	speed up check in SymmetricFunctions

[videlec]

comment:13

Implementation changed to go via action of PermutationGroupElement on ETuple.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

1470131

don't call exponents + force conversion to SymmetricGroup(n)

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 4e3cffe to 1470131

[dimpase]

comment:15

Does it work for Laurent polynomials?

[dimpase]

comment:16

result._data = <int*> sig_malloc(sizeof(int)*result._nonzero*2)

looks like a memory leak, as I don't see a matching sig_free() call.

[tscrim]

comment:17

Replying to @dimpase:

result._data = <int*> sig_malloc(sizeof(int)*result._nonzero*2)

looks like a memory leak, as I don't see a matching sig_free() call.

Why would there be a sig_free()? result is the return value of the function, and it would be the job of the ETuple to handle freeing that in its deallocation I believe.

I also have a few comments:

-        - ``group`` (optional) - if set, test whether the polynomial is
-          symmetric with respect to the given permutation group.
+        - ``group`` (default: symmetric group) -- if set, test whether the
+          polynomial is invariant with respect to the given permutation group

This error message doesn't make sense to me raise ValueError("wrong argument 'group'").

-        for e, coeff in coeffs.items():
-            for g in gens:
-                if coeffs.get(g._act_on_etuple_on_position(e, True), zero) != coeff:
-                    return False
-        return True
+        return all(coeffs.get(g._act_on_etuple_on_position(e, True), zero) == coeff
+                   for e, coeff in coeffs.items() for g in gens)

In monomial.py:

-from sage.combinat.partition import Partition
+from sage.combinat.partition import Partition, _Partitions

-        out = self.sum_of_terms((Partition(e), c)
-                                for (e,c) in f.dict().items()
-                                if all(e[i+1] <= e[i] for i in range(len(e) - 1)))
+        R = self.base_ring()
+        out = self._from_dict({_Partitions(e): R(c)
+                                for (e,c) in f.dict().items()
+                                if all(e[i+1] <= e[i] for i in range(len(e)-1))})

[dimpase]

comment:18

sig_malloc is a thin wrapper around C malloc(). There is no garbage collection in Cython that would magically call sig_free() on what result._data points to - unless I seriously misunderstand something.

[tscrim]

comment:19

In ETuple:

    def __dealloc__(self):
        if self._data != <int*>0:
            sig_free(self._data)

So whenever result (which is an ETuple) is destroyed, then that memory is freed.

[dimpase]

comment:20

OK, I see, thanks for an explanation.

[videlec]

comment:21

Replying to @tscrim:

Replying to @dimpase:

result._data = <int*> sig_malloc(sizeof(int)*result._nonzero*2)

looks like a memory leak, as I don't see a matching sig_free() call.

Why would there be a sig_free()? result is the return value of the function, and it would be the job of the ETuple to handle freeing that in its deallocation I believe.

Indeed. That is exactly what is happening.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

0d72df3

some details

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 1470131 to 0d72df3

[fchapoton]

comment:23

beware that monomial.py has been modified recently, in #29540

[mwageringel]

comment:24

To me the naming "on position" suggests a certain way a permutation acts on lists, independently from the notion of left and right (although this action does satisfy the properties of a left action). So I think either the _on_position or the self_on_left should be removed, as these terms are at odds with each other.

I would probably drop the self_on_left as it is simple enough to call the method on the inverse permutation instead when necessary, but your current implementation does the opposite of what _act_on_list_on_position does. I had not considered this before - sorry. Since you mentioned the action on polynomials, I can see why the current implementation might be preferable. In that case, you could just rename the method to _act_on_etuple.

In any case, left and right seem to be mixed up currently:

sage: S = SymmetricGroup(6)
sage: p, q = S('(1,2,3,4,5,6)'), S('(1,2)(3,4)(5,6)')
sage: from sage.rings.polynomial.polydict import ETuple
sage: e = ETuple([10..15])
sage: right = lambda x, p: p._act_on_etuple_on_position(x, self_on_left=False)
sage: right(e, p * q) == right(right(e, p), q)  # should be True
False

[mwageringel]

comment:25

Replying to @tscrim:

As a more broader question, it might be worthwhile to consider reimplementing generic multivariate polynomials in Cython and only convert to (lib)singular when wanting a Gröbner basis.

This is exactly what the polydict implementation of polynomials does, no? You can construct a polynomial ring via PolynomialRing(..., implementation='generic') to use it, and when you want a Gröbner basis, only then it will convert to Singular (the conversion to libsingular does not seem to be supported).

On the other hand, I cannot think of many operations on polynomials for which random access to arbitrary coefficients is important.

[videlec]

comment:26

Replying to @mwageringel:

To me the naming "on position" suggests a certain way a permutation acts on lists, independently from the notion of left and right (although this action does satisfy the properties of a left action). So I think either the _on_position or the self_on_left should be removed, as these terms are at odds with each other.

I would probably drop the self_on_left as it is simple enough to call the method on the inverse permutation instead when necessary, but your current implementation does the opposite of what _act_on_list_on_position does. I had not considered this before - sorry. Since you mentioned the action on polynomials, I can see why the current implementation might be preferable. In that case, you could just rename the method to _act_on_etuple.

In any case, left and right seem to be mixed up currently:

sage: S = SymmetricGroup(6)
sage: p, q = S('(1,2,3,4,5,6)'), S('(1,2)(3,4)(5,6)')
sage: from sage.rings.polynomial.polydict import ETuple
sage: e = ETuple([10..15])
sage: right = lambda x, p: p._act_on_etuple_on_position(x, self_on_left=False)
sage: right(e, p * q) == right(right(e, p), q)  # should be True
False

I just copied what is in _act_on_list_on_position which claims to be a right action. Same behaviour

sage: S = SymmetricGroup(6)
sage: p, q = S('(1,2,3,4,5,6)'), S('(1,2)(3,4)(5,6)')
sage: right = lambda x,p: p._act_on_list_on_position(x)
sage: e = [10..15]
sage: right(e, p * q) == right(right(e, p), q)
False
sage: right(e, p * q) == right(right(e, q), p)
True

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 0d72df3 to 989716c

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

6f299fa	drop self_on_left
e119938	spring cleanup in sf/monomial (pep, code, doc details)
9570e17	trac 29540 still better code in monomial.py, other details there
f7381ea	trac 29540 fix import of _Partitions
989716c	Merge changes from ticket #29540

[mwageringel]

comment:28

Replying to @videlec:

I just copied what is in _act_on_list_on_position which claims to be a right action. Same behaviour

The action of permutations on matrices implemented in _act_on_ also has it the other way around, but for polynomials it works as expected:

sage: S = SymmetricGroup(6)
sage: p, q = S('(1,2,3,4,5,6)'), S('(1,2)(3,4)(5,6)')
sage: M = matrix.diagonal([1..6])
sage: (M * p) * q == M * (p * q)
False
sage: R = PolynomialRing(QQ, 'x', 6)
sage: (R.0 * p) * q == R.0 * (p * q)
True

Am I misunderstanding something about left and right actions in Sage?

[tscrim]

comment:29

Replying to @mwageringel:

Replying to @tscrim:

As a more broader question, it might be worthwhile to consider reimplementing generic multivariate polynomials in Cython and only convert to (lib)singular when wanting a Gröbner basis.

This is exactly what the polydict implementation of polynomials does, no? You can construct a polynomial ring via PolynomialRing(..., implementation='generic') to use it, and when you want a Gröbner basis, only then it will convert to Singular (the conversion to libsingular does not seem to be supported).

Sorry, I forgot generic is an overloaded word here and probably not the best word. What I meant was more universally shall we do this for rings that can be converted to (lib)singular, like those over Z.

On the other hand, I cannot think of many operations on polynomials for which random access to arbitrary coefficients is important.

Perhaps you're right. I can think of a number of things where you want to iterate over the pairs of coefficients and exponents. I don't think we have a method to do that (the default iterator is quite bad, getting the list of coefficients and list of monomials and zipping them together). I will open a ticket tomorrow or the next day to try and improve the iteration.

[tscrim]

comment:30

Replying to @mwageringel:

The action of permutations on matrices implemented in _act_on_ also has it the other way around, but for polynomials it works as expected:

sage: S = SymmetricGroup(6)
sage: p, q = S('(1,2,3,4,5,6)'), S('(1,2)(3,4)(5,6)')
sage: M = matrix.diagonal([1..6])
sage: (M * p) * q == M * (p * q)
False
sage: R = PolynomialRing(QQ, 'x', 6)
sage: (R.0 * p) * q == R.0 * (p * q)
True

Am I misunderstanding something about left and right actions in Sage?

Here is where it comes from I think:

sage: (M * p.matrix()) * q.matrix() == M * (p.matrix() * q.matrix())
True
sage: (M * p.matrix()) * q.matrix() == M * (p*q).matrix()
True
sage: p.matrix() * q.matrix() == (p*q).matrix()
True
sage: p.matrix() * M == M * p
True

(An oddity that needs fixing: p.matrix() works but p.to_matrix() is a NotImplemented.) Also Sage only knows it has a right action:

sage: (q * p) * M == q * (p * M)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<snip>
TypeError: unsupported operand parent(s) for *: 'Symmetric group of order 6! as a permutation group' and 'Full MatrixSpace of 6 by 6 sparse matrices over Integer Ring'

So I think the thing that needs to change in the perm group element _act_on_ is

             elif is_Matrix(left):
+               return left.with_permuted_columns(~self)
+        else:
+            if is_Matrix(left):
                return left.with_permuted_rows(self)

Addendum: Because of this:

sage: M * p.matrix()
[0 1 0 0 0 0]
[0 0 2 0 0 0]
[0 0 0 3 0 0]
[0 0 0 0 4 0]
[0 0 0 0 0 5]
[6 0 0 0 0 0]
sage: M.with_permuted_columns(~p)
[0 1 0 0 0 0]
[0 0 2 0 0 0]
[0 0 0 3 0 0]
[0 0 0 0 4 0]
[0 0 0 0 0 5]
[6 0 0 0 0 0]

[mwageringel]

Reviewer: Travis Scrimshaw, Markus Wageringel

[mwageringel]

comment:31

Replying to @tscrim:

Sorry, I forgot generic is an overloaded word here and probably not the best word. What I meant was more universally shall we do this for rings that can be converted to (lib)singular, like those over Z.

So you are suggesting to make the generic implementation the default? This implementation is not usually faster than the Singular backend I think, but it depends on the use case of course.

Perhaps you're right. I can think of a number of things where you want to iterate over the pairs of coefficients and exponents. I don't think we have a method to do that (the default iterator is quite bad, getting the list of coefficients and list of monomials and zipping them together). I will open a ticket tomorrow or the next day to try and improve the iteration.

That would indeed be nice to have, as it is such a common operation. I did not know the __iter__ method was implemented for polynomials and have always been zipping coefficients and monomials (or exponents), but it always felt odd to me, especially since this pattern fails for CombinatorialFreeModule elements if one does not pay attention to the sorting.

Regarding the action on matrices, we could handle that on a new ticket, as it alters existing behavior and is not really related to the aim of this ticket.

I am happy with this ticket as it is now. Except maybe there is one little detail from Travis' suggestion:

        out = self._from_dict({_Partitions.element_class(_Partitions, list(e)): R(c)
                               for (e,c) in f.dict().items()

The conversion R(c) should not be necessary here, as the coefficients should already be elements of the base ring.

[mwageringel]

comment:32

Replying to @tscrim:

(An oddity that needs fixing: p.matrix() works but p.to_matrix() is a NotImplemented.)

Do you know where this to_matrix comes from? I stumbled upon this a few days ago, but could not figure out where it was defined.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 989716c to c53ac9b

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

c53ac9b

don't convert since base ring is already correct

[tscrim]

comment:34

Replying to @mwageringel:

Replying to @tscrim:

Sorry, I forgot generic is an overloaded word here and probably not the best word. What I meant was more universally shall we do this for rings that can be converted to (lib)singular, like those over Z.

So you are suggesting to make the generic implementation the default? This implementation is not usually faster than the Singular backend I think, but it depends on the use case of course.

Yes, although I am not sure if this will be a good option. However, IIRC things like multiplying polynomials is really slow and could use another library to speed that up.

Perhaps you're right. I can think of a number of things where you want to iterate over the pairs of coefficients and exponents. I don't think we have a method to do that (the default iterator is quite bad, getting the list of coefficients and list of monomials and zipping them together). I will open a ticket tomorrow or the next day to try and improve the iteration.

That would indeed be nice to have, as it is such a common operation. I did not know the __iter__ method was implemented for polynomials and have always been zipping coefficients and monomials (or exponents), but it always felt odd to me, especially since this pattern fails for CombinatorialFreeModule elements if one does not pay attention to the sorting.

This is now #29595.

Regarding the action on matrices, we could handle that on a new ticket, as it alters existing behavior and is not really related to the aim of this ticket.

I agree that it should be a separate ticket.

I am happy with this ticket as it is now. Except maybe there is one little detail from Travis' suggestion:

        out = self._from_dict({_Partitions.element_class(_Partitions, list(e)): R(c)
                               for (e,c) in f.dict().items()

The conversion R(c) should not be necessary here, as the coefficients should already be elements of the base ring.

If we add that conversion, then we can remove the

assert self.base_ring() == f.base_ring()

which I think is not a good thing to enforce.

[tscrim]

comment:35

Replying to @mwageringel:

Replying to @tscrim:

(An oddity that needs fixing: p.matrix() works but p.to_matrix() is a NotImplemented.)

Do you know where this to_matrix comes from? I stumbled upon this a few days ago, but could not figure out where it was defined.

It comes from the finite Complex reflection group category. This can be easily fixed with an alias in the (finite) Coxeter group category (using the method canonical_matrix) and/or for the specific implementation of permutation groups.

[videlec]

comment:36

Replying to @tscrim:

Replying to @mwageringel:

Replying to @tscrim: I am happy with this ticket as it is now. Except maybe there is one little detail from Travis' suggestion:

        out = self._from_dict({_Partitions.element_class(_Partitions, list(e)): R(c)
                               for (e,c) in f.dict().items()

The conversion R(c) should not be necessary here, as the coefficients should already be elements of the base ring.

If we add that conversion, then we can remove the

assert self.base_ring() == f.base_ring()

which I think is not a good thing to enforce.

Indeed, I found this line a bit weird. Though I did not touch since it was beyond the scope of the ticket.

[mwageringel]

comment:37

Ok, it seems this is ready to be merged then? Let me set this ticket to positive, but please undo if you disagree.

[vbraun]

Changed branch from u/vdelecroix/29553 to c53ac9b