求和 einsum 中多个向量的外积

summing outer product of multiple vectors in einsum

我已经通读了 einsum manual and ajcr's basic introduction

我在非编码环境中对爱因斯坦求和的经验为零,尽管我试图通过一些互联网研究来弥补这一点(会提供链接,但还没有两个以上的声誉)。我还尝试在 python 中使用 einsum 进行试验,看看我是否可以更好地处理事情。

但我仍然不清楚是否可以并有效地执行以下操作:

在等长 (3) 等高 (n) 的两个数组 (a 和 b) 上,逐行产生 ( row i: a on b) 的外积加上 (row 的外积i: b on a),然后将所有外积矩阵求和输出一个最终矩阵。

我知道 'i,j->ij' 会产生一个向量对另一个向量的外积——接下来的步骤让我迷路了。 ('ijk,jik->ij'绝对不是)

我的另一个可用选项是遍历数组并从我用 cython 编写的函数中调用基本函数(双外积和矩阵加法)(使用内置的 numpy 外和求和函数是不是一个选择,它太慢了)。很可能我最终也会将循环本身也移动到 cython。

所以:

  1. 我怎样才能完整地表达我上面描述的过程?

  2. 与在 cython 中做所有事情相比,它会提供真正的收益吗?或者还有其他我不知道的选择吗? (包括我使用 numpy 的效率低于我可能的可能性......)

谢谢。


用例子编辑:

A=np.zeros((3,3))
arrays_1=np.array([[1,0,0],[1,2,3],[0,1,0],[3,2,1]])
arrays_2=np.array([[1,2,3],[0,1,0],[1,0,0],[3,2,1]])
for i in range(len(arrays_1)):
  A=A+(np.outer(arrays_1[i], arrays_2[i])+np.outer(arrays_2[i],arrays_1[i]))

(但请注意,实际上我们处理的数组长度要大得多(即每个内部成员的长度仍为 3,但最多有几千个这样的成员), 这部分代码(不可避免地)被调用了很多次)

如果它有任何帮助,这里是用于求和两个外部产品的 cython:

def outer_product_sum(np.ndarray[DTYPE_t, ndim=1] a_in, np.ndarray[DTYPE_t, ndim=1] b_in):
    cdef double *a = <double *>a_in.data
    cdef double *b = <double *>b_in.data
    return np.array([
[a[0]*b[0]+a[0]*b[0], a[0]*b[1]+a[1]*b[0], a[0] * b[2]+a[2] * b[0]],
[a[1]*b[0]+a[0]*b[1], a[1]*b[1]+a[1]*b[1], a[1] * b[2]+a[2] * b[1]],
[a[2]*b[0]+a[0]*b[2], a[2]*b[1]+a[1]*b[2], a[2] * b[2]+a[2] * b[2]]])

现在,我从 'i in range(len(array))' 循环中调用,如上所示。

爱因斯坦求和只能用于问题的乘法部分(即外积)。幸运的是,求和不必按元素执行,但您可以在 reduce 矩阵上执行。使用示例中的数组:

arrays_1 = np.array([[1,0,0],[1,2,3],[0,1,0],[3,2,1]])
arrays_2 = np.array([[1,2,3],[0,1,0],[1,0,0],[3,2,1]])
A = np.einsum('ki,kj->ij', arrays_1, arrays_2) + np.einsum('ki,kj->ij', arrays_2, arrays_1)

输入数组的形状为 (4,3),求和发生在第一个索引(名为 'k')上。如果应在第二个索引上进行求和,请将下标字符串更改为 'ik,jk->ij'.

无论您可以使用 np.einsum 做什么,使用 np.dot 通常可以更快地完成。在这种情况下,A 是两个点积的总和:

arrays_1 = np.array([[1,0,0],[1,2,3],[0,1,0],[3,2,1]])
arrays_2 = np.array([[1,2,3],[0,1,0],[1,0,0],[3,2,1]])

A1 = (np.einsum('ki,kj->ij', arrays_1, arrays_2) +
      np.einsum('ki,kj->ij', arrays_2, arrays_1))

A2 = arrays_1.T.dot(arrays_2) + arrays_2.T.dot(arrays_1)

print(np.allclose(A1, A2))
# True

%timeit (np.einsum('ki,kj->ij', arrays_1, arrays_2) +
         np.einsum('ki,kj->ij', arrays_2, arrays_1))
# 100000 loops, best of 3: 7.51 µs per loop

%timeit arrays_1.T.dot(arrays_2) + arrays_2.T.dot(arrays_1)
# 100000 loops, best of 3: 4.51 µs per loop