Benchmarking matrix inversion in R

matrix_inversion.R

library(tidyverse)

inv <- function(m) {
  n <- nrow(m)
  stopifnot(n == ncol(m))
  
  cofactor <- function(m, i, j) (-1)^(i+j) * det(m[-i, -j, drop = FALSE])
  adjoint <- function(m) {
    cofactor_matrix <- outer(seq_len(n), seq_len(n), 
                             Vectorize(function(i, j) cofactor(m, i, j)))
    t(cofactor_matrix)
  }
  
  1 / det(m) * adjoint(m)
}

bmark <- function(m, n) {
  df <- microbenchmark::microbenchmark("hand-rolled" = inv(m), 
                                       "built-in" = solve(m))
  df <- as.data.frame(df)
  df$n <- n
  df
}

n <- 2:30
mats <- n %>% map(~matrix(runif(.^2), ., .))
result <- map2_df(mats, n, bmark)

ggplot(result, aes(x = n, y = time / 1e6, color = expr)) + 
  geom_jitter(alpha = 0.25) +
  geom_smooth(aes(group = expr), color = "black") +
  scale_y_log10() +
  theme_bw() +
  theme(legend.position = c(0, 1), legend.justification = c(0, 1),
        legend.background = element_blank()) +
  annotation_logticks(sides = "l") +
  labs(x = "\nMatrix size", y = NULL, color = NULL,
       title = "Benchmarking built-in and hand-rolled matrix inversion in R", 
       subtitle = "Microseconds to calculate inverse")