使用parallel优化数组R的循环

Question

我有一个数组data = array[1:50,1:50,1:50]里面的值是介于-1和1之间的实数

"Data"可以看成是一个50x50x50的立方体。

我需要根据这个方程创建一个相关矩阵(删除所有零)=>

value = (x+y)-|x-y|矩阵大小是2倍的可能组合(50x50x50)*((50x50x50)-1)/2 = 7.812.437.500这2倍=相关矩阵

我这样做了:

假设我们有3x3x3:

arr = array(rnorm(10), dim=c(3,3,3))
data = data.frame(array(arr))

data$voxel <- rownames(data) 
#remove zeros
data<-data[!(data[,1]==0),]
rownames(data) = data$voxel
data$voxel = NULL

#######################################################################################
#Create cluster
no_cores <- detectCores() #- 1
clus <- makeCluster(no_cores)
clusterExport(clus, list("data") , envir=environment())
clusterEvalQ(clus,
             compare_strings <- function(j,i) {
               value <- (data[i,]+data[j,])-abs(data[i,]- data[j,])   
               pair <- rbind(rownames(data)[j],rownames(data)[i],value)
               return(pair)
             })
i = 0 # start 0
kk = 1
table <- data.frame()
ptm <- proc.time()
while(kk<nrow(data)) {
  out <-NULL  
  i = i+1 # fix row
  j = c((kk+1):nrow(data)) # rows to be compared
  #Apply the declared function  
  out = matrix(unlist(parRapply(clus,expand.grid(i,j), function(x,y) compare_strings(x[1],x[2]))),ncol=3, byrow = T)
  table <- rbind(table,out)
  kk = kk +1
}
proc.time() - ptm

结果是data.frame:

v1  v2  v3
1   2   2.70430114250358
1   3   0.199941717684129
... up to 351 rows

但这需要几天…

我还想为这个相关性创建一个矩阵:

   1                         2              3...
1  1                  2.70430114250358 
2  2.70430114250358          1
3...

有更快的方法吗?

谢谢

Answer 1

您的代码中存在许多性能错误:

1. 当你应该依赖矢量化时你循环。
2. 在循环中生成对象。
3. 你并行化循环的每一次迭代，而不是并行化外部循环。

如果你避免了第一个问题，你就可以避免所有这些问题。

显然，您想比较每个行组合。为此，您应该首先获得行索引的所有组合:

combs <- t(combn(1:27, 2))

那么你可以对这些应用比较函数:

compare <- function(j,i, data) {
  as.vector((data[i,]+data[j,])-abs(data[i,]- data[j,]))
}
res <- data.frame(V1 = combs[,1], V2 = combs[,2], 
                  V3 = compare(combs[,1], combs[,2], data))

现在，如果我们想要检查这是否与您的代码给出相同的结果，我们首先需要修复您的输出。通过在矩阵中组合字符(行名)和数字，您得到一个字符矩阵，并且最终数据帧的列都是字符。我们可以在之后使用type.convert来修复这个问题(尽管从一开始就应该避免):

table[] <- lapply(table, function(x) type.convert(as.character(x)))

现在我们可以看到结果是一样的:

all.equal(res, table)
#[1] TRUE

如果你愿意，你可以把结果变成一个稀疏矩阵:

library(Matrix)
m <- sparseMatrix(i = res$V1, j = res$V2, x = res$V3, 
                  dims = c(27, 27), symmetric = TRUE)
diag(m) <- 1