切片NumPy 2d数组,或者如何从nxn数组(n> m)中提取mxm子矩阵？

Question

我想切片NumPy nxn数组.我想提取该数组的m行和列的任意选择(即行数/列数没有任何模式),使其成为一个新的mxm数组.对于这个例子,让我们说数组是4x4,我想从中提取一个2x2数组.

这是我们的数组:

from numpy import *
x = range(16)
x = reshape(x,(4,4))

print x
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]
 [12 13 14 15]]

要删除的行和列是相同的.最简单的情况是当我想要提取一个位于开头或结尾的2x2子矩阵时,即:

In [33]: x[0:2,0:2]
Out[33]: 
array([[0, 1],
       [4, 5]])

In [34]: x[2:,2:]
Out[34]: 
array([[10, 11],
       [14, 15]])

但是,如果我需要删除另一行/列的混合物怎么办？如果我需要删除第一行和第三行/行,从而提取子矩阵[[5,7],[13,15]]怎么办？行/行可以有任何组合.我读到某个地方,我只需要使用行和列的索引数组/列表索引我的数组,但这似乎不起作用:

In [35]: x[[1,3],[1,3]]
Out[35]: array([ 5, 15])

我发现了一种方法,即:

    In [61]: x[[1,3]][:,[1,3]]
Out[61]: 
array([[ 5,  7],
       [13, 15]])

第一个问题是它几乎不可读,尽管我可以忍受它.如果有人有更好的解决方案,我当然希望听到它.

另一件事是我在一个论坛上读到,使用数组索引数组强制NumPy制作所需数组的副本,因此当处理大型数组时,这可能会成为一个问题.为什么这样/这个机制如何运作？

Answer 1

要回答这个问题,我们必须看看如何在Numpy中对多维数组进行索引.我们先说你有x问题中的数组.分配给的缓冲区x将包含从0到15的16个升序整数.如果访问一个元素,比如说x[i,j],NumPy必须计算出该元素相对于缓冲区开头的内存位置.这是通过有效计算i*x.shape[1]+j(并乘以int的大小来获得实际的内存偏移量)来完成的.

如果通过基本切片提取子y = x[0:2,0:2]数组,则生成的对象将与之共享底层缓冲区x.但如果你访问会发生什么y[i,j]？NumPy不能i*y.shape[1]+j用来计算数组中的偏移量,因为属于的数据y在内存中不是连续的.

NumPy通过引入步幅解决了这个问题.在计算用于访问的内存偏移量时x[i,j],实际计算的是i*x.strides[0]+j*x.strides[1](并且这已经包括int的大小因子):

x.strides
(16, 4)

当y像上面那样提取时,NumPy不会创建一个新的缓冲区,但它会创建一个引用相同缓冲区的新数组对象(否则y就等于x.)新的数组对象将具有不同的形状然后x可能是一个不同的开始偏移到缓冲区,但将与x(至少在这种情况下)共享步幅:

y.shape
(2,2)
y.strides
(16, 4)

这样,计算内存偏移量y[i,j]将产生正确的结果.

但NumPy应该做些什么z=x[[1,3]]呢？如果使用原始缓冲区,则strides机制将不允许正确索引z.理论上, NumPy 可以添加一些比步幅更复杂的机制,但是这会使元素访问相对昂贵,从某种程度上违背了数组的整个想法.此外,视图不再是一个非常轻量级的对象.

有关索引的NumPy文档中对此进行了深入介绍.

哦,几乎忘记了你的实际问题:以下是如何使多个列表的索引按预期工作:

x[[[1],[3]],[1,3]]

这是因为索引数组被广播到一个共同的形状.当然,对于这个特定示例,您还可以使用基本切片:

x[1::2, 1::2]

Answer 2

正如Sven所提到的,x[[[0],[2]],[1,3]]将返回与1和3列匹配的0和2行,同时x[[0,2],[1,3]]将返回数组中的值x [0,1]和x [2,3].

我做的第一个例子有一个有用的功能,numpy.ix_.你可以用我的第一个例子做同样的事情x[numpy.ix_([0,2],[1,3])].这可以使您不必输入所有这些额外的括号.

Answer 3

我认为这x[[1,3]][:,[1,3]]几乎不可读.如果您想更清楚自己的意图,可以这样做:

a[[1,3],:][:,[1,3]]

我不是切片方面的专家,但通常情况下,如果您尝试切片到数组并且值是连续的,那么您将返回一个视图,其中步幅值已更改.

例如,在输入33和34中,虽然你得到一个2x2数组,但是步幅是4.因此,当你索引下一行时,指针移动到内存中的正确位置.

显然,这种机制并不适用于一系列指数.因此,numpy必须制作副本.毕竟,许多其他矩阵数学函数依赖于大小,步幅和连续的内存分配.

Answer 4

如果要跳过每隔一行和每隔一列,那么可以使用基本切片来执行此操作:

In [49]: x=np.arange(16).reshape((4,4))
In [50]: x[1:4:2,1:4:2]
Out[50]: 
array([[ 5,  7],
       [13, 15]])

这将返回一个视图,而不是数组的副本.

In [51]: y=x[1:4:2,1:4:2]

In [52]: y[0,0]=100

In [53]: x   # <---- Notice x[1,1] has changed
Out[53]: 
array([[  0,   1,   2,   3],
       [  4, 100,   6,   7],
       [  8,   9,  10,  11],
       [ 12,  13,  14,  15]])

同时z=x[(1,3),:][:,(1,3)]采用先进的索引,从而返回一个副本:

In [58]: x=np.arange(16).reshape((4,4))
In [59]: z=x[(1,3),:][:,(1,3)]

In [60]: z
Out[60]: 
array([[ 5,  7],
       [13, 15]])

In [61]: z[0,0]=0

请注意,x保持不变:

In [62]: x
Out[62]: 
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11],
       [12, 13, 14, 15]])

如果要选择任意行和列,则不能使用基本切片.你必须使用高级索引,使用类似的x[rows,:][:,columns]地方rows和columns序列.这当然会为您提供原始数组的副本,而不是视图.这是人们应该期待的,因为numpy数组使用连续的内存(具有常量步幅),并且无法生成具有任意行和列的视图(因为这将需要非常量的步幅).

Answer 5

使用numpy,您可以为索引的每个组件传递一个切片 - 因此,x[0:2,0:2]上面的示例可以正常工作.

如果您只想均匀地跳过列或行,则可以传递包含三个组件的切片(即启动,停止,步骤).

再次,对于上面的例子:

>>> x[1:4:2, 1:4:2]
array([[ 5,  7],
       [13, 15]])

基本上是:第一维中的切片,索引为1的开始,当索引等于或大于4时停止,并在每次传递中向索引添加2.第二个维度也是如此.再说一遍:这只适用于不断的步骤.

你必须在内部做一些完全不同的语法 - x[[1,3]][:,[1,3]]实际上做的是创建一个新的数组,只包括原始数组中的第1行和第3行(完成x[[1,3]]部分),然后重新切片 - 创建第三个数组 - 仅包括上一个数组的第1列和第3列.