0

Is it possible to do higher precision matrix exponential in Python? I mean obtain higher precision than double floating numbers.

I have following testing code:

import sympy
from sympy import N
import random

n = 100
#A = sympy.Matrix([[random.random(),random.random()],
#                   [random.random(),random.random()]])
A = sympy.Matrix([[1,2],[3,4]])
dlt = 1000
e1 = A.exp()
e1 = N(e1, n)
ee2 = (A/dlt).exp()
ee2 = N(ee2, n)
e2 = sympy.eye(2)
for i in range(dlt):
    e2 = e2*ee2
print(N(max(e1-e2)))

Theoretically, the final result should be zero. With scipy, the error is about 1e-14.

By sympy, if the matrix is like [[1,2],[3,4]], the output of previous code is about 1e-98. However, for random matrix, the error is around 1e-14. Is it possible to get results like 1e-100 for random matrices?

Speed is not concern.

iuradz
  • 1,128
  • 1
  • 12
  • 25

2 Answers2

0

Once you use N, you are in the realm of floating point operations, as such you can never assume that you will reach absolute zero. This is the case with all floating point arithmetic, as discussed here and in many other places. The only reliable solution is to include a suitably chosen eps variable and a function to check.

So instead of checking result == 0 define isZero = lambda val: abs(val) < eps and check isZero(result).

This is a universal problem in floating point operations. In principle, using sympy, you can find real zeros because it is an algebra library, not a floating point math library. However, in the example you gave, not using N (which is what switches to float arithmetic), makes the computation extremely slow.

Community
  • 1
  • 1
mmdanziger
  • 4,466
  • 2
  • 31
  • 47
  • The function N do not return a regular floating point number. For exapmle, a=Symbols('a');N((1/a).subs(a, 3),100) will return a number with 100 threes . – iuradz Jul 17 '15 at 07:29
  • But still, it is a finite precision floating point number, and all limitations of regular floats apply. – Rob Jul 17 '15 at 09:06
0

I made a mistake when trying mpmath. I have tried mpmath again and it's a perfect solution for this problem.

iuradz
  • 1,128
  • 1
  • 12
  • 25
  • I would caution you not to rely on this. It may work on occasion but as long as you are in a float world, you do not reliably have equality between different operations, even if they are mathematically equivalent . – mmdanziger Jul 17 '15 at 09:08
  • Thank you. I'm aware of this. I think I descripted my question in a wrong way. I'm not expecting totally equality. But the precision of double float is not enough for me. I need a higher precision, like 10 times of double floating. This is exactly what mpmath designed for. – iuradz Jul 17 '15 at 14:30