Compute the matrix exponential of an array.
Input with last two dimensions are square
(..., n, n).
The resulting matrix exponential with the same shape of
Implements the algorithm given in , which is essentially a Pade approximation with a variable order that is decided based on the array data.
For input with size
n, the memory usage is in the worst case in the order of
8*(n**2). If the input data is not of single and double precision of real and complex dtypes, it is copied to a new array.
n >= 400, the exact 1-norm computation cost, breaks even with 1-norm estimation and from that point on the estimation scheme given in  is used to decide on the approximation order.
Awad H. Al-Mohy and Nicholas J. Higham, (2009), “A New Scaling and Squaring Algorithm for the Matrix Exponential”, SIAM J. Matrix Anal. Appl. 31(3):970-989, DOI:10.1137/09074721X
Nicholas J. Higham and Francoise Tisseur (2000), “A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra.” SIAM J. Matrix Anal. Appl. 21(4):1185-1201, DOI:10.1137/S0895479899356080
>>> from scipy.linalg import expm, sinm, cosm
Matrix version of the formula exp(0) = 1:
>>> expm(np.zeros((3, 2, 2))) array([[[1., 0.], [0., 1.]], [[1., 0.], [0., 1.]], [[1., 0.], [0., 1.]]])
Euler’s identity (exp(i*theta) = cos(theta) + i*sin(theta)) applied to a matrix:
>>> a = np.array([[1.0, 2.0], [-1.0, 3.0]]) >>> expm(1j*a) array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j], [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]]) >>> cosm(a) + 1j*sinm(a) array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j], [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])