Faster implementation of is_invertible() by checking full rank

[HenryWu2101]

#31274
As mentioned in the issue, a faster invertibility check for matrix in Matrix_gf2e_dense by checking full rank.

📝 Checklist

• The title is concise and informative.
• The description explains in detail what this PR is about.
• I have linked a relevant issue or discussion.
• I have created tests covering the changes.
• I have updated the documentation and checked the documentation preview.

⌛ Dependencies

[github-actions]

Documentation preview for this PR (built with commit cce4350; changes) is ready! 🎉
This preview will update shortly after each push to this PR.

[user202729]

Could probably be implemented for all matrix over exact fields (such as QQ)

[user202729]

On another note, the check was originally added in #30161, since the bug has been fixed in m4rie the check in __invert__ should be redundant anyway.

malb/m4rie@ed0f25e

(actually it suffices to check if element at [n-1][n-1] is 1, but that's out of scope here.)

[HenryWu2101]

On another note, the check was originally added in #30161, since the bug has been fixed in m4rie the check in __invert__ should be redundant anyway.

malb/m4rie@ed0f25e

Hi so I've checked this note, totally agree on the checks part and adjusted accordingly. This should be reflected in the latest commit and should be ok if all checks pan out.

Could probably be implemented for all matrix over exact fields (such as QQ)

I've checked matrix_rational_dense and matrix_integer_dense, both of which integrated the "check full-rank" and "check square" in their _invert_ meth::, so I'm not sure if separately creating a new invertible check and change the original invert function is a good idea. Could still be done though.

[vneiger]

This is changing behaviour: until now, for this type of matrices, one could call inverse on a non-square matrix without issue:

sage: mat = matrix.random(GF(4), 2, 4)
sage: mat.inverse()
[ 0 z2  0  0]
[ 1  0  0  0]
sage: mat
[     0      1 z2 + 1      0]
[z2 + 1      0 z2 + 1 z2 + 1]

(I would say that this behaviour change is actually a good thing, but this means at least checking that this is not breaking anything elsewhere in Sage.)

[vneiger]

Suggestion: use the same description as the general function.

Suggested change

	Return if invertible by checking full rank.
	Return "True" if this matrix is invertible.

[user202729]

use ​`​`True`​`​ instead.

[vneiger]

Could probably be implemented for all matrix over exact fields (such as QQ)

For some exact fields it might help, but why would it help for all of them? Here it seems to help mostly because there is no determinant method for this type of matrices, so sage is calling a slow determinant to determine if the matrix is invertible (this is my understanding after a quick review, please correct me if you have better insight). For example over prime fields or the rationals, why would this approach (computing the rank) be better than computing the determinant?

[vneiger]

Hi so I've checked this note, totally agree on the checks part and adjusted accordingly. This should be reflected in the latest commit and should be ok if all checks pan out.

If I follow the comment of @user202729 correctly, I think that the check that can be removed is the is_invertible part. From what I see in the mentioned commit I am wondering whether you would rather have removed the other check (the one about the number of rows being not zero).

[HenryWu2101]

Hi so I've checked this note, totally agree on the checks part and adjusted accordingly. This should be reflected in the latest commit and should be ok if all checks pan out.

If I follow the comment of @user202729 correctly, I think that the check that can be removed is the is_invertible part. From what I see in the mentioned commit I am wondering whether you would rather have removed the other check (the one about the number of rows being not zero).

Yes I'm afraid that's right. At first I thought the check in newton_john is checking for 0 row, but on 2nd look it's checking for full rank. Will revert that change.

Could probably be implemented for all matrix over exact fields (such as QQ)

For some exact fields it might help, but why would it help for all of them? Here it seems to help mostly because there is no determinant method for this type of matrices, so sage is calling a slow determinant to determine if the matrix is invertible (this is my understanding after a quick review, please correct me if you have better insight). For example over prime fields or the rationals, why would this approach (computing the rank) be better than computing the determinant?

From what I see in either rational_dense or integer_dense, it checks for invertible still via checking square + checking full rank (though integrated somehow in different functions), even both have the determinant method. I guess that was initially designed so; no class that I know of checks for det not 0.

[vneiger]

From what I see in either rational_dense or integer_dense, it checks for invertible still via checking square + checking full rank (though integrated somehow in different functions), even both have the determinant method. I guess that was initially designed so; no class that I know of checks for det not 0.

The general method in matrix0.pyx defines is_invertible as the result of the test self.is_square() and self.determinant().is_unit() . So any subclass that does not override this method is going to use the check for determinant not zero (or more precisely, determinant invertible, in case the matrix is not over a field). That is the case for matrices modn_dense, i.e. over Z/nZ for an integer n (including prime fields when n is prime).

In fact, for matrices in rational_dense, I'm not sure to follow. I don't see an implementation of is_invertible there. There are some checks at the beginning of inversion methods, but I don't see a check that computes the rank. Could you please clarify what part of the code makes you write that "it checks for invertible still via checking square + checking full rank"? I'm asking because if some inversion function checks the rank first, this would hint at potential suboptimal code that would be (hopefully) easily fixed.

[user202729]

For some exact fields it might help, but why would it help for all of them?

The naive algorithm for determinant over (exact) field is to compute the row echelon form then multiply entries on the diagonal, and for rank is to compute the row echelon form then count number of nonzero entries on the diagonal. Both are O(n^3) but you eliminate some multiplication?

Or is it because there happen to have more efficient rank method for this particular class of matrix? Or some other reason?

[vneiger]

For some exact fields it might help, but why would it help for all of them?

The naive algorithm for determinant over (exact) field is to compute the row echelon form then multiply entries on the diagonal, and for rank is to compute the row echelon form then count number of nonzero entries on the diagonal. Both are O(n^3) but you eliminate some multiplication?

Ok, I see the reasoning, thanks for the explanation. Yes we can expect several types of matrices over exact fields to use this strategy: the bulk of the work is done by the computation of some echelon form or triangularization (including "non-naive" fast PLE / PLUQ decompositions), and then the rank or determinant is deduced rather directly. For the determinant this multiplication step is only O(n) multiplications, so it should not influence performance much, except possibly in rare corner cases.

For matrices over finite fields the different types in Sage (based on FFLAS-FFPACK or on M4RI or on M4RIE from what I see) all use some fast triangularization, so indeed it should be a bit cleaner to rely on the rank. On the other hand, for other types of matrices especially when coefficient growth is possible, things are less clear to me. And for matrices not over a field, one should be extra careful (e.g. matrix(Integers(6), [[2]]).rank() raises an error, but there's no issue with determinants).

Or is it because there happen to have more efficient rank method for this particular class of matrix? Or some other reason?

In the present case of this issue (matrix over GF(2**e)), the rank method is significantly more efficient than the determinant one, but not for a good (algorithmic) reason: this only seems to be due to the absence of a dedicated determinant method. Yesterday I checked the M4RIE code and could not find a determinant procedure there (but several fast triangularization procedures, so one can expect that adding the determinant will not be difficult).

[user202729]

In the present case of this issue (matrix over GF(2**e)), the rank method is significantly more efficient than the determinant one, but not for a good (algorithmic) reason: this only seems to be due to the absence of a dedicated determinant method.

Okay that makes sense.

My complaint is that a lot of boilerplate is added, for example the docstring need to be mostly copied. I'd rather not doing that unless absolutely necessary.

So… I guess you can add a timing thing and make the test fail once determinant get fast? That or add the explanation why the method need to exist (link here).

If the benchmark sounds like too much work/too flaky, the latter would be good for me also.

[vneiger]

My complaint is that a lot of boilerplate is added, for example the docstring need to be mostly copied.

Yes, that is not nice. I'm not sure how to avoid this.

add the explanation why the method need to exist (link here).

Maybe that would be a good option (the other one you propose looks good as well but actually doing it looks trickier). I have seen such notes sometimes in docstrings.

[user202729]

Probably better to do #39912 first. Then we decide whether rank is still faster than determinant(algorithm=rref) later. (though I doubt)

Then we need to benchmark whether it makes sense to switch to algorithm=rref by default on other fields.

If that works, the docstring duplication etc. wouldn't be needed.

[whoami730]

This is changing behaviour: until now, for this type of matrices, one could call inverse on a non-square matrix without issue:

sage: mat = matrix.random(GF(4), 2, 4)
sage: mat.inverse()
[ 0 z2  0  0]
[ 1  0  0  0]
sage: mat
[     0      1 z2 + 1      0]
[z2 + 1      0 z2 + 1 z2 + 1]

(I would say that this behaviour change is actually a good thing, but this means at least checking that this is not breaking anything elsewhere in Sage.)

Honestly I would like to understand what does mat.inverse() actually compute here? It does not seem to be possible to compute mat * ~mat or ~mat * mat due to size issues already. Neither does it seem to be the moore penrose inverse.

[vneiger]

Honestly I would like to understand what does mat.inverse() actually compute here? It does not seem to be possible to compute mat * ~mat or ~mat * mat due to size issues already. Neither does it seem to be the moore penrose inverse.

Following the discussion above, you will end up at this question where this was discussed, and a related issue and fix (#30161). The fix with the rank check was not complete as it still "inverts" full row rank matrices... While these links explain well what is being computed, this behaviour should probably be changed.

[user202729]

might as well get this in, I don't see any other issue.

(apart from #39911 , but not sure what can be done about that.)