高效计算每个非零元素为 1 的稀疏数组的列总和

Question

我有一堆 SciPy compressed sparse row (CSR) 格式的数据。当然，大多数元素都是零，而且我进一步知道所有非零元素的值为 1。我想计算矩阵行的不同子集的总和。目前我正在做以下事情：

import numpy as np
import scipy as sp
import scipy.sparse

# create some data with sparsely distributed ones
data = np.random.choice((0, 1), size=(1000, 2000), p=(0.95, 0.05))
data = sp.sparse.csr_matrix(data, dtype='int8')

# generate column-wise sums over random subsets of rows
nrand = 1000
for k in range(nrand):
    inds = np.random.choice(data.shape[0], size=100, replace=False)

    # 60% of time is spent here
    extracted_rows = data[inds]

    # 20% of time is spent here
    row_sum = extracted_rows.sum(axis=0)

最后几行是较大计算管道中的瓶颈。正如我在代码中注释的那样，60% 的时间花在从随机索引中切分数据上，20% 的时间花在计算实际总和上。

在我看来，我应该能够使用我对数组中数据的了解（即，稀疏矩阵中的任何非零值都将为 1；不存在其他值）来更有效地计算这些总和.不幸的是，我无法弄清楚如何。也许只处理 data.indices？我尝试过其他稀疏结构（例如 CSC 矩阵），以及先转换为密集数组，但这些方法都比这种 CSR 矩阵方法慢。

Answer 1

众所周知，稀疏矩阵的索引相对较慢。关于通过直接访问数据属性来解决这个问题有很多问题。

但首先是一些时间安排。使用 data 和 ind 如您所见，我得到

In [23]: datad=data.A  # times at 3.76 ms per loop

In [24]: timeit row_sumd=datad[inds].sum(axis=0)
1000 loops, best of 3: 529 µs per loop

In [25]: timeit row_sum=data[inds].sum(axis=0)
1000 loops, best of 3: 890 µs per loop

In [26]: timeit d=datad[inds]
10000 loops, best of 3: 55.9 µs per loop

In [27]: timeit d=data[inds]
1000 loops, best of 3: 617 µs per loop

稀疏版本比密集版本慢，但相差不大。稀疏索引要慢得多，但它的总和要快一些。

稀疏求和是用矩阵乘积完成的

def sparse.spmatrix.sum
     ....
    return np.asmatrix(np.ones((1, m), dtype=res_dtype)) * self

这表明更快的方法 - 将 inds 转换为适当的 1 数组并相乘。

In [49]: %%timeit
   ....: b=np.zeros((1,data.shape[0]),'int8')
   ....: b[:,inds]=1
   ....: rowmul=b*data
   ....: 
1000 loops, best of 3: 587 µs per loop

这使得稀疏运算与等效的密集运算一样快。（但转换为密集要慢得多）

==================

最后一次测试缺少稀疏 sum 中存在的 np.asmatrix。但时代相似，结果相同

In [232]: timeit b=np.zeros((1,data.shape[0]),'int8'); b[:,inds]=1; x1=np.asmatrix(b)*data
1000 loops, best of 3: 661 µs per loop

In [233]: timeit b=np.zeros((1,data.shape[0]),'int8'); b[:,inds]=1; x2=b*data
1000 loops, best of 3: 605 µs per loop

一个生成矩阵，另一个生成数组。但是两者都在做矩阵乘积，B 的第 2 维对 data 的第 1 维。尽管 b 是一个数组，但任务实际上是委托给 data 及其矩阵乘积 - 以一种不太透明的方式。

In [234]: x1
Out[234]: matrix([[9, 9, 5, ..., 9, 5, 3]], dtype=int8)

In [235]: x2
Out[235]: array([[9, 9, 5, ..., 9, 5, 3]], dtype=int8)

b*data.A 是元素乘法并引发错误； np.dot(b,data.A) 有效，但速度较慢。

较新的 numpy/python 有一个 matmul 运算符。我看到了相同的时间模式：

In [280]: timeit b@dataA           # dense product
100 loops, best of 3: 2.64 ms per loop

In [281]: timeit b@data.A           # slower due to `.A` conversion
100 loops, best of 3: 6.44 ms per loop

In [282]: timeit b@data             # sparse product
1000 loops, best of 3: 571 µs per loop

np.dot 也可以将操作委托给 sparse，但您必须小心。我刚刚用 np.dot(csr_matrix(b),data.A).

挂了我的机器

Answer 2

这是将 data 转换为密集数组并使用 -

以矢量化方式获取所有这些 inds 之后的矢量化方法

# Number of selections as a parameter
n = 100

# Get inds across all iterations in a vectorized manner as a 2D array.
inds2D = np.random.rand(nrand,data.shape[0]).argpartition(n)[:,:n]

# Index into data with those 2D array indices. Then, convert to dense NumPy array, 
# reshape and sum reduce to get the final output 
out = np.array(data.todense())[inds2D.ravel()].reshape(nrand,n,-1).sum(1)

运行时测试 -

1）函数定义：

def org_app(nrand,n):
    out = np.zeros((nrand,data.shape[1]),dtype=int)
    for k in range(nrand):
        inds = np.random.choice(data.shape[0], size=n, replace=False)        
        extracted_rows = data[inds]    
        out[k] = extracted_rows.sum(axis=0)    
    return out


def vectorized_app(nrand,n):
    inds2D = np.random.rand(nrand,data.shape[0]).argpartition(n)[:,:n]
    return np.array(data.todense())[inds2D.ravel()].reshape(nrand,n,-1).sum(1)

时间安排：

In [205]: # create some data with sparsely distributed ones
     ...: data = np.random.choice((0, 1), size=(1000, 2000), p=(0.95, 0.05))
     ...: data = sp.sparse.csr_matrix(data, dtype='int8')
     ...: 
     ...: # generate column-wise sums over random subsets of rows
     ...: nrand = 1000
     ...: n = 100
     ...: 

In [206]: %timeit org_app(nrand,n)
1 loops, best of 3: 1.38 s per loop

In [207]: %timeit vectorized_app(nrand,n)
1 loops, best of 3: 826 ms per loop

高效计算每个非零元素为 1 的稀疏数组的列总和

Efficiently compute columnwise sum of sparse array where every non-zero element is 1

python

arrays

numpy

scipy

sparse-matrix