在 Python 中进行排列广播
Broadcasting in Python with permutations
我知道 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 张量。上面的函数是为了矢量化循环的克罗内克产品。我假设 Matlab 为 A 和 B 添加了 2 个单一维度，这就是它能够重新排列它们的原因。我希望将此代码移植到 Python ，但我认为它没有添加这些额外维度的能力。。我发现 this 成功地添加了额外的维度，但是广播与 matlab 的 bsxfun 不一样。我尝试了明显的翻译（是的，我正在为这些 ndarray 和函数使用 numpy）：
A = A[...,None,None]
B = B[...,None,None]
C = transpose(A,[3,1,4,0,2])*transpose(B,[0,5,1,6,2,3,4])
我收到以下错误：
return transpose(axes)
ValueError: axes don't match array
我的第一个猜测是对 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 乘法的问题
MATLAB 和 numpy 之间的巨大区别在于前者对其数组使用列优先格式，而后者使用行优先格式。推论是隐式单例维度的处理方式不同。
具体来说，MATLAB 明确地忽略了尾随的单一维度：rand(3,3,1,1,1,1,1) 实际上是一个 3x3 矩阵。按照这些思路，您可以使用 bsxfun 对两个数组进行操作，如果它们的 前导 维度匹配： NxNxM 隐式 NxNxMx1x1 与 NxNxMxPxP.
Numpy，另一方面，allows implicit singletons up front。您需要 permute 您的数组，使其 尾部 维度匹配，例如形状 (40,40,9,1,9,1,8) 与形状 (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)
请注意，您尝试执行的操作可能可以使用 numpy.einsum 完成。我建议仔细研究一下。我的意思的一个例子：从你的问题中我了解到你想执行这个：获取元素 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] 这样
C[i1,i2,i3,i4,i5,i6,i7] = A[i2,i4,i5]*B[i1,i3,i5,i6,i7]
（您的特定索引顺序可能会有所不同）。如果这是正确的，你确实可以使用 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 时间）。其次，您应该认真考虑在移植代码时重新安排数据模型。广播在两种语言中的工作方式不同的原因与 column-major/row-major 差异错综复杂地联系在一起。这也意味着在 MATLAB 中，您应该首先使用 leading 索引，因为 A(:,i2,i3) 对应于连续的内存块，而 A(i1,i2,:) 则不是。相反，在 numpy 中 A[i1,i2,:] 是连续的，而 A[:,i2,i3] 不是。
这些考虑表明您应该设置数据的逻辑，以便矢量化操作最好使用 MATLAB 中的前导索引和 numpy 中的尾随索引。您仍然可以使用 numpy.einsum 来执行操作本身，但是与 MATLAB 相比，您的维度应该采用不同的（可能是相反的）顺序，至少如果我们假设两个版本的代码都使用最佳设置。
查看您的 MATLAB 代码，您有 -
C = bsxfun(@times, permute(A,[4,2,5,1,3]), permute(B, [1,6,2,7,3,4,5])
所以，本质上-
B : 1 , 6 , 2 , 7 , 3 , 4 , 5
A : 4 , 2 , 5 , 1 , 3
现在，在 MATLAB 中，我们不得不从更高维度借用单一维度，这就是为什么要为 B 引入 dims 6、7 和 dims [=20 带来的所有麻烦=] 5 对于 A.
在 NumPy 中，我们使用 np.newaxis/None 显式引入那些。因此，对于 NumPy，我们可以这样说 -
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 中，我们可以跳过前导新轴和尾随非单例 dims。所以，我们可以像这样简化 -
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 的 post 意味着等效的 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
我知道 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 张量。上面的函数是为了矢量化循环的克罗内克产品。我假设 Matlab 为 A 和 B 添加了 2 个单一维度，这就是它能够重新排列它们的原因。我希望将此代码移植到 Python ，但我认为它没有添加这些额外维度的能力。。我发现 this 成功地添加了额外的维度，但是广播与 matlab 的 bsxfun 不一样。我尝试了明显的翻译（是的，我正在为这些 ndarray 和函数使用 numpy）：
A = A[...,None,None]
B = B[...,None,None]
C = transpose(A,[3,1,4,0,2])*transpose(B,[0,5,1,6,2,3,4])
我收到以下错误：
return transpose(axes)
ValueError: axes don't match array
我的第一个猜测是对 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 乘法的问题MATLAB 和 numpy 之间的巨大区别在于前者对其数组使用列优先格式，而后者使用行优先格式。推论是隐式单例维度的处理方式不同。
具体来说，MATLAB 明确地忽略了尾随的单一维度：rand(3,3,1,1,1,1,1) 实际上是一个 3x3 矩阵。按照这些思路，您可以使用 bsxfun 对两个数组进行操作，如果它们的 前导 维度匹配： NxNxM 隐式 NxNxMx1x1 与 NxNxMxPxP.
Numpy，另一方面，allows implicit singletons up front。您需要 permute 您的数组，使其 尾部 维度匹配，例如形状 (40,40,9,1,9,1,8) 与形状 (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)
请注意，您尝试执行的操作可能可以使用 numpy.einsum 完成。我建议仔细研究一下。我的意思的一个例子：从你的问题中我了解到你想执行这个：获取元素 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] 这样
C[i1,i2,i3,i4,i5,i6,i7] = A[i2,i4,i5]*B[i1,i3,i5,i6,i7]
（您的特定索引顺序可能会有所不同）。如果这是正确的，你确实可以使用 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 时间）。其次，您应该认真考虑在移植代码时重新安排数据模型。广播在两种语言中的工作方式不同的原因与 column-major/row-major 差异错综复杂地联系在一起。这也意味着在 MATLAB 中，您应该首先使用 leading 索引，因为 A(:,i2,i3) 对应于连续的内存块，而 A(i1,i2,:) 则不是。相反，在 numpy 中 A[i1,i2,:] 是连续的，而 A[:,i2,i3] 不是。
这些考虑表明您应该设置数据的逻辑，以便矢量化操作最好使用 MATLAB 中的前导索引和 numpy 中的尾随索引。您仍然可以使用 numpy.einsum 来执行操作本身，但是与 MATLAB 相比，您的维度应该采用不同的（可能是相反的）顺序，至少如果我们假设两个版本的代码都使用最佳设置。
查看您的 MATLAB 代码，您有 -
C = bsxfun(@times, permute(A,[4,2,5,1,3]), permute(B, [1,6,2,7,3,4,5])
所以，本质上-
B : 1 , 6 , 2 , 7 , 3 , 4 , 5
A : 4 , 2 , 5 , 1 , 3
现在，在 MATLAB 中，我们不得不从更高维度借用单一维度，这就是为什么要为 B 引入 dims 6、7 和 dims [=20 带来的所有麻烦=] 5 对于 A.
在 NumPy 中，我们使用 np.newaxis/None 显式引入那些。因此，对于 NumPy，我们可以这样说 -
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 中，我们可以跳过前导新轴和尾随非单例 dims。所以，我们可以像这样简化 -
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 的 post 意味着等效的 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