计算任意维度数组的外积

Question

我有两个数组A和B，想在它们的最后一个维度上进行外积运算，例如：当A和B是二维时，result[:,i,j]=A[:,i]*B[:,j]。如果我不知道它们是二维还是三维，该怎么做呢？

在我的具体问题中，A和B是从一个更大的三维数组Z中切片得到的。有时可能会使用整数索引A=Z[:,1,:]，B=Z[:,2,:]，而其他时候可能会使用切片A=Z[:,1:3,:]，B=Z[:,4:6,:]。由于scipy“压缩”单例维度，我无法知道我的输入将有哪些维度。

我试图定义的数组外积应满足：

array_outer_product( Y[a,b,:], Z[i,j,:] ) == scipy.outer( Y[a,b,:], Z[i,j,:] )
array_outer_product( Y[a:a+N,b,:], Z[i:i+N,j,:])[n,:,:] == scipy.outer( Y[a+n,b,:], Z[i+n,j,:] ) 
array_outer_product( Y[a:a+N,b:b+M,:], Z[i:i+N, j:j+M,:] )[n,m,:,:]==scipy.outer( Y[a+n,b+m,:] , Z[i+n,j+m,:] )

对于任意的三维数组Y,Z和整数a,b,...i,j,k...n,N,...

我所处理的问题涉及到一个二维空间网格，在每个网格点上有一个矢量值函数。我想能够计算这些向量在第一和第二轴上定义的切片区域内的协方差矩阵（外积）。

- Dave

---

你能发布一些输入和输出的示例吗？ - user545424

你看过numpy.multiply.outer吗？我认为这就是它的作用（但我没有仔细阅读这篇文章）。 - David

Answer 1

您可以尝试使用einsum： http://docs.scipy.org/doc/numpy/reference/generated/numpy.einsum.html。 它是一个很好的工具，可用于高效地操作多维数组。

- Nicolas Barbey

---

我不知道这个函数。看起来非常有用。谢谢。 - Michelle Lynn Gill

如果我使用的是numpy >1.6.0（但我没有），我会尝试这个。 - Dave

Answer 2

在发现numpy/scipy数组中省略号的用法后，我最终将其实现为一个递归函数：

def array_outer_product(A, B, result=None):
    ''' Compute the outer-product in the final two dimensions of the given arrays.
    If the result array is provided, the results are written into it.
    '''
    assert(A.shape[:-1] == B.shape[:-1])
    if result is None:
        result=scipy.zeros(A.shape+B.shape[-1:], dtype=A.dtype)
    if A.ndim==1:
        result[:,:]=scipy.outer(A, B)
    else:
        for idx in xrange(A.shape[0]):
            array_outer_product(A[idx,...], B[idx,...], result[idx,...])
    return result

- Dave

Answer 3

假设我理解你的意思正确，几周前在我的研究中遇到了类似的问题。我意识到Kronecker积只是一个保留维度的外积。因此，你可以像这样做：

import numpy as np

# Generate some data
a = np.random.random((3,2,4))
b = np.random.random((2,5))

# Now compute the Kronecker delta function
c = np.kron(a,b)

# Check the shape
np.prod(c.shape) == np.prod(a.shape)*np.prod(b.shape)

我不确定你最终想要什么形状，但你可以使用数组切片与np.rollaxis、np.reshape、np.ravel等函数结合使用，将元素按照你的意愿进行重排。这样做的缺点是会进行一些额外的计算，具体是否影响取决于你的限制条件。

- Michelle Lynn Gill