加速涉及映射和集成的功能

Question

我继承了 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)

}

为了测试，我使用了以下虚拟函数，但在实际代码中，G 复杂得多，涉及 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)
}

但他们花费的时间几乎完全相同。

Answer 1

一般来说，要查看的位置在最内层的循环中，您可以通过减少它花费的时间或减少调用它的次数来加快它的速度。您有一个内部循环运行 mapply，但随后您从中提取了元素 [test]。这是否意味着所有其他元素都被丢弃了？如果是这样，为什么还要花时间计算额外的元素？

Answer 2

您正在执行大量独立集成。您可以通过同时在不同的内核上执行这些集成来加快速度（如果您有可用的多核处理器）。问题在于 R 默认情况下以单线程方式执行其计算。但是，有许多可用的包允许多线程支持。我最近回答了几个类似的问题 and 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 函数。您可以尝试使用这些设置并通过允许不太严格的数值近似来获得性能提升。作为替代方案，您可以实现自己的自适应 Simpson 正交的简单版本，并使用离散步长。您很可能会像这样获得巨大的性能提升（如果您 able/willing 允许在近似值中出现更多误差）。

编辑： 更新代码以使其适用于所有场景：d and/or t valid/invalid 数字或向量

回复评论 @mawir：你是对的。 ifelse(test, yes, no) 将为测试计算为 TRUE 的行 return 对应的 yes 值，它将 return 相应的 no 行值test 的计算结果为 FALSE。但是，它首先必须评估您的 yes 表达式，以便创建 length(test) 的 yes 向量。这段代码演示了这一点：

> 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 和 d 的最后 5 个值与此场景相关。但是，在 F1 函数内部，ifelse 首先对所有 d 和 t 值计算 mapply 以创建 yes 向量。这就是函数执行需要这么长时间的原因。接下来，它选择满足条件的元素，否则为 0。 F4 函数解决了这个问题。

此外，你说你在t和d是非向量的情况下获得了加速。但是，在这种情况下，没有使用并行化。您通常应该在 senario 中获得最大加速，其中 t/d 之一或两者都是向量。

另一次编辑，以回应 Roland 的评论：如果您不想创建单独的函数文件，则可以将 clusterEvalQ(cl,eval(parse("functionG.R"))) 替换为 clusterExport(cl,"G")。

加速涉及映射和集成的功能

Speeding up a function involving mapply and integrate

optimization

r

integrate

mapply