scipy (稀疏)2D numpy数组的每行/列快速非零索引

xzlaal3s  于 2023-03-08  发布在  其他
关注(0)|答案(3)|浏览(176)

我正在寻找一种最快的方法来获取一个二维数组的每行和每列的非零索引列表。下面是一段代码:

preds = [matrix[:,v].nonzero()[0] for v in range(matrix.shape[1])]
descs = [matrix[v].nonzero()[0] for v in range(matrix.shape[0])]
    • 示例输入:**
matrix = np.array([[0,0,0,0],[1,0,0,0],[1,1,0,0],[1,1,1,0]])
    • 示例输出**
preds = [array([1, 2, 3]), array([2, 3]), array([3]), array([], dtype=int64)]
descs = [array([], dtype=int64), array([0]), array([0, 1]), array([0, 1, 2])]

(The列表被称为预测和描述,因为当矩阵被解释为邻接矩阵时,它们指的是DAG中的前趋和后代,但这对问题来说不是必需的。)

    • 时间示例:**出于时间目的,以下矩阵是一个很好的代表:
test_matrix = np.zeros(shape=(4096,4096),dtype=np.float32)
for k in range(16):
    test_matrix[256*(k+1):256*(k+2),256*k:256*(k+1)]=1
    • 背景:**在我的代码中,对于一个4000x4000的矩阵,这两行占用了75%的时间,而随后的拓扑排序和DP算法只占用了剩下的四分之一的时间。矩阵中大约5%的值是非零的,所以稀疏矩阵的解决方案可能是适用的。

谢谢你。
(On建议也张贴在这里:https://scicomp.stackexchange.com/questions/35242/fast-nonzero-indices-per-row-column-for-sparse-2d-numpy-array还有一些答案,我将在评论中提供时间。此链接包含一个可接受的答案,速度是原来的两倍。

vwoqyblh

vwoqyblh1#

如果你有足够的动力,Numba可以做出惊人的事情。下面是你所需要的逻辑的快速实现。简单地说,它计算np.nonzero()的等价物,但它同时包含了稍后将索引分派成你所需要的格式的信息。这些信息是受sparse.csr.indptrsparse.csc.indptr的启发而产生的。

import numpy as np
import numba as nb

@nb.jit
def cumsum(arr):
    result = np.empty_like(arr)
    cumsum = result[0] = arr[0]
    for i in range(1, len(arr)):
        cumsum += arr[i]
        result[i] = cumsum
    return result

@nb.jit
def count_nonzero(arr):
    arr = arr.ravel()
    n = 0
    for x in arr:
        if x != 0:
            n += 1
    return n

@nb.jit
def row_col_nonzero_nb(arr):
    n, m = arr.shape
    max_k = count_nonzero(arr)
    indices = np.empty((2, max_k), dtype=np.uint32)
    i_offset = np.zeros(n + 1, dtype=np.uint32)
    j_offset = np.zeros(m + 1, dtype=np.uint32)
    n, m = arr.shape
    k = 0
    for i in range(n):
        for j in range(m):
            if arr[i, j] != 0:
                indices[:, k] = i, j
                i_offset[i + 1] += 1
                j_offset[j + 1] += 1
                k += 1
    return indices, cumsum(i_offset), cumsum(j_offset)

def row_col_idx_nonzero_nb(arr):
    (ii, jj), jj_split, ii_split = row_col_nonzero_nb(arr)
    ii_ = np.argsort(jj)
    ii = ii[ii_]
    return np.split(ii, ii_split[1:-1]), np.split(jj, jj_split[1:-1])

根据www.example.com上的@hpaulj答案(row_col_idx_sparse_lil())和@knl答案(row_col_idx_sparse_coo()),与您的方法(下面的row_col_idx_sep())以及其他一些方法相比:scicomp.stackexchange.com ( row_col_idx_sparse_coo() ):
一个一个一个一个一个x一个一个二个一个x一个一个三个一个x一个一个x一个四个一个
对于使用生成的输入:

def gen_input(n, density=0.1, dtype=np.float32):
    arr = np.zeros(shape=(n, n), dtype=dtype)
    indices = tuple(np.random.randint(0, n, (2, int(n * n * density))).tolist())
    arr[indices] = 1.0
    return arr

我们会得到(您的test_matrix具有大约0.06的非零密度):

m = gen_input(4096, density=0.06)
%timeit row_col_idx_sep(m)
# 1 loop, best of 3: 767 ms per loop
%timeit row_col_idx_zip(m)
# 1 loop, best of 3: 660 ms per loop
%timeit row_col_idx_sparse_coo(m)
# 1 loop, best of 3: 205 ms per loop
%timeit row_col_idx_sparse_lil(m)
# 1 loop, best of 3: 498 ms per loop
%timeit row_col_idx_nonzero_nb(m)
# 10 loops, best of 3: 130 ms per loop

这表明它的速度接近最快的基于scipy.sparse的方法的两倍。

ufj5ltwl

ufj5ltwl2#

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

数据存在于整个数组nonzero中,只是没有分解为每行/列数组:

In [183]: np.nonzero(arr)                                                                
Out[183]: (array([1, 2, 2, 3, 3, 3]), array([0, 0, 1, 0, 1, 2]))
In [184]: np.argwhere(arr)                                                               
Out[184]: 
array([[1, 0],
       [2, 0],
       [2, 1],
       [3, 0],
       [3, 1],
       [3, 2]])

可以将array([1, 2, 2, 3, 3, 3])分解为基于另一个数组的子列表[1,2,3],[2,3],[3],[],但这可能需要一些时间来解决逻辑问题,而且不能保证它会比行/列迭代更快。
逻辑运算可以将布尔数组缩减为列或行,给出出现非零值的行或列,但也不是不规则的:

In [185]: arr!=0                                                                         
Out[185]: 
array([[False, False, False, False],
       [ True, False, False, False],
       [ True,  True, False, False],
       [ True,  True,  True, False]])
In [186]: (arr!=0).any(axis=0)                                                           
Out[186]: array([ True,  True,  True, False])
In [187]: np.nonzero((arr!=0).any(axis=0))                                               
Out[187]: (array([0, 1, 2]),)
In [188]: np.nonzero((arr!=0).any(axis=1))                                               
Out[188]: (array([1, 2, 3]),)
In [189]: arr                                                                            
Out[189]: 
array([[0, 0, 0, 0],
       [1, 0, 0, 0],
       [1, 1, 0, 0],
       [1, 1, 1, 0]])

scipy.sparselil确实会生成所需的数据:

In [190]: sparse                                                                         
Out[190]: <module 'scipy.sparse' from '/usr/local/lib/python3.6/dist-packages/scipy/sparse/__init__.py'>
In [191]: M = sparse.lil_matrix(arr)                                                     
In [192]: M                                                                              
Out[192]: 
<4x4 sparse matrix of type '<class 'numpy.longlong'>'
    with 6 stored elements in List of Lists format>
In [193]: M.rows                                                                         
Out[193]: array([list([]), list([0]), list([0, 1]), list([0, 1, 2])], dtype=object)
In [194]: M.T                                                                            
Out[194]: 
<4x4 sparse matrix of type '<class 'numpy.longlong'>'
    with 6 stored elements in List of Lists format>
In [195]: M.T.rows                                                                       
Out[195]: array([list([1, 2, 3]), list([2, 3]), list([3]), list([])], dtype=object)

但是定时可能并不比行或列迭代更好。

js4nwp54

js4nwp543#

例如,也可以使用np.bincount()来执行此操作

# Suppose X is a sparse matrix
In [56]: X.shape
Out[56]: (5830, 93761)

# Non-zero elements per row
In [57]: np.bincount(X.nonzero()[0])
Out[57]: array([316, 305,  98, ..., 284, 915, 721])

# Non-zero elements per column
In [58]: np.bincount(X.nonzero()[1])
Out[58]: array([14, 10,  1, ...,  1,  1,  1])

# Its also reasonably quick
In [59]: %timeit np.bincount(X.nonzero()[0])
24.1 ms ± 393 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

需要注意的一点是,矩阵中的行/列是否全为零:

In [60]: A = np.array([[0, 1, 0], [0, 0, 0]])
In [61]: np.bincount(A.nonzero()[0])
Out[61]: array([1]) # <-- Second row is all zeros, so its omitted here

In [62]: B = np.array([[0, 0, 0], [0, 0, 1]])
In [63]: np.bincount(B.nonzero()[0])
Out[63]: array([0, 1])

相关问题