numpy 点积和矩阵积

Question

我正在使用形状为 (N,)、(N,3) 和 (N,3,3) 的 numpy 数组，它们表示 3D 中的标量、向量和矩阵序列 space。我已经实现了逐点点积、矩阵乘法和 matrix/vector 乘法，如下所示：

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)

如您所见，我使用 einsum 是因为没有更好的解决方案。令我惊讶的是，我无法使用 np.dot... 这似乎不适合这种需要。有没有更numpythonic的方式来实现这些功能？

特别是，如果这些函数也可以通过广播第一个缺失的轴来处理形状 (3,) 和 (3,3)，那就太好了。我想我需要省略号，但我不太明白如何实现结果。

Answer 1

无法将这些操作重塑为一般的 BLAS 调用，并且对于这种大小的数组，循环 BLAS 调用会非常慢。因此，einsum 可能是此类操作的最佳选择。

您的函数可以用省略号概括如下：

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)

Answer 2

作为工作笔记，这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])

或者用 Ophion 的概括

np.einsum(A,[Ellipsis,1,2], B, ...)

根据输入数组的维度生成 [0,1,..] 列表应该不难。

通过专注于概括 einsum 表达式，我忽略了一个事实，即您要重现的是 N 小点积。

np.array([np.dot(i,j) for i,j in zip(a,b)])

值得记住的是 np.dot 使用快速编译代码，并专注于数组较大的计算。你的问题是计算许多小点积之一。

并且没有定义轴的额外参数，np.dot 仅执行 2 种可能的组合，可以表示为：

np.einsum('i,i', v1, v2)
np.einsum('...ij,...jk->...ik', m1, m2)

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)

它可以执行 dot 对多个轴求和。但是在转置和整形完成后，计算就变成了二维数组的标准dot。

这可能已被标记为重复问题。一段时间以来，人们一直在询问关于做多个点积的问题。

Matrix vector multiplication along array axes 建议使用 numpy.core.umath_tests.matrix_multiply

等于：

matrix_multiply(matrices, vectors[..., None])
np.einsum('ijk,ik->ij', matrices, vectors)

matrix_multiply 注释的 C 文档：

* This implements the function
* out[k, m, p] = sum_n { in1[k, m, n] * in2[k, n, p] }.

来自同一目录的

inner1d 对 (N,n) 向量

做同样的事情

inner1d(vector, vector)  
np.einsum('ij,ij->i', vector, vector)
# out[n] = sum_i { in1[n, i] * in2[n, i] }

两者都是 UFunc，并且可以处理最右侧维度的广播。在numpy/core/test/test_ufunc.py中，这些函数用于执行UFunc机制。

matrix_multiply(np.ones((4,5,6,2,3)),np.ones((3,2)))

补充说这种计算可以用*求和，eg

(w*v).sum(-1)
(M*v[...,None]).sum(-1)
(A*B.swapaxes(...)).sum(-1)

在进一步测试中，我认为 inner1d 和 matrix_multiply 匹配您的 dot 和 matrix-matrix 产品案例，如果您添加 matrix-vector 案例[...,None]。看起来它们比 einsum 版本快 2 倍（在我的机器和测试阵列上）。

https://github.com/numpy/numpy/blob/master/doc/neps/return-of-revenge-of-matmul-pep.rst 是 numpy 上对 @ 中缀运算符的讨论。我认为 numpy 开发人员对这个 PEP 的热情不如 Python 开发人员。

numpy 点积和矩阵积

numpy dot product and matrix product

python

numpy

matrix