加速python中的元素数组乘法

Posted

技术标签:

【中文标题】加速python中的元素数组乘法【英文标题】:Speeding up element-wise array multiplication in python 【发布时间】:2013-10-16 07:56:31 【问题描述】:

我一直在尝试使用 numba 和 numexpr 来加快简单的逐元素矩阵乘法。我一直没能得到更好的结果,它们基本上(速度方面)都相当于 numpys 乘法函数。有没有人在这方面有运气?我是否使用了 numba 和 numexpr 错误(我对此很陌生),或者这完全是一种尝试加快速度的坏方法。这是一个可重现的代码,在此先感谢您:

import numpy as np
from numba import autojit
import numexpr as ne

a=np.random.rand(10,5000000)

# numpy
multiplication1 = np.multiply(a,a)

# numba
def multiplix(X,Y):
    M = X.shape[0]
    N = X.shape[1]
    D = np.empty((M, N), dtype=np.float)
    for i in range(M):
        for j in range(N):
            D[i,j] = X[i, j] * Y[i, j]
    return D

mul = autojit(multiplix)
multiplication2 = mul(a,a)

# numexpr
def numexprmult(X,Y):
    M = X.shape[0]
    N = X.shape[1]
    return ne.evaluate("X * Y")

multiplication3 = numexprmult(a,a) 

【问题讨论】:

numexpr 可以胜过 numpy 这样的类似 ufunc 的操作,尤其是将几个串在一起。此外,如果您有多个内核,请尝试设置ne.set_num_cores(N),其中N 是您机器的内核数。 在我的机器上,基于 numexpr 的函数比在单核上运行的 np.multiply() 慢约 15%,但是当我将内核数设置为8. 请记住,您可能会发现您必须重置 Python 进程的核心亲和性才能使用多个核心 - see my answer here。 您可以尝试使用 Theano 来使用您的 GPU。我真的不知道它是否会有所帮助,结果将取决于您的确切硬件,但它可能值得一试。 Here 您将找到一个如何使用 Theano 进行元素矩阵乘法的示例。 如果可以,请将您的 numpy 更新到 1.8。 (在编写它时,即将发布),这应该会提供一个简单的加速。否则,您将不得不使用其他可以使用 SIMD 指令或可以优化您的处理器的东西。 【参考方案1】:

使用fortran 和ctypes 怎么样?

elementwise.F90:

subroutine elementwise( a, b, c, M, N ) bind(c, name='elementwise')
  use iso_c_binding, only: c_float, c_int

  integer(c_int),intent(in) :: M, N
  real(c_float), intent(in) :: a(M, N), b(M, N)
  real(c_float), intent(out):: c(M, N)

  integer :: i,j

  forall (i=1:M,j=1:N)
    c(i,j) = a(i,j) * b(i,j)
  end forall

end subroutine 

elementwise.py:

from ctypes import CDLL, POINTER, c_int, c_float
import numpy as np
import time

fortran = CDLL('./elementwise.so')
fortran.elementwise.argtypes = [ POINTER(c_float), 
                                 POINTER(c_float), 
                                 POINTER(c_float),
                                 POINTER(c_int),
                                 POINTER(c_int) ]

# Setup    
M=10
N=5000000

a = np.empty((M,N), dtype=c_float)
b = np.empty((M,N), dtype=c_float)
c = np.empty((M,N), dtype=c_float)

a[:] = np.random.rand(M,N)
b[:] = np.random.rand(M,N)


# Fortran call
start = time.time()
fortran.elementwise( a.ctypes.data_as(POINTER(c_float)), 
                     b.ctypes.data_as(POINTER(c_float)), 
                     c.ctypes.data_as(POINTER(c_float)), 
                     c_int(M), c_int(N) )
stop = time.time()
print 'Fortran took ',stop - start,'seconds'

# Numpy
start = time.time()
c = np.multiply(a,b)
stop = time.time()
print 'Numpy took ',stop - start,'seconds'

我编译了 Fortran 文件使用

gfortran -O3 -funroll-loops -ffast-math -floop-strip-mine -shared -fPIC \
         -o elementwise.so elementwise.F90

输出产生约 10% 的加速:

 $ python elementwise.py 
Fortran took  0.213667869568 seconds
Numpy took  0.230120897293 seconds
 $ python elementwise.py 
Fortran took  0.209784984589 seconds
Numpy took  0.231616973877 seconds
 $ python elementwise.py 
Fortran took  0.214708089828 seconds
Numpy took  0.25369310379 seconds

【讨论】:

可爱的答案。加速并不令人印象深刻,但我有兴趣玩这个,谢谢。 如杰奇华所说的可爱答案。但是,要获得准确的答案,必须执行第一次 fortran 调用来初始化共享库。第二个电话将给出最准确的答案。加速应该在 50% 左右。另一种最准确的方法是使用循环(假设 100 次调用同一函数)并取平均时间。 为什么加速会在 50% 左右?如何? @innoSPG @JEquihua,我忘了说 50% 是根据我自己的本地测试。谢谢你指出。这可能取决于您的系统配置。【参考方案2】:

你的时间安排如何?

您的随机数组的创建占用了您计算的全部部分,如果您将其包含在您的计时中,您几乎不会看到结果有任何真正的差异, 但是,如果您预先创建它,您实际上可以比较这些方法。

这是我的结果,我始终如一地看到您所看到的。 numpy 和 numba 给出的结果大致相同(numba 更快一点。)

(我没有可用的 numexpr)

In [1]: import numpy as np
In [2]: from numba import autojit
In [3]: a=np.random.rand(10,5000000)

In [4]: %timeit multiplication1 = np.multiply(a,a)
10 loops, best of 3: 90 ms per loop

In [5]: # numba

In [6]: def multiplix(X,Y):
   ...:         M = X.shape[0]
   ...:         N = X.shape[1]
   ...:         D = np.empty((M, N), dtype=np.float)
   ...:         for i in range(M):
   ...:                 for j in range(N):
   ...:                         D[i,j] = X[i, j] * Y[i, j]
   ...:         return D
   ...:         

In [7]: mul = autojit(multiplix)

In [26]: %timeit multiplication1 = np.multiply(a,a)
10 loops, best of 3: 182 ms per loop

In [27]: %timeit multiplication1 = np.multiply(a,a)
10 loops, best of 3: 185 ms per loop

In [28]: %timeit multiplication1 = np.multiply(a,a)
10 loops, best of 3: 181 ms per loop

In [29]: %timeit multiplication2 = mul(a,a)
10 loops, best of 3: 179 ms per loop

In [30]: %timeit multiplication2 = mul(a,a)
10 loops, best of 3: 180 ms per loop

In [31]: %timeit multiplication2 = mul(a,a)
10 loops, best of 3: 178 ms per loop

更新: 我用的是最新版的numba,就compiled it from source: '0.11.0-3-gea20d11-dirty'

我在 Fedora 19 中使用默认的 numpy 进行了测试,“1.7.1” numpy '1.6.1' 从源代码编译,链接到:

更新3 我之前的结果当然是不正确的,我在内循环中返回了 D,所以跳过了 90% 的计算。

这为 ali_m 的假设提供了更多证据,即真的很难比已经非常优化的 c 代码做得更好。

但是,如果您尝试do something more complicated,例如,

np.sqrt(((X[:, None, :] - X) ** 2).sum(-1))

我可以重现 Jake Vanderplas 得到的数字:

In [14]: %timeit pairwise_numba(X)
10000 loops, best of 3: 92.6 us per loop

In [15]: %timeit pairwise_numpy(X)
1000 loops, best of 3: 662 us per loop

看来你正在做的事情到目前为止已经被 numpy 优化了,很难做得更好。

【讨论】:

我正在使用 %% a = np.random.rand(10,5000000) \ mul(a,a) 进行计时 - 数组的创建不包括在计时计算中。您使用的是哪个版本的 numba 和 numpy? @ali_m 我在帖子中回答了。 有趣...我开始怀疑我当前的 numba/pyllvm/llvm 设置可能存在一些微妙的问题(一方面,我遇到了高于 v0.10.2 的 numba 版本的编译器错误)。我会深入研究它 - 也许它可能与 OP 正在经历的事情有关。 我也排除了时序中的数组创建。有趣的。我不知道为什么你会看到 numba 有如此巨大的改进。任何人都可以帮我弄清楚这件事吗? @ali_m 我只是复制粘贴了 ipython 中的原始代码,它把 return D 放在了 i 循环中,因此跳过了 90% 的计算,现在这更有意义了。【参考方案3】:

编辑:不要介意这个答案,我错了(见下面的评论)。


恐怕在 python 中实现比使用 numpy 更快的矩阵乘法会非常非常困难。 NumPy 通常使用内部的 fortran 库,例如 ATLAS/LAPACK,它们经过了非常好的优化。

要检查您的 NumPy 版本是否支持 LAPACK:打开终端,转到 Python 安装目录并输入:

for f in `find lib/python2.7/site-packages/numpy/* -name \*.so`; do echo $f; ldd $f;echo "\n";done | grep lapack

请注意,路径可能因您的 python 版本而异。 如果你打印了一些行,那么你肯定有 LAPACK 支持......所以在单核上实现更快的矩阵乘法将很难实现。

现在我不知道使用多核来执行矩阵乘法,所以你可能想研究一下(见 ali_m 的评论)。

【讨论】:

外部 BLAS/LAPACK 库仅与线性代数运算相关,例如 matrix 乘法。 Elementwise 乘法,就像在 OP 的示例中一样,使用用 C 代码编写的 ufunc,它是 numpy 的固有组件。话虽如此,但我的感觉是,对于像元素乘法这样简单的事情,这两种方法中的任何一种都需要超过手写 C 代码的速度。【参考方案4】:

使用 GPU。使用以下包。

gnumpy

【讨论】:

【参考方案5】:

np.multiply 的速度很大程度上依赖于大小完全相同的数组。

a = np.random.rand(80000,1)
b = np.random.rand(80000,1)

c = np.multiply(a, b)

速度快得要命,而下面的代码需要一分钟多的时间,并用光了我所有的 16 GB 内存:

a = np.squeeze(np.random.rand(80000,1))
b = np.random.rand(80000,1)

c = np.multiply(a, b)

所以我的建议是使用完全相同维度的数组。希望这对寻找如何加速元素乘法的人有用。

【讨论】:

这是因为第二个代码计算外积,而第一个代码进行元素乘法。两种截然不同的操作。第一个产生一个大小为 (80000,) 的数组,第二个产生一个大小为 (80000,80000) 的数组。

以上是关于加速python中的元素数组乘法的主要内容,如果未能解决你的问题,请参考以下文章

矩阵乘法中每个单元格的平均值,而不仅仅是 python 中的总和

为啥 FFT 的计算时间比 intel MKL 中的元素到元素乘法要短?

APL:数组的元素替换和乘法

python如何挑选矩阵中的不相领的列组成新的矩阵

在 Numpy 数组上使用 Pycuda 的 GPU 数组乘法

Python Numpy中的几个矩阵乘法