什么是极度稀疏矩阵的最快乘法方法？

Question

什么是极度稀疏矩阵的最快乘法方法？

3

我有一个非常稀疏的结构化矩阵。我的矩阵每列恰好有一个非零条目。但它很大（10k * 1M），并以以下形式给出（使用随机值作为示例）

rows = np.random.randint(0, 10000, 1000000)
values = np.random.randint(0,10,1000000)

这里的rows指的是每列中非零元素所在的行号。我希望通过S来进行快速矩阵乘法，目前我的做法是将其转换为稀疏矩阵(S)，然后使用S.dot(X)进行和矩阵X(可以是稀疏或密集矩阵)的乘法。

S=scipy.sparse.csr_matrix( (values, (rows, scipy.arange(1000000))), shape = (10000,1000000))

对于大小为1M * 2500且nnz（X）= 8M的X，创建S需要178ms，应用它需要255ms。所以我的问题是，如果给定描述的S，最好的做SX的方法是什么（其中X可以是稀疏或密集的）。由于创建S本身非常耗时，我考虑使用一些特定的方法。我尝试使用循环创建了一些东西，但效果不佳。
简单的循环过程如下：

SX = np.zeros((rows.size,X.shape[1]))
 for i in range(X.shape[0]):
     SX[rows[i],:]+=values[i]*X[i,:]
 return SX

我们能使这个过程更有效吗？
非常感谢任何建议。谢谢

- user1131274

考虑到您的矩阵有多大，这些时间非常不错。我认为除了使用完全不同的框架之外，您无法再挤出更多性能了。 - rayryeng

@rayryeng 我觉得有点慢，因为例如在我给出的例子中，做X^TX需要450毫秒，这与做SA所需的时间相等。考虑到当X是m*n时，做SA是O(mn)，而做X^TX是O(mn^2)，我发现做SA不够快。此外，我知道我没有考虑X的稀疏性，但对于密集的X，我有类似的数字。谢谢。 - user1131274

每列只有一个值，您可以通过使用“行”索引“X”，然后使用密集的“点”和“值”来更好地处理。 - hpaulj

3个回答

1

受稀疏矩阵行求和的最快方法这篇文章的启发，我发现最好的方法是编写循环并将其转换为numba。以下是：

`

@njit
def sparse_mul(SX,row,col,data,values,row_map):
    N = len(data)
    for idx in range(N):
        SX[row_map[row[idx]],col[idx]]+=data[idx]*values[row[idx]]
    return SX
X_coo=X.tocoo()
s=row_map.max()+1
SX = np.zeros((s,X.shape[1]))
sparse_mul(SX,X_coo.row,X_coo.col,X_coo.data,values,row_map)`

这里的row_map是问题中的行。对于大小为(1M* 1K)、1%稀疏度且s=10K的X，使用此方法比从row_map形成稀疏矩阵并进行S.dot(A)要快两倍。

- user1131274

0

据我回忆，Knuth的TAOP（《计算机程序设计艺术》）提到了将稀疏矩阵表示为非零值的链表（适用于您的应用程序）。也许可以尝试这种方法？然后通过遍历链表而不是每个维度上的整个数组来实现。

- Mark Diaz

1

那可能不会有所帮助。scipy 中的稀疏矩阵表示和稀疏矩阵乘法已经非常优化，因此使用未经优化的稀疏矩阵表示进行乘法运算可能比当前基准测试所需的时间更长。（顺便说一下，我没有投反对票）。 - rayryeng

不必担心任何负评 - 如果回答值得的话。任何成功的程序员都不会认为自我比现实更重要。 - Mark Diaz

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Divakar · Accepted Answer

方案 #1

考虑到第一个输入中每一列都只有一个条目，我们可以使用np.bincount函数处理输入的rows、values和X，从而避免创建稀疏矩阵S。

def sparse_matrix_mult(rows, values, X):
    nrows = rows.max()+1
    ncols = X.shape[1]
    nelem = nrows * ncols

    ids = rows + nrows*np.arange(ncols)[:,None]
    sums = np.bincount(ids.ravel(), (X.T*values).ravel(), minlength=nelem)
    out = sums.reshape(ncols,-1).T
    return out

样例运行 -

In [746]: import numpy as np
     ...: from scipy.sparse import csr_matrix
     ...: import scipy as sp
     ...: 

In [747]: np.random.seed(1234)
     ...: m,n = 3,4
     ...: rows = np.random.randint(0, m, n)
     ...: values = np.random.randint(2,10,n)
     ...: X = np.random.randint(2, 10, (n,5))
     ...: 

In [748]: S = csr_matrix( (values, (rows, sp.arange(n))), shape = (m,n))

In [749]: S.dot(X)
Out[749]: 
array([[42, 27, 45, 78, 87],
       [24, 18, 18, 12, 24],
       [18,  6,  8, 16, 10]])

In [750]: sparse_matrix_mult(rows, values, X)
Out[750]: 
array([[ 42.,  27.,  45.,  78.,  87.],
       [ 24.,  18.,  18.,  12.,  24.],
       [ 18.,   6.,   8.,  16.,  10.]])

方法二

使用np.add.reduceat替换np.bincount-

def sparse_matrix_mult_v2(rows, values, X):
    nrows = rows.max()+1
    ncols = X.shape[1]

    scaled_ar = X*values[:,None]
    sidx = rows.argsort()
    rows_s = rows[sidx]
    cut_idx = np.concatenate(([0],np.flatnonzero(rows_s[1:] != rows_s[:-1])+1))
    sums = np.add.reduceat(scaled_ar[sidx],cut_idx,axis=0)

    out = np.empty((nrows, ncols),dtype=sums.dtype)
    row_idx = rows_s[cut_idx]
    out[row_idx] = sums
    return out

运行测试

我无法按照问题中提到的那些大小运行它，因为那对我来说太大了。因此，在缩小的数据集上，这是我得到的结果 -

In [149]: m,n = 1000, 100000
     ...: rows = np.random.randint(0, m, n)
     ...: values = np.random.randint(2,10,n)
     ...: X = np.random.randint(2, 10, (n,2500))
     ...: 

In [150]: S = csr_matrix( (values, (rows, sp.arange(n))), shape = (m,n))

In [151]: %timeit csr_matrix( (values, (rows, sp.arange(n))), shape = (m,n))
100 loops, best of 3: 16.1 ms per loop

In [152]: %timeit S.dot(X)
1 loop, best of 3: 193 ms per loop

In [153]: %timeit sparse_matrix_mult(rows, values, X)
1 loop, best of 3: 4.4 s per loop

In [154]: %timeit sparse_matrix_mult_v2(rows, values, X)
1 loop, best of 3: 2.81 s per loop

因此，提出的方法似乎在性能方面不如numpy.dot厉害，但它们在内存效率方面应该很好。

对于稀疏的X

对于稀疏的X，我们需要进行一些修改，具体列在下面的修改方法中 -

from scipy.sparse import find
def sparse_matrix_mult_sparseX(rows, values, Xs): # Xs is sparse    
    nrows = rows.max()+1
    ncols = Xs.shape[1]
    nelem = nrows * ncols

    scaled_vals = Xs.multiply(values[:,None])
    r,c,v = find(scaled_vals)
    ids = rows[r] + c*nrows
    sums = np.bincount(ids, v, minlength=nelem)
    out = sums.reshape(ncols,-1).T
    return out