We propose Deep Neural Coregionalization, a scalable framework for uncertainty-aware multivariate geostatistics. DNC models multivariate spatial effects through spatially varying latent factors and loadings, assigning deep Gaussian process (DGP) priors to both the factors and the entries of the loading matrix. This joint construction learns shared latent spatial structure together with response-specific, location-dependent mixing weights, enabling flexible nonlinear and space-dependent associations within and across variables. A key contribution is a variational formulation that makes the DGP to deep neural network (DNN) correspondence explicit: maximizing the DGP evidence lower bound (ELBO) is equivalent to training DNNs with weight decay and Monte Carlo (MC) dropout. This yields fast mini-batch stochastic optimization without Markov Chain Monte Carlo (MCMC), while providing principled uncertainty quantification through MC-dropout forward passes as approximate posterior draws, producing calibrated credible surfaces for prediction and spatial effect estimation. Across simulations, DNC is competitive with existing spatial factor models, particularly under strong nonstationarity and complex cross-dependence, while delivering substantial computational gains. In a multivariate environmental case study, DNC captures spatially varying cross-variable interactions, produces interpretable maps of multivariate outcomes, and scales uncertainty quantification to large datasets with orders-of-magnitude reductions in runtime.
Technological and computational advances have created data-rich settings that offer unprecedented opportunities to probe the complexity of large, spatially indexed datasets. Over the past decade, spatial models have played a central role in addressing increasingly complex questions in the environmental sciences. A particularly active frontier is the analysis of multivariate spatial data, where multiple outcomes are recorded at each location. In such settings, it is natural to posit dependence both within locations (among the co-measured variables) and across locations (spatial autocorrelation within each variable). As a motivating example in this article, consider a suite of spatially indexed spectral measures of vegetation activity which typically exhibits strong spatial structure, readily seen in empirical variograms and exploratory maps, while the variables themselves are mutually associated through shared biophysical processes. Analyzing each variable in isolation may recover its marginal spatial pattern yet discards cross-variable information, often degrading interpolation and prediction performance (see, e.g., [Wackernagel, 2003, Chiles and Delfiner, 2012, Cressie and Wikle, 2011]). These considerations strongly motivate joint modeling of multivariate spatial processes. Moreover, in many applications, including ours, the strength and form of the associations among variables vary over space. Addressing this reality requires models that allow spatially varying association structure between variables and deliver computationally efficient inference at scale, suitable for large numbers of locations.
A common entry point to multivariate spatial analysis is to posit a vector-valued latent spatial process, typically a multivariate Gaussian process, equipped with a matrixvalued cross-covariance that encodes both within-variable and cross-variable dependence across space. The literature on such cross-covariances is vast, so we highlight the principal construction paradigms that balance interpretability and flexibility. The most widely used is the linear model of coregionalization (LMC) [Schmidt andGelfand, 2003, Wackernagel, 2003], which represents a multivariate field as a linear combination of independent univariate latent processes, producing cross-covariances as sums of separable components. A second family are convolution methods, which build cross-covariances by convolving shared and variable-specific kernels with latent white-noise or parent processes [Gaspari andCohn, 1999, Majumdar andGelfand, 2007]. A third idea is the latent-dimension approach, which assigns each variable a coordinate in an auxiliary space and then applies a univariate stationary kernel in the augmented domain, yielding valid stationary cross-covariances by construction [Apanasovich and Genton, 2010]. Finally, multivariate Matérn families provide interpretable, practice-friendly specifications wherein each margin enjoys a Matérn form and cross-parameters are constrained to ensure positive definiteness [Gneiting et al., 2010, Genton andKleiber, 2015].
While much of this work targets stationary dependence, where cross-covariances depend only on inter-location distances, modern applications frequently exhibit spatially varying associations tied to geography, ecology, or physics. Empirical evidence from ecology and environmental science (e.g., [Diez and Pulliam, 2007, Ovaskainen et al., 2010, Waddle et al., 2010, Coombes et al., 2015, Guha et al., 2024]) underscores that nonstationary cross-correlation can reveal latent drivers, sharpen scientific understanding, and materially improve prediction at new locations. Methodologically, non-stationarity has been introduced by allowing Matérn parameters (variance, range, smoothness) to vary over space [Kleiber and Nychka, 2012], by extending kernel convolutions to yield spatially varying matrix-valued covariances [Calder, 2008, Majumdar andGelfand, 2007], and by adapting the latent-dimension framework to nonstationary settings [Bornn et al., 2012]. Related hierarchical formulations, such as the matrix-variate Wishart process [Gelfand et al., 2004] and Lagrangian frameworks that induce directionally stronger cross-covariances [Salvaña et al., 2023], offer additional flexibility for modeling spatially varying covariance structures.
Scaling these multivariate spatial models to large n (locations) remains challenging. Recent advances introduce computationally efficient Gaussian-process surrogates-low-rank/predictiveprocess, sparse/nearest-neighbor, and multi-resolution schemes-to deliver nonstationary cross-covariances at scale [Guhaniyogi et al., 2013a, Guhaniyogi, 2017, Zhang and Banerjee, 2022a]. These approaches substantially broaden the reach of multivariate spatial modeling, yet truly massive sample sizes can still strain computation, especially when flexible spatially-dependent association is required across many variables and locations.
Expanding upon the nonstationary Linear Model of C
This content is AI-processed based on open access ArXiv data.