Opened 4 months ago

Closed 3 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) 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 4 months ago by tscrim

  • Authors set to Travis Scrimshaw
  • Branch set to public/linear_algebra/fix_sparse_inverse_inexact_fields-28402
  • Commit set to 7d72279fb6595c4ad3aa9397b62dbb1866391987
  • Status changed from new to needs_review

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:

7d72279Use the special sparse construction only over exact fields.

comment:2 Changed 4 months ago by charpent

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 ?

Last edited 4 months ago by charpent (previous) (diff)

comment:3 follow-up: Changed 4 months ago by tscrim

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: Changed 4 months ago by charpent

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,

Last edited 4 months ago by charpent (previous) (diff)

comment:5 in reply to: ↑ 4 ; follow-up: Changed 4 months ago by charpent

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 4 months ago by tscrim

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 4 months ago by tmonteil

  • 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 value 1.0 is not in the interval, which might be confusing. Perhaps could you replace it with 2e-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 replacing RR with QQ.

comment:8 Changed 4 months ago by git

  • Commit changed from 7d72279fb6595c4ad3aa9397b62dbb1866391987 to df1cbceb4157dcb5099a13bcbeaf7ca12aa4139c

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

df1cbceUse the special sparse construction only over exact fields.

comment:9 Changed 4 months ago by tscrim

  • Status changed from needs_work to needs_review

Good ideas. Done.

comment:10 in reply to: ↑ 5 Changed 3 months ago by gh-mwageringel

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 3 months ago by gh-mwageringel

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 3 months ago by tscrim

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 months ago by tscrim

  • 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 3 months ago by tmonteil

  • Status changed from needs_review to positive_review

comment:15 Changed 3 months ago by tscrim

Thank you.

comment:16 Changed 3 months ago by vbraun

  • 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
Note: See TracTickets for help on using tickets.