首页 文章

SciPy稀疏矩阵(COO,CSR):清除行

提问于
浏览
1

为了创建scipy sparse matrix,我有一个数组或行和列索引 IJ 以及数据数组 V . 我使用它们在COO format中构造一个矩阵,然后将其转换为CSR

matrix = sparse.coo_matrix((V, (I, J)), shape=(n, n))
matrix = matrix.tocsr()

我有一组行索引,对角线上唯一的条目应该是 1.0 . 到目前为止,我通过 I ,查找需要擦除的所有索引,并执行以下操作:

def find(lst, a):
    # From <http://stackoverflow.com/a/16685428/353337>
    return [i for i, x in enumerate(lst) if x in a]

# wipe_rows = [1, 55, 32, ...]  # something something

indices = find(I, wipe_rows)  # takes too long
I = numpy.delete(I, indices).tolist()
J = numpy.delete(J, indices).tolist()
V = numpy.delete(V, indices).tolist()

# Add entry 1.0 to the diagonal for each wipe row
I.extend(wipe_rows)
J.extend(wipe_rows)
V.extend(numpy.ones(len(wipe_rows)))

# construct matrix via coo

这没关系,但 find 往往需要一段时间 .

关于如何提高速度的任何提示? (也许以COO或CSR格式擦除行是一个更好的主意 . )

2 回答

  • 2

    如果您打算一次清除多行,请执行此操作

    def _wipe_rows_csr(matrix, rows):
        assert isinstance(matrix, sparse.csr_matrix)
    
        # delete rows
        for i in rows:
            matrix.data[matrix.indptr[i]:matrix.indptr[i+1]] = 0.0
    
        # Set the diagonal
        d = matrix.diagonal()
        d[rows] = 1.0
        matrix.setdiag(d)
    
        return
    

    是迄今为止最快的方法 . 它并没有真正删除线条,而是将所有条目设置为零,然后用对角线填充 .

    如果实际要删除条目,则必须进行一些数组操作 . 这可能是非常昂贵的,但如果速度没有问题:这

    def _wipe_row_csr(A, i):
        '''Wipes a row of a matrix in CSR format and puts 1.0 on the diagonal.
        '''
        assert isinstance(A, sparse.csr_matrix)
    
        n = A.indptr[i+1] - A.indptr[i]
    
        assert n > 0
    
        A.data[A.indptr[i]+1:-n+1] = A.data[A.indptr[i+1]:]
        A.data[A.indptr[i]] = 1.0
        A.data = A.data[:-n+1]
    
        A.indices[A.indptr[i]+1:-n+1] = A.indices[A.indptr[i+1]:]
        A.indices[A.indptr[i]] = i
        A.indices = A.indices[:-n+1]
    
        A.indptr[i+1:] -= n-1
    
        return
    

    通过对角线上的条目 1.0 替换矩阵 matrix 的给定行 i .

  • 0

    np.in1d 应该是找到 indices 的更快捷方式:

    In [322]: I   # from a np.arange(12).reshape(4,3) matrix
    Out[322]: array([0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3], dtype=int32)
    
    In [323]: indices=[i for i, x in enumerate(I) if x in [1,2]]
    
    In [324]: indices
    Out[324]: [2, 3, 4, 5, 6, 7]
    
    In [325]: ind1=np.in1d(I,[1,2])
    
    In [326]: ind1
    Out[326]: 
    array([False, False,  True,  True,  True,  True,  True,  True, False,
           False, False], dtype=bool)
    
    In [327]: np.where(ind1)   # same as indices
    Out[327]: (array([2, 3, 4, 5, 6, 7], dtype=int32),)
    
    In [328]: I[~ind1]  # same as the delete
    Out[328]: array([0, 0, 3, 3, 3], dtype=int32)
    

    像这样直接操纵 coo 输入通常是一种好方法 . 但另一个是利用 csr 数学能力 . 您应该能够构造一个对齐矩阵,将正确的行归零,然后将其重新添加 .

    这就是我的想法:

    In [357]: A=np.arange(16).reshape(4,4)
    In [358]: M=sparse.coo_matrix(A)
    In [359]: M.A
    Out[359]: 
    array([[ 0,  1,  2,  3],
           [ 4,  5,  6,  7],
           [ 8,  9, 10, 11],
           [12, 13, 14, 15]])
    
    In [360]: d1=sparse.diags([(1,0,0,1)],[0],(4,4))
    In [361]: d2=sparse.diags([(0,1,1,0)],[0],(4,4))
    
    In [362]: (d1*M+d2).A
    Out[362]: 
    array([[  0.,   1.,   2.,   3.],
           [  0.,   1.,   0.,   0.],
           [  0.,   0.,   1.,   0.],
           [ 12.,  13.,  14.,  15.]])
    
    In [376]: x=np.ones((4,),bool);x[[1,2]]=False
    In [378]: d1=sparse.diags([x],[0],(4,4),dtype=int)
    In [379]: d2=sparse.diags([~x],[0],(4,4),dtype=int)
    

    使用 lil 格式执行此操作看起来很简单:

    In [593]: Ml=M.tolil()
    In [594]: Ml.data[wipe]=[[1]]*len(wipe)
    In [595]: Ml.rows[wipe]=[[i] for i in wipe]
    
    In [596]: Ml.A
    Out[596]: 
    array([[ 0,  1,  2,  3],
           [ 0,  1,  0,  0],
           [ 0,  0,  1,  0],
           [12, 13, 14, 15]], dtype=int32)
    

    这是你用 csr 格式做的事情,但是用适当的[1]和[i]列表替换每个行列表很容易 . 但转换时间( tolil 等)可能会影响运行时间 .

相关问题