具有稀疏矩阵的 numpy 元素外积
numpy elementwise outer product with sparse matrices
我想对 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)
感谢 ，我得到了一段有效的代码：
# 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 s ± 3.06 ms（优化 = True 时没有加速）
最初的分解方法每个循环花费 3.01 s ± 16.5 ms
方法 #2：
绝对好，谢谢！多么令人印象深刻的加速！
每个循环 62.6 毫秒 ± 233 微秒
关于numexpr，我尽量避免对外部模块的要求，也不打算使用multicores/threads。这是一个 "embarrassingly" 可并行任务，有数十万个对象需要分析，我将在生产期间将列表分散到可用的 CPU 上。我会尝试一下内存优化。
作为 numexpr 的简短尝试，限制为 1 个线程，执行 1 次乘法，我得到 运行没有 numexpr 的时间为 40 毫秒，有 numexpr 的时间为 52 毫秒。
再次感谢！！
方法 #1
我们可以用np.einsum一次性完成sum-reductions-
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()
我们也可以使用支持multi-core处理的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)
将此解决方案扩展到更多输入时，这应该很重要。
我想对 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)
感谢
# 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 s ± 3.06 ms（优化 = True 时没有加速）
最初的分解方法每个循环花费 3.01 s ± 16.5 ms
方法 #2：
绝对好，谢谢！多么令人印象深刻的加速！
每个循环 62.6 毫秒 ± 233 微秒
关于numexpr，我尽量避免对外部模块的要求，也不打算使用multicores/threads。这是一个 "embarrassingly" 可并行任务，有数十万个对象需要分析，我将在生产期间将列表分散到可用的 CPU 上。我会尝试一下内存优化。
作为 numexpr 的简短尝试，限制为 1 个线程，执行 1 次乘法，我得到 运行没有 numexpr 的时间为 40 毫秒，有 numexpr 的时间为 52 毫秒。
再次感谢！！
方法 #1
我们可以用np.einsum一次性完成sum-reductions-
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()
我们也可以使用支持multi-core处理的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)
将此解决方案扩展到更多输入时，这应该很重要。