Faster ways of creating matrices in R

user_2766780 0
Viewed 582

I have the following example of something that i've been doing, which is formally simple, but I wanted to check what are the potential alternatives to my code -- in order to get faster, if possible. Here it is the example:

Time1=Sys.time()
v=rep(c("A","B","C","D","E","F","G","H","I","J","K","L","M","N","O","P"),
     each=1000)
m=matrix(0,ncol=length(v),nrow=length(v))
for (j in 1:length(v)) {
  for(i in 1:length(v)) {
    if (v[j]==v[i]) {
      m[j,i]=1
    } else {
      next
    }
  }
}
Time2=Sys.time()
Time2-Time1
# Time difference of 1.405404 mins

I am creating a simple relational matrix -- where the vector v1 could be interpreted as being placed as lines and columns and the matrix maps where the results are equal. If they're equal, we get m[j,i]=1; if not equal, m[j,i]=0. As I stated, I would like to make this code go faster. I was trying think of a way to code it as an apply function, but I haven't figured that out for now. Still, I would like to know if there are other options besides what I've said.

EDIT: I made some corrections on the text and I tried to clarify the question.

r matrix
5 Answers

Sure. With the assumption that your sample data is not representative of real data, then this works about 6x faster:

m2 <- +(outer(v, v, `==`))
all.equal(m, m2)
# [1] TRUE

If, however, your real data has significant number of duplications, then @Sathish's method of removing duplications before comparison and propagating that through a matrix is likely much faster.

user_33582720

I think outer approach by @r2evans is the easiest way to construct the matrix. Below is another base R option, using expand.grid

m2 <- matrix(+do.call("==",expand.grid(v,v)),length(v))
user_121587570

I'll throw my hat in the ring with CJ from data.table.

libary(data.table)
m2 <- matrix(+(CJ(v1 = v,v2 = v,sorted=FALSE)[,ans := v1==v2][,ans]),length(v))
all.equal(m,m2)
#[1] TRUE
user_130953260

If you have a lot of zeros, then you can use a sparseMatrix. You fill in the positions which are matches, and leaving the rest zeros. This consumes less memory but you can only use it with certain functions. It's used for example in glmnet lasso.

I think @r2evans solution is the most concise and would work for most cases.

Below I have a few codes from the answers, most of them are indeed faster than the OP's

library(microbenchmark)
library(Matrix)
library(data.table)
setDTthreads(threads =1)

f_sw = function(v){
N = length(v)
i = lapply(v,function(i)which(v==i))
j = rep(1:N,times=sapply(i,length))
as.matrix(sparseMatrix(i=unlist(i),j=j,dims=list(N,N)))
}

f_r2evans = function(v){
  m2 <- +(outer(v, v, `==`))
  return(m2)
}

f_IanCampbell = function(v){ 
  matrix(+(CJ(v1 = v,v2 = v,sorted=FALSE)[,ans := v1==v2][,ans]),length(v))
}

microbenchmark(f_IanCampbell(v),f_sw(v),f_r2evans(v),times=5)

Unit: seconds
             expr       min        lq      mean    median        uq      max
 f_IanCampbell(v) 10.820616 11.325422 12.544146 12.983926 13.126655 14.46411
          f_sw(v)  7.014364  7.228585  8.206858  7.745741  8.877425 10.16818
     f_r2evans(v)  9.117405  9.519443  9.996789  9.896823 10.288586 11.16169
 neval cld
     5   b
     5  a 
     5  a
user_122584590

may be this:

Instead of checking for equality of all possible values, we can do it on the vector and replicate the matrix from it to 1000 times by row and by column. It will give the same result. The column and row order is not maintained by this code. But, using the row names and column names, we can verify the answer to be correct or not.

I used t(), because column bind is faster than row bind.

system.time({
  v1 <- c("A","B","C","D","E","F","G","H","I","J","K","L","M","N","O","P")
  m1 <- sapply(v1, function(x) as.integer(v1 == x))
  rownames(m1) <- colnames(m1)
  m1 <- do.call('cbind', mget(rep('m1', 1000)))
  m1 <- t(m1)
  m1 <- do.call('cbind', mget(rep('m1', 1000)))
  m1 <- t(m1)
}) 
# user  system elapsed 
# 9.32    0.50    9.84 


dim(m1)
# [1] 16000 16000

Another method:

This one (hard-coded) will not do any comparison, but we create values based on what can happen by comparing vector with its own values.

It can be improved by using eval(parse(text=paste())) construct. 

system.time({
  m4 <- matrix(data = 
           c(
             c(rep(1, 1000), rep(0, 15000)),
             c(rep(0, 1000), rep(1, 1000), rep(0, 14000)),
             c(rep(0, 2000), rep(1, 1000), rep(0, 13000)),
             c(rep(0, 3000), rep(1, 1000), rep(0, 12000)),
             c(rep(0, 4000), rep(1, 1000), rep(0, 11000)),
             c(rep(0, 5000), rep(1, 1000), rep(0, 10000)),
             c(rep(0, 6000), rep(1, 1000), rep(0, 9000)),
             c(rep(0, 7000), rep(1, 1000), rep(0, 8000)),
             c(rep(0, 8000), rep(1, 1000), rep(0, 7000)),
             c(rep(0, 9000), rep(1, 1000), rep(0, 6000)),
             c(rep(0, 10000), rep(1, 1000), rep(0, 5000)),
             c(rep(0, 11000), rep(1, 1000), rep(0, 4000)),
             c(rep(0, 12000), rep(1, 1000), rep(0, 3000)),
             c(rep(0, 13000), rep(1, 1000), rep(0, 2000)),
             c(rep(0, 14000), rep(1, 1000), rep(0, 1000)),
             c(rep(0, 15000), rep(1, 1000))), nrow = 16000, ncol = 16000)
})
# user  system elapsed 
# 0.72    0.93    1.82

Note: As @r2evans said, this will not work if the OP's sample data is non-representative of the real data

user_16917230
Related