Opened 9 months ago

Closed 9 months ago

Last modified 9 months ago

#28570 closed defect (fixed)

fix matrix inversion over ZZ

Reported by: chapoton Owned by:
Priority: critical Milestone: sage-9.0
Component: linear algebra Keywords:
Cc: vdelecroix, jhpalmieri Merged in:
Authors: Frédéric Chapoton Reviewers: Vincent Delecroix, John Palmieri
Report Upstream: N/A Work issues:
Branch: 56ff2bb (Commits) Commit:
Dependencies: Stopgaps:

Description

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()

Change History (22)

comment:1 Changed 9 months ago by 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.

comment:2 Changed 9 months ago by 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.

comment:3 follow-up: Changed 9 months ago by 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'
Last edited 9 months ago by chapoton (previous) (diff)

comment:4 in reply to: ↑ 3 ; follow-up: Changed 9 months ago by 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.

comment:5 Changed 9 months ago by chapoton

  • Branch set to public/ticket/28570
  • Commit set to 8ef9966171213db2b23bc87bd7e9e1a824947566

Trying to fix one of the problems.


New commits:

8ef9966trac 28570 fix one of the problems

comment:6 follow-up: Changed 9 months ago by 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

comment:7 in reply to: ↑ 6 Changed 9 months ago by 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.

comment:8 Changed 9 months ago by chapoton

_adjugate does that :

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

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

comment:9 Changed 9 months ago by git

  • Commit changed from 8ef9966171213db2b23bc87bd7e9e1a824947566 to 8c5962392dc769db108d4e0c3d29400ee42286a4

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

8c59623trac 28570 fixing the slowness problem

comment:10 Changed 9 months ago by chapoton

  • Authors set to Frédéric Chapoton
  • Status changed from new to needs_review

Here is a proposal. Please review

comment:11 Changed 9 months ago by chapoton

green bot

comment:12 Changed 9 months ago by 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

comment:13 Changed 9 months ago by chapoton

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

comment:14 Changed 9 months ago by chapoton

  • Authors Frédéric Chapoton deleted
  • Status changed from needs_review to needs_work

comment:15 Changed 9 months ago by git

  • Commit changed from 8c5962392dc769db108d4e0c3d29400ee42286a4 to 56ff2bbdc8e0e83f6128aee294ee2e07a334932f

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

56ff2bbtrac 28570 proper handling of Mat(ZZ)

comment:16 Changed 9 months ago by chapoton

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

comment:17 Changed 9 months ago by chapoton

  • Authors set to Frédéric Chapoton
  • Status changed from needs_work to needs_review

comment:18 Changed 9 months ago by 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.

comment:19 Changed 9 months ago by vdelecroix

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

Last edited 9 months ago by vdelecroix (previous) (diff)

comment:20 Changed 9 months ago by jhpalmieri

  • Reviewers set to Vincent Delecroix, John Palmieri
  • Status changed from needs_review to positive_review

comment:21 Changed 9 months ago by vbraun

  • Branch changed from public/ticket/28570 to 56ff2bbdc8e0e83f6128aee294ee2e07a334932f
  • Resolution set to fixed
  • Status changed from positive_review to closed

comment:22 in reply to: ↑ 4 Changed 9 months ago by jhpalmieri

  • Commit 56ff2bbdc8e0e83f6128aee294ee2e07a334932f deleted

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.)

Note: See TracTickets for help on using tickets.