提高在大型矩阵中计算加权 Jaccard 的性能
Improve performance for computing Weighted Jaccard in a large matrix
R 输入：矩阵（测量 x 样本）(2291 x 265)（矩阵 [i,j]=0 和 1 之间的值）
输出：在所有样本对之间计算的加权 jaccard 的对称相似度矩阵
问题：找到产生输出的最快方法。我找到了使用 "doParallel" 和 "foreach" 的好方法，但这还不够，因为它仍然太慢。我没有找到任何具有能够计算加权 jaccard 功能的包，但也许我错过了它。无论如何，你可以用你喜欢的解决方案和方法回复。感谢大家会回答。
这是我现在的脚本：
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)
我没有可以为您 运行 完成的版本，因为我不完全确定输入的内容以及所需的输出应该是什么。但是，我确实有一些提示可以显着加快您的代码速度。
步骤 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 很慢等等。顺便说一句，我遇到了同样的问题（计算矩阵中所有列的所有成对乘积）。您可以在 MESS 包中以 pairwise_combination_indices 的形式访问该函数（并且您需要 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]
试试这些，看看是否有帮助？
顺便说一句，加权 Jaccard 相似度可以在 Rcpp 中非常快速地计算出所有对。
我找到了一个非常好的解决方案，它替换了所有原始代码并在几行代码中解决了问题。
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)
使用你的回答和@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
主要区别在于我首先计算我们需要计算值的唯一索引组合
R 输入：矩阵（测量 x 样本）(2291 x 265)（矩阵 [i,j]=0 和 1 之间的值）
输出：在所有样本对之间计算的加权 jaccard 的对称相似度矩阵
问题：找到产生输出的最快方法。我找到了使用 "doParallel" 和 "foreach" 的好方法，但这还不够，因为它仍然太慢。我没有找到任何具有能够计算加权 jaccard 功能的包，但也许我错过了它。无论如何，你可以用你喜欢的解决方案和方法回复。感谢大家会回答。 这是我现在的脚本：
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)
我没有可以为您 运行 完成的版本，因为我不完全确定输入的内容以及所需的输出应该是什么。但是，我确实有一些提示可以显着加快您的代码速度。
步骤 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 很慢等等。顺便说一句，我遇到了同样的问题（计算矩阵中所有列的所有成对乘积）。您可以在 MESS 包中以 pairwise_combination_indices 的形式访问该函数（并且您需要 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]
试试这些，看看是否有帮助？
顺便说一句，加权 Jaccard 相似度可以在 Rcpp 中非常快速地计算出所有对。
我找到了一个非常好的解决方案，它替换了所有原始代码并在几行代码中解决了问题。
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)
使用你的回答和@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
主要区别在于我首先计算我们需要计算值的唯一索引组合