{'fraction_non_alphanumeric': 0.04329307742074617, 'fraction_numerical': 0.012346344619365969, 'mean_word_length': 4.495851851851852, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fused Gaussian Process for Very Large Spatial Data

With the development of new remote sensing technology, large or even massive spatial datasets covering the globe becomes available. Statistical analysis of such data is challenging. This article proposes a semiparametric approach to model large or massive spatial datasets. In particular, a Gaussian process with additive components is proposed, with its covariance structure consisting of two components: one component is flexible without assuming a specific parametric covariance function but is able to achieve dimension reduction; the other is parametric and simultaneously induces sparsity. The inference algorithm for parameter estimation and spatial prediction is devised. The resulting spatial prediction method that we call fused Gaussian process (FGP), is applied to simulated data and a massive satellite dataset. The results demonstrate the computational and inferential benefits of the FGP over competing methods and show that the FGP is more flexible and robust against model misspecification.


Introduction
With the advances of remote sensing technologies, massive scientific data can be collected over (a proportion of) the globe. Such spatially correlated data sets allow researchers to investigate various issues in environmental and atmospheric sciences. Classic statistical methods such as kriging have been widely used to model spatial data (Cressie 1993;Cressie and Wikle 2011;Banerjee et al. 2014). However, with large or massive data, direct implementation of these statistical methods becomes computationally prohibitive, since solving the kriging equations involves inversion of an n × n covariance matrix with data of size n, which requires computational cost O(n 3 ) and memory cost O(n 2 ) in general.
To tackle these issues, many recent developments in spatial statistics have focused on modeling large or massive spatial datasets. Most of them assume a specific form for the spatial covariance function known up to several parameters, e.g., the Matérn covariance function, and then use different approaches to represent or approximate this target function, the resulting covariance or precision matrix, or the likelihood function. Methods in this paradigm include approximate likelihood (Stein et al. 2004), covariance tapering (e.g., Furrer et al. 2006;Kaufman et al. 2008), predictive process (Banerjee et al. 2008) and its variants such as Sang and Huang (2012) and Katzfuss (2013), composite likelihood (e.g., Lindsay 1988; Eidsvik et al. 2014), nearest neighbor Gaussian process (Datta et al. 2016), Gaussian Markov random field representation (Lindgren et al. 2011), multiscale approximation for Gaussian process (Katzfuss 2016), and spectral methods (Duan et al. 2015; Guinness and Fuentes 2016;Stroud et al. 2016).
While the richness and flexibility of the methods mentioned above are indisputable, their implementation and performance generally rely on the assumption of a particular parametric form for the spatial covariance function. One of the main difficulties in using these methods to analyze massive data observed on a very large spatial domain such as the globe is to choose a specific covariance function that can represent various spatial structures in data, since misspecification of the spatial covariance function can have a large impact on inferential efficacy. Therefore, a family of covariance functions that are more flexible and robust are more desirable. Some extensions of the aforementioned methods to construct complex spatial processes have been proposed as well. By specifying a nonstationary Matérn covariance function in Paciorek and Schervish (2006), random spatial basis functions are obtained through reversible-jump Markov chain Monte Carlo algorithm, which can increase the flex-ibility of the resulting spatial process constructed from the predictive process in Katzfuss (2013). Another flexible and nonstationary covariance function is developed by adaptively partitioning the spatial domain through a Bayesian treed Gaussian process in Konomi et al. (2014). However, these extensions in general result in much more complicated algorithms, and the resulting computational cost can be too demanding for massive spatial data sets.
A second avenue of recent research has been focused on semiparametric modeling for analyzing large or massive spatial data sets. Methods in this paradigm represent the spatial process as a linear combination of multiresolutional basis functions and random coefficients (see, e.g., Fixed Rank Kriging in Cressie and Johannesson (2008), Gaussian process via local Karhunen-Loève expansion (Chu et al. 2014), and multiresolution Gaussian process (Nychka et al. 2015)). Specifically, a flexible family of nonstationary spatial covariance functions is developed in Cressie and Johannesson (2008) based on a pre-specified multiresolutional and compactly supported basis functions and a general lower-rank covariance matrix. Based on a novel local Karhunen-Loève expansion for a underlying spatial process, another type of flexible and unspecified spatial covariance functions is proposed and the consistency conditions for parameter estimators in the covariance function are established in Chu et al. (2014). Nychka et al. (2015) proposed to represent a spatial process as a linear combination of radial basis functions and random coefficients with precision matrix defined through a spatial autoregressive (SAR) model independently at multiple resolutions. These methods do not assume a parametric form for the spatial covariance function to attain more model flexibility and tend to be more robust against misspecification of the spatial covariance function.
However, these semiparametric models have their own limitations. The model in Cressie and Johannesson (2008) incorporates a relatively small number of basis functions to model the spatial field. The resulting low-rank covariance component allows for fast computation but may also incur sacrifices in capturing the spatial structure at various scales presented in data. By partitioning spatial domain into identically shaped subdomains and assuming independence among different subdomains in Chu et al. (2014), the resulting loss of accuracy needs to be recovered by other approaches in practice such as tapering approximation. Independence structure has been assumed among different resolutions for regularly spaced basis functions in Nychka et al. (2015), and hence can result in lack of flexibility in the resulting covariance function.
Our work combines the advantages of methods from both parametric and semiparametric paradigms. In particular, a Gaussian process is constructed with additive components, which have different basis representations based on a set of inducing variables. Thus, the resulting covariance structure of the latent Gaussian process is decomposed into two components.
The first component is characterized using a semiparametric representation with a relatively small number of basis functions. This low-rank component will not only enable dimension reduction and fast computation but also attain model flexibility. The second component in the covariance structure is defined through an undirected Gaussian graphical model (GGM), also called Gaussian Markov random field (GMRF). We will demonstrate that several stateof-the-art methods for parametric covariance approximation can be viewed as special cases of the GGM under various assumptions. The resulting process is called fused Gaussian process (FGP) as it blends two components that induce a low-rank covariance matrix and a sparse precision matrix, respectively. By taking advantage of the properties from both components, computationally efficient algorithms have been developed for parameter estimation and spatial prediction. Motivated by the small-block idea in Chu et al. (2014), the FGP is generalized to block fused Gaussian process (Block-FGP) by partitioning the spatial domain into blocks. The distributed parallel computing environments can be utilized to further accelerate statistical inferences for massive data in Block-FGP. The superior performance of the proposed approach is demonstrated through simulation studies and a real data example.
The remainder of the paper is organized as follows: Section 2 presents the semiparametric statistical model and discusses about relevant model specification and properties of the FGP. In Section 3, likelihood-based inference including parameter estimation and spatial prediction is devised for the FGP. Extension for the FGP is given and distributed computing algorithm for the Block-FGP is also discussed. Section 4 presents simulation examples as well as analyses with sea surface temperature data from NASA's Terra and Aqua satellites to illustrate the robust performance and inferential benefits of the proposed method. Section 5 concludes the paper with a brief summary and discussion on possible extensions and other applications.

The Fused Gaussian Process
This section starts with the definition the fused Gaussian Process and then discusses its properties and relationship with other state-of-the-art methods.

Low-rank Component S(·) η
The r-dimensional vector η is assumed to be a Gaussian random vector with zero mean and an unknown r × r covariance matrix K. The associated r basis functions S(·) are fixed and known. Such a low-rank component via a basis expansion has been widely used in analyzing large spatial data sets, and many forms of basis functions have been suggested, including local bisquare functions (Cressie and Johannesson 2008), wavelets (Shi and Cressie 2007), cubic B-splines (Chu et al. 2014), and basis functions resulted from a prespecified parametric covariance function and a set of prespecified locations or knots (Banerjee et al. 2008), among the others. In addition to prespecify basis functions, numerical methods for Karhunen-Loève (K-L) expansion has been developed to calculate eigenfunctions, which can be used as basis functions as suggested in Hu (2013). Recently, it has been pointed out in Bradley et al. (2016) that any class of basis functions can be re-weighted and then viewed as eigenfunctions within a K-L expansion, although sensitivity analysis is recommended to choose basis functions. As demonstrated in previous work, incorporating a low-rank basis expansion enables dimension reduction, and thus is able to facilitate computationally feasible inferences. Further discussion of the relationship and differences between our model and other low-rank models is given in the end of this section.

Gaussian-graphical-model component A(·) ξ
The M -dimensional vector ξ is assumed to be a Gaussian random vector with zero mean and nonsingular covariance matrix Σ. Suppose that Q ≡ Σ −1 = (q ij ) denote the corresponding precision matrix with (i, j)-th element q ij . Then the Gaussian random vector ξ can be represented by an undirected Gaussian graphical model (uGGM) via an undirected graph G = (V, E), where V contains the M vertices corresponding to the variables in ξ and the edges E = (e ij ) 1≤i<j≤M indicate whether the variables ξ i and ξ j (i = j) are conditionally independent given all other variables in ξ. Therefore, the variables in ξ are Markov with respect to G: where p(·) represents the probability density function; N i ≡ {j|j ∈ V, and {i, j} ∈ E}, and elements in the precision matrix Q are non-zero only for neighbors and diagonal elements: In spatial statistics, an uGGM is also called Gaussian Markov random field (e.g., Rue and Held 2005). Note that it is also possible to define a directed GGM (dGGM) for ξ. However, since a dGGM can be converted to an uGGM via moralization (i.e., "marriage" of a child node's parent nodes) (e.g., Jordan 2003), an uGGM for ξ is assumed in all numerical examples in this article, while brief discussion of extension with a dGGM is included in Section 5.
The precision matrix Q plays an important role in determining the dependence structure in ξ and is usually modeled as a large structured sparse matrix, known up to a few parameters. For example, an explicit parametric form for Q is provided in Lindgren et al.
(2011) when a Matérn covariance function is assumed. Another commonly used approach to construct Q is to assume a conditional autoregressive (CAR) model (Cressie and Wikle 2011, sec. 4.2). Specifically, the variance-covariance matrix of ξ in the CAR model takes the following form: is symmetric and positive-definite, the parameter γ needs to be between the reciprocal of smallest and largest eigenvalues of H (e.g., Besag (1974)). This CAR model on ξ implies the following conditional distributions, To specify the basis functions A(·) ≡ (A 1 (·), . . . , A M (·)) , suppose that the random vector ξ is defined on a generic lattice over the domain of interest, and zero otherwise. The matrices H and ∆ are also specified according to the neighborhood structure in the lattice. In the following numerical examples, H is constructed from first order neighborhood structure, and ∆ is chosen to be the identity matrix. In practice, {R i } can be determined by the finest resolution for which spatial predictions will be made, and choice of neighborhood structure can be made based on model selection criteria such as cross-validation or Bayesian information criterion (BIC). We assume that the observed data are obtained at the same or a coarser spatial resolution compared to the lattice {R i : i = 1, . . . , M }, and thus M can be much larger than the size of the data n. Such a lattice structure has been introduced and utilized to analyze remote-sensing data in previous studies (Nguyen et al. 2012).
The resulting model for the hidden process Y (s) is given by: which is called fused Gaussian process, since it combines the low-rank and graphical-model components. The following property is derived and its proof is given in Appendix A.
Proposition 1: The spatial process {Z(s) : s ∈ D} obtained from (2.6) is a valid Gaussian process with a positive semidefinite covariance function C Z (·, ·); and satisfies Kolmogorov's consistency condition, i.e., for every finite subset S ⊂ D, let f Z (S) denote the probability density function. Then the probability density function f Z satisfies Kolmogorov's consistency condition: (i) For every finite set S = {s 1 , . . . , s p } ⊂ D with arbitrary positive integer p ∈ N, and for every permutation π(1), . . . , π(p) of 1, . . . , p, we have that Remark: Note that for any p > M , the covariance matrix of (Y (s 1 ), . . . , Y (s p )) only has rank M . Therefore, the process Y (·) is a degenerate Gaussian process. However, since M is very large (M ≈ n or M > n) and is determined by the lattice defined by the finest resolution of interest, as demonstrated in numerical illustrations in Section 4, the FGP outperforms low-rank processes such as Fixed Rank Kriging in terms of spatial prediction. Meanwhile, unlike the Gaussian process whose predictive performance deteriorates with a misspecified covariance function, the proposed FGP is more flexible and able to provide robust spatial predictions.
This section ends with comparison between FGP with several models in literature for analyzing large or massive spatial data. The model in (2.6) contains a low-rank component similar to that in Cressie and Johannesson (2008). As pointed out in later work (e.g., Stein 2014), the performance of low-rank methods can be sensitive to the number of basis functions.
In Section 4.2, we compare the predictive performance of the FGP with Fixed Rank Kriging and show that by adding a Gaussian-Graphical-model component, the FGP is able to provide more accurate spatial predictions even with a small number of basis functions in the low-rank component. Meanwhile, it is worth noting that another two recently suggested methods for massive spatial data, the nearest-neighbor Gaussian process (NNGP) in Datta et al. (2016) and the multi-resoution approximation (M-RA) in Katzfuss (2016), can both be considered as models induced from a GGM component. In particular, the NNGP is induced by a directed Gaussian Graphical model with vertices related to locations in a pre-specified reference set, while the M-RA is built based on a multi-resolution Gaussian graphical model but has assumed independence between vertices from different resolutions and vertices from different clusters (or regions) within the same resolution. In the FGP, it is also possible to include such a GGM component as those in NNGP or M-RA, especially when a specific target parametric covariance function is desired for this component. Therefore, the NNGP and M-RA can be viewed as special cases of parameterization for the GGM component in the FGP, and the inference procedure and related computational advantages presented in Section 3 will still hold. Furthermore, note that many methods for large or massive spatial data, including predictive process and its variant ( rely on the assumption of a specific parametric covariance function C(·, ·). These methods are all designed to approximate or represent the corresponding Gaussian process with this assumed target covariance function C(·, ·). The simulation studies in Section 4.1 illustrate that misspecifying this covariance function in the Gaussian process can have a substantial impact on prediction efficacy, while the FGP is more flexible and robust.

Inference
This section presents how to obtain inferences about model parameters and spatial predictions with the FGP (Section 3.1) and gives the corresponding computational complexity calculations that demonstrate the effect of dimension reduction and utilization of sparse matrices (Section 3.2). Furthermore, the proposed FGP is extended to Block-FGP, which can be implemented in a distributed computing environment to further improve computational efficiency (Section 3.3).

Parameter Estimation and Spatial Prediction
Let θ denote the vector consisting of parameters in {β, K, τ 2 , γ}. Recall that the observed data is Z ≡ (Z(s 1 ), . . . , Z(s n )) . By combining (2.1) and (2.6) and assembling vectors into matrices, the spatial linear mixed effects model for Z can be written as the following matrix form: Up to an additive constant, the negative log-likelihood function is written as: . Evaluation of the negative log-likelihood function requires calculation of the inverse and log-determinant of the n × n matrix C.
Proposition 2: Recall that C ≡ var(Z) = SKS + AQ −1 A + V , it can be proved that: Proof of Proposition 2 can be found in Appendix A. Note that the right-hand sides of (3.1) and (3.2) involve only inversion and determinant of r×r matrices and M ×M sparse matrices, which enable fast evaluation of the negative log-likelihood.
To minimize the negative log-likelihood function l Z (θ), iterative algorithms are devised.
For example, as suggested in Chu et al. (2014), the matrix K is parametrized via its eigendecomposition and then a two-step iterative algorithm can be carried out. First, l Z (θ) is minimized with respect to the eigenvectors of K for a fixed (β, τ 2 , γ) using a Newton-Raphson-type algorithm on a Stiefel manifold (Peng and Paul 2009). Second, with fixed eigenvectors of K, l Z (θ) is minimized with respect to the remaining parameters. Another technique to minimize l Z (θ) is to treat random effects as "missing data", and the expectationmaximization (EM) algorithm (Dempster et al. 1977) can be implemented. To devise the EM algorithm for the FGP, the random vector η is treated as "missing data". Let θ t denote the parameters at the t-th iteration. In the expectation step (E-step), the conditional expectations and covariance matrix for η given the data Z and parameter estimates θ t are derived, respectively: In the maximization step (M-step), θ t+1 is updated by maximizing the so-called Q function obtained in the E-step.
In particular, closed-form updates can be derived for K and β, and numerical optimization procedures such as interior point method and active-set method, are implemented to update τ 2 and γ. To accelerate optimization procedure, parallel optimization algorithms (e.g., alternating direction method of multipliers method in Boyd et al. 2011) can be used in the M-step. Parallel algorithms can also be used to solve a large-scale sparse linear system For spatial prediction, suppose that we are interested in making prediction of Y (·) at a set of locations Conditioning on the parameter vector θ, the predictive distribution, Y P |Z, is derived as:

Computational Complexity
The main computational effort for inferences described in Section 3.1 is devoted to calculating the inverse and log-determinant of the n × n matrix, C ≡ SKS + AQ −1 A + V , in which K is only an r ×r matrix with r n, and Q is an M ×M sparse matrix with M ≈ n or M > n. and thus this matrix will never has memory cost more than O(M r). Unlike the full Gaussian process, the overall memory cost in FGP will never exceed O(M r) since log M is much smaller than r.

Extensions
Note that the computational complexity of the FGP is dependent on M , the size of the GGM component. To further improve computational efficiency with extremely large M (i.e., lattice at extremely fine spatial resolution), the small block method can be used to model ξ (see, for example, Vecchia 1988; Caragea and Smith 2006;Heaton et al. 2017). That is, the spatial domain is partitioned into blocks, and the corresponding graph is paritioned into J subgraphs. Then independence structure is assumed between these subgraphs for ξ. Therefore, the precision matrix Q becomes a block-diagonal sparse matrix, blockdiag(Q 1 , . . . , Q J ) i are the spatial association parameter and conditional marginal variance, respectively; H i is the proximity matrix, and ∆ i is the diagonal matrix with known positive elements for the ith subgraph, for i = 1, . . . , J. The resulting fused Gaussian process is called block fused Gaussian process. Although the elements in ξ are independent between blocks, it is worth noting that the Block-FGP maintains spatial dependence between blocks due to its low-rank component.
The inferences for the Block-FGP can be implemented in a distributed computing environment. According to the partition of blocks, the matrices X, S, A, and V can be decomposed as X = (X 1 , . . . , X J ) , S = (S 1 , . . . , S J ) , A = blockdiag(A 1 , . . . , A J ), and . . , J. The inverse of covariance matrix C of the data can be computed as follows: Then, up to an additive constant, the negative log-likelihood function is written as: where the quantitiesZ j D jZj , S j D jZj , S j D j S j , and log |D −1 j | can be computed in parallel for j = 1, . . . , J, which only require inversion and determinant of M j × M j sparse matrice Q j + A j V −1 j A j and r × r matrix K and thus requires computational cost O(M 1.5 j ) and O(r 3 ), respectively. In a distributed computing system with J worker nodes and one central node, to evaluate the negative log-likelihood function, for j = 1, . . . , J, the jth worker node will need to calculateZ j D jZj , S j D jZj , S j D j S j , and log |D −1 j |, and then store and transfer these quantities to the central node. The associated communication cost between the worker nodes and the central node is only for two scalers, one r × 1 vector, and one r × r matrix per worker node. The corresponding spatial predictions and associated standard errors can also be obtained through distributed computing system.

Numerical Examples
In this section, several simulation examples are provided to demonstrate the robust performance of the FGP and its computational efficiency. A massive sea surface temperature (SST) dataset is used to illustrate the inferential benefits of the FGP. All methods were implemented in MATLAB R2015b, and the MATLAB function fmincon was used to carry out numerical optimization. The simulation studies were carried out on a 4-core HP system with Intel Xeon x5650 CPU and 12 Gigabytes memory. The analysis of the SST data was carried out using 16 cores, and computation for the variance of predictions over entire region of interest is accomplished through parallel computing.

Simulation Examples
The robust performance of the FGP is illustrated in this section. Two different scenarios with two different covariance functions are considered. For comparison, two alternative methods are implemented in addition to the FGP. The first alternative is to perform kriging using the true covariance function, which is called EK. The second alternative is to perform kriging assuming the exponential covariance function regardless of the true underlying covariance structure, which is called MK.
In and ∆ is chosen to be the identity matrix. The EM algorithm described in Section 3.1 is utilized to estimate parameters and obtain spatial predictions based on formulas (3.3).
To implement EK, the true covariance parameters φ, σ 2 , and σ 2 are used to perform kriging with full covariance structure, which should give the best predictions and can be used as a baseline in the simulation experiment. To implement MK, an exponential covariance function model is assumed, but the parameters are estimated based on observed data and predictions are obtained through kriging.
To evaluate the predictive performance, the mean-squared prediction error (MSPE) is

adopted, which is defined as MSPE
where Y (s) is the simulated value andŶ (s) is the predicted one using procedure A = FGP, EK, or MK. We carry out L = 50 runs of simulations and summarize the mean and standard deviation of the In practice, we do not know the covariance structure and need to specify a model for it. Although several parametric covariance functions including exponential or Matérn functions have been widely used, they are sometimes chosen practically, and it may misspecify spatial dependence structure. In Scenario 2 of the simulation study, the potential impact of misspecified covariance structure is investigated, and the robust performance of the FGP is illustrated. Specifically, the same set-up as in Scenario 1 is used here, but the underlying process Y (·) is simulated from a Gaussian process with the sinusoidal covariance function c(h) = σ 2 sin(h/φ)φ/h with σ 2 = 16, and φ = 0.5. The same FGP is fitted using the observed data at n locations. For the procedure EK, the true sinusoidal covariance function with true parameters is used to perform kriging. For the procedure MK, an exponential covariance function is still assumed to carry out kriging, whose parameters are estimated based on observed data.
Based on L = 50 simulation runs, the mean and standard deviations of MSPE A are calculated for all three methods, A = FGP, EK, and MK. The results are summarized in Table 1.
As expected, EK performs the best, because it uses the true covariance structure as well as true parameters. Meanwhile, owing to the misspecified exponential covariance function, MK experiences severe deterioration of predictive performance; the relative efficiency for MK with respect to EK is only 29%. On the other hand, the method FGP outperforms MK by providing substantially better spatial prediction with smaller MSPE. Comparing FGP to the best procedure, namely EK, Table 1 shows that the performance of the FGP is satisfactory, and its relative efficiency with respect to EK is 86%. Meanwhile, the standard deviation of MSPE for FGP over the L = 50 simulation runs is also similar to that of EK and only about one sixth of that of MK, demonstrating the superiority of FGP compared to MK in this scenario.
Comparing the results under these two scenarios in Table 1, we notice that, as expected, performing kriging with the true covariance model and true parameters, namely EK, always gives the best prediction, but EK is certainly impractical since it requires true model and true parameters. Comparing the performance of MK in these two scenarios, it is clear that the specification of the covariance function can be influential: MK experiences a large deterioration of performance in Scenario 2 (where the covariance function is misspecified as an exponential covariance function) compared to Scenario 1 (where the covariance function is correctly specified as an exponential one); the average of MSPE increases substantially with the relative efficiency of MK with respective to EK reduces from beyond 95% to merely 29%, and the standard deviation of MSPE increases by a factor of seven. Other methods are expected to perform similarly as MK when the target covariance function is misspecified such as predictive process (Banerjee et al. 2008), GMRF (Lindgren et al. 2011), NNGP (Datta et al. 2016) or M-RA (Katzfuss 2016), because these methods are designed to approximate or represent the Gaussian process with assumed target covariance function (but with computational advantages). Now consider FGP: the same model is used under these two scenarios.
In both scenarios, the FGP gives satisfactory performance. It only experiences a small deterioration of performance in Scenario 2 compared to Scenario 1. Its relative efficiencies with respect to EK are 98% (Scenario 1) and 86% (Scenario 2), and the standard deviation of MSPE increases less than 10% in Scenario 2. In what follows, the computational advantages of the proposed method FGP is illus- In conclusion, the FGP provides superior predictive performance against misspecified covariance structure and can also achieve efficient computation. In the next subsection, the FGP is carried out to analyze a massive sea surface temperature dataset to further demonstrate its inferential benefits.

Application with Sea Surface Temperature Data
Sea surface temperature plays a vital role in the Earth's atmosphere and climate systems. A complete and accurate map of SST is essential in oceanographic sciences, weather forecasts, and in studying global and regional climate changes (Donlon et al. 2002;Zhang et al. 2009;Hirota et al. 2011;Lau and Waliser 2012). In this section, the performance of the FGP is When implementing the FGP to analyze this SST dataset, the trend term is modeled as a quadratic function in terms of the latitude. That is, the covariates X(·) = [1, latitude(·), latitude 2 (·)] are used in (2.6). For the low-rank component, as suggested by Cressie and Johannesson (2008), the multi-resolutional local bisquared basis functions are employed over the globe. These basis functions are defined as: seconds in CAR. It can be seen that with a more complicated model, the FGP does require more computing time compared to the other two methods, but it is capable of providing better and more reliable predictions. Table 2. Results from the validation study using the SST data: The M SP E A S are given for A = FRK, FGP, and CAR, and S = S 1 (missing by design), S 2 (missing at random) and S 1 ∪ S 2 (overall).  Figure 2 presents the predictions and associate prediction standard errors over this large region (upper panels) as well as zoomed-in maps for a subregion in Indian Ocean (lower panels), which clearly shows that the standard errors reflect the pattern of the missing data, and that the FGP captures spatial variation of SST very well.

Conclusion and Discussion
This article presents the fused Gaussian process that combines the low-rank component and an undirected Gaussian graphical model. The inference algorithms are devised for The larger rectangular region delineated shows the region where spatial predictions will be made. The smaller rectangular region delineated marks the area for which zoomed-in views are shown. Uppermiddle and upper-right panels: the spatial predictions and associated standard errors from the FGP, respectively. Lower panels: zoomed-in views of the data, spatial predictions and associated prediction standard errors over a small region, respectively, from left to right.
where S(R) ≡ (S 1 (R), . . . , S r (R)) ; A(R) ≡ (A 1 (R), . . . , A M (R)) ; S i (R) ≡ R S i (s) ds/|R|; and A i (R) ≡ R A i (s) ds/|R| for R ⊂ R d . Thus, the basis functions can be integrated offline and the formulas for spatial prediction and standard error will be of the same form. The current FGP model can be easily generalized for the space-time framework. For example, one may assume a spatio-temporal random-effect model ν(s, t) = S t (s) η(t), where the basis function S t (·) depends on time t = 0, 1, . . ., and {η(t) : t = 0, 1, . . .} follows an r dimensional vector autoregressive (VAR) model  Ren and Sun 2012) can be assigned to τ 2 and γ. When using the FGP for data fusion, we can generalize the work in Nguyen et al. (2012) and also include a mutlivariate CAR model in the FGP. These topics are currently under investigation. density function. Then the probability density function f Z satisfies Kolmogorov's consistency conditions: (i) For every finite set S = {s 1 , . . . , s p } ⊂ D with arbitrary positive integer p ∈ N, and for every permutation π(1), . . . , π(p) of 1, . . . , p, we have that f Z (s 1 , . . . , s p ) = f Z (s π(1) , . . . , s π(p) ).
Proof of Proposition 1 : Let {s 1 , . . . , s p } be every subset of the spatial domain D for arbitrary integer p. The p-dimensional random vector of Z p ≡ (Z(s 1 ), . . . , Z(s p )) follows multivariate normal distribution since a linear combination of Gaussian random vectors also follows multivariate normal distribution given the fact that Z p = µ p + S p η + A p ξ + p , and thus Z(·) is a Gaussian process. The covariance function of Z(s) and Z(u) is given by for any s, u ∈ D. For any real numbers k i ∈ R, i=1, . . . ,p, it follows from the positive definiteness of r × r matrix K and M × M matrix Q that where S p is p × r basis function matrix; A p is p × M basis matrix; and k p ≡ (k 1 , . . . , k p ) .
So, Z(·) is a valid Gaussian process.
To check the consistency condition for Gaussian process Z(·), it suffices to verify that the covariance matrix C Zp of random vector Z p is positive definite. Notice that C Zp ≡ S p KS p + A p QA p + V p , which is clearly positive definite.

Remark:
The process Y (·) is a degenerate process, since for any p > M , the covariance matrix of (Y (s 1 ), . . . , Y (s p )) only has rank M .
Proposition 2: Recall that C ≡ var(Z) = SKS + AQ −1 A + V , we have: log |C| = log |K −1 + S DS| + log |K| + log |D −1 |, Proof of Proposition 2 : For any n × n invertible matrix E and r × r invertible matrix F, we have Sherman-Morrison-Woodbury formula (Henderson and Searle 1981), Then the first equation for matrix inversion can be derived by letting E = D, U = S, V = S , and E can be inverted using Sherman-Morrison-Woodbury formula as well. According to the Matrix Determinant Lemma, the determinant of matrix E + UFV can be calculated as follows: |E + UFV| = |F −1 + VE −1 U||F||E|.
Thus, the log-determinant equation can also be easily verified after taking the logarithm of both sides.

B The EM Algorithm for Fused Gaussian Process
This section gives details for EM algorithm to estimate parameters. Recall that the variance of measurement error is not estimated in EM algorithm, since this quantity is usually known from the experiment in advance. If it is unknown, a straight line can be fitted near origin for the empirical semivariograms as suggested in Kang et al. (2010). We treat the random vector η as "missing data". The EM algorithm attempts to maximize the complete log-likelihood function ln L(θ|η, Z) iteratively, by replacing it with its conditional expectation of η given the observed data Z. The complete twice negative log-likelihood function is given by −2 ln L(θ|η, Z) = (n + r) ln(2π) + ln |AQ −1 A + V | + ln |K| Given the parameter estimates θ t , the EM algorithm consists of an E-step followed by an M-step defined as follows, for t = 0, 1, . . .: where µ η|Z,θt = E(η|Z, θ t ) = K t S C −1 t (Z − Xβ t ), Σ η|Z,θt = var(η|Z, θ t ) = K t − K t S C −1 t SK t , with C t ≡ SK t S + AQ −1 t A + V , Q t ≡ ∆ −1 (I − γ t H)/τ 2 t , and D ≡ (AQ −1 A + V ) −1 M-step: Find θ t+1 in parameter space Θ such that θ t+1 = Arg sup θ∈Θ Q(θ; θ t ).
The minimization of the function f (τ 2 , γ) is carried out in each iteration of EM algorithm, which requires more time for EM algorithm to converge. In order to make the EM algorithm converge faster, the SQUAREM algorithm can be used for parameter estimation, which accelerates EM algorithm through Akein's acceleration (for details, see Berlinet and Roland 2007;Varadhan and Roland 2008).
The initial value for parameter τ 2 in EM algorithm can be set to 0.1σ 2 Z withσ 2 Z representing the empirical variance in the data Z. The initial value for γ is restricted to the fixed interval (1/λ 1 , 1/λ M ), where λ 1 , λ M are smallest and largest eigenvalues for the proximity matrix H. Initial value for matrix K can be set to the r × r diagonal matrix 0.9σ 2 Z I r , which can ensure that the covariance matrix K is updated with a positive definite matrix in every iteration of EM algorithm. Given initial values for the parameters, the convergence of EM algorithm can be determined by monitoring the negative log-likelihood function.