Scipy 稀疏矩阵循环永远持续 - 需要提高效率

Question

用稀疏矩阵编写这个循环的最有效的时间和内存方式是什么（目前使用 csc_matrix）

for j in range(0, reducedsize):
    xs = sum(X[:, j])
    X[:, j] = X[:, j] / xs.data[0]

例子：

缩小尺寸 (int) - 2500
X (csc_matrix) - 908x2500

循环确实会迭代，但与仅使用 numpy 相比需要很长时间。

Answer 1

In [388]: from scipy import sparse&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp;制作样本矩阵：In [390]: M = sparse.random(10,8,.2, 'csc')&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;我是一个矩阵：In [393]: M.sum(axis=0)&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;Out[393]:&nbsp;matrix([[1.95018736, 0.90924629, 1.93427113, 2.38816133, 1.08713479,&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;0.&nbsp; &nbsp; &nbsp; &nbsp; , 2.45435481, 0.&nbsp; &nbsp; &nbsp; &nbsp; ]])那些 0 在划分时会产生警告 -nan在结果中：In [394]: M/_&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;/usr/local/lib/python3.6/dist-packages/scipy/sparse/base.py:599: RuntimeWarning: invalid value encountered in true_divide&nbsp; return np.true_divide(self.todense(), other)Out[394]:&nbsp;matrix([[0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.27079623,&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; nan, 0.13752665,&nbsp; &nbsp; &nbsp; &nbsp; nan],&nbsp; &nbsp; &nbsp; &nbsp; [0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; ,&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; nan, 0.32825122,&nbsp; &nbsp; &nbsp; &nbsp; nan],&nbsp; &nbsp; &nbsp; &nbsp; [0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; ,&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; nan, 0.&nbsp; &nbsp; &nbsp; &nbsp; ,&nbsp; &nbsp; &nbsp; &nbsp; nan],&nbsp;...&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; nan, 0.&nbsp; &nbsp; &nbsp; &nbsp; ,&nbsp; &nbsp; &nbsp; &nbsp; nan]])0 也会给您的方法带来问题：In [395]: for i in range(8):&nbsp;&nbsp; &nbsp; &nbsp;...:&nbsp; &nbsp; &nbsp;xs = sum(M[:,i])&nbsp;&nbsp; &nbsp; &nbsp;...:&nbsp; &nbsp; &nbsp;M[:,i] = M[:,i]/xs.data[0]&nbsp;&nbsp; &nbsp; &nbsp;...:&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;---------------------------------------------------------------------------IndexError&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; Traceback (most recent call last)<ipython-input-395-0195298ead19> in <module>&nbsp; &nbsp; &nbsp; 1 for i in range(8):&nbsp; &nbsp; &nbsp; 2&nbsp; &nbsp; &nbsp;xs = sum(M[:,i])----> 3&nbsp; &nbsp; &nbsp;M[:,i] = M[:,i]/xs.data[0]&nbsp; &nbsp; &nbsp; 4&nbsp;IndexError: index 0 is out of bounds for axis 0 with size 0但是，如果我们比较没有 0 和的列，则值匹配：In [401]: Out[394][:,:5]&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp;Out[401]:&nbsp;matrix([[0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.27079623],&nbsp; &nbsp; &nbsp; &nbsp; [0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; ],&nbsp; &nbsp; &nbsp; &nbsp; [0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; ],&nbsp; &nbsp; &nbsp; &nbsp; [0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; ],&nbsp; &nbsp; &nbsp; &nbsp; [0.49648886, 0.25626608, 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.19162678, 0.72920377],&nbsp; &nbsp; &nbsp; &nbsp; [0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.30200765, 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; ],&nbsp; &nbsp; &nbsp; &nbsp; [0.50351114, 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.30445113, 0.41129367, 0.&nbsp; &nbsp; &nbsp; &nbsp; ],&nbsp; &nbsp; &nbsp; &nbsp; [0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.74373392, 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; ],&nbsp; &nbsp; &nbsp; &nbsp; [0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.39354122, 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; ],&nbsp; &nbsp; &nbsp; &nbsp; [0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.39707955, 0.&nbsp; &nbsp; &nbsp; &nbsp; ]])In [402]: M.A[:,:5]&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;Out[402]:&nbsp;array([[0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.27079623],&nbsp; &nbsp; &nbsp; &nbsp;[0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; ],&nbsp; &nbsp; &nbsp; &nbsp;[0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; ],&nbsp; &nbsp; &nbsp; &nbsp;[0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; ],&nbsp; &nbsp; &nbsp; &nbsp;[0.49648886, 0.25626608, 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.19162678, 0.72920377],&nbsp; &nbsp; &nbsp; &nbsp;[0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.30200765, 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; ],&nbsp; &nbsp; &nbsp; &nbsp;[0.50351114, 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.30445113, 0.41129367, 0.&nbsp; &nbsp; &nbsp; &nbsp; ],&nbsp; &nbsp; &nbsp; &nbsp;[0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.74373392, 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; ],&nbsp; &nbsp; &nbsp; &nbsp;[0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.39354122, 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; ],&nbsp; &nbsp; &nbsp; &nbsp;[0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.&nbsp; &nbsp; &nbsp; &nbsp; , 0.39707955, 0.&nbsp; &nbsp; &nbsp; &nbsp; ]])回到 [394] 我应该首先将矩阵和转换为稀疏，所以结果也将是稀疏的。稀疏没有逐元素除法，所以我必须先求密集矩阵的逆。0 仍然很麻烦。In [409]: M.multiply(sparse.csr_matrix(1/Out[393]))&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;...Out[409]:&nbsp;<10x8 sparse matrix of type '<class 'numpy.float64'>'&nbsp; &nbsp; with 16 stored elements in Compressed Sparse Column format>

Answer 2

如果您想在没有任何内存开销的情况下执行此操作（就地）始终考虑数据的实际存储方式。csc 矩阵的一个小例子。shape=(5,5)X=sparse.random(shape[0], shape[1], density=0.5, format='csc')print(X.todense())[[0.12146814 0.&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;0.&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;0.04075121 0.28749552]&nbsp;[0.&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;0.92208639 0.&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;0.44279661 0.&nbsp; &nbsp; &nbsp; &nbsp; ]&nbsp;[0.63509196 0.42334964 0.&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;0.&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;0.99160443]&nbsp;[0.&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;0.&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;0.25941113 0.44669367 0.00389409]&nbsp;[0.&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;0.&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;0.&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;0.&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;0.83226886]]i=0 #first columnprint(X.data[X.indptr[i]:X.indptr[i+1]])[0.12146814 0.63509196]一个 Numpy 解决方案所以我们在这里唯一要做的就是逐列修改非零条目。这可以使用部分矢量化的 numpy 解决方案轻松完成。data只是包含所有非零值的数组，indptr存储每列开始和结束的信息。def Numpy_csc_norm(data,indptr):&nbsp; &nbsp; for i in range(indptr.shape[0]-1):&nbsp; &nbsp; &nbsp; &nbsp; xs = np.sum(data[indptr[i]:indptr[i+1]])&nbsp; &nbsp; &nbsp; &nbsp; #Modify the view in place&nbsp; &nbsp; &nbsp; &nbsp; data[indptr[i]:indptr[i+1]]/=xs关于性能，这个就地解决方案已经不错了。如果你想进一步提高性能，你可以使用 Cython/Numba/ 或其他一些可以或多或少容易地用 Python 包装的编译代码。一个 Numba 解决方案@nb.njit(fastmath=True,error_model="numpy",parallel=True)def Numba_csc_norm(data,indptr):&nbsp; &nbsp; for i in nb.prange(indptr.shape[0]-1):&nbsp; &nbsp; &nbsp; &nbsp; acc=0&nbsp; &nbsp; &nbsp; &nbsp; for j in range(indptr[i],indptr[i+1]):&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; acc+=data[j]&nbsp; &nbsp; &nbsp; &nbsp; for j in range(indptr[i],indptr[i+1]):&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; data[j]/=acc表现#Create a not to small example matrixshape=(50_000,10_000)X=sparse.random(shape[0], shape[1], density=0.001, format='csc')#Not in-place from hpauljdef hpaulj(X):&nbsp; &nbsp; acc=X.sum(axis=0)&nbsp; &nbsp; return X.multiply(sparse.csr_matrix(1./acc))%timeit X2=hpaulj(X)#6.54 ms ± 67.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)#All 2 variants are in-place,&nbsp;#but this shouldn't have a influence on the timings%timeit Numpy_csc_norm(X.data,X.indptr)#79.2 ms ± 914 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)#parallel=False -> faster on tiny matrices%timeit Numba_csc_norm(X.data,X.indptr)#626 µs ± 30.6 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)#parallel=True -> faster on larger matrices%timeit Numba_csc_norm(X.data,X.indptr)#185 µs ± 5.39 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)