Python / Scipy / Numpy 中高效的增量稀疏矩阵

Question

在Python中，是否有一种有效的增量更新稀疏矩阵的方法？

 H = lil_matrix((n,m))
 for (i,j) in zip(A,B):
   h(i,j) += compute_something  

看起来这种方法构建稀疏矩阵速度相对较慢（lil_matrix 是最适合的稀疏矩阵类型）。

是否有一种高效构建稀疏矩阵 H 的方法（如使用字典或其他方式）？

- francois rousseau

---

H、A 和 B 的大小是多少？ - Cleb

稀疏矩阵的大小大约为100000 x 100000。因此，dim(A)可能接近于10^6。 - francois rousseau

Answer 1

在https://dev59.com/Zl4c5IYBdhLWcg3wj7F6#27771335中，我介绍了增量矩阵赋值。

lol和dok是如果您想更改值的推荐格式。 csr将为您提供效率警告，并且coo不允许索引。

但我也发现与常规字典索引相比，dok索引较慢。因此，对于许多更改，最好构建一个普通字典（具有相同的元组索引），并从该字典构建dok矩阵。

但是，如果您可以使用快速的numpy向量操作计算H数据值，而不是迭代，那么最好这样做，并从其构建稀疏矩阵（例如，使用coo格式）。实际上，即使使用迭代，这也将更快：

 h = np.zeros(A.shape)
 for k, (i,j) in enumerate(zip(A,B)):
    h[k] = compute_something 
 H = sparse.coo_matrix((h, (A, B)), shape=(n,m))

e.g.

In [780]: A=np.array([0,1,1,2]); B=np.array([0,2,2,1])
In [781]: h=np.zeros(A.shape)
In [782]: for k, (i,j) in enumerate(zip(A,B)):
    h[k] = i+j+k
   .....:     
In [783]: h
Out[783]: array([ 0.,  4.,  5.,  6.])
In [784]: M=sparse.coo_matrix((h,(A,B)),shape=(4,4))
In [785]: M
Out[785]: 
<4x4 sparse matrix of type '<class 'numpy.float64'>'
    with 4 stored elements in COOrdinate format>
In [786]: M.A
Out[786]: 
array([[ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  9.,  0.],
       [ 0.,  6.,  0.,  0.],
       [ 0.,  0.,  0.,  0.]])

请注意，(1,2)的值是4+5的总和。这是从coo到csr的转换的一部分。

在这种情况下，我可以使用以下方法计算h：

In [791]: A+B+np.arange(A.shape[0])
Out[791]: array([0, 4, 5, 6])

因此，无需迭代。

- hpaulj

Answer 2

不要使用csr_matrix或csc_matrix，因为如果你逐步构建它们，它们会比lil_matrix更慢。基于键的稀疏矩阵字典正是你需要的。

from scipy.sparse import dok_matrix
S = dok_matrix((5, 5), dtype=np.float32)
for i in range(5):
    for j in range(5):
        S[i,j] = i+j    # Update elements

- romeric

---

我在其他的SO问题中发现，对于具有相同键值的字典来说，索引普通字典比索引dok（Dictionary of Keys）字典更快。 - hpaulj

我已经测试了lil、csc、csr和dok矩阵。其中，lil_matrix是迄今为止最快的。 - francois rousseau

我相信，如果你看一下scipy的dok格式的实现，它本质上是dict的一个包装器。此外，dok的__getitem__和__setitem__方法是用纯Python编写的，因此索引它们会引入更多的开销。然而，与其他稀疏格式相比，dok在迭代时相对较便宜，因为项访问在算法上仍然是O(1)/常数级别的。老实说，我认为稀疏矩阵并不是真正意义上的迭代对象。 - romeric

Answer 3

一个更快的方法是：

H_ij = compute_something_vectorized()
H = coo_matrix((H_ij, (A, B))).tocsr()

重复坐标的数据将被求和，详见coo_matrix文档。

- user2379410

---

如果能制作出向量化函数，那将是最好的。 - hpaulj

我想将这个函数向量化，但不确定能否实现！ - francois rousseau