python-2.7 - Numpy - 考虑偏移量分割矩阵

Question

给定一个 m x n矩阵我想把它分成正方形a x a (a = 3 或 a = 4)任意偏移量矩阵(最小偏移量 = 1，最大偏移量 = block 大小)，如 Mathematica 的 Partition 功能做:

例如，给定一个 4 x 4 矩阵 A像

1  2  3  4 
5  6  7  8
9  10 11 12
13 14 15 16

如果我给出 3 x 3 block 并且偏移量 = 1，我想得到 4 个矩阵:

1  2  3 
5  6  7 
9  10 11

2  3  4 
6  7  8
10 11 12

5  6  7 
9  10 11
13 14 15

6  7  8
10 11 12
14 15 16

如果矩阵A是 A = np.arange(1, 37).reshape((6,6))我使用偏移量 = 3 的 3 x 3 block ，我想输出这些 block :

1  2  3
7  8  9
3 14 15

 4  5  6
10 11 12
16 17 18

19 20 21
25 26 27
31 32 33

22 23 24
28 29 30
34 35 36

我可以将矩阵 A 作为列表列表，我认为我不需要 NumPy 的功能。我很惊讶array_split也不是 numpy.split开箱即用地提供这个偏移选项，在纯 Python 中使用切片进行编码是否更直接，还是我应该研究 NumPy 的进步？我希望代码高度易读。

Answer 1

正如您所暗示的，有一种方法可以使用 strides

In [900]: M = np.lib.stride_tricks.as_strided(A, shape=(2,2,3,3), strides=(16,4,16,4))
In [901]: M
Out[901]: 
array([[[[ 1,  2,  3],
         [ 5,  6,  7],
         [ 9, 10, 11]],

        [[ 2,  3,  4],
         [ 6,  7,  8],
         [10, 11, 12]]],


       [[[ 5,  6,  7],
         [ 9, 10, 11],
         [13, 14, 15]],

        [[ 6,  7,  8],
         [10, 11, 12],
         [14, 15, 16]]]])
In [902]: M.reshape(4,3,3)  # to get it in form you list
Out[902]: 
array([[[ 1,  2,  3],
        [ 5,  6,  7],
        [ 9, 10, 11]],

       [[ 2,  3,  4],
        [ 6,  7,  8],
        [10, 11, 12]],

       [[ 5,  6,  7],
        [ 9, 10, 11],
        [13, 14, 15]],

       [[ 6,  7,  8],
        [10, 11, 12],
        [14, 15, 16]]])

strides 的一个问题是它是先进的，并且很难向没有太多 numpy 经验的人解释。我在没有太多尝试和错误的情况下想出了这个表格，但我已经在这里闲逛了太久了。 :) )。

但是这个迭代解决方案更容易解释:

In [909]: alist=[]
In [910]: for i in range(2):
     ...:     for j in range(2):
     ...:         alist.append(A[np.ix_(range(i,i+3),range(j,j+3))])
     ...:         
In [911]: alist
Out[911]: 
[array([[ 1,  2,  3],
        [ 5,  6,  7],
        [ 9, 10, 11]]), 
 array([[ 2,  3,  4],
        [ 6,  7,  8],
        [10, 11, 12]]), 
 array([[ 5,  6,  7],
        [ 9, 10, 11],
        [13, 14, 15]]), 
 array([[ 6,  7,  8],
        [10, 11, 12],
        [14, 15, 16]])]

可以用 np.array(alist) 转换成数组.如果它更清晰，使用它没有任何问题。

关于 as_strided 要记住的一件事方法是它是一个 View ，并更改为 M可能会改变A ，并且在 M 中的一个地方发生了变化可以修改M中的几个地方.但那 reshape M可能会变成副本。所以总的来说，从 M 读取值更安全。 , 并将它们用于 sum 等计算和 mean .就地变化可能是不可预测的。

迭代解决方案会产生周围的副本。
np.ogrid 的迭代解决方案而不是 np.ix_ (否则相同的想法):

np.array([A[np.ogrid[i:i+3, j:j+3]] for i in range(2) for j in range(2)])

两个ix_和 ogrid只是构建用于索引 block 的向量对的简单方法:

In [970]: np.ogrid[0:3, 0:3]
Out[970]: 
[array([[0],
        [1],
        [2]]), array([[0, 1, 2]])]

相同的东西，但使用 slice对象:

np.array([A[slice(i,i+3), slice(j,j+3)] for i in range(2) for j in range(2)])

此列表版本将具有类似的 view行为为 as_strided解决方案(列表的元素是 View )。

对于具有非重叠 block 的 6x6，请尝试:

In [1016]: np.array([A[slice(i,i+3), slice(j,j+3)] for i in range(0,6,3) for j i
      ...: n range(0,6,3)])
Out[1016]: 
array([[[ 1,  2,  3],
        [ 7,  8,  9],
        [13, 14, 15]],

       [[ 4,  5,  6],
        [10, 11, 12],
        [16, 17, 18]],

       [[19, 20, 21],
        [25, 26, 27],
        [31, 32, 33]],

       [[22, 23, 24],
        [28, 29, 30],
        [34, 35, 36]]])

假设您想要连续的 block ，内部切片/范围不会改变，只是外部 i 的步进和 j

In [1017]: np.arange(0,6,3)
Out[1017]: array([0, 3])