如何有效地对两个稀疏 SciPy 矩阵进行矩阵乘积并生成密集的 Numpy 数组？

Question

我想执行两个稀疏 SciPy 矩阵的矩阵乘积。然而，就我而言，结果并不稀疏，所以我想将其存储为密集的 NumPy 数组。

是否可以有效地做到这一点，即无需先创建“稀疏”矩阵然后对其进行转换？我可以自由选择任何输入格式（以效率更高者为准）。

---

一个例子：两个具有随机分布零的 10000x10000 99% 稀疏矩阵的乘积将是稠密的：

n = 10_000
a = np.random.randn(n, n) * (np.random.randint(100, size=(n, n)) == 0)
b = np.random.randn(n, n) * (np.random.randint(100, size=(n, n)) == 0)
c = a.dot(b)
np.count_nonzero(c) / c.size # should be 0.63

Answer 1

import numpy as np
from scipy import sparse

n = 10_000
a = sparse.csr_matrix(np.random.randn(n, n) * (np.random.randint(100, size=(n, n)) == 0))
b = sparse.csr_matrix(np.random.randn(n, n) * (np.random.randint(100, size=(n, n)) == 0))
c = a.dot(b)

>>> c
<10000x10000 sparse matrix of type '<class 'numpy.float64'>'
    with 63132806 stored elements in Compressed Sparse Row format>

是的，这是存储这个矩阵的一种非常低效的方法。不过 scipy 中没有办法直接进入密集状态。您可以使用直接进入密集的稀疏BLAS 函数（该函数存在于您所描述的用例中）。

我为此使用了一个用于 MKL的 python 包装器，它包装了 mkl_sparse_spmmd:

from sparse_dot_mkl import dot_product_mkl

c_dense = dot_product_mkl(a, b, dense=True)

>>>> np.sum(c_dense != 0)
63132806

它也是线程化的，所以它比 scipy 快得多。安装 MKL 留给读者（

conda install -c intel mkl

可能是最简单的）