如何在 numpy、scipy 中进行多维矩阵乘法

Question

我有一个形状为 (8,9) 的 scipy 稀疏矩阵和另一个形状为 (9, 12 , 17) 的数组。我想将这些相乘，得到一个大小为 (8,12,17) 的矩阵/数组，其中 (8,9) 矩阵仅有效地乘以第一个维度。我是否必须使用 Kronecker 产品才能做到这一点，或者在 numpy 中有简单的方法吗？

Answer 1

如果m1 是二维稀疏矩阵，m1.A 是它的密集数组形式。维度实际上写了einsum 表达式。

np.einsum('ij,jkl->ikl', m1.A, m2)

例如：

In [506]: M = sparse.random(8, 9, 0.1)
In [507]: A = np.ones((9, 12, 17))
In [508]: np.einsum('ij,jkl->ikl', M.A, A).shape
Out[508]: (8, 12, 17)

Answer 2

正如 hpaulj 在 cmets 中建议的那样，最简单的方法是使用密集矩阵的 np.einsum：

>>> a = np.random.randn(8, 9)
>>> b = np.random.randn(9, 12, 17)
>>> c = np.einsum('ij,jkl->ikl', a, b)
>>> c.shape
(8, 12, 17)

Answer 3

@Divakar 推荐np.tensordot，@hpaulj 和@Praveen 推荐np.einsum。还有一种方法是转轴：

(a @ b.transpose((2, 0, 1))).transpose((1, 2, 0))

对于您引用的小尺寸，np.einsum 和转置似乎更快。但是，一旦您开始放大要沿其相乘的轴的维度，np.tensordot 就会胜过其他两个。

import numpy as np

m, n, k, l = 8, 9, 12, 17
a = np.random.random((m, n))
b = np.random.random((n, k, l))

%timeit np.tensordot(a, b, axes=([1], [0]))
# => 10000 loops, best of 3: 22 µs per loop
%timeit np.einsum("ij,jkl->ikl", a, b)
# => 100000 loops, best of 3: 10.1 µs per loop
%timeit (a @ b.transpose((2, 0, 1))).transpose((1, 2, 0))
# => 100000 loops, best of 3: 11.1 µs per loop

m, n, k, l = 8, 900, 12, 17
a = np.random.random((m, n))
b = np.random.random((n, k, l))

%timeit np.tensordot(a, b, axes=([1], [0]))
# => 1000 loops, best of 3: 198 µs per loop
%timeit np.einsum("ij,jkl->ikl", a, b)
# => 1000 loops, best of 3: 868 µs per loop
%timeit (a @ b.transpose((2, 0, 1))).transpose((1, 2, 0))
# => 1000 loops, best of 3: 907 µs per loop

m, n, k, l = 8, 90000, 12, 17
a = np.random.random((m, n))
b = np.random.random((n, k, l))

%timeit np.tensordot(a, b, axes=([1], [0]))
# => 10 loops, best of 3: 21.7 ms per loop
%timeit np.einsum("ij,jkl->ikl", a, b)
# => 10 loops, best of 3: 164 ms per loop
%timeit (a @ b.transpose((2, 0, 1))).transpose((1, 2, 0))
# => 10 loops, best of 3: 166 ms per loop

Answer 4

这里有几种方法可以让它发挥作用。第二个似乎更好，并且在我测试时快了大约 12 倍。

def multiply_3D_dim_zero_slow(matrix, array):
shape = array.shape
final_shape = (matrix.shape[0], array.shape[1], array.shape[2])
result = np.zeros(final_shape)
for i in xrange(shape[1]):
    for j in xrange(shape[2]):
        result[:, i, j] = matrix * array[:, i, j]
return result.reshape(final_shape)

这是一个更快的版本，它使用 reshape 将多维数组变成二维数组。

def multiply_3D_dim_zero(matrix, array):
shape = array.shape
final_shape = (matrix.shape[0], array.shape[1], array.shape[2])
array_reshaped = array.reshape(shape[0], shape[1] * shape[2])
return (matrix * array_reshaped).reshape(final_shape)ode here

这仅适用于我需要的第一个维度，但可以概括。