问题：在NumPy数组的每个单元中高效评估函数
给定一个NumPy数组A，将相同的函数f应用于每个单元的最快/最有效的方法是什么？
假设我们将分配给A（I，J）的F（A（I，J）） 。
函数f没有二进制输出，因此mask（ing）操作将无济于事。
“显而易见的”双循环迭代（通过每个单元）是否是最佳解决方案？
回答 0
您可以对函数进行矢量化处理，然后在每次需要时将其直接应用于Numpy数组：
import numpy as np
def f(x):
return x * x + 3 * x - 2 if x > 0 else x * 5 + 8
f = np.vectorize(f) # or use a different name if you want to keep the original f
result_array = f(A) # if A is your Numpy array向量化时最好直接指定一个显式输出类型：
f = np.vectorize(f, otypes=[np.float])回答 1
一个类似的问题是：将NumPy数组映射到适当的位置。如果可以为f（）找到一个ufunc，则应使用out参数。
回答 2
如果您使用数字和f(A(i,j)) = f(A(j,i))，则可以使用scipy.spatial.distance.cdist将f定义为A(i)和之间的距离A(j)。
回答 3
我相信我找到了更好的解决方案。将函数更改为python通用函数的想法（请参阅文档），可以在后台进行并行计算。
一个人可以用ufuncC 编写自己的自定义脚本，这肯定会更有效，也可以通过调用np.frompyfunc内置工厂方法来编写。经过测试，此方法比np.vectorize：
f = lambda x, y: x * y
f_arr = np.frompyfunc(f, 2, 1)
vf = np.vectorize(f)
arr = np.linspace(0, 1, 10000)
%timeit f_arr(arr, arr) # 307ms
%timeit f_arr(arr, arr) # 450ms我还测试了较大的样本，并且改进成比例。有关其他方法的性能比较，请参阅这篇文章
回答 4
当2d数组（或nd数组）为C或F连续时，将函数映射到2d数组的任务实际上与将函数映射到1d数组的任务相同-我们只是必须以这种方式查看它，例如通过np.ravel(A,'K')。
例如，这里讨论了1d阵列的可能解决方案。
但是，当2d数组的内存不连续时，情况会稍微复杂一些，因为如果轴处理顺序错误，则希望避免可能发生的高速缓存未命中。
Numpy已经有机器以最佳顺序加工轴。使用这种机械的一种可能性是np.vectorize。但是，numpy的文档np.vectorize指出“主要是为了方便而不是为了性能而提供”-慢速python函数保持慢速python函数以及所有相关的开销！另一个问题是其巨大的内存消耗-例如，参见此SO-post。
当想要具有C函数的性能但要使用numpy的机器时，一个好的解决方案是使用numba创建ufunc，例如：
# runtime generated C-function as ufunc
import numba as nb
@nb.vectorize(target="cpu")
def nb_vf(x):
return x+2*x*x+4*x*x*x它很容易击败，np.vectorize但是当执行相同的功能作为numpy-array乘法/加法时，即
# numpy-functionality
def f(x):
return x+2*x*x+4*x*x*x
# python-function as ufunc
import numpy as np
vf=np.vectorize(f)
vf.__name__="vf"有关时间测量代码，请参见此答案的附录：
Numba的版本（绿色）比python函数（即np.vectorize）快约100倍，这并不奇怪。但这也比numpy功能快约10倍，因为numbas版本不需要中间数组，因此可以更有效地使用缓存。
尽管numba的ufunc方法是可用性和性能之间的良好折衷，但它仍然不是我们能做的最好的选择。然而，没有灵丹妙药或最适合任何任务的方法-人们必须了解什么是局限性以及如何减轻它们。
例如，对于超越函数（例如exp，sin，cos）numba不提供超过任何优点numpy的的np.exp（有没有创建临时数组-高速化的主要来源）。但是，我的Anaconda安装使用Intel的VML来处理大于8192的矢量-如果内存不连续，则无法执行。因此，最好将元素复制到连续内存中，以便能够使用英特尔的VML：
import numba as nb
@nb.vectorize(target="cpu")
def nb_vexp(x):
return np.exp(x)
def np_copy_exp(x):
copy = np.ravel(x, 'K')
return np.exp(copy).reshape(x.shape) 为了公平起见，我关闭了VML的并行化功能（请参见附录中的代码）：
可以看到，一旦VML开始运行，复制的开销就远远超过了补偿。但是，一旦数据对于L3高速缓存而言变得太大，则优势就变得微不足道了，因为任务再次变得与内存带宽绑定。
在另一方面，numba可以使用英特尔的SVML为好，在解释这个帖子：
from llvmlite import binding
# set before import
binding.set_option('SVML', '-vector-library=SVML')
import numba as nb
@nb.vectorize(target="cpu")
def nb_vexp_svml(x):
return np.exp(x)并使用具有并行化功能的VML生成：
numba的版本开销较小，但是对于某些大小，VML甚至比SVML都要高，尽管有额外的复制开销-这并不奇怪，因为numba的ufunc没有并行化。
清单：
A.多项式函数的比较：
import perfplot
perfplot.show(
setup=lambda n: np.random.rand(n,n)[::2,::2],
n_range=[2**k for k in range(0,12)],
kernels=[
f,
vf,
nb_vf
],
logx=True,
logy=True,
xlabel='len(x)'
) B.比较exp：
import perfplot
import numexpr as ne # using ne is the easiest way to set vml_num_threads
ne.set_vml_num_threads(1)
perfplot.show(
setup=lambda n: np.random.rand(n,n)[::2,::2],
n_range=[2**k for k in range(0,12)],
kernels=[
nb_vexp,
np.exp,
np_copy_exp,
],
logx=True,
logy=True,
xlabel='len(x)',
)When the 2d-array (or nd-array) is C- or F-contiguous, then this task of mapping a function onto a 2d-array is practically the same as the task of mapping a function onto a 1d-array – we just have to view it that way, e.g. via np.ravel(A,'K').
Possible solution for 1d-array have been discussed for example here.
However, when the memory of the 2d-array isn’t contiguous, then the situation a little bit more complicated, because one would like to avoid possible cache misses if axis are handled in wrong order.
Numpy has already a machinery in place to process axes in the best possible order. One possibility to use this machinery is np.vectorize. However, numpy’s documentation on np.vectorize states that it is “provided primarily for convenience, not for performance” – a slow python function stays a slow python function with the whole associated overhead! Another issue is its huge memory-consumption – see for example this SO-post.
When one wants to have a performance of a C-function but to use numpy’s machinery, a good solution is to use numba for creation of ufuncs, for example:
# runtime generated C-function as ufunc
import numba as nb
@nb.vectorize(target="cpu")
def nb_vf(x):
return x+2*x*x+4*x*x*x
It easily beats np.vectorize but also when the same function would be performed as numpy-array multiplication/addition, i.e.
# numpy-functionality
def f(x):
return x+2*x*x+4*x*x*x
# python-function as ufunc
import numpy as np
vf=np.vectorize(f)
vf.__name__="vf"
See appendix of this answer for time-measurement-code:
Numba’s version (green) is about 100 times faster than the python-function (i.e. np.vectorize), which is not surprising. But it is also about 10 times faster than the numpy-functionality, because numbas version doesn’t need intermediate arrays and thus uses cache more efficiently.
While numba’s ufunc approach is a good trade-off between usability and performance, it is still not the best we can do. Yet there is no silver bullet or an approach best for any task – one has to understand what are the limitation and how they can be mitigated.
For example, for transcendental functions (e.g. exp, sin, cos) numba doesn’t provide any advantages over numpy’s np.exp (there are no temporary arrays created – the main source of the speed-up). However, my Anaconda installation utilizes Intel’s VML for vectors bigger than 8192 – it just cannot do it if memory is not contiguous. So it might be better to copy the elements to a contiguous memory in order to be able to use Intel’s VML:
import numba as nb
@nb.vectorize(target="cpu")
def nb_vexp(x):
return np.exp(x)
def np_copy_exp(x):
copy = np.ravel(x, 'K')
return np.exp(copy).reshape(x.shape)
For the fairness of the comparison, I have switched off VML’s parallelization (see code in the appendix):
As one can see, once VML kicks in, the overhead of copying is more than compensated. Yet once data becomes too big for L3 cache, the advantage is minimal as task becomes once again memory-bandwidth-bound.
On the other hand, numba could use Intel’s SVML as well, as explained in this post:
from llvmlite import binding
# set before import
binding.set_option('SVML', '-vector-library=SVML')
import numba as nb
@nb.vectorize(target="cpu")
def nb_vexp_svml(x):
return np.exp(x)
and using VML with parallelization yields:
numba’s version has less overhead, but for some sizes VML beats SVML even despite of the additional copying overhead – which isn’t a bit surprise as numba’s ufuncs aren’t parallelized.
Listings:
A. comparison of polynomial function:
import perfplot
perfplot.show(
setup=lambda n: np.random.rand(n,n)[::2,::2],
n_range=[2**k for k in range(0,12)],
kernels=[
f,
vf,
nb_vf
],
logx=True,
logy=True,
xlabel='len(x)'
)
B. comparison of exp:
import perfplot
import numexpr as ne # using ne is the easiest way to set vml_num_threads
ne.set_vml_num_threads(1)
perfplot.show(
setup=lambda n: np.random.rand(n,n)[::2,::2],
n_range=[2**k for k in range(0,12)],
kernels=[
nb_vexp,
np.exp,
np_copy_exp,
],
logx=True,
logy=True,
xlabel='len(x)',
)
回答 5
以上所有答案比较起来都不错，但是如果您需要使用自定义函数进行映射，并且拥有numpy.ndarray，则需要保留数组的形状。
我只比较了两个，但它将保留的形状ndarray。我已将具有100万个条目的数组用于比较。在这里我使用平方函数。我正在介绍n维数组的一般情况。对于二维，只需制作iter2D。
import numpy, time
def A(e):
return e * e
def timeit():
y = numpy.arange(1000000)
now = time.time()
numpy.array([A(x) for x in y.reshape(-1)]).reshape(y.shape)
print(time.time() - now)
now = time.time()
numpy.fromiter((A(x) for x in y.reshape(-1)), y.dtype).reshape(y.shape)
print(time.time() - now)
now = time.time()
numpy.square(y)
print(time.time() - now)输出量
>>> timeit()
1.162431240081787 # list comprehension and then building numpy array
1.0775556564331055 # from numpy.fromiter
0.002948284149169922 # using inbuilt function在这里你可以清楚地看到 numpy.fromiter用户平方函数，可以使用任何选择。如果你的功能是依赖于i, j 那就是数组的索引，迭代上数组的大小一样for ind in range(arr.size)，用numpy.unravel_index得到i, j, ..基于阵列的您1D指数和形状numpy.unravel_index
这个答案是受到我对其他问题的回答的启发 这里
