sagemath / sage

Main repository of SageMath
https://www.sagemath.org
Other
1.33k stars 453 forks source link

Action of MatrixGroup on a MPolynomialRing #4513

Open simon-king-jena opened 15 years ago

simon-king-jena commented 15 years ago

A group of n by n matrices over a field K acts on a polynomial ring with n variables over K. However, this is not implemented yet.

The following should work:

sage: M = Matrix(GF(3),[[1,2],[1,1]])
sage: N = Matrix(GF(3),[[2,2],[2,1]])
sage: G = MatrixGroup([M,N])
sage: m = G.0
sage: n = G.1
sage: R.<x,y> = GF(3)[]
sage: m*x
x + y
sage: x*m
x - y
sage: (n*m)*x == n*(m*x)
True
sage: x*(n*m) == (x*n)*m
True

On the other hand, we still want to have the usual action on vectors or matrices:

sage: x = vector([1,1])
sage: x*m
(2, 0)
sage: m*x
(0, 2)
sage: (n*m)*x == n*(m*x)
True
sage: x*(n*m) == (x*n)*m
True
sage: x = matrix([[1,2],[1,1]])
sage: x*m
[0 1]
[2 0]
sage: m*x
[0 1]
[2 0]
sage: (n*m)*x == n*(m*x)
True
sage: x*(n*m) == (x*n)*m
True

CC: @wdjoyner @malb

Component: commutative algebra

Keywords: matrix group, action, polynomial ring

Author: Simon King

Reviewer: David Loeffler, William Stein

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

simon-king-jena commented 15 years ago
comment:1

Sorry, in the above code I forgot to copy/paste the line

sage: R.<x,y> = GF(3)[]

Moreover, for the reasons above, the ticket should belong to commutative algebra, not just algebra (I was clicking on the wrong button).

simon-king-jena commented 15 years ago
comment:2

The patch matrixgroupCall.patch is without doctests, and I am not sure if it couldn't be done better. So, it needs more work.

For example, Martin mentioned the possibility (off list) to create a pyx file with a Cython function, and then the call method would use that function. It might pay off here, since there are tight loops and since the method has to deal with tuples or lists. So Cdefining might speed things up.

Opinions?

At least, the following now works:

sage: R.<x,y>=GF(3)[]
sage: M=Matrix(GF(3),[[1,2],[1,1]])
sage: G=MatrixGroup([M])
sage: g=G.0
sage: p=x*y^2
sage: g(p)
x^3 + x^2*y - x*y^2 - y^3
simon-king-jena commented 15 years ago

call method for MatrixGroupelement (this time with doc test) and left_matrix_action for MPolynomial

simon-king-jena commented 15 years ago
comment:4

Attachment: matrixgroupCall.patch.gz

Some changes:

Nicolas Thiéry suggested to implement the action of matrix groups by a method for polynomials (I call the method left_matrix_action), and, for convenience, also provide a __call__ method for MatrixGroupElement that refers to left_matrix_action.

This has several advantages:

I have two Cython concerns left:

  1. The innerst loop is

    prod([Im[k]**X[k] for k in xrange(n)])

    where k is c'defined as int. Should this better be done in a for-loop, rather then creating a list and calling prod?

  2. The variable X is of type polydict.ETuple, so I can not directly c'define X. One could do

    cdef tuple X
    for i from 0<i<l:
        X = tuple(Expo[i])

    But would this be faster?

simon-king-jena commented 15 years ago

Attachment: matrixgroupCallNew.patch.gz

Another version of left_matrix_action

simon-king-jena commented 15 years ago
comment:5

In matrixgroupCallNew.patch (to be applied after the first patch), I modified the method according to my above concerns. In the example from my original post, the average running time improves from ~240 microseconds to 164 microseconds, and in a larger example it improved from 6.5s to 5.4s

Nevertheless, I made two separate patches, so that the reviewer (if there is any...) can compare by him- or herself.

Cheers Simon

simon-king-jena commented 15 years ago
comment:6

One observation: Reverse the outer loop

        for i from l>i>=0:
            X = tuple(Expo[i])
            c = Coef[i]
            for k from 0<=k<n:
                if X[k]:
                    c *= Im[k]**X[k]
            q += c

It results in a further improvement of computation time. Is this coincidence? Or is it since summation of polynomials should better start with the smallest summands?

Anyway, I didn't change the patch yet.

simon-king-jena commented 15 years ago

Slight improvement; extended functionality

simon-king-jena commented 15 years ago
comment:7

Attachment: matrixgroupCallNew2.patch.gz

Replying to @simon-king-jena:

One observation: Reverse the outer loop

        for i from l>i>=0:
            X = tuple(Expo[i])
            c = Coef[i]
            for k from 0<=k<n:
                if X[k]:
                    c *= Im[k]**X[k]
            q += c

It results in a further improvement of computation time. Is this coincidence? Or is it since summation of polynomials should better start with the smallest summands?

I made a couple of tests, and there was a small but consistent improvement. So, in the third patch (to be applied after the other two) I did it in that way.

The left_matrix_action shall eventually be used for computing the Reynolds operator of a group action; moreover, the Reynolds operator should be applicable on a list of polynomials. Then, the function would repeatedly compute the image of the ring variables under the action of some group element. But then it would be better to compute that image only once and pass it to left_matrix_action. The new patch provides this functionality. Example (continuing the original example):

sage: L=[X.left_matrix_action(g) for X in R.gens()]
sage: p.left_matrix_action(L)
x^3 + x^2*y - x*y^2 - y^3
wdjoyner commented 15 years ago
comment:8

I did confirm that the patches apply cleanly, that

sage: M = Matrix(GF(3),[[1,2],[1,1]])
sage: G = MatrixGroup([M])
sage: g = G.0
sage: g

[1 2]
[1 1]
sage: P.<x,y> = PolynomialRing(GF(3),2)
sage: p = x*y^2
sage: g(p)
x^3 + x^2*y - x*y^2 - y^3
sage: (x+2*y)*(x+y)^2
x^3 + x^2*y - x*y^2 - y^3

works, and that the code seems well-documented.

However, I can't do testing on this machine (intrepid ubuntu) and some of the code is written in Cython, which I can't really 100% vouch for. Seems okay though and simple enough. Since speed was a topic of the comments above, my only question is that the segment

    396         for i from 0<=i<l: 
    397             X = Expo[i] 
    398             c = Coef[i] 
    399             q += c*prod([Im[k]**X[k] for k in xrange(n)]) 

could probably be rewritten as a one-line sum, which might (or might not) be faster.

Maybe Martin Albrecht could comment on the Cython code?

If Martin (for example) passes the Cython code, and the docstrings pass sage -testall, I would give it a positive review.

malb commented 15 years ago
comment:9
cdef list Im 
if isinstance(M,list): 
  Im = M 

shouldn't Im = M take care of the type checking anyway, so that a try-except block is sufficient? Also, I think maybe the type of p should be checked in the __call__ method and a friendly error message raised? Not sure though.

malb commented 15 years ago
comment:10

Cython code looks good (just read it).

williamstein commented 15 years ago
comment:11

REFEREE REPORT:

Check this out:

sage: R.<x,y> = GF(3)[]
sage: M=Matrix(GF(3),[[1,2],[1,1]])
sage: M2=Matrix(GF(3),[[1,2],[1,0]])
sage: G=MatrixGroup([M, M2])
sage: (G.0*G.1)(p)
-x^2*y + x*y^2 - y^3
sage: G.0(G.1(p))
x^2*y + x*y^2 + y^3

Oops, your left action -- which it better be if you use that notation -- ain't a left action! Oops

-- William

simon-king-jena commented 15 years ago
comment:12

Really Oops. Sorry.

I implemented it analogous to what is done in Singular. Perhaps I am mistaken in the sense that it is supposed to be a right action (which then would deserve another notation).

sage: (G.0(G.1((p))))
-x^2*y + x*y^2 - y^3
sage: (G.1*G.0)(p)
-x^2*y + x*y^2 - y^3

However, I think it doesn't matter what Singular does. I will look up the literature whether one really wants a right or left action, and how the left action is supposed to be.

simon-king-jena commented 15 years ago
comment:13

Replying to @simon-king-jena:

However, I think it doesn't matter what Singular does. I will look up the literature whether one really wants a right or left action...

I mean, something like "one has a left action on a variety, which gives rise to a right action on the coordinate ring". I have to sort it out.

If this is the case, then it should be better implemented in the __mul__ method of polynomials, isn't it? Such as

sage: p*G.1*G.0==p*(G.1*G.0)
True
williamstein commented 15 years ago
comment:14

Left actions should use __call__, right actions should use exponentiation.

Substitution is a right action. Substitution of the inverse is a left action.

simon-king-jena commented 14 years ago

Description changed:

--- 
+++ 
@@ -1,20 +1,48 @@
 A group of n by n matrices over a field K acts on a polynomial ring with n variables over K. However, this is not implemented yet.

-Off list, David Joyner suggested to implement it with a `__call__` method in `matrix_group_element.py`. Then, the following should work:
+The following should work:

-sage: M=Matrix(GF(3),[[1,2],[1,1]]) -sage: G=MatrixGroup([M]) -sage: g=G.0 -sage: p=xy^2 -sage: g(p) -x^3 + x^2y - xy^2 - y^3 -sage: _==(x+2y)(x+y)^2 +sage: M = Matrix(GF(3),[[1,2],[1,1]]) +sage: N = Matrix(GF(3),[[2,2],[2,1]]) +sage: G = MatrixGroup([M,N]) +sage: m = G.0 +sage: n = G.1 +sage: R.<x,y> = GF(3)[] +sage: mx +x + y +sage: xm +x - y +sage: (nm)x == n(mx) +True +sage: x(nm) == (xn)*m True


-Although it concerns `matrix_group_element.py`, I believe this ticket belongs to Commutative Algebra, for two reasons:
-1. An efficient implementation probably requires knowledge of the guts of MPolynomialElement.
-2. My long-term goal is to re-implement my algorithms for the computation of non-modular invariant rings. The current implementation is in the `finvar.lib` library of Singular -- the slow Singular interpreter sometimes is a bottle necks.
+On the other hand, we still want to have the usual action on vectors or matrices:

-One more general technical question: It is `matrix_group_element.py`, hence seems to be pure python. Is it possible to define an additional method in some `.pyx` file using Cython? I don't know if this would be reasonable to do here, but perhaps this could come in handy at some point...
+```
+sage: x = vector([1,1])
+sage: x*m
+(2, 0)
+sage: m*x
+(0, 2)
+sage: (n*m)*x == n*(m*x)
+True
+sage: x*(n*m) == (x*n)*m
+True
+```
+
+```
+sage: x = matrix([[1,2],[1,1]])
+sage: x*m
+[0 1]
+[2 0]
+sage: m*x
+[0 1]
+[2 0]
+sage: (n*m)*x == n*(m*x)
+True
+sage: x*(n*m) == (x*n)*m
+True
+```
simon-king-jena commented 14 years ago

Attachment: trac-4513_matrix_action_on_polynomials.patch.gz

Replaces the other patches

simon-king-jena commented 14 years ago
comment:16

The new patch solves the problems addressed in the new ticket description.

I worked on top of several other tickets, since I somehow cared about number fields. To be precise, I did

hg_sage.import_patch('https://github.com/sagemath/sage-prod/files/10649509/trac_9205-discrete_log.patch.gz')
hg_sage.import_patch('https://github.com/sagemath/sage-prod/files/10649510/trac_9205-doctest.patch.gz')
hg_sage.import_patch('https://github.com/sagemath/sage-prod/files/10649981/trac_9438_IntegerMod_log_vs_PARI.patch.gz')
hg_sage.import_patch('https://github.com/sagemath/sage-prod/files/10649919/trac_9423_gap_for_numberfields.patch.gz')
hg_sage.import_patch('https://github.com/sagemath/sage-prod/files/10649053/8909_gap2cyclotomic.patch.gz')
hg_sage.import_patch('https://github.com/sagemath/sage-prod/files/10649054/trac_8909_catch_exception.patch.gz')
hg_sage.import_patch('https://github.com/sagemath/sage-prod/files/10644116/trac_5618_gap_for_cyclotomic_fields.patch.gz')

before creating my pacth. But this shouldn't matter, I guess the patch applies cleanly to sage-4.4.

simon-king-jena commented 14 years ago

Author: Simon King

simon-king-jena commented 14 years ago
comment:18

I made Cc to malb and wdj, since it both concerns polynomials and groups.

d883e2d8-b666-438e-bb30-c4f09d46f55a commented 13 years ago
comment:19

Sorry to be a nag, but matrix_group_element.py doesn't pass the doctests:

----------------------------------------------------------------------
matrix_group_element.py
SCORE matrix_group_element.py: 87% (14 of 16)

Missing documentation:
     * is_MatrixGroupElement(x):
     * __invert__(self):

----------------------------------------------------------------------

This is the same doctest score that the old unpatched file gets too.

simon-king-jena commented 13 years ago
comment:20

Replying to @sagetrac-nharris:

Sorry to be a nag, but matrix_group_element.py doesn't pass the doctests:

Which one fails?

----------------------------------------------------------------------
matrix_group_element.py
SCORE matrix_group_element.py: 87% (14 of 16)

Missing documentation:
   * is_MatrixGroupElement(x):
   * __invert__(self):

----------------------------------------------------------------------

This is the same doctest score that the old unpatched file gets too.

Of course there is no change. I added code to one already existing method, extending its functionality, and added tests for the new functionality. But this patch is not about inversion of matrix group elements, and I think the patch is not supposed to add documentation to methods that are not in its scope.

So, unless a doc test fails because of my patch, the criticism about not raising the doc test coverage is invalid.

d883e2d8-b666-438e-bb30-c4f09d46f55a commented 13 years ago
comment:21

I'm sorry if I did something that I shouldn't have. I was just following this guideline for reviewing patches (found here http://www.sagemath.org/doc/developer/trac.html#section-review-patches):

Is it documented sufficiently, including both explanation and doctests? This is very important: all code in Sage must have doctests, so even if the patch is for code which did not have a doctest before, the new version must include one. In particular, all new code must be 100% doctested. Use the command sage -coverage to see the coverage percentage of .

simon-king-jena commented 13 years ago
comment:22

Replying to @sagetrac-nharris:

Is it documented sufficiently, including both explanation and doctests? This is very important: all code in Sage must have doctests, so even if the patch is for code which did not have a doctest before, the new version must include one. In particular, all new code must be 100% doctested. Use the command sage -coverage to see the coverage percentage of .

If you would read the patch, you would find:

  1. The patch adds one case to an existing method, namely _act_on_.

  2. The patch adds several doc tests to _act_on_, covering the new functionality.

  3. The original version of _act_on_ already had a doc test.

In particular, it is impossible to detect the doc test change without reading the patch, by just using sage -coverage: The coverage script detects whether _act_on_ has any test (this is the case with or without my patch), but it does not detect whether the patch extends existing documentation to cover a new case.

11d1fc49-71a1-44e1-869f-76be013245a0 commented 13 years ago
comment:23

Apply trac-4513_matrix_action_on_polynomials.patch

(for patchbot -- it's trying to apply all the patches at once)

11d1fc49-71a1-44e1-869f-76be013245a0 commented 13 years ago
comment:24

Apply trac-4513_matrix_action_on_polynomials.patch

(maybe it will work this time?)

11d1fc49-71a1-44e1-869f-76be013245a0 commented 13 years ago
comment:25
Apply trac-4513_matrix_action_on_polynomials.patch

(For some reason patchbot's not picking this up -- I apologise to all human beings reading this for the spam!)

11d1fc49-71a1-44e1-869f-76be013245a0 commented 13 years ago
comment:26

I've had a look at the patch, and I don't think you've addressed William's comment #14 from two years back. The following makes me extremely uneasy:

sage: G = GL(3, 7)
sage: R.<a, b> = GF(7)[]
sage: G.0 * a
[3*a   0   0]
[  0   a   0]
[  0   0   a]
sage: R.<a,b,c> = GF(7)[]
sage: G.0 * a
3*a

It looks like there's some pre-existing coercion mechanism which returns elements of the matrix space over R, and you're overriding it in one case with an alternative coercion that returns completely different answers; this violates a Sage coercion axiom (where there are multiple paths in the coercion diagram, all must give the same answer up to numerical precision issues). Moreover, if you look at the patchbot logs it seems to have found an example where the preexisting coercion gets picked up instead of the new one.

Sorry, that's a thumbs down from me.

David

11d1fc49-71a1-44e1-869f-76be013245a0 commented 13 years ago

Reviewer: David Loeffler, William Stein