优化从各种数组到数组的赋值 - NumPy答案

【问题标题】：Optimizing assignment into an array from various arrays - NumPy优化从各种数组到数组的赋值 - NumPy
【发布时间】：2017-06-09 13:15:46
【问题描述】：

我有四个维度为 3Nx3N 的方阵，分别称为 A、B、C 和 D。

我想将它们组合成一个矩阵。 for 循环的代码如下：

import numpy
N = 3
A = numpy.random.random((3*N, 3*N))
B = numpy.random.random((3*N, 3*N))
C = numpy.random.random((3*N, 3*N))
D = numpy.random.random((3*N, 3*N))

final = numpy.zeros((6*N, 6*N))

for i in range(N):
    for j in range(N):
        for k in range(3):
            for l in range(3):
                final[6*i + k][6*j + l] = A[3*i+k][3*j+l]
                final[6*i + k + 3][6*j + l + 3] = B[3*i+k][3*j+l]
                final[6*i + k + 3][6*j + l] = C[3*i+k][3*j+l]
                final[6*i + k][6*j + l + 3] = D[3*i+k][3*j+l]

是否可以用numpythonic的方式编写前面的for循环？

【问题讨论】：

（离题）|我在想自己，我唯一一次听到有人说“numpythonic”是SO上的其他人。但后来我查了一下，你也是 :)
@Mitch 几乎是一种补品。
是的，我正在以 Numpythonic 的方式使用不同的方式来制作 for 循环。眼睛是错误的，应该是零。
@Guillaume Jacquenot N 不是 3，是任何数字！但现在我明白了，只是做一个可执行代码！酷

标签： python numpy matrix optimization vectorization

【解决方案1】：

将array-slicing 练习成多维张量/数组的大问题！

我们将输出数组初始化为一个多维6D 数组，然后简单地对其进行切片并将被重新整形的四个数组分配为4D 数组。这样做的目的是避免任何堆叠/连接，因为在处理大型数组时，这些操作会很昂贵，而是通过对输入数组进行整形，这只是视图。

这是实现 -

out = np.zeros((N,2,3,N,2,3),dtype=A.dtype)
out[:,0,:,:,0,:] = A.reshape(N,3,N,3)
out[:,0,:,:,1,:] = D.reshape(N,3,N,3)
out[:,1,:,:,0,:] = C.reshape(N,3,N,3)
out[:,1,:,:,1,:] = B.reshape(N,3,N,3)
out.shape = (6*N,6*N)

只是为了解释一下，我们有：

            |------------------------ Axes for selecting A, B, C, D
np.zeros((N,2,3,N,2,3),dtype=A.dtype)
                  |------------------------- Axes for selecting A, B, C, D

因此，长度为(2x2) = 4 的这两个轴（第二和第五）用于在四个输入之间进行选择。

运行时测试

方法-

def original_app(A, B, C, D):
    final = np.zeros((6*N,6*N),dtype=A.dtype)
    for i in range(N):
        for j in range(N):
            for k in range(3):
                for l in range(3):
                    final[6*i + k][6*j + l] = A[3*i+k][3*j+l]
                    final[6*i + k + 3][6*j + l + 3] = B[3*i+k][3*j+l]
                    final[6*i + k + 3][6*j + l] = C[3*i+k][3*j+l]
                    final[6*i + k][6*j + l + 3] = D[3*i+k][3*j+l]
    return final

def slicing_app(A, B, C, D):
    out = np.zeros((N,2,3,N,2,3),dtype=A.dtype)
    out[:,0,:,:,0,:] = A.reshape(N,3,N,3)
    out[:,0,:,:,1,:] = D.reshape(N,3,N,3)
    out[:,1,:,:,0,:] = C.reshape(N,3,N,3)
    out[:,1,:,:,1,:] = B.reshape(N,3,N,3)
    return out.reshape(6*N,6*N)

时间和验证 -

In [147]: # Setup input arrays
     ...: N = 200
     ...: A = np.random.randint(11,99,(3*N,3*N))
     ...: B = np.random.randint(11,99,(3*N,3*N))
     ...: C = np.random.randint(11,99,(3*N,3*N))
     ...: D = np.random.randint(11,99,(3*N,3*N))
     ...: 

In [148]: np.allclose(slicing_app(A, B, C, D), original_app(A, B, C, D))
Out[148]: True

In [149]: %timeit original_app(A, B, C, D)
1 loops, best of 3: 1.63 s per loop

In [150]: %timeit slicing_app(A, B, C, D)
100 loops, best of 3: 9.26 ms per loop

【讨论】：

命令numpy.allclose(out, final) 返回True：这确保了我们想要执行的操作
@GuillaumeJacquenot 感谢您的确认！还添加了运行时测试以及验证部分。

【解决方案2】：

我将从几个一般性观察开始

对于numpy 数组，我们通常使用 [ , ] 语法而不是 [][] 最终[6*i + k][6*j + l] 最终[6*i + k, 6*j + l]

对于从其他人构建的新数组，我们经常使用诸如重塑和切片之类的东西，以便我们可以将它们作为块添加在一起，而不是使用迭代循环

举个简单的例子，取连续的差异：

 y = x[1:] - x[:-1]

关于标题，“矩阵创建”更清晰。 'load' 更多的是从文件中读取数据的意思，如np.loadtxt。

==================

所以N=1，

In [171]: A=np.arange(0,9).reshape(3,3)
In [172]: B=np.arange(10,19).reshape(3,3)
In [173]: C=np.arange(20,29).reshape(3,3)
In [174]: D=np.arange(30,39).reshape(3,3)

In [178]: final
Out[178]: 
array([[ 0,  1,  2, 30, 31, 32],
       [ 3,  4,  5, 33, 34, 35],
       [ 6,  7,  8, 36, 37, 38],
       [20, 21, 22, 10, 11, 12],
       [23, 24, 25, 13, 14, 15],
       [26, 27, 28, 16, 17, 18]])

只需调用bmat即可创建：

In [183]: np.bmat([[A,D],[C,B]]).A
Out[183]: 
array([[ 0,  1,  2, 30, 31, 32],
       [ 3,  4,  5, 33, 34, 35],
       [ 6,  7,  8, 36, 37, 38],
       [20, 21, 22, 10, 11, 12],
       [23, 24, 25, 13, 14, 15],
       [26, 27, 28, 16, 17, 18]])

bmat 混合使用hstack 和vstack。它还产生np.matrix，因此需要.A。 @Divakar's 解决方案一定会更快。

这与 N=3 不匹配。 3x3 块出现故障。但是将数组扩展为 6d（如 Divakar 所做的那样），并交换一些轴，将子块放入正确的顺序。

对于N=3：

In [57]: block=np.bmat([[A,D],[C,B]])
In [58]: b1=block.A.reshape(2,3,3,2,3,3)
In [59]: b2=b1.transpose(1,0,2,4,3,5)
In [60]: b3=b2.reshape(18,18)
In [61]: np.allclose(b3,final)
Out[61]: True

在快速时间测试 (N=3) 中，我的方法大约是 slicing_app 的一半。

出于好奇，bmat 使用字符串输入：np.bmat('A,D;C,B')。那是因为 np.matrix 在几年前就试图给人一种 MATLAB 的感觉。

【讨论】：

【解决方案3】：

你可以只连接它们

水平连接 A 和 B 横向连接 C 和 D

AB的连词和CD的连词垂直连接

示例：

AB = numpy.concatenate([A,B],1)
CD = numpy.concatenate([C,D],1)
ABCD = numpy.concatenate([AB,CD],0)

希望对你有帮助:)

【讨论】：

结果不是预期的：numpy.allclose(ABCD, final) 返回 false
您是否尝试仅打印 ABCD？它应该是所需的矩阵。您可以（应该）检查手动打印矩阵的值，看看它是否符合预期。尝试使用较小的 2x2 矩阵作为示例/测试
你在正确的轨道上。我对bmat 的回答就是这种concatenate。但是 OP 的布局更复杂，需要一些重塑和轴交换。

【解决方案4】：

使用view_as_blocks 的另一种方法：

from skimage.util import view_as_blocks

def by_blocks():
    final = numpy.empty((6*N,6*N))
    a,b,c,d,f= [view_as_blocks(X,(3,3)) for X in [A,B,C,D,final]]
    f[0::2,0::2]=a
    f[1::2,1::2]=b
    f[1::2,0::2]=c
    f[0::2,1::2]=d
    return final

您只需逐块思考，让 view_as_blocks 为您管理步幅和形状。它与其他 numpy 解决方案一样快。

【讨论】：