c++ - R 的 sum() 和 Armadillo 的 accu() 之间的区别

Question

在给定相同输入时，R 的 sum() 函数和 RcppArmadillo 的 accu() 函数的结果存在细微差别。例如以下代码:

R:

vec <- runif(100, 0, 0.00001)
accu(vec)
sum(vec)

C++:

// [[Rcpp::depends("RcppArmadillo")]]
// [[Rcpp::export]]
double accu(arma::vec& obj)
{
    return arma::accu(obj);
}

给出结果:

0.00047941851844312633 (C++)

0.00047941851844312628 (R)

根据http://keisan.casio.com/calculator真正的答案是:

4.79418518443126270948E-4

这些微小的差异在我的算法中加起来并显着影响它的执行方式。有没有办法更准确地总结 C++ 中的 vector ？或者至少获得与 R 相同的结果而无需调用 R 代码？

Answer 1

更新:根据其他人在源代码中找到的内容，我错了 - sum() 不排序。我在下面发现的一致性模式源于这样一个事实，即排序(如下面的某些情况下所做的那样)和使用扩展精度中间值(如 sum() 中所做的那样)会对精度产生类似的影响。 ..

@user2357112 评论如下:

src/main/summary.c ... doesn't do any sorting. (That'd be a lot of expense to add to a summation operation.) It's not even using pairwise or compensated summation; it just naively adds everything up left to right in an LDOUBLE (either long double or double, depending on HAVE_LONG_DOUBLE).

我已经筋疲力尽地在 R 源代码中寻找这个(没有成功 - sum 很难搜索)，但是 我可以通过实验证明当执行 sum ()，R将输入 vector 从小到大排序，以最大化准确率；下面的 sum() 和 Reduce() 结果之间的差异是由于使用了扩展精度。我不知道 accu 是做什么的...

 set.seed(101)
 vec <- runif(100, 0, 0.00001)
 options(digits=20)
 (s1 <- sum(vec))
 ## [1] 0.00052502325481269514554

使用 Reduce("+",...) 只是将元素按顺序添加。

 (s2 <- Reduce("+",sort(vec)))
 ## [1] 0.00052502325481269514554
 (s3 <- Reduce("+",vec))
 ## [1] 0.00052502325481269503712
 identical(s1,s2)  ## TRUE

?sum() 也说

Where possible extended-precision accumulators are used, but this is platform-dependent.

在 RcppArmadillo 中对已排序的 vector 执行此操作会给出与 R 中相同的答案；以原始顺序对 vector 执行此操作会给出不同的答案(我不知道为什么；我的猜测是前面提到的扩展精度累加器，当数据未排序时，它会更多地影响数值结果)。

suppressMessages(require(inline))
code <- '
   arma::vec ax = Rcpp::as<arma::vec>(x);
   return Rcpp::wrap(arma::accu(ax));
 '
## create the compiled function
armasum <- cxxfunction(signature(x="numeric"),
                        code,plugin="RcppArmadillo")
(s4 <- armasum(vec))
## [1] 0.00052502325481269525396
(s5 <- armasum(sort(vec)))
## [1] 0.00052502325481269514554
identical(s1,s5)  ## TRUE

但正如评论中指出的，这并不适用于所有种子:在这种情况下，Reduce() 结果更接近 sum 的结果()

set.seed(123)
vec2 <- runif(50000,0,0.000001)
s4 <- sum(vec2); s5 <- Reduce("+",sort(vec2))
s6 <- Reduce("+",vec2); s7 <- armasum(sort(vec2))
rbind(s4,s5,s6,s7)
##                       [,1]
## s4 0.024869900535651481843
## s5 0.024869900535651658785
## s6 0.024869900535651523477
## s7 0.024869900535651343065

我被难住了。我希望至少 s6 和 s7 是相同的...

我要指出的是，一般来说，当您的算法依赖于这些微小的数值差异时，您可能会感到非常沮丧，因为结果可能会因许多小的和可能超出的情况而有所不同-您的控制因素，例如您使用的特定操作系统、编译器等。