Opened 19 months ago
Last modified 5 weeks ago
#29886 needs_work enhancement
Add generator of symmetric matrices
Reported by:  ghIvoMaffei  Owned by:  

Priority:  major  Milestone:  sage9.6 
Component:  linear algebra  Keywords:  symmetric_matrices matrix_spaces 
Cc:  mkoeppe, tscrim  Merged in:  
Authors:  Ivo Maffei  Reviewers:  Samuel Lelièvre, ... 
Report Upstream:  N/A  Work issues:  
Branch:  u/ghIvoMaffei/symmetric_matrices (Commits, GitHub, GitLab)  Commit:  08fc8931f350c764b29aa1f76ac2a504de7cac5b 
Dependencies:  Stopgaps: 
Description (last modified by )
We add a method to the MatrixSpace
class that returns a generator over
the symmetric matrices in the space. The "symmetry" is defined by the user.
Attachments (1)
Change History (38)
comment:1 followup: ↓ 5 Changed 19 months ago by
comment:2 Changed 19 months ago by
The PEP0008 Python style guide and PEP0257 Docstring conventions are recommended reading.
I should read them again too, I realise I need a refresher on some points there.
comment:3 Changed 19 months ago by
 Commit changed from a24e6a07971f44715896523d15b8c61633ef9e52 to 9b5145b551a552ece2f31c0cf6a42b7c10d08f89
comment:4 Changed 19 months ago by
 Status changed from new to needs_review
comment:5 in reply to: ↑ 1 Changed 19 months ago by
Thanks for the detailed review.
I think I fixed all formatting issues now.
I also took the liberty to change a bit the __iter__
method.
comment:6 Changed 19 months ago by
 Description modified (diff)
 Summary changed from Generator of symmetric matrices. to Add generator of symmetric matrices
Use is None
to test if something is None
:
 if g == None: + if g is None:
I would also use f=None
to mean f
is the identity 
giving the standard "symmetric" matrices:
if f is None: f = lambda x: x
I would write the input as:
 ``f``  function or ``None`` (optional; default: ``None``); if ``None``, use the identity function  ``g``  function or ``None`` (optional; default: ``None``); if ``None``, use the same function as for ``f``
Need to decide what symmetric should mean
 Among possible definitions:
 subdiagonal entries are images by
f
of the corresponding superdiagonal entries, i.e. for any
i < j
we haveA[j,i] = f(A[i,j])
 nondiagonal entries are images by
f
of their mirror entries, i.e. for any
i != j
we haveA[j,i] = f(A[i,j])
 subdiagonal entries are images by
 These are the same if
f
is an involution; if that is assumed it needs to be said  Accordingly, do all or only half the checks
I would avoid formulas in the docstring's first sentence and clarify in a second sentence:
Return a generator of all matrices in this matrix space that are "symmetric" as defined by `f` and `g`. This means subdiagonal entries are images by ``f`` of the corresponding superdiagonal entries, i.e. for any ``i < j`` we have ``A[j,i] = f(A[i,j])``, and diagonal elements are stable by ``g`` , i.e. for any ``i`` we have ``A[i,i] = g(A[i,i])``. The default value of ``None`` for ``f`` means ``f`` is the identity, and the default value of ``None`` for ``g`` means ``g`` is ``f``.
Need ::
here to start the doctest block:
Skewsymmetric matrices::
I would use v
rather than iv
for the loop variable.
Comments could be more concise:
 # If the number of entries is zero, then just  # yield the empty matrix in that case and return + # If no entries, yield empty matrix and return if number_of_entries == 0: yield self(0) return
 # When the base ring is not finite, then we should go  # through and yield the matrices by "weight", which is  # the total number of iterations that need to be done  # on the base ring to reach the matrix. + # For infinite rings, yield matrices by "weight", + # i.e. by number of iterations needed to reach them.
 # In the finite case, we do a similar thing except that  # instead of checking if the diagonal is correct after creating the vector  # we can select all possible diagonal elements a priori + For finite rings, allow all ring elements on the diagonal a priori.
comment:7 Changed 19 months ago by
 Commit changed from 9b5145b551a552ece2f31c0cf6a42b7c10d08f89 to 08fc8931f350c764b29aa1f76ac2a504de7cac5b
Branch pushed to git repo; I updated commit sha1. New commits:
08fc893  fixed more code formatting; allow f=None for nomal symmetric matrices

comment:8 Changed 19 months ago by
Thanks for the review!! I apologise for the delay. I got caught in other tickets.
comment:9 Changed 18 months ago by
reviewer name?
comment:10 Changed 18 months ago by
@slelievre (Samuel Lelièvre) did some reviewing, but never green lighted this. They probably never felt sure about the code.
comment:11 Changed 18 months ago by
 Cc mkoeppe added
 Reviewers set to Samuel Lelièvre
Matthias, does it fit into the tensors stuff you are working on?
comment:12 Changed 18 months ago by
 Reviewers changed from Samuel Lelièvre to Samuel Lelièvre, ...
comment:13 Changed 18 months ago by
Yes, I think it does. In #30229 I define free modules of tensors with prescribed index symmetries. This covers the case of symmetric and antisymmetric matrices, of course. There is a basis and so if the base ring is an enumerated set, so is the module.
comment:14 Changed 18 months ago by
I think that skewsymmetric matrices are requested often enough for them to warrant a named function (which would just forward naturally to what you altready implemented).
comment:15 Changed 18 months ago by
I agree that it would be nice to have facilities in matrix space for symmetric and skewsymmetric matrices, but I think it would be better to define the subspaces/modules of these matrices than to define enumeration methods. Enumeration comes for free from our general facilities for enumerated sets.
comment:16 Changed 18 months ago by
So you suggest having a general SymmetricMatrixSpace
class which extends MatrixSpace
or having the "standard" SymmetricMatrixSpace
and SkewSymmetricMatrixSpace
classes and leaving everything else to an enumeration method?
comment:17 Changed 18 months ago by
 Cc tscrim added
I don't think new classes are really needed. Instead, just methods that construct submodules.
Unfortunately, currently MatrixSpace
does not fully implement the category of modules with bases (see (sage.categories.modules_with_basis
) that it claims to belong to.
sage: M = MatrixSpace(QQ, 3, 3) sage: M in ModulesWithBasis(QQ) True sage: A = M([[0, 0, 1], [0, 0, 0], [0, 0, 0]]) sage: A [0 0 1] [0 0 0] [0 0 0] sage: M.submodule([A]) AttributeError: 'MatrixSpace_with_category' object has no attribute 'get_order'
I would suggest to start with fixing the implementation of this functionality of MatrixSpace
.
Then in the next step, you can define the methods that construct the specific modules that you need.
comment:18 Changed 18 months ago by
I am inclined to wave this through, as it's needed for Ivo's GSoC project, and the time is tight.
comment:19 Changed 18 months ago by
I don't object to adding these specific enumeration methods,
but I think adding this interface with the generality of this f
and g
thing is not a good idea, because a future implementation along the lines that we just discussed will not be able to support general f and g.
comment:20 Changed 18 months ago by
See #30276 for the type of symmetry that we hope to support.
comment:21 Changed 18 months ago by
 Status changed from needs_review to needs_work
I usually thing of the sub/superdiagonal matrices as just being the entries immediately off the diagonal, not in all of the lower/upper triangular part.
A lot of what this method does is really confusing to me. The algorithm should be explained in the docstring, in particular the order which it is iterating through everything.
I also somewhat agree with Matthias here. I would first expect this to return a subspace, which doesn't work for generic f
and g
. Plus the word "symmetric" seems misleading to me for generic f
and g
. Perhaps as a middle ground, we call this symmetric_matrices_iterator()
.
Comments from a quick look at the code:
I don't see why you need to do v2.clone()
since you never modify it. You simply get a sublist that to reassign to v
.
self(0)
> self.zero()
.
This is very strange: for _ in range(nrows):
. I don't understand this code block.
comment:22 followups: ↓ 24 ↓ 32 Changed 18 months ago by
I suggest to implement here only the minimum needed for the GSoC tickets to be done (on a different branch, to preserve this, I suppose).
comment:23 Changed 18 months ago by
I spent some time looking at the submodule
method to see if there was some easy fix. I think I got 90% of the way.
There were some issues in the method sage.categories.finite_dimensional_modules_with_basis.echelon_form
.

src/sage/categories/finite_dimensional_modules_with_basis.py
diff git a/src/sage/categories/finite_dimensional_modules_with_basis.py b/src/sage/categories/finite_dimensional_modules_with_basis.py
a b class FiniteDimensionalModulesWithBasis(CategoryWithAxiom_over_base_ring): 337 337 sage: C.echelon_form([x[0]  x[1], 2*x[1]  2*x[2], x[0]  x[2]]) 338 338 [x[0]  x[2], x[1]  x[2]] 339 339 """ 340 340 if order is not None: 341 order = self._compute_support_order(order) 341 order = self._compute_support_order(elements, order) 342 343 # order may be an ordering of the support of elements 344 # hence not an ordering of the whole basis as needed below 345 # so we extend order to an ordering of the whole basis 346 orderList = list(order) 347 for k in self.basis().keys(): 348 if k not in orderList: 349 orderList.append(k) 350 351 full_order = tuple(orderList) 352 353 342 354 from sage.matrix.constructor import matrix 343 mat = matrix(self.base_ring(), [g._vector_(order= order) for g in elements])355 mat = matrix(self.base_ring(), [g._vector_(order=full_order) for g in elements]) 344 356 # Echelonizing a matrix over a field returned the rref 345 357 if row_reduced and self.base_ring() not in Fields(): 346 358 try:
After applying that patch there is an issue with mutable matrices not being hashable. The issue seems to lie in sage.misc.cachefunc
and I don't have enough knowledge of python or the inner working of Sage to proceed.
Honestly, I had a huge amount of issues with matrices/vectors not being hashable, so I wonder why aren't they immutable by default.
comment:24 in reply to: ↑ 22 Changed 18 months ago by
comment:25 Changed 18 months ago by
#30334 is not urgent.
Given that submodule
does not work, my suggestion for the present ticket would be the following. First implement a method named perhaps MatrixSpace.symmetric_submodule_basis(phase=1)
. As in #30276, for phase=1
, return a basis of the submodule of the symmetric matrices, and for phase=1
, skewsymmetric matrices. Error for all other values of phase
.
comment:26 Changed 18 months ago by
I don't fully understand #30334. What I'm using this method for, is to generate skewsymmetric matrices with diagonal entries 0 (the standard definition A^T = A
would allow diagonal entries 1 in GF(2)) and Hermitian matrices.
Can these be achieved for some phase
?
comment:27 followup: ↓ 28 Changed 18 months ago by
That's important information  from looking at the code it was not clear to me what generality you need.
For Hermitian matrices  over what field?
comment:28 in reply to: ↑ 27 Changed 18 months ago by
Replying to mkoeppe:
For Hermitian matrices  over what field?
Finite fields GF(r^2)
. Practically only GF(4)
and GF(9)
.
comment:29 Changed 18 months ago by
Unfortunately the Hermitian matrices over a field k do not form a subspace (over k)...
comment:30 Changed 18 months ago by
They form a subgroup and they seem quite common... If you think there is no place for a generator of some sort here, then I'll just have the hermitian matrices created where I use them, so that nothing is exposed to the user.
comment:31 Changed 18 months ago by
well, Hermitian matrices form a subspace over the fixed (by the automorphism) subfield. Probably our matrices aren't flexible enough for this, though.
comment:32 in reply to: ↑ 22 Changed 17 months ago by
Replying to dimpase:
I suggest to implement here only the minimum needed for the GSoC tickets to be done (on a different branch, to preserve this, I suppose).
To move forward, I'd suggest to move this methods to the code where they are needed. It can be be refactored later when we figure out the right framework.
comment:33 Changed 15 months ago by
 Milestone changed from sage9.2 to sage9.3
comment:34 Changed 11 months ago by
 Milestone changed from sage9.3 to sage9.4
Setting new milestone based on a cursory review of ticket status, priority, and last modification date.
comment:35 Changed 7 months ago by
comment:36 Changed 6 months ago by
 Milestone changed from sage9.4 to sage9.5
Setting a new milestone for this ticket based on a cursory review.
comment:37 Changed 5 weeks ago by
 Milestone changed from sage9.5 to sage9.6
Remove space before closing parenthesis:
Typo: doesn't contains > doesn't contain.
Use a space after each comma in code, and around
+
signs.Use double backticks, e.g.
``f``
or``ValueError``
for code formatting in docstrings.Typo "funcition" > function.
Use code formatting for
None
:if it is None
>if it is ``None``
Typo:
MetricSpace
>MatrixSpace
.Hint: after making changes, rebuild Sage and test with
before pushing. That would catch things such as
MetricSpace
.Replace this doctest whose output is hard to parse:
by this one whose output takes up less space and is easier to parse:
Likewise, replace this:
with this more concise:
Remove spaces around
M
indef make_symmetric( M ):
.Simplify "can't have symmetric matrices if they are not square" to "symmetric matrices must be square" or "only square matrices can be symmetric".
Use a space after
#
for comments, for instance:For inline comments, double space before
#
and single space after. For instance:No need to comment
#make M symmetric
beforemake_symmetric(M)
, function name says it all.Rather than importing
and using
sage.combinat.integer_vector.IntegerVectors
, instead importand then use
IntegerVectors
.