MatrixObj: basic operations for row/column reduction

[fingolfin]

To be able to write all kinds of algorithms that work at least somewhat efficiently with proper MatrixObj implementations, we need to allow people to write code that avoids extracting and manipulating rows. That means adding a bunch of operating facilitating this. Perhaps we already have some of these, but I couldn't find them, and they are not mentioned in chapter 26 of the ref manual... The names are of course tentative, the functionality is what is important

(CC @mohamed-barakat @ThomasBreuer @user141)

• MultMatrixRow( mat, row, scalar ) does mat[row] :=mat[row] * scalar;

• MultMatrixColumn( mat, col, scalar ) same for columns

• to support non-commutative base domains, also need an operation which does mat[row] :=scalar * mat[row] and same for columns; unfortunately we can't just overload the same function, as a positive integer could denote both a scalar and a row/column index. Perhaps MultMarixRowLeft and MultMarixRowRight ? And then have MultMatrixRow be a synonym for MultMatrixRowRight (@mohamed-barakat suggestions welcome :-)

• AddMatrixRows(mat, dstrow, srcow, scalar) does mat[dstrow] := mat[dstrow] + mat[srcrow] * scalar

• same for columns; and for scalars coming from the left (I think we assume that addition is commutative in too many places for it to make sense to allow variants there, too)

• SwapMatrixRows(mat, row1, row2)

• SwapMatrixColumns(mat, col1, col2)

We might want more (e.g. one could have PermuteMatrixRows(mat, perm)), but those listed above are an absolute minimum. If you can think of others, chime up.

[mohamed-barakat]

Invertible manipulation of rows/columns corresponds to multiplication with an invertible matrix from the left/right. So I would suggest that

MultMatrixRow( mat, row, scalar ) does mat[row] := scalar * mat[row]

while

MultMatrixColumn( mat, column, scalar ) does mat[column] := mat[column] * scalar.

The same applies to AddMatrixRows and AddMatrixColumns.

This will also cover many use-cases for noncommutative rings in a natural way :)

[fingolfin]

That's an excellent point, I like the idea very much. (To support Mohamed's point, this way the operations are symmetric under transposition - a natural way to represent a column-major matrix is to store its transpose as a row-major matrix, or vice versa).

But I wonder if we still should offer operations that perform "the other" transformation; i.e. have AddMatrixRowsLeft and AddMatrixRowsRight (but perhaps with better names)

[ssiccha]

@AnnaKDS @wucas @danielrademacher Also see the issue about parametrized tests. We could try to make some for MultMatrixRow.

[ssiccha]

Also, I think we can collect the documentation for all these functions in a (new?) section basic operations for row/column reduction.

[fingolfin]

A remark regarding performance: operations like SwapMatrixRows should ideally of course also work for "old style matrices" so that one can write code support both them and MatrixObj. But then they should be fast for such matrices, otherwise people don't want to use them. But here's the thing: being operations, they will incur method dispatch overhead, which is one of the things we meant to avoid by using MatrixObj.

Here's some naive benchmarking:

gap> Swap:=function(mat,i,j) mat{[i,j]} := mat{[j,i]}; end;;
gap> DeclareOperation("SwapMatrixRows",[IsMatrixObj,IsPosInt,IsPosInt]);
gap> InstallMethod(SwapMatrixRows,[IsMatrix,IsPosInt,IsPosInt],function(mat,i,j) mat{[i,j]} := mat{[j,i]}; end);
gap> mat:=IdentityMat(20);;
gap> for i in [1..1000000] do mat{[1,2]} := mat{[2,1]}; od; time;
138
gap> for i in [1..1000000] do Swap(mat,1,2); od; time;
193
gap> for i in [1..1000000] do SwapMatrixRows(mat,1,2); od; time;
697
gap> mat:=IdentityMat(100);;
gap> for i in [1..1000000] do mat{[1,2]} := mat{[2,1]}; od; time;
146
gap> for i in [1..1000000] do Swap(mat,1,2); od; time;
188
gap> for i in [1..1000000] do SwapMatrixRows(mat,1,2); od; time;
2385

So to counteract this, we could use the same trick as for many other operations: let SwapMatrixRows be a global function which contains code to deal with "old style matrices" (as above), and which otherwise dispatches to an operation SwapMatrixRowsOp which has all the "real" MatrixObj methods (see also issue #4432 for more on Something vs. SomethingOp).

The only thing is: it's really annoying to have to do this. Which leads me back to the idea proposed in issue #955 by @ChrisJefferson ... hrm

[ssiccha]

InstallEarlyMethod sounds like a great idea!