c++ - 对超长 vector (长度为 ~1e+7 至 ~1e+8)的快速 pnorm() 计算

Question

有没有办法优化pnorm？我的代码遇到了一些瓶颈，经过大量优化和基准测试后，我意识到它来自对非常大的 vector 调用 pnorm。

使用 microbenchmarking 我在我的机器上发现如果 length(u) ~ 5e+7 那么 pnorm(u) 需要 11 秒。

有没有办法在这里使用Rcpp，或者内置的pnorm已经优化过了？

欢迎任何想法。

我在 SO 上找到了这些帖子:Use pnorm from Rmath.h with Rcpp和 How can I use qnorm on Rcpp? .但据我了解，它们的目的是将 R 函数用于 Cpp 代码。

Answer 1

在本节课中，我将展示对 pnorm() 的快速而准确的近似。

在我们开始之前，我们需要明确:我们想要通过使用近似来实现什么？效率/速度/性能，对吧？但这种效率从何而来？

如上所述，pnorm() 计算受内存限制；即使我们进行近似计算，它仍然是内存限制的(因此我们不考虑进一步的并行化)。内存限制问题有

number of floating point operations : memory access = O(1)

换句话说，这个比率是某个常数C。所以我们的目标应该是减少这个常量，即我们要减少浮点运算。

浮点运算的次数通常报告为浮点加法和乘法的次数。其他类型的浮点运算被“转换”为此类度量。下面，我们来比较一下几种常见的浮点运算的开销。

u <- sample(1:10/10, 5e+7, replace = TRUE)

system.time(u + u)
#  user  system elapsed 
# 0.468   0.168   0.639 
system.time(u * u)
#  user  system elapsed 
# 0.424   0.212   0.638 
system.time(u / u)
#  user  system elapsed 
# 0.504   0.204   0.710 
system.time(u ^ 1.1)
#  user  system elapsed 
# 7.240   0.212   7.458 
system.time(sqrt(u))
#  user  system elapsed 
# 2.044   0.176   2.224 
system.time(exp(u))
#  user  system elapsed 
# 4.336   0.208   4.550 
system.time(log(u))
#  user  system elapsed 
# 2.748   0.172   2.925 
system.time(round(u))
#  user  system elapsed 
# 6.836   0.188   7.034 

请注意，加法和乘法成本较低，根和对数成本较高，而某些运算成本非常高，包括幂、指数和舍入。

现在让我们回到 pnorm()，甚至 dnorm() 等，我们有一个指数项要计算。鉴于:

system.time(pnorm(u))
#   user  system elapsed 
# 11.016   0.160  11.193 
system.time(dnorm(u))
#  user  system elapsed 
# 8.844   0.164   9.022 

我们看到计算 pnorm() 和 dnorm() 的大部分时间都用于计算指数。 pnorm() 比 dnorm() 需要更长的时间，因为它还有一个积分!

现在，我们的目标相当明确:我们想用真正便宜的东西取代昂贵的 pnorm() 评估，理想情况下只涉及加法/乘法。我们可以吗？？

历史上出现过很多逼近方法。 @Ben 提到了逻辑近似。 R 函数 plogis() 执行此操作。但是仔细阅读 ?plogis 就会发现它仍然基于指数。

现在，我们可以使用非参数近似来代替那些参数近似吗？但是我们不应该在这里进行回归；相反，我们想使用一些高分辨率精确数据的插值函数来预测 pnorm()。好吧，线性插值是最佳选择，因为它非常高效(由于稀疏性:线性预测矩阵是三对角矩阵)。在 R 中，approx 执行此操作。我将不熟悉此内容的读者转至 ?approx，我将继续进行。

OP 说他只需要标准正态分布，所以我们专注于此。考虑以下近似函数(我没有使用 approxfun 因为我想要自定义 h):

approx.pnorm <- function(u, h = 0.2) {
  x <- seq(from = -4, to = 4, by = h)
  approx(x, pnorm(x), yleft = 0, yright = 1, xout = u)$y
  }

准确的数据取自 [-4, 4] 之间分辨率为 h 的网格。 -4以下的预测为0，4以上的预测为1，满足CDF的要求。给定新值 u，我们根据已知的准确数据通过线性插值来近似 pnorm(u)。

显然，分辨率 h 控制着准确性。考虑以下函数来计算 RMSE 并显示近似曲线:

RMSEh <- function(h) {
  x <- sort(rnorm(1000))
  y <- pnorm(x)
  y1 <- approx.pnorm(x, h)
  plot(x, y, type = "l", lwd = 2); lines(x, y1, col = 2, lwd = 2)
  mean((y - y1) ^ 2)^0.5
  }

par(mfrow = c(1, 3))
RMSEh(1)  # 0.01570339
RMSEh(0.5)  # 0.003968882
RMSEh(0.2)  # 0.000639888

实际上，当h = 0.2时，近似值已经相当不错了。所以我们将在下面使用h = 0.2。

---

基准测试

这应该是最精彩的部分了。在上面我们已经看到 pnorm(u) 的准确计算需要 11 秒。现在

system.time(approx.pnorm(u, h = 0.2))
#  user  system elapsed 
# 2.656   0.172   2.833 

哇，我们快了将近 4 倍!!