incorrect inverse of sparse matrix over inexact rings

[sagetrac-tmonteil]

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.

CC: @tscrim @videlec

Component: linear algebra

Author: Travis Scrimshaw

Branch/Commit: df1cbce

Reviewer: Thierry Monteil

Issue created by migration from https://trac.sagemath.org/ticket/28402

[tscrim]

Branch: public/linear_algebra/fix_sparse_inverse_inexact_fields-28402

[tscrim]

Commit: 7d72279

[tscrim]

comment:1

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:

7d72279

Use the special sparse construction only over exact fields.

[tscrim]

Author: Travis Scrimshaw

[EmmanuelCharpentier]

comment:2

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 ?

[tscrim]

comment:3

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)

[EmmanuelCharpentier]

comment:4

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 : ... ]

HTH,

[EmmanuelCharpentier]

comment:5

Replying to @EmmanuelCharpentier:

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...

[tscrim]

comment:6

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().

[sagetrac-tmonteil]

comment:7

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.

[sagetrac-tmonteil]

Reviewer: Thierry Monteil

[sagetrac-git]

Changed commit from 7d72279 to df1cbce

[sagetrac-git]

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.

[tscrim]

comment:9

Good ideas. Done.

[mwageringel]

comment:10

Replying to @EmmanuelCharpentier:

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?

[mwageringel]

comment:11

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.

[tscrim]

comment:12

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.

[tscrim]

comment:13

Also, I think this should be a blocker because it is a fundamental computation that can silently returns a wrong answer.

[tscrim]

comment:15

Thank you.

[vbraun]

Changed branch from public/linear_algebra/fix_sparse_inverse_inexact_fields-28402 to df1cbce