真正长向量的快速pnorm()计算(长度为~1e + 7到~1e + 8)

Question

有没有办法优化pnorm？我在我的代码中遇到了一些瓶颈,经过大量的优化和基准测试后,我意识到它来自pnorm对真正大的向量的调用.

随着microbenchmarking我的机器上了,如果length(u) ~ 5e+7再pnorm(u)需要11秒.

有没有办法在这里使用Rcpp,或内置pnorm已经优化？

欢迎任何想法.

我在SO上发现了这些帖子:使用Rcth.h中的pnorm和Rcpp,如何在Rcpp上使用qnorm？.但据我所知,他们的目的是将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
  }

准确的数据是在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倍!