在 Python 中实现大型 ndarray 乘法的最快方法

Question

我目前正在研究 Python 耦合 HMM 的实现，它涉及元素乘法、点积和维数 (50,50,40,40,40) 的 ndarrays 和维度 (40,40,40,40) 的 ndarrays。也就是说，非常大......第二个ndarray或多或少是稀疏的，大约有70％的零。

我一直在使用 numpy.einsum 直到现在，它给了我相当慢的结果（算法需要大约 3 小时到运行）。现在的问题是我需要为我的 HMM 优化一些参数，这意味着我将至少需要 10000 运行s，因此我需要将速度提高至少 1000 倍为了让事情保持合理。

我一直在四处寻找在 Python 中执行大型数组操作的最佳方法，现在我对所读的所有内容感到迷茫。所以我有几个问题。在询问他们之前，我只想说明我在 OSX、intel 处理器和 nvidia GPU 的机器上使用带有 Python 3 的最新 anaconda 发行版。另外，我相信我可以展平我的 ndarrays 以将我的问题简化为一个简单的矩阵-矩阵问题。

1)看来BLAS/MKL 图书馆可以提供相当不错的增长。当将 Ananaconda 与 OSX 一起使用时，MKL 似乎也原生链接到 Python。因此，直到现在我一直在使用 MKL，而不自知吗？ np.config.show() 给我这样的东西：

openblas_lapack_info:
  NOT AVAILABLE
blas_opt_info:
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['//anaconda/include']
    libraries = ['mkl_intel_lp64', 'mkl_intel_thread', 'mkl_core', 'iomp5', 'pthread']
    library_dirs = ['//anaconda/lib']
lapack_opt_info:
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['//anaconda/include']
    libraries = ['mkl_lapack95_lp64', 'mkl_intel_lp64', 'mkl_intel_thread', 'mkl_core', 'iomp5', 'pthread']
    library_dirs = ['//anaconda/lib']
lapack_mkl_info:
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['//anaconda/include']
    libraries = ['mkl_lapack95_lp64', 'mkl_intel_lp64', 'mkl_intel_thread', 'mkl_core', 'iomp5', 'pthread']
    library_dirs = ['//anaconda/lib']
blas_mkl_info:
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['//anaconda/include']
    libraries = ['mkl_intel_lp64', 'mkl_intel_thread', 'mkl_core', 'iomp5', 'pthread']
    library_dirs = ['//anaconda/lib']
mkl_info:
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['//anaconda/include']
    libraries = ['mkl_intel_lp64', 'mkl_intel_thread', 'mkl_core', 'iomp5', 'pthread']
    library_dirs = ['//anaconda/lib']

einsum 函数真的可以提高速度吗？如果没有，我如何设法从 Python?

访问 MKL 的 gemm 函数？

2) 这些 Python 库中的任何一个对我的问题有用吗：Theanos、Numba、Cupy？或者只是 Anaconda Accelerate？来自 Anaconda Accelerate 的 cuBlas 是否是我的最佳选择（即我愿意使用单精度或半浮点精度）？

3)用 C 或 C++ 等方式重新编码我的算法是否有用？

4)我尝试使用稀疏矩阵 (scipy.sparse.csc_matrix) 实现我的算法，但它使我的程序变慢了很多。这是意料之中还是我犯了错误？

我知道这会引起很多问题，但这是我第一次遇到这种问题，而且互联网上对此并不是很清楚...

非常感谢！

Answer 1

好吧，由于最近几天我花了很多时间来调查我的问题，我也许可以为像我一样迷茫的人提供一些答案。所以：

1) 是的，默认情况下 MKL 本机安装在 anaconda 中（尽管情况并非总是如此）。但是，np.einsum 并没有从中获益，因为它没有使用 BLAS 优化的 dot numpy 函数。因此，尝试从 mkl 手动访问 gemm 函数是没有用的（即使由于库 anaconda accelerate 这实际上很容易）。使用 np.tensordot、np.multiply 和其他方法重写 np.einsum 的操作要容易得多。

2) 我没有时间探索我询问的所有库，但是张量点 Theanos 操作显然并不比简单的 np.tensordot 快。同样，几年前情况并非如此，因为 np.tensordot 操作并未将 BLAS 用于多维数组。关于所有 cuda 库，与 MKL 相比，它们似乎实际上并没有那么好，除非你有一个非常强大的 CGU（例如参见 [=10=]）

3) 将一些使用 MKL 优化的 Python 代码重写到 C++ 中并没有带来太大的提升（参见 Benchmarking (python vs. c++ using BLAS) and (numpy) ）

4) 仍然不知道为什么我的稀疏矩阵实现比密集矩阵慢。很有可能我做错了。

在 Python 中实现大型 ndarray 乘法的最快方法

Fastest way to implement large ndarray multiplication in Python

python

matrix-multiplication

intel-mkl

anaconda