python - 使用排列在 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 张量。上面的函数是为了矢量化循环的克罗内克产品。我假设 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 乘法的问题。

Answer 1

MATLAB 和 numpy 之间的巨大区别在于前者对其数组使用列优先格式，而后者使用行优先格式。推论是隐式单例维度的处理方式不同。

具体来说，MATLAB 明确忽略了尾随的单一维度:rand(3,3,1,1,1,1,1) 实际上是一个 3x3 矩阵。按照这些思路，您可以使用 bsxfun 对两个数组进行操作，前提是它们的leading 维度匹配:NxNxM 隐含地是 NxNxMx1x1 与 NxNxMxPxP 兼容。

另一方面，Numpy allows implicit singletons up front .您需要以尾随维度匹配的方式置换您的数组，例如形状(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 时间)。其次，您应该认真考虑在移植代码时重新安排数据模型。广播在两种语言中的工作方式不同的原因与列优先/行优先的差异错综复杂地联系在一起。这也意味着在 MATLAB 中你应该首先使用 leading 索引，因为 A(:,i2,i3) 对应于一个连续的内存块，而 A (i1,i2,:) 没有。相反，在 numpy 中 A[i1,i2,:] 是连续的，而 A[:,i2,i3] 不是。

这些考虑表明您应该设置数据的逻辑，以便矢量化操作最好与 MATLAB 中的前导索引和 numpy 中的尾随索引一起使用。您仍然可以使用 numpy.einsum 来执行操作本身，但是与 MATLAB 相比，您的尺寸应该采用不同的(可能是相反的)顺序，至少如果我们假设代码的两个版本都使用最佳设置。