Performance & Scalability Benchmarks

Executive Summary

Small Area Estimation often involves large administrative registries or extensive spatial networks with hundreds or thousands of domains. Traditional implementations in R often rely on interpreted loops, pure R Fisher-scoring routines, or large memory allocations for $D \times D$ dense covariance matrices.

fastsae re-architects core estimation algorithms in compiled C++ using RcppArmadillo and native OpenMP multi-threading, delivering:

Up to 364x faster than sae and 12,300x faster than emdi for Fay-Herriot models at $n = 1,000$ .
Up to 80x faster than sae for Spatial Fay-Herriot models at $n = 1,000$ .
Up to 31x faster than sae for Spatio-Temporal models at $n = 1,000$ .
Massive RAM reduction: Peak memory footprint stays under 25 MB where existing packages require hundreds of megabytes or several gigabytes.

Benchmark Results Table

The following benchmarks were conducted on simulated datasets across domain sizes ranging from $n = 30$ to $n = 1,000$ (with 5 auxiliary covariates):

Metric	`fastsae`	`sae` (Molina & Rao)	`emdi` (Kreutzmann et al.)
Mean Time (EBLUP FH)	0.0015 s	0.291 s	9.64 s
Mean Time (Spatial FH)	0.165 s	12.60 s	8.69 s
Mean Time (Spatio-Temporal FH)	12.5 s	353.0 s	-
Peak Memory (EBLUP FH)	0.055 MB	16.3 MB	824 MB
Peak Memory (Spatial FH)	5.15 MB	408 MB	824 MB
Peak Memory (Spatio-Temporal FH)	0.289 MB	7,822 MB	-
Speedup at n = 1,000 (FH)	Baseline	~364x slower	~12,300x slower

Interactive Benchmark Explorer

Use the interactive controls below to compare execution time, RAM consumption, and iterations per second across sample sizes:

Execution time (median, seconds) — FH · log scale

Speedup Factor (fastsae vs competitors) — FH

Architectural Insights: Why is fastsae so Fast?

Compiled C++ Linear Solvers: Instead of interpreting nested loops in R, fastsae implements Fisher-scoring parameter search and Woodbury identity matrix inversions in Armadillo C++, directly leveraging optimized BLAS/LAPACK routines.
Zero-Copy Matrix Operations: Memory allocations for large intermediate structures ( $V$ , $V^{-1}$ , $P$ ) are avoided or reused across Fisher iterations rather than reallocated on the heap.
OpenMP Multi-Threaded Bootstrap: Bootstrap resampling runs natively in parallel across CPU cores using #pragma omp parallel for, avoiding the serialization overhead of R worker processes (parallel::makeCluster / foreach).
Woodbury Identity for Panel Data: In eblup_stfh, inversion of the block $(DT \times DT)$ covariance matrix $V$ is reduced to operations on individual $D \times D$ and $T \times T$ blocks via Kronecker and Woodbury decomposition, transforming an $O((DT)^3)$ bottleneck into scalable $O(D^3 + T^3)$ steps.
Integrated Nested Laplace Approximations (INLA): For complex hierarchical generalized and spatio-temporal models (hb_area), fastsae utilizes INLA to compute analytical posterior margins directly from sparse Gaussian Markov Random Fields, bypassing Monte Carlo sampling entirely.

Bayesian Spatio-Temporal Benchmark: `fastsae` (INLA) vs `tipsae` (Stan MCMC)

To assess Bayesian small area estimation performance, fastsae::hb_area() was benchmarked against tipsae::fit_sae(), the state-of-the-art Stan MCMC package for spatio-temporal Beta small area models.

The test used the official tipsae panel dataset: emilia (38 health districts in Emilia-Romagna over 5 years, $N = 190$ domains $\times$ years) with spatial polygon contiguity matrix $W$ from emilia_shp.

Metric	`tipsae::fit_sae` (Stan MCMC)	`fastsae::hb_area` (INLA)	Advantage
Model Specification	Besag ICAR + Domain RW(1)	`spatial = "besag"`, `temporal = "rw1"`, `st_interaction = "domain-specific"`	Exact structural equivalence
Computational Engine	Hamiltonian Monte Carlo (NUTS Stan)	Integrated Nested Laplace Approximation	Analytical & deterministic
Execution Time	15.8 s (1 chain, 200 iter) / ~120 s (4 chains)	1.93 s (`simplified.laplace`)	~9x to 50x+ faster
Pearson Correlation ( $r$ )	Baseline	0.9893	Near-identical point estimates
Mean Absolute Error (MAE)	Baseline	0.00304	Negligible numerical error
Convergence Overhead	Requires $\hat{R} < 1.05$ checks, warmup, and tuning	None (closed-form Laplace expansions)	Instant convergence

Interactive Beta SAE Benchmark Explorer ( $n = 30$ to $n = 1,000$ )

Use the interactive controls below to compare runtime, memory footprint, and speedup factors between fastsae::hb_area and tipsae::fit_sae across Standard Beta, Spatial Beta (Besag ICAR), and Spatio-Temporal Beta (Besag + RW1) models:

Execution time (median, seconds) — Standard Beta · log scale

Speedup Factor & Efficiency Multiplier (fastsae vs tipsae) — Standard Beta

Inferential Equivalence & Validation Metrics ( $n = 1,000$ )

To establish that the substantial computational accelerations achieved by fastsae::hb_area preserve full inferential validity relative to Stan MCMC (tipsae::fit_sae), summary concordance and error metrics were evaluated across $n = 1,000$ domains/units for Standard Beta, Spatial Beta (Besag ICAR), and Spatio-Temporal Beta models:

Dimension	Metric	Beta	Spatial Beta	Spatio-Temporal Beta	Statistical Implication
Point Estimates (EBP)	Pearson Correlation ( $r$ )	0.99934	0.99891	0.99614	Near-perfect linear agreement between INLA and MCMC
	Spearman Rank Correlation ( $\rho$ )	0.99921	0.99875	0.99517	Consistent domain priority and rank ordering
	Mean Absolute Error (MAE)	0.01969	0.01456	0.01200	Average estimation discrepancy $< 0.02$
	Root Mean Squared Difference (RMSD)	0.02369	0.01713	0.01458	Negligible $L_2$ deviation across domains
Uncertainty / MSE	MSE Pearson Correlation ( $r_{\text{MSE}}$ )	0.96315	0.89254	0.67957	Preserved area-level precision ordering
	Mean Absolute Difference (MSE)	0.00018	0.00012	0.00025	Virtually identical error variances

Finite Population Simulation & Ground Truth Validation

For comprehensive design-based simulation studies evaluating parameter recovery and empirical properties under repeated survey sampling from finite populations ( $D = 60$ and national-scale $D = 500$ ), see the dedicated article: Finite Population Simulation & Validation.