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

Question

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

16

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

这些微小差异在我的算法中累加，并严重影响其执行方式。有没有一种更准确地对向量求和的方法？或者至少能不调用R代码就得到R相同的结果吗？

- Matthew Lueder

2

至少在这种特定情况下，C++代码似乎比R代码更准确（假设在线工具的结果是可靠的黄金标准，我并不完全相信），所以您是否确切知道在此处使用C++会对准确性产生负面影响？ - Konrad Rudolph

你可能会发现 R 的 .Machine 文档很有用：http://stat.ethz.ch/R-manual/R-patched/library/base/html/zMachine.html - patrickmdnet

假设在线工具是正确的，那么R代码更准确。第一个数字来自C ++，第二个数字来自R。最终，我的算法以任何方式都会收敛于相同的结果。不同之处在于达到目标所需的迭代次数。有时迭代次数的差异在于100个左右。 - Matthew Lueder

2

我认为你应该在调用runif之前使用set.seed来使你的示例更具可重复性，这样其他人可以根据相同的测试用例进行检查。你还可以包含你的log函数的代码吗？ - nrussell

3

为了配合你最后的编辑，sum 确实似乎将总和存储在此处定义的类型中这里。 - alexis_laz

1

我认为你应该将你的编辑作为答案发布（如果它是最好的解决方案，那么请接受它），而不是编辑你的问题。 - Ben Bolker

3个回答

10

我的发现：

我成功地编写了一个函数，可以模仿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()函数快，但它明显比armadillo的accu()函数慢。

- Matthew Lueder

为了基准测试，我通常不会打印这么多位数。我之所以在这里这样做，是因为我之前设置了这个选项（用于测试解决方案）。 - Matthew Lueder

做完了。我认为你评论我打上Rcpp标签的所有问题是很酷的，即使是唠叨，哈哈。 - Matthew Lueder

2

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

4

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

- ArunK

2

我认为OP更喜欢在C++中执行操作。幸运的是，高精度算术库，特别是MPFR，也存在于C++中。 - Konrad Rudolph

Konrad说得对，我更喜欢在C++中进行操作。很酷这个库在C++中存在。我会去看看的。另一个值得考虑的问题是执行操作所需的时间。因此，如果我不需要执行任何昂贵的转换操作，远离armadillo对象，那就太好了。 - Matthew Lueder

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Ben Bolker · Accepted Answer

更新: 根据其他人在源代码中发现的内容，我之前的说法是错误的 - sum() 不会排序。下面我发现的一致性模式源于排序（如下面某些情况所做的）和使用扩展精度中间值（如sum()所做的）对精度具有相似的影响...

@user2357112 在下面进行了评论:

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

但是评论中指出，这种方法并不适用于所有的seed：在这种情况下，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是相同的。

值得注意的是，通常情况下，当您的算法依赖于这些微小的数字差异时，您很可能会感到非常沮丧，因为结果很可能会因许多小的、可能超出您控制范围的因素（如特定操作系统、编译器等）而有所不同。