添加切片时,在 scipy 稀疏 csc 矩阵中,存储元素的数量无形地增长

Number of stored elements grows invisibly, in scipy sparse csc matrix, when adding slices

我正在尝试将 scipy 稀疏 csc 矩阵的切片加在一起。我现在这样做的方式是,输出了具有正确值的矩阵,但是由于一个神秘的原因,存储的元素数量增加了。它比数学上矩阵中非零元素的数量增长得更高。怎么样?

示例:

import numpy as np
from scipy.sparse import *

#define two dense matrices
M1_d = np.array([[0, 0, 1, 2, 3],
                 [4, 5, 0, 0, 0],
                 [6, 0, 0, 7, 8],
                 [9,10, 0,11,12],
                 [0, 0, 0, 0, 0]])
M2_d = 2*M1_d

#define their sparse counterparts (in CSC form)
M1_s = csc_matrix(M1_d)
M2_s = csc_matrix(M2_d)

#doing addition with both forms works fine
Mp_d = M1_d + M2_d
Mp_s = M1_s + M2_s

在此正常加法中,Mp_d 和 Mp_s 具有相同的结果,并且正确数量的元素存储在稀疏矩阵 Mp_s

>>>Mp_s
<5x5 sparse matrix of type '<type 'numpy.int32'>'
with 12 stored elements in Compressed Sparse Column format>

另一方面,当我尝试仅将矩阵的切片加在一起并将其存储在相同大小的矩阵中时。正确的值出来了,但是稀疏矩阵中存储了一些新元素

#initialise some matrices as copies
Mslicep_d = M1_d.copy()
Mslicep_s = M1_s.copy()

#now I want to add only slices of the matrices
clm = [0,2,3] # I want to add these colums
Mslicep_d[:,clm] = M1_d[:,clm] + M2_d[:,clm]
Mslicep_s[:,clm] = M1_s[:,clm] + M2_s[:,clm]

矩阵 Mslicep_s 应再次具有 12 个存储元素。看起来确实如此,但是:

>>>Mslicep_s
<5x5 sparse matrix of type '<type 'numpy.int32'>'
with 20 stored elements in Compressed Sparse Column format>
>>>Mslicep_s.toarray()
array([[ 0,  0,  3,  6,  3],
       [12,  5,  0,  0,  0],
       [18,  0,  0, 21,  8],
       [27, 10,  0, 33, 12],
       [ 0,  0,  0,  0,  0]])

如果你计算这里的值,你会发现有 12 个。但是由于某些原因,稀疏矩阵有 20 个存储元素。

当我在一个循环中执行多个这样的加法时,矩阵变得越来越稀疏,这使得我的代码需要很长时间。谁能告诉我这些新存储的元素是从哪里来的?如何防止这种情况?或者至少如何在每次迭代中快速修复它?

提前致谢。

编辑:让我说清楚。我想执行矩阵稀疏结构不变的加法。

我第一次做问题添加时收到警告:

In [344]: Mslicep_s[:,clm] = M1_s[:,clm] + M2_s[:,clm]
/usr/local/lib/python2.7/site-packages/scipy/sparse/compressed.py:728: SparseEfficiencyWarning: Changing the sparsity structure of a csc_matrix is expensive. lil_matrix is more efficient.
  SparseEfficiencyWarning)

那多出的8个元素全为0,可见Mslicep_s.data。显然它正在插入你的总和的密集版本:

In [355]: (M1_s[:,clm] + M2_s[:,clm]).A
Out[355]: 
array([[ 0,  3,  6],
       [12,  0,  0],
       [18,  0, 21],
       [27,  0, 33],
       [ 0,  0,  0]])

我猜它更喜欢速度而不是存储。

我清理矩阵:

In [347]: Mslicep_s.eliminate_zeros()    
In [348]: Mslicep_s
Out[348]: 
<5x5 sparse matrix of type '<type 'numpy.int32'>'
    with 12 stored elements in Compressed Sparse Column format>

因此,作为即时改进,您可以在每次迭代时执行此步骤。

对于密集数组,select使用 clm 等列表的列是 'advanced indexing',并生成副本,而不是更快的视图。 csc 格式是像这样 selecting 列的最佳稀疏格式,但 M1_s[:,clm] 仍然是一个副本。如您所见,将值插入 csc 矩阵既耗时又 space。我不确定更改 Mslicep_s 的格式是否有帮助。 lil 应该更适合修改稀疏性,尽管更改行比更改列更好。

请记住,csccsr 格式针对矩阵乘法计算进行了优化。他们甚至用乘法做 .sum() 之类的事情。一般来说,如果您可以用这些术语进行计算,您将获得更好的性能,无论是时间还是内存。索引不是稀疏矩阵的强项。

如果您所做的计算不改变稀疏性(不添加非零值),那么您可以直接对 .data 属性数组进行更改。但确实需要了解存储布局。

<todo: code to translate [:,clm] to idx array>
Mslice_s.data[idx] = M1_s.data[idx] + M2_s.data[idx]

这似乎是正确的 idx:

idx = np.concatenate([np.arange(M1_s.indptr[i], M1_s.indptr[i+1]) for i in clm])
# array([0, 1, 2, 5, 6, 7, 8])

可能有更直接的计算方法,但我认为它有正确的想法 - select 基于 selected .indptr 值的范围。