用Python进行置换广播

Question

我知道ndarray上的transpose是为了等同于matlab的permute函数，但是我有一个具体的用例，它不能简单地工作。在matlab中，我有以下内容：

C = @bsxfun(@times, permute(A,[4,2,5,1,3]), permute(B, [1,6,2,7,3,4,5])

其中A是NxNxM形状的3D张量，B是NxNxMxPxP形状的5D张量。上面的函数用于将循环的kronecker产品矢量化。我假设Matlab为A和B都增加了2个单变量维度，这就是为什么它能够重新排列它们。我希望将这些代码移植到Python上，但我不认为它有能力添加这些额外的维度。我发现了this，它成功地增加了额外的维度，但是广播的工作方式与matlab的bsxfun不同。我已经尝试了显而易见的翻译（是的，我正在为这些ndarray和函数使用numpy）：

^{pr2}$

我得到以下错误：

^{3}$

我的第一个猜测是对a和B做一个reshape来添加那些单例维度？在

我现在得到以下错误：

mults = transpose(rho_in,[3,1,4,0,2])*transpose(proj,[0,5,1,6,2,3,4])
ValueError: operands could not be broadcast together with shapes (1,9,1,9,8) (9,1,9,1,8,40,40)

编辑：修改了我的问题，使其不再是增加单例维度，而是在python中正确地广播matlab的乘法运算。在

Answer 1

看看你的MATLAB代码，你有-

C = bsxfun(@times, permute(A,[4,2,5,1,3]), permute(B, [1,6,2,7,3,4,5])

所以，本质上-

^{pr2}$

现在，在MATLAB中，我们不得不从更高的维度中借用单例维度，这就是为什么要为6，7为{}引入{}5为A带来那么多麻烦。在

在NumPy中，我们显式地引入带有^{}/None的元素。因此，对于纽比，我们可以这样说-

B : 1 , N , 2 , N , 3 , 4 , 5 
A : N , 2 , N , 1 , 3 , N , N

，其中N表示新轴。请注意，我们需要在末尾添加新的轴，A以推进其他维度的对齐。相反，这在MATLAB中默认发生。在

使B看起来很容易，因为维度似乎是有序的，我们只需要在适当的位置添加新的轴-B[:,None,:,None,:,:,:]。在

创建这样的A并不是向前看的。忽略A中的N's，我们将得到-A : 2 , 1 , 3。所以，起点应该是置换维度，然后加入被忽略的两个新轴-A.transpose(1,0,2)[None,:,None,:,:,None,None]。在

到目前为止，我们已经-

B (new): B[:,None,:,None,:,:,:]
A (new): A.transpose(1,0,2)[None,:,None,:,:,None,None]

在NumPy中，我们可以跳过前导的新轴和尾随的非单体dim。所以，我们可以这样简化-

B (new): B[:,None,:,None]
A (new): A.transpose(1,0,2)[:,None,...,None,None]

最终的结果是这两个扩展版本之间的乘法-

C = A.transpose(1,0,2)[:,None,...,None,None]*B[:,None,:,None]

运行时测试

我相信@Andras的帖子意味着等效的np.einsum实现类似于：np.einsum('ijk,lmkno->ljmikno',A,B)。在

In [24]: A = np.random.randint(0,9,(10,10,10))
    ...: B = np.random.randint(0,9,(10,10,10,10,10))
    ...: 

In [25]: C1 = np.einsum('ijk,lmkno->ljmikno',A,B)

In [26]: C2 = A.transpose(1,0,2)[:,None,...,None,None]*B[:,None,:,None]

In [27]: np.allclose(C1,C2)
Out[27]: True

In [28]: %timeit np.einsum('ijk,lmkno->ljmikno',A,B)
10 loops, best of 3: 102 ms per loop

In [29]: %timeit A.transpose(1,0,2)[:,None,...,None,None]*B[:,None,:,None]
10 loops, best of 3: 78.4 ms per loop

In [30]: A = np.random.randint(0,9,(15,15,15))
    ...: B = np.random.randint(0,9,(15,15,15,15,15))
    ...: 

In [31]: %timeit np.einsum('ijk,lmkno->ljmikno',A,B)
1 loop, best of 3: 1.76 s per loop

In [32]: %timeit A.transpose(1,0,2)[:,None,...,None,None]*B[:,None,:,None]
1 loop, best of 3: 1.36 s per loop

Answer 2

MATLAB和numpy之间的巨大差异在于前者使用列主格式作为数组，而后者使用row major格式。其推论是隐式单变量维度的处理方式不同。在

具体地说，MATLAB显式地忽略了后面的单变量维数：rand(3,3,1,1,1,1,1)实际上是一个3x3矩阵。沿着这些思路，如果两个数组的前导维度匹配，bsxfun可以对两个数组进行操作：NxNxM是隐式的NxNxMx1x1，它与NxNxMxPxP兼容。在

另一方面，allows implicit singletons up front。您需要permute数组的尾随维度匹配，例如shape (40,40,9,1,9,1,8)与shape (1,9,1,9,8)，结果应该是(40,40,9,9,9,9,8)。在

虚拟示例：

>>> import numpy as np
>>> (np.random.rand(40,40,9,1,9,1,8)+np.random.rand(1,9,1,9,8)).shape
(40, 40, 9, 9, 9, 9, 8)

请注意，您尝试的操作可能可以使用^{}完成。我建议你仔细看看。我的意思是：根据你的问题，我收集到你想要执行的一个例子：取元素A[1:N,1:N,1:M]和B[1:N,1:N,1:M,1:P,1:P]并构造一个新的数组C[1:N,1:N,1:N,1:N,1:M,1:P,1:P]，这样

^{pr2}$

（具体的索引顺序可能会有所不同）。如果这是正确的，您确实可以使用numpy.einsum()：

>>> a = np.random.rand(3,3,2)
>>> b = np.random.rand(3,3,2,4,4)
>>> np.einsum('ijk,lmkno->limjkno',a,b).shape
(3, 3, 3, 3, 2, 4, 4)

不过，有两点需要注意。首先，上面的操作将是非常占用内存的，对于向量化的情况（通常以牺牲内存需求为代价来赢得CPU时间）来说，这是意料之中的事情。其次，在移植代码时，应该认真考虑重新安排数据模型。广播在两种语言中工作方式不同的原因与列主音/行主音的差异有着错综复杂的联系。这也意味着在MATLAB中，您应该首先使用前导索引，因为A(:,i2,i3)对应于一个连续的内存块，而{}则不是。相反，在numpy中A[i1,i2,:]是连续的，而{}不是连续的。在

这些考虑建议您设置数据的逻辑，以便向量化操作最好使用MATLAB中的前导索引和numpy中的尾随索引。您仍然可以使用numpy.einsum来执行操作本身，但是与MATLAB相比，您的维度应该以不同（可能相反）的顺序排列，至少如果我们假设两个版本的代码都使用最佳设置的话。在