Разница между суммой R() и accu() Армадилло

обновление: судя по тому, что другие нашли в исходном коде, я ошибался — 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

Но, как указано в комментариях, это работает не для всех семян: в этом случае результат 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 будут идентичны...

Я укажу, что в целом, когда ваш алгоритм зависит от такого рода крошечных числовых различий, вы, вероятно, будете очень разочарованы, так как результаты, вероятно, будут отличаться на основе множества мелких и, возможно, не в ваших силах. управляющие факторы, такие как конкретная операционная система, компилятор и т. д., с которыми вы работаете.

Что я нашел:

Мне удалось написать функцию, которая может имитировать функцию суммы 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

Итак, мое решение на С++ по-прежнему быстрее, чем сумма R, однако оно значительно медленнее, чем accu() броненосца.

вы можете использовать пакет 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