从SciPy 1.9.0版（发布于2012年7月28日）开始，您可以传递scipy.linalg.expm一个最后两维为正方形的数组--即形状为(..., n, n)的数组--它将以矢量化的方式计算矩阵指数。
文档为here。

浏览代码https://github.com/scipy/scipy/blob/v1.4.1/scipy/sparse/linalg/matfuncs.py#L550-L595，看起来它实际上只假设了一个矩阵。因此，如果你想避免python循环，你需要提取相关的部分，然后自己对它们进行cythonize或numpy-vectorize。只有满足密集矩阵的要求才能让它更容易。

下面是pade近似的代码。这只是针对1个数组，你必须对它进行numpy矢量化/ cythonise。或者只使用numba中的@njit

def matrixexponen(n,a):
    q = 6
    a2 = a.copy ( )
    a_norm = np.linalg.norm ( a2, np.inf )
    ee = ( int ) ( np.log2 ( a_norm ) ) + 1
    s = max ( 0, ee + 1 )
    a2 = a2 / ( 2.0 ** s )
    x = a2.copy ( )
    c = 0.5
    e = np.eye ( n, dtype = np.complex64 ) + c * a2
    d = np.eye ( n, dtype = np.complex64 ) - c * a2
    p = True

    for k in range ( 2, q + 1 ):
        c = c * float ( q - k + 1 ) / float ( k * ( 2 * q - k + 1 ) )
        x = np.dot ( a2, x )
        e = e + c * x
        
        if ( p ):
            d = d + c * x
        else:
            d = d - c * x

        p = not p
    #  E -> inverse(D) * E
    e = np.linalg.solve ( d, e )
    #  E -> E^(2*S)
    for k in range ( 0, s ):
        e = np.dot ( e, e )
    return e