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

Posted 2023-02-16

【中文标题】R 的 sum() 和 Armadillo 的 accu() 之间的区别

【英文标题】：Difference between R's sum() and Armadillo's accu()

【发布时间】：2016-07-26 12:17:18

【问题描述】：

在给定相同输入时，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 代码？

【问题讨论】：

好吧，至少在这种特殊情况下，C++ 代码似乎比 R 代码更准确（假设在线工具的结果是可靠的黄金标准，我'不相信）所以你知道在这里使用 C++ 会对你的准确性产生负面影响吗？

You may find the R .Machine doc useful

我认为你应该通过在runif 之前调用set.seed 来提高你的示例的可重复性，以便其他人可以检查相同的测试用例。您还可以包含log 函数的代码吗？

伴随您的最后一次编辑，sum 实际上似乎将总和存储在定义的类型中 here

我认为您应该将您的编辑作为答案发布（如果它最好地解决您的问题，请接受它）而不是编辑您的问题

【参考方案1】：

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

@user2357112 cmets 下面：

src/main/summary.c ... 不进行任何排序。 （添加到求和操作中会花费很多。）它甚至没有使用成对或补偿求和；它只是天真地将所有内容从左到右添加到 LDOUBLE 中（long double 或 double，取决于 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() 也说

在可能的情况下使用扩展精度累加器，但这取决于平台。

在 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

但正如 cmets 中所指出的，这并不适用于所有种子：在这种情况下，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 是相同的......

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

【讨论】：

我尝试过排序和运行 accu，但是我没有得到和你一样的结果。示例：set.seed(123) vec

抱歉格式化。我不知道如何让代码在评论中看起来不错

你必须使用相同的随机数种子...尝试set.seed(101)而不是set.seed(123)

但如果它只适用于某些种子，那么我们怎么知道相似性不是偶然的呢？

我找到了来源。它在src/main/summary.c 中，它不做任何排序。 （添加到求和操作中会花费很多。）它甚至没有使用成对或补偿求和；它只是天真地将所有内容从左到右添加到 LDOUBLE（long double 或 double，取决于 HAVE_LONG_DOUBLE）。

【参考方案2】：

我发现了什么：

我成功地编写了一个能够模仿 R 的 sum 函数的函数。看起来 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() 慢得多

【讨论】：

对于基准测试，我通常不会打印出这么多数字。我在这里这样做的唯一原因是因为我之前可以选择这样做（用于测试我的问题的解决方案）

完成。我觉得你评论我用 Rcpp 标记的所有问题很酷，即使是唠叨，哈哈

请注意，sum 会进行一些检查以及明确的NA 处理，这就解释了为什么它更慢。

【参考方案3】：

您可以使用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

【讨论】：

我认为 OP 更喜欢在 C++ 中执行操作。幸运的是，C++ 还存在高精度算术库，尤其是 MPFR。

Konrad 关于我更喜欢​​ C++ 中的操作是正确的。这个库存在于 c++ 中，这很酷。我会检查一下。值得考虑的另一点是执行操作所需的时间。因此，如果我不需要对犰狳对象进行任何昂贵的转换操作，那就太好了