Benchmarking estimation speed against other packages

Comparing run times

The {logitr} package includes a runtimes data frame with the benchmark results. The tables below summarize the run times (in seconds) for each package and how many times slower each is relative to single-threaded {logitr} (logitr (1 cores)), which is used as the reference point for a fair, hardware-independent comparison.

library(logitr)
library(dplyr)
library(tidyr)
library(kableExtra) # For tables

times <- runtimes %>% select(-version)
logitr_time <- times %>%
    filter(package == "logitr (1 cores)") %>%
    rename(time_logitr = time_sec)
time_compare <- times %>%
    left_join(select(logitr_time, -package), by = "numDraws") %>%
    mutate(mult = round(time_sec / time_logitr, 1)) 
# Compare raw times (seconds)
time_compare %>%
    select(package, numDraws, time_sec) %>%
    mutate(time_sec = round(time_sec, 1)) %>%
    pivot_wider(names_from = numDraws, values_from = time_sec) %>%
    kbl()

Now compare how many times slower each is relative to single-threaded {logitr} (logitr (1 cores)), the fair, hardware-independent reference point:

time_compare %>%
    select(package, numDraws, mult) %>%
    pivot_wider(names_from = numDraws, values_from = mult) %>%
    kbl()

The figure below visualizes the run times. All of the figures in this article are generated by the data-raw/figures.R script in the package repository.

Note that the estimation times reflect both the cost per iteration (which grows with the number of draws) and the number of iterations each package’s optimizer takes to converge, which can itself vary with the number of draws. This is why some packages – particularly {apollo} – show run times that are not strictly increasing in the number of draws. The overall comparison is robust to this: {logitr} is dramatically faster across every draw count, and its advantage grows with more draws and more cores.

Scaling to large draw counts

The comparison above stops at 1,500 draws because that is roughly where the other packages become impractically slow. {logitr}, by contrast, scales comfortably to 10,000 draws, which can be useful when many random parameters demand a large simulation sample for a stable simulated log-likelihood. The figure below isolates {logitr} and shows how its full-fit estimation time grows with the number of draws for each of its three backends:

• cpu – the native R path (memory-safe, single-threaded),
• cpp (1 core) – the compiled backend on a single core, and
• cpp (10 cores) – the compiled backend on all 10 cores.

These results are stored in the runtimes_draws data frame and generated by the same data-raw/runtimes.R script, run on the same 10-core Apple M-series machine as the cross-package benchmark above.

runtimes_draws %>%
    mutate(time_sec = round(time_sec, 1)) %>%
    pivot_wider(names_from = numDraws, values_from = time_sec) %>%
    kbl()

The compiled backend is several times faster than the native R path at every draw count, and adding cores multiplies that advantage further – so even at 10,000 draws, a model that is essentially out of reach for the other packages estimates in a fraction of the time.

What the speed buys you: a stable log-likelihood

Speed is not the end goal – it is what makes enough draws affordable. The mixed logit log-likelihood is approximated by simulation, so any estimate carries simulation error that shrinks as the number of draws grows. Historical defaults are small: {mlogit} and {gmnl} default to 40 draws, Stata’s mixlogit to 50, and {logitr} itself defaulted to 50 prior to version 1.2.0. Those defaults date from when estimation was slow enough that more draws felt like a luxury, but they leave enough simulation noise in the log-likelihood to matter.

To make that noise visible, the same mixed logit model used in the benchmarks above was estimated repeatedly at increasing draw counts. At each draw count the model was estimated 10 times using randomized MLHS draws with a different seed each time, with every run starting from the same converged reference solution – so at a given draw count, the across-seed spread in the converged log-likelihood reflects simulation error: the noise in the simulated log-likelihood shifts the converged value, and at low draw counts it can distort the surface enough that some runs converge to different local optima entirely. The results ship with the package in the loglik_draws data frame, generated by the same data-raw/runtimes.R script (run on the same 10-core machine as above).

loglik_summary <- loglik_draws %>%
    group_by(numDraws) %>%
    summarise(
        min = min(logLik),
        max = max(logLik),
        range = max(logLik) - min(logLik),
        sd = sd(logLik),
        median_time_sec = median(time_sec)
    ) %>%
    mutate(across(min:median_time_sec, \(x) round(x, 1)))

kbl(loglik_summary)

At 50 draws the converged log-likelihood ranges over about 89 units across seeds – easily enough to flip a likelihood-based model comparison – while at 500 draws the spread is roughly 3.3 times smaller and keeps narrowing from there. (The narrowing is not perfectly monotonic – with few draws, an occasional run wanders to a different local optimum, which is itself a symptom of the simulation noise, but the trend is unmistakable.) And because of the compiled backend, those 500 draws cost only about 0.7 seconds on this model. This is why {logitr} defaults to numDraws = 500: at these speeds there is no longer a good reason to accept the simulation noise of a 40-50 draw default.