4

我无法理解为什么以下内容不起作用：

我有一个可以是三维或六维的数组前置因子。我有一个具有四个维度的偶极子阵列。dipoles 的前三个维度与 prefactor的后三个维度匹配。

由于我不知道prefactor的形状，我使用省略号来说明prefactor中的三个可选维度：

numpy.einsum('...lmn,lmno->...o', prefactor, dipoles)

（在这里的例子中，prefactor.shape 是 (1, 1, 1, 160, 160, 128) 和 dipoles.shape 是 (160, 160, 128, 3)。执行时，我得到错误：

操作数 1 没有足够的维度来匹配广播，并且无法扩展，因为在开始和结束时都指定了爱因斯坦和下标

但是，当我在第二项中添加省略号时，它确实有效：

numpy.einsum('...lmn,...lmno->...o', prefactor, dipoles)

只是我不明白为什么，因为那里不需要省略号。有人知道这里发生了什么吗？

在http://comments.gmane.org/gmane.comp.python.numeric.general/53705上提出了同样的问题，但还没有令人满意的答案。

4

1 回答 1

5

这个问题有一个github问题：

https://github.com/numpy/numpy/issues/2455 einsum 中索引符号的改进（Trac #1862）

错误案例：

einsum('ij...,j->ij...',A,B)

当前的解决方法需要（空）省略号：

einsum('ij...,j...->ij...',A,B)

它看起来像einsum循环遍历字符串参数和操作数次，识别索引和广播类型（右、左、中间、无）和操作维度。有了这个，它构造了一个numpy.nditer. 正是在构建引发此错误op_axes的 nditer 时。einsum我不知道测试标准是否太严格（ibroadcast >= ndim），或者它是否需要采取额外的步骤来构造op_axes这个参数的正确性。

https://github.com/numpy/numpy/issues/2619展示了如何nditer用于复制einsum行为。以此为基础，我可以复制您的计算：

prefactor = np.random.random((1, 1, 1, 160, 160, 128))
dipoles = np.random.random((160, 160, 128, 3))
x = numpy.einsum('...lmn,...lmno->...o', prefactor, dipoles)
#numpy.einsum('...lmn,lmno->...o', prefactor, dipoles)  # not work

op_axes = [[0,1,2,3,4,5,-1], [-1,-1,-1,0,1,2,3], [0,1,2,-1,-1,-1,3]]
flags = ['reduce_ok','buffered', 'external_loop', 'delay_bufalloc', 'grow_inner']
op_flags = [['readonly']]*nops + [['allocate','readwrite']]
it = np.nditer([prefactor,dipoles,None], flags, op_flags, op_axes=op_axes)
it.operands[nops][...] = 0
it.reset()
#it.debug_print()
for (x,y,w) in it:
    w[...] += x*y
print "\nnditer usage:"
print it.operands[nops] # == x
print it.operands[nops].shape # (1, 1, 1, 3)

该op_axes线表示einsum从 推导出的内容'...lmn,...lmno->...o'。

我正在https://github.com/hpaulj/numpy-einsum上探索这个问题。

我有一个用 Python 代码einsum_py.py模拟的。np.einsum与此问题相关的部分是parse_subscripts()，特别是prepare_op_axes()。似乎只需要BROADCAST_RIGHT迭代（从末尾开始）即可正确创建op_axes，而不管椭圆在下标中的位置。它还删除了此问题的核心错误消息。

该einsum.c.src存储库上的文件具有此更改，并且可以使用当前的主发行版正确编译（只需替换文件并构建）。它可以针对test_einsum.py以及此问题的示例进行测试。

我已经为此更改提交了拉取请求。

于 2013-10-05T22:39:57.170 回答