R 的 sum() 和犰狳的 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++ 中的向量？或者至少无需调用 R 代码即可获得与 R 相同的结果？

Answer 1

您可以使用 mpfr 程序包（可靠的多精度浮点数）并指定小数点

 library("Rmpfr")
 set.seed(1)
 vec <- runif(100, 0, 0.00001)
#      [1] 2.655087e-06 3.721239e-06 5.728534e-06 9.082078e-06 2.016819e-06 8.983897e-06 9.446753e-06 6.607978e-06 6.291140e-06 6.178627e-07 2.059746e-06
#     [12] 1.765568e-06 6.870228e-06 3.841037e-06 7.698414e-06 4.976992e-06 7.176185e-06 9.919061e-06 3.800352e-06 7.774452e-06 9.347052e-06 2.121425e-06
#     [23] 6.516738e-06 1.255551e-06 2.672207e-06 3.861141e-06 1.339033e-07 3.823880e-06 8.696908e-06 3.403490e-06 4.820801e-06 5.995658e-06 4.935413e-06
#    [34] 1.862176e-06 8.273733e-06 6.684667e-06 7.942399e-06 1.079436e-06 7.237109e-06 4.112744e-06 8.209463e-06 6.470602e-06 7.829328e-06 5.530363e-06
#     [45] 5.297196e-06 7.893562e-06 2.333120e-07 4.772301e-06 7.323137e-06 6.927316e-06 4.776196e-06 8.612095e-06 4.380971e-06 2.447973e-06 7.067905e-07
#     [56] 9.946616e-07 3.162717e-06 5.186343e-06 6.620051e-06 4.068302e-06 9.128759e-06 2.936034e-06 4.590657e-06 3.323947e-06 6.508705e-06 2.580168e-06
#     [67] 4.785452e-06 7.663107e-06 8.424691e-07 8.753213e-06 3.390729e-06 8.394404e-06 3.466835e-06 3.337749e-06 4.763512e-06 8.921983e-06 8.643395e-06
#     [78] 3.899895e-06 7.773207e-06 9.606180e-06 4.346595e-06 7.125147e-06 3.999944e-06 3.253522e-06 7.570871e-06 2.026923e-06 7.111212e-06 1.216919e-06
#     [89] 2.454885e-06 1.433044e-06 2.396294e-06 5.893438e-07 6.422883e-06 8.762692e-06 7.789147e-06 7.973088e-06 4.552745e-06 4.100841e-06 8.108702e-06
#     [100] 6.049333e-06


sum(mpfr(vec,10))
#    1 'mpfr' number of precision  53   bits 
#    [1] 0.00051783234812319279

Answer 2

更新：根据其他人在源代码中找到的内容，我错了 - 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 将输入向量从小到大排序以最大化精度~~;以下 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 中对排序向量执行此操作给出与 R 中相同的答案；按原始顺序对向量执行此操作给出了不同的答案（我不知道为什么；我的猜测是前面提到的扩展精度累加器，当数据未排序时，它会更多地影响数值结果）。

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 是相同的...

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

Answer 3

我发现了什么：

我成功地编写了一个能够模仿 R 的求和函数的函数。看起来 R 使用更高精度的变量来存储每个加法运算的结果。

我写的：

// [[Rcpp::depends("RcppArmadillo")]]
// [[Rcpp::export]]
double accu2(arma::vec& obj)
{
    long double result = 0;
    for (auto iter = obj.begin(); iter != obj.end(); ++iter)
    {
        result += *iter;
    }
    return result;
}

速度比较：

set.seed(123)
vec <- runif(50000, 0, 0.000001)
microbenchmark(
  sum(vec),
  accu(vec),
  accu2(vec)
)


       expr    min     lq     mean  median      uq    max neval
   sum(vec) 72.155 72.351 72.61018 72.6755 72.7485 75.068   100
  accu(vec) 48.275 48.545 48.84046 48.7675 48.9975 52.128   100
 accu2(vec) 69.087 69.409 70.80095 69.6275 69.8275 182.955  100

因此，我的 C++ 解决方案仍然比 R 的 sum 快，但是它比犰狳的 accu()

慢得多

R 的 sum() 和犰狳的 accu() 之间的区别

Difference between R's sum() and Armadillo's accu()

c++

precision

r

armadillo

rcpp