关键词：
【中文标题】加速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) 的数组。加速 NeDB 中的数组插入
】加速NeDB中的数组插入【英文标题】：SpeedinguparrayinsertsinNeDB【发布时间】：2020-09-0104:51:11【问题描述】：我正在制作一个应用程序，该应用程序从redditAPI中获取用户的cmets并将它们加载到本地数据库中。我正在使用NeDB作为数据... 查看详情
APL：数组的元素替换和乘法
...、3、4）分别替换为0、10、100、1000。所以我想把整个数组中的1映射到0、2映射到10、3映射到100、4映射到10 查看详情
在 Numpy 数组上使用 Pycuda 的 GPU 数组乘法
...原始的numpy逐点乘法要慢得多。我希望使用GPU获得良好的加速。zz0是complex128类型，(64,256,1 查看详情
矩阵乘法中每个单元格的平均值，而不仅仅是 python 中的总和
】矩阵乘法中每个单元格的平均值，而不仅仅是python中的总和【英文标题】：meanofeverycellinmatrixmultiplcationinsteadofjustthesuminpython【发布时间】：2021-12-1801:09:43【问题描述】：我知道使用2二维数组的numpy，它会给我矩阵乘法。以下是... 查看详情
python元素的元素乘法(代码片段)
numpy数组中元素的乘法和除法给出整数结果
】numpy数组中元素的乘法和除法给出整数结果【英文标题】：MultiplicationandDivisionofelementsinanumpyarraygivesintegerresults【发布时间】：2017-11-0402:06:55【问题描述】：importnumpyasnpA=np.array([[2,1,-1,8],[-3,-1,2,-11],[-2,1,2,-3]])B=A[1]+A[0]*(-A[1][0]/A[0] 查看详情
如何使用流从两个列表或数组乘法中查找元素对
】如何使用流从两个列表或数组乘法中查找元素对【英文标题】：Howtousestreamstofindpairsofelementsfromtwolistsorarraymultiplication【发布时间】：2017-07-0209:10:48【问题描述】：我有两个数字列表，我想找到所有可能的数字对。例如，给定列... 查看详情
Python中的KMeans：ValueError：使用序列设置数组元素
】Python中的KMeans：ValueError：使用序列设置数组元素【英文标题】：KMeansinPython:ValueError:settinganarrayelementwithasequence【发布时间】：2016-06-1205:11:51【问题描述】：我正在尝试使用numpy和sklearn在Python中执行kmeans聚类。我有一个包含45列... 查看详情
加速必须遍历整个列表的 Python 代码
】加速必须遍历整个列表的Python代码【英文标题】：SpeedingupPythoncodethathastogothroughentirelist【发布时间】：2017-08-1403:56:30【问题描述】：我有一个问题需要（至少可以肯定）遍历整个列表来解决。问题是找出列表中最多的连续数字... 查看详情
为啥 FFT 的计算时间比 intel MKL 中的元素到元素乘法要短？
】为啥FFT的计算时间比intelMKL中的元素到元素乘法要短？【英文标题】：WhycalculatingtimeofFFTisshorterthanelement-to-elementmultiplicationinintelMKL?为什么FFT的计算时间比intelMKL中的元素到元素乘法要短？【发布时间】：2019-03-1902:51:14【问题描... 查看详情
AVX中的6元素双精度向量矩阵向量乘法
】AVX中的6元素双精度向量矩阵向量乘法【英文标题】：6elementdoubleprecisionvectormatrixvectormultiplyinAVX【发布时间】：2013-05-1919:47:25【问题描述】：我需要以双精度进行以下操作：数字表示值在内存中的存储方式。我想用AVX来实现这... 查看详情
如何从 Python 中的数组中选择随机元素？ [复制]
】如何从Python中的数组中选择随机元素？[复制]【英文标题】：HowdoIselectarandomelementfromanarrayinPython?[duplicate]【发布时间】：2010-11-0617:17:19【问题描述】：我用谷歌搜索的第一个例子没有用。这应该是微不足道的，对吧？【问题讨... 查看详情
检查字符串是不是包含Python中数组中的多个元素
】检查字符串是不是包含Python中数组中的多个元素【英文标题】：CheckifstringcontainsmorethanoneelementfromarrayinPython检查字符串是否包含Python中数组中的多个元素【发布时间】：2016-06-1706:09:14【问题描述】：我在我的项目中使用正则表... 查看详情
什么是数组的维度，python的ndim的使用
...就是几维。numpy中直接用*即可表示数与向量的乘法，参考python2.7的一个例子：inportnumpyasnp a=np.array([1,2,3,4])#向量 b=5#数 printa*b ++++++++++++ [5,10,15,20] NumPy数组的下标从0开始。 同一个NumPy数组中... 查看详情
将文件逐行读入Python中的数组元素[重复]
】将文件逐行读入Python中的数组元素[重复]【英文标题】：ReadingafilelinebylineintoelementsofanarrayinPython[duplicate]【发布时间】：2013-04-1919:58:03【问题描述】：所以在Ruby中我可以做到以下几点：testsite_array=Array.newy=0File.open(\'topsites.txt\').... 查看详情
Armadillo C ++中的元素方式向量或矩阵乘法
】ArmadilloC++中的元素方式向量或矩阵乘法【英文标题】：elementwisevectorormatrixmultiplicationinArmadilloC++【发布时间】：2014-11-0622:01:19【问题描述】：#include<iostream>#include<armadillo>usingnamespacestd;usingnamespacearma;intmain()vecx=(1 查看详情
Python中的高效数组替换
】Python中的高效数组替换【英文标题】：EfficientArrayreplacementinPython【发布时间】：2011-11-1305:58:48【问题描述】：我想知道在给定一些标准的情况下，用数组中的其他随机元素替换数组中的元素的最有效方法是什么。更具体地说，... 查看详情
如何将所有数组的元素添加到python中的一个列表中[重复]
】如何将所有数组的元素添加到python中的一个列表中[重复]【英文标题】：howtoaddallarray\'selementstoonelistinpython[duplicate]【发布时间】：2013-04-2508:44:36【问题描述】：有一个看起来像这样的二维数组：myarray=[[\'jacob\',\'mary\'],[\'jack\',\'... 查看详情