在处理完这个和其他回复之后,这里的优化策略(和近似加速)似乎是
- (30x) 选择适当的数据表示形式——矩阵,而不是 data.frame
- (1.5x) 减少不必要的数据副本——列的差异,而不是 rowMeans
- 循环结构作为
*apply
函数(强调代码结构,简化内存管理,并提供类型一致性)
- (2x) 循环外的提升向量操作——列上的 abs 和 sum 变成矩阵上的 abs 和 colSums
整体加速约 100 倍。对于这种大小和复杂性的代码,使用编译器或并行包将无效。
我把你的代码放到一个函数中
f0 <- function(x) {
y <- rowMeans(x)
totaldiff <- numeric()
for (i in 1:ncol(x)){
x.after <- x
x.after[,i] <- sample(x[,i])
diff <- abs(y-rowMeans(x.after))
totaldiff[i] <- sum(diff)
}
which.max(totaldiff)
}
在这里我们有
x <- data.frame(matrix(runif(50*100),nrow=50,ncol=100)) #larger example
set.seed(123)
system.time(res0 <- f0(x))
## user system elapsed
## 1.065 0.000 1.066
您的数据可以表示为矩阵,对 R 矩阵的操作比对 data.frames 的操作要快。
m <- matrix(runif(50*100),nrow=50,ncol=100)
set.seed(123)
system.time(res0.m <- f0(m))
## user system elapsed
## 0.036 0.000 0.037
identical(res0, res0.m)
##[1] TRUE
这可能是最大的加速。但是对于这里的具体操作,我们不需要计算更新矩阵的行均值,只需要计算改组一列后均值的变化
f1 <- function(x) {
y <- rowMeans(x)
totaldiff <- numeric()
for (i in 1:ncol(x)){
diff <- abs(sample(x[,i]) - x[,i]) / ncol(x)
totaldiff[i] <- sum(diff)
}
which.max(totaldiff)
}
for
循环没有遵循正确的模式来填充结果向量(totaldiff
你想“预分配和填充”,所以totaldiff <- numeric(ncol(x))
)但我们可以使用 ansapply
并让 R 担心(这种内存管理是使用 apply 系列函数)
f2 <- function(x) {
totaldiff <- sapply(seq_len(ncol(x)), function(i, x) {
sum(abs(sample(x[,i]) - x[,i]) / ncol(x))
}, x)
which.max(totaldiff)
}
set.seed(123); identical(res0, f1(m))
set.seed(123); identical(res0, f2(m))
时间是
> library(microbenchmark)
> microbenchmark(f0(m), f1(m), f2(m))
Unit: milliseconds
expr min lq median uq max neval
f0(m) 32.45073 33.07804 33.16851 33.26364 33.81924 100
f1(m) 22.20913 23.87784 23.96915 24.06216 24.66042 100
f2(m) 21.02474 22.60745 22.70042 22.80080 23.19030 100
@flodel 指出vapply
可以更快(并提供类型安全)
f3 <- function(x) {
totaldiff <- vapply(seq_len(ncol(x)), function(i, x) {
sum(abs(sample(x[,i]) - x[,i]) / ncol(x))
}, numeric(1), x)
which.max(totaldiff)
}
然后
f4 <- function(x)
which.max(colSums(abs((apply(x, 2, sample) - x))))
仍然更快(ncol(x)
是一个常数因素,因此被删除)- abs
andsum
被提升到 之外sapply
,可能会以额外的内存使用为代价。注释中关于编译函数的建议总体上是好的;这里还有一些时间
> microbenchmark(f0(m), f1(m), f1.c(m), f2(m), f2.c(m), f3(m), f4(m))
Unit: milliseconds
expr min lq median uq max neval
f0(m) 32.35600 32.88326 33.12274 33.25946 34.49003 100
f1(m) 22.21964 23.41500 23.96087 24.06587 24.49663 100
f1.c(m) 20.69856 21.20862 22.20771 22.32653 213.26667 100
f2(m) 20.76128 21.52786 22.66352 22.79101 69.49891 100
f2.c(m) 21.16423 21.57205 22.94157 23.06497 23.35764 100
f3(m) 20.17755 21.41369 21.99292 22.10814 22.36987 100
f4(m) 10.10816 10.47535 10.56790 10.61938 10.83338 100
其中“.c”是编译版本和
编译在用 for 循环编写的代码中特别有用,但对矢量化代码没有多大作用;这里显示了编译 f1 的 for 循环的一个小但一致的改进,但不是 f2 的 sapply。