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

Question

我想在 python 中做三个（或四个）大型二维数组的元素外积（值是 float32 舍入到小数点后 2 位）。它们都有相同的行数“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)

感谢ruankesi and divakar，我得到了一段有效的代码：

# 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 值来加快计算速度。 真的吗？如果是，有人可以帮助我吗？
当然，如果你对更好的实现有什么建议，我也很感兴趣。我不需要任何中间结果，只需要最终的一维结果。
我被困在这部分代码已经好几个星期了，我快疯了！

谢谢！

在 Divakar 回答后编辑

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

方法二：
绝对很棒，谢谢！多么令人印象深刻的加速！
每个循环 62.6 毫秒 ± 233 微秒

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

Answer 1

方法#1

我们可以使用np.einsum一次性做减和-

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

另外，通过将np.einsum 中的optimize 标志设置为True 来使用BLAS。

方法#2

我们可以使用broadcasting 来做第一步，正如发布的代码中提到的那样，然后利用 tensor-matrix-multiplcation 和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，它支持多核处理，内存效率也更好，得到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)

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