NumPy点积和矩阵乘积

Question

我正在处理形状为(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)

如您所见，我使用einsum作为最佳解决方案。 令我惊讶的是，我无法使用np.dot……似乎不适合这个需要。 是否有更多numpythonic的方法来实现这些函数？

特别是如果函数可以通过广播第一个缺失轴来处理形状（3，）和（3,3），那将是很好的。 我认为我需要ellipsis，但我不太理解如何实现结果。

- Emanuele Paolini

---

根据两个输入的维度生成'ij,ij...'字符串不应该很难。还有一种子列表方法可以指定这些索引。 - hpaulj

如果N=3，规则是什么？即两个(3,3)的数组。是'ij,ij->i'还是'jk,ik->ij'还是'jk,kl->jl'? - hpaulj

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)

- Daniel

---

这些适用于 v=np.ones((3,)) 和 A=np.ones((3,3))。只是没有明显的方法来组合这三种情况。 - hpaulj

对我来说，看起来很奇怪的是，在这个问题 http://stackoverflow.com/questions/28230296/numpy-row-order-in-the-representation-of-vector-funcions 中，有人指向了 http://legacy.python.org/dev/peps/pep-0465/ 该提案建议在Python中添加一个新的运算符“@”来完全执行我需要的操作（如果我理解正确的话）。因此，似乎NumPy社区正在提议一个新的Python运算符来执行目前没有任何专用函数执行的操作... - Emanuele Paolini

Python社区可以添加@运算符，但是给ndarray赋予意义的将是numpy开发人员。它是否包括您想要的广播类型是一个未决问题。 - hpaulj

@EmanuelePaolini 将矩阵乘法推广到更高的维度存在问题，因为有两种方法可以实现。Numpy将其视为组合内积（文档），这与您在此处提出的方法不同。 - Daniel

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 只执行可能的组合中的两个，可以表示为：

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)

它可以在多个轴上执行带求和的“点乘”运算。但在转置和重塑完成后，计算就变成了标准的二维数组“点乘”。

---

这可能被标记为重复问题。人们一直在询问如何进行多个点积。 沿数组轴进行矩阵向量乘法建议使用numpy.core.umath_tests.matrix_multiply https://stackoverflow.com/a/24174347/901925中等于：

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

< p > 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)))

这句话的意思是：“https://dev59.com/rXLYa4cB1Zd3GeqPXnD6#16704079 表示，这种计算可以使用 * 和 sum 函数来完成，例如：”

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

---

经过进一步测试，我认为inner1d和matrix_multiply与您的dot和matrix-matrix乘积情况相匹配，并且如果添加[...,None]到matrix-vector情况中，则也匹配。看起来它们比einsum版本快2倍（在我的机器和测试数组上）。 https://github.com/numpy/numpy/blob/master/doc/neps/return-of-revenge-of-matmul-pep.rst是关于numpy中@中缀运算符的讨论。我认为numpy开发人员对这个PEP的热情不如Python开发人员高。

- hpaulj

---

这里有一个很好的解释，说明为什么BLAS速度会变慢。我正在通过使用中间张量和BLAS调用来优化einsum表达式。我想知道这是否是我们应该考虑的情况。在此处查看PR。 - Daniel

去年在解决einsum（涉及省略号）bug时，我编写了一个纯Python模拟器。它可以帮助跟踪索引字符串的解析过程。代码位于 https://github.com/hpaulj/numpy-einsum。 - hpaulj