如何从nxn数组中提取mxm子矩阵 (n>m)？

Question

我想对一个NumPy的nxn数组进行切片。我想从这个数组中提取一些任意的m行和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

我觉得 x[[1,3]][:,[1,3]] 这个写法并不难懂。如果你想让你的意图更清晰，可以这样写：

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

我对切片不是特别专业，但通常来说，如果你想从一个数组中切出一部分，并且这些值是连续的，你会得到一个视图，这个视图的步长值会改变。

举个例子，在你的输入33和34中，虽然你得到的是一个2x2的数组，但它的步长是4。因此，当你索引下一行时，指针会移动到内存中的正确位置。

显然，这种机制在处理索引数组时就不太适用了。所以，numpy不得不进行复制。毕竟，很多其他的矩阵运算函数都依赖于大小、步长和连续的内存分配。

Answer 2

要回答这个问题，我们需要了解在Numpy中如何对多维数组进行索引。首先，假设你有一个数组x，这个数组里面会包含从0到15的16个递增的整数。如果你想访问其中的一个元素，比如x[i,j]，Numpy需要计算这个元素在内存中的位置，也就是它距离数组开头有多远。这是通过计算i*x.shape[1]+j来实现的（然后再乘以一个整数的大小，得到实际的内存偏移量）。

如果你通过基本切片提取一个子数组，比如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]（这已经包括了整数大小的因素）：

x.strides
(16, 4)

当像上面那样提取y时，Numpy并不会创建一个新的缓冲区，而是会创建一个新的数组对象来引用同一个缓冲区（否则y就会和x完全相同）。新的数组对象会有不同的形状，可能在缓冲区中的起始偏移量也不同，但会和x共享步幅（至少在这种情况下是这样的）：

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

这样，计算y[i,j]的内存偏移量就能得到正确的结果。

但是，对于像z=x[[1,3]]这样的情况，Numpy该怎么做呢？如果使用原始缓冲区，步幅机制就无法正确索引。理论上，Numpy可以添加一些比步幅更复杂的机制，但这样会使得访问元素变得相对昂贵，这有点违背了数组的初衷。此外，视图也不再是一个轻量级的对象了。

这个问题在Numpy的索引文档中有详细的讨论。

哦，差点忘了你的实际问题：这是如何让多个列表的索引按预期工作的方法：

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

这是因为索引数组会被广播到一个共同的形状。当然，对于这个特定的例子，你也可以使用基本切片来实现：

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

Answer 3

正如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])]来做和我第一个例子一样的事情。这样可以省去输入那么多额外的括号。