我继承了 R 的一些代码,但它运行得非常慢。大部分时间花在评估该形式的函数上(大约有 15 个具有不同被积函数 G 的此类函数):
TMin <- 0.5
F <- function (t, d) {
result <- ifelse(((d > 0) & (t > TMin)),
mapply(function(t, d) integrate(G, lower=0, upper=t, t, d)$value, t, d),
0)
return(result)
}
为了测试,我使用以下虚拟函数,但在实际代码中,Gs 要复杂得多,涉及 exp()、log()、dlnorm()、plnorm() 等。
G <- function(x, t, d) {
mean(rnorm(1e5))
x + t - d
}
在最坏的情况下,F 将被计算大约 200 万次。
该函数可以通过 3 种不同的方式调用:
t 是单个数字,d 是数值向量,或者,
t 是数值向量,d 是单个数字,或者,
t是一个数值向量并且是一个数值向量
有没有一种(简单的)方法来加速这个功能?
到目前为止,我已经尝试过一些变体(以摆脱 ifelse 循环):
F2 <- function (t,d) {
TempRes <- mapply(function(t, d) integrate(G, lower=0, upper=t, t, d)$value, t, d)
TempRes[(d <= 0) | (t <= TMin)] <- 0
result <- TempRes
return(result)
}
和
F3 <- function (t,d) {
result <- rep(0, max(length(t),length(d)))
test <- ((d > 0) & (t > TMin))
result[test] <- mapply(function(t, d) integrate(G, lower=0, upper=t, t, d)$value, t, d)[test]
return(result)
}
但它们花费的时间几乎完全相同。
最佳答案
您正在执行大量独立集成。您可以通过同时在单独的内核上执行这些集成来加快速度(如果您有可用的多核处理器)。问题是 R 默认情况下以单线程方式执行计算。但是,有许多可用的软件包允许多线程支持。我最近回答过几个类似的问题here和 here ,以及有关相关包和功能的一些附加信息。
此外,正如 @Mike Dunlavey 已经提到的,您应该避免对 t
的值执行积分。和d
不符合您的标准。 (您当前正在对这些值执行不需要的函数计算,然后用 0 覆盖结果)。
我在下面添加了可能的改进。请注意,您必须使用函数 G
创建一个单独的文件。包括在内以便在集群节点上对其进行评估。在下面的代码中,假设该文件名为 functionG.R
代码片段:
library(doParallel)
F4 <- function(t,d) {
results = vector(mode="numeric",max(length=length(t),length(d))) # Zero vector
logicalVector <- ((d > 0) & (t > TMin))
relevantT <- t[logicalVector]
relevantD <- d[logicalVector] # when d is single element, NA values created
if(length(relevantT) > 1 | length(relevantD) > 1)
{
if(length(d)==1) # d is only one element instead of vector --> replicate it
relevantD <- rep(d,length(relevantT))
if(length(t)==1) # t is only one element instead of vector --> replicate it
relevantT <- rep(t,length(relevantD))
cl <- makeCluster(detectCores());
registerDoParallel(cl)
clusterEvalQ(cl,eval(parse("functionG.R")))
integrationResults <- foreach(i=1:length(relevantT),.combine="c") %dopar%
{
integrate(G,lower=0,upper=relevantT[i],relevantT[i],relevantD[i])$value;
}
stopCluster(cl)
results[logicalVector] <- integrationResults
}
else if(length(relevantT==1)) # Cluster overhead not needd
{
results[logicalVector] = integrate(G,lower=0,upper=relevantT,relevantT,relevantD)$value;
}
return(results)
}
我的 CPU 包含 6 个启用超线程的物理核心 (x2)。结果如下:
> t = -5000:20000
> d = -5000:20000
>
> start = Sys.time()
> testF3 = F3(t,d)
> timeNeededF3 = Sys.time()-start
>
> start = Sys.time()
> testF4 = F4(t,d)
> timeNeededF4 = Sys.time()-start;
> timeNeededF3
Time difference of 3.452825 mins
> timeNeededF4
Time difference of 29.52558 secs
> identical(testF3,testF4)
[1] TRUE
运行此代码时,内核似乎一直在使用。但是,您可以通过更有效地围绕核心预分割数据来进一步优化此代码,然后在单独的核心上使用应用类型函数。
如果需要更多优化,您还可以更深入地查看 integrate
功能。您可以尝试这些设置,并通过允许不太严格的数值近似来获得性能增益。作为替代方案,您可以实现自己的简单版本的自适应辛普森求积并使用离散步长。您很可能会像这样获得巨大的性能提升(如果您能够/愿意在近似值中允许更多错误)。
编辑:
更新了代码以使其能够在所有场景中工作:d
和/或 t
有效/无效数字或向量
回复评论
@mawir:你是对的。 ifelse(test, yes, no)
会返回对应的yes
测试评估结果为 TRUE
的行的值,它将返回相应的 no
test
的行的值评估为FALSE
。但是,它首先必须评估您的 yes
表达式以创建 yes
向量length(test)
。这段代码演示了这一点:
> t = -5000:5
> d = -5000:5
>
> start = Sys.time()
> testF1 = F(t,d)
> timeNeededF1 = Sys.time()-start
> timeNeededF1
Time difference of 43.31346 secs
>
> start = Sys.time()
> testF4 = F4(t,d)
> timeNeededF4 = Sys.time()-start
> timeNeededF4
Time difference of 2.284134 secs
仅 t
的最后 5 个值和d
与此场景相关。
然而,在 F1
里面功能ifelse
评估 mapply
全部d
和t
首先值以创建 yes
向量。这就是函数执行时间如此长的原因。接下来,它选择满足条件的元素,否则选择 0。 F4
函数可以解决这个问题。
此外,您说您在 t
的情况下获得了加速和d
是非向量。然而,在这种情况下,不使用并行化。通常,您应该在满足 t
之一或两者的情况下获得最大加速比/d
是向量。
另一个编辑,回应罗兰的评论:
您可以替换 clusterEvalQ(cl,eval(parse("functionG.R")))
与 clusterExport(cl,"G")
如果您不想创建单独的函数文件。
关于r - 加速涉及映射和集成的函数,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/30951515/