I'm working with numpy arrays of shape (N,), (N,3) and (N,3,3) which represent sequences of scalars, vectors and matrices in 3D space. I have implemented pointwise dot product, matrix multiplication, and matrix/vector multiplication as follows:
我使用的是形状(N，)， (N,3)和(N,3,3)的numpy数组，它们表示三维空间中的标量、向量和矩阵的序列。我已经实现了点积，矩阵乘法，矩阵/向量乘法如下:
def dot_product(v, w):
return np.einsum('ij, ij -> i', v, w)
def matrix_vector_product(M, v):
return np.einsum('ijk, ik -> ij', M, v)
def matrix_matrix_product(A, B):
return np.einsum('ijk, ikl -> ijl', A, B)
As you can see I use einsum for lack of a better solution. To my surprise I was not able to use np.dot... which seems not suitable for this need. Is there a more numpythonic way to implement these function?
如您所见，我使用einsum是因为没有更好的解决方案。令我惊讶的是我不能使用np.dot…这似乎不适合这种需要。是否有更好的方法来实现这些函数?
In particular it would be nice if the functions could work also on the shapes (3,) and (3,3) by broadcasting the first missing axis. I think I need ellipsis, but I don't quite understand how to achieve the result.
特别是，如果函数可以通过广播第一个缺失的轴在形状(3，)和(3,3)上工作，那就更好了。我想我需要省略，但我不太明白如何才能达到目的。
2 个解决方案
#1
4
These operations cannot be reshaped into general BLAS calls and looping BLAS calls would be quite slow for arrays of this size. As such, einsum is likely optimal for this kind of operation.
这些操作不能被重构为一般的BLAS调用，对于这种大小的数组来说，循环的BLAS调用将非常缓慢。因此，einsum可能是这种操作的最佳选择。
Your functions can be generalized with ellipses as follows:
你的函数可以概括为省略号如下:
def dot_product(v, w):
return np.einsum('...j,...j->...', v, w)
def matrix_vector_product(M, v):
return np.einsum('...jk,...k->...j', M, v)
def matrix_matrix_product(A, B):
return np.einsum('...jk,...kl->...jl', A, B)
#2
2
Just as working notes, these 3 calculations can also be written as:
正如工作笔记一样，这3项计算也可以写成:
np.einsum(A,[0,1,2],B,[0,2,3],[0,1,3])
np.einsum(M,[0,1,2],v,[0,2],[0,1])
np.einsum(w,[0,1],v,[0,1],[0])
Or with Ophion's generalization
或与Ophion泛化
np.einsum(A,[Ellipsis,1,2], B, ...)
It shouldn't be hard to generate the [0,1,..] lists based on the dimensions of the inputs arrays.
生成[0,1]应该不难。基于输入数组的维度的列表。
By focusing on generalizing the einsum expressions, I missed the fact that what you are trying to reproduce is N small dot products.
通过对einsum表达式进行归纳，我忽略了这样一个事实，即你要复制的是N个小点积。
np.array([np.dot(i,j) for i,j in zip(a,b)])
It's worth keeping mind that np.dot uses fast compiled code, and focuses on calculations where the arrays are large. Where as your problem is one of calculating many small dot products.
值得注意的是np。dot使用快速编译的代码，并专注于数组较大的计算。当你的问题是计算许多小点积的时候。
And without extra arguments that define axes, np.dot performs just 2 of the possible combinations, ones which can be expressed as:
没有定义坐标轴的额外参数，np。dot只执行两种可能的组合，这些组合可以表示为:
np.einsum('i,i', v1, v2)
np.einsum('...ij,...jk->...ik', m1, m2)
An operator version of dot would face the same limitation - no extra parameters to specify how the axes are to be combined.
一个操作符版本的dot将面临同样的限制——没有额外的参数来指定如何组合这些轴。
It may also be instructive to note what tensordot does to generalize dot:
注意tensordot是如何概括dot的:
def tensordot(a, b, axes=2):
....
newshape_a = (-1, N2)
...
newshape_b = (N2, -1)
....
at = a.transpose(newaxes_a).reshape(newshape_a)
bt = b.transpose(newaxes_b).reshape(newshape_b)
res = dot(at, bt)
return res.reshape(olda + oldb)
It can perform a dot with summation over several axes. But after the transposing and reshaping is done, the calculation becomes the standard dot with 2d arrays.
它可以对几个轴做一个点求和。但在置换和重组完成后，计算变成了二维数组的标准点。
This could have been flagged as a duplicate issue. People have asking about doing multiple dot products for some time.
这可能被标记为一个重复的问题。一段时间以来，人们一直在问做多个点积。
Matrix vector multiplication along array axes suggests using numpy.core.umath_tests.matrix_multiply
矩阵向量沿阵列轴相乘建议使用numpi .core.umath_tests.matrix_multiply
https://*.com/a/24174347/901925 equates:
https://*.com/a/24174347/901925相当于:
matrix_multiply(matrices, vectors[..., None])
np.einsum('ijk,ik->ij', matrices, vectors)
The C documentation for matrix_multiply notes:
矩阵乘法的C文档说明:
* This implements the function
* out[k, m, p] = sum_n { in1[k, m, n] * in2[k, n, p] }.
inner1d from the same directory does the same same for (N,n) vectors
来自同一目录的inner1d对(N, N)向量执行相同的操作。
inner1d(vector, vector)
np.einsum('ij,ij->i', vector, vector)
# out[n] = sum_i { in1[n, i] * in2[n, i] }
Both are UFunc, and can handle broadcasting on the right most dimensions. In numpy/core/test/test_ufunc.py these functions are used to exercise the UFunc mechanism.
两者都是UFunc，并且可以处理在最正确的维度上的广播。核心/测试/ test_ufunc numpy /。py这些函数用于执行UFunc机制。
matrix_multiply(np.ones((4,5,6,2,3)),np.ones((3,2)))
https://*.com/a/16704079/901925 adds that this kind of calculation can be done with * and sum, eg
https://*.com/a/16704079/901925补充说这种计算可以用*和sum来完成。
(w*v).sum(-1)
(M*v[...,None]).sum(-1)
(A*B.swapaxes(...)).sum(-1)
On further testing, I think inner1d and matrix_multiply match your dot and matrix-matrix product cases, and the matrix-vector case if you add the [...,None]. Looks like they are 2x faster than the einsum versions (on my machine and test arrays).
在进一步的测试中，我认为inner1d和matrix_multiply与你的点和矩阵-矩阵乘积的情况相匹配，而矩阵-向量的情况如果你加上[…，None]。看起来它们比einsum版本(在我的机器和测试数组上)快了2x。
https://github.com/numpy/numpy/blob/master/doc/neps/return-of-revenge-of-matmul-pep.rst is the discussion of the @ infix operator on numpy. I think the numpy developers are less enthused about this PEP than the Python ones.
https://github.com/numpy/numpy/blob/master/doc/neps/return of matmu -pep.rst是对numpy上的@ infix操作符的讨论。我认为numpy开发人员对这个PEP的热情不如Python开发人员。
#1
4
These operations cannot be reshaped into general BLAS calls and looping BLAS calls would be quite slow for arrays of this size. As such, einsum is likely optimal for this kind of operation.
这些操作不能被重构为一般的BLAS调用，对于这种大小的数组来说，循环的BLAS调用将非常缓慢。因此，einsum可能是这种操作的最佳选择。
Your functions can be generalized with ellipses as follows:
你的函数可以概括为省略号如下:
def dot_product(v, w):
return np.einsum('...j,...j->...', v, w)
def matrix_vector_product(M, v):
return np.einsum('...jk,...k->...j', M, v)
def matrix_matrix_product(A, B):
return np.einsum('...jk,...kl->...jl', A, B)
#2
2
Just as working notes, these 3 calculations can also be written as:
正如工作笔记一样，这3项计算也可以写成:
np.einsum(A,[0,1,2],B,[0,2,3],[0,1,3])
np.einsum(M,[0,1,2],v,[0,2],[0,1])
np.einsum(w,[0,1],v,[0,1],[0])
Or with Ophion's generalization
或与Ophion泛化
np.einsum(A,[Ellipsis,1,2], B, ...)
It shouldn't be hard to generate the [0,1,..] lists based on the dimensions of the inputs arrays.
生成[0,1]应该不难。基于输入数组的维度的列表。
By focusing on generalizing the einsum expressions, I missed the fact that what you are trying to reproduce is N small dot products.
通过对einsum表达式进行归纳，我忽略了这样一个事实，即你要复制的是N个小点积。
np.array([np.dot(i,j) for i,j in zip(a,b)])
It's worth keeping mind that np.dot uses fast compiled code, and focuses on calculations where the arrays are large. Where as your problem is one of calculating many small dot products.
值得注意的是np。dot使用快速编译的代码，并专注于数组较大的计算。当你的问题是计算许多小点积的时候。
And without extra arguments that define axes, np.dot performs just 2 of the possible combinations, ones which can be expressed as:
没有定义坐标轴的额外参数，np。dot只执行两种可能的组合，这些组合可以表示为:
np.einsum('i,i', v1, v2)
np.einsum('...ij,...jk->...ik', m1, m2)
An operator version of dot would face the same limitation - no extra parameters to specify how the axes are to be combined.
一个操作符版本的dot将面临同样的限制——没有额外的参数来指定如何组合这些轴。
It may also be instructive to note what tensordot does to generalize dot:
注意tensordot是如何概括dot的:
def tensordot(a, b, axes=2):
....
newshape_a = (-1, N2)
...
newshape_b = (N2, -1)
....
at = a.transpose(newaxes_a).reshape(newshape_a)
bt = b.transpose(newaxes_b).reshape(newshape_b)
res = dot(at, bt)
return res.reshape(olda + oldb)
It can perform a dot with summation over several axes. But after the transposing and reshaping is done, the calculation becomes the standard dot with 2d arrays.
它可以对几个轴做一个点求和。但在置换和重组完成后，计算变成了二维数组的标准点。
This could have been flagged as a duplicate issue. People have asking about doing multiple dot products for some time.
这可能被标记为一个重复的问题。一段时间以来，人们一直在问做多个点积。
Matrix vector multiplication along array axes suggests using numpy.core.umath_tests.matrix_multiply
矩阵向量沿阵列轴相乘建议使用numpi .core.umath_tests.matrix_multiply
https://*.com/a/24174347/901925 equates:
https://*.com/a/24174347/901925相当于:
matrix_multiply(matrices, vectors[..., None])
np.einsum('ijk,ik->ij', matrices, vectors)
The C documentation for matrix_multiply notes:
矩阵乘法的C文档说明:
* This implements the function
* out[k, m, p] = sum_n { in1[k, m, n] * in2[k, n, p] }.
inner1d from the same directory does the same same for (N,n) vectors
来自同一目录的inner1d对(N, N)向量执行相同的操作。
inner1d(vector, vector)
np.einsum('ij,ij->i', vector, vector)
# out[n] = sum_i { in1[n, i] * in2[n, i] }
Both are UFunc, and can handle broadcasting on the right most dimensions. In numpy/core/test/test_ufunc.py these functions are used to exercise the UFunc mechanism.
两者都是UFunc，并且可以处理在最正确的维度上的广播。核心/测试/ test_ufunc numpy /。py这些函数用于执行UFunc机制。
matrix_multiply(np.ones((4,5,6,2,3)),np.ones((3,2)))
https://*.com/a/16704079/901925 adds that this kind of calculation can be done with * and sum, eg
https://*.com/a/16704079/901925补充说这种计算可以用*和sum来完成。
(w*v).sum(-1)
(M*v[...,None]).sum(-1)
(A*B.swapaxes(...)).sum(-1)
On further testing, I think inner1d and matrix_multiply match your dot and matrix-matrix product cases, and the matrix-vector case if you add the [...,None]. Looks like they are 2x faster than the einsum versions (on my machine and test arrays).
在进一步的测试中，我认为inner1d和matrix_multiply与你的点和矩阵-矩阵乘积的情况相匹配，而矩阵-向量的情况如果你加上[…，None]。看起来它们比einsum版本(在我的机器和测试数组上)快了2x。
https://github.com/numpy/numpy/blob/master/doc/neps/return-of-revenge-of-matmul-pep.rst is the discussion of the @ infix operator on numpy. I think the numpy developers are less enthused about this PEP than the Python ones.
https://github.com/numpy/numpy/blob/master/doc/neps/return of matmu -pep.rst是对numpy上的@ infix操作符的讨论。我认为numpy开发人员对这个PEP的热情不如Python开发人员。