具有稀疏矩阵的 numpy 元素外积

我想在 python 中对三个（或四个）大型二维数组进行逐元素外积（值是 float32 四舍五入到小数点后两位）。它们都具有相同的行数“n”，但具有不同的列数“i”、“j”、“k”。
所得数组的形状应为 (n, i*j*k)。然后，我想对结果的每一列求和，最终得到形状 (i*j*k) 的一维数组。

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)

谢谢阮克斯和迪瓦卡 https://stackoverflow.com/questions/42378936/numpy-elementwise-outer-product，我得到了一段有效的代码：

# 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稀疏矩阵中不受支持）

问题3： 内存使用情况

如果我尝试使用第四个矩阵，我会遇到内存问题。


我想将此代码转换为稀疏矩阵会节省大量内存，并通过忽略大量 0 值来加快计算速度。 真的吗？如果是，有人可以帮助我吗？
当然，如果您有什么更好的实施建议，我也很感兴趣。我不需要任何中间结果，只需要最终的一维结果。
我已经被困在这部分代码上几个星期了，我快疯了！

谢谢你！

在 Divakar 的回答后编辑

方法#1：
非常好的一个班轮，但令人惊讶地比原来的方法慢（？）。
在我的测试数据集上，方法 #1 每个循环需要 4.98 s ± 3.06 ms（优化 = True 时没有加速）
最初的分解方法每个循环花费 3.01 s ± 16.5 ms

方法#2：
绝对棒极了，谢谢！多么令人印象深刻的加速！
每个循环 62.6 ms ± 233 µs

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

方法#1

我们可以用np.einsum https://docs.scipy.org/doc/numpy/reference/generated/numpy.einsum.html一次性进行求和减少 -

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

另外，玩一下optimize标记在np.einsum通过将其设置为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 module https://github.com/pydata/numexpr/blob/master/doc/user_guide.rst#enabling-intel-vml-support支持多核处理并实现更好的内存效率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)

当将此解决方案扩展到更多数量的输入时，这应该是很重要的。