Two-way frequency table followed by matrix multiplication - high running time

user_18893821 0
Viewed 34

I'm new to R, and trying to calculate the product between a fixed matrix to a 2-way frequency table for any combinations of columns in a dataframe or matrix and divide it by the sequence length (aka number of rows which is 15), the problem is that the running time increases dramatically when performing it on 1K sequences (1K columns). the goal is to use it with as much as possible sequences (more than 10 minutes, for 10K could be more than 1hr)

mat1 <- matrix(sample(LETTERS),ncol = 100,nrow = 15)
mat2 <- matrix(sample(abs(rnorm(26,0,3))),ncol=26,nrow=26)
rownames(mat2) <- LETTERS
colnames(mat2) <- LETTERS
diag(mat2) <- 0

test_vec <- c()
for (i in seq(ncol(mat1)-1)){  
  for(j in seq(i+1,ncol(mat1))){
    
    s2 <- table(mat1[,i],mat1[,j]) # create 2-way frequency table
    mat2_1 <- mat2
    mat2_1 <- mat2_1[rownames(mat2_1) %in% rownames(s2), 
                               colnames(mat2_1) %in% colnames(s2)]
    calc <- ((1/nrow(mat1))*sum(mat2_1*s2))
    test_vec <- append(test_vec,calc)

  }}

Thanks for the help.

matrix-multiplication r matrix performance
1 Answers

Here is an approach that converts mat1 to a data.table, and converts all the columns to factors, and uses table(..., exclude=NULL)

library(data.table)
m=as.data.table(mat1)[,lapply(.SD, factor, levels=LETTERS)]
g = combn(colnames(m),2, simplify = F)
result = sapply(g, function(x) sum(table(m[[x[1]]], m[[x[2]]], exclude=NULL)*mat2)/nrow(m))

Check equality:

sum(result-test_vec>1e-10)
[1] 0

Here there are 4950 combinations (100*99/2), but the number of combinations will increase quickly as nrow(mat1) increases (as you point out). You might find in that case that a parallelized version works well.

library(doParallel)
library(data.table)
registerDoParallel()

m=as.data.table(mat1)[,lapply(.SD, factor, levels=LETTERS)]
g = combn(colnames(m),2, simplify = F)
result = foreach(i=1:length(g), .combine=c) %dopar%
  sum(table(m[[g[[i]][1]]], m[[g[[i]][2]]], exclude=NULL)*mat2)
result = result/nrow(m)
user_44475400
Related