**TL;DR：**问题中提供的所有方法都非常低效。事实上，Numpy显然是次优的，没有办法只使用Numpy有效地计算它。不过，还有比问题中提供的更快的解决方案。

讲解及快速实现

Numpy代码使用强大的通用 * 迭代器 * 将给定的计算（如矩阵乘法）应用于多个数组切片。这样的迭代器对于实现 broadcasting 以及生成einsum的相对简单的实现非常有用。问题是，当迭代次数很大而数组很小时，它们的开销也很大。这正是在您的用例中发生的情况。虽然Numpy迭代器的开销可以通过优化Numpy代码来减轻，但在这个特定用例中，没有办法将开销减少到可以忽略不计的时间。实际上，每个矩阵只有12个浮点运算要执行。相对较新的x86-64处理器能够使用标量FMA单元在小于10纳秒内计算每个矩阵。事实上，可以使用SIMD指令，以便在仅几纳秒内计算每个矩阵。
首先，我们可以通过自己对第一个轴沿着的向量进行矩阵乘法来消除内部Numpy迭代器的开销。这正是explicit_2x2_matrices_multiplication所做的！
虽然explicit_2x2_matrices_multiplication应该快得多，但它仍然是次优的：它执行非连续内存读取，创建几个无用的临时数组，每个Numpy调用都会引入少量开销。更快的解决方案是编写Numba/Cython代码，以便底层编译器可以生成针对2x2矩阵优化的非常好的指令序列。好的编译器甚至可以在这种情况下生成SIMD指令，这对于Numpy来说是不可能的。下面是生成的代码：

import numba as nb
@nb.njit('(float64[:,:,::1], float64[:,:,::1])')
def compute_fastest(matrices_a, matrices_b):
    assert matrices_a.shape == matrices_b.shape
    assert matrices_a.shape[1] == 2 and matrices_a.shape[2] == 2

    n = matrices_a.shape[0]
    result_matrices = np.empty((n, 2, 2))
    for k in range(n):
        for i in range(2):
            for j in range(2):
                result_matrices[k,i,j] = matrices_a[k,i,0] * matrices_b[k,0,j] + matrices_a[k,i,1] * matrices_b[k,1,j]

    return result_matrices

字符串

性能结果

下面是我的机器上的性能结果，使用i5- 9600 KF CPU，用于1000 x2 x2矩阵：

Naive einsum:                           214   µs
matrices_a @ matrices_b:                102   µs
explicit_2x2_matrices_multiplication:    24   µs
compute_fastest:                          2.7 µs   <-----

型

讨论

Numba实现达到4.5 GFlops。每个矩阵的计算仅需12个CPU周期（2.7 ns）！我的机器在实践中能够达到300 GFlops（理论上为432 GFlops），但只有50 GFlops与1核心和12.5 GFlops使用标量代码（理论上为18 GFlops）。操作的粒度太小，多线程无法使用（生成线程的开销至少为几微秒）。此外，SIMD代码几乎不能使de FMA单元饱和，因为输入数据布局需要SIMD混洗，因此50 GFlops实际上是乐观的上限。因此，我们可以有把握地说Numba实现是一个非常有效的实现。尽管如此，由于SIMD指令，可以编写更快的代码（我希望在实践中加速约x2）。话虽如此，使用SIMD内部函数编写本机代码以帮助编译器生成快速SIMD代码确实远非易事（更不用说最终代码将是丑陋且难以维护的）。因此，SIMD实现可能不值得付出努力。

第一个“批次”维度上matmul的粗略线性验证：
你的（1000，2，2）数组：

In [353]: timeit matrices_a@matrices_b
251 µs ± 767 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

字符串
和一半和十分之一的大小：

In [354]: timeit matrices_a[:500]@matrices_b[:500]
129 µs ± 783 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
In [355]: timeit matrices_a[:100]@matrices_b[:100]
28.7 µs ± 532 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

型
你的露骨

In [360]: explicit_2x2_matrices_multiplication(matrices_a, matrices_b).shape
Out[360]: (1000, 2, 2)
In [361]: timeit explicit_2x2_matrices_multiplication(matrices_a, matrices_b)
59.9 µs ± 1.91 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

型
np.einsum不会尝试任何重新排序或其他优化：

In [362]: print(np.einsum_path('ijk,ikl->ijl',matrices_a, matrices_b, optimize='greedy
     ...: ')[1])
  Complete contraction:  ijk,ikl->ijl
         Naive scaling:  4
     Optimized scaling:  4
      Naive FLOP count:  1.600e+04
  Optimized FLOP count:  1.600e+04
   Theoretical speedup:  1.000
  Largest intermediate:  4.000e+03 elements
--------------------------------------------------------------------------
scaling                  current                                remaining
--------------------------------------------------------------------------
   4                ikl,ijk->ijl                                 ijl->ijl

型