Add generator of symmetric matrices

[gh-Ivo-Maffei]

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.

[slelievre]

Remove space before closing parenthesis:

- def symmetric_matrices(self, f, g=None ):
+ def symmetric_matrices(self, f, g=None):

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``

- skew-symmetric matrices:
+ Skew-symmetric matrices:

Typo: MetricSpace -> MatrixSpace.

Hint: after making changes, rebuild Sage and test with

sage -t src/sage/matrix/matrix_space.py

before pushing. That would catch things such as MetricSpace.

Replace this doctest whose output is hard to parse:

sage: gen = MS.symmetric_matrices(f)
sage: for M in gen:
....:     print(M)

by this one whose output takes up less space and is easier to parse:

sage: list(MS.symmetric_matrices(f))

Likewise, replace this:

sage: gen = MS.symmetric_matrices(f)
sage: i = 0
sage: for M in gen:
....:     i+= 1
....:     print(M)
....:     if i == 3: break

with this more concise:

sage: gen = MS.symmetric_matrices(f)
sage: [next(gen) for _ in range(3)]

Remove spaces around M in def 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:

- #When the base ring is not finite,
+ # When the base ring is not finite,

For inline comments, double space before # and single space after. For instance:

matrix_entries = []  # entries of matrix with lower half 0

No need to comment #make M symmetric before make_symmetric(M), function name says it all.

Rather than importing

import sage.combinat.integer_vector

and using sage.combinat.integer_vector.IntegerVectors, instead import

from sage.combinat.integer_vector import IntegerVectors

and then use IntegerVectors.

[slelievre]

The ​PEP-0008 Python style guide and ​PEP-0257 Docstring conventions are recommended reading.

I should read them again too, I realise I need a refresher on some points there.

[git]

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

​1b38de3	fixed formatting
​9b5145b	removed trailing whitespaces

[gh-Ivo-Maffei]

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.

[slelievre]

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 have A[j,i] = f(A[i,j])
  • nondiagonal entries are images by f of their mirror entries, i.e. for any i != j we have A[j,i] = f(A[i,j])
• 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:

Skew-symmetric 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.

[git]

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

​08fc893

fixed more code formatting; allow f=None for nomal symmetric matrices

[gh-Ivo-Maffei]

Thanks for the review!! I apologise for the delay. I got caught in other tickets.

[dimpase]

reviewer name?

[gh-Ivo-Maffei]

@slelievre (Samuel Lelièvre) did some reviewing, but never green lighted this. They probably never felt sure about the code.

[dimpase]

Matthias, does it fit into the tensors stuff you are working on?

[mkoeppe]

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.

[dimpase]

I think that skew-symmetric matrices are requested often enough for them to warrant a named function (which would just forward naturally to what you altready implemented).

[mkoeppe]

I agree that it would be nice to have facilities in matrix space for symmetric and skew-symmetric 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.

[gh-Ivo-Maffei]

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?

[mkoeppe]

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.

[dimpase]

I am inclined to wave this through, as it's needed for Ivo's GSoC project, and the time is tight.

[mkoeppe]

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.

[mkoeppe]

See #30276 for the type of symmetry that we hope to support.

[tscrim]

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.

[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).

[gh-Ivo-Maffei]

mutable matrices issue full traceback

[gh-Ivo-Maffei]

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

@@ -337 +337 @@
                 sage: C.echelon_form([x[0] - x[1], 2*x[1] - 2*x[2], x[0] - x[2]])
                 [x[0] - x[2], x[1] - x[2]]
             """
             if order is not None:
-                order = self._compute_support_order(order)
+                order = self._compute_support_order(elements, order)
+
+            # order may be an ordering of the support of elements
+            # hence not an ordering of the whole basis as needed below
+            # so we extend order to an ordering of the whole basis
+            orderList = list(order)
+            for k in self.basis().keys():
+                if k not in orderList:
+                    orderList.append(k)
+
+            full_order = tuple(orderList)
+
+            
             from sage.matrix.constructor import matrix
-            mat = matrix(self.base_ring(), [g._vector_(order=order) for g in elements])
+            mat = matrix(self.base_ring(), [g._vector_(order=full_order) for g in elements])
             # Echelonizing a matrix over a field returned the rref
             if row_reduced and self.base_ring() not in Fields():
                 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.

[mkoeppe]

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).

I have opened #30334 for the MatrixSpace.submodule issue.

[mkoeppe]

#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, skew-symmetric matrices. Error for all other values of phase.

[gh-Ivo-Maffei]

I don't fully understand #30334. What I'm using this method for, is to generate skew-symmetric 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?

[mkoeppe]

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?

[gh-Ivo-Maffei]

Replying to mkoeppe:

For Hermitian matrices - over what field?

Finite fields GF(r^2). Practically only GF(4) and GF(9).

[mkoeppe]

Unfortunately the Hermitian matrices over a field k do not form a subspace (over k)...

[gh-Ivo-Maffei]

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.

[dimpase]

well, Hermitian matrices form a subspace over the fixed (by the automorphism) subfield. Probably our matrices aren't flexible enough for this, though.

[mkoeppe]

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.

[mkoeppe]

Setting new milestone based on a cursory review of ticket status, priority, and last modification date.

[slelievre]

Submodules of matrix spaces are taken care of at #31995 (superseding #30334).

[mkoeppe]

Setting a new milestone for this ticket based on a cursory review.