# Why expm(2*A) != expm(A) @ expm(A)

1049 views
7

According to Matrix exponential, if `XY = YX`, then `exp(X)exp(Y) = exp(X+Y)`. However when I run the following in Python:

``````import numpy as np
from scipy.linalg import expm

A = np.arange(1,17).reshape(4,4)

print(expm(2*A))
[[ 306.63168024  344.81465009  380.01335176  432.47730444]
[ 172.59336774  195.36562731  214.19453937  243.76985501]
[ -35.40485583  -39.87705598  -42.94545895  -50.01324379]
[-168.44316833 -190.32607875 -209.76427134 -237.72069322]]

print(expm(A) @ expm(A))
[[1.87271814e+30 2.12068332e+30 2.36864850e+30 2.61661368e+30]
[4.32685652e+30 4.89977229e+30 5.47268806e+30 6.04560383e+30]
[6.78099490e+30 7.67886126e+30 8.57672762e+30 9.47459398e+30]
[9.23513328e+30 1.04579502e+31 1.16807672e+31 1.29035841e+31]]
``````

I get two very different results. Note that `@` is just the dot product.

I also tried it in Matlab and the two results are the same as expected. What am I missing here?

Maybe this is related to the rounding errors?

What version of `numpy`? on `'1.15.0'`, `(expm(2*A) == expm(A)@expm(A)).all()` returns `TRUE`

I have Numpy '1.15.3'

11

To me the root cause of this issue should be similar to the one described in the following NumPy bug report.

posted this