I am trying to compute a matrix exponential using scipy.sparse.linalg.expm but i realize it does not provide correct result when the numbers in the matrix are getting moderately high.
As a sanity check, i created a simple matrix
X = np.array([[06.1, 52.0],[016., 01.6 ]])
and then compute the product of exp(X)*exp(-X) to recover the identity matrix
Id = np.dot(expm(X),expm(-X))
I am expecting to get Id = ([[1, 0],[0, 1]]) but i get:
[[ 2.39599500e+08 -2.12377558e+08],[-6.89818242e+07 4.93089944e+08]]
Here is the full code:
import numpy as np
from scipy.sparse.linalg import expm
X = np.array([[06.1, 52.0],
[016., 01.6 ]])
Id = np.dot(expm(X),expm(-X))
print(Id)
To be noted that if I used smaller numbers in the matrix i will eventually get
Id = ([[1, 0],[0, 1]])
Another puzzling fact, a colleague of mine managed to get a consistent result when typing the exact same code as above on his python software.
Any thought? Thanks in advance
I think you triggered a known issue of scipy expm for 2x2 matrices, see https://github.com/scipy/scipy/issues/19584 and https://github.com/scipy/scipy/pull/20261. The fix is recent and probably not yet in the scipy version you are using. Also see the work in progess https://github.com/scipy/scipy/issues/20262
If you really need to compute expm of this matrix, one solution may be do normalize it by some integer and call np.linalg.matrix_power. Not sure this really improves numerical precision.