New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

matrices and groups: Smith Normal Form and an infinity test for fp groups #12705

Merged
merged 10 commits into from Jun 19, 2017

Conversation

Projects
None yet
5 participants
@valglad
Contributor

valglad commented Jun 2, 2017

This PR implements a function for finding a Smith Normal Form of a matrix over a principal ideal domain. It works by creating a diagonal matrix out of the tuple of invariant factors returned by smith_normal_factors and appending zero rows/columns to preserve the matrix shape.

smith_normal_factors can only work if the user manually adds an attribute "ring" to the matrix which should be equal to the domain over which the entries of the matrix are to be considered.
( @jksuom You mentioned I might want to use RawMatrix for this. I didn't because I didn't encounter any problems with the regular Matrix. Is it going to run into trouble in certain special cases?)

smith_normal_factors can be used to test if a given FpGroup is infinite by considering its abelianisation. The relators of the abelianisation can be written as an integer matrix, and if smith_normal_factors of that matrix contains 0, the abelianisation is infinite and hence the whole group is. This test is added as a check in the .order() method of the FpGroup class.

I noticed that SymPy's Integers can't be handled by the default integer ring which is PythonIntegerRing - it's meant only for ints. I've changed it to work with Integers as well, but maybe a SympyIntegerRing should be implemented instead?

Additionally, is_PID attribute for domains is introduced and has_Field and has_Ring attributes are replaced with is_Field and is_Ring.

@jksuom

This comment has been minimized.

Show comment
Hide comment
@jksuom

jksuom Jun 3, 2017

Member

I've changed it to work with Integers as well, but maybe a SympyIntegerRing should be implemented instead?

There are many different implementations of gcdex depending of the domain. There is one in polytools for SymPy's objects like Integers and polynomial expressions (yes, integers and polynomials are currently handled by the same function, which may not be the best solution), and the two (non-SymPy) integer classes PythonInteger and GMPYInteger each have their own methods. There should be no need to modify these methods.

Member

jksuom commented Jun 3, 2017

I've changed it to work with Integers as well, but maybe a SympyIntegerRing should be implemented instead?

There are many different implementations of gcdex depending of the domain. There is one in polytools for SymPy's objects like Integers and polynomial expressions (yes, integers and polynomials are currently handled by the same function, which may not be the best solution), and the two (non-SymPy) integer classes PythonInteger and GMPYInteger each have their own methods. There should be no need to modify these methods.

@valglad

This comment has been minimized.

Show comment
Hide comment
@valglad

valglad Jun 3, 2017

Contributor

There are many different implementations of gcdex depending of the domain. There is one in polytools for SymPy's objects like Integers

That means I can't use something general like domain.gcdex() in the code for the smith normal form. Should I have a variable function _gcdex instead that is set at the beginning to gcdex from polytools if the given domain is ZZ and the entries are Integers (and to domain.gcdex otherwise)? Alternatively, if we are willing to restrict our attention to only Integers and polynomials, I could just use gcdex from polytools for all cases.

Contributor

valglad commented Jun 3, 2017

There are many different implementations of gcdex depending of the domain. There is one in polytools for SymPy's objects like Integers

That means I can't use something general like domain.gcdex() in the code for the smith normal form. Should I have a variable function _gcdex instead that is set at the beginning to gcdex from polytools if the given domain is ZZ and the entries are Integers (and to domain.gcdex otherwise)? Alternatively, if we are willing to restrict our attention to only Integers and polynomials, I could just use gcdex from polytools for all cases.

@jksuom

This comment has been minimized.

Show comment
Hide comment
@jksuom

jksuom Jun 3, 2017

Member

That means I can't use something general like domain.gcdex() in the code for the smith normal form.

I think you can. The important domains seem to have gcdex(). There should be no need to use the polytools version when the matrix entries are domain elements.

Member

jksuom commented Jun 3, 2017

That means I can't use something general like domain.gcdex() in the code for the smith normal form.

I think you can. The important domains seem to have gcdex(). There should be no need to use the polytools version when the matrix entries are domain elements.

Functions returning normal forms of matrices
'''
def smith_normal_form(m, domain = None):

This comment has been minimized.

@kshitij10496

kshitij10496 Jun 3, 2017

Member

Can you add a docstring for this function along with information regarding type of inputs and returned values?
Few examples to go along with this would be ideal.

@kshitij10496

kshitij10496 Jun 3, 2017

Member

Can you add a docstring for this function along with information regarding type of inputs and returned values?
Few examples to go along with this would be ideal.

This comment has been minimized.

@jksuom

jksuom Jun 8, 2017

Member

There should also be some reference though there are not many apart from Wikipedia. That will contain further references to the theoretical background (Lang's Algebra, 3rd ed.) and applications to modules over PIDs (invariant factors, elementary divisors).

@jksuom

jksuom Jun 8, 2017

Member

There should also be some reference though there are not many apart from Wikipedia. That will contain further references to the theoretical background (Lang's Algebra, 3rd ed.) and applications to modules over PIDs (invariant factors, elementary divisors).

ind = [j for j in range(m.cols) if m[0,j] != 0]
if ind:
m = m.permute_cols([[0, ind[0]]])

This comment has been minimized.

@jksuom

jksuom Jun 4, 2017

Member

Else both the first column and the first row consist of zeros and we have to scan the whole matrix.

@jksuom

jksuom Jun 4, 2017

Member

Else both the first column and the first row consist of zeros and we have to scan the whole matrix.

This comment has been minimized.

@valglad

valglad Jun 5, 2017

Contributor

If both the first column and row are zero, the function will just add 0 to the invariants list which it should do in this case anyway.

@valglad

valglad Jun 5, 2017

Contributor

If both the first column and row are zero, the function will just add 0 to the invariants list which it should do in this case anyway.

This comment has been minimized.

@jksuom

jksuom Jun 6, 2017

Member

I think we should find the invariants in (divisibility) order starting with the gcd of all matrix entries and leaving the zeros, if any, to the end. Otherwise expressed, we should use the inverted inclusion order of the ideals generated by the invariants, first the largest one, generated by all matrix entries, and the smallest ideals, {0} in particular, coming last.

This means that we should add 0 to the invariant list only when there are no nonzero entries left. Therefore all entries should be scanned at this point.

@jksuom

jksuom Jun 6, 2017

Member

I think we should find the invariants in (divisibility) order starting with the gcd of all matrix entries and leaving the zeros, if any, to the end. Otherwise expressed, we should use the inverted inclusion order of the ideals generated by the invariants, first the largest one, generated by all matrix entries, and the smallest ideals, {0} in particular, coming last.

This means that we should add 0 to the invariant list only when there are no nonzero entries left. Therefore all entries should be scanned at this point.

"The matrix entries must be over a principal ideal domain")
else:
domain = m.ring

This comment has been minimized.

@jksuom

jksuom Jun 6, 2017

Member

We should check here that the matrix is not empty. If it is (i.e., len(m) == 0), then the empty list is returned.

@jksuom

jksuom Jun 6, 2017

Member

We should check here that the matrix is not empty. If it is (i.e., len(m) == 0), then the empty list is returned.

@debugger22

This comment has been minimized.

Show comment
Hide comment
@debugger22

debugger22 Jun 7, 2017

Member

Few caveats:

  1. has_Ring and has_Field, publicly accessible fields, are being changed without deprecation. IMO the formers make more sense.
  2. setattr isn't a very good way to set something on an instance from a user point of view. Initializing Matrix.ring with None and then later setting the value would be a better API.
Member

debugger22 commented Jun 7, 2017

Few caveats:

  1. has_Ring and has_Field, publicly accessible fields, are being changed without deprecation. IMO the formers make more sense.
  2. setattr isn't a very good way to set something on an instance from a user point of view. Initializing Matrix.ring with None and then later setting the value would be a better API.
@debugger22

This comment has been minimized.

Show comment
Hide comment
@debugger22

debugger22 Jun 7, 2017

Member

I see you've written this in your blog post

Hopefully, at some point in the future matrices will have this attribute by default.

What's the issue in adding it to Matrix now?

Member

debugger22 commented Jun 7, 2017

I see you've written this in your blog post

Hopefully, at some point in the future matrices will have this attribute by default.

What's the issue in adding it to Matrix now?

@jksuom

This comment has been minimized.

Show comment
Hide comment
@jksuom

jksuom Jun 7, 2017

Member

Matrix modules are currently being restructured, so it is better not add anything there now. Since this code is still experimental, I think it would be possible to add the attribute ring locally to simplify interfacing with smith_normal_form. It need not remain like that in the future. In principle, it could be given as a keyword argument, but eventually it should be a proper attribute set when the matrix is initialized, and smith_normal_form can be written now in a form anticipating that.

Member

jksuom commented Jun 7, 2017

Matrix modules are currently being restructured, so it is better not add anything there now. Since this code is still experimental, I think it would be possible to add the attribute ring locally to simplify interfacing with smith_normal_form. It need not remain like that in the future. In principle, it could be given as a keyword argument, but eventually it should be a proper attribute set when the matrix is initialized, and smith_normal_form can be written now in a form anticipating that.

@debugger22

This comment has been minimized.

Show comment
Hide comment
@debugger22

debugger22 Jun 8, 2017

Member

@jksuom Alright, that makes sense. Thanks!

Member

debugger22 commented Jun 8, 2017

@jksuom Alright, that makes sense. Thanks!

@jksuom

This comment has been minimized.

Show comment
Hide comment
@jksuom

jksuom Jun 19, 2017

Member

I think this is technically ready for merging. Only I would prefer the name 'invariant factors' instead of 'Smith normal factors' which does not seem to be used in the literature.

Member

jksuom commented Jun 19, 2017

I think this is technically ready for merging. Only I would prefer the name 'invariant factors' instead of 'Smith normal factors' which does not seem to be used in the literature.

@jksuom jksuom merged commit b3326ff into sympy:master Jun 19, 2017

1 check passed

continuous-integration/travis-ci/pr The Travis CI build passed
Details

@valglad valglad deleted the valglad:smith_normal branch Jun 19, 2017

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment