Imperative vs. Functional Factorials

Functions

These functions are pretty self explanatory - the only interesting ones are using the purrr package comparing map, walk, and reduce.

recurser_acc <- function(x, acc = 1) {
  if (x == 0) acc
  else recurser_acc(x-1, acc = acc * x)
}
  
recurser <- function(x) {
  if (x == 0) 1
  else x * recurser(x-1)
}

forlooper <- function(x) {
  begin <- 1
  for (num in x:1) begin <- num * begin
  begin
}

maplooper <- function(x) {
  begin <- 1
  purrr::map(x:1, function(num) begin <<- num * begin)
  begin
}

walker <- function(x) {
  begin <- 1
  purrr::walk(x:1, function(num) begin <<- num * begin)
  begin
}

reducer <- function(x) purrr::reduce(x:1, prod)

benchmark <- function(n) {
  microbenchmark::microbenchmark(
    factorial = factorial(n),
    forlooper = forlooper(n),
    recurser = recurser(n),
    recurser_acc = recurser_acc(n),
    maplooper = maplooper(n),
    walker = walker(n),
    reducer = reducer(n),
    times = 1000,
    unit = "us"
  )
}

And here are the preliminary tests. Box plots aren’t always the answer, clearly.

results <- purrr::map(c(5, 100, 1000), function(n) {
  benchmark(n) |>
    dplyr::mutate(n_eval = n)
}) |> dplyr::bind_rows()

library(ggplot2)
results |> 
  dplyr::mutate(norm_time = time / min(time), .by = n_eval) |>
  ggplot(aes(
    expr, norm_time,
    color = as.character(n_eval)
  )) +
  geom_boxplot() +
  # scale_y_continuous(limits = c(0, 6000)) +
  theme_bw()

So maybe we should plot it in in a different way. Let’s also plot the geometric mean to mitigate the effect of outliers.

results |> 
  dplyr::summarize(
    norm_time = exp(mean(log(time))),
    .by = c(n_eval, expr)
  ) |>
  dplyr::mutate(norm_time = norm_time / max(norm_time), .by = n_eval) |>
  ggplot(aes(
    expr, norm_time, 
    color = as.character(n_eval), 
    shape = as.character(n_eval) 
  )) +
  geom_point() +
  facet_grid(rows = vars(n_eval))

I am surprised to see that the length of the factorial actually didn’t affect the hierarchy of the results! Of course the base function would be the fastest, but the speed of the for loop here flies in the face of the point of R. However, we are comparing basically un-vectorized functions to a for loop, which is kind of against the spirit. Really we would need to vectorize these functions and then calculate the factorials of a vector (which we will do below). I will leave out recurser_acc and walker since they aren’t really informative when comparing whole paradigms.

# Factorial obviously is unchanged (vectorized by default)

forlooper_vec <- function(x) {
  result <- numeric(length(x))
  for (i in seq(1, length(x))) result[i] <- forlooper(x[i])
  result
}

maplooper_vec <- function(x) {
  purrr::map_dbl(x, maplooper)
}

reducer_vec <- function(x) {
  purrr::map_dbl(x, reducer)
}

Ok, and now let’s do the same benchmark again but with these guys using a vector. This time we’ll change the length of the vector to see if that makes a difference more than the actual number we use.

smaller_benchmark <- function(n) {
  microbenchmark::microbenchmark(
    factorial = factorial(n),
    forlooper = forlooper_vec(n),
    maplooper = maplooper_vec(n),
    reducer = reducer_vec(n),
    times = 500,
    unit = "us"
  )
}

results_vec <-  c(1, 5, 20) |>
  purrr::map(function(x) rep( c(5, 100, 1000), x)) |>
  purrr::map(function(n) {
    smaller_benchmark(n) |>
      dplyr::mutate(n_eval = length(n))
  }) |>
  dplyr::bind_rows()

results_vec |> 
  dplyr::summarize(
    norm_time = exp(mean(log(time))),
    .by = c(n_eval, expr)
  ) |>
  dplyr::mutate(norm_time = norm_time / max(norm_time), .by = n_eval) |>
  ggplot(aes(
    expr, norm_time, 
    color = as.character(n_eval), 
    shape = as.character(n_eval) 
  )) +
  geom_point() +
  facet_grid(rows = vars(n_eval))

It looks like even in a vectorized scenario we still see the forlooper beating out the map and mapreduce algorithms! I’m honestly disappointed, I find the reduce function to be very satisfying. I’m also surprised by the consistent trend we see across the length of the input vector (i.e. the performance hit is proportional to the size of the input vector). However, people that like to do things themselves might be pleased about the performance of the for loop functions relative to the base function, meaning that simple for loops are comparable in performance to just straight C function calls - pretty cool!