基于一维阵列的高效Numpy二维阵列构建

Question

我有一个这样的数组：

A = array([1,2,3,4,5,6,7,8,9,10])

我试着得到一个这样的数组：

B = array([[1,2,3],
          [2,3,4],
          [3,4,5],
          [4,5,6]])

其中每行（具有固定的任意宽度）移动一个。A的数组有10万个记录，我正在努力寻找一种有效的方法在纽比做这件事。目前我使用的vstack和for循环是缓慢的。有更快的方法吗？

编辑：

width = 3 # fixed arbitrary width
length = 10000 # length of A which I wish to use
B = A[0:length + 1]
for i in range (1, length):
    B = np.vstack((B, A[i, i + width + 1]))

Answer 1

事实上，有一种更有效的方法。。。使用vstack等的缺点是，您正在复制数组。

顺便说一下，这实际上和@Paul的答案是一样的，但是我发布这个只是为了更详细地解释一些事情。。。

有一种方法可以只使用视图，这样就可以复制no内存。

我直接从Erik Rigtorp's post to numpy-discussion那里借来的，他又从基思·古德曼的Bottleneck（非常有用）那里借来的。

基本技巧是直接操作strides of the array（对于一维数组）：

import numpy as np

def rolling(a, window):
    shape = (a.size - window + 1, window)
    strides = (a.itemsize, a.itemsize)
    return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)

a = np.arange(10)
print rolling(a, 3)

其中，a是您的输入数组，window是您想要的窗口长度（在您的例子中是3）。

这将产生：

[[0 1 2]
 [1 2 3]
 [2 3 4]
 [3 4 5]
 [4 5 6]
 [5 6 7]
 [6 7 8]
 [7 8 9]]

但是，在原始a和返回的数组之间绝对没有内存重复。这意味着它的速度和规模远优于其他选择。

例如（使用a = np.arange(100000)和window=3）：

%timeit np.vstack([a[i:i-window] for i in xrange(window)]).T
1000 loops, best of 3: 256 us per loop

%timeit rolling(a, window)
100000 loops, best of 3: 12 us per loop

如果我们将其推广到N维数组的最后一个轴上的“滚动窗口”，我们将得到Erik Rigtorp的“滚动窗口”函数：

import numpy as np

def rolling_window(a, window):
   """
   Make an ndarray with a rolling window of the last dimension

   Parameters
   ----------
   a : array_like
       Array to add rolling window to
   window : int
       Size of rolling window

   Returns
   -------
   Array that is a view of the original array with a added dimension
   of size w.

   Examples
   --------
   >>> x=np.arange(10).reshape((2,5))
   >>> rolling_window(x, 3)
   array([[[0, 1, 2], [1, 2, 3], [2, 3, 4]],
          [[5, 6, 7], [6, 7, 8], [7, 8, 9]]])

   Calculate rolling mean of last dimension:
   >>> np.mean(rolling_window(x, 3), -1)
   array([[ 1.,  2.,  3.],
          [ 6.,  7.,  8.]])

   """
   if window < 1:
       raise ValueError, "`window` must be at least 1."
   if window > a.shape[-1]:
       raise ValueError, "`window` is too long."
   shape = a.shape[:-1] + (a.shape[-1] - window + 1, window)
   strides = a.strides + (a.strides[-1],)
   return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)

所以，让我们看看这里发生了什么。。。操纵数组的strides看起来有点神奇，但是一旦你了解了发生的事情，就一点也不神奇了。numpy数组的步长描述了沿给定轴递增一个值所必须采取的步骤的字节大小。因此，在64位浮点的一维数组中，每个项的长度是8字节，x.strides是(8,)。

x = np.arange(9)
print x.strides

现在，如果我们将其重塑为一个2D，3x3数组，那么跨步将是(3 * 8, 8)，因为我们必须跳24字节以沿第一个轴递增一步，并跳8字节以沿第二个轴递增一步。

y = x.reshape(3,3)
print y.strides

类似地，转置与反转数组的步伐是一样的：

print y
y.strides = y.strides[::-1]
print y

很明显，数组的跨步和数组的形状是紧密相连的。如果我们改变其中一个，我们必须相应地改变另一个，否则我们就没有一个有效的内存缓冲区描述，它实际上保存了数组的值。

因此，如果要同时更改数组的形状和大小，就不能只通过设置x.strides和x.shape来完成，即使新的跨距和形状是兼容的。

这就是numpy.lib.as_strided进来的地方。它实际上是一个非常简单的函数，只需同时设置数组的跨步和形状。

它检查两者是否兼容，但不检查旧的跨步和新的形状是否兼容，如果将两者单独设置，则会发生这种情况。（实际上是通过numpy's ^{}来实现的，它允许任意类将内存缓冲区描述为numpy数组。）

所以，我们所做的就是使一个项沿着一个轴向前移动（64位数组的情况下为8字节），但是也只沿着另一个轴向前移动8字节。

换言之，在“窗口”大小为3的情况下，数组的形状为(whatever, 3)，但它没有对第二维度执行完整的3 * x.itemsize，而是只向前推进一个项目，有效地使新数组的行成为原始数组的“移动窗口”视图。

（这也意味着x.shape[0] * x.shape[1]与新数组的x.size不同。）

无论如何，希望这能让事情变得更清楚。。

Answer 2

你用的是哪种方法？

import numpy as np
A = np.array([1,2,3,4,5,6,7,8,9,10])
width = 3

np.vstack([A[i:i-len(A)+width] for i in xrange(len(A)-width)])
# needs 26.3µs

np.vstack([A[i:i-width] for i in xrange(width)]).T
# needs 13.2µs

如果你的宽度相对较低（3），并且你有一个很大的A（10000个元素），那么差别就更重要：第一个是32.4ms，第二个是44μs。

Answer 3

这个解决方案不是由python循环有效地实现的，因为在使用numpy数组时最好避免使用各种类型的检查。如果您的阵列非常高，您会注意到：

newshape = (4,3)
newstrides = (A.itemsize, A.itemsize)
B = numpy.lib.stride_tricks.as_strided(A, shape=newshape, strides=newstrides)

这将提供数组a的视图。如果您想要一个新的数组，可以进行编辑，也可以这样做，但是在最后使用.copy()。

跨步细节：

本例中的newstrides元组将是（4,4），因为数组有4字节的项，并且您希望继续在i维中的单个项步骤中单步遍历数据。第二个值“4”是指j维中的跨步（在正常的4x4数组中，它是16）。因为在这种情况下，您还需要在j维中以4字节的步骤递增从缓冲区的读取。

乔给出了一个很好的，详细的描述，并使事情非常清楚，当他说，所有这些伎俩所做的是改变步伐和形状的同时。