Add IsGenericMatrixRep matrix object type

[fingolfin]

Store dense matrices as plain row lists for faster entry access.

Partially addresses issue #2148 (the unresolved part is that of a "vector object that is not in IsList)
Closes #2973 (an older attempt of something similar)

This PR was created with help of Codex by OpenAI: I planned the work with it, then it did a full first implementation, which was improved by me at first via some rounds of prompts requesting tweaks, then a bunch of manual work to implement details and speedups the AI was not aware.

Since this is fresh, I am sure there are some bugs, and definitely opportunities for things to be improved. And of course plenty of places in the code where using these MatrixObj will blow up because they are not row lists... but that's for future work.

---

Some microbenchmarks show the benefit, namely that performance for basic arithmetic in many cases matches that of ist-of-list matrices, while IsPlistMatrixRep is basically always slower, sometimes quite substantially so.

I am deliberately doing this over a finite field that is not covered by our compressed matrices, as that is IMHO going to be one of the central use cases in GAP; though for integer/rational matrices, I don't think things will be substantially difference, at least relative to each other. There is no point comparing e.g. over GF(2) to our compressed GF(2) matrices, which of course will be much faster than any of the three implementations being tested here.

Setup:

R:=GF(257);;n:=20;;
M1:=IdentityMat(n,R);;
M2:=IdentityMatrix(IsPlistMatrixRep,R,n);;
M3:=IdentityMatrix(IsFlatPlistMatrixRep,R,n);;

Addition:

gap> for i in [1..1000] do x:=M1+M1; od; time;
4
gap> for i in [1..1000] do x:=M2+M2; od; time;
12
gap> for i in [1..1000] do x:=M3+M3; od; time;
4

Multiplication:

gap> for i in [1..1000] do x:=M1*M1; od; time;
66
gap> for i in [1..1000] do x:=M2*M2; od; time;
1448
gap> for i in [1..1000] do x:=M3*M3; od; time;
60

Inversion:

gap> for i in [1..1000] do x:=Inverse(M1); od; time;
7
gap> for i in [1..1000] do x:=Inverse(M2); od; time;
83
gap> for i in [1..1000] do x:=Inverse(M3); od; time;
4

Powers:

gap>
gap> for i in [1..1000] do x:=M1^5; od; time;
187
gap> for i in [1..1000] do x:=M2^5; od; time;
4350
gap> for i in [1..1000] do x:=M3^5; od; time;
182

Element access (plist-of-plist benefits from a direct kernel implementation; we could add that for both IsPlistMatrixRep and IsFlatPlistMatrixRep, if so desired)

gap> for i in [1..1000000] do x:=M1[1,1]; od; time;
19
gap> for i in [1..1000000] do x:=M2[1,1]; od; time;
78
gap> for i in [1..1000000] do x:=M3[1,1]; od; time;
66

Example of a typical matrix property: IsDiagonalMat

gap> for i in [1..10000] do x:=IsDiagonalMat(M1); od; time;
132
gap> for i in [1..10000] do x:=IsDiagonalMat(M2); od; time;
191
gap> for i in [1..10000] do x:=IsDiagonalMat(M3); od; time;
146

IsOne (surprisingly is slower for IsFlatPlistMatrixRep, I do not yet know why)

gap> for i in [1..10000] do x:=IsOne(M1); od; time;
144
gap> for i in [1..10000] do x:=IsOne(M2); od; time;
145
gap> for i in [1..10000] do x:=IsOne(M3); od; time;
172

PositionNonZeroInRow is among the core functionality for row reduction:

gap> for i in [1..1000000] do x:=PositionNonZeroInRow(M1,1); od; time;
334
gap> for i in [1..1000000] do x:=PositionNonZeroInRow(M2,1); od; time;
486
gap> for i in [1..1000000] do x:=PositionNonZeroInRow(M3,1); od; time;
376

[ThomasBreuer]

just a few minor remarks

if we want to describe this representation (in comparison with others) then my impression is:

• a matrix type that can be multiplied with any type of vector object, but formally defines IsPlistVectorRep as its CompatibleVectorFilter,
• a matrix type that deliberately does not support row access,
• matrix type that supports matrices with zero rows or columns (something which should hold for all matrix types except IsPlistRep.

for curiosity:
Where does the name IsFlatPlistMatrixRep come from?
GAP has a function Flat, thus one could get the idea that an m times n "flat" matrix is internally represented by a "flat" list of length mn.

[fingolfin]

The name IsFlatPlistMatrixRep comes from some recent discussion thread involving you and James, someone brought it up. But yeah, I do see the confusion with Flat.

In PR #6227 I suggested that we could rename IsPlistMatrixRep to IsRowPlistMatrixRep; and then this new type could become the new IsPlistMatrixRep. That's at least in principle viable since no GAP package uses IsPlistMatrixRep for their implementation; in fact only forms uses it in two .tst files, no other code has active code involving it.

However, it might break private user code, so I am not keen on that.

I am open to alternatives. Some braingstorming with AI yields these, but perhaps someone has a better idea:

1. IsDirectPlistMatrixRep
2. IsNestedPlistMatrixRep
3. IsPlainPlistMatrixRep
4. IsGenericMatrixRep

[ThomasBreuer]

I like the idea IsGenericMatrixRep, from the viewpoint that there are no conditions of the base domain and thus this representation can be used for all computation in the context of "matrices over some base domain".

When one wants to "upgrade" code that was written for list-of-lists matrices to matrix objects, switching to IsGenericMatrixRep could be suggested as a first step: Replace the input matrices by matrices in IsGenericMatrixRep, run the code, and fix the errors that happen due to the (mis)use of features of list-of-lists matrices.

[james-d-mitchell]

for curiosity: Where does the name IsFlatPlistMatrixRep come from? GAP has a function Flat, thus one could get the idea that an m times n "flat" matrix is internally represented by a "flat" list of length mn.

I have the same question as @ThomasBreuer, I also expected this. Could you possibly indicate what the key difference between matrices in IsFlatPlistMatrixRep and those in IsPlistMatrixRep is @fingolfin? Is it that the rows stored in matrix aren't matrices of the same type as the matrix itself but rather "just" plain lists? Is this also the reason for not providing access to the rows?

[ThomasBreuer]

Concerning the difference between IsFlatPlistMatrixRep and IsPlistMatrixRep:
The idea of IsPlistMatrixRep was that the internal data are a plain list of vector objects in IsPlistVectorRep. In this situation, it is natural to define that these vector objects are the rows.
For IsFlatPlistMatrixRep, the idea is not to support access to rows; there are no natural vector objects for that purpose, and the advantage is that the matrix can handle the multiplication with any kind of vector object, by creating an object of the same representation as the result.

[fingolfin]

Everything @ThomasBreuer wrote is true, but the main advantage of this approach is speed: by not having every row be a "rich" and thus expensive object, we avoid a lot overhead in both memory usage and speed. Also, we can directly call into GAP kernel routines for e.g. multiplying two "list of list matrices". As a result, the new type can be almost as fast as the classic "list of list matrices", and in some cases even faster (as it avoids expensive method dispatch overhead)

[fingolfin]

I've now renamed this to IsGenericMatrixRep, as discussed.

[fingolfin]

@ThomasBreuer @james-d-mitchell if there are no further objections (also to the new name), I'll merge this later today

[james-d-mitchell]

@ThomasBreuer @james-d-mitchell if there are no further objections (also to the new name), I'll merge this later today

I have no objections but also haven't looked in detail at this. I can do that if desired but not until later today.

[ThomasBreuer]

As far as I see, the changes refer only to the renaming. Since I had voted for the new name, I do not object.

[fingolfin]

@james-d-mitchell here is my view: if you would like to have a look at this, I'd be happy to wait for you, also beyond today! But don't feel compelled to look at it -- in the end @ThomasBreuer had a look, and even if I messed up a lot, we can still revert or fix things :-)

[ThomasBreuer]

One more comment:

Since I do now understand the implementation of IsZmodnZMatrixRep better, I think it would be appropriate to

• change IsZmodnZMatrixRep to store a list of lists of reduced integers, not a list of IsZmodnZVectorRep objects,
• support the multiplication of such a IsZmodnZMatrixRep object with a vector (from both sides) only for vectors in IsZmodnZVectorRep, and
• do not allow the access to rows of such a IsZmodnZMatrixRep object until we see where this turns out to be useful (at the cost of creating a new vector object).

[james-d-mitchell]

It'd be good to clarify what the difference is between IsPlistMatrixRep and IsGenericMatrixRep is here, as described in one of the comments elsewhere.

[james-d-mitchell]

For example, although I know what is intended "a dense plain list of dense plain row lists" and "is not a row list matrix" seem contradictory. More words might indicate more clearly what is meant here.

[james-d-mitchell]

Does the leading "F" stand for "Flat"? If so it should probably be "G" for generic no?

[james-d-mitchell]

Is there a reason that this exists? It seems to only be used in NumberColumns couldn't is just be Length(m![FROWSPOS][1]) instead? Or I suppose

if Length(m![FROWSPOS]) = 0 then
return 0;
else
return Length(m![FROWSPOS][1]);
fi;

If this is for performance reasons, then maybe a comment indicating that?

[fingolfin]

Then we'd loose support for 0 x m matrices, which is really annoying in lots of situations in my practical experience.

[james-d-mitchell]

This isn't uniformly formatted see line 76 for example

[james-d-mitchell]

Suggested change

	Error( "<row> is out of bounds" );
	Error( "<row> ", row, " is out of bounds, expected value in [", 1, ", ", NrRows(M),"]");

[james-d-mitchell]

Suggested change

	Error( "<M> and <N> are not compatible" );
	Error( "the 1st and 2nd arguments <M> and <N> are not compatible, their BaseDomains must be identical objects" );

[fingolfin]

The user will be on a break loop, where M and N can be accessed, and if they so care, they can also see that they are argument 1 and 2 of the function in which the error is raised; so I really don't like adding these position information to error message.

[james-d-mitchell]

This could be more nuanced too.

[fingolfin]

Yes it could be. But this is a first iteration. We might decided to discard it again in a couple months. I am not sure it is worth it polishing it to this degree before we even have a single place (other than the tests) using it?

In fact I am tempted to also remove it from the documentation, I think it is a mistake to put it there right now, until we have actually used it a bit and decide to keep it.

[james-d-mitchell]

Suggested change

	ErrorNoReturn( "\\*: Matrices do not fit together" );
	ErrorNoReturn( "\\*: the factors (matrices) do not fit together, the 1st factor has ", colsA, "columns, and the 2nd factor has ", rowsB, " rows, expected these to be equal");

or something similar.

[james-d-mitchell]

These error messages could be more descriptive or the if-statements could be combined into one.

[james-d-mitchell]

All of these 1 line functions could be lambdas to save "vertical screen real estate"

[fingolfin]

Some but not all (just those with a return)

[james-d-mitchell]

Looks slightly weird to me that the opening bracket is in the same line as the first row, but the closing bracket is essentially in a line by itself.

[james-d-mitchell]

There are some remnants of the original name "Flat" that should be changed to "Generic" to avoid confusion in the future. The documentation should also be updated to clearly explain what the difference between IsGenericMatrixRep and IsPlistMatrixRep. Otherwise, I've made some suggestions that I think might improve the error messages/usability of the new code, which are less important.