使用排列在Python中进行广播

Question

我理解,transpose在一个ndarray意图是相当于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,但我不认为它具有添加这些额外维度的能力..我发现这成功地增加了额外的尺寸,但是广播不是同样的matlab bsxfun.我试过了明显的翻译(是的,我正在使用numpy来实现这些ndarray和功能):

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

Answer 1

MATLAB和numpy之间的巨大差异在于前者使用列主要格式作为其数组,而后者使用行主要格式.推论是隐式单例维度的处理方式不同.

具体来说,MATLAB明确忽略尾随单例维度:rand(3,3,1,1,1,1,1)实际上是一个3x3矩阵.沿着这些方向,bsxfun如果它们的前导尺寸匹配,则可以使用两个阵列进行操作:NxNxM隐式地NxNxMx1x1 兼容NxNxMxPxP.

另一方面,Numpy 允许隐含的单身人士在前面.您需要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时间).其次,您应该认真考虑在移植代码时重新安排数据模型.广播在两种语言中的工作方式不同的原因是与列主要/行主要区别错综复杂地联系在一起.这也意味着在MATLAB中你应该首先使用前导索引,因为它A(:,i2,i3)对应于一个连续的内存块,而A(i1,i2,:)不是.相反,numpy A[i1,i2,:]是连续的,而A[:,i2,i3]不是.

这些考虑建议您应该设置数据的后勤,使得向量化操作最好与MATLAB中的前导索引和numpy中的尾随索引一起使用.您仍然可以使用numpy.einsum自己执行操作,但是与MATLAB相比,您的尺寸应该处于不同的(可能是反向的)顺序,至少如果我们假设两个版本的代码都使用最佳设置.