使用并行优化数组 R 的循环
Optimize loops for arrays R using parallel
我有一个数组 data = array[1:50,1:50,1:50] 里面的值是 -1, 1 之间的实数。
"Data" 可以视为 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...
有没有更快的方法?
谢谢
您的代码中存在一些性能错误:
- 当你应该依赖向量化时你循环了。
- 你在循环中增长一个对象。
- 您并行化了循环的每次迭代,而不是并行化了外部循环。
如果你避免了第一个问题,你就可以避免所有这些问题。
显然,您想要比较行的每个组合。为此,您应该首先获取行索引的所有组合:
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
我有一个数组 data = array[1:50,1:50,1:50] 里面的值是 -1, 1 之间的实数。
"Data" 可以视为 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...
有没有更快的方法?
谢谢
您的代码中存在一些性能错误:
- 当你应该依赖向量化时你循环了。
- 你在循环中增长一个对象。
- 您并行化了循环的每次迭代,而不是并行化了外部循环。
如果你避免了第一个问题,你就可以避免所有这些问题。
显然,您想要比较行的每个组合。为此,您应该首先获取行索引的所有组合:
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