Improving stLMM runtime on multicore systems

Initial practical guidance for OpenMP, BLAS/LAPACK, and CHOLMOD

Andrew O. Finley

1 Purpose

This vignette gives practical starting points for improving stLMM runtime on multicore systems. The recommendations are intentionally cautious. Runtime depends on the model, data structure, sparse-matrix fill-in, hardware, memory bandwidth, BLAS/LAPACK implementation, and the Matrix/CHOLMOD build used by R.

More testing is needed across machines and model classes. Treat these notes as initial guidance rather than fixed rules. Users are encouraged to run small timing experiments on their own analyses, and feedback about configurations that help or hurt is very welcome.

2 Runtime Layers

There are three computational layers to keep separate.

First, stLMM has OpenMP-parallel sections in selected compiled kernels. These are controlled by the n_omp_threads argument. This is not a general “make everything parallel” switch; it only applies to code regions that were written with OpenMP.

Second, dense linear algebra is handled by BLAS and LAPACK (Blackford et al. 2002; Anderson et al. 1999). R supplies these libraries to packages. If R is linked against an optimized implementation such as OpenBLAS (OpenBLAS Project 2025), MKL, or BLIS, dense linear algebra can be faster. If that implementation is threaded, BLAS/LAPACK can also create its own worker threads.

Third, sparse Cholesky factorization is handled by CHOLMOD (Chen et al. 2008) through the R Matrix package (Bates et al. 2025). This layer is important for NNGP, CAR, CAR-time, and autoregressive terms because these terms produce sparse precision matrices. Sparse Cholesky runtime depends strongly on ordering, fill-in, supernode sizes, symbolic structure, and memory traffic.

These layers have different controls. n_omp_threads does not set OpenBLAS, MKL, or BLIS threads. BLAS thread controls do not necessarily control CHOLMOD’s own OpenMP threading. Increasing several thread counts at once can oversubscribe the machine and make a run slower.

3 Practical Observations

During development benchmarks, the most reliable improvement came from stLMM’s own OpenMP regions: runtime improved as n_omp_threads increased up to a point, then flattened or worsened. The useful number of threads depended on the machine and workload.

Threaded BLAS did not improve the benchmark workloads we tested. In those runs, single-threaded BLAS combined with stLMM OpenMP was fastest. When a verified threaded BLAS was used, CHOLMOD-related time became slower, and prediction-heavy timings did not improve. This is an observation from specific development benchmarks, not a guarantee for every model or machine.

CHOLMOD threading needs special care. Many R Matrix installations use CHOLMOD without CHOLMOD OpenMP support. In that case, CHOLMOD’s own OpenMP controls cannot help. Supernodal CHOLMOD can still call BLAS for dense frontal-matrix operations, so a threaded BLAS may be used inside CHOLMOD, but whether this helps depends on the sparse structure, fill-in, and size of dense updates. For some workloads, threaded BLAS overhead can outweigh any benefit.

The practical starting point is therefore:

• start with BLAS threads set to one;
• sweep n_omp_threads over a small grid appropriate for the machine;
• only then test BLAS threads if dense linear algebra is suspected to matter;
• change one layer at a time and keep the model, data, seed, and sample count fixed.

4 Checking the Current Setup

Before timing, record the R and numerical-library setup.

Code

sessionInfo()
extSoftVersion()["BLAS"]
packageVersion("Matrix")
packageVersion("stLMM")

Do not assume that an OpenBLAS library is threaded just because its path contains openblas. Some systems provide both single-threaded and threaded OpenBLAS builds. On Linux, ldd can help show which libraries R and packages load.

ldd "$(R RHOME)/lib/libR.so"
ldd "$(R RHOME)/lib/libRblas.so"

Rscript -e 'so <- list.files(system.file("libs", package = "Matrix"), full.names = TRUE); system2("ldd", so)'
Rscript -e 'so <- list.files(system.file("libs", package = "stLMM"), full.names = TRUE); system2("ldd", so)'

If you are testing threaded BLAS, verify actual thread use rather than relying only on environment variables. On Linux, one simple check during a running benchmark is:

ps -o pid,nlwp,pcpu,comm,args -p <R-process-id>

NLWP is the number of lightweight process threads. It should increase when a threaded BLAS or OpenMP region is active.

5 Thread Controls

For interactive work, RhpcBLASctl is often convenient if the linked BLAS supports runtime control.

Code

library(RhpcBLASctl)

RhpcBLASctl::blas_set_num_threads(1)

For reproducible benchmarks, prefer setting thread counts before R starts. For example, with OpenBLAS:

OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1 Rscript timing-script.R
OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=4 Rscript timing-script.R
OPENBLAS_NUM_THREADS=4 OMP_NUM_THREADS=1 Rscript timing-script.R

For stLMM’s own OpenMP regions, use the explicit model argument.

Code

fit <- stLMM(
  y ~ x + nngp(lon, lat, m = 15),
  data = dat,
  n_samples = 1000,
  n_omp_threads = 4,
  priors = list(
    resid = list(tau_sq = ig(2, 1)),
    nngp_1 = list(sigma_sq = ig(2, 1), phi = uniform(0.1, 10))
  )
)

Choose a thread grid that fits the machine. A laptop, a shared server, and a large workstation should not use the same grid. A reasonable first sweep is a small set of powers of two, capped at a sensible value for the machine.

Code

physical_cores <- parallel::detectCores(logical = FALSE)
if (is.na(physical_cores)) {
  physical_cores <- parallel::detectCores(logical = TRUE)
}

omp_grid <- unique(pmin(c(1, 2, 4, 8, 16, physical_cores), physical_cores))
omp_grid

On larger machines, it can also be useful to include the number of logical CPUs. For compute-heavy workloads, logical CPUs from simultaneous multithreading do not usually double throughput, but they are worth testing if the run is memory or latency limited.

6 Matrix, CHOLMOD, and Ordering

stLMM uses the CHOLMOD interface exposed by Matrix. Installing SuiteSparse or METIS system packages does not necessarily change an already installed Matrix package. METIS-based orderings are available only if the CHOLMOD library used by Matrix was built with that support.

The cholmod_control argument can be passed as a string, for example cholmod_control = "metis", or as a control list with an ordering entry. The supported ordering values are "auto", "best", "natural", "amd", "metis", "nesdis", and "colamd". The default "auto" uses CHOLMOD’s default strategy. The "best" option asks CHOLMOD to try its full built-in method set, which can make symbolic analysis slower.

For large sparse models, ordering is worth testing. The best ordering can differ by model structure, spatial/temporal support, neighbor count, and Matrix/CHOLMOD build.

Code

fit <- stLMM(
  y ~ x + nngp(lon, lat, m = 15),
  data = dat,
  n_samples = 1000,
  n_omp_threads = 4,
  cholmod_control = "amd",
  priors = list(
    resid = list(tau_sq = ig(2, 1)),
    nngp_1 = list(sigma_sq = ig(2, 1), phi = uniform(0.1, 10))
  )
)

7 Benchmarking Your Own Runs

The most useful benchmark is the target analysis, or a shortened version of it. Keep the model, data, seed, and number of MCMC samples fixed while changing one runtime setting at a time.

A first experiment is to hold BLAS threads at one and vary n_omp_threads.

run_fit <- function(omp_threads) {
  stLMM(
    y ~ x + nngp(lon, lat, m = 15),
    data = dat,
    n_samples = 200,
    n_omp_threads = omp_threads,
    priors = list(
      resid = list(tau_sq = ig(2, 1)),
      nngp_1 = list(sigma_sq = ig(2, 1), phi = uniform(0.1, 10))
    )
  )
}

system.time(run_fit(1))
system.time(run_fit(4))
system.time(run_fit(8))

For more reproducible comparisons, use separate R processes so thread environment variables are fixed before R starts.

OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1 Rscript timing-script.R
OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=4 Rscript timing-script.R
OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=8 Rscript timing-script.R

Then, if you want to test threaded BLAS, hold n_omp_threads fixed and vary the BLAS layer. Monitor actual thread use and elapsed time.

OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=4 Rscript timing-script.R
OPENBLAS_NUM_THREADS=2 OMP_NUM_THREADS=4 Rscript timing-script.R
OPENBLAS_NUM_THREADS=4 OMP_NUM_THREADS=4 Rscript timing-script.R

Before sharing benchmark results, record:

• CPU model, number of physical cores, number of logical CPUs, and memory;
• sessionInfo() and extSoftVersion()["BLAS"];
• Matrix and stLMM versions;
• loaded BLAS/LAPACK libraries, when available;
• thread settings such as OPENBLAS_NUM_THREADS, MKL_NUM_THREADS, BLIS_NUM_THREADS, and OMP_NUM_THREADS;
• model family, sparse term type, neighbor count, data size, and cholmod_control;
• elapsed time summaries across repeated runs when possible.

8 Practical Checklist

• Use an optimized BLAS/LAPACK implementation, but do not assume threaded BLAS will help.
• Start benchmarking with BLAS threads set to one.
• Sweep n_omp_threads over a machine-appropriate grid.
• Avoid setting BLAS and OpenMP thread counts high at the same time unless oversubscription is the experiment.
• For large sparse models, test CHOLMOD ordering choices.
• Verify actual thread use when testing threaded libraries.
• Repeat timing experiments for each major model/data class.
• Share feedback when alternative BLAS/LAPACK, Matrix/CHOLMOD builds, ordering choices, or thread settings materially improve performance.

Anderson, Edward, Zhaojun Bai, Christian Bischof, et al. 1999. LAPACK Users’ Guide. 3rd ed. Society for Industrial; Applied Mathematics.

Bates, Douglas, Martin Maechler, and Mikael Jagan. 2025. Matrix: Sparse and Dense Matrix Classes and Methods. https://doi.org/10.32614/CRAN.package.Matrix.

Blackford, L. Susan, James Demmel, Jack Dongarra, et al. 2002. “An Updated Set of Basic Linear Algebra Subprograms (BLAS).” ACM Transactions on Mathematical Software 28 (2): 135–51. https://doi.org/10.1145/567806.567807.

Chen, Yanqing, Timothy A. Davis, William W. Hager, and Sivasankaran Rajamanickam. 2008. “Algorithm 887: CHOLMOD, Supernodal Sparse Cholesky Factorization and Update/Downdate.” ACM Transactions on Mathematical Software 35 (3): 1–14. https://doi.org/10.1145/1391989.1391995.

OpenBLAS Project. 2025. OpenBLAS: An Optimized BLAS Library. https://www.openblas.net/.