Cách nhanh nhất để đặt chéo hai vectơ lôgic lớn trong R

Question

Đối với hai vectơ hợp lý, x và y, dài> 1E8, cách nhanh nhất để tính toán các bảng chéo 2x2 là gì?

Tôi nghi ngờ câu trả lời là viết nó trong C/C++, nhưng tôi tự hỏi nếu có điều gì đó trong R đã khá thông minh về vấn đề này, vì nó không phải là không phổ biến.

Ví dụ mã, cho 300M mục (cảm thấy tự do để N = 1E8 nếu 3E8 quá lớn; tôi đã chọn tổng kích thước dưới 2,5 GB (2,4 GB) .Tôi nhắm mục tiêu mật độ 0,02, chỉ để làm cho nó nhiều hơn thú vị (người ta có thể sử dụng một vector thưa thớt, nếu điều đó giúp, nhưng kiểu chuyển đổi có thể mất thời gian)

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) 

một số phương pháp rõ ràng:.

1. table
2. bigtabulate
3. Hoạt động lôgic đơn giản (ví dụ: sum(x & y))
4. Vector nhân (boo)
5. data.table
6. Một số những điều trên, với parallel từ gói multicore (hoặc parallel gói mới)

Tôi đã lấy một đâm vào người đầu tiên ba lựa chọn (xem câu trả lời của tôi), nhưng tôi cảm thấy rằng phải có một cái gì đó tốt hơn và nhanh hơn.

Tôi thấy rằng table hoạt động rất chậm. bigtabulate có vẻ như quá mức cần thiết cho một cặp vectơ logic. Cuối cùng, thực hiện các hoạt động hợp lý vani có vẻ giống như một kludge, và nó nhìn vào mỗi vectơ quá nhiều lần (3X? 7X?), Chưa kể rằng nó lấp đầy nhiều bộ nhớ bổ sung trong quá trình xử lý, đó là một thời gian thải lớn.

Phép nhân Vector thường là một ý tưởng tồi, nhưng khi véc tơ thưa thớt, người ta có thể có được lợi thế từ việc lưu trữ nó như vậy, và sau đó sử dụng phép nhân vectơ.

Hãy thay đổi N và p, nếu điều đó sẽ chứng minh bất kỳ hành vi thú vị nào của các chức năng tabulation. :)

---

Update 1. Câu trả lời đầu tiên của tôi cho timings trên ba phương pháp ngây thơ, đó là cơ sở để tin table là chậm. Tuy nhiên, điều quan trọng để nhận ra là phương pháp "hợp lý" là không hiệu quả. Nhìn vào những gì nó làm:

• 4 hoạt động vector logic
• 4 chuyển đổi loại hình (logic để nguyên hoặc FP - cho sum)
• 4 summations vector
• 8 bài tập (1 cho các hoạt động hợp lý, 1 để tổng kết)

Không chỉ vậy, nhưng nó thậm chí không được biên soạn hoặc song song. Tuy nhiên, nó vẫn đánh bại chiếc quần của table.Lưu ý rằng bigtabulate, với chuyển đổi loại bổ sung (1 * cbind...) vẫn bị đánh bại table.

Cập nhật 2. Vì vậy, bất kỳ ai cũng chỉ ra rằng các vectơ logic trong hỗ trợ R NA và đó sẽ là một cờ lê trong hệ thống cho các bảng chéo này (điều này đúng trong hầu hết các trường hợp). từ is.na() hoặc is.finite(). :) Tôi đã gỡ lỗi NA và các giá trị không hữu hạn khác - they've been a headache for me recently. Nếu bạn không biết có hay không phải tất cả các mục của bạn là NA, bạn có thể thử nghiệm với any(is.na(yourVector)) - đây sẽ là khôn ngoan trước khi bạn áp dụng một số ý tưởng phát sinh trong Q này & A.

---

Cập nhật 3. Brandon Bertelsen hỏi một câu hỏi rất hợp lý trong các ý kiến: tại sao lại sử dụng quá nhiều dữ liệu khi một mẫu phụ (tập đầu tiên, sau khi tất cả, là một mẫu ;-)) có thể là thích hợp cho mục đích tạo một bảng chéo? Không được trôi quá xa vào số liệu thống kê, nhưng dữ liệu phát sinh từ các trường hợp mà các quan sát TRUE là rất hiếm, đối với cả hai biến. Một là kết quả của một dị thường dữ liệu, khác do một lỗi có thể trong mã (có thể lỗi vì chúng ta chỉ thấy kết quả tính toán - suy nghĩ của biến x là "Garbage In", và y là "Garbage Out". , câu hỏi đặt ra là liệu các vấn đề trong đầu ra gây ra bởi mã chỉ là những trường hợp dữ liệu bất thường, hoặc có một số trường hợp khác có dữ liệu tốt không? (Đây là lý do tại sao tôi hỏi một câu hỏi về stopping when a NaN, NA, or Inf is encountered.)

Điều đó cũng giải thích lý do vì sao ví dụ của tôi có xác suất thấp cho các giá trị TRUE, những giá trị này thực sự xảy ra ít hơn 0,1% thời gian. (I E. các vị trí của TRUE trong mỗi bộ) và đếm các giao điểm đã đặt. Tôi tránh thiết lập giao lộ vì tôi đã bị đốt cháy một thời gian trước bởi Matlab (vâng, đây là R, nhưng chịu với tôi), mà trước tiên sẽ sắp xếp các yếu tố của một bộ trước khi nó thực hiện giao lộ. (Tôi mơ hồ nhớ lại sự phức tạp thậm chí còn lúng túng hơn:. Như O(n^2) thay vì O(n log n))

Answer 1

Dưới đây là kết quả cho các phương pháp logic, table, và bigtabulate, cho 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 

Trong trường hợp này , table là một thảm họa.

Để so sánh, đây là 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 

Tại thời điểm này, có vẻ như viết chức năng logic của chính mình là tốt nhất, mặc dù rằng lạm dụng sum, và kiểm tra từng vector logic nhiều lần. Tôi chưa thử biên dịch các chức năng, nhưng điều đó sẽ mang lại kết quả tốt hơn.

Cập nhật 1 Nếu chúng tôi đưa ra bigtabulate giá trị mà đã là các số nguyên, ví dụ: nếu chúng ta làm loại chuyển đổi 1 * cbind(v1,v2) ngoài bigtabulate thì nhiều N = 3E6 là 1.80, thay vì 2.4. N = 3E8 nhiều liên quan đến phương pháp "hợp lý" chỉ là 1,21, thay vì 1,53.

---

Cập nhật 2

Như Joshua Ulrich đã chỉ ra, chuyển sang vector bit là một cải tiến đáng kể - chúng tôi đang phân bổ và di chuyển xung quanh dữ liệu ít hơn rất nhiều: vectơ logic R của tiêu thụ 4 byte cho mỗi mục nhập ("Tại sao?", bạn có thể yêu cầu ... Well, I don't know, but an answer may turn up here.), trong khi một bit bit tiêu thụ, tốt, một bit, mỗi mục - tức là 1/32 càng nhiều dữ liệu. Vì vậy, x tiêu thụ 1,2e9 byte, trong khi xb (phiên bản bit trong mã bên dưới) chỉ tiêu thụ 3,75e7 byte.

Tôi đã giảm table và các biến thể bigtabulate từ các điểm chuẩn được cập nhật (N = 3e8). Lưu ý rằng logicalB1 giả định rằng dữ liệu đã là một bit bit, trong khi logicalB2 là hoạt động tương tự với hình phạt cho chuyển đổi loại. Khi các vectơ logic của tôi là kết quả của các phép toán trên các dữ liệu khác, tôi không có lợi ích của việc bắt đầu với một vectơ bit. Tuy nhiên, hình phạt phải trả là tương đối nhỏ. [Chuỗi "logical3" chỉ thực hiện 3 phép toán logic, và sau đó thực hiện phép trừ. Vì nó là cross-lập bảng, chúng ta biết tổng số, như DWin đã nhận xét.]

 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 

Bây giờ chúng ta đã tăng tốc này lên đến chụp chỉ 1,8-2,8 giây, thậm chí với nhiều sự thiếu hiệu quả tổng. Có không nghi ngờ gì là có thể thực hiện điều này trong vòng chưa đến 1 giây, với các thay đổi bao gồm một hoặc nhiều: mã C, biên dịch và xử lý đa lõi. Sau khi tất cả 3 (hoặc 4) hoạt động hợp lý khác nhau có thể được thực hiện một cách độc lập, mặc dù đó vẫn là một sự lãng phí của chu kỳ tính toán.

Tương tự nhất của những người thách thức tốt nhất, logical3B2, nhanh hơn khoảng 80X so với table. Đó là khoảng 10 X nhanh hơn so với hoạt động hợp lý ngây thơ. Và nó vẫn còn rất nhiều chỗ để cải thiện.

---

Đây là mã để sản xuất ở trên. CHÚ Ý Tôi khuyên bạn nên bình luận một số hoạt động hoặc vectơ, trừ khi bạn có nhiều RAM - việc tạo ra x, x1 và xb, cùng với các đối tượng y tương ứng sẽ chiếm một chút dung lượng bộ nhớ.

Ngoài ra, lưu ý: Tôi phải sử dụng 1L làm hệ số nguyên cho bigtabulate, thay vì chỉ 1. Tại một thời điểm nào đó, tôi sẽ chạy lại với thay đổi này và sẽ đề xuất thay đổi đó cho bất kỳ ai sử dụng phương thức 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

Nếu bạn đang thực hiện nhiều thao tác trên vectơ lôgic rất lớn, hãy xem gói bit. Nó tiết kiệm được một tấn bộ nhớ bằng cách lưu trữ các boolean như là các boolean 1 bit thực sự.

Điều này không giúp ích với table; nó thực sự làm cho nó tồi tệ hơn vì có nhiều giá trị duy nhất trong vector bit do cách nó được xây dựng. Nhưng nó thực sự sẽ giúp so sánh logic.

# 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

Một chiến thuật khác nhau là để xem xét chỉ cần đặt nút giao thông, sử dụng các chỉ số của TRUE giá trị, lợi dụng mà các mẫu được rất thiên vị (ví dụ: chủ yếu FALSE).

Để kết thúc, tôi giới thiệu func_find01 và bản dịch sử dụng gói bit (func_find01B); tất cả các mã không xuất hiện trong câu trả lời ở trên được dán bên dưới.

Tôi đã chạy lại toàn bộ đánh giá N = 3e8, ngoại trừ quên sử dụng func_find01B; Tôi đọc lại các phương pháp nhanh hơn chống lại nó, trong một lần thứ hai.

  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 

Chỉ cần "nhanh" phương pháp:

 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 

Vì vậy, find01B là nhanh nhất trong số các phương pháp không sử dụng vectơ chút tiền được chuyển đổi, bằng lãi mỏng (2,099 giây so với 2,327 giây). Trường hợp đã làm find02 đến từ đâu? Sau đó tôi đã viết một phiên bản sử dụng các bit bit được tính toán trước. Đây là nhanh nhất.

Nói chung, thời gian chạy của phương pháp "phương pháp chỉ số" có thể bị ảnh hưởng bởi xác suất chung & cận biên. Tôi nghi ngờ rằng nó sẽ đặc biệt cạnh tranh khi xác suất thậm chí còn thấp hơn, nhưng người ta phải biết rằng một ưu tiên, hoặc thông qua một mẫu phụ.

---

Update 1. Tôi cũng đã đề nghị đúng lúc Josh O'Brien của, sử dụng tabulate() thay vì table(). Kết quả, sau 12 giây trôi qua, khoảng 2X find01 và khoảng một nửa số bigtabulate2. Bây giờ mà các phương pháp tốt nhất đang tiếp cận 1 giây, điều này cũng là tương đối chậm:

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 4

Dưới đây là một câu trả lời với Rcpp, lập bảng chỉ những mục mà không phải là cả 0 . Tôi nghi ngờ phải có một số cách để cải thiện điều này, vì điều này là bất thường chậm; nó cũng là nỗ lực đầu tiên của tôi với Rcpp, vì vậy có thể có một số sự thiếu hiệu quả rõ ràng liên quan đến việc di chuyển dữ liệu xung quanh. Tôi đã viết một ví dụ đó là vanilla có mục đích rõ ràng, mà nên cho người khác chứng minh làm thế nào điều này có thể được cải thiện.

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

kết quả Timing cho N = 3E8:

user system elapsed 
10.930 1.620 12.586 

này có hơn 6X càng lâu càng func_find01B trong câu trả lời thứ 2 của tôi.

Answer 5

Đây là câu trả lời sử dụng đường 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) 
}) 

mà kết quả trong:

> 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 

Vì vậy, chúng ta có thể có được một chút tốc độ lên trên gói chút, mặc dù tôi ngạc nhiên khi thấy mức độ cạnh tranh trong thời gian được tính.

Cập nhật: Trong danh dự của Iterator, đây là một giải pháp 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