大型Numpy Scipy CSR矩阵，按行操作

Question

我想遍历CSR矩阵的行，并将每个元素除以该行总和，类似于这里所示的操作：numpy divide row by row sum。

我的问题是我正在处理一个大矩阵：(96582, 350138)。当应用链接帖子中的操作时，由于返回的矩阵是密集的，它会膨胀我的内存。

因此，这是我的第一次尝试：

for row in counts:
    row = row / row.sum()

不幸的是，这并不会对矩阵产生影响，因此我想到了第二个方法，创建一个新的csr矩阵，并使用vstack连接行：

from scipy import sparse
import time

start_time = curr_time = time.time()
mtx = sparse.csr_matrix((0, counts.shape[1]))
for i, row in enumerate(counts):
   prob_row = row / row.sum()
   mtx = sparse.vstack([mtx, prob_row])
   if i % 1000 == 0:
      delta_time = time.time() - curr_time
      total_time = time.time() - start_time
      curr_time = time.time()
      print('step: %i, total time: %i, delta_time: %i' % (i, total_time, delta_time))

这个功能很好，但经过一些迭代后，它变得越来越慢：

step: 0, total time: 0, delta_time: 0
step: 1000, total time: 1, delta_time: 1
step: 2000, total time: 5, delta_time: 4
step: 3000, total time: 12, delta_time: 6
step: 4000, total time: 23, delta_time: 11
step: 5000, total time: 38, delta_time: 14
step: 6000, total time: 55, delta_time: 17
step: 7000, total time: 88, delta_time: 32
step: 8000, total time: 136, delta_time: 47
step: 9000, total time: 190, delta_time: 53
step: 10000, total time: 250, delta_time: 59
step: 11000, total time: 315, delta_time: 65
step: 12000, total time: 386, delta_time: 70
step: 13000, total time: 462, delta_time: 76
step: 14000, total time: 543, delta_time: 81
step: 15000, total time: 630, delta_time: 86
step: 16000, total time: 722, delta_time: 92
step: 17000, total time: 820, delta_time: 97

有什么建议吗？vstack为什么会越来越慢？

- nadre

---

请参考以下链接：https://stackoverflow.com/a/45339754 和 https://stackoverflow.com/q/44080315。 - hpaulj

与密集数组一样，循环中的重复连接速度较慢。在列表中累积结果并执行一个 vstack 更快。 - hpaulj

Answer 1

vstack 是一个 O(n) 操作，因为它需要为结果分配内存，然后将传递给它的所有数组的内容复制到结果数组中。

您可以简单地使用multiply 来执行此操作：

>>> res = counts.multiply(1 / counts.sum(1))  # multiply with inverse
>>> res.todense()
matrix([[ 0.33333333,  0.        ,  0.66666667],
        [ 0.        ,  0.        ,  1.        ],
        [ 0.26666667,  0.33333333,  0.4       ]])

---

但是使用np.lib.stride_tricks.as_strided也可以很容易地完成所需操作（相对高效）。此as_strided函数还允许在数组上执行更复杂的操作（如果没有您的情况的方法或函数）。

例如，使用scipy文档中的csr示例：

>>> from scipy.sparse import csr_matrix
>>> import numpy as np
>>> row = np.array([0,0,1,2,2,2])
>>> col = np.array([0,2,2,0,1,2])
>>> data = np.array([1.,2,3,4,5,6])
>>> counts = csr_matrix( (data,(row,col)), shape=(3,3) )
>>> counts.todense()
matrix([[ 1.,  0.,  2.],
        [ 0.,  0.,  3.],
        [ 4.,  5.,  6.]])

您可以按以下方式将每行除以其总和：

您可以按以下方式将每行除以其总和：

>>> row_start_stop = np.lib.stride_tricks.as_strided(counts.indptr, 
                                                     shape=(counts.shape[0], 2),
                                                     strides=2*counts.indptr.strides)
>>> for start, stop in row_start_stop:   
...    row = counts.data[start:stop]
...    row /= row.sum()
>>> counts.todense()
matrix([[ 0.33333333,  0.        ,  0.66666667],
        [ 0.        ,  0.        ,  1.        ],
        [ 0.26666667,  0.33333333,  0.4       ]])

- MSeifert

---

非常好的答案，你更新行的方式比我的高效得多。像一百倍！这应该是被采纳的答案。 - Anis

Answer 2

编辑

@MSeifert的答案更加高效，应该是正确的方法。我认为写counts[i, :] 意味着进行了某种列切片操作，而我之前没有意识到这一点。文档明确表示，在csr_matrix上执行这些操作非常低效。Way确实有一个很好的例子。

回答

文档说行切片是高效的，我认为你应该这样做

for i in range(counts.shape[0]):
    counts[i,:] /= counts[i,:].sum()

这种方式可以直接编辑您的矩阵，保持其稀疏性，并且无需使用vstack。我不确定这是否是最有效的操作，但至少您不应该有内存问题，并且在计算行时没有减速效应：

import time()
s = time.time()
for i in range(counts.shape[0]):
    counts[i, :] /= (counts[i, :].sum() + 1)
    if i % 1000 == 0:
        e = time.time()
        if i > 0:
            print i, e-s
        s = time.time()

性能表现:

1000 6.00199794769 2000 6.02894091606 3000 7.44459486008 4000 7.10011601448 5000 6.16998195648 6000 7.79510307312 7000 7.00139117241 8000 7.37821507454 9000 7.28075814247 ...

MSeifert 答案的性能表现:

row_start_stop = np.lib.stride_tricks.as_strided(counts.indptr, shape=(counts.shape[0], 2),
                                                 strides=2*counts.indptr.strides)
for i, (start, stop) in enumerate(row_start_stop):
    row = counts.data[start:stop]
    row /= row.sum()
    if i % 1000 == 0:
        e = time.time()
        if i > 0:
            print i,e-s
        s = time.time()

1000 0.00735783576965 2000 0.0108380317688 3000 0.0102109909058 4000 0.0131571292877 5000 0.00670218467712 6000 0.00608897209167 7000 0.00663685798645 8000 0.0164499282837 9000 0.0061981678009 ... 关于为什么使用vstack会变慢，@MSeifert的答案非常好。

- Anis

---

只是出于好奇，for item in counts: item /= item.sum() 怎么样？我已经相信被接受的答案应该是 Mseifert 的，这只是对隐式和显式切片性能差异的个人好奇。 - Uvar

问题在于OP的问题。当执行for item in counts时，item并不是对counts实际行的引用。这就是为什么修改不会影响矩阵计数的原因。所以我可以测量性能，但我真的看不出有什么意义。 - Anis

也许我对稀疏矩阵的理解还不够扎实。:/ 不管怎样，正如您在我的回答中所看到的，对于密集矩阵，您可以通过这种方式引用实际行。我以为它会以同样的方式起作用.. :$ - Uvar

不正确。OP在他的问题中尝试的方法也不适用于密集矩阵。 - Uvar

让我们在聊天室里继续这个讨论。 - Uvar