Inverse of a non-invertible matrix over the integers is a (non-integer) integer matrix

[james-d-mitchell]

In the current GAP master branch I get the following, inverting a non-invertible integer MatrixObj results in a matrix that believes it is an integer matrix but whose entries are not integers.

Observed behaviour

gap> mat := Matrix(Integers, [[0, -3, 0, -2], [-1, 1, -1, 0],
> [0, 1, 0, 1], [0, 0, 2, 0]]);;
gap> InverseOp(mat);
<4x4-matrix over Integers>
gap> Display(last);
<4x4-matrix over Integers:
[[ -1, -1, -2, -1/2 ]
 [ -1, 0, -2, 0 ]
 [ 0, 0, 0, 1/2 ]
 [ 1, 0, 3, 0 ]
]>

Expected behaviour

gap> mat := Matrix(Integers, [[0, -3, 0, -2], [-1, 1, -1, 0],
> [0, 1, 0, 1], [0, 0, 2, 0]]);;
gap> InverseOp(mat);
fail

or perhaps at a stretch:

gap> mat := Matrix(Integers, [[0, -3, 0, -2], [-1, 1, -1, 0],
> [0, 1, 0, 1], [0, 0, 2, 0]]);;
gap> InverseOp(mat);
<4x4-matrix over Rationals>
gap> Display(last);
<4x4-matrix over Rationals:
[[ -1, -1, -2, -1/2 ]
 [ -1, 0, -2, 0 ]
 [ 0, 0, 0, 1/2 ]
 [ 1, 0, 3, 0 ]
]>

[james-d-mitchell]

Here's another example of an integer matrix with non-integer entries:

gap> Matrix(Integers, [[1, 1], [2, E(8)]]);
<2x2-matrix over Integers>
gap> IsInt(E(8));
false

Again maybe this is the intended behaviour, but this is unexpected, I'd have expected this to throw an error.

[hulpke]

This is an interesting question: What should Inverse do for a matrix declared over a ring, which is invertible over the fraction field? And should Matrix check whether the given entries actually lie over the ring? @ThomasBreuer might know.

[ThomasBreuer]

Concerning the example Matrix(Integers, [[1, 1], [2, E(8)]]), this input must yield either an error or fail, and the documentation of Matrix must state what happens in such cases.
Currently the documentation of Matrix just states that the returned matrix is defined over the given base domain, and the documentation of BaseDomain refers to an introductory section that states that the matrix entries must either lie in the base domain or naturally embed into the base domain.
Thus Matrix must check this property, and perhaps something like MatrixNC should be provided, which omits these checks.

The first example is more subtle. Again, the documentation should be made precise about what shall happen. The role of the BaseDomain value of matrix objects seems to be not defined very restrictively. There is the following statement about the idea behind vector and matrix objects.

We want to be able to write (efficient) code that is independent of the actual representation (in the sense of GAP's representation filters, see Section 13.4) and preserves it.

If both Integers and Rationals are supported for a given representation filter of matrix objects then one could argue that it would be allowed to create inverses over different base domains.
On the other hand, this would be against the idea of getting efficient code. And there is already code (for example a KroneckerProduct method) that accepts only arguments over the same base domain.
The documentation must state what is allowed in general, for example what the sum and the product of two matrix objects over different base domains are.
My understanding is that operations such as Inverse should not try to create matrix objects with new (guessed) base domains but should signal an error or return fail; for the sake of consistency (since we can think of Inverse calling Matrix for creating its result), I think that the behaviour should be the same as for Matrix(Integers, [[1, 1], [2, E(8)]]).

[james-d-mitchell]

I agree with @ThomasBreuer, I think the "correct" behaviour would be for Inverse and Matrix(Integers, [[1, 1], [2, E(8)]]) to return fail or give an error, and that there should be a "no check" version called MatrixNC.

[hulpke]

That seems to indicate that we should introduce MatrixNC, rename the existing method to MatrixNC, and add a wrapper function Matrix that tests for membership in the domain. I think as for using it the extra cost is relevant, as long as matrix operations call MatrixNC, the overhead only occurs onve when creating the initial objects. Does this seem reasonable?

[james-d-mitchell]

Seems completely reasonable to me @hulpke, it's what we did in Semigroups and it's the test failures there that indicated this issue to me when converting from our own matrices to MatrixObjs.

[fingolfin]

I agree with @ThomasBreuer and @james-d-mitchell regarding what should happen. I like the idea of MatrixNC. We may also need NewMatrixNC (note that NewMatrix is intended to be the low-level way to construct a matrix, while Matrix methods are intended to be implemented in terms of NewMatrix -- however, I don't think we make that clear in the docs right now (because it wasn't always like this).

[fingolfin]

Likewise IdentityMatrix is implemented in terms of NewIdentityMatrix etc. -- the idea here being that not every MatrixObj implementation should have to implement all variants of IdentityMatrix etc. , they only need one; and then we also offer defaults for e.g. NewIdentityMatrix which delegate to NewMatrix. So that a basic MatrixObj implementation should be able to get away with only implementing a single NewMatrix constructor (modulo bugs and omissions in our code for ensuring that, of course). All of this ought to be documented, though sigh

[fingolfin]

I submitted issue #4891 to remind us of this