1D阵列的高效Numpy 2D阵列构造

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的数组是10k记录长,我试图在Numpy中找到一种有效的方法.目前我正在使用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的答案完全相同,但我发布的内容只是为了更详细地解释一下......

有一种方法可以只使用视图来执行此操作,以便不会重复内存.

我直接从Erik Rigtorp的帖子中借用这个来讨论numpy讨论,后者又从Keith Goodman的瓶颈中借用它(这非常有用!).

基本技巧是直接操纵数组的步幅(对于一维数组):

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来实现的__array_interface__,它允许任意类将内存缓冲区描述为numpy数组.)

因此,我们所做的就是沿着一个轴向前移动一个项目(在64位阵列的情况下为8个字节),但也只沿另一个轴向前移动8个字节.

换句话说,在"窗口"大小为3的情况下,阵列的形状为(whatever, 3),但不是3 * x.itemsize为第二个维度踩一个完整的,它只向前迈出一个项目,有效地使新数组的行"移动"窗口"查看原始数组.

(这也意味着x.shape[0] * x.shape[1]与x.size新阵列不同.)

无论如何,希望这会让事情变得更加清晰......

Answer 2

这个解决方案没有通过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字节步长增加缓冲区的读取.

乔给出了一个很好的,详细的描述,当他说所有这些诀窍都是同时改变步幅和形状时,事情变得清晰.