迈向外存 ND 数组 -- 矩阵乘法基准测试

这项工作得到了 Continuum Analytics 和 XDATA 项目（作为 Blaze 项目的一部分）的支持

太长不看 我们在外存点积上对 dask 进行了基准测试。我们还比较并说明了使用优化 BLAS 的动机。

结果如下

免责声明：此文基于实验性的、可能包含 bug 的代码。它尚未准备好供公众使用。

引言

这是使用 NumPy、Blaze 和 dask 构建外存 nd 数组系列文章的第四篇。您可以在这里查看这些文章

1. 简单任务调度,
2. 前端可用性
3. 多线程调度器

我们现在给出外存矩阵乘法的性能数据。

矩阵乘法

稠密矩阵乘法是计算密集型的，而非 I/O 密集型的。我们大部分时间都在进行算术运算，相对较少的时间用于数据传输。因此，我们可能能够在不损失性能的情况下从磁盘读取大型数据。

当两个 $n\times n$ 矩阵相乘时，我们读取 $n^2$ 字节，但执行 $n^3$ 次计算。每字节有 $n$ 次计算，所以相对而言，I/O 是廉价的。

我们通常使用每秒千兆浮点运算次数（GFLOPS）来衡量单 CPU 的速度。我们来看看我的笔记本电脑使用 NumPy 在单线程内存中执行矩阵乘法的表现。

>>> import numpy as np
>>> x = np.ones((1000, 1000), dtype='f8')
>>> %timeit x.dot(x)  # Matrix-matrix multiplication
10 loops, best of 3: 176 ms per loop

>>> 1000 ** 3 / 0.176 / 1e9  # n cubed computations / seconds / billion
>>> 5.681818181818182

好的，NumPy 的矩阵乘法速度大约在 6 GFLOPS。在我这台机器上，np.dot 函数调用的是 BLAS 库中的 GEMM 操作。目前我的 numpy 只使用了参考 BLAS。（您可以使用 np.show_config() 来检查。）

从磁盘进行矩阵乘法

对于大到无法完全载入内存的矩阵，我们分部分计算结果，必要时从磁盘载入块。我们使用多线程来并行化此过程。我们上一篇文章展示了 NumPy+Blaze+Dask 如何为我们自动完成这一切。

我们进行一个简单的数值实验，使用 HDF5 作为磁盘存储。

我们安装所需库

conda install -c blaze blaze
pip install git+https://github.com/ContinuumIO/dask

我们在磁盘上设置一个模拟数据集

>>> import h5py
>>> f = h5py.File('myfile.hdf5')
>>> A = f.create_dataset(name='A', shape=(200000, 4000), dtype='f8',
...                                chunks=(250, 250), fillvalue=1.0)
>>> B = f.create_dataset(name='B', shape=(4000, 4000), dtype='f8',
...                                chunks=(250, 250), fillvalue=1.0)

我们告诉 Dask+Blaze 如何与该数据集交互

>>> from dask.obj import Array
>>> from blaze import Data, into

>>> a = into(Array, 'myfile.hdf5::/A', blockshape=(1000, 1000))  # dask things
>>> b = into(Array, 'myfile.hdf5::/B', blockshape=(1000, 1000))
>>> A = Data(a)  # Blaze things
>>> B = Data(b)

我们计算所需结果，并将其存回磁盘

>>> %time into('myfile.hdf5::/result', A.dot(B))
2min 49s

>>> 200000 * 4000 * 4000 / 169 / 1e9
18.934911242

18.9 GFLOPS，大约比内存中计算快 3 倍。乍一看这令人困惑——从磁盘读取不是应该更慢吗？我们的加速得益于并行使用了四个核心。这很好，我们没有经历从磁盘读取带来的太多减速。

就好像我们的整个硬盘都变成了内存一样。

OpenBLAS

参考 BLAS 很慢；它是很久以前写的。OpenBLAS 是一个现代实现。我使用系统安装程序（apt-get）安装了 OpenBLAS，然后重新配置并重建了 numpy。OpenBLAS 支持多核心。我们将展示使用一个线程和四个线程的计时。

export OPENBLAS_NUM_THREADS=1
ipython

>>> import numpy as np
>>> x = np.ones((1000, 1000), dtype='f8')
>>> %timeit x.dot(x)
10 loops, best of 3: 89.8 ms per loop

>>> 1000 ** 3 / 0.0898 / 1e9  # compute GFLOPS
11.135857461024498

export OPENBLAS_NUM_THREADS=4
ipython

>>> %timeit x.dot(x)
10 loops, best of 3: 46.3 ms per loop

>>> 1000 ** 3 / 0.0463 / 1e9  # compute GFLOPS
21.598272138228943

这比参考实现快了大约四倍。如果您还没有通过其他方式（比如 dask）进行并行化，那么您应该使用像 OpenBLAS 或 MKL 这样的现代 BLAS。

OpenBLAS + dask

最后，我们再次运行我们的磁盘实验，这次使用了 OpenBLAS。我们分别测试了 OpenBLAS 在单线程和多线程下的表现。

我们将跳过代码（它与上面相同），并在下面提供一个综合的结果表格。

遗憾的是，外存解决方案在使用 OpenBLAS 后并没有显著提升。实际上，当 OpenBLAS 和 dask 都尝试并行化时，性能反而会下降。

结果

太长不看 在进行计算密集型工作时，不必担心使用磁盘，但不要同时使用两种并行机制。

主要收获

1. 在计算密集型任务中，从磁盘操作不会损失太多性能
2. 实际上，当没有优化 BLAS 可用时，我们可以提升性能。
3. Dask 并未从优化的 BLAS 中获得多少益处。这令人沮丧且惊讶。我曾期望性能能够随着单核内存性能而扩展。也许这表明存在其他限制因素
4. 不应过分地用这些数字进行推断。它们仅与矩阵乘法等高度计算密集型操作相关。

另外，感谢 Wesley Emeneker 找到了我们的内存泄漏之处，使得这些结果成为可能。