为了更明确，这里是如何（如果我正确理解你的最后一个轴是什么，那只是执行N次线性代数运算的方法;这是numpy通常处理的方式，相反，第一个轴），您可以使用np.moveaxis。
例如，指数的预期结果（同样，如果我理解正确的话），可以通过这种方式获得

np.moveaxis(scipy.linalg.expm(np.moveaxis(M,2,0)),0,2)

也就是说，将最后一个轴移动到第一个位置：np.moveaxis(M,2,0)是一个N×3×3阵列。expm执行N次指数运算，对N个3x 3矩阵中的每个矩阵执行一次，并返回N个3x 3数组，对应于N个结果。由于您想要的是一个3x 3xN的数组，所以您可以将结果的第一个轴移回最后一个位置。
请注意，np.moveaxis几乎不做任何事情。它只是改变了阵列的strides，这会影响数据的访问方式。它不分配任何新的内存（它是相同的数据），它不移动任何数据。np.moveaxis(M,2,0)[0,1,2]是与M[1,2,0]相同的元素，但有另一个名称（我的意思是不相同的“值”，内存中的相同数据; np.moveaxis(M,2,0)[0,1,2]=15实际上改变了M[1,2,0]的值。因此，简而言之，它不需要花费任何费用，或者几乎不需要花费任何费用。
但这并不意味着你可以移动轴，然后对N×3×3矩阵和N×3向量执行任何线性代数，就像你对3x 3矩阵和3-向量一样。大多数numpy操作都是这样工作的，但不是全部。有时候这是不可能的。
以@矩阵向量乘法为例。它不能工作，因为它是模棱两可的。A@B是矩阵-矩阵乘法，如果A和B是二维阵列。如果A是一个2D数组，B是一个1D数组，则它是一个矩阵向量乘法。
但是如果A是一个3D数组，V是一个2D数组呢？在numpy的精神中，我们使用额外的轴进行“矢量化”。但这必须被理解为N个2D-2D阵列（总是具有相同的V）吗？还是N二维-一维乘法？几乎是一样的，你会说。毕竟，根据定义，矩阵-矩阵乘法只是M矩阵-向量乘法。但在这种情况下，“额外”维度是最后一个，而不是第一个。
因此，简而言之，它是模糊的。并且numpy中每个操作数超过2个轴的@选择的定义不是您想要的。
话虽如此，您可以使用其他方法。例如，你可以在你的vector中添加另一个1大小的维度，这样它就被@视为2D（只有1列，在数学上，这与vector相同，但在numpy中并不完全相同，例如，在这种情况下就不是这样），

# move axis of M to have a N×3×3 array (aka N 3x3 matrix)   moveaxis
# move axis of V to have a N×3 array (aka N vectors)        moveaxis
# add a dimension V to have N×3×1 array (aka N 3x1 matrix)  [:,:,None]
# perform N matrix-matrix multiplication  (N (3×3)@(3×1))   @
# result is N×3×1 array
# remove the last axis to get a N×3 array                   [:,:,0]
# move back axis to get a 3×N array
MV = np.moveaxis((np.moveaxis(M,2,0) @ np.moveaxis(V,1,0)[:,:,None])[:,:,0], 0, 1)

您还可以使用einsum（至少更短）以更简单的方式完成此操作。但要理解完整的einsum需要一些挠头）

MV = np.einsum('ij...,j...->i...', M,V)

另一方面，矩阵-矩阵乘法的工作原理类似于expm示例。也就是说，如果你有两个N×3×3矩阵，那么只使用@就可以简单地在3x 3矩阵上执行N次乘法。由于您拥有的是3x3×N矩阵，因此在简单的情况下，只需移动轴返回，然后将结果的轴移动回来，就像np.expm一样

np.moveaxis(np.moveaxis(A,2,0) @ np.moveaxis(B,2,0), 0,2)

或者，同样是einsum（但在这种情况下，这是不必要的，因为一般的moveaxis(op(moveaxis(...)))技巧是有效的。

np.einsum('ij...,jk...->ik...',A,B)

你需要做一些修改。
首先，你的数组应该是统一的：

import numpy as np

A = np.array([[[1,2], [0,0], [0,0]],
              [[0,0], [2,4], [0,0]],
              [[0,0], [0,0], [3,6]]])

>>> A
array([[[1, 2],
        [0, 0],
        [0, 0]],

       [[0, 0],
        [2, 4],
        [0, 0]],

       [[0, 0],
        [0, 0],
        [3, 6]]])

第二，您要查找的矩阵运算是scipy.linalg.expm：

import scipy

B = scipy.linalg.expm(np.array([[1, 0, 0], [0, 2, 0], [0, 0, 3]]))

>>> B
array([[ 2.71828183,  0.        ,  0.        ],
       [ 0.        ,  7.3890561 ,  0.        ],
       [ 0.        ,  0.        , 20.08553692]])

最后，expm在最后两个维度上需要一个正方形数组，所以你需要稍微操作一下你的数组：

A1 = np.transpose(A, (2, 0, 1))
A2 = scipy.linalg.expm(A1)
A3 = np.transpose(A, (1, 2, 0))

>>> A1
array([[[1, 0, 0],
        [0, 2, 0],
        [0, 0, 3]],

       [[2, 0, 0],
        [0, 4, 0],
        [0, 0, 6]]])

>>> A2
array([[[  2.71828183,   0.        ,   0.        ],
        [  0.        ,   7.3890561 ,   0.        ],
        [  0.        ,   0.        ,  20.08553692]],

       [[  7.3890561 ,   0.        ,   0.        ],
        [  0.        ,  54.59815003,   0.        ],
        [  0.        ,   0.        , 403.42879349]]])

>>> A3
array([[[  2.71828183,   7.3890561 ],
        [  0.        ,   0.        ],
        [  0.        ,   0.        ]],

       [[  0.        ,   0.        ],
        [  7.3890561 ,  54.59815003],
        [  0.        ,   0.        ]],

       [[  0.        ,   0.        ],
        [  0.        ,   0.        ],
        [ 20.08553692, 403.42879349]]])

创建数组：

In [101]: A = np.empty((3,3),object); A[:]=0    
In [102]: A[[0,1,2],[0,1,2]] = [np.arange(3) for _ in range(3)]    
In [103]: A
Out[103]: 
array([[array([0, 1, 2]), 0, 0],
       [0, array([0, 1, 2]), 0],
       [0, 0, array([0, 1, 2])]], dtype=object)

不能全部使用：

In [104]: np.exp(A)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
AttributeError: 'numpy.ndarray' object has no attribute 'exp'

The above exception was the direct cause of the following exception:

TypeError                                 Traceback (most recent call last)
Cell In[104], line 1
----> 1 np.exp(A)

TypeError: loop of ufunc does not support argument 0 of type numpy.ndarray which has no callable exp method

np.vectorize或np.frompyfunc可以单独应用函数。但请注意非对角线的变化：

In [105]: np.frompyfunc(np.exp, 1,1)(A)
Out[105]: 
array([[array([1.        , 2.71828183, 7.3890561 ]), 1.0, 1.0],
       [1.0, array([1.        , 2.71828183, 7.3890561 ]), 1.0],
       [1.0, 1.0, array([1.        , 2.71828183, 7.3890561 ])]],
      dtype=object)

In [106]: np.frompyfunc(lambda x:2*x+1, 1,1)(A)
Out[106]: 
array([[array([1, 3, 5]), 1, 1],
       [1, array([1, 3, 5]), 1],
       [1, 1, array([1, 3, 5])]], dtype=object)

但是对于这样的操作员来说不需要这个

In [107]: 2*A+1
Out[107]: 
array([[array([1, 3, 5]), 1, 1],
       [1, array([1, 3, 5]), 1],
       [1, 1, array([1, 3, 5])]], dtype=object)

忽略offdiagonals将需要一些进一步的覆盖实用程序，例如数组测试：

In [110]: np.frompyfunc(lambda x:np.exp(x) if isinstance(x,np.ndarray) else x, 1,1)(A)
Out[110]: 
array([[array([1.        , 2.71828183, 7.3890561 ]), 0, 0],
       [0, array([1.        , 2.71828183, 7.3890561 ]), 0],
       [0, 0, array([1.        , 2.71828183, 7.3890561 ])]], dtype=object)

In [112]: np.frompyfunc(lambda x:2*x+1.1 if isinstance(x,np.ndarray) else x, 1,1)(A)
Out[112]: 
array([[array([1.1, 3.1, 5.1]), 0, 0],
       [0, array([1.1, 3.1, 5.1]), 0],
       [0, 0, array([1.1, 3.1, 5.1])]], dtype=object)

这种事情可以概括。只是不要期望像编译的数字numpy代码那样的性能。