8

np.einsum用来乘以概率表,例如:

np.einsum('ijk,jklm->ijklm', A, B)

问题是我正在处理超过 26 个随机变量(轴),所以如果我为每个随机变量分配一个字母,我的字母就会用完。有没有另一种方法可以指定上述操作来避免这个问题,而无需使用混乱的np.sum操作np.dot

4

3 回答 3

6

简短的回答是,您可以使用 52 个字母(大写和小写)中的任何一个。这就是英语中的所有字母。任何更高级的轴名称都必须映射到这 52 个或一组等效的数字上。实际上,您将希望在任何一次einsum通话中使用这 52 个中的一小部分。


@kennytm建议使用替代输入语法。一些示例运行表明这不是解决方案。26 仍然是实际限制(尽管有可疑的错误消息)。

In [258]: np.einsum(np.ones((2,3)),[0,20],np.ones((3,4)),[20,2],[0,2])
Out[258]: 
array([[ 3.,  3.,  3.,  3.],
       [ 3.,  3.,  3.,  3.]])

In [259]: np.einsum(np.ones((2,3)),[0,27],np.ones((3,4)),[27,2],[0,2])
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-259-ea61c9e50d6a> in <module>()
----> 1 np.einsum(np.ones((2,3)),[0,27],np.ones((3,4)),[27,2],[0,2])

ValueError: invalid subscript '|' in einstein sum subscripts string, subscripts must be letters

In [260]: np.einsum(np.ones((2,3)),[0,100],np.ones((3,4)),[100,2],[0,2])
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-260-ebd9b4889388> in <module>()
----> 1 np.einsum(np.ones((2,3)),[0,100],np.ones((3,4)),[100,2],[0,2])

ValueError: subscript is not within the valid range [0, 52]

我不完全确定为什么您需要超过 52 个字母(大写和小写),但我确定您需要进行某种映射。您不想一次einsum使用超过 52 个轴编写字符串。生成的迭代器将太大(对于内存或时间)。

我正在描绘某种可以用作的映射函数:

 astr = foo(A.names, B.names)
 # foo(['i','j','k'],['j','k','l','m'])
 # foo(['a1','a2','a3'],['a2','a3','b4','b5'])
 np.einsum(astr, A, B)

https://github.com/hpaulj/numpy-einsum/blob/master/einsum_py.py

是 Python 版本的einsum. 粗略地说einsum,解析下标字符串,创建一个op_axes列表,可用于np.nditer设置所需的积和计算。使用此代码,我可以查看翻译是如何完成的:

__name__块中的一个例子:

    label_str, op_axes = parse_subscripts('ik,kj->ij',    Labels([A.ndim,B.ndim]))
    print op_axes
    # [[0, -1, 1], [-1, 1, 0], [0, 1, -1]]  fine
    # map (4,newaxis,3)(newaxis,3,2)->(4,2,newaxis)
    print sum_of_prod([A,B],op_axes)

您的示例,具有完整的诊断输出是

In [275]:  einsum_py.parse_subscripts('ijk,jklm->ijklm',einsum_py.Labels([3,4])) 
jklm
{'counts': {105: 1, 106: 2, 107: 2, 108: 1, 109: 1}, 
 'strides': [], 
 'num_labels': 5, 
 'min_label': 105, 
 'nop': 2, 
 'ndims': [3, 4], 
 'ndim_broadcast': 0, 
 'shapes': [], 
 'max_label': 109}
[('ijk', [105, 106, 107], 'NONE'), 
 ('jklm', [106, 107, 108, 109], 'NONE')]
 ('ijklm', [105, 106, 107, 108, 109], 'NONE')
iter labels: [105, 106, 107, 108, 109],'ijklm'
op_axes [[0, 1, 2, -1, -1], [-1, 0, 1, 2, 3], [0, 1, 2, 3, 4]]

Out[275]: 
(<einsum_py.Labels at 0xb4f80cac>,
 [[0, 1, 2, -1, -1], [-1, 0, 1, 2, 3], [0, 1, 2, 3, 4]])

使用'ajk,jkzZ->ajkzZ'更改标签,但结果相同op_axes


这是翻译功能的初稿。它应该适用于任何列表列表(可散列项):

def translate(ll):
    mset=set()
    for i in ll: 
        mset.update(i)
    dd={k:v for v,k in enumerate(mset)}
    x=[''.join([chr(dd[i]+97) for i in l]) for l in ll]
    #  ['cdb', 'dbea', 'cdbea']
    y=','.join(x[:-1])+'->'+x[-1]
    # 'cdb,dbea->cdbea'

In [377]: A=np.ones((3,1,2),int)
In [378]: B=np.ones((1,2,4,3),int)
In [380]: ll=[list(i) for i in ['ijk','jklm','ijklm']]
In [381]: y=translate(ll)
In [382]: y
Out[382]: 'cdb,dbea->cdbea'

In [383]: np.einsum(y,A,B).shape
Out[383]: (3, 1, 2, 4, 3)

使用set映射索引对象意味着最终的索引字符是无序的。只要您指定不应该成为问题的 RHS。我也忽略了ellipsis

==================

输入的列表版本einsum被转换为einsum_list_to_subscripts()(in numpy/core/src/multiarray/multiarraymodule.c) 中的下标字符串版本。它替换ELLIPSIS为“...”。如果( s < 0 || s > 2*26)其中s一个子列表中的数字在哪里,则会引发 [0,52] 错误消息。并转换s为字符串

        if (s < 26) {
            subscripts[subindex++] = 'A' + s;
        }
        else {
            subscripts[subindex++] = 'a' + s;

但看起来第二种情况不起作用;我收到类似 26 的错误。

ValueError: invalid subscript '{' in einstein sum subscripts string, subscripts must be letters

'a'+s是错误的,如果s>26

In [424]: ''.join([chr(ord('A')+i) for i in range(0,26)])
Out[424]: 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'

In [425]: ''.join([chr(ord('a')+i) for i in range(0,26)])
Out[425]: 'abcdefghijklmnopqrstuvwxyz'

In [435]: ''.join([chr(ord('a')+i) for i in range(26,52)])
Out[435]: '{|}~\x7f\x80\x81\x82\x83\x84\x85\x86\x87\x88\x89\x8a\x8b\x8c\x8d\x8e\x8f\x90\x91\x92\x93\x94'

'a'+s是错误的;应该是:

In [436]: ''.join([chr(ord('a')+i-26) for i in range(26,52)])
Out[436]: 'abcdefghijklmnopqrstuvwxyz'

我提交了https://github.com/numpy/numpy/issues/7741

这么长时间后这个bug的存在表明子列表格式并不常见,并且在该列表中使用大数字的频率更低。

于 2016-06-13T18:08:42.087 回答
4

您可以使用einsum(op0, sublist0, op1, sublist1, ..., [sublistout])表单代替i,j,ik->ijk,API 不限于 52 个轴*。文档中显示了这种详细形式如何对应于 ijk 形式。

OP的

np.einsum('ijk,jklm->ijklm', A, B)

将被写为

np.einsum(A, [0,1,2], B, [1,2,3,4], [0,1,2,3,4])

(* 注意:实现仍然限于 26 轴。请参阅@hpaulj 的答案他的错误报告以获取解释)


来自 numpy 示例的等价物:

>>> np.einsum('ii', a)
>>> np.einsum(a, [0,0])

>>> np.einsum('ii->i', a)
>>> np.einsum(a, [0,0], [0])

>>> np.einsum('ij,j', a, b)
>>> np.einsum(a, [0,1], b, [1])

>>> np.einsum('ji', c)
>>> np.einsum(c, [1,0])

>>> np.einsum('..., ...', 3, c)
>>> np.einsum(3, [...], c, [...])

>>> np.einsum('i,i', b, b)
>>> np.einsum(b, [0], b, [0])

>>> np.einsum('i,j', np.arange(2)+1, b)
>>> np.einsum(np.arange(2)+1, [0], b, [1])

>>> np.einsum('i...->...', a)
>>> np.einsum(a, [0, ...], [...])

>>> np.einsum('ijk,jil->kl', a, b)
>>> np.einsum(a, [0,1,2], b, [1,0,3], [2,3])
于 2016-06-13T16:57:26.117 回答
2

如果您在谈论ijk示例中的字母并且有更多可用的字母字符,那么不,您不能。

在此处此处的 einsum numpy 代码中, numpy 正在逐个检查每个字符,isalpha似乎没有办法创建超过 1 个字符的名称。

也许您可以使用大写字母,但问题的主要答案是您不能为轴命名超过 1 个字符。

于 2016-06-13T16:50:07.097 回答