High Performance Benchmarks

Joseph Wood

03/29/2022

This document serves as an overview for measuring the performance of RcppAlgos against other tools for generating combinations, permutations, and partitions. This stackoverflow post: Permutations and combinations with/without replacement and for distinct/non-distinct items/multiset has some benchmarks. You will note that the examples in that post were relatively small. The benchmarks below will focus on larger examples where performance really matters and for this reason we only consider the packages arrangements, partitions, and RcppAlgos.

Setup Information

For the benchmarks below, we used a Macbook Pro i7 16Gb machine. We also tested on a Windows and Linux machine with similar specs and obtained similar results.

library(RcppAlgos)
library(partitions)
library(arrangements)
library(microbenchmark)
pertinent_output <- capture.output(sessionInfo())

options(digits = 4)
cat(paste(pertinent_output[1:3], collapse = "\n"))
#> R version 4.1.3 (2022-03-10)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS Monterey 12.3

cat(pertinent_output[which(pertinent_output == "other attached packages:") + 1L])
#> [1] microbenchmark_1.4-7 arrangements_1.1.9   partitions_1.10-2    RcppAlgos_2.5.1

numThreads <- as.integer(RcppAlgos::stdThreadMax() / 2)
print(numThreads)
#> [1] 4

Combinations

Combinations - Distinct

set.seed(13)
v1 <- sort(sample(100, 30))
m <- 21
t1 <- comboGeneral(v1, m, Parallel = T)
t2 <- combinations(v1, m)
identical(t1, t2)
#> [1] TRUE
dim(t1)
#> [1] 14307150       21
rm(t1, t2)
invisible(gc())
microbenchmark(cbRcppAlgosPar = comboGeneral(v1, m, nThreads = numThreads),
               cbRcppAlgosSer = comboGeneral(v1, m),
               cbArrangements = combinations(v1, m),
               times = 15, unit = "relative")
#> Unit: relative
#>            expr   min    lq  mean median    uq   max neval cld
#>  cbRcppAlgosPar 1.000 1.000 1.000  1.000 1.000 1.000    15  a
#>  cbRcppAlgosSer 3.581 2.518 2.516  2.536 2.387 2.336    15   b
#>  cbArrangements 3.564 2.502 2.510  2.510 2.398 2.372    15   b

Combinations - Repetition

v2 <- v1[1:10]
m <- 20
t1 <- comboGeneral(v2, m, repetition = TRUE, nThreads = numThreads)
t2 <- combinations(v2, m, replace = TRUE)
identical(t1, t2)
#> [1] TRUE
dim(t1)
#> [1] 10015005       20
rm(t1, t2)
invisible(gc())
microbenchmark(cbRcppAlgosPar = comboGeneral(v2, m, TRUE, nThreads = numThreads),
               cbRcppAlgosSer = comboGeneral(v2, m, TRUE),
               cbArrangements = combinations(v2, m, replace = TRUE),
               times = 15, unit = "relative")
#> Unit: relative
#>            expr   min    lq  mean median    uq   max neval cld
#>  cbRcppAlgosPar 1.000 1.000 1.000  1.000 1.000 1.000    15 a
#>  cbRcppAlgosSer 2.350 2.554 2.459  2.538 2.462 2.186    15   c
#>  cbArrangements 2.421 2.420 2.369  2.420 2.344 2.341    15  b

Combinations - Multisets

myFreqs <- c(2, 4, 4, 5, 3, 2, 2, 2, 3, 4, 1, 4, 2, 5)
v3 <- as.integer(c(1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610))
t1 <- comboGeneral(v3, 20, freqs = myFreqs, nThreads = numThreads)
t2 <- combinations(freq = myFreqs, k = 20, x = v3)
identical(t1, t2)
#> [1] TRUE
dim(t1)
#> [1] 14594082       20
rm(t1, t2)
invisible(gc())
microbenchmark(cbRcppAlgosPar = comboGeneral(v3, 20, freqs = myFreqs, nThreads = numThreads),
               cbRcppAlgosSer = comboGeneral(v3, 20, freqs = myFreqs),
               cbArrangements = combinations(freq = myFreqs, k = 20, x = v3),
               times = 10, unit = "relative")
#> Unit: relative
#>            expr   min    lq  mean median    uq   max neval cld
#>  cbRcppAlgosPar 1.000 1.000 1.000  1.000 1.000 1.000    10 a
#>  cbRcppAlgosSer 2.608 2.589 2.378  2.505 2.107 2.160    10  b
#>  cbArrangements 3.776 3.931 3.634  3.829 3.220 3.397    10   c

Permutations

Permutations - Distinct

v4 <- as.integer(c(2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59))
t1 <- permuteGeneral(v4, 6, nThreads = numThreads)
t2 <- permutations(v4, 6)
identical(t1, t2)
#> [1] TRUE
dim(t1)
#> [1] 8910720       6
rm(t1, t2)
invisible(gc())
microbenchmark(cbRcppAlgosPar = permuteGeneral(v4, 6, nThreads = numThreads),
               cbRcppAlgosSer = permuteGeneral(v4, 6),
               cbArrangements = permutations(v4, 6),
               times = 15, unit = "relative")
#> Unit: relative
#>            expr   min    lq  mean median    uq   max neval cld
#>  cbRcppAlgosPar 1.000 1.000 1.000  1.000 1.000 1.000    15 a
#>  cbRcppAlgosSer 2.324 2.340 2.178  2.380 2.436 1.527    15  b
#>  cbArrangements 3.050 3.035 2.893  3.033 3.524 1.950    15   c


## Indexing permutation example with the partitions package
t1 <- permuteGeneral(11, nThreads = 4)
t2 <- permutations(11)
t3 <- perms(11)

dim(t1)
#> [1] 39916800       11

identical(t1, t2)
#> [1] TRUE
identical(t1, t(as.matrix(t3)))
#> [1] TRUE
rm(t1, t2, t3)
invisible(gc())

microbenchmark(cbRcppAlgosPar = permuteGeneral(11, nThreads = 4),
               cbRcppAlgosSer = permuteGeneral(11),
               cbArrangements = permutations(11),
               cbPartitions   = perms(11),
               times = 5, unit = "relative")
#> Unit: relative
#>            expr    min    lq  mean median    uq   max neval cld
#>  cbRcppAlgosPar  1.000 1.000 1.000  1.000 1.000 1.000     5 a
#>  cbRcppAlgosSer  3.129 3.299 2.621  3.284 1.947 2.289     5  b
#>  cbArrangements  3.664 3.630 2.791  3.639 2.201 2.032     5  b
#>    cbPartitions 10.040 9.932 7.974 10.901 6.497 5.797     5   c

Permutations - Repetition

v5 <- v3[1:5]
t1 <- permuteGeneral(v5, 10, repetition = TRUE, nThreads = numThreads)
t2 <- permutations(v5, 10, replace = TRUE)
identical(t1, t2)
#> [1] TRUE
dim(t1)
#> [1] 9765625      10
rm(t1, t2)
invisible(gc())
microbenchmark(cbRcppAlgosPar = permuteGeneral(v5, 10, TRUE, nThreads = numThreads),
               cbRcppAlgosSer = permuteGeneral(v5, 10, TRUE),
               cbArrangements = permutations(x = v5, k = 10, replace = TRUE),
               times = 10, unit = "relative")
#> Unit: relative
#>            expr   min    lq  mean median    uq   max neval cld
#>  cbRcppAlgosPar 1.000 1.000 1.000  1.000 1.000 1.000    10  a
#>  cbRcppAlgosSer 3.261 3.300 4.096  3.296 3.412 6.233    10   b
#>  cbArrangements 3.185 3.155 3.067  3.121 3.109 2.655    10   b

Permutations - Multisets

v6 <- sort(runif(12))
t1 <- permuteGeneral(v6, 7, freqs = rep(1:3, 4), nThreads = numThreads)
t2 <- permutations(freq = rep(1:3, 4), k = 7, x = v6)
identical(t1, t2)
#> [1] TRUE
dim(t1)
#> [1] 19520760        7
rm(t1, t2)
invisible(gc())
microbenchmark(cbRcppAlgosPar = permuteGeneral(v6, 7, freqs = rep(1:3, 4), nThreads = numThreads),
               cbRcppAlgosSer = permuteGeneral(v6, 7, freqs = rep(1:3, 4)),
               cbArrangements = permutations(freq = rep(1:3, 4), k = 7, x = v6),
               times = 10, unit = "relative")
#> Unit: relative
#>            expr   min    lq  mean median    uq   max neval cld
#>  cbRcppAlgosPar 1.000 1.000 1.000  1.000 1.000 1.000    10  a
#>  cbRcppAlgosSer 3.396 3.412 2.370  2.930 1.600 1.971    10   b
#>  cbArrangements 3.539 3.552 2.521  3.021 2.074 1.990    10   b

Partitions

Partitions - Distinct

All Distinct Partitions

t1 <- comboGeneral(0:140, freqs=c(140, rep(1, 140)),
                   constraintFun = "sum", comparisonFun = "==",
                   limitConstraints = 140)
t2 <- partitions(140, distinct = TRUE)
t3 <- diffparts(140)

# Each package has different output formats... we only examine dimensions
# and that each result is a partition of 140
identical(dim(t1), dim(t2))
#> [1] TRUE
identical(dim(t1), dim(t(t3)))
#> [1] TRUE
all(rowSums(t1) == 140)
#> [1] TRUE
all(rowSums(t2) == 140)
#> [1] TRUE
all(colSums(t3) == 140)
#> [1] TRUE

dim(t1)
#> [1] 9617150      16
rm(t1, t2, t3)
invisible(gc())
microbenchmark(cbRcppAlgosPar = partitionsGeneral(0:140, freqs=c(140, rep(1, 140)), nThreads = numThreads),
               cbRcppAlgosSer = partitionsGeneral(0:140, freqs=c(140, rep(1, 140))),
               cbArrangements = partitions(140, distinct = TRUE),
               cbPartitions   = diffparts(140),
               times = 10, unit = "relative")
#> Unit: relative
#>            expr    min     lq   mean median     uq    max neval cld
#>  cbRcppAlgosPar  1.000  1.000  1.000  1.000  1.000  1.000    10 a
#>  cbRcppAlgosSer  3.033  3.055  3.030  3.105  3.613  2.345    10  b
#>  cbArrangements  2.843  2.875  2.939  3.184  3.860  2.202    10  b
#>    cbPartitions 19.742 20.600 18.420 20.926 20.458 12.175    10   c

Restricted Distinct Partitions

t1 <- comboGeneral(160, 10,
                   constraintFun = "sum", comparisonFun = "==",
                   limitConstraints = 160)
t2 <- partitions(160, 10, distinct = TRUE)
identical(t1, t2)
#> [1] TRUE
dim(t1)
#> [1] 8942920      10
rm(t1, t2)
invisible(gc())
microbenchmark(cbRcppAlgosPar = partitionsGeneral(160, 10, nThreads = numThreads),
               cbRcppAlgosSer = partitionsGeneral(160, 10),
               cbArrangements = partitions(160, 10, distinct = TRUE),
               times = 10, unit = "relative")
#> Unit: relative
#>            expr   min    lq  mean median    uq   max neval cld
#>  cbRcppAlgosPar 1.000 1.000 1.000  1.000 1.000 1.000    10  a
#>  cbRcppAlgosSer 3.134 3.141 3.002  3.143 3.342 1.838    10   b
#>  cbArrangements 3.430 3.426 3.066  3.394 3.379 1.925    10   b

Partitions - Repetition

All Partitions

t1 <- comboGeneral(0:65, repetition = TRUE, constraintFun = "sum",
                   comparisonFun = "==", limitConstraints = 65)
t2 <- partitions(65)
t3 <- parts(65)

# Each package has different output formats... we only dimensions
# and that each result is a partition of 65
identical(dim(t1), dim(t2))
#> [1] TRUE
identical(dim(t1), dim(t(t3)))
#> [1] TRUE
all(rowSums(t1) == 65)
#> [1] TRUE
all(rowSums(t2) == 65)
#> [1] TRUE
all(colSums(t3) == 65)
#> [1] TRUE
dim(t1)
#> [1] 2012558      65
rm(t1, t2, t3)
invisible(gc())
microbenchmark(cbRcppAlgosPar = partitionsGeneral(0:65, repetition = TRUE,
                                                  nThreads = numThreads),
               cbRcppAlgosSer = partitionsGeneral(0:65, repetition = TRUE),
               cbArrangements = partitions(65),
               cbPartitions   = parts(65),
               times = 20, unit = "relative")
#> Unit: relative
#>            expr    min     lq  mean median    uq   max neval cld
#>  cbRcppAlgosPar  1.000  1.000 1.000  1.000 1.000 1.000    20 a
#>  cbRcppAlgosSer  2.996  2.985 2.362  2.691 2.102 1.688    20  b
#>  cbArrangements  2.901  2.844 2.297  2.608 2.042 1.831    20  b
#>    cbPartitions 10.852 10.695 7.695  9.581 6.311 4.514    20   c

Restricted Partitions

t1 <- comboGeneral(100, 15, TRUE, constraintFun = "sum",
                   comparisonFun = "==", limitConstraints = 100)
t2 <- partitions(100, 15)
identical(t1, t2)
#> [1] TRUE
dim(t1)
#> [1] 9921212      15
rm(t1, t2)

# This takes a really long time... not because of restrictedparts,
# but because apply is not that fast. This transformation is
# needed for proper comparisons. As a result, we will compare
# a smaller example
# t3 <- t(apply(as.matrix(restrictedparts(100, 15, include.zero = F)), 2, sort))
t3 <- t(apply(as.matrix(restrictedparts(50, 15, include.zero = F)), 2, sort))
identical(partitions(50, 15), t3)
#> [1] TRUE
rm(t3)
invisible(gc())
microbenchmark(cbRcppAlgosPar = partitionsGeneral(100, 15, TRUE,
                                                  nThreads = numThreads),
               cbRcppAlgosSer = partitionsGeneral(100, 15, TRUE),
               cbArrangements = partitions(100, 15),
               cbPartitions   = restrictedparts(100, 15,
                                                include.zero = FALSE),
               times = 10, unit = "relative")
#> Unit: relative
#>            expr    min     lq   mean median    uq   max neval cld
#>  cbRcppAlgosPar  1.000  1.000  1.000  1.000 1.000 1.000    10 a
#>  cbRcppAlgosSer  3.054  3.049  2.222  1.857 2.132 2.115    10  b
#>  cbArrangements  3.035  3.047  2.191  1.853 2.088 2.096    10  b
#>    cbPartitions 15.445 15.745 10.670  9.610 9.221 9.002    10   c

Partitions - Multisets

Currenlty, RcppAlgos is the only package capable of efficiently generating partitions of multisets. Therefore, we will only time RcppAlgos and use this as a reference for future improvements.

t1 <- comboGeneral(120, 10, freqs=rep(1:8, 15),
                   constraintFun = "sum", comparisonFun = "==",
                   limitConstraints = 120)
dim(t1)
#> [1] 7340225      10
all(rowSums(t1) == 120)
#> [1] TRUE
microbenchmark(cbRcppAlgos = partitionsGeneral(120, 10, freqs=rep(1:8, 15)),
               times = 10)
#> Unit: milliseconds
#>         expr   min      lq   mean median     uq    max neval
#>  cbRcppAlgos 715.2   724.8  750.5  735.3  791.5  804.5    10

### In RcppAlgos 2.3.6 - 2.4.3
#> Unit: seconds
#>         expr   min      lq   mean median     uq    max neval
#>  cbRcppAlgos 1172.9 1175.1 1223.5 1193.3 1200.5 1482.0    10    10

Iterators

We will show one example from each category to demonstrate the efficiency of the iterators in RcppAlgos. The results are similar for the rest of the cases not shown.

Combinations

pkg_arrangements <- function(n, total) {
    a <- icombinations(n, as.integer(n / 2))
    for (i in 1:total) a$getnext()
}

pkg_RcppAlgos <- function(n, total) {
    a <- comboIter(n, as.integer(n / 2))
    for (i in 1:total) a@nextIter()
}

total <- comboCount(18, 9)
total
#> [1] 48620

microbenchmark(cbRcppAlgos    = pkg_RcppAlgos(18, total),
               cbArrangements = pkg_arrangements(18, total),
               times = 15, unit = "relative")
#> Unit: relative
#>            expr   min    lq  mean median uq   max neval cld
#>     cbRcppAlgos  1.00  1.00  1.00   1.00  1  1.00    15  a 
#>  cbArrangements 26.56 25.94 24.94  25.39 25 19.61    15   b

Permutations

pkg_arrangements <- function(n, total) {
    a <- ipermutations(n)
    for (i in 1:total) a$getnext()
}

pkg_RcppAlgos <- function(n, total) {
    a <- permuteIter(n)
    for (i in 1:total) a@nextIter()
}

total <- permuteCount(8)
total
#> [1] 40320

microbenchmark(cbRcppAlgos    = pkg_RcppAlgos(8, total),
               cbArrangements = pkg_arrangements(8, total),
               times = 15, unit = "relative")
#> Unit: relative
#>            expr   min    lq  mean median    uq   max neval cld
#>     cbRcppAlgos  1.00  1.00  1.00   1.00  1.00  1.00    15  a 
#>  cbArrangements 26.12 25.83 24.53  24.95 24.03 22.15    15   b

Partitions

pkg_partitions <- function(n, total) {
    a <- firstpart(n)
    for (i in 1:(total - 1)) a <- nextpart(a)
}

pkg_arrangements <- function(n, total) {
    a <- ipartitions(n)
    for (i in 1:total) a$getnext()
}

pkg_RcppAlgos <- function(n, total) {
    a <- partitionsIter(0:n, repetition = TRUE)
    for (i in 1:total) a@nextIter()
}

total <- partitionsCount(0:40, repetition = TRUE)
total
#> [1] 37338

microbenchmark(cbRcppAlgos    = pkg_RcppAlgos(40, total),
               cbArrangements = pkg_arrangements(40, total),
               cbPartitions   = pkg_partitions(40, total),
               times = 15, unit = "relative")
#> Unit: relative
#>            expr   min    lq  mean median    uq   max neval cld
#>     cbRcppAlgos  1.00  1.00  1.00   1.00  1.00  1.00    15 a  
#>  cbArrangements 20.84 20.97 20.53  20.51 20.69 18.32    15  b 
#>    cbPartitions 38.10 38.29 37.01  37.66 37.01 30.80    15   c