通过平均多个单元简化矩阵

Question

我有一个大的二维numpy矩阵，需要缩小（例如：从100x100转换为10x10）。你知道吗

我的目标基本上是：将nxn矩阵分解成更小的mxm矩阵，平均这些mxm切片中的细胞，然后用这些mxm切片构造一个新的（更小的）矩阵。你知道吗

我在考虑使用matrix[a::b, c::d]之类的方法来提取较小的矩阵，然后对这些值求平均值，但这似乎过于复杂了。有没有更好的方法来实现这一点？你知道吗

Answer 1

我来晚了，但我推荐scipy.ndimage.zoom()作为一个现成的解决方案。它使用从0到5的任意顺序的样条插值来减小大小（或增大大小）。从你的问题来看，订单0对你来说已经足够了。你知道吗

from scipy import ndimage as ndi
import numpy as np

M=np.arange(1000000).reshape(1000,1000)

shrinkby=10

Mfilt = ndi.filters.uniform_filter(input=M, size=shrinkby)
Msmall = ndi.interpolation.zoom(input=Mfilt, zoom=1./shrinkby, order=0)

这就是你所需要的。指定缩放比指定所需的输出大小稍微不方便，但至少对于order=0来说，这种方法非常快。你知道吗

输出大小是每个维度输入的10%，即

print M.shape, Msmall.shape

给出(1000, 1000) (100, 100)和您可以从中获得的速度

%timeit Mfilt = ndi.filters.uniform_filter(input=M, size=shrinkby)
%timeit Msmall = ndi.interpolation.zoom(input=Mfilt, zoom=1./shrinkby, order=0)

在我的机器上给10 loops, best of 3: 20.5 ms per loop的uniform_filter调用和1000 loops, best of 3: 1.67 ms per loop的zoom调用。你知道吗

Answer 2

如果没有ski-learn，您可以简单地重塑，并采取适当的平均值。你知道吗

M=np.arange(10000).reshape(100,100)
M1=M.reshape(10,10,10,10)
M2=M1.mean(axis=(1,3))

快来看看我的斧子对不对

In [127]: M2[0,0]
Out[127]: 454.5

In [128]: M[:10,:10].mean()
Out[128]: 454.5

In [131]: M[-10:,-10:].mean()
Out[131]: 9544.5

In [132]: M2[-1,-1]
Out[132]: 9544.5

添加.transpose([0,2,1,3])将2个平均维度放在末尾，就像view_as_blocks一样。你知道吗

对于这种(100,100)情况，重塑方法比as_strided方法快2倍，但两者都相当快。你知道吗

然而，直接的解决方案并不比重塑慢多少。你知道吗

as_strided(M,shape=(10,10,10,10),strides=(8000,80,800,8)).mean((2,3))
as_strided(M,shape=(10,10,10,10),strides=(8000,800,80,8)).mean((1,3))

Answer 3

可以使用^{}函数（在scikit图像中）将数组拆分为块。你知道吗

对于2D数组，这将返回一个4D数组，其块按行排序：

>>> import skimage.util as ski
>>> import numpy as np
>>> a = np.arange(16).reshape(4,4) # 4x4 array
>>> ski.view_as_blocks(a, (2,2))
array([[[[ 0,  1],
         [ 4,  5]],

        [[ 2,  3],
         [ 6,  7]]],


       [[[ 8,  9],
         [12, 13]],

        [[10, 11],
         [14, 15]]]])

沿最后两个轴取平均值返回一个2D数组，每个块中都有平均值：

>>> ski.view_as_blocks(a, (2,2)).mean(axis=(2,3))
array([[  2.5,   4.5],
       [ 10.5,  12.5]])

注意：view_as_blocks通过修改步长返回数组视图（它也适用于二维以上的数组）。它是使用as_strided在NumPy中实现的，因此如果您没有访问scikit映像库的权限，您可以copy the code from here。你知道吗