Il modo più veloce di cross-tabulare due vettori logici enormi in R

Question

Per due vettori logici, x e y, di lunghezza> 1E8, qual è il modo più veloce per calcolare le tabelle a croce 2x2?

Sospetto che la risposta sia scrivere in C/C++, ma mi chiedo se c'è qualcosa in R che è già abbastanza intelligente su questo problema, in quanto non è raro.

codice di esempio, per 300M voci (sentitevi liberi di lasciare N = 1E8 se 3E8 è troppo grande;. Ho scelto una dimensione totale poco meno di 2,5 GB (2.4GB) ho preso di mira una densità di 0,02, solo per renderlo più interessante (si potrebbe usare un vettore scarsa, se questo aiuta, ma il tipo di conversione può richiedere tempo)

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) 

Alcuni metodi evidenti:.

1. table
2. bigtabulate
3. Operazioni logiche semplici (ad es. sum(x & y))
4. Vector moltiplicazione (boo)
5. data.table
6. Alcuni di quanto sopra, con parallel dal pacchetto multicore (o il nuovo pacchetto parallel)

ho preso una pugnalata alla prima tre opzioni (vedi la mia risposta), ma sento che ci deve essere qualcosa di meglio e più veloce.

Trovo che table lavori molto lentamente. bigtabulate sembra eccessivo per una coppia di vettori logici. Infine, eseguire le operazioni logiche alla vaniglia sembra un kludge, e guarda ogni vettore troppe volte (3X? 7X?), Senza contare che riempie molta memoria aggiuntiva durante l'elaborazione, che è un enorme perdita di tempo.

La moltiplicazione di vettore è in genere una cattiva idea, ma quando il vettore è scarso, si può ottenere un vantaggio dall'archiviazione come tale e quindi dall'uso della moltiplicazione vettoriale.

Sentitevi liberi di variare N e p, se questo dimostrerà qualcosa di interessante comportamento delle funzioni di tabulazione. :)

---

Update 1. La mia prima risposta dà tempi su tre metodi ingenue, che è la base per credere table è lento. Tuttavia, la cosa fondamentale da comprendere è che il metodo "logico" è grossolanamente inefficiente. Guardare cosa sta facendo:

• 4 operazioni logiche vettore
• 4 conversioni di tipo (logici a intero o FP - per sum)
• 4 sommatorie vettore
• 8 assegnazioni (1 per l'operazione logica, 1 per la sommatoria)

Non solo, ma non è nemmeno compilato o parallelizzato. Eppure, batte ancora i pantaloni di table.Si noti che bigtabulate, con un tipo di conversione in più (1 * cbind...) batte ancora table.

Update 2. Per timore che qualcuno sottolineare che i vettori logici a sostegno R NA, e che questa sarà una chiave nel sistema per questi tabulati croce (che è vero nella maggior parte dei casi), vorrei sottolineare che i miei vettori vengono da is.na() o is.finite(). :) Ho debug di NA e altri valori non finiti - they've been a headache for me recently. Se non si sa se o non tutte le voci sono NA, si potrebbe verificare con any(is.na(yourVector)) - questo sarebbe saggio prima di adottare alcune delle idee che sorgono in questo Q & A.

---

Update 3. Brandon Bertelsen ha fatto una domanda molto ragionevole nei commenti: perché usare così tanti dati quando un sottocampione (l'insieme iniziale, dopo tutto, è un campione ;-)) potrebbe essere adeguato ai fini della creazione di un cross-tabulazione? Non andare troppo lontano nelle statistiche, ma i dati derivano da casi in cui le osservazioni TRUE sono molto rare, per entrambe le variabili. Uno è il risultato di un'anomalia di dati, l'altro a causa di un possibile bug nel codice (possibile errore, perché vediamo solo il risultato computazionale -. Pensare variabile x come "Garbage In", e y come "Garbage Out" come risultato , la questione è se i problemi in uscita causati dal codice sono gli unici casi in cui i dati sono anomali, o ci sono alcuni altri casi in cui buoni dati va male? (Questo è il motivo per cui ho chiesto a una domanda su stopping when a NaN, NA, or Inf is encountered.)

questo spiega anche il motivo per cui il mio esempio ha una bassa probabilità di TRUE valori, questi si verificano in realtà molto meno dello 0,1% del tempo

questo suggerisce un percorso soluzione diversa Sì: si suggerisce che possiamo usare due indici.?(Cioè le posizioni di TRUE in ciascun set) e il conteggio delle intersezioni. Evitai impostato intersezioni perché ero bruciato un po 'indietro da Matlab (sì, questo è R, ma portare con me), che prima ordinare elementi di un insieme prima che fa un incrocio. (I vagamente la complessità era ancora più imbarazzante:. Come O(n^2) anziché O(n log n))

Answer 1

qui risultati del metodo logico, table e bigtabulate, per N = 3E8:

  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 

In questo caso , table è un disastro.

Per confronto, qui è 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 

A questo punto, sembra che scrivere le proprie funzioni logiche è la cosa migliore, anche se che gli abusi sum, ed esamina ogni vettore logica più volte. Non ho ancora provato a compilare le funzioni, ma questo dovrebbe produrre risultati migliori.

Update 1 Se diamo bigtabulate valori che sono già interi, cioè se facciamo la conversione di tipo 1 * cbind(v1,v2) fuori di bigtabulate, quindi il multiplo N = 3E6 è 1,80, invece di 2.4. Il multiplo N = 3E8 relativo al metodo "logico" è solo 1.21, anziché 1.53.

---

Update 2

Come Joshua Ulrich ha fatto notare, la conversione in vettori di bit è un miglioramento significativo - stiamo ripartizione e muoversi molto meno di dati: vettori logiche di R consumano 4 byte per ogni voce ("Perché?", potresti chiedere ... Well, I don't know, but an answer may turn up here.), mentre un vettore bit consuma, beh, un bit, per voce, cioè 1/32 di dati. Quindi, x consuma 1,2e9 byte, mentre xb (la versione bit nel codice di seguito) consuma solo 3,75e7 byte.

Sono caduto table e le variazioni bigtabulate dai benchmark aggiornati (N = 3e8). Si noti che logicalB1 presuppone che i dati siano già un vettore bit, mentre logicalB2 è la stessa operazione con la penalità per la conversione del tipo. Poiché i miei vettori logici sono i risultati di operazioni su altri dati, non ho il vantaggio di iniziare con un vettore bit. Tuttavia, la sanzione da pagare è relativamente piccola. [La serie "logical3" esegue solo 3 operazioni logiche e quindi esegue una sottrazione. Dal momento che è tabulazione incrociata, sappiamo che il totale, come DWin ha osservato.]

 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 

ora abbiamo accelerato questo fino a prendere solo 1,8-2,8 secondi, anche con molte inefficienze lordo. C'è senza dubbio dovrebbe essere possibile farlo in meno di 1 secondo, con le modifiche che includono uno o più dei seguenti: codice C, compilazione e elaborazione multicore. Dopo tutte le 3 (o 4) operazioni logiche differenti potrebbero essere fatte indipendentemente, anche se questo è ancora uno spreco di cicli di calcolo.

Il più simile dei migliori sfidanti, logical3B2, è circa 80 volte più veloce di table. È circa 10 volte più veloce dell'ingenua operazione logica. E ha ancora molto spazio per migliorare.

---

Ecco il codice per produrre quanto sopra. NOTA Si consiglia di commentare alcune delle operazioni o dei vettori, a meno che non si disponga di molta RAM: la creazione di x, x1 e xb, insieme agli oggetti corrispondenti, occuperanno un bel po 'di memoria.

Inoltre, nota: avrei dovuto usare 1L come moltiplicatore intero per bigtabulate, invece di 1. A un certo punto, eseguirò nuovamente questa modifica e consiglierei di passare a chiunque utilizzi l'approccio bigtabulate.

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)} 
) 

Answer 2

Se stai facendo un sacco di operazioni su grandi vettori logiche, dare un'occhiata al pacchetto bit. Salva una tonnellata di memoria memorizzando i booleani come veri booleani a 1 bit.

Questo non aiuta con table; in realtà peggiora perché ci sono più valori unici nel vettore bit a causa di come è costruito. Ma lo in realtà aiuta con i confronti logici.

# 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 

Answer 3

Una tattica differente è considerare appena impostati intersezioni, utilizzando gli indici dei valori TRUE, approfittando che i campioni sono molto influenzati (cioè principalmente FALSE).

A tal fine, introdurre func_find01 e una traduzione che utilizza il pacchetto bit (func_find01B); tutto il codice che non appare nella risposta sopra è incollato sotto.

Ho eseguito di nuovo l'intera valutazione N = 3e8, tranne che ho dimenticato di utilizzare func_find01B; Riavviai i metodi più veloci contro di esso, in un secondo passaggio.

  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 

Solo i metodi "veloci":

 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 

Quindi, find01B è più veloce tra i metodi che non utilizzano vettori di bit pre-convertiti, con un margine sottile (2.099 secondi rispetto a 2.327 secondi). Da dove proviene lo find02? Successivamente ho scritto una versione che utilizza vettori bit pre-calcolati. Questo è ora il più veloce.

In generale, il tempo di esecuzione dell'approccio "metodo indici" può essere influenzato dalle probabilità marginali & congiunte. Sospetto che sarebbe particolarmente competitivo quando le probabilità sono ancora più basse, ma bisogna sapere che a priori, o attraverso un sottocampione.

---

Update 1. Ho anche cronometrato il suggerimento di Josh O'Brien, utilizzando tabulate() invece di table(). I risultati, a 12 secondi trascorsi, sono circa 2X find01 e circa la metà di bigtabulate2. Ora che i migliori metodi si stanno avvicinando 1 secondo, questo è anche relativamente lento:

user system elapsed 
7.670 5.140 12.815 

---

Codice:

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)} 
) 

Answer 4

Ecco una risposta con Rcpp, tabulazione solo quelle voci che non sono entrambi 0 . Sospetto che ci debbano essere diversi modi per migliorarlo, dato che questo è insolitamente lento; è anche il mio primo tentativo con Rcpp, quindi potrebbero esserci alcune evidenti inefficienze associate allo spostamento dei dati. Ho scritto un esempio di vaniglia volutamente esplicita, che dovrebbe consentire agli altri di dimostrare come questo possa essere migliorato.

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") 

risultati Timing per N = 3E8:

user system elapsed 
10.930 1.620 12.586 

questo richiede più di 6X finché func_find01B nella mia seconda risposta.

Answer 5

Ecco una risposta utilizzando zucchero 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) 
}) 

che si traduce in:

> 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 

Quindi, siamo in grado di ottenere un po 'di velocità fino oltre il pacchetto po', anche se sono sorpreso di come i tempi sono competitivi.

Aggiornamento: In onore di Iterator, ecco una soluzione iteratore 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