csr_matrix加法时如何计算索引？答案

【问题标题】：how to calculate indices when doing csr_matrix addition?csr_matrix加法时如何计算索引？
【发布时间】：2019-05-24 12:22:49
【问题描述】：

我正在学习稀疏矩阵运算的 scipy 代码。但是，我对在进行加法运算时如何计算 csr_matrix 索引感到困惑。我的代码很简单：

B = A + A.T
A.toarray()
array([[1., 1., 1., 1., 0., 1., 1., 0., 1., 1., 1., 1.],
   [0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.],
   [0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.],
   [0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.],
   [0., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
   [0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.],
   [0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.],
   [0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.],
   [0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.],
   [0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.],
   [0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.],
   [0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.]])
A.indices:
array([ 0,  2,  1,  5, 10,  3,  6, 11,  9,  8,  
    1,  2,  9, 10,  5,  8,  6, 11,  3,  7,  
    2, 10,  1,  5,  9,  6, 11,  7,  8,  3,  
    3,  8,  6,  7,  5, 11,  9, 10,  2,  1,  
    4,  9,  1,  7,  8, 11,  3,  6, 10,  5,  
    5,  6, 11, 10,  8,  9,  3,  2,  7,  1,  
    6,  5,  8, 11,  9, 10,  3,  7,  2,  1,  
    7,  9,  3,  6,  8,  5, 11, 10,  2,  1,  
    8,  6, 11,  9,  5,  3, 10,  7,  1,  2,  
    9, 11, 10,  8,  6,  5,  7,  3,  2,  1, 
    10, 11, 5,  9,  6,  8,  2,  3,  1,  7, 
    11, 10,  5,  6,  8,  9,  3,  7,  2, 1], dtype=int32)

我调试了深入 scipy 的代码以获取更多详细信息，发现 .T 操作使矩阵 A 变为 CSC 格式矩阵。然后，除了重载操作，A.T 将再次转换为 CSR 格式矩阵，但其索引更改如下：

A.T.toarray():
array([[1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
   [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
   [1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.],
   [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
   [0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
   [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
   [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
   [0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
   [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
   [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
   [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
   [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]])
A.T.indices:
array([ 0,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11,  0,  1,  2,  3,
    5,  6,  7,  8,  9, 10, 11,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9,
   10, 11,  4,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11,  0,  1,
    2,  3,  4,  5,  6,  7,  8,  9, 10, 11,  1,  2,  3,  4,  5,  6,  7,
    8,  9, 10, 11,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11,  0,
    1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11,  0,  1,  2,  3,  4,  5,
    6,  7,  8,  9, 10, 11,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10,
   11], dtype=int32)

我无法理解的是 B.indices，如下所示：

B.toarray()
array([[2., 1., 1., 1., 0., 1., 1., 0., 1., 1., 1., 1.],
   [1., 2., 2., 2., 1., 2., 2., 2., 2., 2., 2., 2.],
   [1., 2., 2., 2., 0., 2., 2., 2., 2., 2., 2., 2.],
   [1., 2., 2., 2., 1., 2., 2., 2., 2., 2., 2., 2.],
   [0., 1., 0., 1., 2., 1., 1., 1., 1., 1., 1., 1.],
   [1., 2., 2., 2., 1., 2., 2., 2., 2., 2., 2., 2.],
   [1., 2., 2., 2., 1., 2., 2., 2., 2., 2., 2., 2.],
   [0., 2., 2., 2., 1., 2., 2., 2., 2., 2., 2., 2.],
   [1., 2., 2., 2., 1., 2., 2., 2., 2., 2., 2., 2.],
   [1., 2., 2., 2., 1., 2., 2., 2., 2., 2., 2., 2.],
   [1., 2., 2., 2., 1., 2., 2., 2., 2., 2., 2., 2.],
   [1., 2., 2., 2., 1., 2., 2., 2., 2., 2., 2., 2.]])
B.indices
array([ 8,  9, 11,  6,  3, 10,  5,  1,  2,  0,  
    4,  0,  7,  3, 11,  6,  8,  5, 10,  9,  2,  1,  
    0,  3,  8,  7, 11,  6,  9,  5,  1, 10,  2,  
    4,  0,  1,  2, 10,  9, 11,  5,  7,  6,  8,  3,  
    5, 10,  6,  3, 11,  8,  7,  1,  9,  4,
    4,  0,  1,  7,  2,  3,  9,  8, 10, 11,  6,  5,  
    4,  0,  1,  2,  7,  3, 10,  9, 11,  8,  5,  6,  
    4,  1,  2, 10, 11,  5,  8,  6,  3,  9,  7,  
    4,  0,  2,  1,  7, 10,  3,  5,  9, 11,  6,  8,
    4,  0,  1,  2,  3,  7,  5,  6,  8, 10, 11,  9,  
    4,  0,  7,  1,  3,  2,  8,  6,  9,  5, 11, 10,  
    4,  0,  1,  2,  7,  3,  9,  8,  6,  5, 10, 11], dtype=int32)

结果 B 是通过 C++ 计算的，所以库 _sparsetools.cpython-35m-x86_64-linux-gnu.so 在 scipy 中，因此我无法获取它的代码。

这个问题几乎把我逼疯了。希望你们中的一些人可以帮助我。

【问题讨论】：

我认为你把事情弄得太复杂了。您的A 比测试所需的更大更密集。 A+A.T 的加法与两个不相关矩阵的加法（只要它们在形状上匹配）B = A1 + A2 没有任何不同。一个逻辑方法，假设两者都是csr，是逐行迭代的。 B 行的索引将是 A1 和 A2 的并集。现在可以删除等于 0 的总和，并在最后的 B.eliminate_zeros() 调用中删除。
您的A.T.indices 实际上是A.T.tocsr().indices，不是吗？使csr 副本对索引进行排序。

标签： python numpy scipy sparse-matrix

【解决方案1】：

如果我从您的A.todense() 显示中创建A：

In [156]: A = np.array([[1., 1., 1., 1., 0., 1., 1., 0., 1., 1., 1., 1.],
     ...:    [0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.],
     ...:    [0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.],
     ...:    [0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.],
     ...:    [0., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
     ...:    [0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.],
     ...:    [0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.],
     ...:    [0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.],
     ...:    [0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.],
     ...:    [0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.],
     ...:    [0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.],
     ...:    [0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.]])

还有一个稀疏的：

In [159]: M = sparse.csr_matrix(A)
In [160]: M
Out[160]: 
<12x12 sparse matrix of type '<class 'numpy.float64'>'
    with 120 stored elements in Compressed Sparse Row format>

In [162]: M.indices
Out[162]: 
array([ 0,  1,  2,  3,  5,  6,  8,  9, 10, 11,  1,  2,  3,  5,  6,  7,  8,
        9, 10, 11,  1,  2,  3,  5,  6,  7,  8,  9, 10, 11,  1,  2,  3,  5,
        6,  7,  8,  9, 10, 11,  1,  3,  4,  5,  6,  7,  8,  9, 10, 11,  1,
        2,  3,  5,  6,  7,  8,  9, 10, 11,  1,  2,  3,  5,  6,  7,  8,  9,
       10, 11,  1,  2,  3,  5,  6,  7,  8,  9, 10, 11,  1,  2,  3,  5,  6,
        7,  8,  9, 10, 11,  1,  2,  3,  5,  6,  7,  8,  9, 10, 11,  1,  2,
        3,  5,  6,  7,  8,  9, 10, 11,  1,  2,  3,  5,  6,  7,  8,  9, 10,
       11], dtype=int32)

这些索引与您的不同 - 或者它们相同但按行排序：

In [164]: for i in range(12):
     ...:     print(M.indices[M.indptr[i]:M.indptr[i+1]])
     ...:     
[ 0  1  2  3  5  6  8  9 10 11]
[ 1  2  3  5  6  7  8  9 10 11]
[ 1  2  3  5  6  7  8  9 10 11]
[ 1  2  3  5  6  7  8  9 10 11]
[ 1  3  4  5  6  7  8  9 10 11]
[ 1  2  3  5  6  7  8  9 10 11]
[ 1  2  3  5  6  7  8  9 10 11]
[ 1  2  3  5  6  7  8  9 10 11]
[ 1  2  3  5  6  7  8  9 10 11]
[ 1  2  3  5  6  7  8  9 10 11]
[ 1  2  3  5  6  7  8  9 10 11]
[ 1  2  3  5  6  7  8  9 10 11]

对于其转置的csr 版本：

In [165]: M1 = M.T.tocsr()
In [166]: for i in range(12):
     ...:     print(M1.indices[M1.indptr[i]:M1.indptr[i+1]])

[0]
[ 0  1  2  3  4  5  6  7  8  9 10 11]
[ 0  1  2  3  5  6  7  8  9 10 11]
[ 0  1  2  3  4  5  6  7  8  9 10 11]
[4]
[ 0  1  2  3  4  5  6  7  8  9 10 11]
[ 0  1  2  3  4  5  6  7  8  9 10 11]
[ 1  2  3  4  5  6  7  8  9 10 11]
[ 0  1  2  3  4  5  6  7  8  9 10 11]
[ 0  1  2  3  4  5  6  7  8  9 10 11]
[ 0  1  2  3  4  5  6  7  8  9 10 11]
[ 0  1  2  3  4  5  6  7  8  9 10 11]

通过切换到具有相同属性的csc 格式进行转置速度很快，但正如您注意到的那样，它已转换回csr。

然后取总和：

In [167]: B = M+M.T
In [168]: for i in range(12):
     ...:     print(B.indices[B.indptr[i]:B.indptr[i+1]])

[ 0  1  2  3  5  6  8  9 10 11]
[ 0  1  2  3  4  5  6  7  8  9 10 11]
[ 0  1  2  3  5  6  7  8  9 10 11]
[ 0  1  2  3  4  5  6  7  8  9 10 11]
[ 1  3  4  5  6  7  8  9 10 11]
[ 0  1  2  3  4  5  6  7  8  9 10 11]
[ 0  1  2  3  4  5  6  7  8  9 10 11]
[ 1  2  3  4  5  6  7  8  9 10 11]
[ 0  1  2  3  4  5  6  7  8  9 10 11]
[ 0  1  2  3  4  5  6  7  8  9 10 11]
[ 0  1  2  3  4  5  6  7  8  9 10 11]
[ 0  1  2  3  4  5  6  7  8  9 10 11]

由于M 索引按行排序，因此模式相当明显。

In [172]: M.has_sorted_indices
Out[172]: True

我的版本也是canonical。

A 的所有行都有 10 个非零项。 B 有 10、11、12；用 0-2 个非零填充的转置。

您的A 具有未排序的索引。 A.T.indices 已排序。我想人们可以通过比较无序的B.indices 和A.indices 来推断出一些关于评估方法的信息。但是这值得吗？

===

我重新创建了您的非规范矩阵：

In [186]: altindices = np.array([ 0,  2,  1,  5, 10,  3,  6, 11,  9,  8,  
     ...:     1,  2,  9, 10,  5,  8,  6, 11,  3,  7,  
     ...:     2, 10,  1,  5,  9,  6, 11,  7,  8,  3,  
     ...:     3,  8,  6,  7,  5, 11,  9, 10,  2,  1,  
     ...:     4,  9,  1,  7,  8, 11,  3,  6, 10,  5,  
     ...:     5,  6, 11, 10,  8,  9,  3,  2,  7,  1,  
     ...:     6,  5,  8, 11,  9, 10,  3,  7,  2,  1,  
     ...:     7,  9,  3,  6,  8,  5, 11, 10,  2,  1,  
     ...:     8,  6, 11,  9,  5,  3, 10,  7,  1,  2,  
     ...:     9, 11, 10,  8,  6,  5,  7,  3,  2,  1, 
     ...:     10, 11, 5,  9,  6,  8,  2,  3,  1,  7, 
     ...:     11, 10,  5,  6,  8,  9,  3,  7,  2, 1], dtype='int32')
    In [188]: M2 = sparse.csr_matrix((M.data, altindices, M.indptr))

M2.T.indices 与 M2.indices 相同（csc 等效项），但 M2.T.tocsr() 在显示时进行排序。 B2 = M2+M2.T 有你的索引。

如果我只评估一行，我会得到相同的索引，这证实了我的猜测，它逐行执行求和（使用indptr 进行迭代，如上所述）：

In [194]: M2[0,:]+(M2.T)[0,:]
Out[194]: 
<1x12 sparse matrix of type '<class 'numpy.float64'>'
    with 10 stored elements in Compressed Sparse Row format>
In [195]: _.indices
Out[195]: array([ 8,  9, 11,  6,  3, 10,  5,  1,  2,  0], dtype=int32)

【讨论】：