Recursive Nearest Neighbor Co-Kriging Models for Big Multiple Fidelity Spatial Data Sets

TL;DR

This paper introduces the recursive nearest neighbor co-kriging (RNNC) model for fast, scalable Bayesian inference in large, multi-fidelity spatial datasets, significantly reducing computation time while maintaining accuracy.

Contribution

The paper develops the RNNC model and two efficient inference methods, including a novel MCMC-free approach, enabling rapid analysis of massive spatial datasets.

Findings

Reduced computational time from hours to minutes.

Maintained high prediction accuracy with new methods.

Demonstrated effectiveness on satellite data.

Abstract

Big datasets are gathered daily from different remote sensing platforms. Recently, statistical co-kriging models, with the help of scalable techniques, have been able to combine such datasets by using spatially varying bias corrections. The associated Bayesian inference for these models is usually facilitated via Markov chain Monte Carlo (MCMC) methods which present (sometimes prohibitively) slow mixing and convergence because they require the simulation of high-dimensional random effect vectors from their posteriors given large datasets. To enable fast inference in big data spatial problems, we propose the recursive nearest neighbor co-kriging (RNNC) model. Based on this model, we develop two computationally efficient inferential procedures: a) the collapsed RNNC which reduces the posterior sampling space by integrating out the latent processes, and b) the conjugate RNNC, an MCMC free inference which significantly reduces the computational time without sacrificing prediction accuracy. An important highlight of conjugate RNNC is that it enables fast inference in massive multifidelity data sets by avoiding expensive integration algorithms. The efficient computational and good predictive performances of our proposed algorithms are demonstrated on benchmark examples and the analysis of the High-resolution Infrared Radiation Sounder data gathered from two NOAA polar orbiting satellites in which we managed to reduce the computational time from multiple hours to just a few minutes.

Tables

Table 1. Table 1: The estimation of parameters in nested input dataset, using sequential NNCGP, collapsed RNNC and conjugate RNNC models.

	True	Nested data-set
	values	Sequential NNCGP		Collapsed RNNC		Conjugate RNNC
${\beta }_1$	10	10.29	(9.93,10.57)	9.96	(9.60,10.32)	10.02
${\beta }_2$	1	0.77	(0.59,1.04)	0.87	(0.59,1.13)	0.82
${\sigma }_1^2$	4	3.55	(2.77,4.38)	3.46	(2.96,4.27)	3.15
${\sigma }_2^2$	1	0.81	(0.27, 2.05)	0.98	(0.43, 1.88)	0.79
$1/{\varphi }_1$	10	10.42	(8.15,13.47)	10.50	(8.59,13.90)	12.1
$1/{\varphi }_2$	10	14.96	(3.37, 20.29)	15.69	(5.92, 19.98)	19.6
${\gamma }_1$	1	0.99	(0.98,1.00)	0.99	(0.98,1.00)	0.99
${\tau }_1^2$	0.1	0.12	(0.10,0.14)	0.15	(0.10,0.19)	0.12
${\tau }_2^2$	0.05	0.07	(0.03,0.11)	0.10	(0.04,0.19)	0.16
$m$	10	-	-	-	-	-

Table 2. Table 2: Performance measures for the predictive ability of the sequential NNCGP model, collapsed RNNC model and conjugate RNNC model.

	Nested data-set
	Sequential NNCGP	Collapsed RNNC	Conjugate RNNC
RMSPE	0.63	0.69	0.72
NSME	0.77	0.75	0.71
CRPS	0.45	0.44	0.41
CVG(95%)	0.91	0.88	0.93
ALCI(95%)	1.93	1.92	2.54
Time(Hour)	4.5	4.7	0.08

Table 3. Table 3: The estimation of parameters in non-nested input dataset using sequential NNCGP, collapsed RNNC and conjugate RNNC models.

	True	Non-nested data-set
	values	Sequential NNCGP		Collapsed RNNC		Conjugate RNNC
${\beta }_1$	10	9.71	(9.36, 10.16)	9.97	(9.52,10.41)	9.71
${\beta }_2$	1	0.87	(0.39,1.36)	1.23	(0.24,2.19)	1.27
${\sigma }_1^2$	4	3.51	(2.71,4.52)	3.28	(3.02,3.72)	3.84
${\sigma }_2^2$	1	1.05	(0.18,2.31)	1.00	(0.64, 1.49)	1.19
$1/{\varphi }_1$	10	10.77	(8.07,13.91)	13.23	(9.93,15.85)	6.50
$1/{\varphi }_2$	10	12.61	(3.93,24.07)	16.01	(12.84, 19.92)	20.00
${\gamma }_1$	1	0.99	(0.98,1.05)	0.97	(0.94,0.99)	0.96
${\tau }_1^2$	0.1	0.13	(0.10,0.15)	0.10	(0.07,0.15)	0.11
${\tau }_2^2$	0.05	0.16	(0.04,0.23)	0.18	(0.03,0.29)	0.10
$m$	10	-	-	-	-	-

Table 4. Table 4: Performance measures for the predictive ability of the Sequential NNCGP model, collapsed RNNC model and conjugate RNNC model, in non-nested input.

	Non-nested data-set
	Sequential NNCGP	Collapsed RNNC	Conjugate RNNC
RMSPE	0.93	0.92	1.07
NSME	0.75	0.78	0.75
CRPS	0.61	0.61	0.65
CVG(95%)	0.93	0.97	0.95
ALCI(95%)	1.49	1.76	2.49
Time(Hour)	4.1	4.3	0.05

Equations

$$z_t(\mathbf{s})=y_t(\mathbf{s})+{ϵ}_t,$$

$$y_t(\mathbf{s})={\zeta }_{t-1}(\mathbf{s})y_{t-1}(\mathbf{s})+{\delta }_t(\mathbf{s}),$$

$${\delta }_t(\mathbf{s})={\mathbf{h}}_t^T(\mathbf{s}){\mathbf{\beta }}_t+w_t(\mathbf{s}),$$

$$z_t(\mathbf{s})$$

$$y_t(\mathbf{s})$$

$${\delta }_t(\mathbf{s})$$

$${\hat{y}}_{t-1}(\mathbf{s})|{\mathbf{Z}}_{t-1},{\hat{y}}_{t-2}(\mathbf{s}),{\mathbf{\theta }}_{t-1},{\mathbf{\beta }}_{t-1}\sim$$

$$L({\mathbf{Z}}_{1:T}|\cdot )$$

$$=N({\mathbf{Z}}_1|{\mathbf{h}}_1({\mathbf{S}}_1){\mathbf{\beta }}_1+{\mathbf{w}}_1,{\tau }_1\mathbf{I})\prod _{t=2}^{T}N({\mathbf{Z}}_t|{\zeta }_{t-1}({\mathbf{S}}_t)\circ {\hat{y}}_{t-1}({\mathbf{S}}_t)+{\mathbf{h}}_t({\mathbf{S}}_t){\mathbf{\beta }}_t+{\mathbf{w}}_t,{\tau }_t\mathbf{I}),$$

$$w_t({\mathbf{s}}_{t,1})$$

$$w_t({\mathbf{s}}_{t,i})$$

$$z_1({\mathbf{S}}_1)|{\mathbf{\beta }}_1,{\mathbf{\theta }}_1,{\tau }_1$$

$$z_t({\mathbf{S}}_t)|{\mathbf{\beta }}_t,{\mathbf{\theta }}_t,{\tau }_t,{\zeta }_t({\mathbf{S}}_t),{\hat{y}}_{t-1}({\mathbf{S}}_t)$$

$${\tilde{\mathbf{\Lambda }}}_t^{-1}={\tau }_t^{-2}\mathbf{I}-{\tau }_t^{-4}({\tilde{\mathbf{C}}}_t({\mathbf{\theta }}_t{)}^{-1}+{\tau }_t^{-2}\mathbf{I}{)}^{-1},$$

$$\text{det}({\tilde{\mathbf{\Lambda }}}_t)={\tau }_t^{2n}\text{det}({\tilde{\mathbf{C}}}_t({\mathbf{\theta }}_t))\text{det}({\tilde{\mathbf{C}}}_t({\mathbf{\theta }}_t{)}^{-1}+{\tau }_t^{-2}\mathbf{I}).$$

$$p({\mathbf{\Theta }}_{1:T},{\hat{\mathbf{y}}}_{1:T-1}({\mathbf{S}}_{2:T}^{\ast })|{\mathbf{Z}}_{1:T})$$

$$\propto p({\mathbf{\Theta }}_1)\tilde{L}({\mathbf{Z}}_1|{\mathbf{\Theta }}_1)\prod _{t=2}^{T}p({\mathbf{\Theta }}_t)\tilde{L}({\mathbf{Z}}_t|{\mathbf{\Theta }}_t,{\hat{\mathbf{y}}}_{t-1}({\mathbf{S}}_t))\tilde{p}({\hat{\mathbf{y}}}_{t-1}({\mathbf{S}}_t^{\ast })|\cdot ),$$

$${\mathbf{\beta }}_t|{\mathbf{Z}}_t,{\hat{\mathbf{y}}}_{t-1}({\mathbf{S}}_t),{\mathbf{\theta }}_t,{\mathbf{\gamma }}_t{\tau }_t^2\sim N({\mathbf{V}}_{{\beta }_t}^{\ast }{\mathbf{\mu }}_{{\beta }_t}^{\ast },{\mathbf{V}}_{{\beta }_t}^{\ast }),$$

$${\mathbf{\gamma }}_t|{\mathbf{Z}}_t,{\hat{\mathbf{y}}}_{t-1}({\mathbf{S}}_t),{\mathbf{\theta }}_t,{\mathbf{\beta }}_t,{\tau }_t^2\sim N({\mathbf{V}}_{{\gamma }_t}^{\ast }{\mathbf{\mu }}_{{\gamma }_t}^{\ast },{\mathbf{V}}_{{\gamma }_t}^{\ast }),$$

$$p({\mathbf{\beta }}_t,{\mathbf{\gamma }}_{t-1},{\sigma }_t^2$$

$$\times N({\mathbf{Z}}_t|{\zeta }_{t-1}({\mathbf{S}}_t)\circ {\hat{y}}_{t-1}({\mathbf{S}}_t)+{\mathbf{h}}_t^T{\mathbf{\beta }}_t,{\sigma }_t^2{\tilde{\mathbf{\Sigma }}}_t).$$

$${\mathbf{\gamma }}_{t-1}|{\mathbf{\beta }}_t,{\sigma }_t^2,{\mathbf{Z}}_t,{\hat{y}}_{t-1}({\mathbf{S}}_t))$$

$${\mathbf{\beta }}_t|{\sigma }_t^2,{\mathbf{Z}}_t,{\hat{y}}_{t-1}({\mathbf{S}}_t)$$

$${\sigma }_t^2|{\mathbf{Z}}_t,{\hat{y}}_{t-1}({\mathbf{S}}_t)$$

$${\hat{z}}_t({\mathbf{S}}_{t,k})=1_{t>1}(t){\mathbf{g}}_{t-1}^T({\mathbf{S}}_{t,k}){\hat{\mathbf{\gamma }}}_{t-1}{\hat{y}}_{t-1}({\mathbf{S}}_{t,k})+{\mathbf{h}}_t^T({\mathbf{S}}_{t,k}){\hat{\mathbf{\beta }}}_t+V_{t,{\mathbf{S}}_{t,k}}{\mu }_{t,{\mathbf{S}}_{t,k}}.$$

$$\text{RMSPE}=\frac{\sum_{k=1}^{K}\biggl{[}\sum_{\mathbf{s}=\mathbf{S}_{t,k}}(z_{t}(\mathbf{s})-\hat{z}_{t}(\mathbf{s}))^{2}/n_{k}\biggr{]}}{K}.$$

$${\bar{\hat{y}}}_t({\mathbf{S}}_t^{\ast })=1_{t>1}(t){\mathbf{g}}_{t-1}^T({\mathbf{S}}_t^{\ast }){\hat{\mathbf{\gamma }}}_{t-1}{\hat{y}}_{t-1}({\mathbf{S}}_t^{\ast })+{\mathbf{h}}_t^T({\mathbf{S}}_t^{\ast }){\hat{\mathbf{\beta }}}_t+V_{t,{\mathbf{S}}_t^{\ast }}{\mu }_{t,{\mathbf{S}}_t^{\ast }}$$

$${\hat{y}}_{t-1}(\mathbf{s})|{\mathbf{Z}}_{t-1},{\hat{y}}_{t-2}(\mathbf{s}))\sim N({\zeta }_{t-2}{\hat{y}}_{t-2}(\mathbf{s})+{\mathbf{h}}_{t-1}^T(\mathbf{s}){\mathbf{\beta }}_{t-1}+V_{t-1,\mathbf{s}}{\mu }_{t-1,\mathbf{s}},V_{t-1,\mathbf{s}})$$

$$\tilde{p}({\mathbf{w}}_t|\cdot )$$

$$=\text{exp}\left(-\frac{1}{2}{\mathbf{w}}_t^T{\mathbf{B}}_t^T{\mathbf{F}}_t^{-1}{\mathbf{B}}_t{\mathbf{w}}_t\right),$$

$${\mathbf{B}}_{t,s_{t,i},j}=\{\begin{array}{l}1,\text{if}i=j,\\ -{\mathbf{B}}_{t,s_{t,i}}[,k],\text{if}{s}_{t,j}\text{is the}{k}^{th}\text{element in}{N}_t(s_{t,i}),\\ 0,\text{Others}.\end{array}$$

$${\mu }_t({\mathbf{s}}_{t,k})=$$

$$+{\mathbf{h}}_t^T({\mathbf{s}}_{t,k}){\mathbf{\beta }}_t+w_t({\mathbf{s}}_{t,k}),$$

$${\Lambda }_{k,l}^{(t,t)}=\text{cov}(z_t({\mathbf{s}}_{t,k}),z_t({\mathbf{s}}_{t,l})|\cdot )=\sum _{i=1}^{t-1}{1}_{\{{\mathbf{s}}_{t,k},{\mathbf{s}}_{t,l}\notin {\mathbf{S}}_i\}}({\mathbf{s}}_{t,k},{\mathbf{s}}_{t,l})\left\{\prod _{j=i}^{t-1}{\zeta }_j({\mathbf{s}}_{t,k}{)}^T{\zeta }_j({\mathbf{s}}_{t,l})\right\}{C}_i({\mathbf{s}}_{t,k},{\mathbf{s}}_{t,l}|{\mathbf{\theta }}_i)$$

$$+1_{{\mathbf{s}}_{t,k}={\mathbf{s}}_{t,l}}({\mathbf{s}}_{t,k},{\mathbf{s}}_{t,l}){\tau }_t^2,$$

$${\Lambda }_{k,l}^{(t,t^{\prime })}=\text{cov}(z_t(s_{t,k}),z_{t^{\prime }}({\mathbf{s}}_{t^{\prime },l})|\cdot )=\sum _{i=1}^{\text{min}(t,t^{\prime })-1}{1}_{\{{\mathbf{s}}_{t,k},{\mathbf{s}}_{t^{\prime },l}\notin {\mathbf{S}}_i\}}({\mathbf{s}}_{t,k},{\mathbf{s}}_{t^{\prime },l})\left\{\prod _{j=i}^{\text{min}(t,t^{\prime })-1}{\zeta }_j({\mathbf{s}}_{t,k}{)}^T{\zeta }_j({\mathbf{s}}_{t^{\prime },l})\right\}$$

$$\times C_i({\mathbf{s}}_{t,k{\mathbf{s}}_{t^{\prime },l}}|{\mathbf{\theta }}_i)+1_{\{{\mathbf{s}}_{t,k},{\mathbf{s}}_{t^{\prime },l}\notin {\mathbf{S}}_{\text{min}(t,t^{\prime })}\}}({\mathbf{s}}_{t,k},{\mathbf{s}}_{t^{\prime },l})C_{\text{min}(t,t^{\prime })}({\mathbf{s}}_{t,k},{\mathbf{s}}_{t^{\prime },l}|{\mathbf{\theta }}_{\text{min}(t,t^{\prime })}),$$

$${\mathbf{V}}_{{\beta }_t}^{\ast }$$

$${\mathbf{\mu }}_{{\beta }_t}^{\ast }$$

Full text

Keywords: Recursive co-kriging; Nearest neighbor Gaussian process; Remote sensing

1 Introduction

Global geophysical information is measured daily by numerous satellite sensors. Due to aging and exposure to the harsh environment of space the satellite sensors degrade over time, resulting in decreased performance reliability. Decreased performance may affect data measurement accuracy(Goldberg, 2011). In addition, newer satellites with technologically more advanced sensors provide information of higher fidelity than older sensors. These discrepancies in sensor performance have created the need to develop efficient methods to analyse daily global remote sensing measurements with varying fidelity. Here, our work is motivated by a data set produced from the high-resolution infrared radiation sounder (HIRS), which provides hundred of thousands of measurements from multiple satellite platforms daily.

Multiple methods in remote sensing have been developed to assess satellite sensor performance and consistency (Chander et al., 2013; Xiong et al., 2010; National Research Council, 2004). These methods do not account for spatial correlation and oversimplify the relationship between sensors. Furthermore, statistical methods to analyse these data sets which account for spatial correlation pose challenges due to the multifideltiy presence as well as the size and computationally intensive procedures. Nguyen et al. (2012, 2017) have proposed data fusion techniques to model multivariate spatial data at potentially different spatial resolutions based on fixed ranked kriging (Cressie and Johannesson, 2008). The accuracy of this approach relies on the number of basis functions and can only capture large scale variation of the covariance function. When the data sets are dense, strongly correlated, and the noise effect is sufficiently small, the low rank kriging techniques have difficulty accounting for small scale variation (Stein, 2014).

Many statistical methods for large spatially correlated data sets have been developed over the past two decade. For instance, we distinguish, the low-rank approximation methods (Banerjee et al., 2008; Cressie and Johannesson, 2008), approximate likelihood methods (Stein et al., 2004; Gramacy and Apley, 2015), covariance tapering methods (Furrer et al., 2006; Kaufman et al., 2008; Du et al., 2009), sparse structures (Lindgren et al., 2011; Nychka et al., 2015; Datta et al., 2016; Ma and Kang, 2020; Michele Peruzzi and Finley, 2022), lower dimensional conditional distributions (Vecchia, 1988; Stein et al., 2004; Datta et al., 2016; Katzfuss and Guinness, 2021), and multiple-scale approximation (Sang and Huang, 2012; Katzfuss, 2017; Abdulah et al., 2023; Shirota et al., 2023). All these methods have been developed for data sets obtained from the same source or instrument which translates to a single fidelity data source. However, their extension to multi-fidelity data sets is not straightforward.

Autoregressive co-kriging models (Kennedy and O’Hagan, 2000; Qian et al., 2005; Le Gratiet, 2013), originally built for computer simulation problems, can be used for the analysis of multiple fidelity remote sensing observations with spatially nested structure and no random error. Nested design in the multifidelity setting means the design points at the higher fidelity levels are subsets of the lower fidelity ones. Konomi and Karagiannis (2021) and Ma et al. (2022) relaxed the nested design requirements by properly introducing an imputation mechanism. However, the aforesaid methods rely on Gaussian process models and are computationally impossible for big data problems. For cases when the observed space can be expressed as a tensor product Konomi et al. (2023) uses a separable covariance function within the co-kriging model to improve the computational efficiency. For large data sets which are irregularly positioned over space and contaminated with random error, Cheng et al. (2021) proposed the nearest neighbour co-kriging Gaussian process (NNCGP) to embed nearest neighbor Gaussian process (NNGP; Datta et al., 2016) into an autoregressive co-kriging model to make computations possible. NNCGP achieves this by using imputation ideas into the latent variables to construct a nested reference set of multiple NNGP levels. Although NNCGP makes the analysis of big multi-fidelity data sets computationally possible, its computational speed depends on an expensive iterative MCMC procedure which makes it impractical for analysing daily large data sets.

To overcome the iterative MCMC procedure, we propose a recursive formulation based on the latent variable of the NNCGP model following similar ideas with Le Gratiet and Garnier (2014), who proposed the recursive formulation directly into the observations. Based on this new formulation, which we call recursive nearest neighbors co-kriging (RNNC), we are able to build a nearest neighbors co-kriging model with $T$ levels by building $T$ conditionally independent NNGPs. This enables the development of two alternative inferential procedures which aim to reduce high-dimensional parametric space, improve convergence, and reduce computational time in comparison to the NNCGP. Both proposed procedures are able to address applications for large non-nested and irregular spatial data sets from different platforms and with varying quality. The first proposed procedure, called Collapsed RNNC, reduces the MCMC posterior sampling space by integrating out the spatial latent variables. Based on the collapsed RNNC, we propose an MCMC free procedure to speed up Bayesian inference. We build an algorithm which sequentially decomposes the parametric space into conditionally independent parts for each fidelity level where we can apply a K-fold optimization method. Each sequential step can be viewed as collapsed NNGP (Finley et al., 2019) where the bases functions of the Gaussian process mean are determined at the previous step. We name this second inferential procedure conjugate RNNC. We note that the MCMC free procedure proposed in (Finley et al., 2019) cannot be applied directly in the NNCGP model because the computational complexity of the $K$-fold cross-validation method depends on the dimension of the parametric space. Our simulation study and our analysis of the HERS data set shows that the proposed conjugate RNNC procedure reduces the computational time notably without significantly sacrificing prediction accuracy over the existing NNCGP approach.

The layout of the paper is as follows. In Section 2, we introduce the high-resolution infrared radiation sounder data studied in this work. In Section 3, we review the NNCGP model. In Section 4, we introduce the proposed RNNC model. In Section 4.1, we integrate out the latent variables from the model and design an MCMC algorithm for this model. In Section 4.2, we design an MCMC free approach tailored to the proposed RNNC model that facilitates parametric and predictive inference. In Section 5, we investigate the performance of the proposed procedure. Specifically in Section 5.1 we introduce two simulation studies and in Section 5.2 we implement the proposed method for the analysis of data sets from two satellites, NOAA-14 and NOAA-15. Finally,in Section 6 we give a summary and conclusion.

2 High-resolution Infrared Radiation Sounder Data

Satellite soundings have been providing measurements of the Earth’s atmosphere, oceans, land, and ice since the 1970s to support the study of global climate system dynamics. Long term observations from past and current environmental satellites are widely used in developing climate data records (CDR) (National Research Council, 2004). HIRS mission objectives include observations of atmospheric temperature, water vapor, specific humidity, sea surface temperature, cloud cover, and total column ozone. The HIRS instrument is comprised of twenty channels, including twelve longwave channels, seven shortwave channels, and one visible channel. The dataset being considered in this study is limb-corrected HIRS swath data as brightness temperatures (Jackson et al., 2003). The data is stored as daily files, where each daily file records approximately 120,000 geolocated observations. The current archive includes data from NOAA-5 through NOAA-17 along with Metop-02, covering the time period of 1978-2017. In all, this data archive is more than 2 TB, with an average daily file size of about 82 MB. The HIRS CRD faces some common challenges regarding the consistency and accuracy over time, due to degradation of sensors and intersatellite discrepancies. Furthermore, there is missing information caused by atmospheric conditions such as thick cloud cover.

We examine HIRS Channel 5 observations from a single day, March 1, 2001, as illustrated in Figure 1. On this day, we may exploit a period of temporal overlap in the NOAA POES series where two satellites captured measurements: NOAA-14 and NOAA-15. The HIRS sensors on these two satellites have similar technical designs which allow us to ignore the spectral and spatial footprint differences. NOAA-14 became operational in December 1994 while NOAA-15 became operational in October 1998. The spatial resolution footprint for both satellites is approximately 10 km at nadir. Given the sensor age difference, it is reasonable to consider that the instruments on-board NOAA-15 are in better condition than those of NOAA-14 and hence provide more accurate data. Therefore, we treat observations from NOAA-14 as a low fidelity dataset, and those from NOAA-15 as a high fidelity dataset.

3 Nearest Neighbor Co-kriging Gaussian Process

Let $y_{t}(\mathbf{s})$ denote the output function at the spatial location $s$ at fidelity level $t=1,...,T$ in a system with fidelity $T$ levels. The fidelity level index $t$ runs from the least accurate to the most accurate one. Let $z_{t}(\mathbf{s})$ denote the observed output at location $s$. We specify the co-kriging model as:

[TABLE]

where $z_{t}(\mathbf{s})$ is contaminated by additive random noise $\epsilon_{t}\sim N(0,\tau_{t}^{2})$ for $t=2,\ldots,T$, and $y_{1}(\mathbf{s})=\mathbf{h}^{T}_{1}(\mathbf{s})\bm{\beta_{1}}+w_{1}(\mathbf{s})$ is the noiseless output. Here, $\zeta_{t-1}(\mathbf{s})$ and $\delta_{t}(\mathbf{s})$ represent the scale and additive discrepancies between systems with fidelity levels $t$ and $t-1$, $\mathbf{h}_{t}(\cdot)$ is a vector of preselected bases functions, and $\bm{\beta}_{t}$ is a vector of coefficients at fidelity level $t$. The latent random function $w_{t}(\mathbf{s})$ is modeled as a Gaussian process, mutually independent for different $t$; i.e. $w_{t}(\cdot)\sim GP(0,C_{t}(\cdot,\cdot;\bm{\theta}_{t}))$ where $C_{t}(\cdot,\cdot;\bm{\theta}_{t})$ is a covariance function with covariance parameters $\bm{\theta}_{t}$ at fidelity level $t$. Any well defined covariance function can be used $C_{t}(\mathbf{s},\mathbf{s}^{\prime}|\bm{\theta}_{t})=\sigma_{t}^{2}R(\mathbf{s},\mathbf{s}^{\prime}|\bm{\phi}_{t})$, where $\bm{\theta}_{t}=\{\sigma_{t}^{2},\bm{\phi}_{t}\}$. This indicates that discrepancy term $\delta_{t}(\mathbf{s})$ is a Gaussian process. Finally, the unknown scale discrepancy function $\zeta_{t-1}(\mathbf{s})$ is modeled as a basis expansion $\zeta_{t-1}(\mathbf{s}|\bm{\gamma}_{t-1})=\mathbf{g}_{t-1}(\mathbf{s})^{T}\bm{\gamma}_{t-1}$ (usually low degree), where $\mathbf{g}_{t}(\mathbf{s})$ is a vector of polynomial basis functions and $\{\bm{\gamma}_{t-1}\}$ is a vector of random coefficients, for $t=2,\dots,T$.

Let us assume the system is observed at $n_{t}$ locations at fidelity level $t$. Let $\mathbf{S}_{t}=\{\mathbf{s}_{t,1},\dots,\mathbf{s}_{t,n_{t}}\}$ be the set of $n_{t}$ observed locations, let $\mathbf{w}_{t}={w}_{t}(\mathbf{S}_{t})=\{w_{t}(\mathbf{s}_{t,1}),\dots,w_{t}(\mathbf{s}_{t,n_{t}})\}$ the latent spatial random effect vector at fidelity level $t$, and let $\mathbf{Z}_{t}={z}_{t}(\mathbf{S}_{t})=\{z_{t}(\mathbf{s}_{t,1}),\dots,z_{t}(\mathbf{s}_{t,n_{t}})\}$ represent the observed output at fidelity level $t$. If data $\{\mathbf{Z}_{t}\}$ are observed in non-nested locations across the fidelity levels, the calculation of the likelihood requires $\mathcal{O}((\sum_{t=1}^{T}n_{t})^{3})$ flops to invert the covariance matrix of the observations (denoted by $\bm{\Lambda}$) and additional $\mathcal{O}((\sum_{t=1}^{T}n_{t})^{2})$ memory to store it as explained in (Konomi and Karagiannis, 2021). To reduce the computational complexity, Cheng et al. (2021) proposed NNCGP which assigns conditionally independent NNGP models within a nested reference set. For Bayesian inference, Cheng et al. (2021) proposed a Gibbs sampler taking advantage of the nearest neighbour structure at each fidelity level. However, this sampler is based on updating a conditionally independent high dimensional latent variable which could cause slow convergence and high autocorrelation (Liu et al., 1994). The slow convergence can significantly increase the number of the Gibbs sampler iterations $I$. The overall computational cost of NNCGP for $m$ neighbours and non-nested spatial locations is $\mathcal{O}(I\times(\sum_{t=1}^{T}n_{t})m^{3})$ floating point operations (flops). Also, for a fixed computational budget, the produced Monte Carlo estimates may be sensitive to the initial values of the MCMC sampler. Despite reducing the computational complexity to linear for every MCMC iteration, the number of the iterations can significantly increase computational cost. This simple observation makes the existing NNCGP too expensive for the vast majority of real remote sensing applications.

4 Recursive Nearest Neighbor Co-kriging Model

Improvement in the convergence of MCMC can be achieved by integrating out the latent variable $\mathbf{w}=(\mathbf{w}_{1},\dots,\mathbf{w}_{S})$ from the Bayesian hierarchical NNCGP model which allows dimension reduction in the sampling space and the involved posterior distributions. However, integrating out the latent variables $\mathbf{w}$ in the NNCGP model is not feasible under non-nested designs. This is because the posterior distribution of latent variable of the lower fidelity is affected by the likelihood of higher fidelity. To make possible the integration of latent variables $\mathbf{w}$, we propose a recursive formulation for the co-kriging model by using ideas similar to (Le Gratiet and Garnier, 2014). Precisely, our proposed recursive nearest neighbors co-kriging (RNNC) has the following hierarchical structure:

[TABLE]

where $\delta_{t}(\mathbf{s})$ is a Gaussian process as before and $\hat{y}_{t-1}(\mathbf{s})$ is a Gaussian process with distribution $[\hat{y}_{t-1}(\mathbf{s})|\mathbf{Z}_{t-1},\hat{y}_{t-2}(\mathbf{s}),\bm{\theta}_{t-1},\bm{\beta}_{t-1}]$. Essentially, we express $y_{t}(\mathbf{s})$ (the Gaussian process response at level $t$) as a function of the Gaussian process $y_{t-1}(\mathbf{s})$ conditioned by the values $\mathbf{Z}^{(t-1)}=(\mathbf{Z}_{1},\dots,\mathbf{Z}_{t-1})$. For computational efficiency, we assume NNGP independent priors for $w_{t}(\mathbf{s})$, $t=1,\dots,T$. Based on the NNGP priors, the conditional distribution can be computed for all types of reference sets. So, based on the recursive representation we can relax the nested condition on the NNCGP nested reference set. Specifically,

[TABLE]

with, \mu_{t-1,\mathbf{s}}=V_{t-1,\mathbf{s}}^{-1}\mathbf{B}_{t-1,\mathbf{s}}\bigl{[}z_{t-1}(N_{t-1}(\mathbf{s}))-\mathbf{h}^{T}_{t-1}(N_{t-1}(\mathbf{s}))\bm{\beta}_{t-1}-\zeta_{t-2}(N_{t-1}(\mathbf{s}))\circ\hat{y}_{t-2}(N_{t-1}(\mathbf{s}))\bigr{]}, $\mathbf{B}_{t,\mathbf{s}}=C_{\mathbf{s},N_{t}(\mathbf{s})}^{T}C_{N_{t}{\mathbf{s}}}^{-1}$, and $V_{t,\mathbf{s}}=C(\mathbf{s},\mathbf{s})-C_{\mathbf{s},N_{t}(\mathbf{s})}^{T}C_{N_{t}{\mathbf{s}}}^{-1}C_{\mathbf{s},N_{t}(\mathbf{s})}$. The $\circ$ represents the Hadamard product between two matrices.

Using the Markovian property of the co-kriging model (O’Hagan, 1998), the joint likelihood of the proposed model in (4.1) can be factorized as a product of likelihoods at different fidelity levels conditional on $\hat{y}_{t-1}(\mathbf{S}_{t})=\{\hat{y}_{t-1}(\mathbf{s}_{t,1}),\dots,\hat{y}_{t-1}(\mathbf{s}_{t,n_{t}})\}$ for $t=2,\dots,T$ and prior $\mathbf{w}_{t}$ for $t=1,...,T$, i.e.:

[TABLE]

where $(\cdot)$ denotes all the parameters associated with the model. This representation makes it possible to integrate out the latent variable $\mathbf{w}_{t}$ independently for each fidelity level $t=1,\dots,T$ .

4.1 Collapsed Recursive Nearest Neighbor Co-kriging Model

We represent the multivariate Gaussian latent variable $\mathbf{w}_{t}(\mathbf{S}_{t})$ as a linear model:

[TABLE]

for $t=1,\dots,T$. We set $\eta_{t,i}\sim N(0,d_{t,i,i})$ independently for all $t,i$, $d_{t,1,1}=\variation(w_{t,1})$ and $d_{t,i,i}=\variation(w_{t,i}|\{w_{t,j};j<i\})$ for $i=2,\dots,n_{t}$ and $t=1,\dots,T$. In a matrix form we can write $\mathbf{w}_{t}(\mathbf{S}_{t})=\mathbf{A}_{t}\mathbf{w}_{t}(\mathbf{S}_{t})+\bm{\eta}_{t}$, where $\mathbf{A}_{t}$ is an $n\times n$ strictly lower-triangular matrix and $\bm{\eta}_{t}\sim N(0,\mathbf{D})$ and $\mathbf{D}$ is diagonal. Based on the structure of $\mathbf{A}_{t}$, we can write the covariance of each level as $\mathbf{C}_{t}(\bm{\theta}_{t})=(\mathbf{I}_{t}-\mathbf{A}_{t})^{-1}\mathbf{D}_{t}(\mathbf{I}_{t}-\mathbf{A}_{t})^{-T}$. The NNGP prior constructs a sparse strictly lower triangular matrix $\mathbf{A}$ with no more than $m\ (\text{where }m\ll n)$ non-zero entries in each row resulting in an approximation of the covariance matrix $\mathbf{C}_{t}$. So the approximated inverse $\tilde{\mathbf{C}}_{t}^{-1}(\bm{\theta}_{t})=(\mathbf{I}_{t}-\mathbf{A}_{t})\mathbf{D}_{t}^{-1}(\mathbf{I}_{t}-\mathbf{A}_{t})^{T}$ is a sparse matrix and can be computed based on $\mathcal{O}(n_{t}m^{3})$ operations.

We call the integrated version of the above model collapsed RNNC model. Specifically, after integrating out ${\bf w}_{t}$ the proposed RNNC model can be written as:

[TABLE]

for $t=2,\ldots,T$, where $\tilde{\Lambda}_{t}(\bm{\theta}_{t},\tau_{t})=\tilde{C}_{t}(\bm{\theta}_{t})+\tau_{t}^{2}\mathbf{I}=\sigma_{t}^{2}\tilde{\mathbf{R}}_{t}(\bm{\phi}_{t})+\tau_{t}^{2}\mathbf{I}$ is the covariance matrix of the observations, $\tilde{C}_{t}(\bm{\theta}_{t})$ is the sparse covariance matrix with parameters $\bm{\theta}_{t}=\{\sigma_{t}^{2},\bm{\phi}_{t}\}$ and $\tau_{t}^{2}$ is the variance of the error $\epsilon_{t}$ at level $t$. By applying Sherman-Morrison-Woodbury formula, the inverse and determinant of $\tilde{\bm{\Lambda}}$ get the computationally convenient form

[TABLE]

For simplicity let us denote $\bm{\Theta}_{t}=(\bm{\beta}_{t},\gamma_{t},\bm{\theta}_{t},\tau_{t})$. Based on this representation, the joint posterior approximation of all the unknowns is:

[TABLE]

where $\tilde{L}(\mathbf{Z}_{t}|\bm{\Theta}_{t},\hat{\mathbf{y}}_{t-1}(\mathbf{S}_{t}))$ is the approximated likelihood using the sparse representation and $\tilde{p}(\hat{\mathbf{y}}_{t-1}(\mathbf{S}_{t}^{*})|\cdot)$ the nearest neighbor Gaussian process prediction at locations $\mathbf{S}_{t}^{*}=\bigcup\limits_{i=t+1}^{T}\mathbf{S}_{i}\backslash\mathbf{S}_{t}=\{s_{t,1}^{*},\ldots,s_{t,n_{t}^{*}}^{*}\}$ as a set of knots of fidelity level $t$. This contains the observed locations that are not in the $t^{th}$ level but in the higher fidelity levels. The $(\cdot)$ represents the parameters and data necessary to produce the prediction distortion at level $t-1$. Note that the prediction probability can be excluded for cases with hierarchically nested structure for the spatial locations. For each level, a Gibbs sampler can be employed to fascilitate inference based on the conditional representation $p(\bm{\Theta}_{T}|\hat{\mathbf{y}}_{t-1}(\mathbf{S}_{t}),\mathbf{Z}_{T})$ and $p(\hat{\mathbf{y}}_{t-1}(\mathbf{S}_{t}^{*})|\bm{\Theta}_{T},\mathbf{Z}_{T})$ which is given in Eq. (4.2).

For $\bm{\Theta}_{t}=(\bm{\beta}_{t},\gamma_{t},\bm{\theta}_{t},\tau_{t})$, by assigning independent conjugate prior $\bm{\beta}_{t}\sim N(\bm{\mu}_{\bm{\beta}_{t}},\mathbf{V}_{\bm{\beta}_{t}})$ and $\gamma_{t}\sim N(\bm{\mu}_{\gamma_{t}},\mathbf{V}_{\gamma_{t}})$, we achieve explicit forms of the conditional distribution for parameters $\bm{\beta}_{t}$ and $\gamma_{t}$ as:

[TABLE]

where $\mathbf{V}_{\beta_{t}}^{*}\bm{\mu}_{\beta_{t}}^{*},\mathbf{V}_{\gamma_{t}}^{*},\bm{\mu}_{\gamma_{t}}^{*}$ are given in Appendix C. Finally, for each fidelity level $t$, we use a Metropolis-Hastings (MH) algorithm targeting the distribution $p(\bm{\theta}_{t},\tau_{t}^{2}|\mathbf{Z}_{t},\hat{\mathbf{y}}_{t-1}(\mathbf{S}_{t}),\bm{\beta}_{t})$ to carry out the inference.

In the special case of hierarchical nested structure for the spatial locations, we can avoid sampling from $p(\hat{\mathbf{y}}_{t-1}(\mathbf{S}_{t}^{*})|\bm{\Theta}_{T},\mathbf{Z}_{T})$ using the observed locations $z_{t-1}(\mathbf{S}_{t}^{*})$. We can prove that the mean and variance of the predictive distribution at level $T$ of the collapsed RNNC is the same as the mean and variance of the predictive distribution of the NNCGP. The proof is very similar to Le Gratiet and Garnier (2014) in the sense that we just need to substitute the GP priors with the NNGP priors and add the nugget effect in each level.

4.2 Conjugate Recursive Nearest Neighbor Co-kriging Model

Both NNCGP and collapsed RNNC models rely on the MCMC inference which can be practically prohibitive when analyzing thousands or millions of spatial data sets. Following recent work by Finley et al. (2019), we propose a MCMC free procedure to achieve exact Bayesian inference at a more practical time. Because the computational efficiency of the estimation procedure in Finley et al. (2019) is sensitive to the number of parameters, it cannot be applied directly to our model. We utilize the RNNC model conditionally independent posterior representation to decompose the parametric space into smaller different groups based on the fidelity levels. To make MCMC free inference possible, we re-parameterize the covariance function of the collapsed recursive co-kriging model as $\tilde{\bm{\Lambda}}_{t}(\bm{\theta}_{t},\tilde{\tau}_{t}^{2})=\sigma_{t}^{2}\tilde{\bm{\Sigma}}_{t}$, where $\tilde{\bm{\Sigma}}_{t}=\tilde{\mathbf{R}}_{t}+\tilde{\tau}_{t}^{2}\mathbf{I}$ , $\tilde{\mathbf{R}}_{t}$ is the nearest-neighbor approximation correlation matrix, and $\tilde{\tau}_{t}^{2}=\frac{\tau_{t}^{2}}{\sigma_{t}^{2}}$. To avoid the computational bottleneck due to the MCMC, we propose to make fast estimation $(\bm{\phi}_{t},\tilde{\tau}_{t}^{2})$ through a cross-validation approach for each level as well as use the prediction means of $\mathbf{y}_{t-1}(\mathbf{S}_{t})$ based on the estimated values. We estimate $\hat{y}_{t}(\mathbf{S}_{t}^{*})$ by the posterior mean $\bar{\hat{y}}_{t}(\mathbf{S}_{t}^{*})=\mathbf{1}_{t>1}(t)\mathbf{g}_{t-1}^{T}(\mathbf{S}_{t}^{*})\hat{\bm{\gamma}}_{t-1}\hat{y}_{t-1}(\mathbf{S}_{t}^{*})+\mathbf{h}_{t}^{T}(\mathbf{S}_{t}^{*})\hat{\bm{\beta}}_{t}+V_{t,\mathbf{S}_{t}^{*}}\mu_{t,\mathbf{S}_{t}^{*}}$. In the case that we have nested locations, $y_{t}(\mathbf{s}_{u})$ for a location $\mathbf{s}_{u}\in\mathbf{S}_{t-1}$ is estimated with an empirical approach $\bar{\hat{y}}_{t-1}(\mathbf{s}_{u})$ as $z_{t-1}(\mathbf{s}_{u})$ and its variance is equal to the variance of the nugget effect. Given $\bm{\phi}_{t},\tilde{\tau}_{t}^{2}$ and $\hat{y}_{t}(\mathbf{S}_{t})$, the convariance matrix $\tilde{\bm{\Sigma}}_{t}$ can be calculated analytically.

For computational convenience, we assign an independent conjugate prior for the parameters of each level such as $p(\bm{\beta}_{1},\dots,\bm{\beta}_{T},\sigma_{1}^{2},\dots,\sigma_{T}^{2},\bm{\gamma}_{1},\dots,\bm{\gamma}_{T})=\prod_{t=1}^{T}p(\bm{\beta}_{t})p(\sigma_{t}^{2})p(\bm{\gamma}_{t})$ such as $\bm{\beta}_{t}\sim N(\bm{\mu}_{\bm{\beta}_{t}},\sigma_{t}^{2}\mathbf{V}_{\bm{\beta}_{t}})$, $\sigma_{t}^{2}\sim IG(a_{t},b_{t})$, and $\bm{\gamma}_{t}\sim N(\bm{\mu}_{\bm{\gamma}_{t}},\sigma_{t+1}^{2}\mathbf{V}_{\bm{\gamma}_{t}})$. Based on this specifications, the posterior density function can be separated for each level $t$ such as:

[TABLE]

We can compute the full conditional density function of $\bm{\gamma}_{t-1},\bm{\beta}_{t}$, and $\sigma_{t}^{2}$ as

[TABLE]

were $\tilde{\mathbf{V}}_{\bm{\gamma}_{t-1}},\tilde{\bm{\mu}}_{\bm{\gamma}_{t-1}},\tilde{\mathbf{V}}_{\beta_{t}}\tilde{\bm{\mu}}_{\beta_{t}},\sigma_{t}^{2}\tilde{\mathbf{V}}_{\beta_{t}},a_{t}^{*}$, and $b_{t}^{*}$ are given analytically in Appendix D. Note that for $t=1$, $\bm{\gamma}_{0}$ and $y_{0}(\mathbf{S}_{t})$ do not exist. Also the conditional posterior density function of $\bm{\beta}_{1}$ and $\sigma_{1}^{2}$ are slightly different as explained in Appendix D.

A $K$-fold cross-validation method is used for the selection of optimal values for the parameters $\bm{\phi}_{t}$ and $\tilde{\tau}_{t}^{2}$ at level $t$ that provide best prediction performance for the model, from a group of candidates. The criteria for choosing $\bm{\phi}_{t}$ and $\tilde{\tau}_{t}^{2}$ can be the root mean square prediction error (RMSPE) over the $K$ folds of data set. The geolocated observations of the $t$ fidelity are partitioned into $K$ equal size subsets. Then, one of the subsets is used as a test set and the others are used for training. The procedure is repeated $K$ times such that each subset is used once as a test set. The computational complexity of these procedures is reduced significantly from the use of the NNGP priors in the recursive co-kriging model. The estimation, tuning and prediction procedure of conjugate RNNC model are given in Algorithm 1. Similar to NNCGP model, the conjugate RNNC model analyzes the data set of each fidelity level sequentially from the lowest level to the highest. For each single fidelity level $t$, the conjugate RNNC model is able to run in parallel for tuning the parameter $\bm{\phi}_{t}$ and $\tilde{\tau}_{t}^{2}$ using a $K$ fold cross validation procedure. Step $3$, for given values of $(\bm{\phi}_{t},\tilde{\tau}_{t}^{2})$, to estimate $(\hat{\bm{\beta}}_{t},\hat{\sigma}_{t}^{2})$ requires $\mathcal{O}(n_{t}m^{3}+n_{t}mp_{t}^{2})$ floating point operations (flops) where $p_{t}$ is the dimension of $\bm{\beta}_{t}$. Step 6, to predict at new locations $\mathbf{S}_{t}^{*}$ requires $\mathcal{O}(n_{t}^{*}m^{3})$ flops where $n_{t}^{*}$ is the dimension of $\mathbf{S}_{t}^{*}$. Step 7, to predict at a new location requires $\mathcal{O}(m^{3})$ flops. When we use parallel computing within a fidelity level for parameters $(\bm{\phi}_{t},\tilde{\tau}_{t}^{2})$ the computation in step $3$ becomes extremely fast. The conjugate RNNC model provides an empirical estimation of spatial effect parameter $\bm{\phi}_{t}$ and noise parameter $\tilde{\tau_{t}}^{2}$ within a given resolution. We note that the proposed MCMC free inference can be viewed as a sequential optimization technique which splits the parametric space into several lower dimension components where we can apply conditional independent conjugate NNGP models.

5 Synthetic Data Example and Real Data Analysis

We study the performance of our proposed procedures, the conjugate RNNC and the collapsed RNNC, as well as compare their performances with that of the sequential NNCGP. The empirical study is based on one synthetic data set example with nested and one with non-nested input data sets. Also we use a real satellite data set application. As measures of performance we use the root mean squared prediction errors (RMSPE), coverage probability of the $95\%$ equal tail credible interval (CVG($95\%$)), average length of the $95\%$ equal tail credible interval (ALCI($95\%$)), and continuous rank probability score (CRPS) (Gneiting and Raftery, 2007). Details of these measurements are given in Appendix E. The simulations were performed in MATLAB R2018a, on a personal computer with specifications (intelR i7-3770 3.4GHz Processor, RAM 8.00GB, MS Windows 64bit).We have also included a simulation study with four levels of fidelity in the supplementary materials.

5.1 Simulation Study

We consider a two-fidelity level system represented by the hierarchical statistical model (3.1) defined on a two dimensional unit square domain with univariate observation data sets for both $\mathbf{Z}_{1}$ and $\mathbf{Z}_{2}$. Let the design matrix be $\mathbf{h}(\mathbf{S}_{t})=\mathbbm{1}$, the autoregressive coefficient function be an unknown constant $\zeta_{1}(\mathbf{s})=\gamma_{1}$, and exponential covariance functions. We generate two synthetic data sets for the above statistical model. The true values of the parameters are listed in Table 1 and Table 3. The data sets on the nested spatial locations consists of observations $\mathbf{Z}_{1}$ and $\mathbf{Z}_{2}$ from $100\times 100$ grids $\mathbf{S}_{1}$ and $\mathbf{S}_{2}$, respectively. The data sets, shown in Figures 2a and 2b are based on a fully non-nested input where the low fidelity observations $\mathbf{Z}_{1}$ and the high fidelity observations $\mathbf{Z}_{2}$ are generated at irregularly located at point in sets $\mathbf{S}_{1}$ and $\mathbf{S}_{2}$ of size $5000$, while $\mathbf{S}_{1}\cap\mathbf{S}_{2}=\emptyset$. In all data sets, a few small square regions from $\mathbf{Z}_{2}$ are treated as a testing data-set, and the rest of $\mathbf{Z}_{2}$ and $\mathbf{Z}_{1}$ are treated as training data sets. The testing regions for the non-nested input can be seen as white boxes in Figure 2(b).

Regarding the Bayesian inference, we compared the sequential NNCGP model, with the proposed collapsed RNNC model and that with the conjugate RNNC model, on both nested and non-nested data sets. We assigned similar non-informative priors for all the four models. We assign independent conjugate prior on parameters $\beta_{1}\sim N(0,1000)$, $\beta_{2}\sim N(0,1000)$, and scale parameter $\gamma_{1}$. We assign independent inverse gamma prior on spatial variance parameters $\sigma_{1}^{2}\sim IG(2,1)$, and $\sigma_{2}^{2}\sim IG(2,1)$ and on the noise parameters $\tau_{1}^{2}\sim IG(2,1)$, $\tau_{2}^{2}\sim IG(2,1)$. We also assign uniform prior on the range parameters $\phi_{1}\sim U(0,25)$, and $\phi_{2}\sim U(0,25)$. For the collapsed RNNC model, we run Markov chain Monte Carlo (MCMC) samplers for $35000$ iterations where the first $5000$ iterations are discarded as a burn-in, and convergence of the MCMC sampler was diagnosed from the individual trace plots. The RMSPE with a 5-fold cross-validation was used for the conjugate RNNC model. We select $(\phi_{t},\tilde{\tau}_{t}^{2})$ on a grid, for $\phi_{t}$ from the range $[0.1,25]$, and for $\tilde{\tau}_{t}^{2}$ from the range $[0.0005,0.4]$. No significant differences were observed when we used 3-fold cross-validation and 7-fold cross-validation approach.

In Tables 1 and 3, we report the Monte Carlo estimates of the posterior means and the associated $95\%$ marginal credible intervals of the unknown parameters using the two different NNCGP based procedures: sequential NNCGP, collapsed RNNC, along with the posterior mean and tuned values of parameters using conjugate RNNC, with $m=10$. There is no significant difference in the estimation of parameters for all MCMC based models (NNCGP and collapsed RNNC) and the true values of the parameters are successfully included in the $95\%$ marginal credible intervals. The introduction of latent interpolants may have caused a small overestimation of $\tau_{2}^{2}$ for all models. Instead, the conjugate RNNC is underestimating the variance of the nugget for the second fidelity level. The parameter estimations can be improved with a semi-nested or nested structure between the observed locations of the fidelity levels, and it is also shown for the auto-regressive co-kriging model in Konomi and Karagiannis (2021). The conjugate RNNC model has similar performance on estimating the mean of the parameters compared to the NNCGP and RNNC models. However, it does not provide uncertainties regarding these estimations.

In Table 2 and Table 4, we report standard performance measures (defined in Appendix E) for the sequential NNCGP, collapsed RNNC and conjugate RNNC with $m=10$ number of neighbours. All performance measures indicate that the collapsed RNNC model has similar predictive ability with the sequential NNCGP model. The conjugate RNNC model produced RMSPE value that is 10% larger than other NNCGP and collapsed RNNC models, but it is still significantly smaller than the RMSPE values from single level NNGP and combined NNGP. The tables also show that the running time of the collapsed RNNC model is not different from the sequential NNCGP model, this is consistent with our previous discussion since the two procedures have the same computational complexity. We observe that the conjugate RNNC model has extremely smaller running time compared to sequential NNCGP models, since the inference for conjugate RNNC model requires the same amount of running time as one iteration in collapsed RNNC model. It is worth pointing out that the cross validation process and tuning process in Algorithm 1 are independent from each other, which makes conjugate RNNC model benefit from parallel computation environments and greatly reduce computational time.

Figure 2 and 3 provide the non-nested synthetic observations and the prediction plots from sequential NNCGP, collapsed RNNC, conjugate RNNC, combined NNGP and single level NNGP models. We observed that for the testing regions the NNCGP models provides similar prediction surfaces and all NNCGP models has the better presentation of patterns in prediction surface comparing to single level NNGP model.

5.2 Application to High-resolution Infrared Radiation Sounder data

We model our data based on the two-fidelity level conjugate RNNC model and on the two-fidelity level sequential NNCGP model. Moreover, we provide comparisons with the single level NNGP model and combined NNGP model. We consider a linear model for the mean of the Gaussian processes, in $y_{1}(\cdot)$ and $\delta_{2}(\cdot)$, with linear basis function representation $\mathbf{h}(s_{t})$ and coefficients $\bm{\beta}_{t}=\{\beta_{0,t},\beta_{1,t},\beta_{2,t}\}^{T}$. We consider the scalar discrepancy $\zeta(\mathbf{s})$ to be unknown constant and equal to $\gamma$. The number of nearest neighbors $m$ is set to 10, and the spatial process $\mathbf{w}_{t}$ is considered to have a diagonal anisotropic exponential covariance function.

We assign independent normal distribution priors with zero mean and large variances for $\beta_{0,t},\beta_{1,t},\beta_{2,t}$ and $\gamma$. We assign independent uniform prior distributions $U(0,1000)$ to the range correlation parameters $({\phi}_{t,1},{\phi}_{t,2})$ for $t=1,2$. Also, we assign independent $IG(2,1)$ prior distributions for the variance parameters $\sigma_{t}^{2}$ and $\tau_{t}^{2}$. For the Bayesian inference of the sequential NNCGP, we run the MCMC sampler with of $35,000$ iterations where the first $5,000$ iterations are discarded as a burn-in. For the Bayesian inference of conjugate RNNC, we consider using posterior means as the estimated values for parameters $\beta_{0,t},\beta_{1,t},\beta_{2,t}$, $\sigma_{t}^{2}$ and $\gamma$; we also use posterior means as the imputation values for latent process $\tilde{\mathbf{y}}_{t}$ and for the prediction values of $z(s_{p})$ at location $s_{p}\not\in\tilde{\mathbf{S}}_{t}$.

The prediction performance metrics of the four different methods are given in Table 5. Compared to the single level NNGP model and combined NNGP model, the sequential NNCGP model and conjugate RNNC model produced a 20-30% smaller RMSPE and their NSME is closer to 1. The sequential NNCGP model and collapsed RNNC model also produced larger CVG and smaller ALCI than the single level NNGP model and combined NNGP model. The result suggests that the NNCGP and RNNC models have a substantial improvement in terms of predictive accuracy in real data analysis too. In the prediction plots (Figure 4) of the testing data of NOAA-15, we observe that RNNC models are more capable of capturing the pattern of the testing data than single level NNGP model and combined NNGP model. This is reasonable because the observations from NOAA-14 have provided information of the testing region, and comparing to combined NNGP model, the NNCGP and RNNC models are capable of modeling the discrepancy of observations from different satellites. In the non-nested structure, the computational complexity of the single level NNGP model is $\mathcal{O}(n_{2}m^{3})$ and that of NNCGP model is $\mathcal{O}((n_{1}+n_{2})m^{3})$, for an MCMC iteration. However, the whole computational complexity of the conjugate RNNC model is $\mathcal{O}((n_{1}+n_{2})m^{3})$ with parallel computational environment, which makes it remarkably computationally efficient without losing significant prediction accuracy. This is consistent with the running times of the models shown in Table 5.

We apply the MCMC free conjugate RNNC model for gap-filling predictions based upon a discrete global grid. We chose to use $1^{\circ}$ latitude by $1.25^{\circ}$ longitude ($1^{\circ}\times 1.25^{\circ}$) pixels as grids with global spatial coverage from $-70^{\circ}$ to $70^{\circ}$N. By applying the NNCGP model, we predict gridded NOAA-15 brightness temperature data on the center of the grids, based on the NOAA-14 and NOAA-15 swath-based spatial support. The prediction plot (Figure 5) illustrates the ability of the MCMC free conjugate RNNC model to handle large irregularly spaced data sets and produce a gap-filled composite gridded dataset. The resulting global image of the brightness temperature is practically the same as the sequential NNCGP.

6 Summary and conclusions

We have proposed a new computationally efficient co-kriging method, the recursive nearest neighbor Autoregressive Co-Kriging (RNNC) model, for the analysis of large and multi-fidelity spatial data sets. In particular, we proposed two computationally efficient inferential procedures: a) the collapsed RNNC, and b) the conjugate RNNC. Regarding the collapsed RNNC, we integrate out the latent variables of the RNNC model which enables the factorisation of the likelihood into terms involving smaller and sparse covariance matrices within each level. Then, a prediction focused approximation is applied to the aforesaid model to further speed up the computation. The cross-validation using grid search on a two or three dimensional space is a computationally feasible method to estimate the hyperparameters. Regarding the proposed conjugate RNNC, it is MCMC free and at most computationally linear in the total number of all spatial locations of all fidelity levels. We compared the proposed collapsed RNNC and conjugate RNNC with NNCGP in a simulation study and a real data application of intersatellite calibration. We observed that similar to NNCGP, the collapsed and conjugate RNNC were also able to improve the accuracy of the prediction for the HIRS brightness temperatures from the NOAA-15 polar-orbiting satellite by incorporating information from an older version of the same HIRS sensor on board the polar orbiting satellite NOAA-14. The RNNC can be viewed as a modularization approach to NNCGP model in Bayesian statistics Bayarri et al. (2009) where the analysis is done in steps rather than jointly.

The proposed procedures can be used for a variety of large multi-fidelity data sets in remote sensing with overlapping areas of observed locations. A natural extension of our model can be done based on a recently proposed sparse plus Low-rank Gaussian Process (SPLGP) (Shirota et al., 2023) who used a combination of Gaussian predictive process and NNGP in an MCMC-free framework. A natural choice for introducing the non-stationarity in the conjugate RNNC is to use non-dynamic partition methods such as (Matthew J. Heaton and Terres, 2017; Konomi et al., 2019). Moreover, we can use more complex Vechias approximations (Vecchia, 1988; Stein et al., 2004; Guinness, 2018; Katzfuss et al., 2020), similar to the NNGP, where the ordering of the data is more complicated but results in a better approximation. These Vechias approximation techniques of ordering can be applied naturally in the proposed RNNC model, however, they are out of the scope of this paper and will be investigated in future work. Next steps will include extending the proposed method in the multivariate setting by using ideas from parallel partial autoregressive co-kriging (Ma et al., 2022) and NNGP spatial factor models (Taylor-Rodriguez et al., 2018). Spherical covariance function can be used for global data set analysis (Guinness and Fuentes, 2016), however their extension to anisotropic representation is not straightforward. Still, work needs to be done in developing new strategies for tuning the hyperparameters in more complex covariance functions with multiple parameters within a fidelity level.

Acknowledgements

The research of Konomi and Kang was supported in part by National Science Foundation grant NSF DMS-2053668 and the Taft Research Center at the University of Cincinnati. Kang was also supported in part by Simons Foundation’s Collaboration Award (#317298 and #712755).

Appendix

Appendix A NNGP specifications

The posterior distribution of

[TABLE]

where $\mathbf{F}_{t}=\text{diag}(F_{t,\mathbf{s}_{t,1}},F_{t,\mathbf{s}_{t,2}},\ldots,F_{t,\mathbf{s}_{t,n_{t}}})$, \mathbf{B}_{t}=\Big{(}\mathbf{B}_{t,1}^{T},\mathbf{B}_{t,2}^{T},\ldots,\mathbf{B}_{t,n_{t}}^{T}\Big{)}^{T}, and for each element in $\mathbf{B}_{t}$, we have \mathbf{B}_{t,i}=\Big{(}\mathbf{B}_{t,s_{t,i},1}^{T},\mathbf{B}_{t,s_{t,i},2}^{T},\ldots,\mathbf{B}_{t,s_{t,i},n_{t}}^{T}\Big{)}^{T} and

[TABLE]

Appendix B Mean and Variance Specifications

The mean vector $\bm{\mu}=\left(\mu_{1}(\mathbf{s}_{1,1}),\ldots,\mu_{1}(\mathbf{s}_{1,n_{1}}),\ldots,\mu_{T}(\mathbf{s}_{T,n_{T}})\right)$ is

[TABLE]

for $t=1,\ldots T$, $i=1,\ldots,n_{t}$. $\mathbf{1}_{\{\cdot\}}(\cdot)$ is the indicator function, and covariance matrix $\bm{\Lambda}$ is a block matrix with blocks $\Lambda^{(1,1)},\ldots,\Lambda^{(1,T)},\ldots,\Lambda^{(T,T)}$, and the size of $\bm{\Lambda}$ is $\sum_{t=1}^{T}n_{t}\times\sum_{t=1}^{T}n_{t}$. The $\Lambda^{(t,t)}$ components are calculated as:

[TABLE]

for $t\ \text{and}\ t^{\prime}=1,\ldots,T$; $k=1,\ldots,n_{t}$; $l=1,\ldots,n_{t^{\prime}}$, and

[TABLE]

for $t\neq t^{\prime}$, $\Lambda^{(t,t^{\prime})}$.

Appendix C Gibbs Sampler

[TABLE]

Appendix D Conjugate Conditional Probabilities

we derive the posterior distribution as

[TABLE]

The full conditional density function of $\bm{\gamma}_{t-1}$ is

[TABLE]

After integrate $\bm{\gamma}_{t-1}$ out, the conditional posterior density function of $\bm{\beta}_{t}$ is

[TABLE]

The marginalized posterior density function of $\sigma_{t}$ is

[TABLE]

the conditional posterior density function of $\bm{\beta}_{1}$ is

[TABLE]

and the marginal posterior density function of $\sigma_{1}^{2}$ is

[TABLE]

Appendix E Performance Metrics

In the empirical comparisons, we used the following performance metrics:

Root mean square prediction error (RMSPE) is defined as

[TABLE]

where $y^{\text{obs}}$ is the observed value in test data-set and $y_{i}^{\text{pred}}$ is the predicted value from the model. It measures the accuracy of the prediction from model. Smaller values of RMSPE indicate more a accurate model. 2. 2.

Nash-Sutcliffe model efficiency coefficient (NSME) is defined as:

[TABLE]

where $y^{\text{obs}}$ is the observed value in test data-set and $y_{i}^{\text{pred}}$ is the predicted value from the model. NSME gives the relative magnitude of the residual variance from data and the model variance. NSME values closer to $1$ indicate that the model has a better predictive performance. 3. 3.

95% CVG is the coverage probability of 95% equal tail prediction interval. 95% CVG values closer to $0.95$ indicate better prediction performance for the model. 4. 4.

95% ALCI is average length of 95% equal tail prediction interval. Smaller 95% ALCI values indicate better prediction performance for the model. 5. 5.

Continuous Ranked Probability Score (CRPS), which is defined as

[TABLE]

where $y$ is a scalar quantity that needs to forecast, and y admits an underlying distribution that is described by CDF $F$; $G$ is a CDF that is chosen by the forecaster in order to predict $F$. $H$ is a unit step function and $H_{y}(x)$ indicates a centered Heaviside function $H(x-y)$. CRPS is negatively oriented and the lower scores imply better performance.

Appendix F Simulation Study with Four Fidelity Levels

We consider a system with four levels of fidelity represented by the hierarchical statistical model (3.1) defined on a two dimensional unit square domain with univariate observation data sets $\mathbf{Z}=(\mathbf{Z}_{1},\mathbf{Z}_{2},\mathbf{Z}_{3},\mathbf{Z}_{4})$ with corresponding spatial locations $\mathbf{S}=(\mathbf{S}_{1},\mathbf{S}_{2},\mathbf{S}_{3},\mathbf{S}_{4})$. We assume a constant mean component for each level with design matrix $\mathbf{h}(\mathbf{S}_{t})=\mathbbm{1}$ and autoregressive coefficient function to be constant as $\gamma_{1}=1.1$, $\gamma_{2}=0.9$, $\gamma_{3}=1$. We chose isotropic exponential covariance function for the Gaussian process latent variables with parameters $(\sigma_{1}^{2}=2,\phi_{1}=12)$, $(\sigma_{2}^{2}=1,\phi_{2}=6)$, $(\sigma_{3}^{2}=0.8,\phi_{3}=8)$, and $(\sigma_{4}^{2}=0.5,\phi_{4}=3)$. Each of the levels are contaminated with white noise with variance $\tau_{1}=0.1$, $\tau_{2}=0.05$, $\tau_{1}=0.05$, and $\tau_{4}=0.01$. To make possible simulation of a Gaussian process with the specifications given in model (3.1), we constrain our overall sample size to $n=12,000$ and each level of fidelity $n_{1}=n_{2}=n_{3}=n_{4}=3,000$.

We generate a synthetic data set for the above statistical model. The data sets, shown in Figures 6a, 6b, 6c, 6d show the training data set for each of the four fidelity levels and 6e show the testing data from the fourth fidelity level. We apply the conjugate RNNC model to the generated data sets. We assign independent conjugate prior on parameters $\beta_{t}\sim N(0,1000)$ for $t=1,\dots,4$, and scale parameter $\gamma_{j}\sim N(0,1000)$ for $j=1,2,3$. We assign independent inverse gamma prior on spatial variance parameters $\sigma_{t}^{2}\sim IG(2,1)$, and on the noise parameters $\tau_{t}^{2}\sim IG(2,1)$ for $t=1,\dots,4$. We also assign uniform prior on the range parameters $\phi_{1,t}\equiv\phi_{2,t}\sim U(0,30)$ for $t=1,\dots,4$. The RMSPE with a 5-fold cross-validation was used for the conjugate RNNC model. We select $(\phi_{t},\tilde{\tau}_{t}^{2})$ on a grid such that $\phi_{t}$ ranges at $[0.1,25]$ and $\tilde{\tau}_{t}^{2}$ ranges at $[0.0005,0.4]$. As in the case of two fidelity levels, no significant differences were observed when we used 3-fold cross-validation and 7-fold cross-validation approach.

We compare our conjugate RNNC with the conjugate NNGP proposed in (Finley et al., 2019) using only the fourth single level training data. The prior specification for the conjugate NNGP are the same as above. Figure 7 gives the prediction mean and the prediction standard deviations using both methods. Figure 7a gives the single layer conjugate NNGP and compared to our proposed method Figure 7 is less accurate. Figure 7b shows the single layer conjugate NNGP prediction standard deviations which is considerably larger than the prediction standard deviations using conjugate RNNC. Moreover, the computed conjugate RNNC RMSPE is $0.64$ with 95% CVG $0.94$ and ALCI $2.37$. Instead, the computed conjugate NNGP RMSPE is $1.68$ with 95% CVG $0.86$ and ALCI $3.98$. These comparisons shows that accounting for the multi-fidelity dependencies can improve our prediction accuracy and its variation. The computational time for our proposed method is around four times slower than the conjugate NNGP in a single level. More precisely our method take around $117$ sec and the single layer conjugate NNGP $27$ seconds). This is very normal to expect since the computational complexity is linear to the observed data.