从一维数组高效构建二维Numpy数组

38 投票
7 回答
14600 浏览
提问于 2025-04-16 11:23

我有一个这样的数组:

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有一万条记录,我想在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]))

7 个回答

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.4毫秒,而第二种方法只需要44微秒。

11

这个解决方案在用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字节的步长从缓冲区读取数据。

Joe给出了一个很好的详细描述,他清楚地说明了这个技巧的作用就是同时改变步长和形状。

61

其实,还有一种更高效的方法来实现这个功能……使用 vstack 等方法的缺点是,你会复制一份数组。

顺便提一下,这个方法和 @Paul 的回答基本上是一样的,但我发这个是想更详细地解释一下……

有一种方法可以仅通过视图来实现,这样就不会重复占用内存。

我直接借用自 Erik Rigtorp 在 numpy-discussion 的帖子,他又是从 Keith Goodman 的 Bottleneck 中借来的(这个库非常有用!)。

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

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.stridesx.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,对于你的新数组来说。)

总之,希望这能让事情稍微清楚一些……

撰写回答