使用并行优化数组 R 的循环

Question

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

“数据”可以被视为 50x50x50 的立方体。

我需要根据这个方程创建一个相关矩阵（去除所有零）=>

值 = (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))

现在，如果我们想检查这是否与您的代码给出相同的结果，我们首先需要修复您的输出。通过将字符（行名）与矩阵中的数字组合，您将获得一个字符矩阵，并且最终 data.frame 的列都是字符。之后我们可以使用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