Rcppの欠損値の高速検査

Question

この質問はNA values in Rcpp conditionalにリンクしています。

私は基本的に、複数の（ダブル）要素をループするRcppコードをいくつか持っています。そして、要素ごとに欠損値があるかどうかを確認する必要があります（ベクトル化は使用できません）。のは、ベクターに欠損値の数をカウントしてみましょう、同じように最小限の再現性の例：

#include <Rcpp.h> 
using namespace Rcpp; 

// [[Rcpp::export]] 
int nb_na(const NumericVector& x) { 
    int n = x.size(); 
    int c = 0; 
    for (int i = 0; i < n; i++) if (R_IsNA(x[i])) c++; 
    return c; 
} 

// [[Rcpp::export]] 
int nb_na3(const NumericVector& x) { 
    int n = x.size(); 
    int c = 0; 
    for (int i = 0; i < n; i++) if (x[i] == 3) c++; 
    return c; 
} 

// [[Rcpp::export]] 
LogicalVector na_real(NumericVector x) { 
    return x == NA_REAL; 
} 

その後、Rで、我々が得る：リンク質問で述べたように

> x <- rep(c(1, 2, NA), 1e4) 

> x2 <- replace(x, is.na(x), 3) 

> microbenchmark::microbenchmark(
+ nb_na(x), 
+ nb_na3(x2) 
+) 
Unit: microseconds 
     expr  min  lq  mean median  uq  max neval 
    nb_na(x) 135.633 135.982 153.08586 139.753 140.3115 1294.928 100 
nb_na3(x2) 22.490 22.908 30.14005 23.188 23.5025 684.026 100 

> all.equal(nb_na(x), nb_na3(x2)) 
[1] TRUE 

> na_real(x[1:3]) 
[1] NA NA NA 

、あなたがすることができませんx[i] == NA_REALにチェックすると、常に欠損値が返されるためです。ただし、R_IsNA(x[i])を使用すると数値との等価性を確認する速度がはるかに遅くなります（例：3）。

基本的には、単一の値が欠損値であることを確認できる解決策が必要です。この解決法は、数値と等しいかどうかをチェックするほど速くなければなりません。

Answer 1

浮動小数点値が欠落しているかどうかを確認する作業は高価で、周囲には道がありません。これは、Rとは無関係です。

これは、Rの今後のALTREPに興奮している理由の1つです。たとえば、二重/実ベクトルに欠損値が含まれているかどうかを追跡することができます。時間を無駄にする必要はありません。 ALTREPについて言及するために更新されていませんが、https://github.com/HenrikBengtsson/Wishlist-for-R/issues/12

Answer 2

の欠落値をチェックするか、NaN固有のバリアントをチェックすることで、特定の値をチェックするよりも常に高価になります。それは単なる浮動小数点演算です。

ただし、コードにまだ改善の余地があります。 R_IsNAの代わりにNumericVector::is_naを使用することをお勧めしますが、これは主に化粧品です。

分岐が高価な場合があります。つまり、if (R_IsNA(x[i])) c++;をc += NumericVector::is_na(x[i])に置き換えます。

// [[Rcpp::export]] 
int nb_na4(const NumericVector& x) { 
    int n = x.size(); 
    int c = 0; 
    for (int i = 0; i < n; i++) c += NumericVector::is_na(x[i]) ; 
    return c; 
} 

が続いintに反復してx[i]はstd::count_if algorithmを使用することによって置き換えることができますアクセス：これは、このバージョンを提供します。これはそれが理由だ。このバージョンにつながる：今、パフォーマンスはまだ十分ではない場合

// [[Rcpp::export]] 
int nb_na5(const NumericVector& x) { 
    return std::count_if(x.begin(), x.end(), NumericVector::is_na) ; 
} 

、私は一般的にRcppParallelパッケージからtbbライブラリを使用し、このため、並列化をしようとする場合があります。

// [[Rcpp::export]] 
int nb_na6(const NumericVector& x) { 
    return tbb::parallel_reduce( 
    tbb::blocked_range<const double*>(x.begin(), x.end()), 
    0, 
    [](const tbb::blocked_range<const double*>& r, int init) -> int { 
     return init + std::count_if(r.begin(), r.end(), NumericVector::is_na); 
    }, 
    [](int x, int y){ return x+y; } 
) ; 
} 

この機能を備えたベンチマーキング：私のマシン上で

library(microbenchmark) 

bench <- function(n){ 
    x <- rep(c(1, 2, NA), n) 
    microbenchmark(
    nb_na = nb_na(x), 
    nb_na4 = nb_na4(x), 
    nb_na5 = nb_na5(x), 
    nb_na6 = nb_na6(x) 
) 
} 
bench(1e5) 

私が取得：

> bench(1e4) 
Unit: microseconds 
expr min  lq  mean median  uq  max neval cld 
nb_na 84.358 94.6500 107.41957 110.482 118.9580 137.393 100 d 
nb_na4 59.984 69.4925 79.42195 82.442 85.9175 106.567 100 b 
nb_na5 65.047 75.2625 85.17134 87.501 93.0315 116.993 100 c 
nb_na6 39.205 51.0785 59.20582 54.457 68.9625 97.225 100 a 

> bench(1e5) 
Unit: microseconds 
expr  min  lq  mean median  uq  max neval cld 
nb_na 730.416 732.2660 829.8440 797.4350 872.3335 1410.467 100 d 
nb_na4 520.800 521.6215 598.8783 562.7200 657.1755 1059.991 100 b 
nb_na5 578.527 579.3805 664.8795 626.5530 710.5925 1166.365 100 c 
nb_na6 294.486 345.2050 368.6664 353.6945 372.6205 897.552 100 a 

もう一つの方法は完全に浮動小数点演算を回避し、ベクトルがlong longのベクトルであるふりをすることです、別名64ビットの整数値をビットパターンNA_REALと比較する：

このハックを使用して0

> devtools::install_github("ThinkR-open/seven31") 
    > seven31::reveal(NA, NaN, +Inf, -Inf) 
    0 11111111111 (NaN) 0000000000000000000000000000000000000000011110100010 : NA 
    0 11111111111 (NaN) 1000000000000000000000000000000000000000000000000000 : NaN 
    0 11111111111 (NaN) 0000000000000000000000000000000000000000000000000000 : +Inf 
    1 11111111111 (NaN) 0000000000000000000000000000000000000000000000000000 : -Inf 

シリアル溶液：

// [[Rcpp::export]] 
int nb_na7(const NumericVector& x){ 
    const long long* p = reinterpret_cast<const long long*>(x.begin()) ; 
    long long na = *reinterpret_cast<long long*>(&NA_REAL) ; 

    return std::count(p, p + x.size(), na) ; 

} 

そして平行バージョン：

// [[Rcpp::export]] 
int nb_na8(const NumericVector& x){ 
    const long long* p = reinterpret_cast<const long long*>(x.begin()) ; 
    long long na = *reinterpret_cast<long long*>(&NA_REAL) ; 

    auto count_chunk = [=](const tbb::blocked_range<const long long*>& r, int init) -> int { 
    return init + std::count(r.begin(), r.end(), na); 
    } ; 

    return tbb::parallel_reduce( 
    tbb::blocked_range<const long long*>(p, p + x.size()), 
    0, 
    count_chunk, 
    [](int x, int y){ return x+y; } 
) ; 

} 

    > bench(1e5) 
    Unit: microseconds 
    expr  min  lq  mean median  uq  max neval cld 
    nb_na 730.346 762.5720 839.9479 857.5865 881.8635 1045.048 100  f 
    nb_na4 520.946 521.6850 589.0911 578.2825 653.4950 832.449 100 d 
    nb_na5 578.621 579.3245 640.9772 616.8645 701.8125 890.736 100  e 
    nb_na6 291.115 307.4300 340.1626 344.7955 360.7030 484.261 100 c 
    nb_na7 122.156 123.4990 141.1954 132.6385 149.7895 253.988 100 b  
    nb_na8 69.356 86.9980 109.6427 115.2865 126.2775 182.184 100 a  

    > bench(1e6) 
    Unit: microseconds 
    expr  min  lq  mean median  uq  max neval cld 
    nb_na 7342.984 7956.3375 10261.583 9227.7450 10869.605 79757.09 100 d 
    nb_na4 5286.970 5721.9150 7659.009 6660.2390 9234.646 31141.47 100 c 
    nb_na5 5840.946 6272.7050 7307.055 6883.2430 8205.117 10420.48 100 c 
    nb_na6 2833.378 2895.7160 3891.745 3049.4160 4054.022 18242.26 100 b 
    nb_na7 1661.421 1791.1085 2708.992 1916.6055 2232.720 60827.63 100 ab 
    nb_na8 650.639 869.6685 1289.373 939.0045 1291.025 10223.29 100 a 

これは、NAを表すために1ビットのみのパターンがあると想定しています。

は、ここに参考のために私の全体のファイルです：

#include <Rcpp.h> 
#include <RcppParallel.h> 

// [[Rcpp::depends(RcppParallel)]] 
// [[Rcpp::plugins(cpp11)]] 
using namespace Rcpp; 

// [[Rcpp::export]] 
int nb_na(const NumericVector& x) { 
    int n = x.size(); 
    int c = 0; 
    for (int i = 0; i < n; i++) if (R_IsNA(x[i])) c++; 
    return c; 
} 

// [[Rcpp::export]] 
int nb_na4(const NumericVector& x) { 
    int n = x.size(); 
    int c = 0; 
    for (int i = 0; i < n; i++) c += NumericVector::is_na(x[i]) ; 
    return c; 
} 

// [[Rcpp::export]] 
int nb_na5(const NumericVector& x) { 
    return std::count_if(x.begin(), x.end(), NumericVector::is_na) ; 
} 

// [[Rcpp::export]] 
int nb_na6(const NumericVector& x) { 
    return tbb::parallel_reduce( 
    tbb::blocked_range<const double*>(x.begin(), x.end()), 
    0, 
    [](const tbb::blocked_range<const double*>& r, int init) -> int { 
     return init + std::count_if(r.begin(), r.end(), NumericVector::is_na); 
    }, 
    [](int x, int y){ return x+y; } 
) ; 
} 

// [[Rcpp::export]] 
int nb_na7(const NumericVector& x){ 
    const long long* p = reinterpret_cast<const long long*>(x.begin()) ; 
    long long na = *reinterpret_cast<long long*>(&NA_REAL) ; 

    return std::count(p, p + x.size(), na) ; 

} 

// [[Rcpp::export]] 
int nb_na8(const NumericVector& x){ 
    const long long* p = reinterpret_cast<const long long*>(x.begin()) ; 
    long long na = *reinterpret_cast<long long*>(&NA_REAL) ; 

    auto count_chunk = [=](const tbb::blocked_range<const long long*>& r, int init) -> int { 
    return init + std::count(r.begin(), r.end(), na); 
    } ; 

    return tbb::parallel_reduce( 
    tbb::blocked_range<const long long*>(p, p + x.size()), 
    0, 
    count_chunk, 
    [](int x, int y){ return x+y; } 
) ; 

} 

/*** R 
library(microbenchmark) 

bench <- function(n){ 
    x <- rep(c(1, 2, NA), n) 
    microbenchmark(
    nb_na = nb_na(x), 
    nb_na4 = nb_na4(x), 
    nb_na5 = nb_na5(x), 
    nb_na6 = nb_na6(x), 
    nb_na7 = nb_na7(x), 
    nb_na8 = nb_na8(x) 
) 
} 
bench(1e5) 
bench(1e6) 
*/