在 einsum 中对多个向量的外积求和

Question

我已经阅读了einsum manual 和ajcr 的basic introduction

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

但我仍然不清楚这样做是否既可行又有效：

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

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

我的另一个可用选项是遍历数组并从我用 cython 编写的函数中调用基本函数（双外积和矩阵加法）（使用内置的 numpy 外部和 sum 函数不是选项，它太慢了）。我很可能最终也会将循环本身移动到 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))”循环中调用。

Answer 1

爱因斯坦求和只能用于问题的乘法部分（即外积）。幸运的是，求和不必逐元素执行，但您可以在归约矩阵上执行此操作。使用示例中的数组：

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'。

Answer 2

无论您可以使用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