Parallelism: who owns the cores

Two layers of work

A T-Rex run has two layers that could use several cores. The outer layer is the \(K\) random experiments, which are independent of each other and share the design. The inner layer is the work inside one T-solver: the passes over the design that refresh the correlations, the search for the next knot of the path, the standardization. The rule of the package is that one layer owns the cores at a time. Either the experiments run in parallel and every solver runs serially, or the experiments run in sequence and a solver may split its passes. Nested teams of threads are never formed.

The default is serial. A selector runs its experiments one after the other, every solver runs on one thread, and a T-solver used on its own runs serially. Parallel work is opt-in through two options of trex_control(): max_outer_threads for the experiments and max_inner_threads for the solvers’ passes.

The default thread budget of the package

The package is single-threaded by default on every machine. Every parallel region of the backend runs on an explicit budget, and the methods without an option of their own (the scalers, the clustering, the ML regressions, trex_normalize_design(), the dummy and data generators) take the package’s default budget, which starts at 1. trex_get_num_threads() reads that default, trex_set_num_threads() sets it for the session and returns the previous value, and the n_threads argument of those methods overrides it per object. Inside a selector the same kernels take max_inner_threads. trex_get_max_threads() reports the threads OpenMP can run on this machine, an upper bound for a budget. In a build without OpenMP every setting is a no-op and the machine count is one.

trex_get_num_threads()
#> [1] 2
trex_get_max_threads()
#> [1] 12

The selectors and the solvers do not read the default budget. Their budget comes from the control.

The solvers’ passes

trex_control(max_inner_threads = m) gives every solver inside the random experiments a budget of m threads for its passes over the design (the correlation refreshes, the step-length search, the standardization). The default is 1. A pass is split only when it touches at least a million entries of the design, because the cost of forming a team of threads exceeds the gain below that, so small problems run serially whatever the budget says. The seeded generators (the dummies, the synthetic designs) count their work in the same unit with the measured cost of a random draw (about 80 entries) and of seeding an engine (up to 32,000), so they split from a few thousand entries. The budget is live when the experiments run in sequence. Inside parallel experiments the solvers stay serial. A T-solver used on its own takes the same budget through tsolver_options(n_threads = ), and $set_num_threads() changes it on a live solver.

Parallel random experiments

trex_control(parallel_rnd_experiments = TRUE, max_outer_threads = m) runs the random experiments on m threads (at most \(K\)). Inside the parallel experiments every solver runs serially, so the number of threads in use is max_outer_threads and nothing more. Both options are whole numbers of at least 1. The lean L-loop strategies ("SEEDED", "PERMUTATION", "PERMUTATION_SEEDED") run their experiments in sequence in every case and give the cores to the solvers.

set.seed(12)
n <- 100
p <- 150
beta <- c(rep(1, 5), rep(0, p - 5))
X <- matrix(rnorm(n * p), n, p)
y <- as.numeric(X %*% beta + rnorm(n))
serial <- TRexSelector$new(X, y, tFDR = 0.1, seed = 42, verbose = FALSE,
                           control = trex_control(K = 10))
serial$select()
parallel <- TRexSelector$new(X, y, tFDR = 0.1, seed = 42, verbose = FALSE,
                             control = trex_control(K = 10, parallel_rnd_experiments = TRUE,
                                                    max_outer_threads = 2))
parallel$select()
budget <- TRexSelector$new(X, y, tFDR = 0.1, seed = 42, verbose = FALSE,
                           control = trex_control(K = 10, max_inner_threads = 2))
budget$select()
identical(serial$selected_indices, parallel$selected_indices)
#> [1] TRUE
identical(serial$selected_indices, budget$selected_indices)
#> [1] TRUE

Determinism

A seeded selector derives one seed per random experiment from seed and the experiment index, so the dummies of experiment \(k\) do not depend on which thread draws them or on how many threads there are. The selection of a seeded run is therefore the same with experiments in sequence and in parallel, and the same for any thread count of the inner passes. Every reduction inside the solvers is summed in a fixed order for the same reason.

Combining with your own parallel loops

A Monte Carlo study usually parallelizes over trials. Then the trials own the cores, and the package should run serially inside every worker, which is the default. Two details matter.

trials <- parallel::mclapply(seq_len(40), function(mc) {
  trex_set_num_threads(1)
  set.seed(mc)
  X <- matrix(rnorm(n * p), n, p)
  y <- as.numeric(X %*% beta + rnorm(n))
  sel <- TRexSelector$new(X, y, tFDR = 0.1, verbose = FALSE)   # seed = -1 per trial
  sel$select()
  trex_fdp(sel$selected_indices, which(beta != 0))
}, mc.cores = 4)
mean(unlist(trials))

Builds without OpenMP

On Linux and Windows R supplies the OpenMP flags. CRAN’s R for macOS has none, and the package then links the OpenMP runtime R ships with the header of Homebrew’s libomp, or compiles single-threaded when no header is found. TREX_OPENMP=no at installation forces the single-threaded build. A single-threaded build selects the same variables, it only runs the experiments and the passes in sequence.

The CRAN rule

CRAN checks packages on at most two cores. The vignettes and the tests of this package follow the rule with the idiom of their setup chunk, which caps the count for the duration of the page and restores it at the end.

old_threads <- trex_set_num_threads(min(2L, trex_get_max_threads()))
# ... the work ...
trex_set_num_threads(old_threads)

References