python - 与极稀疏矩阵相乘的最快方法是什么？

Question

我有一个极其稀疏的结构化矩阵。我的矩阵每列只有一个非零条目。但是它很大(10k * 1M)并以以下形式给出(例如uisng随机值)

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 需要 178 毫秒，应用它需要 255 毫秒。所以我的问题是，鉴于我的 S 如上所述，做 SX(其中 X 可能稀疏或密集)的最佳方法是什么。由于创建 S 本身非常耗时，所以我在考虑一些临时的东西。我确实尝试使用循环创建一些东西，但它甚至没有关闭。
简单的循环过程看起来像这样

SX = np.zeros((rows.size,X.shape[1]))
对于范围内的我(X.shape[0]):
SX[行[i],:]+=值[i]*X[i,:]
返回 SX
我们能让它变得高效吗？

非常感谢任何建议。谢谢

Answer 1

方法 #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

sample 运行-

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.]])

方法 #2

使用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