加速numpy.dot函数

Question

我有一个使用numpy的脚本，其中约50％的运行时间花费在以下代码中：

s = numpy.dot(v1, v1)

其中

v1 = v[1:]

v是一个包含4000个元素的连续内存的1D ndarray，其数据类型为float64 （v.strides为(8,)）。

有什么建议可以加速这个过程吗？

编辑 这是在英特尔硬件上。这是我numpy.show_config()的输出：

atlas_threads_info:
    libraries = ['lapack', 'ptf77blas', 'ptcblas', 'atlas']
    library_dirs = ['/usr/local/atlas-3.9.16/lib']
    language = f77
    include_dirs = ['/usr/local/atlas-3.9.16/include']

blas_opt_info:
    libraries = ['ptf77blas', 'ptcblas', 'atlas']
    library_dirs = ['/usr/local/atlas-3.9.16/lib']
    define_macros = [('ATLAS_INFO', '"\\"3.9.16\\""')]
    language = c
    include_dirs = ['/usr/local/atlas-3.9.16/include']

atlas_blas_threads_info:
    libraries = ['ptf77blas', 'ptcblas', 'atlas']
    library_dirs = ['/usr/local/atlas-3.9.16/lib']
    language = c
    include_dirs = ['/usr/local/atlas-3.9.16/include']

lapack_opt_info:
    libraries = ['lapack', 'ptf77blas', 'ptcblas', 'atlas']
    library_dirs = ['/usr/local/atlas-3.9.16/lib']
    define_macros = [('ATLAS_INFO', '"\\"3.9.16\\""')]
    language = f77
    include_dirs = ['/usr/local/atlas-3.9.16/include']

lapack_mkl_info:
  NOT AVAILABLE

blas_mkl_info:
  NOT AVAILABLE

mkl_info:
  NOT AVAILABLE

- NPE

---

可以提供一下时间结果吗？顺便说一句，在我的简陋机器上，形状为（4000，）的随机向量使用点乘（dot(., .)）需要大约6微秒。谢谢。 - eat

@eat：在我的机器上，相同的操作只需要不到5微秒。虽然我正在做很多这样的操作，但它们会累加起来。 - NPE

好的，单个“点”似乎是相当高效的。然而，如果您想展示更多的代码，有人可能会找到如何优化计算的方法。谢谢。 - eat

@tillsten：是的，谢谢，它慢了好几倍（可能是因为它使用了一个临时数组）。 - NPE

@OferHelman: 这是相当久以前的事情了（3.5年），我很抱歉我不太记得了。:-（但我似乎记得我最终使用了Intel MKL。我不确定它是否对这个特定问题有所帮助。 - NPE

Answer 1

也许罪魁祸首是复制传递给 dot 的数组。正如Sven所说，点积依赖于BLAS操作。这些操作需要按照连续的C顺序存储的数组。如果传递给 dot 的两个数组都是以C_CONTIGUOUS方式存储的，则应该会看到更好的性能。

当然，如果您传递给dot的两个数组确实是1D（8,），则应将C_CONTIGUOUS和F_CONTIGUOUS标志都设置为True；但如果它们是（1,8），则可能会出现混合顺序。

>>> w = NP.random.randint(0, 10, 100).reshape(100, 1)
>>> w.flags
   C_CONTIGUOUS : True
   F_CONTIGUOUS : False
   OWNDATA : False
   WRITEABLE : True
   ALIGNED : True
   UPDATEIFCOPY : False

一种替代方法：使用BLAS中的_GEMM，该函数通过模块scipy.linalg.fblas公开。（显然，由于使用了fblas，数组A和B都是Fortran顺序。）

from scipy.linalg import fblas as FB
X = FB.dgemm(alpha=1., a=A, b=B, trans_b=True)

- doug

---

只是为了明确，输入数组是一维的（输出是标量）。 - NPE

Answer 2

你的数组不是很大，所以ATLAS可能没有发挥太大作用。以下是Fortran程序的计时结果？假设ATLAS没有发挥太大作用，这应该可以让你了解如果没有Python开销，dot()函数的速度有多快。使用gfortran -O3编译，我得到了5 +/- 0.5微秒的速度。

    program test

    real*8 :: x(4000), start, finish, s
    integer :: i, j
    integer,parameter :: jmax = 100000

    x(:) = 4.65
    s = 0.
    call cpu_time(start)
    do j=1,jmax
        s = s + dot_product(x, x)
    enddo
    call cpu_time(finish)
    print *, (finish-start)/jmax * 1.e6, s

    end program test

- matt

Answer 3

我能想到的唯一加速方法是确保你的NumPy安装已经编译了优化的BLAS库（比如ATLAS）。numpy.dot()是仅有的几个使用BLAS的NumPy函数之一。

- Sven Marnach

---

好建议（+1）。我已经更新了我的问题，并附上了我的numpy配置。它似乎是针对ATLAS构建的。 - NPE

@aix：你的配置看起来很好（虽然我不太确定如何解释它 :) ）。但仔细想想，当你的代码大部分时间都在乘以中等大小的向量时，你可能已经接近内存带宽的极限了，因此任何优化都只会使处理器更长时间地等待新数据。 - Sven Marnach

Answer 4

如果编译正确，numpy.dot方法将使用多线程处理。确保使用top检查。我知道有些人在使用atlas的numpy库时无法实现多线程处理。

此外，建议尝试使用针对intel mkl库进行编译的numpy版本。它们包含的blas例程应该比intel硬件上的atlas更快。

你可以尝试使用enthought的python发行版。它包含了所有这些内容，并且对于拥有edu邮件帐户的人来说是免费的。

- Lucas

---

你能提供一个链接，让我们获取一个针对Intel MKL库编译的NumPy版本吗？ - Ofer Helman

Anaconda通常会这样做。 - dermen