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:  sage8.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_fields28402
 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
comment:3 followup: ↓ 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@SMP36PQ8T2:~/sagebuild$ ./sage optional  grep v ot_install [package]...............................[latest version] ([version]) 4ti2....................................1.6.7.p0 (1.6.7.p0) bliss...................................0.73+debian1+sage20160802.p0 (0.73+debian1+sage20160802.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+lattepatches20190502 (2.3.0+lattepatches20190502) lrslib..................................062+autotools20170303.p1 (062+autotools20170303.p1) meataxe.................................1.0.p0 (1.0.p0) mpir....................................3.0.0644faf502c56f97d9accd301965fc57d6ec70868.p0 (3.0.0644faf502c56f97d9accd301965fc57d6ec70868.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...................................2015src7.2.0.p1 (2015src7.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 ; followup: ↓ 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
EDIT :
sage: P=sage.misc.package.list_packages('optional') sage: P.update(sage.misc.package.list_packages('experimental')) sage: [P.get(u).get('name') for u in P.keys() if P.get(u).get('installed')] ['fricas', 'dot2tex', 'giacpy_sage', 'python2']
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 ; followup: ↓ 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 CtrlD. :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 3based 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 [IAinv]
and there are no nonzero 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
1e14
, the expected value1.0
is not in the interval, which might be confusing. Perhaps could you replace it with2e12
(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 3based 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 nonzero 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_fields28402 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 nonzero 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.