Recent variational Bayes methods for geospatial regression, proposed as an alternative to computationally expensive Markov chain Monte Carlo (MCMC) sampling, have leveraged Nearest Neighbor Gaussian processes (NNGP) to achieve scalability. Yet, these variational methods remain inferior in accuracy and speed compared to spNNGP, the state-of-the-art MCMC-based software for NNGP. We introduce spVarBayes, a suite of fast variational Bayesian approaches for large-scale geospatial data analysis using NNGP. Our contributions are primarily computational. We replace auto-differentiation with a combination of calculus of variations, closed-form gradient updates, and linear response corrections for improved variance estimation. We also accommodate covariates (fixed effects) in the model and offer inference on the variance parameters. Simulation experiments demonstrate that we achieve comparable accuracy to spNNGP but with reduced computational costs, and considerably outperform existing variational inference methods in terms of both accuracy and speed. Analysis of a large forest canopy height dataset illustrates the practical implementation of proposed methods and shows that the inference results are consistent with those obtained from the MCMC approach. The proposed methods are implemented in publicly available Github R-package spVarBayes.

Small area estimation models are essential for estimating population characteristics in regions with limited sample sizes, thereby supporting policy decisions, demographic studies, and resource allocation, among other use cases. The spatial Fay-Herriot model is one such approach that incorporates spatial dependence to improve estimation by borrowing strength from neighboring regions. However, this approach often requires substantial computational resources, limiting its scalability for high-dimensional datasets, especially when considering multiple (multivariate) responses. This paper proposes two methods that integrate the multivariate spatial Fay-Herriot model with spatial random effects, learned through variational autoencoders, to efficiently leverage spatial structure. Importantly, after training the variational autoencoder to represent spatial dependence for a given set of geographies, it may be used again in future modeling efforts, without the need for retraining. Additionally, the use of the variational autoencoder to represent spatial dependence results in extreme improvements in computational efficiency, even for massive datasets. We demonstrate the effectiveness of our approach using 5-year period estimates from the American Community Survey over all census tracts in California.