Multidimensional matrix exponentiation
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.