Самый быстрый способ перекрестного табулирования двух массивных логических векторов в R

Для двух логических векторов x и y, длины > 1E8, что является самым быстрым способом вычисления кросс-таблиц 2x2?

Я подозреваю, что ответ заключается в том, чтобы написать его на C/С++, но мне интересно, есть ли что-то в R, которое уже довольно хорошо разбирается в этой проблеме, так как это не редкость.

Пример кода для записей 300M (не стесняйтесь N = 1E8, если 3E8 слишком велика, я выбрал общий размер чуть менее 2,5 ГБ (2,4 ГБ). Я нацелился на плотность 0,02, чтобы сделать ее более интересной (можно использовать разреженный вектор, если это помогает, но преобразование типа может занять время).

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)

Некоторые очевидные методы:

table
bigtabulate
Простые логические операции (например, sum(x & y))
Векторное умножение (boo)
data.table
Некоторые из вышеперечисленных, с parallel из пакета multicore (или нового пакета parallel)

Я взял удар в первых трех вариантах (см. мой ответ), но я чувствую, что должно быть что-то лучше и быстрее.

Я нахожу, что table работает очень медленно. bigtabulate кажется излишним для пары логических векторов. Наконец, выполнение ванильных логических операций похоже на kludge, и он смотрит на каждый вектор слишком много раз (3X? 7X?), Не говоря уже о том, что он заполняет много дополнительной памяти во время обработки, что является массовым отпуском времени.

Векторное умножение обычно является плохой идеей, но когда вектор разрежен, можно получить преимущество от его хранения как такового, а затем использовать векторное умножение.

Не изменяйте N и p, если это продемонстрирует что-либо интересное поведение функций табуляции.:)

Обновление 1. Мой первый ответ дает тайминги по трем наивным методам, который является основой для веры table медленным. Тем не менее, главное, чтобы понять, что "логический" метод крайне неэффективен. Посмотрите, что он делает:

4 операции логического вектора
4 преобразования типа (от логического до целого или FP - для 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. Брэндон Бертельсен задал очень разумный вопрос в комментариях: зачем использовать так много данных, когда подвыбор (исходный набор, в конце концов, является образцом;-)) может быть адекватным для целей создания кросс-табуляция? Не слишком сильно заходить в статистику, но данные возникают из случаев, когда наблюдения TRUE очень редки для обеих переменных. Один из них является результатом аномалии данных, другой - из-за возможной ошибки в коде (возможная ошибка, потому что мы видим только вычислительный результат - думаем о переменной x как "Garbage In" и y как "Garbage Out", В результате возникает вопрос, являются ли проблемы в выходе, вызванные кодом, исключительно теми случаями, когда данные являются аномальными, или есть некоторые другие случаи, когда хорошие данные ухудшаются? (Вот почему я задал вопрос о останавливается, когда встречаются NaN, NA или Inf.)

Это также объясняет, почему мой пример имеет небольшую вероятность для значений TRUE; они действительно происходят намного меньше, чем 0,1% времени.

Означает ли это другой путь решения? Да: это предполагает, что мы можем использовать два индекса (т.е. Местоположения TRUE в каждом наборе) и подсчетные пересечения. Я избегал заданных перекрестков, потому что я был сожжен некоторое время назад Matlab (да, это R, но медведь со мной), который сначала будет сортировать элементы набора до того, как он пересечет. (Я смутно вспоминаю, что сложность была еще более смущающей: например, O(n^2) вместо O(n log n).)

Ответ 1

Если вы делаете много операций над огромными логическими векторами, посмотрите на bit. Он сохраняет тонну памяти, сохраняя логические значения как истинные 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

Ответ 2

Вот результаты для логического метода, table и bigtabulate, для 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

В этом случае 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 и анализирует каждый логический вектор несколько раз. Я еще не пробовал компилировать функции, но это должно дать лучшие результаты.

Обновление 1 Если мы дадим значения bigtabulate, которые уже являются целыми числами, т.е. если мы делаем преобразование типа 1 * cbind(v1,v2) вне bigtabulate, то N = 3E6 кратно - 1,80 вместо 2.4, Множество N = 3E8 относительно "логического" метода составляет всего 1,21 вместо 1,53.

Обновление 2

Как отметил Джошуа Ульрих, преобразование в битовые векторы является значительным улучшением - мы выделяем и перемещаем LOT меньше данных: R логических векторов потребляют 4 байта на запись ( "Почему?", вы можете спросить.. Ну, я не знаю, но ответ может появиться здесь.), тогда как бит-вектор потребляет, ну, один бит, за запись - то есть 1/32 столько же данных. Итак, x потребляет 1,2e9 байт, а xb (версия бит в коде ниже) потребляет всего 3,75e7 байт.

Я сбросил варианты table и bigtabulate из обновленных тестов (N = 3e8). Обратите внимание, что logicalB1 предполагает, что данные уже являются битовым вектором, а logicalB2 - это одна и та же операция с штрафом за преобразование типа. Поскольку мои логические векторы являются результатом операций над другими данными, у меня нет возможности начать с битового вектора. Тем не менее, штраф, подлежащий выплате, относительно невелик. [Серия "logical3" выполняет только 3 логические операции и затем вычитает. Поскольку это кросс-табуляция, мы знаем общее число, как заметил DWin.]

        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 секунды, даже при большой неэффективности. Существует без сомнения, это должно быть осуществимо для этого в течение менее 1 секунды с изменениями, включая один или несколько из следующих: C-код, компиляция и многоядерная обработка. После того, как все 3 (или 4) различные логические операции могут быть выполнены независимо, хотя это все еще является пустой тратой вычислительных циклов.

Самый похожий из лучших претендентов, logical3B2, примерно на 80 раз быстрее, чем table. Это примерно в 10 раз быстрее, чем наивная логическая операция. И у него все еще есть много возможностей для улучшения.

Вот код для создания выше. ПРИМЕЧАНИЕ. Я рекомендую комментировать некоторые операции или векторы, если у вас много ОЗУ - создание x, x1 и xb, а также соответствующий y объектов, займет довольно немного памяти.

Также обратите внимание: я должен был использовать 1L как целочисленный множитель для bigtabulate, а не только 1. В какой-то момент я буду повторно запускаться с этим изменением и рекомендовал бы это изменение всем, кто использует подход 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)}
)

Ответ 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

Итак, мы можем немного ускорить работу над бит-пакетом, хотя я удивлен, насколько конкурентоспособны времена.

Обновление: в честь Iterator, это решение итератора 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

Ответ 4

Другая тактика состоит в том, чтобы рассматривать только заданные пересечения, используя индексы значений TRUE, пользуясь тем, что сэмплы очень смещены (т.е. в основном FALSE).

С этой целью я представляю func_find01 и перевод, который использует пакет bit (func_find01B); весь код, который не отображается в приведенном выше ответе, вставлен ниже.

Я повторил полную оценку N = 3e8, за исключением того, что забыл использовать func_find01B; Я повторил более быстрые методы против него, на втором проходе.

          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. Я также приурочил предложение Джоша О'Брайена, используя tabulate() вместо table(). Результаты, прошедшие через 12 секунд, составляют около 2X 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)}
)

Ответ 5

Здесь ответ с Rcpp, содержащий только те записи, которые не оба 0. Я подозреваю, что должно быть несколько способов улучшить это, поскольку это необычно медленно; это также моя первая попытка с Rcpp, поэтому могут быть некоторые очевидные неэффективности, связанные с перемещением данных. Я написал пример, который является целенаправленно простой ванилью, что должно позволить другим продемонстрировать, как это можно улучшить.

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

Результаты синхронизации для N = 3E8:

   user  system elapsed 
 10.930   1.620  12.586

В моем втором ответе это занимает более 6X, если func_find01B.