如何在大稀疏矩阵中找到非零元素的索引？

Question

我有两个大小为100000 X 100000的方阵（a和b）。我需要计算这两个矩阵的差（c = a-b）。得到的矩阵'c'是一个稀疏矩阵，也就是说大部分元素都是零。我想找到所有非零元素的索引。我需要多次进行这个操作（超过100次）。

最简单的方法是使用两个循环来处理。但是这样计算量很大。你能告诉我有没有什么算法或者库，最好是用R、Python或C语言，能尽可能快地完成这个操作吗？

Answer 1

你可以使用 c.nonzero() 这个方法：

>>> from scipy.sparse import lil_eye
>>> c = lil_eye((4, 10)) # as an example
>>> c
<4x10 sparse matrix of type '<type 'numpy.float64'>'
        with 4 stored elements in LInked List format>
>>> c.nonzero()
(array([0, 1, 2, 3], dtype=int32), array([0, 1, 2, 3], dtype=int32))
>>> import numpy as np
>>> np.ascontiguousarray(c)
array([  (0, 0) 1.0
  (1, 1)        1.0
  (2, 2)        1.0
  (3, 3)        1.0], dtype=object)

你不需要先计算出 c 矩阵，也能找到 c = a - b 中非零元素的索引；你可以直接用 (a != b).nonzero() 来实现：

>>> a = np.random.random_integers(2, size=(4,4))
>>> b = np.random.random_integers(2, size=(4,4))
>>> (a != b).nonzero()
(array([0, 0, 1, 1, 1, 2, 3]), array([1, 2, 1, 2, 3, 2, 0]))
>>> a - b
array([[ 0,  1,  1,  0],
       [ 0,  1, -1, -1],
       [ 0,  0,  1,  0],
       [-1,  0,  0,  0]])

Answer 2

在R语言中，如果你使用Matrix这个包，并且用sparseMatrix把坐标列表转换成稀疏矩阵，那么你可以通过以下方式把它转换回三列：

TmpX <- as(M, "dgTMatrix")
X3col <- matrix(c(TmpX@i, TmpX@j, TmpX@val), ncol = 3)

这样你就能得到稀疏矩阵中的坐标和数值。

根据矩阵A和B中非零元素的位置，有时候用坐标列表处理会比用稀疏矩阵表示更好，因为这样你可以直接利用向量化操作，而不需要依赖稀疏矩阵包来达到最佳性能。我通常会根据需要的算法性能，在不同的编程语言中交替使用COO格式或稀疏矩阵支持。

---

更新1：我之前不知道你的两个矩阵A和B是稠密的。因此，找到矩阵C中非零元素的最简单方法其实是先不进行减法 - 直接比较A和B的元素。逻辑比较应该比减法更快。首先，找出A和B中不相等的元素A != B，然后只对这些元素进行减法。接下来，你只需要把A和B中索引的向量化形式转换成它们的（行，列）表示。这和Matlab中的ind2sub和sub2ind类似 - 你可以查看这个R参考文档来了解计算方法。

Answer 3

因为你有两个密集矩阵，所以双重循环是你唯一的选择。你根本不需要稀疏矩阵的类，因为你只想知道哪些位置 (i,j) 的值是不同的，也就是 a[i,j] != b[i,j]。

在像R和Python这样的编程语言中，双重循环的效率会比较低。我可能会用更底层的代码来写这个双重循环，然后把不同的索引加到一个列表里。不过，毫无疑问，R、Python等这些解释型语言的高手们肯定有更高效的方法来实现这个，而不需要使用底层代码。