Is there a more efficient way of exponentiating a multidimensional matrix? This is how I am currently doing it:
import numpy as np from scipy import linalg A = np.random.normal(0,0.01, (3, 2,2)) P = np.zeros((3,2,2)) for i in range(3): P[i] = linalg.expm( np.multiply(A[i], 0.5)) Note that my real array is much larger and any effort to speed up the process is greatly appreciated.
Dig through the source code for expm, which leads ultimately to a rather extensive bit of code in scipy.sparse.linalg.matfuncs.py. This uses a well known Pade approximation. It is specifically designed to work with on a square nxn matrix. So generalizing to a (m, n, n) matrix without iterating on m does not work.
Years ago I sped this up by using f2py to calculate an analytical (n,3,3) inverse. But I recall a simpler expm implementation. The current one appears to have a lot of helper functions and special cases, presumably hoping to gain speed.
If your matrix is always (m, 2,2) there might be a special case.
This might be better as a Stackoverflow question, since it has less to do with programming style or code organization, and more with calculation alternatives in the numpy and scipy libraries. But a SO search for expm and numpy/scipy turns up very few questions. The closest suggests looking at the source code just as I suggested.
