Повышение производительности для вычисления взвешенной Jaccard в большой матрице - PullRequest
Новая
       51

Повышение производительности для вычисления взвешенной Jaccard в большой матрице

0 голосов
/ 29 апреля 2018

R вход: матрица (измеряет x выборок) (2291 x 265) (матрица [i, j] = значение от 0 до 1)

Вывод: симметричная матрица подобия взвешенной карты jaccard, вычисленная между всеми парами выборок

Проблема: найти самый быстрый способ получения результата. Я нашел хороший способ, используя «doParallel» и «foreach», но этого недостаточно, потому что он все еще слишком медленный. Я не нашел ни одного пакета с функцией, которая могла бы вычислить взвешенную карту, но, возможно, я пропустил это. В любом случае вы можете ответить с решением и методом, который вам нравится. Спасибо всем, отвечу. Это мой скрипт на данный момент: 

rm(list=ls())

#Load libraries ----
require(doParallel)
require(foreach)
require(doSNOW)
require(doMPI)
#Imported data ----
dim(input_m) #2291 x 265

#Set clusters ----
no_cores <- 3
cl <- makeCluster(as.integer(no_cores))
registerDoParallel(cl)

#I build all the combinations of the pairs of samples ----
samples=seq(1:ncol(input_m))
combs<-as.matrix(expand.grid(samples,samples))
combs<-unique(t(parApply(cl=cl, combs, 1, sort)))

#Prepare the resulting matrix ----
res_m <- matrix(ncol = ncol(input_m), nrow = ncol(input_m))
rownames(res_m)=colnames(input_m)
colnames(res_m)=colnames(input_m)

#Compute Weighted Jaccard similarity btw all pairs of samples ----
sim_m=foreach(s = 1:nrow(combs), .combine=rbind, .noexport=c("pair","num","den"), .inorder=FALSE) %dopar% {
    pair=input_m[,c(combs[s,1],combs[s,2])]
    num=sum(apply(pair,1,min))
    den=sum(apply(pair,1,max))
    return(c(combs[s,1],combs[s,2],num/den))
}

#Fill the prepared matrix with the results in sim_m
for (k in 1:nrow(sim_m)){
    sim=sim_m[k,3]
    idx1=sim_m[k,1]
    idx2=sim_m[k,2]
    res_m[idx1,idx2]=sim
    res_m[idx2,idx1]=sim
}

#Stop clusters
stopCluster(cl)

Ответы [ 3 ]

0 голосов
/ 30 апреля 2018

Я нашел очень хорошее решение, которое заменяет весь исходный код и решает проблему в несколько строк кода. 

rm(list=ls())
load("data.rda")
# dim(input_m) 2291 x 265
res_m=outer(1:ncol(input_m), 1:ncol(input_m) , FUN=Vectorize(function(r,c){
 require(matrixStats);
 sum(rowMins(input_m[,c(r,c)]))/sum(rowMaxs(input_m[,c(r,c)]))}))
rownames(res_m)=colnames(input_m)
colnames(res_m)=colnames(input_m)
0 голосов
/ 02 мая 2018

используя ваш ответ и комментарии @HenrikB Мне удалось написать более быстрый подход: 

## simulate data
nr <- 2291; nc <- 265
set.seed(420)
input_m <- matrix(rnorm(nr * nc), nrow = nr, ncol = nc)
input_m[1:5, 1:5]
#             [,1]       [,2]        [,3]        [,4]        [,5]
# [1,] -0.76774389  1.2623614  2.44166184 -1.86900934  1.61130129
# [2,] -1.44513238 -0.5469383 -0.31919480 -0.03155421  0.09293325
# [3,] -0.71767075 -0.2753542  2.28792301  0.41545393 -0.47370802
# [4,]  0.06410398  1.4956864  0.06859527  2.19689076 -0.96428109
# [5,] -1.85365878  0.1609678 -0.52191522 -0.79557319 -0.33021108

jaccardLuke <- function(input_m) {
  res_m = outer(1:ncol(input_m), 1:ncol(input_m) ,
                FUN = Vectorize(function(r,c) {
                  require(matrixStats)
                  sum(rowMins(input_m[,c(r,c)]))/sum(rowMaxs(input_m[,c(r,c)]))
                  })
                )
  rownames(res_m) = colnames(input_m)
  colnames(res_m) = colnames(input_m)
  res_m
}

jaccardHenrikB <- function(input_m) {
  require(matrixStats)
  res_m = outer(1:ncol(input_m), 1:ncol(input_m) ,
                FUN = Vectorize(function(r, r2) {
                  x <- rowRanges(input_m, cols = c(r, r2))
                  s <- colSums(x)
                  s[1] / s[2]
                })
  )
  rownames(res_m) = colnames(input_m)
  colnames(res_m) = colnames(input_m)
  res_m
}

Моя функция: 

jaccardMinem <- function(input_m) {
  require(data.table)
  require(matrixStats)

  samples <- 1:ncol(input_m)
  comb <- CJ(samples, samples)
  comb[, i := .I]
  comb <- melt(comb, 'i')
  setorder(comb, value)
  v2 <- paste0("V", 1:2)
  comb[, variable2 := v2 , keyby = i]
  comb2 <- dcast(comb, i ~ variable2, value.var = 'value')
  combUnique <- unique(comb2, by = c('V1', 'V2'))

  XX <- apply(combUnique[, -'i'], 1, function(x) {
    x2 <- rowRanges(input_m, cols = x)
    s <- colSums2(x2)
    s[1] / s[2]
  })

  set(combUnique, j = 'xx', value = XX)
  rez2 <- merge(comb2, combUnique[, -'i'], by = c('V1', 'V2'), all.x = T)
  setorder(rez2, i)
  rez2 <- array(rez2$xx, dim = rep(ncol(input_m), 2))
  rownames(rez2) <- colnames(input_m)
  colnames(rez2) <- colnames(input_m)
  rez2
}

Проверьте, все ли равны: 

all.equal(jaccardLuke(input_m), jaccardHenrikB(input_m))
# [1] TRUE
all.equal(jaccardLuke(input_m), jaccardMinem(input_m))
# [1] TRUE

бенчмаркинг: 

system.time(jaccardLuke(input_m)) # 6.05 sek
system.time(jaccardHenrikB(input_m)) # 2.75 sek
system.time(jaccardMinem(input_m)) # 1.74 sek

## for larger data:
nr <- 5000; nc <- 500
set.seed(420)
input_m <- matrix(rnorm(nr * nc), nrow = nr, ncol = nc)
system.time(jaccardLuke(input_m)) # 41.55 sek
system.time(jaccardHenrikB(input_m)) # 19.87 sek
system.time(jaccardMinem(input_m)) # 11.17 sek

главное отличие в том, что я сначала вычисляю уникальные комбинации индексов, для которых нам нужно вычислить значения

0 голосов
/ 29 апреля 2018

У меня нет готовой версии, которая будет работать для вас, потому что я не совсем уверен, как выглядит ввод и каким должен быть желаемый вывод. Однако у меня есть некоторые подсказки, чтобы значительно ускорить ваш код.

Шаг 1

Ваша самая большая свинья это кусок кода 

samples=seq(1:ncol(input_m))
combs<-as.matrix(expand.grid(samples,samples))
combs<-unique(t(parApply(cl=cl, combs, 1, sort)))

expand.grid медленный, sort медленный и т. Д. И т. Д. Кстати, у меня была та же проблема (вычисление всех парных произведений всех столбцов в матрице). Вы можете получить доступ к функции как pairwise_combination_indices в пакете MESS (и вам нужна версия github): 

devtools::install_github("ekstroem/MESS")

Теперь посмотрите на это увеличение скорости. f() соответствует вашим трем строкам выше 

microbenchmark::microbenchmark(f(100), MESS::pairwise_combination_indices(100, self=TRUE))
Unit: microseconds
                                                 expr        min          lq
                                               f(100) 355670.517 386745.3550
 MESS::pairwise_combination_indices(100, self = TRUE)     31.006     44.3855
        mean     median         uq        max neval cld
 414465.6852 409732.726 427356.848 575404.135   100   b
     85.7078     65.962     84.804    679.408   100  a

Теперь вам нужно вычислить матрицу индексов для 265 столбцов, а не только для 100, поэтому прирост скорости должен быть еще больше. Никакое количество ядер не должно с этим конкурировать, поэтому замените три строки на 

combs <- MESS::pairwise_combination_indices(ncols(input_m), self=TRUE)

Шаг 2

Ваш последний цикл должен быть векторизован, и вы можете обойтись (не проверено) 

res_m[cbind(sim_m[k,1], sim_m[k,2])] = sim_m[k,3]
res_m[cbind(sim_m[k,2], sim_m[k,1])] = sim_m[k,3]

Попробуйте это и посмотрите, поможет ли это?

Сходство по весовой Жаккарте, между тем все пары, вероятно, могут быть очень быстро вычислены в Rcpp.

...