How to vectorize comparing each row of matrix with all other rows

Question

I am trying to compare each row with all other rows in a matrix to count the number of differences of each row with all other rows. The result is then stored in the bottom left triangle of a matrix.

So for example when row m[1,] is compared with rows m[2,] and m[3,] the difference counts are stored at positions of mat[c(2:3), 1] in the result matrix.

My problem is that my input matrix can have upto 1e+07 rows and the current implementation (speed and memory) will not scale due to n^2 comparisons. Suggestions and help would be appreciated.

diffMatrix <- function(x) {

  rows <- dim(x)[1] #num of rows
  cols <- dim(x)[2] #num of columns

  if (rows <= 1) stop("'x' must have atleast two rows")

  #potential failure point
  mat <- matrix(, rows, rows)

  # fill bottom left triangle columns ignoring the diagonal
  for (row in 1:(rows-1)) { 
    rRange <- c((1+row):rows)
    m <- matrix(x[row,], nrow=rows-row, ncol=cols, byrow = T)
    mat[rRange, row] <- rowSums(m != x[-1:-row, ])
  }

  return (mat)
}   


m <- matrix(sample(1:12, 12, replace=T), ncol=2, byrow=TRUE)
m
#     [,1] [,2]
#[1,]    8    1
#[2,]    4    1
#[3,]    8    4
#[4,]    4    5
#[5,]    3    1
#[6,]    2    2

x <- diffMatrix(m)
x
#     [,1] [,2] [,3] [,4] [,5] [,6]
#[1,]   NA   NA   NA   NA   NA   NA
#[2,]    1   NA   NA   NA   NA   NA
#[3,]    1    2   NA   NA   NA   NA
#[4,]    2    1    2   NA   NA   NA
#[5,]    1    1    2    2   NA   NA
#[6,]    2    2    2    2    2   NA

m <- matrix(sample(1:5, 50000, replace=T), ncol=10, byrow=TRUE)
# system.time(x <- diffMatrix(m))
#   user  system elapsed 
#  20.39    0.38   21.43

How much memory do you have? You may have to resort to writing out the results to disk. — Scott Ritchie, Jan 01 '14 at 11:39
A custom client-app is front-ending concurrent requests on a pool of R servers. Each request is invoking a R function that the diffMatrix needs to be part of. Hence low memory and super high speed is critical. — user3147662, Jan 02 '14 at 04:31
Any way you look at it, you're going to have to compute, and store, `n^2` comparisons. Unless you can precompute this matrix, and store it to be looked up by the client later, what you're asking is literally impossible. — Scott Ritchie, Jan 02 '14 at 04:34
True, but there should be some way to mix the expand.grid idea (http://stackoverflow.com/questions/19933788/r-compare-all-the-columns-pairwise-in-matrix) and my loop over half rows approach. I will try @alexis_laz's recommendation tomorrow. — user3147662, Jan 02 '14 at 05:55
You may be able to cut corners depending on what type of data the matrix contains: float, integer, integer within a limited range, logical, or something else? Is it sparse? — smci, Mar 31 '14 at 13:35

score 2 · Answer 1 · answered Jan 01 '14 at 01:47

Here is an alternative using .Call (seems valid, but I can't guarantee):

library(inline)

ff = cfunction(sig = c(R_mat = "matrix"), body = '
SEXP mat, dims, ans, dimans;

PROTECT(dims = getAttrib(R_mat, R_DimSymbol));
PROTECT(dimans = allocVector(INTSXP, 2));
R_len_t *pdims = INTEGER(dims), *pdimans = INTEGER(dimans);
PROTECT(ans = allocVector(INTSXP, pdims[0]*pdims[0]));
R_len_t *pans = INTEGER(ans);
pdimans[0] = pdims[0];
pdimans[1] = pdims[0];

for(int ina = 0; ina < LENGTH(ans); ina++) {
   pans[ina] = NA_INTEGER;
}

switch(TYPEOF(R_mat)) {
   case REALSXP:
   {
    PROTECT(mat = coerceVector(R_mat, REALSXP));
    double *pmat = REAL(mat);
    for(int i = 0; i < pdims[0]; i++) {
       R_len_t ir;
       for(ir = i+1; ir < pdims[0]; ir++) {
          R_len_t diffs = 0;
          for(int ic = 0; ic < pdims[1]; ic++) {
             if(pmat[i + ic*pdims[0]] != pmat[ir + ic*pdims[0]]) {
               diffs++;
             }
          }
          pans[ir + i*pdims[0]] = diffs;
       }
    }
    break;
   }

   case INTSXP:
   {
    PROTECT(mat = coerceVector(R_mat, INTSXP));
    R_len_t *pmat = INTEGER(mat);
    for(int i = 0; i < pdims[0]; i++) {
       R_len_t ir; 
       for(ir = i+1; ir < pdims[0]; ir++) {
          R_len_t diffs = 0;
          for(int ic = 0; ic < pdims[1]; ic++) {
             if(pmat[i + ic*pdims[0]] != pmat[ir + ic*pdims[0]]) {
               diffs++;
             }
          }
          pans[ir + i*pdims[0]] = diffs;
       }
    }
    break;
   }
}

setAttrib(ans, R_DimSymbol, dimans);

UNPROTECT(4);

return(ans);
')

m = matrix(c(8,4,8,4,3,2,1,1,4,5,1,2), ncol = 2)
ff(m)
#     [,1] [,2] [,3] [,4] [,5] [,6]
#[1,]   NA   NA   NA   NA   NA   NA
#[2,]    1   NA   NA   NA   NA   NA
#[3,]    1    2   NA   NA   NA   NA
#[4,]    2    1    2   NA   NA   NA
#[5,]    1    1    2    2   NA   NA
#[6,]    2    2    2    2    2   NA
all.equal(diffMatrix(m), ff(m))
#[1] TRUE

m2 <- matrix(sample(1:5, 50000, replace=T), ncol=10, byrow=TRUE)
library(microbenchmark)
microbenchmark(diffMatrix(m2), ff(m2), times = 10)
#Unit: milliseconds
#           expr       min        lq   median        uq        max neval
# diffMatrix(m2) 6972.9778 7049.3455 7427.807 7633.7581 11337.3154    10
#         ff(m2)  422.3195  469.5771  530.470  661.8299   862.3092    10
all.equal(diffMatrix(m2), ff(m2))
#[1] TRUE

If you're going to go the `C` route, I'd highly recommend checking out `Rcpp`. It provides a lot of syntactic sugar that will make the underlying C much easier to write. — Scott Ritchie, Jan 01 '14 at 11:37

How to vectorize comparing each row of matrix with all other rows

1 Answers1

Linked