【问题标题】：Speeding up Rcpp `anyNA` equivalent加速 Rcpp `anyNA` 等效
【发布时间】：2021-03-08 15:23:03
【问题描述】：

这个问题与this old question和this old question有关。

R 具有很好的包装函数anyNA，可以更快地评估any(is.na(x))。在 Rcpp 中工作时，可以通过以下方式给出类似的最小实现：

// CharacterVector example
#include <Rcpp.h>
using namespace Rcpp;
template<typename T, typename S>
bool any_na(S x){
  T xx = as<T>(x);
  for(auto i : xx){
    if(T::is_na(i))
      return true;
  }
  return false;
}

// [[Rcpp::export(rng = false)]]
LogicalVector any_na(SEXP x){
  return any_na<CharacterVector>(x);
}

// [[Rcpp::export(rng = false)]]
SEXP overhead(SEXP x){
  CharacterVector xx = as<CharacterVector>(x);
  return wrap(xx);
}
/***R

library(microbenchmark)
vec <- sample(letters, 1e6, TRUE)
vec[1e6] <- NA_character_
any_na(vec)
# [1] TRUE
*/

但是将其与anyNA 的性能进行比较，我对下面的基准感到惊讶

library(microbenchmark)
microbenchmark(
  Rcpp = any_na(vec), 
  R = anyNA(vec),
  overhead = overhead(vec), 
  unit = "ms"
)
Unit: milliseconds
     expr      min        lq     mean    median       uq      max neval cld
     Rcpp 2.647901 2.8059500 3.243573 3.0435010 3.675051 5.899100   100   c
        R 0.800300 0.8151005 0.952301 0.8577015 0.961201 3.467402   100  b 
 overhead 0.001300 0.0029010 0.011388 0.0122510 0.015751 0.048401   100 a

最后一行是从SEXP 到CharacterVector 的来回转换所产生的“开销”（结果可以忽略不计）。显而易见，Rcpp 版本比 R 版本慢大约 3.5 倍。我很好奇，所以我检查了 Rcpp 的 is_na 的源代码，并没有发现性能缓慢的明显原因，我继续检查 source for anyNA for R's own character vectors's 并使用 R 的 C API 思维重新实现该函数以加快速度

// Added after SEXP overhead(SEXP x){ --- }
inline bool anyNA2(SEXP x){
  R_xlen_t n = Rf_length(x);
  for(R_xlen_t i = 0; i < n; i++){
    if(STRING_ELT(x, i) == NA_STRING)
      return true;
  }
  return false;
}
// [[Rcpp::export(rng = false)]]
SEXP any_na2(SEXP x){
  bool xx = anyNA2(x);
  return wrap(xx);
}
// [[Rcpp::export(rng = false)]]
SEXP any_na3(SEXP x){
  Function anyNA("anyNA");
  return anyNA(x);
}
/***R
microbenchmark(
  Rcpp = any_na(vec), 
  R = anyNA(vec),
  R_C_api = any_na2(vec),
  Rcpp_Function = any_na3(vec),
  overhead = overhead(vec),
  unit = "ms"
)
# Unit: milliseconds
# expr      min        lq       mean    median       uq      max neval  cld
# Rcpp 2.654901 2.8650515 3.54936501 3.2392510 3.997901 8.074201   100    d
# R 0.803701 0.8303015 1.01017200 0.9400015 1.061751 2.019902   100  b
# R_C_api 2.336402 2.4536510 3.01576302 2.7220010 3.314951 6.905101   100   c
# Rcpp_Function 0.844001 0.8862510 1.09259990 0.9597505 1.120701 3.011801   100  b
# overhead 0.001500 0.0071005 0.01459391 0.0146510 0.017651 0.101401   100 a
*/

请注意，我还包含一个简单的包装器，调用 anyNA 到 Rcpp::Function。再一次，anyNA 的这个实现不仅比基本实现慢一点，而且很多慢。

所以问题变成了2折：

为什么 Rcpp 这么慢？
源自 1：如何“更改”以加快代码速度？

这些问题本身并不是很有趣，但如果这会影响 Rcpp 实现的多个部分，而这些部分总体上可能会获得显着的性能提升，那就很有趣了。

SessonInfo()

sessionInfo()
R version 4.0.3 (2020-10-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)

Matrix products: default

locale:
[1] LC_COLLATE=English_Denmark.1252  LC_CTYPE=English_Denmark.1252    LC_MONETARY=English_Denmark.1252 LC_NUMERIC=C                     LC_TIME=English_Denmark.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] microbenchmark_1.4-7    cmdline.arguments_0.0.1 glue_1.4.2              R6_2.5.0                Rcpp_1.0.6             

loaded via a namespace (and not attached):
 [1] codetools_0.2-18 lattice_0.20-41  mvtnorm_1.1-1    zoo_1.8-8        MASS_7.3-53      grid_4.0.3       multcomp_1.4-15  Matrix_1.2-18    sandwich_3.0-0   splines_4.0.3   
[11] TH.data_1.0-10   tools_4.0.3      survival_3.2-7   compiler_4.0.3

编辑（不仅仅是windows问题）：

我想确保这不是“Windows 问题”，所以我在运行 linux 的 Docker 容器中检查并执行了该问题。结果如下图，非常相似

# Unit: milliseconds
#           expr    min      lq     mean  median      uq     max neval
#           Rcpp 2.3399 2.62155 4.093380 3.12495 3.92155 26.2088   100
#              R 0.7635 0.84415 1.459659 1.10350 1.42145 12.1148   100
#        R_C_api 2.3358 2.56500 3.833955 3.11075 3.65925 14.2267   100
#  Rcpp_Function 0.8163 0.96595 1.574403 1.27335 1.56730 11.9240   100
#       overhead 0.0009 0.00530 0.013330 0.01195 0.01660  0.0824   100

会话信息：

sessionInfo()
R version 4.0.2 (2020-06-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04 LTS

Matrix products: default
BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/openblas-openmp/libopenblasp-r0.3.8.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=C
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
 [9] LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] microbenchmark_1.4-7 Rcpp_1.0.5

loaded via a namespace (and not attached):
[1] compiler_4.0.2 tools_4.0.2

【问题讨论】：

（完成：你为什么要使用inline？当我们没有别的东西时它是一个非常有用的包（并且仍然可以正常工作）但是 Rcpp 属性使事情变得更加紧凑和可读。）跨度>
我不完全确定我在哪里使用inline@DirkEddelbuettel？我使用Rcpp::sourceCpp 使代码在这里可重现，而实际上代码被包装在一个包中（评估时间没有区别）
有道理。老实说，我没有考虑到这一点。为了读者，我会编辑。
另外，完全从臀部拍摄，你在你的第一个函数中的第一条语句强制一个转换（已经矢量化并在下面编译代码）R函数不会有.如果你把双腿绑在一起，在比赛中很难击败某人......
请@DirkEddelbuettel 按照我对您的评论和问题本身的回答中的链接进行操作。（他们直接转到 R 源代码中的doANY）。可以理解。我们都有有限的时间。 do_any2 类似于 R C 源中 anyNA 中 STRSXP 的第 2357 - 2360 行，由 do_anyNA 调用（do_anyNA 由 R 中的 anyNA 调用）。代码的其余部分处理错误检查和其他输入类型（例如VECSXP）。所以我在 R-source 中简单地消除了do_anyNA 和anyNA 的“不必要”部分，并在any_na2 中实现了它们。我在示例中找不到梨。

标签： r rcpp

【解决方案1】：

这是一个有趣的问题，但答案非常简单：STRING_ELT 有两个版本，一个由 R 内部使用，或者如果您在 Rinlinedfuns.h 中设置了 USE_RINTERNALS 宏，而在memory.c 中设置了一个用于 plebs 的版本.

对比两个版本，可以看到pleb版本的check比较多，这充分说明了速度上的差异。

如果你真的想要速度并且不关心安全性，你通常可以比 R 至少领先一点。

// [[Rcpp::export(rng = false)]]
bool any_na_unsafe(SEXP x) {
  SEXP* ptr = STRING_PTR(x);
  R_xlen_t n = Rf_xlength(x);
  for(R_xlen_t i=0; i<n; ++i) {
    if(ptr[i] == NA_STRING) return true;
  }
  return false;
}

替补：

> microbenchmark(
+   R = anyNA(vec),
+   R_C_api = any_na2(vec),
+   unsafe = any_na_unsafe(vec),
+   unit = "ms"
+ )
Unit: milliseconds
    expr    min      lq     mean  median      uq     max neval
       R 0.5058 0.52830 0.553696 0.54000 0.55465  0.7758   100
 R_C_api 1.9990 2.05170 2.214136 2.06695 2.10220 12.2183   100
  unsafe 0.3170 0.33135 0.369585 0.35270 0.37730  1.2856   100

虽然这样写是不安全的，但如果你在开始的循环之前添加一些检查就可以了。

【讨论】：

哇很好的答案。我试图找出是否有一些指针细节。我以为我已经检查了所有内容，但没想到我忽略了*_ELT 的第二个定义。我仍然很困惑为什么它如此不同（anyNA 在每次迭代中调用STRING_ELT）但这确实解释了时间差异。感谢您加倍努力。我打算加一个赏金。它似乎最终会得到你的答案。 :-)
啊，是的，脑残。 CHK 或 CHKZLN 会在每个值上调用，检查它们是否为 UNPROTECTed。有道理。
做得很好。顺便说一句，我建议删除unit="ms"，因为microbenchmark 最清楚要显示的内容，可以是nanoseconds。 unit="relative" 也可以发光。
我在想我们在 Rcpp 包中的（旧）基准测试示例：基本上没有什么可以忍受仅仅分配一个指针并增加它以节省（函数调用开销）与子集运算符。但是我忘记了这些宏就是这样做的，很高兴你把它放出来了。

【解决方案2】：

事实证明，这个问题很好地说明了为什么有些人对微基准测试进行抨击和咆哮。

基线是一个内置原语

这里应该被击败的函数实际上是一个原始的，所以它已经有点棘手了

> anyNA
function (x, recursive = FALSE)  .Primitive("anyNA")
>

ALTREP 降低性能水平

接下来，一个小实验表明基线函数anyNA() 永远不会循环。我们定义了一个很短的向量srt 和一个很长的向量lng，它们都包含一个NA 值。结果... R 通过 ALTREP 在数据结构标题中保留匹配位进行优化并且检查的成本与长度无关：

> srt <- c("A",NA_character_); lng <- c(rep("A", 1e6), NA_character_)
> microbenchmark(short=function(srt) { anyNA(srt) }, 
+                long=function(lng) { anyNA(lng) }, times=1000)
Unit: nanoseconds
  expr min lq   mean median uq   max neval cld
 short  48 50 69.324     51 53  5293  1000   a
  long  48 50 92.166     51 52 15494  1000   a
>

注意这里的单位（纳秒）和花费的时间。我们正在测量单个位。

（编辑：擦掉那个。匆忙想到我的，见 cmets。）

Rcpp 函数有一些小的开销

这不是新的和记录的。如果您查看 Rcpp Attributes 生成的代码，方便地给我们一个与我们指定的 C++ 函数同名的 R 函数，您会发现至少涉及一个其他函数调用。加上一个烘焙的try/catch层，RNG设置（这里关闭）等等。这不能为零，并且如果根据任何合理的东西摊销，这两者都没有出现在测量值中。

然而，这里的练习被设置为匹配一个查看一位的原始函数。这是一场谁也赢不了的比赛。所以这是我的决赛桌

> microbenchmark(anyNA = anyNA(vec), Rcpp_plain = rcpp_c_api(vec), 
+     Rcpp_tmpl = rcpp_any_na(vec), Rcpp_altrep = rcpp_altrep(vec), 
+     times = .... [TRUNCATED] 
Unit: microseconds
        expr      min      lq     mean   median      uq      max neval  cld
       anyNA  643.993  658.43  827.773  700.729  819.78  6280.85  5000 a   
  Rcpp_plain 1916.188 1952.55 2168.708 2022.017 2191.64  8506.71  5000    d
   Rcpp_tmpl 1709.380 1743.04 1933.043 1798.788 1947.83  8176.10  5000   c 
 Rcpp_altrep 1501.148 1533.88 1741.465 1590.572 1744.74 10584.93  5000  b

它包含原始的 R 函数，原始的（模板化的）C++ 函数，看起来还是不错的，使用 Rcpp 的东西（及其小开销）仅使用 C API（加上自动包装器输入/输出）有点慢 - - 然后比较 Michel 的 checkmate 包中的一个函数，确实查看了 ALTREP 位。而且速度几乎没有。

所以我们在这里真正看到的是函数调用的开销，它妨碍了测量微操作。所以不，Rcpp 不能比高度优化的原语更快。这个问题看起来很有趣，但归根结底，有点不合时宜。有时值得努力。

我的代码版本如下。

// CharacterVector example
#include <Rcpp.h>
using namespace Rcpp;
template<typename T, typename S>
bool any_na(S x){
    T xx = as<T>(x);
    for (auto i : xx){
        if (T::is_na(i))
            return true;
    }
    return false;
}

// [[Rcpp::export(rng = false)]]
LogicalVector rcpp_any_na(SEXP x){
    return any_na<CharacterVector>(x);
}

// [[Rcpp::export(rng = false)]]
SEXP overhead(SEXP x){
    CharacterVector xx = as<CharacterVector>(x);
    return wrap(xx);
}

// [[Rcpp::export(rng = false)]]
bool rcpp_c_api(SEXP x) {
    R_xlen_t n = Rf_length(x);
    for (R_xlen_t i = 0; i < n; i++) {
        if(STRING_ELT(x, i) == NA_STRING)
            return true;
  }
  return false;
}

// [[Rcpp::export(rng = false)]]
SEXP any_na3(SEXP x){
  Function anyNA("anyNA");
  return anyNA(x);
}

// courtesy of the checkmate package
// [[Rcpp::export(rng=false)]]
R_xlen_t rcpp_altrep(SEXP x) {
#if defined(R_VERSION) && R_VERSION >= R_Version(3, 5, 0)
    if (STRING_NO_NA(x))
        return 0;
#endif
    const R_xlen_t nx = Rf_xlength(x);
    for (R_xlen_t i = 0; i < nx; i++) {
        if (STRING_ELT(x, i) == NA_STRING)
            return i + 1;
    }
    return 0;
}


/***R
library(microbenchmark)

srt <- c("A",NA_character_)
lng <- c(rep("A", 1e6), NA_character_)
microbenchmark(short = function(srt) { anyNA(srt) },
               long = function(lng) { anyNA(lng) },
               times=1000)

N <- 1e6
vec <- sample(letters, N, TRUE)
vec[N] <- NA_character_
anyNA(vec)                      # to check


microbenchmark(
  anyNA       = anyNA(vec),
  Rcpp_plain  = rcpp_c_api(vec),
  Rcpp_tmpl   = rcpp_any_na(vec),
  Rcpp_altrep = rcpp_altrep(vec),
  #Rcpp_Function = any_na3(vec),
  #overhead = overhead(vec),
  times = 5000
#  unit="relative"
)
*/

【讨论】：

通过这个问题，我学习了微基准测试（用于 C 和 RCPP 函数），什么是原语（它们是用汇编编写的这么快吗？），ALTREP 表示替代表示（这句话让我调查）并且 C++ 并不那么令人生畏......（当然这并不容易）。感谢 Dirk 的精彩回答！
Hey Dirk，感谢您的回答（这是一本很棒的书！），我不会在 Windows R-4.0.3 上复制您的结果，所以我将进行更多调查（更具体地说：anyNA 在性能上不是静态的）。附带说明一下，我从来没有想过要击败 anyNA，而只是想发现为什么性能如此不同。此外，在这种情况下，Rcpp 属性所产生的开销非常小，以至于微不足道。
其实三思而后行，很容易找到。 microbenchmark 不会将 function 评估为调用，而只是将 function 评估为表达式。所以你的第一个基准测试只创建了 2 个函数 1e4 次，看起来好像没有“迭代”。（查看microbenchmark(anyNA(srt), anyNA(lng))）
哎呀。你就在那儿。我有一个较早的变体，带有 srt 和 lng 作为参数，我无法向自己解释为什么它们会有所不同。（SEXP 不包含有效载荷。）无论如何，我想现在已经完成了。我仍然认为您以错误的方式看待这个问题，并且正在想象一个没有的大问题。 Rcpp 函数永远不会击败 R 原语，尤其是在（本质上）无法测量的小操作上。所以结果是：性能并不像你声称的那样“如此不同”。它是纳秒。
再次感谢您抽出宝贵时间 Dirk。非常感谢，虽然我不完全同意你的结论。 “一个长 1500 行但最终调用相同的 3 行代码的函数如何比只有 3 行代码的函数更高效？” （可能存在夸大的行数）。我猜这是因为我缺乏将 C++ 代码优化到核心的知识（可能我也需要更深入地研究 R 源代码）。尽管如果我很幸运并且有人通过答案，这个问题仍然存在。 ;-)