Completely confused about simple Sage matrix operation error

[David Ingerman]

The following matrix operation produces wrong answer in online Sage:
M=matrix(RR,[[7,3,10,13],[1,1,2,2],[1,2,3,4],[1,3,5,7]]);det(M);invM=M^(-1);invM*M;det(invM)
If one changes RR to QQ the answers turn correct. Or it is enough to change one corner matrix entry to:
M=matrix(RR,[[8,3,10,13],[1,1,2,2],[1,2,3,4],[1,3,5,7]]);det(M);invM=M^(-1);invM*M;det(invM)
to get right answers.
How can one be sure about Sage linear algebra w/this problem? Am I missing something obvious? It took me several hours to zero in on this bug from a large program. Ant feedback will be appreciated.

[David Ingerman]

'Sage Version 5.4, Release Date: 2012-11-09'

[Harald Schilly]

On Monday, July 1, 2013 10:45:36 AM UTC+2, David Ingerman wrote:

The following matrix operation produces wrong answer in online Sage:
M=matrix(RR,[[7,3,10,13],[1,1,2,2],[1,2,3,4],[1,3,5,7]]);det(M);invM=M^(-1);invM*M;det(invM)

RR stands for the "real numbers" with the usual 53bits of precision, e.g. 5.123957322….
QQ are rational numbers where two large integers build up each number,
e.g. 41000041/333000333000333000333000333000333000333
Therefore, QQ has a much higher precision … but is much slower and uses more memory.

In your case, the lack of precision in RR causes you troubles and you have to find a way to pose the problem you want to solve differently. You cannot rely on QQ, because in bad cases, the expressions blow up and eat all your memory. More generally, this is not a Sage related problem, but related to all calculations your are doing "natively" with your CPU.

To see in advance when this happens, you have to calculate the conditional number of the matrix. I think that's only in numpy (or I haven't found it).

sage: M=matrix(RR,[[7,3,10,13],[1,1,2,2],[1,2,3,4],[1,3,5,7]])
sage: import numpy as np
sage: np.linalg.cond(M)
104.85355762315329

http://en.wikipedia.org/wiki/Condition_number

Here are some decomposition methods that might help:

http://en.wikipedia.org/wiki/Matrix_decomposition

H

[David Ingerman]

But this is a very small matrix and conditional number is not large...

[David Joyner]

Using RDF is better:

sage: M=matrix(RDF,[[7,3,10,13],[1,1,2,2],[1,2,3,4],[1,3,5,7]]);det(M);invM=M^(-1);invM*M;det(invM)
-7.0
[ 1.0 0.0 0.0 -8.881784197e-16]
[ 0.0 1.0 0.0 -8.881784197e-16]
[ 4.4408920985e-16 0.0 1.0 8.881784197e-16]
[-5.55111512313e-17 -2.77555756156e-16 -1.66533453694e-16 1.0]
-0.142857142857

>
> --
> You received this message because you are subscribed to the Google Groups "sage-support" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sage-support...@googlegroups.com.
> To post to this group, send email to sage-s...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sage-support.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

[Mike Hansen]

I think problem is actually due to the inverse using a non-numerically
stable echelon form algorithm for inexact fields. For example, if you
using matrices over RDF:

M=matrix(RDF,[[7,3,10,13],[1,1,2,2],[1,2,3,4],[1,3,5,7]]);det(M);invM=M^(-1);invM*M;det(invM)

you don't get this problem. My guess is that if we implemented
something like partial pivoting for inexact fields.

--Mike
--Mike