从 SciPy 稀疏矩阵获取左、右、上、下非零邻居

Question

假设我有一个 2D SciPy 稀疏矩阵：

import numpy as np
from scipy.sparse import csc_matrix

arr = np.array([[0, 0, 1, 0, 1],
                [1, 0, 0, 1, 0],
                [0, 1, 1, 0, 0],
                [1, 0, 0, 1, 0],
                [0, 1, 0, 0, 0],
               ])

csc = csc_matrix(arr)

对于矩阵中的每个非零元素，我想创建四个新的稀疏矩阵，其中包含与下一个最近的 Left、Right、Up 和 Down 非零邻居对应的索引。末端的元素可以有环绕的邻居（想想水平和垂直方向的圆形双向链表或环形）。如果一个元素是其行/列中唯一的非零元素，则相应的索引将指向它自己。此外，由于索引可以具有零值（当引用第一行或第一列时）并且与自然为零的元素无法区分，我们将这些零索引设置为 -1 以消除实际索引与零元素的歧义。

对于上面的矩阵，密集的 Left 和 Down 矩阵如下所示：

left = np.array([[0, 0, 4,  0, 2],
                 [3, 0, 0, -1, 0],
                 [0, 2, 1,  0, 0],
                 [3, 0, 0, -1, 0],
                 [0, 1, 0,  0, 0],
                ])

down = np.array([[0, 0,  2, 0, -1],
                 [3, 0,  0, 3,  0],
                 [0, 4, -1, 0,  0],
                 [1, 0,  0, 1,  0],
                 [0, 2,  0, 0,  0],
                ])

请记住，索引值为 -1 的元素实际上是对索引零的引用。当然，我需要将这些矩阵设为稀疏矩阵形式，因为我的真实矩阵太大太稀疏，无法放入内存。

Answer 1

这是一种可能的方法来做左边的邻居。 它不是特别有效，但如果整个矩阵中没有很多非零条目，它可能工作得很好。您可以通过获取每行的非零条目并仅计算一次 j[i==row] 来稍微优化它。

请注意，我只是将索引上移 1，而不是将 0 设置为 -1。

i,j = csc.nonzero()
ind = sp.sparse.csc_matrix(csc.shape,dtype='int')
for row in range(csc.shape[0]):
    ind[row,j[i==row]] = np.roll(j[i==row]+1,1)

ind.A = array([[0, 0, 5, 0, 3],
   [4, 0, 0, 1, 0],
   [0, 3, 2, 0, 0],
   [4, 0, 0, 1, 0],
   [0, 2, 0, 0, 0]])

Answer 2

In [183]: arr = np.array([[0, 0, 1, 0, 1],
     ...:                 [1, 0, 0, 1, 0],
     ...:                 [0, 1, 1, 0, 0],
     ...:                 [1, 0, 0, 1, 0],
     ...:                 [0, 1, 0, 0, 0],
     ...:                ])
     ...:                
In [184]: from scipy import sparse
In [185]: M = sparse.lil_matrix(arr)
In [186]: M.rows
Out[186]: 
array([list([2, 4]), list([0, 3]), list([1, 2]), list([0, 3]), list([1])],
      dtype=object)

这与您从密集数组中获得的信息相同：

In [187]: [np.where(row)[0] for row in arr]
Out[187]: [array([2, 4]), array([0, 3]), array([1, 2]), array([0, 3]), array([1])]

我假设您已经知道如何从密集数组生成所需的left（或right），所以我不会深入这些细节（我懒得与您的包装规范搏斗)。

对于列：

 In [189]: M.T.rows
 Out[189]: 
 array([list([1, 3]), list([2, 4]), list([0, 2]), list([1, 3]), list([0])],
  dtype=object)

您可以使用csc 格式：

In [190]: Mc = sparse.csc_matrix(arr)
In [191]: Mc.indptr
Out[191]: array([0, 2, 4, 6, 8, 9], dtype=int32)
In [192]: Mc.indices
Out[192]: array([1, 3, 2, 4, 0, 2, 1, 3, 0], dtype=int32)
In [193]: for i in range(5):
     ...:     print(Mc.indices[Mc.indptr[i]:Mc.indptr[i+1]])
     ...:     
[1 3]
[2 4]
[0 2]
[1 3]
[0]

在此示例中，所有行或列都只有 1 或 2 个非零值。我想在更大更一般的情况下会有很多非零值。同样对于csc（和csr，每个“行”的索引可能没有排序 - 有一种稀疏方法可以解决这个问题。

至于构建返回稀疏矩阵，您可以修改副本的data 属性（它将具有相同的稀疏度）。

In [194]: M.data
Out[194]: 
array([list([1, 1]), list([1, 1]), list([1, 1]), list([1, 1]), list([1])],
      dtype=object)
In [195]: Mc.data
Out[195]: array([1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=int64)

或者从数组构造一个稀疏矩阵（对于coo 样式输入是正常的）。

---

使用我的lil 版本，tch's 解决方案稍微快一些：

ind = sparse.lil_matrix(M.shape,dtype='int')
for i,row in enumerate(M.rows):
    k = np.array(row)
    ind[i,k] = np.roll(k+1,1)

用我替换 data 的想法更好：

ind = M.copy()
for row,dat in zip(ind.rows,ind.data):
    k = np.array(row)
    dat[:] = np.roll(k+1,1).tolist()

或Mr = Mc.tocsr()

ind = Mr.copy()
for i in range(Mr.shape[0]):
    slc = slice(Mr.indptr[i],Mr.indptr[i+1])
    k = Mr.indices[slc]
    ind.data[slc] = np.roll(k+1,1)

Answer 3

更矢量化的方法：

csc = csc_matrix(arr)
inds = (csc.indices,csc.indptr)
irows = np.split(*inds)[1:-1]

down = csc_matrix((np.hstack([np.roll(row,-1) for row in irows]),*inds))
up = csc_matrix((np.hstack([np.roll(row,1) for row in irows]),*inds))

检查：

>>> down.A 
array([[0, 0, 2, 0, 0],
       [3, 0, 0, 3, 0],
       [0, 4, 0, 0, 0],
       [1, 0, 0, 1, 0],
       [0, 2, 0, 0, 0]], dtype=int32)

Left 和 Right 可以通过 CSR 表示获得。

我不认为将 0 编码为 -1 是一个好主意，因为 if 会破坏所有稀疏计算的改进。只有csc.nonzeros()设计的地方必须参观。

Answer 4

一个可能的答案（密集形式）：

ix, iy = csc.nonzero()
w = np.where(np.insert(np.diff(ix), 0,1) != 0)[0]
iy2 = np.concatenate([np.roll(_, 1) for _ in np.split(iy,w)])
iy2[iy2==0] = -1

left = csc_matrix(arr.shape)
left[ix, iy] = iy2

ix, iy = csc.transpose().nonzero()
w = np.where(np.insert(np.diff(ix), 0,1) != 0)[0]
iy2 = np.concatenate([np.roll(_, 1) for _ in np.split(iy,w)])
iy2[iy2==0] = -1

down = csc_matrix(arr.T.shape)
down[ix, iy] = iy2
down = down.transpose()
print(left.todense(), '\n', down.todense())


 >> [[ 0  0  4  0  2]
 [ 3  0  0 -1  0]
 [ 0  2  1  0  0]
 [ 3  0  0 -1  0]
 [ 0  1  0  0  0]]

[[ 0  0  2  0 -1]
 [ 3  0  0  3  0]
 [ 0  4 -1  0  0]
 [ 1  0  0  1  0]
 [ 0  2  0  0  0]]