从 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)

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

对于上面的矩阵，密集的左矩阵和下矩阵看起来像：

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] 来稍微优化它。

请注意，我只是将索引向上移动一位，而不是将 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

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

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]]

Answer 3

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，每个 'row' 的索引可能未排序 - 有一种稀疏方法可以解决这个问题。

至于构建 return 稀疏矩阵，您可以修改副本的 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 4

更向量化的方法：

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 表示得到。

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

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

Get Left, Right, Up, Down Nonzero Neighbor from SciPy Sparse Matrix

python

numpy

matrix

scipy

sparse-matrix