Opened 18 months ago
Closed 18 months ago
#28402 closed defect (fixed)
incorrect inverse of sparse matrix over inexact rings
Reported by: | tmonteil | Owned by: | |
---|---|---|---|
Priority: | blocker | Milestone: | sage-8.9 |
Component: | linear algebra | Keywords: | |
Cc: | tscrim, vdelecroix | Merged in: | |
Authors: | Travis Scrimshaw | Reviewers: | Thierry Monteil |
Report Upstream: | N/A | Work issues: | |
Branch: | df1cbce (Commits, GitHub, GitLab) | Commit: | df1cbceb4157dcb5099a13bcbeaf7ca12aa4139c |
Dependencies: | Stopgaps: |
Description
As reported on this ask question, we have:
sage: B=matrix(RR, [[1/6, -1/24, -1/30, 1/120,1/12, 0, 0, 0, 0], [-1/24,1/60,1/60, 1/420, -1/24, 0, 0, 0, 0], [-1/30,1/60, 2/105, 1/140, -1/20, 0, 0, 0, 0], [1/120, 1/420, 1/140, 13/1260, -1/40, 0, 0, 0, 0 ....: ], [1/12, -1/24, -1/20, -1/40, 1/3, -1/24, -1/30, 1/120,1/12], [0, 0, 0, 0, -1/24,1/60,1/60, 1/420, -1/24], [0, 0, 0, 0, -1/30,1/60, 2/105, 1/140, -1/20], [0, 0, 0, 0, 1/120, 1/420, 1/140, 13/1260, - ....: 1/40], [0, 0, 0, 0,1/12, -1/24, -1/20, -1/40, 1/6]],sparse=true) sage: (B.inverse()*B).norm(1) 138.4999999999923
Note that reverting #24122, by removing the call to build_inverse_from_augmented_sparse
in the __invert__
method of sage/matrix/matrix0.pyx
(line 5410), leads to a correct answer.
Change History (16)
comment:1 Changed 18 months ago by
- Branch set to public/linear_algebra/fix_sparse_inverse_inexact_fields-28402
- Commit set to 7d72279fb6595c4ad3aa9397b62dbb1866391987
- Status changed from new to needs_review
comment:2 Changed 18 months ago by
Can't reproduce:
sage: sage.version.version '8.9.beta7' sage: B=matrix(RR, [[1/6, -1/24, -1/30, 1/120,1/12, 0, 0, 0, 0], [-1/24,1/60,1/60 ....: , 1/420, -1/24, 0, 0, 0, 0], [-1/30,1/60, 2/105, 1/140, -1/20, 0, 0, 0, 0], ....: [1/120, 1/420, 1/140, 13/1260, -1/40, 0, 0, 0, 0], [1/12, -1/24, -1/20, -1 ....: /40, 1/3, -1/24, -1/30, 1/120,1/12], [0, 0, 0, 0, -1/24,1/60,1/60, 1/420, - ....: 1/24], [0, 0, 0, 0, -1/30,1/60, 2/105, 1/140, -1/20], [0, 0, 0, 0, 1/120, 1 ....: /420, 1/140, 13/1260, -1/40], [0, 0, 0, 0,1/12, -1/24, -1/20, -1/40, 1/6]], ....: sparse=true) sage: (B.inverse()*B).norm(1) 1.0000000000019988
No relevant packages installed:
charpent@zen-book-flip:~$ sage -optional | grep -v ot_instal [package]...............................[latest version] ([version]) dot2tex.................................2.11.3.p0 (2.11.3.p0) fricas..................................1.3.5 (1.3.5) gap_packages............................4.10.2.p1 (4.10.2.p1) giacpy_sage.............................0.6.7 (0.6.7) libsemigroups...........................0.6.7 (0.6.7) python2.................................2.7.15.p1 (2.7.15.p1)
Do you use any relevant packages ?
comment:3 follow-up: ↓ 4 Changed 18 months ago by
Maybe it is computer dependent? I am able to reproduce it on 8.9.beta8
(previously on beta7):
sage: (B.inverse()*B).norm(1) 138.4999999999923 sage: B.det() -0.000000000000000
uqtscrim@SMP-36PQ8T2:~/sage-build$ ./sage -optional | grep -v ot_install [package]...............................[latest version] ([version]) 4ti2....................................1.6.7.p0 (1.6.7.p0) bliss...................................0.73+debian-1+sage-2016-08-02.p0 (0.73+debian-1+sage-2016-08-02.p0) cmake...................................3.11.0 (3.11.0) coxeter3................................8ac9c71723c8ca57a836d6381aed125261e44e9e (8ac9c71723c8ca57a836d6381aed125261e44e9e) dot2tex.................................2.11.3.p0 (2.11.3.p0) e_antic.................................0.1.3 (0.1.3) fricas..................................1.3.5 (1.3.5) frobby..................................0.9.0.p2 (0.9.0.p2) gambit..................................15.1.1.p0 (15.1.1.p0) latte_int...............................1.7.5.p0 (1.7.5.p0) libsemigroups...........................0.6.7 (0.6.7) lidia...................................2.3.0+latte-patches-2019-05-02 (2.3.0+latte-patches-2019-05-02) lrslib..................................062+autotools-2017-03-03.p1 (062+autotools-2017-03-03.p1) meataxe.................................1.0.p0 (1.0.p0) mpir....................................3.0.0-644faf502c56f97d9accd301965fc57d6ec70868.p0 (3.0.0-644faf502c56f97d9accd301965fc57d6ec70868.p0) ninja_build.............................1.8.2 (1.8.2) normaliz................................3.7.2 (3.7.2) ore_algebra.............................0.3 (0.3) p_group_cohomology......................3.2 (3.2) pynormaliz..............................2.7 (2.7) python2.................................2.7.15.p1 (2.7.15.p1) qhull...................................2015-src-7.2.0.p1 (2015-src-7.2.0.p1) sirocco.................................2.0.p0 (2.0.p0) tides...................................2.0.p0 (2.0.p0) topcom..................................0.17.7 (0.17.7)
comment:4 in reply to: ↑ 3 ; follow-up: ↓ 5 Changed 18 months ago by
Replying to tscrim:
Maybe it is computer dependent?
Debian testing running on core i7 + 16 GB RAM here.
On another machine, slightly smalle (Debian testing, core i5 + 8 GB RAM):
sage: (B.inverse()*B).norm(1) 1.0000000000019988 sage: B.det() -0.000000000000000
I am able to reproduce it on
8.9.beta8
(previously on beta7):
I'll upgrade to beta8 and report results.
[ Bandwidth savings : <snip>... ]
HTH,
comment:5 in reply to: ↑ 4 ; follow-up: ↓ 10 Changed 18 months ago by
Replying to charpent:
I'll upgrade to beta8 and report results.
Same results :
sage: %cpaste Pasting code; enter '--' alone on the line to stop or use Ctrl-D. :sage: B=matrix(RR, [[1/6, -1/24, -1/30, 1/120,1/12, 0, 0, 0, 0], [-1/24,1/60,1/60, 1/420, -1/24, 0, 0, 0, 0], [-1/30,1/60, 2/105, 1/140, -1/20, 0, 0, 0, 0], [1/120, 1/420, 1/140, 13/1260, -1/40, 0, 0, 0, 0 :....: ], [1/12, -1/24, -1/20, -1/40, 1/3, -1/24, -1/30, 1/120,1/12], [0, 0, 0, 0, -1/24,1/60,1/60, 1/420, -1/24], [0, 0, 0, 0, -1/30,1/60, 2/105, 1/140, -1/20], [0, 0, 0, 0, 1/120, 1/420, 1/140, 13/1260, - :....: 1/40], [0, 0, 0, 0,1/12, -1/24, -1/20, -1/40, 1/6]],sparse=true) :sage: (B.inverse()*B).norm(1) :-- 1.0000000000019988 sage: B.det() -0.000000000000000
I'm using a Python 3-based Sage 8.9.beta8, but I can't see how it would make a difference...
comment:6 Changed 18 months ago by
I would argue this change still makes sense for the reasons I gave in comment:1. It is assuming that after doing the rref we obtain [I|Ainv]
and there are no non-zero entries in I
except for the diagonal, which I would say is dangerous unless we implement a sparse matrix specific echelonization procedure that guarantees such entries are exactly 0 and do not appear in the self.dict()
.
comment:7 Changed 18 months ago by
- Reviewers set to Thierry Monteil
- Status changed from needs_review to needs_work
Thanks for the fix. I would suggest to:
- change the tolerance in the doctest: with the current tolerance of
1e-14
, the expected value1.0
is not in the interval, which might be confusing. Perhaps could you replace it with2e-12
(and maybe even replace the expected answer to 1.0). - to test the
build_inverse_from_augmented_sparse
part, e.g. by adding a similar doctest replacingRR
withQQ
.
comment:8 Changed 18 months ago by
- Commit changed from 7d72279fb6595c4ad3aa9397b62dbb1866391987 to df1cbceb4157dcb5099a13bcbeaf7ca12aa4139c
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
df1cbce | Use the special sparse construction only over exact fields.
|
comment:10 in reply to: ↑ 5 Changed 18 months ago by
Replying to charpent:
I'm using a Python 3-based Sage 8.9.beta8, but I can't see how it would make a difference...
I also cannot replicate the issue with Python 3, whereas I can using Python 2, on those same machines. Any idea what could be causing this?
comment:11 Changed 18 months ago by
Actually, having had a look at the implementation, the reason why this seemingly works with Python 3 is that iterating over dicts is done in an ordered fashion, whereas it is done in arbitrary order with Python 2. Thus, the extraneous non-zero entries in the left half of the augmented matrix are consistently overwritten by entries on the right if they exist.
All in all, I agree this is the correct solution.
comment:12 Changed 18 months ago by
Thank you for investigating it. That would also explain why it works for some people in Python 2 as well.
Patchbot is (morally) green as well.
comment:13 Changed 18 months ago by
- Priority changed from major to blocker
Also, I think this should be a blocker because it is a fundamental computation that can silently returns a wrong answer.
comment:14 Changed 18 months ago by
- Status changed from needs_review to positive_review
comment:15 Changed 18 months ago by
Thank you.
comment:16 Changed 18 months ago by
- Branch changed from public/linear_algebra/fix_sparse_inverse_inexact_fields-28402 to df1cbceb4157dcb5099a13bcbeaf7ca12aa4139c
- Resolution set to fixed
- Status changed from positive_review to closed
Okay, so I figured out the problem. Basically it comes from the inexactness of the arithmetic. There are a lot of non-zero entries in the part of the augmented matrix that are normally 0 but are not when working over inexact fields. So I just moved that special construction to only be for exact fields as the resulting matrix is much more likely to be dense anyways for inexact fields.
New commits:
Use the special sparse construction only over exact fields.