Performance

The sequential Poisson method, and order sampling methods more generally, are simple and consequently not computationally expensive. This makes them suitable for a range of different applications, especially when drawing a sample from a large population. Despite this, there are two optimizations in this package to keep drawing samples fast.

Calculating inclusion probabilities

The first optimization concerns the inclusion probabilities when there are take-all units. As seen in vignette("take-all"), the algorithm finds take-all units one at a time without recomputing the inclusion probabilities many times. This is much faster than the naive approach when drawing a large sample, and is on par with the usual algorithm that finds take-all units in batches.

library(sps)

# Make a population with 200 take-all units.
x <- c(rep(1, 1e6 - 200), rep(1e6, 200))
n <- 1e3


# Naive implementation.
ip <- function(x, n, alpha = 0.001) {
  p <- \(x, n) x * (n / sum(x))
  ta_units <- integer(0)
  pi <- p(x, n)
  max_ts <- which.max(pi)
  while (pi[max_ts] > 1 - alpha) {
    ta_units <- c(ta_units, max_ts)
    pi <- p(replace(pi, max_ts, 0), n - length(ta_units))
    max_ts <- which.max(pi)
  }
  replace(pi, ta_units, 1)
}

cols <- c("expression", "median", "mem_alloc", "n_itr")
 
bench::mark(
  inclusion_prob(x, n),
  ip(x, n),
  sampling::inclusionprobabilities(x, n)
)[cols]

#> # A tibble: 3 × 4
#>   expression                               median mem_alloc n_itr
#>   <bch:expr>                             <bch:tm> <bch:byt> <int>
#> 1 inclusion_prob(x, n)                    24.18ms    61.1MB    21
#> 2 ip(x, n)                                  1.41s       3GB     1
#> 3 sampling::inclusionprobabilities(x, n)   70.1ms     134MB     8

Partial sorting

The second optimization is recognizing that both the computation of inclusion probabilities and the sequential Poisson method can benefit from partial sorting algorithms. In both cases, only the \(n\) largest/smallest elements of a vector are needed and, when installed, the topn() function from the kit package is used to avoid a complete sort. This is a drop-in replacement for order(x)[1:n] that can be faster when n is much smaller than the length of x, and can have a modest impact on performance when drawing a sample from a large population.

options(sps.usekit = TRUE)
bench::mark(sps(x, n))[cols]

#> # A tibble: 1 × 4
#>   expression   median mem_alloc n_itr
#>   <bch:expr> <bch:tm> <bch:byt> <int>
#> 1 sps(x, n)    75.5ms     111MB     6

options(sps.usekit = FALSE)
bench::mark(sps(x, n))[cols]

#> # A tibble: 1 × 4
#>   expression   median mem_alloc n_itr
#>   <bch:expr> <bch:tm> <bch:byt> <int>
#> 1 sps(x, n)     153ms     126MB     4

Partial sorting generally speeds up drawing smaller samples as well, but the effect is not as large because the sequential Poisson method is already quick.

Iteratively drawing a sample

Both optimizations above can be seen when iteratively drawing a sequential Poisson sample as creating an iterator involves finding the point at which each unit enters the take-all stratum. Although this is more expensive than drawing a single sample, it quickly pays off because of the small cost to iterating over units in the sample.

bench::mark(s <- sps_iterator(x, n), min_iterations = 10)[cols]

#> # A tibble: 1 × 4
#>   expression                median mem_alloc n_itr
#>   <bch:expr>              <bch:tm> <bch:byt> <int>
#> 1 s <- sps_iterator(x, n)    931ms     336MB    10

bench::mark(s())[cols]

#> # A tibble: 1 × 4
#>   expression   median mem_alloc n_itr
#>   <bch:expr> <bch:tm> <bch:byt> <int>
#> 1 s()          2.82ms    15.2MB    80

bench::mark(sps(x, n))

#> # A tibble: 1 × 13
#>   expression      min median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc total_time
#>   <bch:expr> <bch:tm> <bch:>     <dbl> <bch:byt>    <dbl> <int> <dbl>   <bch:tm>
#> 1 sps(x, n)     158ms  159ms      6.13     126MB     7.66     4     5      653ms
#> # ℹ 4 more variables: result <list>, memory <list>, time <list>, gc <list>