scipy稀疏矩阵和导致密集矩阵-如何加强结果稀疏性？

Question

在scipy.sparse.csr_matrix的一个轴上求和会产生一个numpy.matrix对象。考虑到我的稀疏矩阵真的很稀疏，我发现这种行为非常恼人。

下面是一个示例：

dense = [[ 0.,  0.,  0.,  0.,  0.],
         [ 1.,  0.,  0.,  0.,  0.],
         [ 0.,  0.,  0.,  0.,  0.],
         [ 0.,  0.,  0.,  0.,  0.],
         [ 2.,  0.,  4.,  0.,  0.]]


from scipy.sparse import csr_matrix
sparse = csr_matrix(dense)

print(sparse.sum(1))

with结果：

matrix([[ 0.],
        [ 1.],
        [ 0.],
        [ 0.],
        [ 6.]])

如何在不将矩阵隐式转换为密集格式的情况下，在列求和操作中强制使用稀疏性？在这个例子中，我只使用了一个小的n矩阵，但是我的矩阵要大得多，也更稀疏，所以通过密集表示法是一种很大的空间浪费。

Answer 1

sparse使用矩阵乘法执行求和：

In [136]: np.matrix(np.ones(M.shape[1]))@M                                      
Out[136]: matrix([[3., 0., 4., 0., 0.]])
In [137]: M@np.matrix(np.ones((M.shape[1],1)))                                  
Out[137]: 
matrix([[0.],
        [1.],
        [0.],
        [0.],
        [6.]])
In [138]: timeit M@np.matrix(np.ones((M.shape[1],1)))                           
91.5 µs ± 268 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [139]: timeit M.sum(1)                                                       
96.6 µs ± 647 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

时代是相似的。这两种方法都会产生np.matrix结果。

如果乘以二维数组，我会得到一个数组结果，有点出人意料的是，时间要好得多：

In [140]: timeit M@np.ones((M.shape[1],1))                                      
24.4 µs ± 1.09 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [141]: M@np.ones((M.shape[1],1))                                             
Out[141]: 
array([[0.],
       [1.],
       [0.],
       [0.],
       [6.]])

我可以把这个数组放回一个稀疏矩阵中--但需要花费一定的时间：

In [142]: csr_matrix(M@np.ones((M.shape[1],1)))                                 
Out[142]: 
<5x1 sparse matrix of type '<class 'numpy.float64'>'
    with 2 stored elements in Compressed Sparse Row format>
In [143]: timeit csr_matrix(M@np.ones((M.shape[1],1)))                          
391 µs ± 17.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

或者我们可以先创建一个稀疏矩阵：

In [144]: M@csr_matrix(np.ones((M.shape[1],1)))                                 
Out[144]: 
<5x1 sparse matrix of type '<class 'numpy.float64'>'
    with 2 stored elements in Compressed Sparse Row format>
In [145]: timeit M@csr_matrix(np.ones((M.shape[1],1)))                          
585 µs ± 5.28 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

即使从循环中删除提取器矩阵的创建也会导致速度变慢：

In [146]: %%timeit m1 = csr_matrix(np.ones((M.shape[1],1))) 
     ...: M@m1                                                                     
227 µs ± 4.72 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

像这样的sum (几乎)总是增加结果的密度。每行至少有一个非零值的矩阵比具有许多纯零行的矩阵更常见。在现实世界中，计时可能会有所不同，但尝试使用相同的内存可能不会为您带来太多好处。

如果我更详细地查看由稀疏矩阵乘法生成的csr矩阵：

In [147]: res = M@csr_matrix(np.ones((M.shape[1],1)))                           
In [148]: res                                                                   
Out[148]: 
<5x1 sparse matrix of type '<class 'numpy.float64'>'
    with 2 stored elements in Compressed Sparse Row format>
In [149]: res.indptr                                                            
Out[149]: array([0, 0, 1, 1, 1, 2], dtype=int32)
In [150]: res.indices                                                           
Out[150]: array([0, 0], dtype=int32)

indptr数组每行有一个值(+1)，因此此列矩阵的内存使用量实际上高于密集等效数组。同样的csc格式的res会更紧凑，只有2个元素的indptr。

还可以直接使用csr矩阵的indptr、indices、data属性，本质上是对行进行迭代并对每个行求和，然后从中创建一个新的稀疏矩阵。在某些情况下，与sparse方法相比，我们已经实现了速度改进。但你必须理解数据存储，并对整个事情保持智能。