在Numpy中做‘A[tuple(B.T)]的更快的方法
在Numpy中做‘A[tuple(B.T)]的更快的方法
提问于 2015-08-15 19:17:32
我有一个M by N数组I，其中的每一行都是一个索引，一个N维数组A。我想要一个向量化表达式从A获得M索引值的一维数组，我发现A[tuple(I.T)]做的是正确的，但是分析表明，尽管被向量化，它还是非常昂贵的。它也不是特别优雅或“自然”，A[I]和A[I.T]做了完全不同的事情。
怎样才是正确的方法？它也应该适用于任务分配，例如
A[tuple(I.T)] = 1回答 2
Stack Overflow用户
回答已采纳
发布于 2015-08-15 20:28:45
我想你说的是：
In [398]: A=np.arange(24).reshape(4,6)
In [401]: I=np.array([[0,1],[1,2],[3,4],[0,0],[2,5]])
In [402]: tuple(I.T)
Out[402]: (array([0, 1, 3, 0, 2]), array([1, 2, 4, 0, 5]))
In [403]: A[tuple(I.T)]
Out[403]: array([ 1, 8, 22, 0, 17])这是http://docs.scipy.org/doc/numpy/reference/arrays.indexing.html#purely-integer-array-indexing -纯整数数组高级索引。
这总是比返回视图的基本索引慢。基本索引选择连续的数据块，或可以大步选择的值。用你的索引是不可能的。
看看一些时间安排：
In [404]: timeit tuple(I.T)
100000 loops, best of 3: 3.4 µs per loop
In [405]: timeit A[tuple(I.T)]
100000 loops, best of 3: 10 µs per loop
In [406]: %%timeit i,j=tuple(I.T)
.....: A[i,j]
.....:
100000 loops, best of 3: 4.86 µs per loop构建元组需要花费大约1/3的时间。i,j=I.T只是稍微快了一点。但索引是最大的一部分。
A[i,j]和A[(i,j)] (和A.__getitem__((i,j))一样)是一样的。因此，将I.T封装在tuple中只会产生2个索引数组，每个维度都有一个。
在该数组的扁平版本上索引速度更快：
In [420]: J= np.ravel_multi_index(tuple(I.T),A.shape)
In [421]: J
Out[421]: array([ 1, 8, 22, 0, 17], dtype=int32)
In [422]: A.flat[J]
Out[422]: array([ 1, 8, 22, 0, 17])
In [425]: timeit A.flat[J]
1000000 loops, best of 3: 1.56 µs per loop
In [426]: %%timeit
.....: J= np.ravel_multi_index(tuple(I.T),A.shape)
.....: A.flat[J]
.....:
100000 loops, best of 3: 11.2 µs per loop因此，能够预先计算和重用索引将节省您的时间，但是从A中选择一组单独的值将需要额外的时间，这是不可能的。
只是为了好玩，比较索引A和每一行I所需的时间。
In [442]: timeit np.array([A[tuple(i)] for i in I])
100000 loops, best of 3: 17.3 µs per loop
In [443]: timeit np.array([A[i,j] for i,j in I])
100000 loops, best of 3: 15.7 µs per loopStack Overflow用户
发布于 2015-08-16 08:19:09
你可以用另一种方式使用linear indexing，就像-
def ravel_einsum(A,I):
# Get A's shape and calculate cummulative dimensions based on it
shp = np.asarray(A.shape)
cumdims = np.append(1,shp[::-1][:-1].cumprod())[::-1]
# Use linear indexing of A to extract elements from A corresponding
# to linear indexing of it with I
return A.ravel()[np.einsum('ij,j->i',I,cumdims)]运行时测试
判例1：
In [84]: # Inputs
...: A = np.random.randint(0,10,(3,5,2,4,5,2,6,8,2,5,3,4,3))
...: I = np.mod(np.random.randint(0,10,(5,A.ndim)),A.shape)
...:
In [85]: %timeit A[tuple(I.T)]
10000 loops, best of 3: 27.7 µs per loop
In [86]: %timeit ravel_einsum(A,I)
10000 loops, best of 3: 48.3 µs per loop案例2：
In [87]: # Inputs
...: A = np.random.randint(0,10,(3,5,4,2))
...: I = np.mod(np.random.randint(0,5,(10000,A.ndim)),A.shape)
...:
In [88]: %timeit A[tuple(I.T)]
1000 loops, best of 3: 357 µs per loop
In [89]: %timeit ravel_einsum(A,I)
1000 loops, best of 3: 240 µs per loop案例3：
In [90]: # Inputs
...: A = np.random.randint(0,10,(30,50,40,20))
...: I = np.mod(np.random.randint(0,50,(5000,A.ndim)),A.shape)
...:
In [91]: %timeit A[tuple(I.T)]
1000 loops, best of 3: 220 µs per loop
In [92]: %timeit ravel_einsum(A,I)
10000 loops, best of 3: 168 µs per loop页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/32028402
复制相关文章
相似问题
腾讯云开发者
