计算具有大量参数组合的函数的最有效方法

Question

计算具有大量参数组合的函数的最有效方法

8

我尝试做的最简单示例：

dX_i <- rnorm(100, 0, 0.0002540362)

p_vec <- seq(0, 1, 0.25)  
gamma_vec <- seq(1, 2, 0.25)     
a_vec <- seq(2, 6, 1)
sigma_hat_vec <- c(0.03201636, 0.05771143, 0.07932116, 0.12262327, 0.15074560)
delta_j_vec <- c(0.0000005850109, 0.0000011700217, 0.0000017550326, 0.0000035100651, 0.0000052650977)

parameters <- expand.grid("p" = p_vec, "gamma" = gamma_vec, "a" = a_vec, "sigma_hat" = sigma_hat_vec, "delta_j" = delta_j_vec)


result <- sapply(1:nrow(parameters), function(x) {
  tmp <- parameters[x,]
  p <- tmp$p
  a <- tmp$a
  gamma <- tmp$gamma
  sigma_hat <- tmp$sigma_hat
  delta_j <- tmp$delta_j

  B <- sum( (abs(dX_i)^p) * ( abs(dX_i) < gamma * a * sigma_hat * delta_j^(1/2) ))

  return(B)
})

目标：根据所有组合的p、a、gamma、sigma_hat和delta_j，在向量dX上计算B。

但实际情况是，网格“parameters”有大约60万行，“dX_i”的长度为大约80k。此外，我有一个包含约1000个“dX_i”的列表。因此，我希望尽可能高效地进行这个计算。其他方法，例如将“parameters”转换为data.table并在该数据表中运行sapply，并没有给出速度提升。

我尝试了并行化函数（受限于在虚拟Windows机器上运行脚本）。

cl <- makePSOCKcluster(numCores)
num.iter <- 1:nrow(parameters)
parSapply(cl, num.iter, function(x, parameters, dX_i) {
  tmp <- parameters[x,]
  p <- tmp$p
  a <- tmp$a
  gamma <- tmp$gamma
  sigma_hat <- tmp$sigma_hat
  delta_j <- tmp$delta_j
  sum( (abs(dX_i)^p) * ( abs(dX_i) < gamma * a * sigma_hat * delta_j^(1/2) ))
}, parameters, dX_i)
stopCluster(cl)

虽然这给了我加速的效果，但我仍然觉得自己没有用最有效的方式解决这个问题，如果有建议，将不胜感激。

- zonfl

也许可以使用贝叶斯搜索？ - Bruno

2

只是好奇，它目前有多快，而“足够快”又有多快？ - Jon Spring

1

你有计算那么多项的充分理由吗？ - user1196549

2

你真的需要每一种组合吗？你的目标是什么？如果你正在寻找最小值或最大值，考虑使用优化器，它将能够实现比网格搜索更智能/更高效的搜索模式。例如，optim或optimx包。 - Gregor Thomas

1

@YalDan 这是回答Jon问题的好开端，但我认为他要求澄清的一个重要部分是你需要它快多少？2倍速度提升可以吗？10倍？100倍？ - Gregor Thomas

显示剩余3条评论

3个回答

10

当我想加速难以向量化的代码时，我经常转向Rcpp。最终，您正在尝试总结abs(dX_i)^p，限制在小于阈值gamma * a * sigma_hat * delta_j^(1/2)的abs(dX_i)值中。您希望对一堆p和阈值的配对执行此操作。您可以通过以下方式实现：

library(Rcpp)
cppFunction(
"NumericVector proc(NumericVector dX_i, NumericVector thresh, NumericVector p) {
  const int n = thresh.size();
  const int m = dX_i.size();
  NumericVector B(n);
  for (int i=0; i < n; ++i) {
    B[i] = 0;
    for (int j=0; j < m; ++j) {
      if (dX_i[j] < thresh[i]) {
        B[i] += pow(dX_i[j], p[i]);
      } else {
        break;
      }
    }
  }
  return B;
}"
)
result2 <- proc(sort(abs(dX_i)), parameters$gamma * parameters$a * parameters$sigma_hat * parameters$delta_j^(1/2), parameters$p)
all.equal(result, result2)
# [1] TRUE

请注意，我的代码对dX_i的绝对值进行排序，这样它就可以在遇到第一个超过阈值的值时停止计算。

在我的电脑上，我看到从您的代码到Rcpp代码的速度提升了20倍（使用system.time测量），时间从0.158秒降至0.007秒。

- josliber

非常感谢！加速效果显著，使用我的代码单次迭代需要几分钟的时间，而使用您的函数几乎可以立即得到结果。我没有想到会有如此快速和简单的解决方案，看来是时候学习一些C++了。 - zonfl

5

有一点观察结果是，您的参数集中实际上有每个p值的大量重复。您可以单独处理每个p值；这样，您只需要对特定p值计算一次dX_i的和。

result4 <- rep(NA, nrow(parameters))
sa_dX_i <- sort(abs(dX_i))
thresh <- parameters$gamma * parameters$a * parameters$sigma_hat * parameters$delta_j^(1/2)
loc <- findInterval(thresh, sa_dX_i)
loc[loc == 0] <- NA  # Handle threshold smaller than everything in dX_i
for (pval in unique(parameters$p)) {
  this.p <- parameters$p == pval
  cs_dX_i_p <- cumsum(sa_dX_i^pval)
  result4[this.p] <- cs_dX_i_p[loc[this.p]]
}
result4[is.na(result4)] <- 0  # Handle threshold smaller than everything in dX_i
all.equal(result, result4)
# [1] TRUE

为了看到这个效果，让我们将原始数据集扩大到问题所描述的规模（约有 600,000 行参数和约有 80,000 个值在 dX_i 中）：

set.seed(144)
dX_i <- rnorm(80000, 0, 0.0002540362)
p_vec <- seq(0, 1, 0.025)  
gamma_vec <- seq(1, 2, 0.025)     
a_vec <- seq(2, 6, 0.3)
sigma_hat_vec <- c(0.03201636, 0.05771143, 0.07932116, 0.12262327, 0.15074560)
delta_j_vec <- c(0.0000005850109, 0.0000011700217, 0.0000017550326, 0.0000035100651, 0.0000052650977)
parameters <- expand.grid("p" = p_vec, "gamma" = gamma_vec, "a" = a_vec, "sigma_hat" = sigma_hat_vec, "delta_j" = delta_j_vec)
dim(parameters)
# [1] 588350      5
length(unique(parameters$p))
# [1] 41

加速效果相当惊人——这段代码在我的计算机上只需要0.27秒，而我在回答这个问题时发布的与 Rcpp 相关的代码需要 655 秒（使用纯 R 实现的加速比达到 2400 倍！）。显然，这种加速方法只适用于 parameters 数据框中相对较少的重复p值。如果每个p值都是唯一的，则这种方法可能比其他提出的方法慢得多。

- josliber

1

这太不可思议了@josliber。我在我的大约9500万行数据集上运行了这段代码，速度从1天降至40分钟，最终仅需12秒！非常感谢您的帮助！顺便说一下，p通常不会超过seq(0, 6, 0.25)，所以这不应该是一个问题。 - zonfl

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- F. Privé · Accepted Answer

@josliber的回答非常好。然而，它让人觉得R不好...并且你必须转向C++来获得更好的性能。

他们的回答中实现了三个技巧:

预计算阈值向量
预计算dX_i的绝对值
对这些值进行排序以尽早停止求和

前两个技巧只是一种称为“向量化”的R技巧->基本上在整个向量上执行操作（例如gamma * a * sigma_hat * delta_j^(1/2)或abs()），而不是在循环中的单个元素上执行。

当使用sum(dX_i^p * vec_boolean)时，这正是您要做的；它被向量化了（*和sum），因此应该非常快。

如果我们仅实现这两个技巧（我们无法以相同的方式执行第三个技巧，因为它会破坏向量化），则结果为：

abs_dX_i <- abs(dX_i)
thresh <- with(parameters, gamma * a * sigma_hat * sqrt(delta_j))
p <- parameters$p
result3 <- sapply(1:nrow(parameters), function(i) {
  in_sum <- (abs_dX_i < thresh[i])
  sum(abs_dX_i[in_sum]^p[i])
})
all.equal(result, result3) # TRUE

如果我们对这三个解决方案进行基准测试：

microbenchmark::microbenchmark(
  OP = {
    result <- sapply(1:nrow(parameters), function(x) {
      tmp <- parameters[x,]
      p <- tmp$p
      a <- tmp$a
      gamma <- tmp$gamma
      sigma_hat <- tmp$sigma_hat
      delta_j <- tmp$delta_j

      B <- sum( (abs(dX_i)^p) * ( abs(dX_i) < gamma * a * sigma_hat * delta_j^(1/2) ))

      return(B)
    })
  },
  RCPP = {
    result2 <- proc(sort(abs(dX_i)), parameters$gamma * parameters$a *
                      parameters$sigma_hat * parameters$delta_j^(1/2), parameters$p)
  },
  R_VEC = {
    abs_dX_i <- abs(dX_i)
    thresh <- with(parameters, gamma * a * sigma_hat * sqrt(delta_j))
    p <- parameters$p
    result3 <- sapply(1:nrow(parameters), function(i) {
      in_sum <- (abs_dX_i < thresh[i])
      sum(abs_dX_i[in_sum]^p[i])
    })
  },
  times = 10
)

我们得到：

Unit: milliseconds
  expr      min       lq      mean   median       uq      max neval
    OP 224.8414 235.4075 289.90096 270.2767 347.1727 399.3262    10
  RCPP  14.8172  15.4691  18.83703  16.3979  20.3829  29.6624    10
 R_VEC  28.3136  29.5964  32.82456  31.4124  33.2542  45.8199    10

通过轻微修改 R 代码，它可以大幅提高速度。这比 Rcpp 代码慢不到两倍，并且可以像以前使用 parSapply() 一样轻松并行化。