Efficient computation of areal mean prediction variance in Gaussian Processes

External question — Cross Validated Stack Exchange Author: Ahmad Awad Original post: https://stats.stackexchange.com/questions/674678 License: CC BY-SA 4.0 — https://creativecommons.org/licenses/by-sa/4.0/ Adaptation: HTML converted to plain text; contact email addresses removed. Block kriging and gaussian process areal averaging In geostatistics, block kriging provides formulas for predicting the average value $$ Z_B = |B|^{-1} \int_B Z(s) \, ds $$ over a domain $B$ and its variance. These involve integrals of the covariance function $C(h)$ over $B$ , but in practice they are computed by discretizing $B$ into $m$ points. In Gaussian Process (GP) regression, we have data $Z$ at locations $\{s_i\}_{i=1}^n$ and a posterior GP with mean function $\mu(\cdot)$ and covariance function $k(\cdot,\cdot)$ . To estimate $Z_B$ , we typically discretize $B$ into $m$ points $\{t_j\}_{j=1}^m$ . The posterior at these points is $$ Z^\ast \sim \mathcal{N}(\hat{\mu}, \Sigma), $$ with $$ \hat{\mu} = K_\ast^T K^{-1} Z, \quad \Sigma = K_{\ast\ast} - K_\ast^T K^{-1} K_\ast. $$ The discrete areal mean is $$ \hat{Z}_B = m^{-1} \sum_{j=1}^m Z(t_j) $$ with variance $$ \operatorname{Var}_{\mathrm{disc}} = m^{-2} \mathbf{1}^T \Sigma \mathbf{1} \quad (1) $$ Direct computation of (1) is $O(m^3)$ due to the need for $\Sigma$ , which is prohibitive for fine grids. I am wondering what would be efficient and statistically sound approximators for large $m$ ? I have a few thoughts that would like your inputs and would appreciate other suggestions as well. 1. Monte Carlo (MC) over the grid Draw $L$ joint samples $$ z^{(l)} \sim \mathcal{N}(\hat{\mu}, \Sigma), \quad l=1,\dots,L, $$ compute each areal mean $$ z_B^{(l)} = m^{-1} \sum_{j=1}^m z_j^{(l)}, $$ and then estimate $$ \widehat{\operatorname{Var}}_{\mathrm{MC}} = \frac{1}{L-1} \sum_{l=1}^L (z_B^{(l)} - \bar{z}_B)^2. $$ Question 1: How efficient and feasible is this approach? What strategies could make it more efficient? and what are good heuristics for $L$ ? 2. Single large MC sample as a direct quadratic-form estimator Draw one joint sample with size $n$ $$ z \sim \mathcal{N}(\hat{\mu}, \Sigma) $$ and compute the single-sample statistic $$ q = n^{-2} (z - \hat{\mu})^T (z - \hat{\mu}) = n^{-2} \sum_{i,j} (z_i - \hat{\mu}_i)(z_j - \hat{\mu}_j) $$ Question 2: Can $q$ serve as an estimator for $\operatorname{Var}_{\mathrm{disc}}$ ? What are its bias and variance? How does its statistical efficiency compare to the Monte Carlo variance estimator $\widehat{\operatorname{Var}}_{\mathrm{MC}}$ that uses multiple samples? I seek understanding of the computational-statistical trade-offs between these strategies when $m$ is large ( $\approx 10^5$ – $10^9$ ). Which approach is most promising in terms of error control and computational cost for practical spatial averaging with Gaussian Processes?