求和 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。
所以:
我怎样才能完整地表达我上面描述的过程?
与在 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
我已经通读了 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。
所以:
我怎样才能完整地表达我上面描述的过程?
与在 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