使用稀疏矩阵进行NumPy逐元素外积

Question

使用稀疏矩阵进行NumPy逐元素外积

pythonnumpysparse-matrixmatrix-multiplicationelementwise-operations

3

我希望在Python中对三个（或四个）大型2D数组进行逐元素外积运算（值为float32，舍入到2位小数）。它们都有相同的行数“n”，但不同的列数“i”、“j”、“k”。
结果数组应该具有形状（n，i*j*k）。然后，我想对结果的每一列求和，以得到形状为（i*j*k）的1D数组。

np.shape(a) = (75466, 10)
np.shape(b) = (75466, 28)
np.shape(c) = (75466, 66)

np.shape(intermediate_result) = (75466, 18480)
np.shape(result) = (18480)

感谢ruankesi和divakar，我得到了一段可行的代码：

链接

关于numpy元素级外积。

# Multiply first two matrices
first_multi = a[...,None] * b[:,None]
# could use np.einsum('ij,ik->ijk',a,b), which is slightly faster
ab_fills = first_multi.reshape(a.shape[0], a.shape[1]*b.shape[1])

# Multiply the result with the third matrix
second_multi = ab_fills[..., None] * c[:,None]
abc_fills = second_multi.reshape(ab_fills.shape[0], ab_fills.shape[1] * c.shape[1])

# Get the result: sum columns and get a 1D array of length 10*28*66 = 18 480
result = np.sum(abc_fills, axis = 0)

问题1: 性能

每次操作大约需要3秒，但我必须重复多次这个操作，而有些矩阵甚至更大（行数）。虽然可以接受，但让它变快会更好。

问题2: 我的矩阵是稀疏的

确实，例如，“a”包含70％的0。我尝试使用scipy csc_matrix，但确实无法得到一个可用的版本。（为了在此处获取逐元素外积，我通过转换为3D矩阵进行，但scipy sparse_matrix不支持）

问题3: 内存使用

如果我尝试处理第4个矩阵，我会遇到内存问题。

我想象中将这段代码转换为sparse_matrix会节省很多内存，并通过忽略众多0值使计算更快。这是真的吗？如果是，有人能帮帮我吗？当然，如果你有更好的实现建议，我也非常感兴趣。我不需要任何中间结果，只需要最终的1D结果。我被卡在这段代码上已经几周了，我要疯了！

谢谢！

Divakar答案后的编辑

方法#1：
非常好的一行代码但比原始方法慢得出乎意料（？）。
在我的测试数据集上，方法#1每次循环需要4.98秒±3.06毫秒（使用optimize = True没有加速）
原始分解方法需要3.01秒±16.5毫秒每次循环

方法#2：
绝对棒，谢谢！真是一个令人印象深刻的加速！
每次循环需要62.6毫秒±233微秒

关于numexpr，我尽可能避免对外部模块的要求，并且不打算使用多核/线程。这是一个“令人尴尬”的可并行任务，有成千上万的对象需要分析，在生产过程中我将把列表分散到可用的CPU上。我会尝试使用它来进行内存优化。
作为使用限制为1个线程，执行1个乘法的numexpr的简短尝试中，我在没有numexpr的情况下获得了40ms的运行时间，而使用numexpr则为52ms。
再次感谢！！

- plancton

发布的解决方案对您有用吗？ - Divakar

太好了，太棒了。谢谢！我正在编辑帖子，加上更长的答案。 - plancton

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Divakar · Accepted Answer

方法一

我们可以使用np.einsum来一次性进行求和缩减 -

result = np.einsum('ij,ik,il->jkl',a,b,c).ravel()

此外，尝试使用np.einsum中的optimize标志，将其设置为True以使用BLAS。

方法#2：

我们可以使用broadcasting来执行第一步，正如发布的代码中所提到的那样，然后利用np.tensordot进行张量矩阵乘法。

def broadcast_dot(a,b,c):
    first_multi = a[...,None] * b[:,None]
    return np.tensordot(first_multi,c, axes=(0,0)).ravel()

我们还可以使用支持多核处理并实现更好的内存效率的 numexpr 模块来获取 first_multi。这为我们提供了一个修改后的解决方案，如下所示 -

import numexpr as ne

def numexpr_broadcast_dot(a,b,c):
    first_multi = ne.evaluate('A*B',{'A':a[...,None],'B':b[:,None]})
    return np.tensordot(first_multi,c, axes=(0,0)).ravel()

给定数据集大小的随机浮点数据的时间 -

In [36]: %timeit np.einsum('ij,ik,il->jkl',a,b,c).ravel()
4.57 s ± 75.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [3]: %timeit broadcast_dot(a,b,c)
270 ms ± 103 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [4]: %timeit numexpr_broadcast_dot(a,b,c)
172 ms ± 63.8 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

仅为了展示使用 numexpr 的改进感受 -

In [7]: %timeit a[...,None] * b[:,None]
80.4 ms ± 2.64 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [8]: %timeit ne.evaluate('A*B',{'A':a[...,None],'B':b[:,None]})
25.9 ms ± 191 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

这在将该解决方案扩展到更多输入时应该是实质性的。