6

我有两个大小为 100000 X 100000 的平方矩阵（a，b）。我必须取这两个矩阵（c = ab）的差异。结果矩阵'c'是一个稀疏矩阵。我想找到所有非零元素的索引。我必须多次执行此操作（> 100）。

最简单的方法是使用两个 for 循环。但这是计算密集型的。你能告诉我任何算法或包/库最好在 R/python/c 中尽快做到这一点吗？

4

6 回答 6

4

由于您有两个密集矩阵，因此双 for 循环是您唯一的选择。您根本不需要稀疏矩阵类，因为您只想知道(i,j).a[i,j] != b[i,j]

在 R 和 Python 等语言中，双重 for 循环的性能会很差。我可能会在本机代码中为双 for 循环编写此代码，并将索引添加到列表对象。但毫无疑问，解释代码（即 R、Python 等）的向导知道有效的方法来做到这一点，而无需求助于本机编码。

于 2011-09-09T13:21:22.190 回答
3

在 R 中，如果您使用Matrix包，并将坐标列表sparseMatrix转换为稀疏矩阵，则可以通过以下方式转换回第 3 列：

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 where 的条目A != B，然后仅减去这些条目。接下来，您只需将 A 和 B 中索引的矢量化转换为它们的 (row, col) 表示。这类似于 Matlab 的 ind2sub 和 sub2ind - 请查看此 R 参考以进行计算。

于 2011-09-09T12:54:00.203 回答
1

看看numpy它有你想要的一切，还有更多！

看到这个稀疏矩阵支持

于 2011-09-09T12:15:50.787 回答
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]])
于 2011-09-09T12:25:57.310 回答
1

这段代码耗时不到 0.1 秒。

m <- matrix(rpois(1000000,0.01),ncol=1000)
m0 <- lapply(seq(NCOL(m)),function(x) which(m[,x] != 0))

编辑：对于任何大小的稀疏矩阵（适合内存）。

数据

library(data.table)

N <- 1e+5
n <- 1e+6

ta <- data.table(r=sample(seq(N), n,replace=TRUE),
                 c=sample(seq(N), n,replace=TRUE),
                 a=sample(1:20,n,replace=TRUE))
tb <- data.table(r=sample(seq(N), n,replace=TRUE),
                 c=sample(seq(N), n,replace=TRUE),
                 b=sample(1:20,n,replace=TRUE))
setkey(ta,r,c)
setkey(tb,r,c)

代码

system.time(tw <- ta[tb][is.na(a)|is.na(b)|(a-b != 0),list(r=r,c=c)])
于 2011-09-09T13:27:09.630 回答
1

我没有计时，但最简单的代码是

all.indices<- which (C>0, arr.ind=T)
于 2011-09-09T15:20:43.507 回答