How to find indices of non zero elements in large sparse matrix?

i have two sq matrix (a, b) of size in order o开发者_运维知识库f 100000 X 100000. I have to take difference of these two matrix (c = a-b). Resultant matrix 'c' is a sparse matrix. I want to find the indices of all non-zero elements. I have to do this operation many times (>100).

Simplest way is to use two for loops. But that's computationally intensive. Can you tell me any algorithm or package/library preferably in R/python/c to do this as quickly as possible?

---

Since you have two dense matrices then the double for loop is the only option you have. You don't need a sparse matrix class at all since you only want to know the list of indices (i,j) for which a[i,j] != b[i,j].

In languages like R and Python the double for loop will perform poorly. I'd probably write this in native code for a double for loop and add the indices to a list object. But no doubt the wizards of interpreted code (i.e. R, Python etc.) know efficient ways to do it without resorting to native coding.

---

In R, if you use the Matrix package, and sparseMatrix for the conversion from the coordinate list to the sparse matrix, then you can convert back to the 3 column via:

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

This will give you the coordinates and values in the sparse matrix.

Depending on the locations of non-zero entries in A and B, you may find it much better to work with the coordinate list than the sparse matrix representation (there are, by the way, dozens of sparse matrix representations), as you can take direct advantage of vectorized operations, rather than rely upon your sparse matrix package to perform optimally. I tend to alternate between using the COO or sparse matrix support in different languages, depending on how I will get the fastest performance for the algorithm of interest.

---

Update 1: I was unaware that your two matrices, A and B, are dense. As such, the easiest solution for finding non-zero entries in C is quite simply to not even subtract at first - just compare the entries of A and B. A logical comparison should be faster than a subtraction. First, find the entries of A and B where A != B, then subtract just those entries. Next, you simply need to convert from the vectorization of indices in A and B to their (row, col) representation. This is similar to ind2sub and sub2ind of Matlab - take a look at this R reference for the calculations.

---

You could use c.nonzero() method:

>>> 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)

You don't need to calculate c matrix to find out indexes of non-zero elements in c = a - b; you could do (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]])

---

I haven't timed it, but the simplest code is

all.indices<- which (C>0, arr.ind=T)

---

This code takes less then 0.1s.

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

EDIT: For sparse matrices of any size (which fits memory).

DATA

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)

CODE

system.time(tw <- ta[tb][is.na(a)|is.na(b)|(a-b != 0),list(r=r,c=c)])

---

have a look at numpy it have everything you ask for and more!

See this for sparse matrix support