15

長さが1E8を超える2つの論理ベクトル、xおよびの場合、2x2クロス集計を計算する最も速い方法は何ですか？y

答えはC/C ++で書くことだと思いますが、Rには、珍しいことではないので、この問題についてすでにかなり賢いものがあるのではないかと思います。

コード例、300Mエントリの場合（3E8が大きすぎる場合は、N = 1E8にしてください。合計サイズは2.5GB（2.4GB）未満を選択しました。密度を0.02に設定しました。これは、より面白くするためです（1つは可能です）。それが役立つ場合は、スパースベクトルを使用しますが、型変換には時間がかかる場合があります）。

set.seed(0)
N = 3E8
p = 0.02
x = sample(c(TRUE, FALSE), N, prob = c(p, 1-p), replace = TRUE)
y = sample(c(TRUE, FALSE), N, prob = c(p, 1-p), replace = TRUE)

いくつかの明白な方法：

  1. table
  2. bigtabulate
  3. 単純な論理演算（例sum(x & y)）
  4. ベクトル乗算（ブーイング）
  5. data.table
  6. 上記の一部parallel、multicoreパッケージ（または新しいparallelパッケージ）から

私は最初の3つのオプション（私の答えを参照）を突き刺しましたが、もっと良いものともっと速いものがあるに違いないと感じています。

table動作が非常に遅い ことがわかりました。bigtabulate論理ベクトルのペアにとってはやり過ぎのようです。最後に、バニラ論理演算を実行することは厄介なことのように見え、処理中に大量の追加メモリをいっぱいにすることは言うまでもなく、各ベクトルを何度も調べます（3X？7X？）。これは膨大な時間の浪費です。

ベクトル乗算は通常は悪い考えですが、ベクトルがスパースである場合は、それをそのまま保存してからベクトル乗算を使用することで利点が得られる場合があります。

自由に変更して、それが集計関数の興味深い動作を示す場合は、自由に変更Nしてください。p:)

table更新1.私の最初の答えは、3つの素朴な方法のタイミングを示しています。これは、遅いと信じる根拠です。しかし、理解すべき重要なことは、「論理的」な方法は非常に非効率的であるということです。それが何をしているのか見てください：

  • 4つの論理ベクトル演算
  • 4つの型変換（論理から整数またはFP-for sum）
  • 4つのベクトルの合計
  • 8つの割り当て（論理演算用に1つ、合計用に1つ）

それだけでなく、コンパイルも並列化もされていません。それでも、それはまだズボンを打ち負かしtableます。、追加の型変換（）を使用しても、がビートになることbigtabulateに注意してください。1 * cbind...table

更新2.Rの論理ベクトルがサポートNAしていること、およびこれらのクロス集計のシステムのレンチになることを誰もが指摘しないように（ほとんどの場合に当てはまります）、私のベクトルはis.na()またはから来ていることを指摘する必要がありis.finite()ます。:)私はデバッグNAやその他の非有限値を使用してきましたが、最近は頭痛の種になっています。すべてのエントリがであるかどうかわからない場合はNA、でテストすることができますany(is.na(yourVector))。これは、このQ＆Aで生じるアイデアのいくつかを採用する前に賢明です。

更新3.BrandonBertelsenはコメントで非常に合理的な質問をしました：サブサンプル（結局のところ、初期セットはサンプルです;-)）がクロスを作成する目的に十分であるのに、なぜこれほど多くのデータを使用するのですか？集計？統計にあまり深く入り込まないでください。ただし、データは、TRUE両方の変数について、観測が非常にまれな場合から発生します。1つはデータ異常の結果であり、もう1つはコードのバグの可能性によるものです（計算結果のみが表示されるため、バグの可能性があります。変数xを「ガベージイン」およびy「ガベージアウト」と考えてください。結果として、質問は、コードによって引き起こされる出力の問題が、データが異常である場合だけであるのか、それとも良いデータが悪くなる他の例があるのか​​ということです（これが私が質問した理由です、、、またはが検出されたときNaNに停止NAInfします。）

TRUEこれは、私の例の値の確率が低い理由も説明しています。これらは実際には0.1％未満の時間で発生します。

これは別のソリューションパスを示唆していますか？TRUEはい：2つのインデックス（つまり、各セット内のの位置）を使用して、セットの共通部分をカウントできることを示しています。交差点を実行する前にセットの要素を最初に並べ替えるMatlab（はい、これはRですが、我慢してください）によってしばらく前に燃やされたため、セットの交差点を避けました。（私は漠然と複雑さがさらに恥ずかしかったことを思い出します：O(n^2)代わりにのようにO(n log n)。）

4

5 に答える 5

11

巨大な論理ベクトルに対して多くの操作を行っている場合は、ビットパッケージを確認してください。ブール値を真の1ビットブール値として格納することにより、大量のメモリを節約します。

これは役に立ちませんtable; それがどのように構築されているかにより、ビットベクトルにはより多くの一意の値があるため、実際にはさらに悪化します。しかし、それは論理的な比較に本当に役立ちます。

# N <- 3e7
require(bit)
xb <- as.bit(x)
yb <- as.bit(y)
benchmark(replications = 1, order = "elapsed", 
    bit = {res <- func_logical(xb,yb)},
    logical = {res <- func_logical(x,y)}
)
#      test replications elapsed relative user.self sys.self user.child sys.child
# 1     bit            1   0.129  1.00000     0.132    0.000          0         0
# 2 logical            1   3.677 28.50388     2.684    0.928          0         0
于 2012-02-07T11:16:51.393 に答える
9

N = 3E8の場合の論理メソッド、、、の結果はtable次のとおりです。bigtabulate

         test replications elapsed relative user.self sys.self
2     logical            1  23.861 1.000000     15.36     8.50
3 bigtabulate            1  36.477 1.528729     28.04     8.43
1       table            1 184.652 7.738653    150.61    33.99

この場合、tableは災害です。

比較のために、ここにN=3E6があります。

         test replications elapsed relative user.self sys.self
2     logical            1   0.220 1.000000      0.14     0.08
3 bigtabulate            1   0.534 2.427273      0.45     0.08
1       table            1   1.956 8.890909      1.87     0.09

この時点では、を悪用しsum、各論理ベクトルを複数回調べますが、独自の論理関数を作成するのが最善のようです。関数のコンパイルはまだ試していませんが、より良い結果が得られるはずです。

更新1bigtabulateすでに整数である値を指定する場合、つまりbigtabulateの外部で型変換を行う場合、1 * cbind(v1,v2)N=3E6の倍数は2.4ではなく1.80になります。「論理的」方法に対するN=3E8の倍数は、1.53ではなく1.21にすぎません。

アップデート2

Joshua Ulrichが指摘しているように、ビットベクトルへの変換は大幅な改善です。つまり、データの割り当てと移動を大幅に減らしています。Rの論理ベクトルはエントリごとに4バイトを消費します（「なぜ？ 」わかりませんが、答えはここに表示される可能性があります。）一方、ビットベクトルは、エントリごとに1ビットを消費します。つまり、1/32のデータを消費します。したがって、x1.2e9バイトを消費しますが、xb（以下のコードのビットバージョン）は3.75e7バイトしか消費しません。

更新されたベンチマーク（N = 3e8）からのバリエーションを削除tableしました。データがすでにビットベクトルであると想定していることにbigtabulate注意してください。これは、型変換のペナルティを伴う同じ操作です。私の論理ベクトルは他のデータに対する演算の結果であるため、ビットベクトルから始める利点はありません。それにもかかわらず、支払われるペナルティは比較的小さいです。[「logical3」シリーズは、3つの論理演算のみを実行してから、減算を実行します。分割表なので、DWinが述べているように、合計がわかります。]logicalB1logicalB2

        test replications elapsed  relative user.self sys.self
4 logical3B1            1   1.276  1.000000      1.11     0.17
2  logicalB1            1   1.768  1.385580      1.56     0.21
5 logical3B2            1   2.297  1.800157      2.15     0.14
3  logicalB2            1   2.782  2.180251      2.53     0.26
1    logical            1  22.953 17.988245     15.14     7.82

多くの重大な非効率性があっても、これをわずか1.8〜2.8秒で高速化できました。Cコード、コンパイル、マルチコア処理の1つ以上を含む変更を加えれば、1秒以内にこれを実行できることは間違いありません。結局のところ、3つ（または4つ）の異なる論理演算は、それでも計算サイクルの無駄ですが、独立して実行できます。

最高の挑戦者の中で最も類似しているのは、logical3B2よりも約80倍高速ですtable。単純な論理演算よりも約10倍高速です。そして、それはまだ改善の余地がたくさんあります。

上記を生成するためのコードは次のとおりです。 注RAMが大量にある場合を除いて、操作またはベクトルの一部をコメントアウトすることをお勧めします。、、、および、を対応するオブジェクトとともに作成すると、xかなりのメモリが消費されます。x1xby

また、注意： 。だけでなく、1Lの整数乗数として使用する必要がありました。ある時点で、この変更を再実行し、このアプローチを使用するすべての人に変更をお勧めします。bigtabulate1bigtabulate

library(rbenchmark)
library(bigtabulate)
library(bit)

set.seed(0)
N <- 3E8
p <- 0.02

x <- sample(c(TRUE, FALSE), N, prob = c(p, 1-p), replace = TRUE)
y <- sample(c(TRUE, FALSE), N, prob = c(p, 1-p), replace = TRUE)

x1 <- 1*x
y1 <- 1*y

xb <- as.bit(x)
yb <- as.bit(y)

func_table  <- function(v1,v2){
    return(table(v1,v2))
}

func_logical  <- function(v1,v2){
    return(c(sum(v1 & v2), sum(v1 & !v2), sum(!v1 & v2), sum(!v1 & !v2)))
}

func_logicalB  <- function(v1,v2){
    v1B <- as.bit(v1)
    v2B <- as.bit(v2)
    return(c(sum(v1B & v2B), sum(v1B & !v2B), sum(!v1B & v2B), sum(!v1B & !v2B)))
}

func_bigtabulate    <- function(v1,v2){
    return(bigtabulate(1*cbind(v1,v2), ccols = c(1,2)))
}

func_bigtabulate2    <- function(v1,v2){
    return(bigtabulate(cbind(v1,v2), ccols = c(1,2)))
}

func_logical3   <- function(v1,v2){
    r1  <- sum(v1 & v2)
    r2  <- sum(v1 & !v2)
    r3  <- sum(!v1 & v2)
    r4  <- length(v1) - sum(c(r1, r2, r3))
    return(c(r1, r2, r3, r4))
}

func_logical3B   <- function(v1,v2){
    v1B <- as.bit(v1)
    v2B <- as.bit(v2)
    r1  <- sum(v1B & v2B)
    r2  <- sum(v1B & !v2B)
    r3  <- sum(!v1B & v2B)
    r4  <- length(v1) - sum(c(r1, r2, r3))
    return(c(r1, r2, r3, r4))
}

benchmark(replications = 1, order = "elapsed", 
    #table = {res <- func_table(x,y)},
    logical = {res <- func_logical(x,y)},
    logicalB1 = {res <- func_logical(xb,yb)},
    logicalB2 = {res <- func_logicalB(x,y)},

    logical3B1 = {res <- func_logical3(xb,yb)},
    logical3B2 = {res <- func_logical3B(x,y)}

    #bigtabulate = {res <- func_bigtabulate(x,y)},
    #bigtabulate2 = {res <- func_bigtabulate2(x1,y1)}
)
于 2012-02-07T06:05:40.090 に答える
3

これがRcpp砂糖を使った答えです。

N <- 1e8
x <- sample(c(T,F),N,replace=T)
y <- sample(c(T,F),N,replace=T)

func_logical  <- function(v1,v2){
    return(c(sum(v1 & v2), sum(v1 & !v2), sum(!v1 & v2), sum(!v1 & !v2)))
}


library(Rcpp)
library(inline)

doCrossTab1 <- cxxfunction(signature(x="integer", y = "integer"), body='
  Rcpp::LogicalVector Vx(x);
  Rcpp::LogicalVector Vy(y);
  Rcpp::IntegerVector V(4);

  V[0] = sum(Vx*Vy);
  V[1] = sum(Vx*!Vy);
  V[2] = sum(!Vx*Vy);
  V[3] = sum(!Vx*!Vy);
  return( wrap(V));
  '
, plugin="Rcpp")

system.time(doCrossTab1(x,y))

require(bit)
system.time(
{
xb <- as.bit(x)
yb <- as.bit(y)
func_logical(xb,yb)
})

その結果：

> system.time(doCrossTab1(x,y))
   user  system elapsed 
  1.067   0.002   1.069 
> system.time(
+ {
+ xb <- as.bit(x)
+ yb <- as.bit(y)
+ func_logical(xb,yb)
+ })
   user  system elapsed 
  1.451   0.001   1.453

ですから、時代の競争の激しさに驚いていますが、ビットパッケージよりも少しスピードアップすることができます。

更新：イテレーターに敬意を表して、ここにRcppイテレーターソリューションがあります：

doCrossTab2 <- cxxfunction(signature(x="integer", y = "integer"), body='
  Rcpp::LogicalVector Vx(x);
  Rcpp::LogicalVector Vy(y);
  Rcpp::IntegerVector V(4);
    V[0]=V[1]=V[2]=V[3]=0;
  LogicalVector::iterator itx = Vx.begin();
  LogicalVector::iterator ity = Vy.begin();
  while(itx!=Vx.end()){
    V[0] += (*itx)*(*ity);
    V[1] += (*itx)*(!*ity);
    V[2] += (!*itx)*(*ity);
    V[3] += (!*itx)*(!*ity);    
    itx++;
    ity++;
  }
  return( wrap(V));
  '
, plugin="Rcpp")

system.time(doCrossTab2(x,y))
#   user  system elapsed 
#  0.780   0.001   0.782
于 2012-02-14T03:12:45.590 に答える
2

TRUE別の戦術は、値のインデックスを使用して、サンプルが非常に偏っている (つまり、ほとんど ) ことを利用して、単にセットの交差を考慮することFALSEです。

そのために、パッケージ ( )func_find01を使用した翻訳を紹介します。上記の回答に表示されていないコードはすべて下に貼り付けられています。bitfunc_find01B

を使用するのを忘れたことを除いて、完全な N=3e8 評価を再実行しましたfunc_find01B。2回目のパスで、それに対してより高速な方法を再実行しました。

          test replications elapsed   relative user.self sys.self
6   logical3B1            1   1.298   1.000000      1.13     0.17
4    logicalB1            1   1.805   1.390601      1.57     0.23
7   logical3B2            1   2.317   1.785054      2.12     0.20
5    logicalB2            1   2.820   2.172573      2.53     0.29
2       find01            1   6.125   4.718798      4.24     1.88
9 bigtabulate2            1  22.823  17.583205     21.00     1.81
3      logical            1  23.800  18.335901     15.51     8.28
8  bigtabulate            1  27.674  21.320493     24.27     3.40
1        table            1 183.467 141.345917    149.01    34.41

「速い」方法だけ：

        test replications elapsed relative user.self sys.self
3     find02            1   1.078 1.000000      1.03     0.04
6 logical3B1            1   1.312 1.217069      1.18     0.13
4  logicalB1            1   1.797 1.666976      1.58     0.22
2    find01B            1   2.104 1.951763      2.03     0.08
7 logical3B2            1   2.319 2.151206      2.13     0.19
5  logicalB2            1   2.817 2.613173      2.50     0.31
1     find01            1   6.143 5.698516      4.21     1.93

そのため、find01B事前変換されたビット ベクトルを使用しない方法の中で、わずかな差 (2.099 秒対 2.327 秒) で最速です。どこfind02から来たの？その後、事前に計算されたビット ベクトルを使用するバージョンを作成しました。これが現在最速です。

一般に、「指数法」アプローチの実行時間は、限界確率と同時確率の影響を受ける可能性があります。確率がさらに低い場合は特に競争力があると思いますが、アプリオリに、またはサブサンプルを介して知っておく必要があります。

更新 1. Josh O'Brien の提案のタイミングを計り、tabulate()代わりにtable(). 12 秒経過した時点での結果は、約 2 倍find01で約半分ですbigtabulate2。最良の方法が 1 秒に近づいているため、これも比較的遅いです。

 user  system elapsed 
7.670   5.140  12.815

コード：

func_find01 <- function(v1, v2){
    ix1 <- which(v1 == TRUE)
    ix2 <- which(v2 == TRUE)

    len_ixJ <- sum(ix1 %in% ix2)
    len1    <- length(ix1)
    len2    <- length(ix2)
    return(c(len_ixJ, len1 - len_ixJ, len2 - len_ixJ,
             length(v1) - len1 - len2 + len_ixJ))
}

func_find01B <- function(v1, v2){
    v1b = as.bit(v1)
    v2b = as.bit(v2)

    len_ixJ <- sum(v1b & v2b)
    len1 <- sum(v1b)
    len2 <- sum(v2b)

    return(c(len_ixJ, len1 - len_ixJ, len2 - len_ixJ,
             length(v1) - len1 - len2 + len_ixJ))
}

func_find02 <- function(v1b, v2b){
    len_ixJ <- sum(v1b & v2b)
    len1 <- sum(v1b)
    len2 <- sum(v2b)

    return(c(len_ixJ, len1 - len_ixJ, len2 - len_ixJ,
             length(v1b) - len1 - len2 + len_ixJ))
}

func_bigtabulate2    <- function(v1,v2){
    return(bigtabulate(cbind(v1,v2), ccols = c(1,2)))
}

func_tabulate01 <- function(v1,v2){
    return(tabulate(1L + 1L*x + 2L*y))
}

benchmark(replications = 1, order = "elapsed", 
    table = {res <- func_table(x,y)},
    find01  = {res <- func_find01(x,y)},
    find01B  = {res <- func_find01B(x,y)},
    find02  = {res <- func_find01B(xb,yb)},
    logical = {res <- func_logical(x,y)},
    logicalB1 = {res <- func_logical(xb,yb)},
    logicalB2 = {res <- func_logicalB(x,y)},

    logical3B1 = {res <- func_logical3(xb,yb)},
    logical3B2 = {res <- func_logical3B(x,y)},

    tabulate    = {res <- func_tabulate(x,y)},
    bigtabulate = {res <- func_bigtabulate(x,y)},
    bigtabulate2 = {res <- func_bigtabulate2(x1,y1)}
)
于 2012-02-08T15:42:27.087 に答える
0

Here's an answer with Rcpp, tabulating only those entries that are not both 0. I suspect there must be several ways to improve this, as this is unusually slow; it's also my first attempt with Rcpp, so there may be some obvious inefficiencies associated with moving the data around. I wrote an example that is purposefully plain vanilla, which should let others demonstrate how this can be improved.

library(Rcpp)
library(inline)
doCrossTab <- cxxfunction(signature(x="integer", y = "integer"), body='
  Rcpp::IntegerVector Vx(x);
  Rcpp::IntegerVector Vy(y);
  Rcpp::IntegerVector V(3);
  for(int i = 0; i < Vx.length(); i++) {
    if( (Vx(i) == 1) & ( Vy(i) == 1) ){ V[0]++; } 
    else if( (Vx(i) == 1) & ( Vy(i) == 0) ){ V[1]++; } 
    else if( (Vx(i) == 0) & ( Vy(i) == 1) ){ V[2]++; } 
 }
  return( wrap(V));
  ', plugin="Rcpp")

Timing results for N = 3E8:

   user  system elapsed 
 10.930   1.620  12.586

This takes more than 6X as long as func_find01B in my 2nd answer.

于 2012-02-13T13:18:53.380 に答える