I think you saw on Twitter that I just released a clone of (hommage to?) margins to CRAN. I’ve been running some benchmarks and finding significant speed gains when unit_ses=TRUE. I’m not 100% sure, but I think there’s a very good chance that you could close that gap by implementing a super simple algebra trick.
The opportunity is that JVJ' requires computing a very wide matrix. But since we only need the diagonal to compute standard errors, we can use a computational trick to get 1000x speedup and much lighter memory footprint.
I’m not sure exactly where this could be implemented because I haven’t looked at your code base in a little while, but I thought I’d leave this here in case you or a contributor feels like looking into this.
To illustrate, I’ll do a similar calculation with model matrices instead of jacobians. First, let’s create a variance-covariance matrix and a large model matrix by stacking mtcars multiple times:
library(bench)
mod <- lm(mpg ~ factor(cyl) + hp + drat, mtcars)
mtbig <- do.call("rbind", lapply(1:500, \(x) mtcars))
mm <- model.matrix(mod, data = mtbig)
vcovmat <- vcov(mod)
Before running the benchmark, let’s convince ourselves that the two strategies yield the same results at 10 digits tolerance:
Hi @leeper,
I think you saw on Twitter that I just released a clone of (hommage to?)
marginsto CRAN. I’ve been running some benchmarks and finding significant speed gains whenunit_ses=TRUE. I’m not 100% sure, but I think there’s a very good chance that you could close that gap by implementing a super simple algebra trick.The opportunity is that
JVJ'requires computing a very wide matrix. But since we only need the diagonal to compute standard errors, we can use a computational trick to get 1000x speedup and much lighter memory footprint.I’m not sure exactly where this could be implemented because I haven’t looked at your code base in a little while, but I thought I’d leave this here in case you or a contributor feels like looking into this.
To illustrate, I’ll do a similar calculation with model matrices instead of jacobians. First, let’s create a variance-covariance matrix and a large model matrix by stacking
mtcarsmultiple times:Before running the benchmark, let’s convince ourselves that the two strategies yield the same results at 10 digits tolerance:
Run benchmark: