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
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.