Benchmarks

Each program runs three ways: as plain R on GNU R, as the same plain R on ry’s runtime, and with inline annotations checked in strict mode on ry’s runtime. All three compute the same result:

program	time			ry vs GNU R	peak memory
program	GNU R, untyped	ry, untyped	ry, typed and strict	ry vs GNU R	GNU R, untyped	ry, untyped	ry, typed and strict
`array_dot`	0.278 s	0.043 s	0.047 s	6.5x	81 MiB	103 MiB	124 MiB
`array_filter`	0.246 s	0.043 s	0.038 s	5.7x	81 MiB	102 MiB	123 MiB
`array_map`	0.385 s	0.054 s	0.067 s	7.1x	94 MiB	123 MiB	141 MiB
`array_sum`	0.175 s	0.033 s	0.033 s	5.3x	81 MiB	102 MiB	123 MiB
`closures`	1.100 s	0.570 s	0.607 s	1.9x	99 MiB	145 MiB	148 MiB
`drawdown`	0.335 s	0.074 s	0.079 s	4.5x	73 MiB	133 MiB	141 MiB
`fib`	1.312 s	0.010 s	0.010 s	131.2x	54 MiB	43 MiB	43 MiB
`log_rank`	0.128 s	0.137 s	0.136 s	0.9x	107 MiB	56 MiB	58 MiB
`loop_sum`	0.312 s	0.069 s	0.068 s	4.5x	50 MiB	82 MiB	106 MiB
`vasicek`	0.055 s	0.040 s	0.040 s	1.4x	106 MiB	64 MiB	70 MiB
`vector_arith`	0.290 s	0.242 s	0.249 s	1.2x	109 MiB	157 MiB	164 MiB

Code that spends its time in R itself gains the most. Recursion (fib) runs over 100 times faster, and loops over the elements of a vector, or carrying scalar state through one, run 4 to 7 times faster. Code that spends its time in vectorized builtins gains little, because GNU R already runs those in C: vector_arith, vasicek and log_rank are within a third of GNU R, and log_rank is slightly slower. The typed column stays within a few percent of the untyped one, so today types make code safer, not faster.

Memory goes both ways. ry’s process peaks higher than GNU R’s in most programs, up to 1.8 times as high in drawdown, and the typed programs peak another 2 to 24 MiB higher, while fib, log_rank and vasicek peak lower on ry.

The programs

Each program exists twice in tools/bench/why/: as plain R, and typed. Open one to compare the two versions; the tab you choose applies to every program.

array_dot 6.5x

Element loops over a million doubles: a dot product of `x[[i]] * y[[i]]`. The typed program stores the same values in Arrays, which hold no NA and no attributes.

Plain R
Typed

array_dot <- function(x, y, times) {
  total <- 0
  n <- length(x)
  for (k in seq_len(times)) {
    for (i in seq_len(n)) total <- total + x[[i]] * y[[i]]
  }
  total
}
x <- seq(0, 1, length.out = 1000000L)
y <- rev(x)
measure("array_dot", function() array_dot(x, y, 5L))

array_dot <- function(x: Array<double>, y: Array<double>, times: integer) -> double {
  total <- 0
  n <- List.len(x)
  for (k in seq_len(times)) {
    for (i in seq_len(n)) total <- total + x[[i]] * y[[i]]
  }
  total
}
values: double[] <- seq(0, 1, length.out = 1000000L)
reversed: double[] <- rev(values)
x <- as.Array(values)
y <- as.Array(reversed)
measure("array_dot", function() array_dot(x, y, 5L))

array_filter 5.7x

Element loops over a million doubles: a filter counting the items above one half. The typed program stores the same values in Arrays, which hold no NA and no attributes.

Plain R
Typed

array_filter <- function(x, y, times) {
  total <- 0
  n <- length(x)
  for (k in seq_len(times)) {
    for (i in seq_len(n)) if (x[[i]] > 0.5) total <- total + 1
  }
  total
}
x <- seq(0, 1, length.out = 1000000L)
y <- rev(x)
measure("array_filter", function() array_filter(x, y, 5L))

array_filter <- function(x: Array<double>, y: Array<double>, times: integer) -> double {
  total <- 0
  n <- List.len(x)
  for (k in seq_len(times)) {
    for (i in seq_len(n)) if (x[[i]] > 0.5) total <- total + 1
  }
  total
}
values: double[] <- seq(0, 1, length.out = 1000000L)
reversed: double[] <- rev(values)
x <- as.Array(values)
y <- as.Array(reversed)
measure("array_filter", function() array_filter(x, y, 5L))

array_map 7.1x

Element loops over a million doubles: a map writing `out[[i]] <- x[[i]] * 2 + 1` into a preallocated vector. The typed program stores the same values in Arrays, which hold no NA and no attributes.

Plain R
Typed

array_map <- function(x, y, times) {
  total <- 0
  n <- length(x)
  out <- numeric(n)
  for (k in seq_len(times)) {
    for (i in seq_len(n)) out[[i]] <- x[[i]] * 2 + 1
    total <- total + sum(out)
  }
  total
}
x <- seq(0, 1, length.out = 1000000L)
y <- rev(x)
measure("array_map", function() array_map(x, y, 5L))

array_map <- function(x: Array<double>, y: Array<double>, times: integer) -> double {
  total <- 0
  n <- List.len(x)
  out <- as.Array(numeric(n))
  for (k in seq_len(times)) {
    for (i in seq_len(n)) out[[i]] <- x[[i]] * 2 + 1
    total <- total + List.sum(out)
  }
  total
}
values: double[] <- seq(0, 1, length.out = 1000000L)
reversed: double[] <- rev(values)
x <- as.Array(values)
y <- as.Array(reversed)
measure("array_map", function() array_map(x, y, 5L))

array_sum 5.3x

Element loops over a million doubles: a sum of `x[[i]]`. The typed program stores the same values in Arrays, which hold no NA and no attributes.

Plain R
Typed

array_sum <- function(x, y, times) {
  total <- 0
  n <- length(x)
  for (k in seq_len(times)) {
    for (i in seq_len(n)) total <- total + x[[i]]
  }
  total
}
x <- seq(0, 1, length.out = 1000000L)
y <- rev(x)
measure("array_sum", function() array_sum(x, y, 5L))

array_sum <- function(x: Array<double>, y: Array<double>, times: integer) -> double {
  total <- 0
  n <- List.len(x)
  for (k in seq_len(times)) {
    for (i in seq_len(n)) total <- total + x[[i]]
  }
  total
}
values: double[] <- seq(0, 1, length.out = 1000000L)
reversed: double[] <- rev(values)
x <- as.Array(values)
y <- as.Array(reversed)
measure("array_sum", function() array_sum(x, y, 5L))

closures 1.9x

Closures and lapply: small closures called in a loop, then lapply and vapply over a sequence.

Plain R
Typed

closures <- function(n) {
  add <- function(a, b) a + b
  scale <- 3
  times_scale <- function(a) a * scale
  s <- 0
  for (i in 1:n) s <- add(s, times_scale(i))
  squares <- lapply(seq_len(n), function(x) x * (x + 0.5))
  halves <- vapply(seq_len(n), function(x) x / 2, numeric(1))
  squared <- unlist(squares)
  s + sum(squared) + sum(halves)
}
measure("closures", function() closures(300000L))

closures <- function(n: integer) -> double {
  add <- function(a: double, b: double) -> double a + b
  scale <- 3
  times_scale <- function(a: double) -> double a * scale
  s <- 0
  for (i in 1:n) s <- add(s, times_scale(i))
  squares <- lapply(seq_len(n), function(x: integer) -> double x * (x + 0.5))
  halves <- vapply(seq_len(n), function(x: integer) -> double x / 2, numeric(1))
  squared: double[] <- unlist(squares)
  s + sum(squared) + sum(halves)
}
measure("closures", function() closures(300000L))

drawdown 4.5x

A hot path from the finance code base: maximum drawdown and EWMA volatility, scalar state carried through one pass over a return series.

Plain R
Typed

risk_report <- function(returns, lambda = 0.94) {
  wealth <- 1
  peak <- 1
  worst <- 0
  variance <- returns[1] * returns[1]
  for (i in seq_along(returns)) {
    r <- returns[i]
    wealth <- wealth * (1 + r)
    if (wealth > peak) peak <- wealth
    drawdown <- (peak - wealth) / peak
    if (drawdown > worst) worst <- drawdown
    variance <- lambda * variance + (1 - lambda) * r * r
  }
  c(max_drawdown = worst, volatility = sqrt(variance * 252))
}
returns <- sin(seq_len(2000000L) * 0.37) * 0.01 + 0.0002
measure("drawdown", function() sum(risk_report(returns)))

risk_report <- function(returns: double[], lambda: double = 0.94) -> double[] {
  wealth <- 1
  peak <- 1
  worst <- 0
  variance <- returns[1] * returns[1]
  for (i in seq_along(returns)) {
    r <- returns[i]
    wealth <- wealth * (1 + r)
    if (wealth > peak) peak <- wealth
    drawdown <- (peak - wealth) / peak
    if (drawdown > worst) worst <- drawdown
    variance <- lambda * variance + (1 - lambda) * r * r
  }
  c(max_drawdown = worst, volatility = sqrt(variance * 252))
}
returns <- sin(seq_len(2000000L) * 0.37) * 0.01 + 0.0002
measure("drawdown", function() sum(risk_report(returns)))

fib 131.2x

Recursion: about 2.7 million calls of a small function.

Plain R
Typed

fib <- function(n) if (n < 2) n else fib(n - 1) + fib(n - 2)
measure("fib", function() fib(30))

fib <- function(n: double) -> double if (n < 2) n else fib(n - 1) + fib(n - 2)
measure("fib", function() fib(30))

log_rank 0.9x

A hot path from the biostatistics code base: the log-rank statistic, a loop over event times with whole-vector comparisons in its body.

Plain R
Typed

log_rank <- function(data) {
  times <- sort(unique(data$time[data$event == 1L]))
  observed <- 0
  expected <- 0
  variance <- 0
  for (t in times) {
    risk <- data$time >= t
    n <- sum(risk)
    n1 <- sum(risk & data$arm == "treated")
    d <- sum(data$time == t & data$event == 1L)
    d1 <- sum(data$time == t & data$event == 1L & data$arm == "treated")
    observed <- observed + d1
    expected <- expected + d * n1 / n
    if (n > 1) variance <- variance + d * (n1 / n) * (1 - n1 / n) * (n - d) / (n - 1)
  }
  (observed - expected)^2 / variance
}
subjects <- 4000L
data <- list(
  time = (seq_len(subjects) * 7919L) %% 1000L + 1L,
  event = as.integer(seq_len(subjects) %% 3L != 0L),
  arm = ifelse(seq_len(subjects) %% 2L == 0L, "treated", "control")
)
measure("log_rank", function() log_rank(data))

log_rank <- function(data: list{time: integer[], event: integer[], arm: character[]}) -> double {
  times <- sort(unique(data$time[data$event == 1L]))
  observed <- 0
  expected <- 0
  variance <- 0
  for (t in times) {
    risk <- data$time >= t
    n <- sum(risk)
    n1 <- sum(risk & data$arm == "treated")
    d <- sum(data$time == t & data$event == 1L)
    d1 <- sum(data$time == t & data$event == 1L & data$arm == "treated")
    observed <- observed + d1
    expected <- expected + d * n1 / n
    if (n > 1) variance <- variance + d * (n1 / n) * (1 - n1 / n) * (n - d) / (n - 1)
  }
  (observed - expected)^2 / variance
}
subjects <- 4000L
data <- list(
  time = (seq_len(subjects) * 7919L) %% 1000L + 1L,
  event = as.integer(seq_len(subjects) %% 3L != 0L),
  arm = ifelse(seq_len(subjects) %% 2L == 0L, "treated", "control")
)
measure("log_rank", function() log_rank(data))

loop_sum 4.5x

A numeric loop: a scalar double accumulator in a counted loop.

Plain R
Typed

loop_sum <- function(n) {
  s <- 0
  for (i in 1:n) s <- s + i * 0.5
  s
}
measure("loop_sum", function() loop_sum(10000000L))

loop_sum <- function(n: integer) -> double {
  s <- 0
  for (i in 1:n) s <- s + i * 0.5
  s
}
measure("loop_sum", function() loop_sum(10000000L))

vasicek 1.4x

A simulation kernel from the finance code base: Vasicek interest-rate paths, one matrix row per step.

Plain R
Typed

vasicek_paths <- function(paths = 200L, steps = 2000L, r0 = 0.03, kappa = 0.5, theta = 0.04, sigma = 0.01, dt = 1 / 12, seed = 9L) {
  set.seed(seed)
  rates <- matrix(r0, nrow = steps + 1L, ncol = paths)
  for (t in seq_len(steps)) {
    shock <- stats::rnorm(paths)
    previous <- rates[t, ]
    rates[t + 1L, ] <- previous + kappa * (theta - previous) * dt + sigma * sqrt(dt) * shock
  }
  rates
}
measure("vasicek", function() sum(vasicek_paths()))

vasicek_paths <- function(paths: integer = 200L, steps: integer = 2000L, r0: double = 0.03, kappa: double = 0.5, theta: double = 0.04, sigma: double = 0.01, dt: double = 1 / 12, seed: integer = 9L) -> matrix {
  set.seed(seed)
  rates <- matrix(r0, nrow = steps + 1L, ncol = paths)
  for (t in seq_len(steps)) {
    shock <- stats::rnorm(paths)
    previous: double[] <- rates[t, ]
    rates[t + 1L, ] <- previous + kappa * (theta - previous) * dt + sigma * sqrt(dt) * shock
  }
  rates
}
measure("vasicek", function() sum(vasicek_paths()))

vector_arith 1.2x

Vector code: whole-vector arithmetic on a million doubles, repeated in a loop.

Plain R
Typed

vector_arith <- function(x, times) {
  total <- 0
  for (k in seq_len(times)) {
    y <- x * 2 + x / 3 - 1
    total <- total + sum(y)
  }
  total
}
x <- seq(0, 1, length.out = 1000000L)
measure("vector_arith", function() vector_arith(x, 20L))

vector_arith <- function(x: double[], times: integer) -> double {
  total <- 0
  for (k in seq_len(times)) {
    y <- x * 2 + x / 3 - 1
    total <- total + sum(y)
  }
  total
}
x: double[] <- seq(0, 1, length.out = 1000000L)
measure("vector_arith", function() vector_arith(x, 20L))

harness measure()

Shared by every benchmark: runs `run()` once to warm up (GNU R's JIT compiles, ry tiers up), then three timed times, and prints the checksum and the best time in seconds.

Plain R
Typed

measure <- function(name, run) {
  value <- run()
  best <- Inf
  for (k in 1:3) {
    started <- proc.time()[[3]]
    again <- run()
    finished <- proc.time()[[3]]
    best <- min(best, finished - started)
    if (!identical(again, value)) stop("the benchmark is not deterministic")
  }
  cat(sprintf("%s %s %.4f\n", name, format(value, digits = 15), best))
}

measure <- function(name: character, run: fn() -> double) -> NULL {
  value <- run()
  best: double <- Inf
  for (k in 1:3) {
    started: double <- proc.time()[[3]]
    again <- run()
    finished: double <- proc.time()[[3]]
    best <- min(best, finished - started)
    if (!identical(again, value)) stop("the benchmark is not deterministic")
  }
  cat(sprintf("%s %s %.4f\n", name, format(value, digits = 15), best))
}

How it was measured

Each script times its own work with the harness above, so a time leaves out startup and, for the typed column, the static check. The harness runs the program once to warm up (GNU R compiles it to bytecode, and ry compiles hot code to machine code), then three more times, and keeps the fastest. just bench-why --write runs every program twice per engine, alternating the engines so they never share the CPU, and stores the fastest of those six runs in tools/bench/why/results.json, which every table on this site reads. The fastest run is reported rather than the median because other load on a machine only ever slows a run down.

Peak memory is the peak resident set size (RSS) of the whole process, R’s startup included, as the operating system reports it when the process exits. Each cell keeps the largest peak of that engine’s runs, so no engine gets credit for a lucky run.

Measured on Intel(R) Xeon(R) Processor @ 2.80GHz, 4 cores, Linux 6.18.44-fc-v185, against GNU R 4.3.3.