Simulation-based Evaluation of Numerical Algorithms forEnhancing Geostatistical Model Computation using MCMC: A Case Study in Disease Mapping

Background Bayesian geostatistical models incorporate spatial dependence through covariance matrices of location-specific random effects. In large spatial datasets, repeated manipulation of these covariance matrices represents a major computational bottleneck. Sparse matrix representations offer a potential strategy for improving computational feasibility, although their impact may vary across numerical algorithms. Methods We evaluated the computational behaviour of three numerical strategies Cholesky decomposition, Conjugate Gradient (CG), and Gibbs-based inversion under varying covariance matrix sizes (300, 600, and 1,000 locations) and sparsity levels (20%, 50%, and 100% zero-valued off-diagonal elements). Performance was assessed in terms of computational time and approximation accuracy relative to a deterministic reference solution. Results The results indicate that computational efficiency is strongly influenced by both matrix sparsity and algorithmic structure Cholesky decomposition serves as a deterministic benchmark, providing exact solutions for symmetric positive definite matrices with machine-level precision. The Conjugate Gradient method demonstrates rapid deterministic convergence, attaining near-zero Root Mean Square Error (RMSE) within a small fraction of the computational time required by direct factorisation, with convergence speed influenced primarily by matrix conditioning. In contrast, Gibbs sampling exhibits stochastic convergence, characterised by gradual RMSE reduction and transient variability due to Monte Carlo fluctuations. Extended iterations reduce sampling variability and stabilise accuracy, consistent with asymptotic sampling theory. Our results highlight the fundamental distinction between deterministic solvers, which optimise numerical precision and efficiency, and stochastic simulation-based approaches, which prioritise posterior uncertainty quantification. Conclusion Sparse covariance structures substantially improve computational feasibility in large-scale Bayesian geostatistical modelling. This comparative assessment provides practical guidance for solver selection in Gaussian geostatistical modelling: CG is the optimal strategy for computationally efficient deterministic recovery, whereas Gibbs sampling retains methodological relevance for fully Bayesian inference.