3

我知道transposeon anndarray旨在等同于 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 ，但我认为它没有添加这些额外维度的能力。. 我发现这成功地添加了额外的维度，但是广播与 matlab 的bsxfun. 我已经尝试了明显的翻译（是的，我正在为这些ndarray's 和函数使用 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

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

我现在收到以下错误：

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 乘法。

4

2 回答 2

6

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

具体来说，MATLAB 明确忽略了尾随单例维度：rand(3,3,1,1,1,1,1)实际上是一个3x3矩阵。沿着这些思路，如果它们的前导维度匹配，您可以使用bsxfun对两个数组进行操作：隐含地 与.NxNxMNxNxMx1x1NxNxMxPxP

另一方面，Numpy允许在前面使用隐式单例。您需要以它们的尾随permute尺寸匹配的方式对您的数组进行匹配，例如 shape与 shape ，结果应该是 shape 。(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 中您应该首先使用前导索引，因为A(:,i2,i3)对应于连续的内存块，而A(i1,i2,:)没有。相反，在 numpyA[i1,i2,:]中是连续的，而A[:,i2,i3]不是。

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

于 2016-09-09T23:03:11.793 回答
2

查看您的 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 中，我们不得不从更高的维度借用单例维度，这就是引入 dims 6、7forB和 dims 4 5for的所有麻烦的原因A。

np.newaxis在 NumPy 中，我们使用/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的看起来并不简单。忽略N's，A我们会有 - 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 中，我们可以跳过前导的新轴和尾随的非单例暗淡。所以，我们可以像这样简化——

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
于 2016-09-10T07:40:14.713 回答