NumPy数组加权部分和

Question

我有一个numpy数组，例如[a,b,c,d,e,...]，想要计算一个看起来像[x*a+y*b, x*b+y*c, x*c+y*d,...]的数组。我的想法是首先将原始数组分成类似[[a,b],[b,c],[c,d],[d,e],...]的东西，然后使用np.average指定轴和权重（在我的情况下为x+y=1），甚至使用np.dot攻击这个数组。不幸的是，我不知道如何创建这样的[a,b],[b,c],...对的数组。非常感谢任何帮助，或者即使是实现主要任务的完全不同的想法 :-)

- Simon Righley

---

我认为这个嵌套生成器是有效的：[(v1*x + v2*y) for v1, v2 in [arr[i:i+2] for i in xrange(len(arr)-1)]] - jon

Answer 1

最快、最简单的方法是手动提取数组的两个切片并将它们相加：

>>> arr = np.arange(5)
>>> x, y = 10, 1
>>> x*arr[:-1] + y*arr[1:]
array([ 1, 12, 23, 34])

如果你想将其推广到三元组、四元组等，这可能会变得很麻烦。但是，你可以使用as_strided以更通用的形式从原始数组创建成对数组：

>>> from numpy.lib.stride_tricks import as_strided

>>> arr_pairs = as_strided(arr, shape=(len(arr)-2+1,2), strides=arr.strides*2)
>>> arr_pairs
array([[0, 1],
       [1, 2],
       [2, 3],
       [3, 4]])

当然，使用as_strided的好处就是，就像数组切片一样，它不涉及数据复制，只是改变了内存的视图方式，因此创建这个数组的成本几乎为零。

现在可能最快的方法是使用np.dot:

>>> xy = [x, y]
>>> np.dot(arr_pairs, xy)
array([ 1, 12, 23, 34])

- Jaime

Answer 2

这似乎是一个相关问题，与 numpy.correlate 相关。

a
Out[61]: array([0, 1, 2, 3, 4, 5, 6, 7])

b
Out[62]: array([1, 2])

np.correlate(a,b,mode='valid')
Out[63]: array([ 2,  5,  8, 11, 14, 17, 20])

根据数组大小和BLAS点积的不同，性能会有很大差异。具体情况因人而异：

arr = np.random.rand(1E6)

b = np.random.rand(2)

np.allclose(jamie_dot(arr,b),np.convolve(arr,b[::-1],mode='valid'))
True

%timeit jamie_dot(arr,b)
100 loops, best of 3: 16.1 ms per loop

%timeit np.correlate(arr,b,mode='valid')
10 loops, best of 3: 28.8 ms per loop

使用英特尔 MKL BLAS 和 8 个核心，对于大多数实现来说，np.correlate 可能会更快。

此外，@Jamie 的帖子中有一个有趣的观察：

%timeit b[0]*arr[:-1] + b[1]*arr[1:]
100 loops, best of 3: 8.43 ms per loop

他的评论还将使用np.convolve(a,b[::-1],mode=valid)的复杂语法简化为更简单的correlate语法。

- Daniel

---

这实际上更像是一个correlate问题，但是没错，+1。实际上，如果你查看源代码，np.convolve(a, b)调用了np.correlate(a, b[::-1])。 - Jaime

@Jamie 哦，真不错，我从没注意到这个函数——知道你可以做到这种事情而不需要卷积技巧是很好的。 - Daniel

@Ophion，整洁，非常有见地！ - Simon Righley

Answer 3

如果您有一个小数组，我会创建一个偏移副本：

shifted_array=numpy.append(original_array[1:],0)
result_array=x*original_array+y*shifted_array

这里你需要把数组存储两次在内存中，因此这种解决方案非常浪费内存，但你可以摆脱for循环。

如果你有大型数组，你确实需要一个循环（但更好的选择是列表推导式）：

result_array=[x*original_array[i]+y*original_array[i+1] for i in xrange(len(original_array)-1)]

它给你与Python列表相同的结果，除了最后一项，无论如何都应该特别处理。

根据一些随机试验，在小于2000个项目的数组中，第一个解决方案似乎比第二个更快，但即使对于相对较小的数组（我的PC上有几万个），也会遇到MemoryError问题。

因此，通常使用列表推导式，但如果您确定只在小型（最大1-2千）数组上运行，则有更好的机会。

创建像[[a,b],[b,c],[c,d],[d,e],...]这样的新列表既浪费内存又浪费时间，因为您还需要for循环（或类似）来创建它，并且必须将每个旧值两次存储在新数组中，因此您最终将存储原始数组三次。

- leeladam

Answer 4

另一种方法是在数组中创建正确的对应关系：a = np.array([a,b,c,d,e,...])，根据数组的大小进行重塑：b = np.array([x, y, ...])，然后利用numpy的广播规则来进行操作。

a = np.arange(8) 
b = np.array([1, 2])

a = a.repeat(2)[1:-1]
ans = a.reshape(-1, b.shape[0]).dot(b)

时间（在我的电脑上）：

@Ophion's solution:
# 100000 loops, best of 3: 4.67 µs per loop

This solution:
# 100000 loops, best of 3: 9.78 µs per loop

所以，它比较慢。@Jaime的解决方案更好，因为它不像这个复制数据。

- Saullo G. P. Castro