Fastest way to cross-tabulate two massive logical vectors in R

Question

For two logical vectors, x and y, of length > 1E8, what is the fastest way to calculate the 2x2 cross tabulations?

I suspect the answer is to write it in C/C++, but I wonder if there is something in R that is already quite smart about this problem, as it's not uncommon.

Example code, for 300M entries (feel free to let N = 1E8 if 3E8 is too big; I chose a total size just under 2.5GB (2.4GB). I targeted a density of 0.02, just to make it more interesting (one could use a sparse vector, if that helps, but type conversion can take time).

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)

Some obvious methods:

1. table
2. bigtabulate
3. Simple logical operations (e.g. sum(x & y))
4. Vector multiplication (boo)
5. data.table
6. Some of the above, with parallel from the multicore package (or the new parallel package)

I've taken a stab at the first three options (see my answer), but I feel that there must be something better and faster.

I find that table works very slowly. bigtabulate seems like overkill for a pair of logical vectors. Finally, doing the vanilla logical operations seems like a kludge, and it looks at each vector too many times (3X? 7X?), not to mention that it fills a lot of additional memory during processing, which is a massive time waster.

Vector multiplication is usually a bad idea, but when the vector is sparse, one may get an advantage out of storing it as such, and then using vector multiplication.

Feel free to vary N and p, if that will demonstrate anything interesting behavior of the tabulation functions. :)

---

Update 1. My first answer gives timings on three naive methods, which is the basis for believing table is slow. Yet, the key thing to realize is that the "logical" method is grossly inefficient. Look at what it's doing:

• 4 logical vector operations
• 4 type conversions (logical to integer or FP - for sum)
• 4 vector summations
• 8 assignments (1 for the logical operation, 1 for the summation)

Not only that, but it's not even compiled or parallelized. Yet, it still beats the pants off of table. Notice that bigtabulate, with an extra type conversion (1 * cbind...) still beats table.

Update 2. Lest anyone point out that logical vectors in R support NA, and that that will be a wrench in the system for these cross tabulations (which is true in most cases), I should point out that my vectors come from is.na() or is.finite(). :) I've been debugging NA and other non-finite values - they've been a headache for me recently. If you don't know whether or not all of your entries are NA, you could test with any(is.na(yourVector)) - this would be wise before you adopt some of the ideas arising in this Q&A.

---

Update 3. Brandon Bertelsen asked a very reasonable question in the comments: why use so much data when a sub-sample (the initial set, after all, is a sample ;-)) might be adequate for the purposes of creating a cross-tabulation? Not to drift too far into statistics, but the data arises from cases where the TRUE observations are very rare, for both variables. One is a result of a data anomaly, the other due to a possible bug in code (possible bug because we only see the computational result - think of variable x as "Garbage In", and y as "Garbage Out". AS a result, the question is whether the issues in the output caused by the code are solely those cases where the data is anomalous, or are there some other instances where good data goes bad? (This is why I asked a question about stopping when a NaN, NA, or Inf is encountered.)

That also explains why my example has a low probability for TRUE values; these really occur much less than 0.1% of the time.

Does this suggest a different solution path? Yes: it suggests that we may use two indices (i.e. the locations of TRUE in each set) and count set intersections. I avoided set intersections because I was burned awhile back by Matlab (yes, this is R, but bear with me), which would first sort elements of a set before it does an intersection. (I vaguely recall the complexity was even more embarrassing: like O(n^2) instead of O(n log n).)

Answer 1

If you're doing a lot of operations on huge logical vectors, take a look at the bit package. It saves a ton of memory by storing the booleans as true 1-bit booleans.

This doesn't help with table; it actually makes it worse because there are more unique values in the bit vector due to how it's constructed. But it really helps with logical comparisons.

# 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 2

A different tactic is to consider just set intersections, using the indices of the TRUE values, taking advantage that the samples are very biased (i.e. mostly FALSE).

To that end, I introduce func_find01 and a translation that uses the bit package (func_find01B); all of the code that doesn't appear in the answer above is pasted below.

I re-ran the full N=3e8 evaluation, except forgot to use func_find01B; I reran the faster methods against it, in a second pass.

          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

Just the "fast" methods:

        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

So, find01B is fastest among methods that do not use pre-converted bit vectors, by a slim margin (2.099 seconds versus 2.327 seconds). Where did find02 come from? I subsequently wrote a version that uses pre-computed bit vectors. This is now the fastest.

In general, the running time of the "indices method" approach may be affected by the marginal & joint probabilities. I suspect that it would be especially competitive when the probabilities are even lower, but one has to know that a priori, or via a sub-sample.

---

Update 1. I've also timed Josh O'Brien's suggestion, using tabulate() instead of table(). The results, at 12 seconds elapsed, are about 2X find01 and about half of bigtabulate2. Now that the best methods are approaching 1 second, this is also relatively slow:

 user  system elapsed 
7.670   5.140  12.815 

---

Code:

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 3

Here is an answer using Rcpp sugar.

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

which results 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 

So, we can get a little speed up over the bit package, though I'm surprised at how competitive the times are.

Update: In honor of Iterator, here is a Rcpp iterator solution:

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