Opened 3 years ago
Closed 3 years ago
#28402 closed defect (fixed)
incorrect inverse of sparse matrix over inexact rings
Reported by:  Thierry Monteil  Owned by:  

Priority:  blocker  Milestone:  sage8.9 
Component:  linear algebra  Keywords:  
Cc:  Travis Scrimshaw, Vincent Delecroix  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 3 years ago by
Authors:  → Travis Scrimshaw 

Branch:  → public/linear_algebra/fix_sparse_inverse_inexact_fields28402 
Commit:  → 7d72279fb6595c4ad3aa9397b62dbb1866391987 
Status:  new → needs_review 
comment:2 Changed 3 years 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@zenbookflip:~$ 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 followup: 4 Changed 3 years 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 followup: 5 Changed 3 years 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 followup: 10 Changed 3 years 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 3 years 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 3 years ago by
Reviewers:  → Thierry Monteil 

Status:  needs_review → 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 3 years ago by
Commit:  7d72279fb6595c4ad3aa9397b62dbb1866391987 → 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 Changed 3 years 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 3 years 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 3 years 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 3 years ago by
Priority:  major → 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 3 years ago by
Status:  needs_review → positive_review 

comment:16 Changed 3 years ago by
Branch:  public/linear_algebra/fix_sparse_inverse_inexact_fields28402 → df1cbceb4157dcb5099a13bcbeaf7ca12aa4139c 

Resolution:  → fixed 
Status:  positive_review → 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.