Scipy.sparse.csr_matrix:如何获得十大价值和指数？

Question

我有一个很大的csr_matrix,我对前十个值及其每一行的指数感兴趣.但我没有找到一种操纵矩阵的好方法.

这是我目前的解决方案,主要思想是逐行处理它们:

row = csr_matrix.getrow(row_number).toarray()[0].ravel()
top_ten_indicies = row.argsort()[-10:]
top_ten_values = row[row.argsort()[-10:]]

通过这样做,csr_matrix没有充分利用其优点.它更像是一个强力解决方案.

Answer 1

csr在这种情况下,我没有看到格式的优点.当然,所有非零值都收集在一个.data数组中,其中包含相应的列索引.indices.但它们的长度各不相同.这意味着它们不能并行处理或使用numpy数组步长处理.

一种解决方案是将这些块填充到公共长度块中.这是什么.toarray().然后你可以用argsort(axis=1) or withargpartition` 找到最大值.

另一种方法是将它们分成行大小的块,然后处理每个块.这就是你正在做的事情.getrow.打破它们的另一种方法是转换为lil格式,并处理.data和.rows数组的子列表.

可能的第三种选择是使用该ufunc reduceat方法.这使您可以将ufunc reduction方法应用于数组的顺序块.有建立ufunc就像np.add利用这一点. argsort不是这样的功能.但是有一种方法可以ufunc从Python函数构造一个,并且比普通的Python迭代获得一些适度的速度.[我需要查看最近的一个SO问题来说明这一点.]

我将用一个更简单的函数来说明其中一些,总结行.

如果A2是csr矩阵.

A2.sum(axis=1)  # the fastest compile csr method
A2.A.sum(axis=1)  # same, but with a dense intermediary
[np.sum(l.data) for l in A2]  # iterate over the rows of A2
[np.sum(A2.getrow(i).data) for i in range(A2.shape[0])]  # iterate with index
[np.sum(l) for l in A2.tolil().data]  # sum the sublists of lil format
np.add.reduceat(A2.data, A2.indptr[:-1])  # with reduceat

A2.sum(axis=1)实现为矩阵乘法.这与排序问题无关,但仍然是查看求和问题的有趣方法.记住csr格式是为了有效的乘法而开发的.

对于我当前的样本矩阵(为另一个SO稀疏问题创建)

<8x47752 sparse matrix of type '<class 'numpy.float32'>'
     with 32 stored elements in Compressed Sparse Row format>

有些比较时间

In [694]: timeit np.add.reduceat(A2.data, A2.indptr[:-1])
100000 loops, best of 3: 7.41 µs per loop

In [695]: timeit A2.sum(axis=1)
10000 loops, best of 3: 71.6 µs per loop

In [696]: timeit [np.sum(l) for l in A2.tolil().data]
1000 loops, best of 3: 280 µs per loop

其他一切都是1ms或更长.

我建议专注于开发你的单行函数,例如:

def max_n(row_data, row_indices, n):
    i = row_data.argsort()[-n:]
    # i = row_data.argpartition(-n)[-n:]
    top_values = row_data[i]
    top_indices = row_indices[i]  # do the sparse indices matter?
    return top_values, top_indices, i

然后看看如何适应这些迭代方法之一. tolil()看起来最有希望.

我没有解决如何收集这些结果的问题.它们应该是列表列表,10列的数组,每行10个值的另一个稀疏矩阵等等吗？

排序大的稀疏和保存的顶部K值和列索引的每一行 - 几年前的类似问题,但没有答案.

scipy稀疏矩阵中每行或每列的Argmax - 最近寻找argmax行的问题csr.我讨论了一些相同的问题.

如何在numpy中加速循环？ - 如何使用np.frompyfunc创建的示例ufunc.我不知道结果函数是否具有该.reduceat方法.

增加稀疏矩阵中前k个元素的值 - 获取csr的前k个元素(不是按行).案例argpartition.

行总和实现np.frompyfunc:

In [741]: def foo(a,b):
    return a+b  
In [742]: vfoo=np.frompyfunc(foo,2,1)
In [743]: timeit vfoo.reduceat(A2.data,A2.indptr[:-1],dtype=object).astype(float)
10000 loops, best of 3: 26.2 µs per loop

这是可观的速度.但是我想不出一种编写二进制函数的方法(需要2个参数)来实现argsort通过简化.所以这可能是这个问题的一个障碍.

Answer 2

只是为了回答最初的问题（对于像我这样发现这个问题寻找 copy-pasta 的人），这里有一个使用多处理的解决方案，该解决方案基于@hpaulj 的转换为lil_matrix，并遍历行的建议

from multiprocessing import Pool

def _top_k(args):
    """
    Helper function to process a single row of top_k
    """
    data, row = args
    data, row = zip(*sorted(zip(data, row), reverse=True)[:k])
    return data, row

def top_k(m, k):
    """
    Keep only the top k elements of each row in a csr_matrix
    """
    ml = m.tolil()
    with Pool() as p:
        ms = p.map(_top_k, zip(ml.data, ml.rows))
    ml.data, ml.rows = zip(*ms)
    return ml.tocsr()

Answer 3

需要迭代行并分别获取每行的顶部索引。但是这个循环可以被jited（和并行化）以获得极快的功能。

@nb.njit(cache=True)
def row_topk_csr(data, indices, indptr, K):
    m = indptr.shape[0] - 1
    max_indices = np.zeros((m, K), dtype=indices.dtype)
    max_values = np.zeros((m, K), dtype=data.dtype)

    for i in nb.prange(m):
        top_inds = np.argsort(data[indptr[i] : indptr[i + 1]])[::-1][:K]
        max_indices[i] = indices[indptr[i] : indptr[i + 1]][top_inds]
        max_values[i] = data[indptr[i] : indptr[i + 1]][top_inds]

    return max_indices, max_values

像这样称呼它：

top_pred_indices, _ = row_topk_csr(csr_mat.data, csr_mat.indices, csr_mat.indptr, K)

我需要频繁执行此操作，并且此函数对我来说足够快，在 1mil x 400k 稀疏矩阵上执行时间小于 1 秒。

HTH。