大型 scipy 稀疏矩阵上的快速列访问

Question

我正在使用 scipy 的 csc 稀疏矩阵，目前代码中的主要瓶颈是类似于以下的一行

for i in range(multiply_cols.shape[0]):
   F = F - factor*values[i]*mat.getcol(multiply_cols[i])

我正在使用的矩阵非常大，通常大于10**6x10**6，我不想将它们转换为密集矩阵。事实上，我有一个限制，即始终使用csc 格式的矩阵。我的尝试表明，转换为 coo_matrix 或 lil_matrix 也没有回报。

这是我使用csc、csr 和coo 的初步尝试：

n=1000
sA = csc_matrix(np.random.rand(n,n))
F = np.random.rand(n,1)
multiply_cols = np.unique(np.random.randint(0,int(0.6*n),size=n))
values = np.random.rand(multiply_cols.shape[0])

def foo1(mat,F,values,multiply_cols):
    factor = 0.75
    for i in range(multiply_cols.shape[0]):
        F = F - factor*values[i]*mat.getcol(multiply_cols[i])

def foo2(mat,F,values,multiply_cols):
    factor = 0.75
    mat = mat.tocsr()
    for i in range(multiply_cols.shape[0]):
        F = F - factor*values[i]*mat.getcol(multiply_cols[i])

def foo3(mat,F,values,multiply_cols):
    factor = 0.75
    mat = mat.tocoo()
    for i in range(multiply_cols.shape[0]):
        F = F - factor*values[i]*mat.getcol(multiply_cols[i])

def foo4(mat,F,values,multiply_cols):
    factor = 0.75
    mat = mat.tolil()
    for i in range(multiply_cols.shape[0]):
        F = F - factor*values[i]*mat.getcol(multiply_cols[i])

我得到它们的时间：

In [41]: %timeit foo1(sA,F,values,multiply_cols)
10 loops, best of 3: 133 ms per loop

In [42]: %timeit foo2(sA,F,values,multiply_cols)
1 loop, best of 3: 999 ms per loop

In [43]: %timeit foo3(sA,F,values,multiply_cols)
1 loop, best of 3: 6.38 s per loop

In [44]: %timeit foo4(sA,F,values,multiply_cols)
1 loop, best of 3: 45.1 s per loop

所以coo_matrix 和lil_matrix 肯定不是一个好的选择。有谁知道这样做的更快方法。检索底层indptr、indices 和data 是否有自定义cython 解决方案是一个不错的选择？

Answer 1

我发现

Sparse matrix slicing using list of int

sparse 矩阵的列（或行）索引本质上是一个矩阵乘法任务 - 构造一个 1 和 0 正确混合的稀疏矩阵，然后相乘。行（和列）总和也是通过乘法完成的。

这个函数实现了这个想法。 M 是一个 1 列稀疏矩阵，values 在 multiply_cols 插槽中：

def wghtsum(sA, values, multiply_cols):
   cols = np.zeros_like(multiply_cols)
   M=sparse.csc_matrix((values,(multiply_cols,cols)),shape=(sA.shape[1],1))
   return (sA*M).A

测试：

In [794]: F1=wghtsum(sA,values,multiply_cols)
In [800]: F2=(sA[:,multiply_cols]*values)[:,None]  # Divaker's
In [802]: np.allclose(F1,F2)
Out[802]: True

与@Divakar's 解决方案相比，它节省了适度的时间：

In [803]: timeit F2=(sA[:,multiply_cols]*values)[:,None]
100 loops, best of 3: 18.3 ms per loop
In [804]: timeit F1=wghtsum(sA,values,multiply_cols)
100 loops, best of 3: 6.57 ms per loop

=======

sA 创建时是密集的 - 它是密集随机数组的稀疏再现。 sparse.rand 可用于创建具有已定义稀疏度的稀疏随机矩阵。

在测试您的foo1 时，我遇到了getcol 的问题：

In [818]: sA.getcol(multiply_cols[0])
...
TypeError: an integer is required
In [819]: sA.getcol(multiply_cols[0].item())
Out[819]: 
<1000x1 sparse matrix of type '<class 'numpy.float64'>'
    with 1000 stored elements in Compressed Sparse Column format>
In [822]: sA[:,multiply_cols[0]]
Out[822]: 
<1000x1 sparse matrix of type '<class 'numpy.float64'>'
    with 1000 stored elements in Compressed Sparse Column format>

我怀疑这是由 scipy 版本差异引起的。

In [821]: scipy.__version__
Out[821]: '0.17.0'

这个问题在 0.18 中确实消失了；但我找不到相关的问题/拉取请求。

Answer 2

你可以使用向量化方法，使用从稀疏矩阵中切出的列与values 的矩阵乘法，就像这样 -

F -= (mat[:,multiply_cols]*values*factor)[:,None]

基准测试

似乎foo1 是问题中列出的最快的。因此，让我们针对该方法对提议的方法进行计时。

函数定义-

def foo1(mat,F,values,multiply_cols):
    factor = 0.75
    outF = F.copy()
    for i in range(multiply_cols.shape[0]):
        outF -= factor*values[i]*mat.getcol(multiply_cols[i])
    return outF

def foo_vectorized(mat,F,values,multiply_cols):
    factor = 0.75 
    return F - (mat[:,multiply_cols]*values*factor)[:,None]

在更大的稀疏集合上进行计时和验证 -

In [242]: # Setup inputs
     ...: n = 3000
     ...: mat = csc_matrix(np.random.randint(0,3,(n,n))) #Sparseness with  0s
     ...: F = np.random.rand(n,1)
     ...: multiply_cols = np.unique(np.random.randint(0,int(0.6*n),size=n))
     ...: values = np.random.rand(multiply_cols.shape[0])
     ...: 

In [243]: out1 = foo1(mat,F,values,multiply_cols)

In [244]: out2 = foo_vectorized(mat,F,values,multiply_cols)

In [245]: np.allclose(out1, out2)
Out[245]: True

In [246]: %timeit foo1(mat,F,values,multiply_cols)
1 loops, best of 3: 641 ms per loop

In [247]: %timeit foo_vectorized(mat,F,values,multiply_cols)
10 loops, best of 3: 40.3 ms per loop

In [248]: 641/40.3
Out[248]: 15.905707196029779

我们有一个 15x+ 加速！