Bixverse R <> Rust interface • bixverse

Why Rust?

Honest answer? The creation of this package was based on the desire of the author to play with Rust. To make it useful, the initial playground was stuff relevant to day-to-day work in bioinformatics and computational biology: correlations, PCA, other matrix factorisation methods, and so on. It quickly became apparent that Rust delivered incredible speed gains over R — sometimes 100x — with millions of hypergeometric tests completing in mere seconds. The results were striking. Memory safety and the fantastic work by the rextendr team made interfacing R and Rust straightforward and, frankly, quite fun.

Correlations and covariance

One of the most typical operations in computational biology is computing correlations and covariances; and these were the first functions implemented, largely to understand faer and the Rust/R interface. What becomes apparent very quickly is the following. For 1000 x 1000 matrices on an M1 Max MacBook Pro (using Apple’s Accelerate framework for R BLAS):

Covariance

Implementation	Median (ms, 20 runs)
base R	515.73
Rust	18.66

A 27x increase in speed.

Pearson correlation

Implementation	Median (ms, 20 runs)
base R	512.83
Rust	20.90

A 24x increase in speed.

Spearman correlation

Implementation	Median (ms, 20 runs)
base R	603.00
Rust	29.47

A 20x increase in speed.

The gains are even more pronounced on larger matrices. For a 1000 x 5000 matrix:

Pearson correlation (1000 x 5000)

Implementation	Median (ms, 20 runs)
base R	12994.37
Rust	419.88

A 30x increase in speed.

Given that matrix algebra underlies a large proportion of computational bioinformatics workflows, this makes a compelling case for Rust under the hood. The rest of this vignette shows how to use the exposed Rust functions directly, with runnable benchmarks you can verify on your own machine.

Using the Rust functions in more general R code

bixverse exposes a number of R/Rust functions that are useful beyond the higher-level package functionality. This vignette covers:

Correlations (via rs_cor()).
Covariance (via rs_covariance()).
Covariance to correlation (via rs_cov2cor()).
Distance calculations (via rs_dist()).
Correlations between two matrices (via rs_cor2()).
Mutual information between columns (via rs_mutual_information()).
Set similarities (via rs_set_similarity_list() or rs_set_similarity()).

These may be useful for your own package development.

Note

The vignette was built on a GitHub runner which is a 2-core machine with ~8 GB of memory. The performance differences tend to be (even) more pronounced on machines with more cores and computational power.

Set-up

library(bixverse)

Any function prefixed with rs_ calls directly into Rust and can panic if misused — check inputs before passing them through. There is a LOT of these functions in the package, and they are all exported. Why? I want other people to be able to integrate the code into their own R packages, or if they feel brave, they can just use the Rust crate directly, see here.

Correlations, covariance, and cosine distance

no_rows <- 500L
no_cols <- 500L

set.seed(10101L)

random_data <- matrix(
  data = rnorm(no_rows * no_cols),
  ncol = no_cols,
  nrow = no_rows
)

Pearson’s correlation

r_pearsons_res <- cor(random_data)

rust_pearsons_res <- rs_cor(random_data, spearman = FALSE)

# tiny numerical differences do exist
all.equal(r_pearsons_res, rust_pearsons_res, tolerance = 1e-15)
#> [1] TRUE

Let’s check the speed.

microbenchmark::microbenchmark(
  r = cor(random_data),
  rust = rs_cor(random_data, spearman = FALSE),
  times = 10L
)
#> Unit: milliseconds
#>  expr      min        lq      mean   median        uq       max neval
#>     r 79.88985 79.894583 80.025159 79.93054 79.951688 80.790340    10
#>  rust  2.59969  2.644867  3.223947  2.76099  3.768263  4.531062    10

Spearman’s correlation

What about ranked correlations?

r_spearman_res <- cor(random_data, method = "spearman")

rust_spearman_res <- rs_cor(random_data, spearman = TRUE)

all.equal(r_spearman_res, rust_spearman_res, tolerance = 1e-15)
#> [1] TRUE

And speed

microbenchmark::microbenchmark(
  r = cor(random_data, method = "spearman"),
  rust = rs_cor(random_data, spearman = TRUE),
  times = 10L
)
#> Unit: milliseconds
#>  expr       min         lq       mean     median         uq        max neval
#>     r 99.618797 100.198322 101.332465 100.360694 102.682619 104.420984    10
#>  rust  3.964636   4.230623   5.448958   4.635082   6.336115   9.317088    10

Covariance

Covariance matrices are also very common…

r_covar_res <- cov(random_data)

rust_covar_res <- rs_covariance(random_data)

all.equal(r_covar_res, rust_covar_res, tolerance = 1e-15)
#> [1] TRUE

And speed comparisons

microbenchmark::microbenchmark(
  r = cov(random_data),
  rust = rs_covariance(random_data),
  times = 10L
)
#> Unit: milliseconds
#>  expr       min        lq      mean    median        uq       max neval
#>     r 79.696258 79.718150 79.823597 79.752462 79.805791 80.462631    10
#>  rust  1.783492  1.808829  1.892554  1.862494  1.956059  2.120354    10

Covariance to correlation

r_cor_from_covar <- cov2cor(r_covar_res)

rust_cor_from_covar <- rs_cov2cor(rust_covar_res)

all.equal(r_cor_from_covar, rust_cor_from_covar, tolerance = 1e-15)
#> [1] TRUE

This is where the Rust implementation shows the least improvement — the operation is already cheap in base R — but the interface remains consistent with the rest of the rs_ family.

microbenchmark::microbenchmark(
  r = cov2cor(r_covar_res),
  rust = rs_cov2cor(rust_covar_res),
  times = 10L
)
#> Unit: microseconds
#>  expr      min       lq     mean    median       uq      max neval
#>     r 1483.614 1500.068 1767.007 1506.3335 1543.784 4002.222    10
#>  rust  788.757  797.950 1150.785  868.2005  920.394 3440.033    10

Correlations between two matrices

So far, all of the functions did pairwise column-based correlations. What about two correlations against each other?

no_rows_2 <- 500L
no_cols_2 <- 400L

set.seed(23L)

random_data_2 <- matrix(
  data = rnorm(no_rows_2 * no_cols_2, mean = 2, sd = 2),
  ncol = no_cols_2,
  nrow = no_rows_2
)

Let’s compare against R:

r_cor_2_matrices <- cor(random_data, random_data_2)

rust_cor_2_matrices <- rs_cor2(random_data, random_data_2, spearman = FALSE)

# small precision differences are further propagated here
all.equal(r_cor_2_matrices, rust_cor_2_matrices, tolerance = 1e-14)
#> [1] TRUE

And speed:

microbenchmark::microbenchmark(
  r = cor(random_data, random_data_2),
  rust = rs_cor2(random_data, random_data_2, spearman = FALSE),
  times = 10L
)
#> Unit: milliseconds
#>  expr       min         lq       mean    median         uq        max neval
#>     r 128.49934 128.509880 128.729562 128.58997 128.597501 129.937791    10
#>  rust   2.07034   2.133304   2.304129   2.15804   2.242196   3.363689    10

Distance metrics

Euclidean distance

Note that dist() computes row-wise distances, while rs_dist() expects columns to represent observations — hence the transpose.

r_euclidean_distance <- as.matrix(
  dist(random_data, method = "euclidean")
)

rs_euclidean_distance <- rs_dist(
  t(random_data),
  distance_type = "euclidean"
)

all.equal(
  r_euclidean_distance,
  rs_euclidean_distance,
  tolerance = 1e-14,
  check.attributes = FALSE
)
#> [1] TRUE

And speed…

microbenchmark::microbenchmark(
  r = dist(random_data, method = "euclidean"),
  rust = rs_dist(random_data, distance_type = "euclidean"),
  times = 10L
)
#> Unit: milliseconds
#>  expr       min        lq      mean    median        uq       max neval
#>     r 55.562845 56.182479 56.963754 56.424785 56.757617 62.293308    10
#>  rust  1.810081  1.822069  1.882675  1.873857  1.905884  2.082638    10

Manhattan distance

r_manhattan_distance <- as.matrix(
  dist(random_data, method = "manhattan")
)

rs_manhattan_distance <- rs_dist(
  t(random_data),
  distance_type = "manhattan"
)

all.equal(
  r_manhattan_distance,
  rs_manhattan_distance,
  tolerance = 1e-14,
  check.attributes = FALSE
)
#> [1] TRUE

And speed… ?

microbenchmark::microbenchmark(
  r = dist(random_data, method = "manhattan"),
  rust = rs_dist(random_data, distance_type = "manhattan"),
  times = 10L
)
#> Unit: milliseconds
#>  expr       min        lq      mean    median        uq      max neval
#>     r 55.373682 55.780619 55.982446 56.078118 56.253184 56.40636    10
#>  rust  9.574217  9.596755  9.739862  9.725167  9.788497 10.18708    10

Canberra distance

r_canberra_distance <- as.matrix(
  dist(random_data, method = "canberra")
)

rs_canberra_distance <- rs_dist(
  t(random_data),
  distance_type = "canberra"
)

all.equal(
  r_canberra_distance,
  rs_canberra_distance,
  tolerance = 1e-14,
  check.attributes = FALSE
)
#> [1] TRUE

Same benchmarks as before…

microbenchmark::microbenchmark(
  r = dist(random_data, method = "canberra"),
  rust = rs_dist(random_data, distance_type = "canberra"),
  times = 10L
)
#> Unit: milliseconds
#>  expr      min       lq     mean   median       uq      max neval
#>     r 60.09060 60.27037 60.99067 60.64449 60.93880 64.54331    10
#>  rust 11.72236 11.76186 11.87293 11.80268 11.83009 12.42257    10

Mutual information

set.seed(246L)

nrows <- 100
ncols <- 100

mat <- matrix(data = rnorm(nrows * ncols), nrow = nrows, ncol = ncols)
rownames(mat) <- sprintf("sample_%i", 1:nrows)
colnames(mat) <- sprintf("feature_%i", 1:ncols)

The Rust version defaults to sqrt(nrow()) bins when n_bins = NULL. Two discretisation strategies are available.

Equal width:

rust_res_mi <- rs_mutual_info(
  mat,
  n_bins = NULL,
  normalise = FALSE,
  strategy = "equal_width"
)
rownames(rust_res_mi) <- colnames(rust_res_mi) <- colnames(mat)

infotheo_res_mi <- infotheo::mutinformation(infotheo::discretize(
  mat,
  disc = "equalwidth",
  nbins = sqrt(nrow(mat))
))

all.equal(rust_res_mi, infotheo_res_mi)
#> [1] TRUE

microbenchmark::microbenchmark(
  infotheo = infotheo::mutinformation(infotheo::discretize(
    mat,
    disc = "equalwidth",
    nbins = sqrt(nrow(mat))
  )),
  rust = rs_mutual_info(
    mat,
    n_bins = NULL,
    normalise = FALSE,
    strategy = "equal_width"
  ),
  times = 10L
)
#> Unit: milliseconds
#>      expr      min       lq      mean    median        uq       max neval
#>  infotheo 57.63905 58.34736 58.525190 58.683416 58.888217 59.083498    10
#>      rust  1.90932  1.91704  1.983103  1.977391  2.001887  2.209116    10

Equal frequency:

rust_res_mi <- rs_mutual_info(
  mat,
  n_bins = NULL,
  normalise = FALSE,
  strategy = "equal_freq"
)
rownames(rust_res_mi) <- colnames(rust_res_mi) <- colnames(mat)

infotheo_res_mi <- infotheo::mutinformation(infotheo::discretize(
  mat,
  disc = "equalfreq",
  nbins = sqrt(nrow(mat))
))

all.equal(rust_res_mi, infotheo_res_mi)
#> [1] TRUE

microbenchmark::microbenchmark(
  infotheo = infotheo::mutinformation(infotheo::discretize(
    mat,
    disc = "equalfreq",
    nbins = sqrt(nrow(mat))
  )),
  rust = rs_mutual_info(
    mat,
    n_bins = NULL,
    normalise = FALSE,
    strategy = "equal_freq"
  ),
  times = 10L
)
#> Unit: milliseconds
#>      expr       min        lq      mean    median        uq       max neval
#>  infotheo 70.271203 70.436259 71.103496 70.612161 70.859521 75.733284    10
#>      rust  2.623295  2.645007  2.709139  2.705528  2.729373  2.883623    10

Set similarities

bixverse also exposes Rust-accelerated set similarity calculations. To illustrate the performance, below is a comparison against a naive R implementation of Jaccard similarity across two lists.

jaccard_sim <- function(x, y) {
  length(intersect(x, y)) / length(union(x, y))
}

set.seed(42L)

random_sets_1 <- purrr::map(
  1:1000L,
  ~ paste(sample(LETTERS, sample(1:20, 1)))
)

random_sets_2 <- purrr::map(
  1:500L,
  ~ paste(sample(LETTERS, sample(1:20, 1)))
)

Starting with a sequential purrr approach as the baseline:

tictoc::tic()

r_results <- purrr::map(
  random_sets_1,
  \(x) {
    purrr::map_dbl(random_sets_2, \(y) jaccard_sim(x, y))
  },
  .progress = TRUE
)
#>  ■■■■■■■■■                         27% |  ETA: 10s
#>  ■■■■■■■■■■■■■■■                   46% |  ETA:  8s
#>  ■■■■■■■■■■■■■■■■■■■■■             68% |  ETA:  5s
#>  ■■■■■■■■■■■■■■■■■■■■■■■■■■■■      89% |  ETA:  2s

similarity_matrix <- matrix(
  data = unlist(r_results),
  nrow = length(random_sets_1),
  ncol = length(random_sets_2),
  byrow = TRUE
)

tictoc::toc()
#> 14.186 sec elapsed

Parallelising via furrr gives a meaningful speedup:

future::plan(strategy = future::multisession(workers = parallel::detectCores()))

tictoc::tic()

r_results <- furrr::future_map(
  random_sets_1,
  \(x) {
    purrr::map_dbl(random_sets_2, \(y) jaccard_sim(x, y))
  },
  .progress = TRUE
)

similarity_matrix <- matrix(
  data = unlist(r_results),
  nrow = length(random_sets_1),
  ncol = length(random_sets_2),
  byrow = TRUE
)

tictoc::toc()
#> 6.608 sec elapsed

future::plan(strategy = future::sequential())

mirai reduces process scheduling overhead (somtimes) compared to future and does a bit better:

mirai::daemons(parallel::detectCores())

tictoc::tic()

r_results_mirai <- mirai::mirai_map(
  random_sets_1,
  \(x, sets_2) {
    purrr::map_dbl(sets_2, \(y) {
      length(intersect(x, y)) / length(union(x, y))
    })
  },
  .args = list(sets_2 = random_sets_2)
)[]

similarity_matrix <- matrix(
  data = unlist(r_results_mirai),
  nrow = length(random_sets_1),
  ncol = length(random_sets_2),
  byrow = TRUE
)

tictoc::toc()
#> 6.598 sec elapsed

mirai::daemons(0)

Depending on your system this optimised parallel R version may be around 4x faster than the naive sequential baseline. Now compare against the Rust implementation:

tictoc::tic()

rust_res <- rs_set_similarity_list2(
  s_1_list = random_sets_1,
  s_2_list = random_sets_2,
  overlap_coefficient = FALSE
)

tictoc::toc()
#> 0.031 sec elapsed

all.equal(similarity_matrix, rust_res, tolerance = 1e-15)
#> [1] TRUE

The margin here is not subtle. The functions shown in this vignette are a representative sample of what is exposed in bixverse: there are many more specialised functions throughout the package worth exploring (if you are brave and can deal with a panic here and there).