#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, GitHub, GitLab) | 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 3 years ago by
comment:2 Changed 3 years ago by
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: ↓ 4 Changed 3 years ago by
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'
comment:4 in reply to: ↑ 3 ; follow-up: ↓ 22 Changed 3 years ago by
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 3 years ago by
- Branch set to public/ticket/28570
- Commit set to 8ef9966171213db2b23bc87bd7e9e1a824947566
comment:6 follow-up: ↓ 7 Changed 3 years ago by
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 3 years ago by
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 3 years ago by
_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 3 years ago by
- Commit changed from 8ef9966171213db2b23bc87bd7e9e1a824947566 to 8c5962392dc769db108d4e0c3d29400ee42286a4
Branch pushed to git repo; I updated commit sha1. New commits:
8c59623 | trac 28570 fixing the slowness problem
|
comment:10 Changed 3 years ago by
- Status changed from new to needs_review
Here is a proposal. Please review
comment:11 Changed 3 years ago by
green bot
comment:12 Changed 3 years ago by
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 3 years ago by
Hello Vincent. Feel free to propose a much better fix.
comment:14 Changed 3 years ago by
- Status changed from needs_review to needs_work
comment:15 Changed 3 years ago by
- Commit changed from 8c5962392dc769db108d4e0c3d29400ee42286a4 to 56ff2bbdc8e0e83f6128aee294ee2e07a334932f
Branch pushed to git repo; I updated commit sha1. New commits:
56ff2bb | trac 28570 proper handling of Mat(ZZ)
|
comment:16 Changed 3 years ago by
The case of fields is in fact handled by the previous part of the code.
comment:17 Changed 3 years ago by
- Status changed from needs_work to needs_review
comment:18 Changed 3 years ago by
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 3 years ago by
Looks good to me. Thanks for fixing my stupid mistake from #28189!
comment:20 Changed 3 years ago by
- Reviewers set to Vincent Delecroix, John Palmieri
- Status changed from needs_review to positive_review
comment:21 Changed 3 years ago by
- 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 3 years ago by
- 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.)
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.