Repository navigation
scmBayesPost v0.3.0
scmBayesPost v0.3.0 — Major performance release
The problem
On datasets with large numbers of treated units (J0 ~ 1879 in the
motivating application), the Gibbs samplers for the Bayesian probit
first stage (selection_probit_bayes) were taking approximately
21 seconds per iteration — around 6 hours for a 2,000-iteration
run. Profiling identified three dominant costs:
- Cholesky factorisation of a
[K*J0 x K*J0]matrix (e.g.1879 x 1879) at every iteration, despite this matrix being block-diagonal - Sparse matrix-vector product
X_block %*% beta_bdover a
[N_stacked x K*J0]matrix at every iteration - Per-iteration
data.tablemerge to alignnu_hatwithX_block
row ordering
What changed
The samplers now exploit the block-diagonal structure of X_block
throughout. Key changes:
- The
[K*J0 x K*J0]Cholesky is replaced byJ0independent
[K x K]factorisations — forK = 1(the common single-treatment
case) this collapses to scalar arithmetic with no matrix operations X_block %*% beta_bdis replaced by a unit-by-unit loop over
precomputed row ranges and cached column slices- The
data.tablemerge is replaced by a one-time integer index map - Truncated normal draws use a fast inverse-CDF sampler with
precomputed treated/untreated splits rhoandsigma2updates use direct vector operations rather than
sparse matrix products
Expected speedup
The theoretical reduction in the Cholesky step alone is from
O((K*J0)^3) to O(J0 * K^3) — for K=1, J0=1879 this is a
factor of 1879^2 ≈ 3.5 million. In practice the overall per-iteration
speedup depends on which operations remain dominant after the fix
(primarily X_fs %*% delta, a 292k × 162 dense matrix-vector product).
Multi-threaded BLAS is recommended for further acceleration:
RhpcBLASctl::blas_set_num_threads(parallel::detectCores(logical = FALSE))API unchanged
No changes to function signatures, arguments, or output. Existing code
runs without modification.
Installation
remotes::install_github("andrebonfrer/scmBayesPost@v0.3.0")