如何计算非常大的SciPy稀疏矩阵的点积

Question

我正在尝试计算一个很大的矩阵和它自己之间的点积。

这个矩阵的形状是 (371744, 36154)，也就是说它有371744行和36154列。里面有577731个非零元素，说明这个矩阵非常稀疏。

mat1是一个稀疏矩阵，使用的是scipy库的csr_matrix格式。如果我用mat1 * mat1.T来计算，就会出现一个值错误。这看起来是因为结果矩阵中的非零元素太多，导致索引指针溢出，具体可以参考这里。

    dp_data = data_m * data_m.T
  File "/usr/lib/python2.7/dist-packages/scipy/sparse/base.py", line 247, in __mul__
    return self._mul_sparse_matrix(other)
  File "/usr/lib/python2.7/dist-packages/scipy/sparse/base.py", line 300, in _mul_sparse_matrix
    return self.tocsr()._mul_sparse_matrix(other)
  File "/usr/lib/python2.7/dist-packages/scipy/sparse/compressed.py", line 290, in _mul_sparse_matrix
    indices = np.empty(nnz, dtype=np.intc)
ValueError: negative dimensions are not allowed

我也尝试过用np.dot来计算。

但是文档上说，
"从NumPy 1.7开始，np.dot并不知道稀疏矩阵，因此使用它会导致意想不到的结果或错误。应该先获取对应的稠密矩阵"

当我尝试用mat1.toarray()或todense()时，出现了内存错误，因为这个矩阵实在太大了！我只有16GB的内存！对于较小的输入，程序似乎运行得很好！

    data_array = data_m.toarray()
  File "/usr/lib/python2.7/dist-packages/scipy/sparse/compressed.py", line 550, in toarray
    return self.tocoo(copy=False).toarray()
  File "/usr/lib/python2.7/dist-packages/scipy/sparse/coo.py", line 219, in toarray
    B = np.zeros(self.shape, dtype=self.dtype)
MemoryError

我使用的版本是：
Numpy 1.8.1
Numpy 0.9.0

那我该如何进行这个乘法运算呢？

Answer 1

首先，稀疏输出矩阵的大小和计算所需的CPU工作量，取决于你的稀疏矩阵的结构。如果有很多空行，事情就简单多了。相反，如果你的数据分布得比较均匀，那就需要进行很多计算。

首先要明白的是，在这个特定的情况下（a * a.T），你实际上是在计算每一行（一个有36154个元素的向量）与每一行的点积。这可以帮助你将计算时间减少一半，因为结果是对称的。（这意味着你大约需要计算50,000,000,000个向量的点积。）

现在，问题是你是否着急。如果你很着急（性能很重要），那么根据你数据中非零元素的分布，可能会有一些方法可以加快计算速度。

一个相对简单的算法如下：

# a is in row-form (n_row lists of non-zeros in each row)

for r in 0..n_row-1
    if all elements in the row are zero
        skip the row
    for c in r..n_row-1
        d := dot(a[c],a[r])
        if d <> 0
            result[r,c] = d
            result[c,r] = d

点积的计算很简单，只需要找到a[c]和a[r]中非零元素集合的交集。大多数情况下，交集是空的，只有在交集不为空时才需要进行计算。

根据你矩阵中空行的数量，这个过程可能会相对快速。另一方面，如果你没有任何空行或列，那么就需要花费时间计算50,000,000,000个集合的交集。（在这种情况下，大多数交集都是1之间的比较，所以比较起来很简单。）

这种方法需要的内存很少，但仍然需要花费大量时间，可能从几个小时到几天，除非有很多空行。

Answer 2

调用稀疏矩阵的点积方法：

dp_data = data_m.dot(data_m)

numpy.dot 是一个叫做 通用函数，它并不知道你的矩阵是稀疏的。而 scipy.sparse.csc_matrix.dot 是一个专门为你的矩阵类型设计的方法，因此它使用了稀疏算法。