用爱因斯坦求和来写求和的正确方法对我来说是个谜,所以我想在我的代码中尝试一下。我在一些案例中取得了成功,但主要是通过反复试验。
现在有一个案例我想不通。首先,一个基本问题。对于分别为Nx1
和1xN
的两个矩阵A
和B
,AB
是 NxN
但 BA
是 1x1
。当我想用 np.einsum
计算 NxN
情况时,我可以这样做:
import numpy as np
a = np.asarray([[1,2]])
b = np.asarray([[2,3]])
print np.einsum('ij,ji->ij', a, b)
最后的数组是2x2
。然而
a = np.asarray([[1,2]])
b = np.asarray([[2,3]])
print np.einsum('ij,ij->ij', a, b)
返回一个 1x2
数组。我不太明白为什么这不能给出正确的结果。
例如,对于上述情况,numpy
的指南说箭头可用于强制求和或阻止求和发生。但这对我来说似乎很模糊;在上述情况下,我不明白 numpy 如何根据索引的顺序(显然会发生变化)决定输出数组的最终大小。
形式上我知道以下内容:当箭头右侧没有任何内容时,可以在数学上将求和写为 $\sum\limits_{i=0}^{N}\sum\limits_{j=0 }^{M} A_{ij}B_{ij}$
对于 np.einsum('ij,ij',A,B)
,但是当有箭头时,我不知道如何用正式的数学表达式来解释它。
最佳答案
In [22]: a
Out[22]: array([[1, 2]])
In [23]: b
Out[23]: array([[2, 3]])
In [24]: np.einsum('ij,ij->ij',a,b)
Out[24]: array([[2, 6]])
In [29]: a*b
Out[29]: array([[2, 6]])
这里所有部分(包括输出)中索引的重复被解释为逐个元素的乘法。什么都没有总结。 a[i,j]*b[i,j] = c[i,j]
对于所有 i,j
.
In [25]: np.einsum('ij,ji->ij',a,b)
Out[25]:
array([[2, 4],
[3, 6]])
In [28]: np.dot(a.T,b).T
Out[28]:
array([[2, 4],
[3, 6]])
In [38]: np.outer(a,b)
Out[38]:
array([[2, 3],
[4, 6]])
同样没有求和,因为相同的索引出现在左侧和右侧。 a[i,j]*b[j,i] = c[i,j]
,换句话说:
[[1*2, 2*2],
[1*3, 2*3]]
实际上是一个外部产品。看看如何a
针对 b.T
广播可能有帮助:
In [69]: np.broadcast_arrays(a,b.T)
Out[69]:
[array([[1, 2],
[1, 2]]),
array([[2, 2],
[3, 3]])]
在语句的左侧,重复的索引表示相乘的维度。左右匹配决定是否相加。
np.einsum('ij,ji->j',a,b) # array([ 5, 10]) sum on i only
np.einsum('ij,ji->i',a,b) # array([ 5, 10]) sum on j only
np.einsum('ij,ji',a,b) # 15 sum on i and j
不久前,我编写了一个等同于 einsum
的纯 Python ,主要关注它如何解析字符串。目标是创建一个 nditer
它用它进行产品总和计算。但这不是一个简单的脚本,即使在 Python 中也是如此:
https://github.com/hpaulj/numpy-einsum/blob/master/einsum_py.py
显示这些求和规则的更简单序列:
In [53]: c=np.array([[1,2],[3,4]])
In [55]: np.einsum('ij',c)
Out[55]:
array([[1, 2],
[3, 4]])
In [56]: np.einsum('ij->i',c)
Out[56]: array([3, 7])
In [57]: np.einsum('ij->j',c)
Out[57]: array([4, 6])
In [58]: np.einsum('ij->',c)
Out[58]: 10
使用没有 1
的数组维度消除了广播并发症:
In [71]: b2=np.arange(1,7).reshape(2,3)
In [72]: np.einsum('ij,ji',a2,b2)
...
ValueError: operands could not be broadcast together with remapped shapes [original->remapped]: (2,3)->(2,3) (2,3)->(3,2)
或者我应该说,它暴露了尝试的广播。
省略号为 einsum 解释增加了一定程度的复杂性。我在解决 ...
的使用中的错误时开发了上述 github 代码.但是我并没有花太多精力来完善文档。
Ellipsis broadcasting in numpy.einsum
当您想要一个可以处理各种大小的数组的表达式时,省略号最有用。如果您的阵列始终是二维的,它不会做任何额外的事情。
例如,考虑 dot
的概括, 乘以 A
的最后一个维度倒数第二个 B
.使用省略号,我们可以编写一个可以处理 2d、3D 和更大数组混合的表达式:
np.einsum('...ij,...jk',np.ones((2,3)),np.ones((3,4))) # (2,4)
np.einsum('...ij,...jk',np.ones((5,2,3)),np.ones((3,4))) # (5,2,4)
np.einsum('...ij,...jk',np.ones((5,2,3)),np.ones((5,3,4))) # (5,2,4)
np.einsum('...ij,...jk',np.ones((5,2,3)),np.ones((7,5,3,4))) # (7,5,2,4)
np.einsum('...ij,...jk->...ik',np.ones((5,2,3)),np.ones((7,5,3,4)) # (7, 5, 2, 4)
最后一个表达式使用默认的右侧索引 ...ik
, 省略号加上非求和指标。
你原来的例子可以写成
np.einsum('...j,j...->...j',a,b)
实际上它填充了 i
(或更多维度)以匹配数组的维度。
如果 a
也可以工作或 b
是 1 天:
np.einsum('...j,j...->...j',a,b[0,:])
np.dot
泛化到更大维度的方式不同
dot(a, b)[i,j,k,m] = sum(a[i,j,:] * b[k,:,m])
用einsum表示为:
np.einsum('ijo,kom->ijkm',np.ones((2,3,4)),np.ones((3,4,2)))
可以概括为
np.einsum('...o,kom->...km',np.ones((4,)),np.ones((3,4,2)))
或
np.einsum('ijo,...om->ij...m',np.ones((2,3,4)),np.ones((3,4,2)))
但我不认为我可以在 einsum
中完全复制它.也就是说,我无法告诉它填写 A
的索引。 , 其次是不同的 B
.
关于python - `numpy.einsum` 是如何工作的?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/26064594/