首先要做的就是.一个人需要understand是怎么优化这样一个代码的，然后是profile，然后是estimate the time，只有到那时才能找到更好的解决方案.
TL;DR:这两个版本效率都很低，并且是串联运行的.Neither BLAS libraries nor Numpy are designed for optimizing this use-case.事实上，当最后一个轴非常小时，即使是基本的块操作也不是有效的(即.4).这可以通过使用Numba编写一个专门为您的矩阵大小而设计的实现来优化.
理解
@是一个Python操作符，但它在内部调用Numpy函数，例如+或*.它对所有矩阵执行循环迭代，并在每个矩阵上调用高度优化的BLAS实现.BLAS是一个数值代数库.有很多现有的BLAS，但Numpy的默认BLAS通常是OpenBLAS，这是非常优化的，特别是对于大型矩阵.还请注意，np.einsum可以调用特定模式的BLAS实现(但如果正确设置了optimize标志)，但这里的情况并非如此.还值得一提的是，np.einsum对于输入中的2个矩阵进行了很好的优化，但对于3个矩阵则不是很好，并且没有针对参数中的更多矩阵进行优化.这是因为可能性的数量呈指数级增长，并且代码手动进行优化.有关np.einsum如何工作的更多信息，请阅读How is numpy.einsum implemented?.
配置文件
问题是你要乘以many very-small matrices and most BLAS implementations are not optimized for that.事实上，Numpy也是如此:与计算相比，通用循环迭代的成本可能会变得很大，更不用说对BLAS的函数调用了.对Numpy的分析表明，np.einsum实现中最慢的函数是PyArray_TransferNDimToStrided.此函数不是主计算函数，而是辅助函数.事实上，主计算功能只需要总时间的20%，这留下了很大的改进空间！对于BLAS实现也是如此:cblas_dgemv仅占用约20%以及dgemv_n_HASWELL(BLAS cblas_dgemv函数的主要计算内核).其余的几乎都是BLAS库或Numpy的开销(两者的时间大约都有一半).此外，both version runs serially.实际上，np.einsum并没有优化为与多线程一起运行，BLAS不能使用多线程，因为矩阵太小，所以多线程可能很有用(因为多线程有很大的开销).这意味着both versions are pretty inefficient.
性能指标
要知道这些版本的效率有多低，人们需要知道要做的计算量和处理器的速度.触发器(浮点运算)的个数由np.einsum_path提供，为5.120e+05(对于优化实现，否则为6.144e+05).主流CPU通常以多线程和数十个GFlop/s的速度并行执行>；=np.einsum_pathGFlop/s.例如，我的i5-9600KF处理器可以实现300-400GFlops/s的并行和50-60GFlopps/s的串行化.由于BLAS版本(BEST)的计算持续0.52ms，这意味着代码以1GFlops/s的速度运行，与最佳结果相比，这是一个很差的结果.
最佳化
加速计算的解决方案是设计一个针对您特定大小的矩阵进行优化的Numba(JIT compiler)或Cython(从Python到C编译器)实现.事实上，对于通用代码来说，最后一个维度太小了，速度不快.在这种情况下，即使是基本的编译代码也不会很快:与实际计算相比，即使是C循环的开销也可能相当大.我们可以告诉编译器某个矩阵轴的大小在编译时很小并且是固定的，因此编译器可以生成更快的代码(这要归功于循环展开、平铺和SIMD指令).这可以通过Numba中的基本断言来完成.此外，如果没有使用特殊的浮点(FP)值，如NaN或不正常的数字，我们可以使用fastmath=True标志来加快计算速度.此标志还会影响结果的准确性，因为它假设FP数学是关联的(这不是真的).简而言之，为了提高性能，它打破了IEEE-754标准.以下是生成的代码:
import numba as nb
# use `fastmath=True` for better performance if there is no
# special value used and the accuracy is not critical.
@nb.njit('(float64[:,:,::1], float64[:,:,::1])', fastmath=True)
def compute_fast_einsum(diff, inv_covs):
ni, nj, nk = diff.shape
nl = inv_covs.shape[2]
assert inv_covs.shape == (nj, nk, nl)
assert nk == 4 and nl == 4
res = np.empty((ni, nj), dtype=np.float64)
for i in range(ni):
for j in range(nj):
s = 0.0
for k in range(nk):
for l in range(nl):
s += diff[i, j, k] * inv_covs[j, k, l] * diff[i, j, l]
res[i, j] = s
return res
%%timeit
diff = xvals - means.squeeze(1)
compute_fast_einsum(diff, inv_covs)
结果
以下是我的计算机上的性能结果(平均值±标准戴夫.共7次运行，每次1000次循环):
@ operator: 602 µs ± 3.33 µs per loop
einsum: 698 µs ± 4.62 µs per loop
Numba code: 193 µs ± 544 ns per loop
Numba + fastmath: 177 µs ± 624 ns per loop
Best Numba: < 100 µs <------ 6x-7x faster !
请注意，diff的计算花费了diff微秒，这是低效的.这一点也可以使用Numba进行优化.实际上，在基于i循环中，可以从其他数组即时计算diff的值.这使得计算对缓存更加友好.这个版本在搜索结果中被称为"BEST Numba".请注意，Numba版本甚至没有使用多线程.尽管如此，多线程的开销通常约为5-500微秒，因此在一些机器上使用多线程可能会更慢(在主流PC上，即.不是计算服务器，开销通常是5-diff微秒，在我的机器上大约是10微秒).