homalg-project / homalg_project

Deposited packages of the homalg project
11 stars 18 forks source link

[Gauss] incorrect internal representation of sparse matrices over an explicitly specified ring #606

Open jdobben opened 4 months ago

jdobben commented 4 months ago

The package gauss doesn't seem to always properly store matrices that are constructed over an explicitly specified ring. For example, I thought that the following would be a valid way of constructing a matrix over GF(3):

LoadPackage("gauss");
SN := SparseMatrix( 2, 3, [ [ 2 ], [ 1, 3 ] ], [ [ 1 ], [ 3, 2 ] ], GF(3) );

The documentation doesn't specify that this is incorrect, so I assumed that this call would automatically convert the integer entries to elements of GF(3). However, when recovering the matrix from gauss, it becomes clear that no conversion has taken place:

gap> N := ConvertSparseMatrixToMatrix(SN);
[ [ 0*Z(3), 1, 0*Z(3) ], [ 3, 0*Z(3), 2 ] ]

On rare occasions, this could lead to problems when doing further computations with these sparse matrices. For instance, I have had a case where computing the rank (using the Rank function from gauss) failed because the algorithm ended up trying to compute -1/3 mod 3, presumably because earlier steps in the Gaussian elimination algorithm had set that pivot to some rational number (like 3 or -1/3) instead of an element of GF(3). However, I have not been able to produce a small example that exhibits this behaviour, and I am currently not at liberty to share the code that led to a failed rank computation.

mohamed-barakat commented 4 months ago

According to the specification, the ring elements of the argument entries must lie in ring. This is why your input is invalid. A valid input would be:

gap> LoadPackage("gauss", false);
gap> SN := SparseMatrix( 2, 3, [ [ 2 ], [ 3 ] ], [ [ 1 ], [ 2 ] ] * One(GF(3)), GF(3) );
gap> Rank(SN);
jdobben commented 4 months ago

By "specification", do you mean the package documentation? Because I don't see it there, and I assumed from the documentation that there would be automatic type casting when specifying the ring (why else would that option be there?).

To be sure, I am aware that type casting is certainly not in the documentation; it's just what I expected to happen. Maybe my problem is that I'm carrying over expectations from other programming languages? 😅 Like, I would expect any constructor to automatically recognize and convert the input to whatever internal representation it uses, or raise an error if it doesn't understand the input data. I do not expect it to quietly go along and cause problems later. 🤔 If performance is an issue, I feel that it would be more appropriate to have a separate SparseMatrixNC function which skips the check (similar to the NC versions of many other GAP functions).

I guess my three recent issues (#604, #605, #606) can be marked as duplicates, as they all boil down to a feature request to:

mohamed-barakat commented 4 months ago

This is in fact what HomalgMatrix does, using SparseMatrix as one of many low-level matrix constructors.

I am happy about PR improving the documentation of Gauss, but am reluctant to convenience in SparseMatrix.