Mpmath determinant calculation bug

[Sameer Grover]

I have a complex 2X2 matrix for which I'd like to calculate the
determinant. Consider the following code:

#========================
from __future__ import division
import numpy as np
import sympy.mpmath as mpmath

m3 = np.matrix([[ 5.18895807e-12 +0.00000000e+00j,
1.92716917e+11 +0.00000000e+00j],
[ 0.00000000e+00 -1.93643371e-11j,
0.00000000e+00 -5.16413237e+10j]])
m = mpmath.matrix(m3)
print np.linalg.det(m3)
print mpmath.det(m)
print m[0,0]*m[1,1]- m[1,0]*m[0,1]
#========================

The output is :
3.4638706823j
0
(0.0 + 3.46387068230212j)

mpmath's det function gives a determinant of zero even though the matrix
is non-singular. We get the correct value by multiplying term-by-term.

Can somebody confirm this bug?

Thank you.

Regards,
Sameer

[Ondřej Čertík]

Seems that I am getting the same as you are:

>>> from __future__ import division
>>> import numpy as np
>>> import sympy.mpmath as mpmath
>>>
>>> m3 = np.matrix([[ 5.18895807e-12 +0.00000000e+00j,

... 1.92716917e+11 +0.00000000e+00j],
... [ 0.00000000e+00 -1.93643371e-11j,
... 0.00000000e+00 -5.16413237e+10j]])

>>> m = mpmath.matrix(m3)
>>> print np.linalg.det(m3)

3.4638706823j
>>> print mpmath.det(m)
0

>>> print m[0,0]*m[1,1]- m[1,0]*m[0,1]

(0.0 + 3.46387068230212j)


I am using SymPy a35203f.

Ondrej

[Matthew Rocklin]

My guess is that it's a precision issue in mpmath. Perhaps it orders the computations in a way that causes cancellation. The range in magnitudes between 5e-12 1.9e11 is greater than the precision of a 64 bit floating point number (roughly 1e16). It's not surprising that software fails or that results vary in this regime.  

SymPy's matrix handles this fine by the way

In [14]: Matrix(m3).det()

Out[14]: 3.46387068230212⋅ⅈ

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To post to this group, send email to sy...@googlegroups.com.
To unsubscribe from this group, send email to sympy+un...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/sympy?hl=en.

[Nicolas Pettiaux]

Wouldn't it be interesting to see how other softwares handle the case
(maxima, mathematica, maple, reduce or numerically only octave, scilab
and matlab) ?

Nicolas

--
Nicolas Pettiaux

[Aaron Meurer]

This bug seems to relate specifically to mpmath, so I would ask about
it on the mpmath mailing list.

Aaron Meurer

[Sameer Grover]

Thanks for the input. I've reported it on the mpmath bug tracker
http://code.google.com/p/mpmath/issues/detail?id=229.

Sameer