# Inverses of matrices over RR

47 views

### Håkan Granath

Nov 5, 2022, 2:59:27 AM11/5/22
to sage-devel
Hi, there seems to be a problem with inverses of matrices with elements in RR. It only occurs very sporadically for me, but here is an example:

a = RR(-4967757600021511 / 2**106)
b = RR(-7769080564883485 / 2**52)
c = RR( 5221315298319565 / 2**53)

m = matrix([[a, b], [c, -a]])

print(m)
print()
print(~m)

On my machines it produces the output

[-6.12323399573677e-17     -1.72508242466029]
[    0.579682446302195  6.12323399573677e-17]

[     4.00000000000000      1.72508242466029]
[   -0.579682446302195 -6.12323399573676e-17]

Clearly the element 4 is wrong (the correct inverse is -m). Is this a known bug?

Some system information:

SageMath version 9.7, using Python 3.10.5
OS: Ubuntu 20.04.5 LTS
CPU: Intel(R) Core(TM) i7-7700 CPU @ 3.60GHz

Best regards,

Håkan Granath

### Emmanuel Charpentier

Nov 5, 2022, 6:20:58 AM11/5/22
to sage-devel

something is definitely unhinged here : On 9.8.beta3 running on Debian testing on core i7 + 16 GB RAM, after running :

``````a = RR(-4967757600021511 / 2**106)
b = RR(-7769080564883485 / 2**52)
c = RR( 5221315298319565 / 2**53)
m = matrix([[a, b], [c, -a]])
M = matrix([[var("p%d%d"%(u, v), latex_name="p_{%s,%d}"%(u, v))
for v in range(2)]
for u in range(2)])
S = dict(zip(M.list(), [a, b, c, -a]))
MN = M.apply_map(lambda u:u.subs(S))
``````

one gets :

``````sage: m.parent()
Full MatrixSpace of 2 by 2 dense matrices over Real Field with 53 bits of precision
sage: m*~m
[     1.00000000000000 -1.23259516440783e-32]
[     2.31872978520878      1.00000000000000]
sage: (~m)*m
[    1.00000000000000    -6.90032969864117]
[6.16297582203915e-33     1.00000000000000]
sage: MN.parent()
Full MatrixSpace of 2 by 2 dense matrices over Symbolic Ring
sage: MN*~MN
[     1.00000000000000 -1.23259516440783e-32]
[     2.31872978520878      1.00000000000000]
sage: (~MN)*MN
[    1.00000000000000    -6.90032969864117]
[6.16297582203915e-33     1.00000000000000]
``````

all being wrong, wrong, wrong

However :

``````sage: (M*~M).apply_map(lambda u:u.subs(S))
[     1.00000000000000                     0]
[-3.54953126192945e-17      1.00000000000000]
sage: ((~M)*M).apply_map(lambda u:u.subs(S))
[    1.00000000000000 1.05630833481279e-16]
[                   0     1.00000000000000]
``````

both being acceptable.

One also notes that the form of :

``````sage: ~M
[1/p00 - p01*p10/(p00^2*(p01*p10/p00 - p11))               p01/(p00*(p01*p10/p00 - p11))]
[              p10/(p00*(p01*p10/p00 - p11))                      -1/(p01*p10/p00 - p11)]
sage: (~M).apply_map(simplify)
[1/p00 - p01*p10/(p00^2*(p01*p10/p00 - p11))               p01/(p00*(p01*p10/p00 - p11))]
[              p10/(p00*(p01*p10/p00 - p11))                      -1/(p01*p10/p00 - p11)]
``````

is somewhat unexpected ; one expects :

``````sage: (~M).apply_map(lambda u:u.simplify_full())
[-p11/(p01*p10 - p00*p11)  p01/(p01*p10 - p00*p11)]
[ p10/(p01*p10 - p00*p11) -p00/(p01*p10 - p00*p11)]
``````

which is also the form returned by `maxima` :

``````sage: maxima_calculus.invert(M).sage()
[-p11/(p01*p10 - p00*p11)  p01/(p01*p10 - p00*p11)]
[ p10/(p01*p10 - p00*p11) -p00/(p01*p10 - p00*p11)]
``````

giac :

``````sage: giac.inverse(giac(M)).sage()
[[-p11/(p01*p10 - p00*p11), p01/(p01*p10 - p00*p11)],
[p10/(p01*p10 - p00*p11), -p00/(p01*p10 - p00*p11)]]
``````

fricas :

``````sage: fricas.inverse(M._fricas_()).sage()
[-p11/(p01*p10 - p00*p11)  p01/(p01*p10 - p00*p11)]
[ p10/(p01*p10 - p00*p11) -p00/(p01*p10 - p00*p11)]
``````

mathematica :

``````sage: mathematica.Inverse(M).sage()
[[-p11/(p01*p10 - p00*p11), p01/(p01*p10 - p00*p11)],
[p10/(p01*p10 - p00*p11), -p00/(p01*p10 - p00*p11)]]
``````

and (somewhat un-backconvertible) :

``````sage: sympy.sympify(M)^-1._sage_()
Matrix([
[ p11/(p00*p11 - p01*p10), -p01/(p00*p11 - p01*p10)],
[-p10/(p00*p11 - p01*p10),  p00/(p00*p11 - p01*p10)]])
``````

This is, IMNSHO, a critical bug. Could you open a tichet for this, and mark it as such ?

### Dima Pasechnik

Nov 5, 2022, 7:08:40 AM11/5/22
Well, applying a naive exact linear algebra routine to inexact data,
and that's what Sage is doing here, is prone to errors.
sage: A=m.augment(identity_matrix(RR,2))
sage: A
[-6.12323399573677e-17 -1.72508242466029 1.00000000000000
0.000000000000000]
[ 0.579682446302195 6.12323399573677e-17 0.000000000000000
1.00000000000000]
sage: A.echelonize(algorithm="classical");A # OOOPS! - that's the default here
[ 1.00000000000000 0.000000000000000 4.00000000000000
1.72508242466029]
[ 0.000000000000000 1.00000000000000 -0.579682446302195
-6.12323399573676e-17]
sage: A=m.augment(identity_matrix(RR,2))
sage: A.echelonize(algorithm='scaled_partial_pivoting');A # that's how
it should be
[ 1.00000000000000 0.000000000000000 6.12323399573677e-17
1.72508242466029]
[ 0.000000000000000 1.00000000000000 -0.579682446302195
-6.12323399573677e-17]
> --
> You received this message because you are subscribed to the Google Groups "sage-devel" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+...@googlegroups.com.

### Dima Pasechnik

Nov 5, 2022, 7:15:26 AM11/5/22