fix matrix inversion over ZZ

[chapoton]

This was broken by #28189.

Probably due to the new "inverse_of_unit" method taking over the existing one.

This cause the following to never stop:

sage: P = posets.TamariLattice(7)
sage: H = P._hasse_diagram
sage: M = H._leq_matrix
sage: M.inverse_of_unit()

[nbruin]

Possibly relevant point: Prior to #28189, this routine returned a sparse matrix. If sparsity isn't taken into account in the new code that could very well explain a drastic slow-down.

[jhpalmieri]

First, the command M.inverse_of_unit() does actually finish for me, eventually. It takes 5 minutes or so, in contrast to ~M which takes under 1 second.

From ​reading the old code, it looks like M.inverse_of_unit() just computed the inverse and then coerced the inverse matrix to have the same base ring. In the new version, it tries to be more clever, in particular calling M.adjugate(). I think M.adjugate() is the problem for sparse integer matrices, or indeed for integer matrices in general. Even in 8.9, running M.adjugate() or M.dense_matrix().adjugate() on the matrix in the ticket description takes a long time.

One way around this would be to add an elif clause to the new inverse_of_unit which handles the case when the base ring is the integers. I don't know if this is the best choice.

[chapoton]

and we also have now the following not-so-funny bug:

sage: m = matrix(Zmod(13**3),2,2,[3,5,7,11],sparse=True)
sage: m.inverse()     
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-12-91aa9265ee1c> in <module>()
----> 1 m.inverse()

/home/chapoton/sage3/local/lib/python3.7/site-packages/sage/matrix/matrix2.pyx in sage.matrix.matrix2.Matrix.inverse (build/cythonized/sage/matrix/matrix2.c:66498)()
   8919 
   8920         """
-> 8921         return ~self
   8922 
   8923     def adjugate(self):

/home/chapoton/sage3/local/lib/python3.7/site-packages/sage/matrix/matrix0.pyx in sage.matrix.matrix0.Matrix.__invert__ (build/cythonized/sage/matrix/matrix0.c:35802)()
   5400                 return ~self.matrix_over_field()
   5401             else:
-> 5402                 return self.inverse_of_unit()
   5403         else:
   5404             A = self.augment(self.parent().identity_matrix())

/home/chapoton/sage3/local/lib/python3.7/site-packages/sage/matrix/matrix0.pyx in sage.matrix.matrix0.Matrix.inverse_of_unit (build/cythonized/sage/matrix/matrix0.c:36801)()
   5501             # CRT or p-adic lifting.
   5502             N = R.characteristic()
-> 5503             m, D = self.lift_centered()._invert_flint()
   5504             if not D.gcd(N).is_one():
   5505                 raise ZeroDivisionError("input matrix must be nonsingular")

/home/chapoton/sage3/local/lib/python3.7/site-packages/sage/structure/element.pyx in sage.structure.element.Element.__getattr__ (build/cythonized/sage/structure/element.c:4607)()
    487             AttributeError: 'LeftZeroSemigroup_with_category.element_class' object has no attribute 'blah_blah'
    488         """
--> 489         return self.getattr_from_category(name)
    490 
    491     cdef getattr_from_category(self, name):

/home/chapoton/sage3/local/lib/python3.7/site-packages/sage/structure/element.pyx in sage.structure.element.Element.getattr_from_category (build/cythonized/sage/structure/element.c:4716)()
    500         else:
    501             cls = P._abstract_element_class
--> 502         return getattr_from_other_class(self, cls, name)
    503 
    504     def __dir__(self):

/home/chapoton/sage3/local/lib/python3.7/site-packages/sage/cpython/getattr.pyx in sage.cpython.getattr.getattr_from_other_class (build/cythonized/sage/cpython/getattr.c:2614)()
    392         dummy_error_message.cls = type(self)
    393         dummy_error_message.name = name
--> 394         raise AttributeError(dummy_error_message)
    395     attribute = <object>attr
    396     # Check for a descriptor (__get__ in Python)

AttributeError: 'sage.matrix.matrix_integer_sparse.Matrix_integer_sparse' object has no attribute '_invert_flint'

[nbruin]

Replying to chapoton:

and we also have now the following not-so-funny bug

One would almost think the reviewer on that ticket could have done a more thorough job :-)

More seriously, I think this shows we are running into weaknesses of our review system. I don't have concrete suggestions on how to improve it, but if this kind of thing happens more often, it might be worth considering if there are procedural changes we can make that reduce the likelihood of such bugs slipping through, without unduly narrowing the bottleneck that review already is.

[chapoton]

Trying to fix one of the problems.

---

New commits:

​8ef9966

trac 28570 fix one of the problems

[chapoton]

Not related, but also quite bad:

sage: m = matrix(Zmod(2**2),1,1,[1],sparse=True)
sage: m.det()
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-10-6b1d86e9fa93> in <module>()
----> 1 m.det()

/home/chapoton/sage3/local/lib/python3.7/site-packages/sage/matrix/matrix2.pyx in sage.matrix.matrix2.Matrix.det (build/cythonized/sage/matrix/matrix2.c:16941)()
   1918             6
   1919         """
-> 1920         return self.determinant(*args, **kwds)
   1921 
   1922     def apply_morphism(self, phi):

/home/chapoton/sage3/local/lib/python3.7/site-packages/sage/matrix/matrix_modn_sparse.pyx in sage.matrix.matrix_modn_sparse.Matrix_modn_sparse.determinant (build/cythonized/sage/matrix/matrix_modn_sparse.cpp:9560)()
    844 
    845         if algorithm is None or algorithm == "linbox":
--> 846             r, d = self._rank_det_linbox()
    847             self.cache('rank', r)
    848             self.cache('det', d)

/home/chapoton/sage3/local/lib/python3.7/site-packages/sage/matrix/matrix_modn_sparse.pyx in sage.matrix.matrix_modn_sparse.Matrix_modn_sparse._rank_det_linbox (build/cythonized/sage/matrix/matrix_modn_sparse.cpp:8655)()
    682 
    683         if not is_prime(self.p):
--> 684             raise TypeError("only GF(p) supported via LinBox")
    685 
    686         cdef givaro.Modular_uint64 * F = new givaro.Modular_uint64(<uint64_t> self.p)

TypeError: only GF(p) supported via LinBox

[jhpalmieri]

Replying to chapoton:

Not related, but also quite bad:

This one is also not new: it's been present at least since Sage 8.6.

[chapoton]

_adjugate does that :

A = self.charpoly().shift(-1)(self)

which seems to be a rather terrible idea for matrices of large size.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​8c59623

trac 28570 fixing the slowness problem

[chapoton]

Here is a proposal. Please review

[chapoton]

green bot

[vdelecroix]

I am not sure this is a smart move. Your proposed solution performs a modular inversion for each coefficient instead of a single one. Here is a use case where one can identify a x2 slowdown.

sage: def inv1(m):
....:     p = m.base_ring().characteristic()           
....:     minv, d = m.lift_centered()._invert_flint()
....:     return d.inverse_mod(p) * minv.change_ring(R)
sage: def inv2(m):
....:     return (~self.lift_centered()).change_ring(R)
sage: p = next_prime(2**50)
sage: R = GF(p)
sage: m = random_matrix(R, 50)   # assuming we get something invertible...
sage: %timeit minv = inv1(m)
10 loops, best of 3: 93.5 ms per loop
sage: %timeit minv = inv2(m)
10 loops, best of 3: 174 ms per loop

[chapoton]

Hello Vincent. Feel free to propose a much better fix.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​56ff2bb

trac 28570 proper handling of Mat(ZZ)

[chapoton]

The case of fields is in fact handled by the previous part of the code.

[jhpalmieri]

I ran an experiment:

sage: p = next_prime(2**50)
sage: R = GF(p)
sage: %timeit random_matrix(R, 50).inverse_of_unit()

Without this branch:

10 loops, best of 5: 186 ms per loop

With this branch:

10 loops, best of 5: 185 ms per loop

As Frédéric says, fields are handled elsewhere, so of course there is no slowdown with this branch.

Are there other concerns? It certainly fixes the integer slowdown and the Z/p**n bug from comment:3.

[vdelecroix]

Looks good to me. Thanks for fixing my stupid mistake from #28189!

[jhpalmieri]

Replying to nbruin:

More seriously, I think this shows we are running into weaknesses of our review system. I don't have concrete suggestions on how to improve it, but if this kind of thing happens more often, it might be worth considering if there are procedural changes we can make that reduce the likelihood of such bugs slipping through, without unduly narrowing the bottleneck that review already is.

Here's an idea, but I don't know how feasible it is: extract every test in the matrix directory and make sure it works with both sparse and dense matrices. (Maybe not every test: maybe you don't want to mess with the ones in *_dense.*, etc.)