Opened 3 years ago

Closed 3 years ago

#19870 closed enhancement (fixed)

Cythonize matrix group elements

Reported by: tscrim Owned by: tscrim
Priority: major Milestone: sage-7.1
Component: cython Keywords: matrix group, elements
Cc: jdemeyer, vbraun Merged in:
Authors: Travis Scrimshaw Reviewers: Frédéric Chapoton
Report Upstream: N/A Work issues:
Branch: ddd203e (Commits) Commit: ddd203eb47d9d221d4c407a5cb9c09700f49a0f1
Dependencies: Stopgaps:

Description

They are used in a lot of places and should be cythonized. With my branch, I get ~10% speedup of my operations. Moreover, the class hierarchy for matrix group elements can and needs to be changed in order to do this. Specifically, I removed the MatrixGroupElement_base (the AffineGroupElement overrode all but 1 method) and GroupElementMixinLibGAP (only used by MatrixGroupElement_gap).

Change History (25)

comment:1 Changed 3 years ago by tscrim

  • Branch set to public/groups/cythonize_matrix_group_element-19870
  • Cc jdemeyer vbraun added
  • Commit set to 6358d539440f548b4fac34cfcaed78a4563f2868
  • Status changed from new to needs_review

While this new structure has some near-duplication by not using the calls to matrix(), it does away with a diamond problem and has fewer function calls overall by accessing the internal structure when necessary. Moreover, the original doctests for MatrixGroupElement_generic were not actually testing that class.


New commits:

989b6eaAdding changes from #19821.
6358d53Cythonized matrix_gp/group_element.py and simplified the class structure.

comment:2 Changed 3 years ago by chapoton

  • Status changed from needs_review to needs_work

many failing doctests, see patchbot

comment:3 Changed 3 years ago by git

  • Commit changed from 6358d539440f548b4fac34cfcaed78a4563f2868 to 11109267de4a9eee93b6a32126ff253e485f1910

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

65ce465Fixing modform_hecketriangle due to changes.
1110926Fixing last doctests.

comment:4 Changed 3 years ago by tscrim

  • Status changed from needs_work to needs_review

Fixed.

comment:5 Changed 3 years ago by vdelecroix

For MatrixGroupElement_generic the __invert__ is very bad since it actually does a copy. If you have a matrix over ZZ its inverse belongs to matrices over QQ.

sage: m = matrix(3, [6, -1, -22, -5, 1, 18, -2, -4, 17])
sage: m.inverse().parent()
Full MatrixSpace of 3 by 3 dense matrices over Rational Field

Hence it would be better to use change_ring before the wrapping.

On the other hand there exists _invert_unit(). Compare

sage: m = matrix(3, [6, -1, -22, -5, 1, 18, -2, -4, 17])
sage: %timeit m.inverse().change_ring(ZZ)
10000 loops, best of 3: 53.7 µs per loop
sage: %timeit m._invert_unit()
100000 loops, best of 3: 6.56 µs per loop

comment:6 Changed 3 years ago by git

  • Commit changed from 11109267de4a9eee93b6a32126ff253e485f1910 to 29f72539754e8efba97deb07a93d735af559a552

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

3efe9b6Merge branch 'develop' into public/groups/cythonize_matrix_group_element-19870
29f7253Special case for dense matrices over ZZ and making sure the inverse is in the base ring.

comment:7 Changed 3 years ago by tscrim

At some point we do have to make a copy of the underlying matrix, and we can't modify self._matrix. It turns out that _invert_unit is only applicable to dense matrices over \ZZ, but it deserves a special case since checking if the base ring (in python) is less than 100 ns on my machine. I also made sure that the result lives in the base ring it should (which covers the sparse \ZZ case).

comment:8 Changed 3 years ago by vdelecroix

Why not

try:
    M = M._invert_unit()
except AttributeError:
    M = ~M
    if M.base_ring() is not parent.base_ring():
        M = M.change_ring(parent.base_ring())

That would let the possibility for other rings/sparse matrices to implement an _invert_unit if appropriate. Though, if somebody changes the name of _invert_unit it will not be discovered. An alternative would be to implement a naive _invert_unit directly as a matrix method (basically following what you did).

comment:9 Changed 3 years ago by tscrim

The exception handling is slower when it has to occur:

def foo(M):
    try:                                    
        return M._invert_unit()
    except AttributeError:
        return ~M

def bar(M):
    if M.base_ring() is ZZ and M.is_dense():
        return M._invert_unit()
    else:                 
        return ~M
sage: M = random_matrix(ZZ, 3, 3, sparse=True)
sage: M.det()
-1
sage: %timeit foo(M)
The slowest run took 15.06 times longer than the fastest. This could mean that an intermediate result is being cached 
1000 loops, best of 3: 643 µs per loop
sage: %timeit bar(M)
1000 loops, best of 3: 639 µs per loop
sage: M = random_matrix(QQ, 3, 3)
sage: M.det()
1
sage: %timeit foo(M)
The slowest run took 20.62 times longer than the fastest. This could mean that an intermediate result is being cached 
100000 loops, best of 3: 10.4 µs per loop
sage: %timeit bar(M)
The slowest run took 17.06 times longer than the fastest. This could mean that an intermediate result is being cached 
100000 loops, best of 3: 7.5 µs per loop

However, it is faster when no exception is generated:

sage: M = random_matrix(ZZ, 3, 3)
sage: M.det()
1
sage: %timeit foo(M)
The slowest run took 16.19 times longer than the fastest. This could mean that an intermediate result is being cached 
100000 loops, best of 3: 4.39 µs per loop
sage: %timeit bar(M)
The slowest run took 11.94 times longer than the fastest. This could mean that an intermediate result is being cached 
100000 loops, best of 3: 4.61 µs per loop

Yet, I don't think that is enough of a speedup to justify the other cases (in particular, the \QQ example).

comment:10 Changed 3 years ago by tscrim

Actually, perhaps directly implementing a default _inverse_unit is best as it might be useful for matrices over other non-fields. Thoughts?

comment:11 Changed 3 years ago by vdelecroix

I like the idea of a _inverse_unit but I would actually be in favor of _inverse_unit_. Moreover, such method makes sense for any ring, not only matrices. We just run into the same kind of trouble while implementing any (multiplicative) subgroup of units. The trivial example being

sage: parent(~ZZ(1))
Rational Field

comment:12 Changed 3 years ago by tscrim

I agree that _inverse_unit_ is better since it is (more like) a special function. The quick solution would be to put a default for this in the Rings category or in Element which just does return self.parent()(~self). Or do you think we should do something a little more comprehensive on a separate ticket?

comment:13 follow-up: Changed 3 years ago by jdemeyer

  1. I think you should add a .pxd file such that other Cython modules can cimport this if needed.
  1. Use the standard copyright template.
  1. Remove __cmp__ = _cmp_ twice (it doesn't work and it's not needed).
  1. You don't actually use this: from sage.structure.element import have_same_parent.
  1. There is not much point in making matrix() a cpdef method: you can just directly access the _matrix attribute from Cython.

comment:14 Changed 3 years ago by git

  • Commit changed from 29f72539754e8efba97deb07a93d735af559a552 to 5333b0f2c69a99623c43c2f660d7441a9de23a6d

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

5333b0fImplementing reviewer comments.

comment:15 in reply to: ↑ 13 Changed 3 years ago by tscrim

  • Milestone changed from sage-7.0 to sage-7.1

Replying to jdemeyer:

  1. I think you should add a .pxd file such that other Cython modules can cimport this if needed.
  1. Use the standard copyright template.
  1. Remove __cmp__ = _cmp_ twice (it doesn't work and it's not needed).
  1. You don't actually use this: from sage.structure.element import have_same_parent.
  1. There is not much point in making matrix() a cpdef method: you can just directly access the _matrix attribute from Cython.

All done. I also made _matrix a public attribute since doing access that way was fastest way to access the defining matrix by my testing.

comment:16 Changed 3 years ago by chapoton

  • Reviewers set to Frédéric Chapoton
  • Status changed from needs_review to positive_review

ok, good enough for me. And patchbots are green.

In a later ticker, we should overload the is_one method.

comment:17 Changed 3 years ago by vbraun

  • Status changed from positive_review to needs_work

Merge conflict, try the next beta...

comment:18 Changed 3 years ago by git

  • Commit changed from 5333b0f2c69a99623c43c2f660d7441a9de23a6d to df18b1300c4f94802f3d8706129846046aa38562

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

df18b13Merge branch 'public/groups/cythonize_matrix_group_element-19870' of trac.sagemath.org:sage into public/groups/cythonize_matrix_group_element-19870

comment:19 Changed 3 years ago by tscrim

  • Status changed from needs_work to positive_review

Trivial rebase because git forgot that the file changed names (even though I had done it with git mv).

comment:20 follow-up: Changed 3 years ago by jdemeyer

  • Status changed from positive_review to needs_work

In the .pxd, you should not repeat methods coming from base classes. So the _cmp_ and _mul_ should not be there.

comment:21 Changed 3 years ago by git

  • Commit changed from df18b1300c4f94802f3d8706129846046aa38562 to ddd203eb47d9d221d4c407a5cb9c09700f49a0f1

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

ddd203eFixing group_element.pxd file by not redefining inherited methods.

comment:22 in reply to: ↑ 20 Changed 3 years ago by tscrim

  • Status changed from needs_work to needs_review

Replying to jdemeyer:

In the .pxd, you should not repeat methods coming from base classes. So the _cmp_ and _mul_ should not be there.

Done.

comment:23 Changed 3 years ago by jdemeyer

  • Status changed from needs_review to positive_review

comment:24 Changed 3 years ago by tscrim

Thank you.

comment:25 Changed 3 years ago by vbraun

  • Branch changed from public/groups/cythonize_matrix_group_element-19870 to ddd203eb47d9d221d4c407a5cb9c09700f49a0f1
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.