c++ - 在很长的 vector 上(长度从〜1e + 7到〜1e + 8)快速进行pnorm()计算

有没有一种方法可以优化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代码。

最佳答案

在本节中，我将演示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
  }

精确数据在h之间的分辨率[-4, 4]的网格上获取。 -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倍!