问Numpy中的矩阵求逆速度

EN

Numpy使用LAPACK的LU分解(正如@Murali在注释中指出的那样)。LU分解通常是由OpenBLAS在大多数平台上完成的(或者是运行类似的Intel )。OpenBLAS使用dgemm调用(矩阵乘法)，从而使LU分解更快。dgemm是大矩阵最优化的BLAS原语之一。它利用多个线程和SIMD指令，所以非常快。尽管如此，在小型矩阵上，使用多线程会引入一个重要的开销(创建线程、分配工作和等待结果需要时间)。对于较大的矩阵，这种开销往往是可以忽略不计的。对于具有更小矩阵的SIMD指令也是如此: SIMD指令对于计算相对较大的数组来说是快速的，但与标量指令相比，它们具有较高的延迟时间(对于小型矩阵可以更好地流水线化)。循环展开也会产生同样的效果(对于非常小的矩阵)以及CPU缓存。在此基础上，Numpy完成了一些类型检查和内部工作(在其他情况下优化内部循环)，这只需几微秒。所有这些因素的结合导致代码对于小矩阵的效率很低，对于较大的矩阵则更快。这种开销在测量的速度上产生了偏差.

在我的机器上，对于时间为100的矩阵，dgemm调用占19%，dtrsm调用占11%，dlaswp调用占7% (其他调用占用可忽略的时间，其余调用在等待或内核中丢失)。对于大小为200的矩阵，dgemm调用占29%，dtrsm占12%，dlaswp占7%。对于大小为1000的矩阵，dgemm调用占57%，dtrsm调用占5%，dlaswp占6%。正如我们所看到的，大部分时间都浪费在大小为100的矩阵上，而大部分时间则被BLAS操作用于1000大小的矩阵。

注意，LU分解很难并行化。这是一个活跃的研究领域(持续数十年)。此外，不幸的是，最好的并行化方法(通常是基于任务的方法)还没有在像OpenBLAS这样的库中使用。AFAIK，MKL使用他们一点。最先进的实现速度更快，特别是对于小矩阵，因为dgemm对大型矩阵进行了适当的优化。

可以在OpenBLAS中使用OMP_NUM_THREADS=1禁用多线程。这样，我得到了151个大小为100的矩阵和800个大小为200的矩阵。我们应该考虑Numpy开销的~4和其他非常量的开销，如页面错误、分配、SIMD指令的延迟、缓存效果、分支预测等，也就是说，因子达到5.44 (没有Numpy开销)。对于200比400，是5.67。对于400比800，它甚至是6.8。除其他开销外，该因子仍不接近8的主要原因是，OpenBLAS没有被优化以运行具有1个线程的LU因式分解:在这种情况下，内部参数的设置是次优化的，而算法不是最优的。