Minimizing Negative Transfer of Knowledge in Multivariate Gaussian Processes: A Scalable and Regularized Approach

Recently there has been an increasing interest in the multivariate Gaussian process (MGP) which extends the Gaussian process (GP) to deal with multiple outputs. One approach to construct the MGP and account for non-trivial commonalities amongst outputs employs a convolution process (CP). The CP is based on the idea of sharing latent functions across several convolutions. Despite the elegance of the CP construction, it provides new challenges that need yet to be tackled. First, even with a moderate number of outputs, model building is extremely prohibitive due to the huge increase in computational demands and number of parameters to be estimated. Second, the negative transfer of knowledge may occur when some outputs do not share commonalities. In this paper we address these issues. We propose a regularized pairwise modeling approach for the MGP established using CP. The key feature of our approach is to distribute the estimation of the full multivariate model into a group of bivariate GPs which are individually built. Interestingly pairwise modeling turns out to possess unique characteristics, which allows us to tackle the challenge of negative transfer through penalizing the latent function that facilitates information sharing in each bivariate model. Predictions are then made through combining predictions from the bivariate models within a Bayesian framework. The proposed method has excellent scalability when the number of outputs is large and minimizes the negative transfer of knowledge between uncorrelated outputs. Statistical guarantees for the proposed method are studied and its advantageous features are demonstrated through numerical studies.


Introduction
Gaussian process regression models are widely used in several fields due to their desirable properties, such as flexibility, ease of implementation, uncertainty quantification and natural Bayesian interpretation (Quiñonero-Candela and Rasmussen, 2005). Recently, there has been increasing interest in extending GP models to deal with multivariate outputs (also known as cokriging) due to their prevalence in many applications. For example, in manufacturing plants, hard to sample performance indicators can be predicted from other correlated and cheap to sample indicators (Osborne et al., 2008). Also, the future evolution of sensor signals from in-service devices can be predicted using previously observed signals from similar devices in the historical database (Kontar et al., 2017a). Other applications arise in geostatistics, wireless networks and computer experiments.
The multivariate Gaussian process draws its roots from multitask learning. When multiple datasets from related outputs exist, integrative analysis can be advantageous relative to learning outputs independently. This integrative analysis, which leverages commonalities among related outputs to improve prediction and learning accuracy is referred to as multitask learning (Caruana, 1998;Yuan et al., 2012). The key feature in multitask learning is to provide a shared representation between training and testing outputs to allow the inductive transfer of knowledge. From an MGP perspective, this transfer of knowledge is achieved through specifying a valid positive semidefinite covariance function that models the dependencies of all data points within an output and across different outputs (Conti et al., 2009).
Traditionally, in MGP models, outputs are jointly modeled with a separable covariance structure, i.e. correlation over the same the input space is separable from between-output correlation (Conti and OHagan, 2010;Qian et al., 2008;Zhou et al., 2011). In such methods, the separable covariance function is of the form t × cov(x, x ) where t is the between-output covariance matrix and cov(., .) is a covariance function over inputs x ∈ R D , the same for all outputs. This assumption is appealing due to the simplified covariance structure and significant reduction of model parameters, however it restricts all marginal GP's for all outputs to share the same set of covariance parameters defined in cov (., .).
On the other hand, nonseparable covariance functions allow outputs to possess both shared and unique features, as different outputs share information through different covariance parameters, therefore, accounting for non-trivial commonalities in the data. Recent work on nonseparable covariance functions are mainly based on convolution processes (CP) (Majumdar and Gelfand, 2007;Fricker et al., 2013;Melkumyan and Ramos, 2011). Earlier work, is known in the geostatistics literature as the linear model of coregionalization (LMC) (Goulard and Voltz, 1992) and can be seen as a special case of the CP framework (Álvarez et al., 2010). The CP is based on the idea that a GP can be constructed by convolving a latent function with a smoothing kernel (Ver Hoef and Barry, 1998). Thus, if each output is expressed as a convolution of a latent function drawn from a GP, and if we share these latent functions across multiple convolutions, then, multiple outputs can be expressed as a jointly distributed GP (Boyle and Frean, 2005;Alvarez and Lawrence, 2009). As referred to byÁlvarez and Lawrence (2011), the key feature of the CP approach is facilitating the non-instantaneous mixing of base processes, where, for instance, each output can be described using its own length-scale, which gives an added flexibility for describing the data.
Despite the elegance of the MGP established using CP, denoted as MGCP, model building is extremely prohibitive and often impractical even with a moderate number of outputs. Further, it is not uncommon for modern engineering systems to have a large number of outputs (Higdon et al., 2008;McFarland et al., 2008). For instance, when considering qualitative factors, the number of outputs in computer experiments increases dramatically based on all potential combinations (Han et al., 2009). This leads to a set of considerable challenges.
• Challenge 1 (Computational complexity): The fact that the full covariance function of the joint GP should be considered results in significant computational burden and numerical issues. This challenge, in fact, is a direct consequence of accounting for multiple outputs and has been recently tackled in some literature (Álvarez and Lawrence, 2011;Kontar et al., 2017b).
However, two other challenges arise with the CP construction and have yet to be tackled.
• Challenge 2 (High dimensional parameter space -Number of parameters): The flexibility of the CP is based on providing different covariance parameters for different output levels, therefore, even for a moderate-scale problem, the number of parameters in the covariance function can easily reach hundreds or even thousands. This will lead to significant difficulties in solving the optimization problem to find the maximum likelihood estimator (MLE) for the parameters, specifically under non-convex and highly nonlinear settings.
• Challenge 3 (Negative transfer of knowledge): The integrative analysis of multiple outputs implicitly assumes that these outputs share some commonalities. However, if this does not hold, negative transfer of knowledge may occur, which leads to decreased performance relative to learning tasks separately (Pan and Yang, 2010). This is specifically important in the CP approach which, unlike separable approaches, implicitly implies that functions have unique features. Even though negative transfer is a very important issue, no research has handled this issue in the context of MGP models.
In the current literature, nonseparable modeling using MGCP is prohibitive even for a moderate number of outputs, due to challenges 1 and 2. Also, no literature has addressed the negative transfer of knowledge (challenge 3 ) in MGCP that results from the integrative analysis of outputs that share no commonalities. This article aims to simultaneously overcome these challenges through proposing a regularized pairwise modeling approach for MGCP models. The proposed method has excellent scalability when the number of outputs is large and minimizes the negative transfer of knowledge between uncorrelated outputs. Our approach is based on breaking down the high dimensional MGCP into a group of bivariate GP models, where the shared latent function parameters that facilitate the sharing of information in each bivariate model are penalized to prevent information sharing between outputs with no commonalities. Consistency in estimation and variable selection are then established. Empirical evidence demonstrates that the proposed method can: (1) achieve similar prediction performance as the full multivariate approach when the output dimension is low, (2) outperform the full multivariate approach, with only a fraction of its computational needs, when the output dimension is high, (3) outperform the full multivariate approach when some functions are uncorrelated even when the output dimension is low.
The rest of the article is organized as follows. Section 2 provides some preliminaries related to the CP construction. In Section 3, we motivate the proposed method through expanding on the proposed challenges. Section 4 introduces our regularized pairwise modeling approach and proves some corresponding statistical properties. The advantageous features of our proposed method are then demonstrated through benchmarking our method with other reference methods in Section 5. Finally, Section 6 concludes this article with discussions. Technical details are deferred to the appendix.

Preliminaries
Consider a set of N output functions y(x) = [y 1 (x), y 2 (x), ..., y N (x)] t from inputs x lying in some input space X ⊂ R D . The MGP is defined as where the function F : , ..., x ip i ] t and p i represents the number of observations for output i. Also, denote X = [X t i , ..., X t N ] t to be the matrix of all input points. We do not restrict X i = X j . The formulation in (1) is a general decomposition of an MGP and will reduce to a univariate GP if N = 1. As shown in Figure 1, the underlying principle of the MGP is to borrow strength from a sample of curves, through the shared representation in (1) which enables integrative analysis, in order to predict individual outputs. For instance, assuming P = p i , at any new input x 0 ∈ X belonging to output i ∈ I the joint distribution of the observed values from all outputs and the target function value at the test location y 0 i = y i (x 0 ) is given by

General Multivariate GP Setting
where y = [y t 1 , y t 2 , ..., y t N ] t are the noisy observed targets corresponding to the latent func- is the covariance matrix relating all input points for all outputs with cov . Following multivariate normal theory, the predictive distribution of y i (x 0 ) denoted as pr(·|y) is given as The mean in (3) is the empirical best linear unbiased estimator (EBLUP) of y i (x 0 ), while the variance is divided into three parts, the first is the variance of the variable under study, σ 2 i represents additive noise and the last part is the variance reduction due to the EBLUP approximation (Stein and Corsten, 1991).
As shown from (3), the sharing of information is achieved through cov f ij (x, x ) which models the variations both within and across different outputs, to capture their relatedness and improve prediction accuracy. This covariance function is typically assumed to belong to a known parametric class (Calder and Cressie, 2007) is a valid positive semidefinite covariance function, such that Θ f is a set that contains the true parameters θ * f for the covariance of f (x). Although in univariate settings there are many well known valid autocovariance functions, however, in the MGP it is extremely challenging to define crosscovariance functions that result in valid covariance matrices (Boyle, 2007).
An alternative to directly parametrizing a covariance function is to construct a GP, f i : R D → R, through convolving a Gaussian white noise process X(x) with a smoothing kernel K i (x) = α i k i (x) where α i ∈ R and k i : R D → R. Since the base process is a GP, and a convolution is a linear operator on a function, then the convolved process is also a GP (Thiébaux and Pedder, 1987;Barry et al., 1996;Higdon, 2002).
where denotes a convolution, cov(X i (u), X i (u )) = δ(u − u ) = δ(u − u) = δ uu and δ is the Dirac delta function. This approach is equivalent to applying a stable linear filter, where the output f i (x) is a weighted integral over the input signal X(x), weighted according to the impulse response K i (x). This requires the filter to be stable, where the output is bounded for all bounded input signals (Haykin, 2008), i.e. for a positive real finite number a, Therefore the only restriction for constructing a valid GP is that the impulse response/kernel is absolutely integrable ∞ −∞ |k i (u)|du < ∞. Some applications and extensions of the CP approach for the single ouput case are presented in Wikle (2002) and Calder and Cressie (2007).
Under the CP construction as shown in (4), if we share the same latent function X(x), across multiple outputs f i (x), i ∈ I, then all outputs can be expressed as a jointly distributed GP, i.e. MGP (Higdon, 2002). The resulting covariance function will then only depend on the displacement vector d = x − x ∈ R D and is given as As shown in (5), instead of directly parametrizing a positive semi-definite covariance function we only need to specify the parameters of a smoothing kernel, where the resulting covariance function must be valid by construction. Therefore, the key advantages is that we can exploit influence of multiple latent functions, X q (x) where q ∈ {1, 2, .., Q}, to share information across different output levels through different covariance parameters encoded in the kernels K qi (x) as shown in (6) (Fricker et al., 2013;Alvarez et al., 2012).

Challenges and Motivation
In this section, we expand on the MGCP challenges and present the motivation for our proposed method.

Computational complexity
As shown in the previous section, the MGCP is fully parametrized through the kernel parameters θ f and measurement noise σ = {σ 1 , ..., σ N }. Now denote θ = {θ t f , σ t } t and all the observed data as D = {D i , ..., D N }, then the likelihood of the joint MGCP is given as As shown in (7), learning from the likelihood requires the inversion of C f ,f + Σ and calculating its determinant at each step/iteration. The covariance matrix of the CP, assuming p observations for each of the N outputs, scales as N p leading to O(P 3 ) = O(N 3 p 3 ) computational complexity and O(N 2 p 2 ) storage. Therefore, CP modeling approaches are plagued by extremely high computational loads and numerical issues. In addition, more input points are usually required for multivariate outputs, resulting in a further increase in the complexity. Although separable approaches consider multiple outputs however their restrictive covariance functions often lead to structured covariances which are easily manipulated (Qian et al., 2008;Zhou et al., 2011). On the other hand, some recent approaches have tried to tackle this computational issue in the context of nonseparable covariances (Álvarez and Lawrence, 2011;Nguyen et al., 2014). For instance,Álvarez and Lawrence (2011) proposed an inducing variable approximation similar to the well-known partially independent training conditional (PITC) approximation (Quiñonero-Candela and Rasmussen, 2005). Their key assumption is that that outputs {y i (x) : i ∈ I} would be conditionally independent if all latent functions X q (u) are observed at some inducing locations, i.e. pr( ). This assumption leads to the inversion of a block diagonal covariance matrix which reduces the complexity to O(N p 3 ) which matches that of modeling with N independent GPs. Further, Kontar et al. (2017b,a) assumed training output as independent and proposed sharing latent functions X q (u) only between test and training outputs. This lead to a block arrowhead covariance matrix which also reduced the complexity of learning from the likelihood to that of modeling with independent GPs. It crucial to note that despite the reduced computational complexity due to the manipulation of the covariance matrix, we are still required to estimate an extremely large number of parameters even with a moderate number of outputs. This challenge is further discussed in the next subsection.

Number of parameters
Despite the disadvantages associated with high computational complexity, the main drawback, which renders nonseparable modeling prohibitive even for a moderate N , is the large number of the parameters that need to be estimated. This large number of parameters is a direct consequence of the CP construction which inherits its flexibility from the ability to share information through different kernels resulting in different covariance parameters for different output levels. It is clear from (3) that prediction accuracy is greatly dependent on the parameter estimates which are obtained from minimizing the negative log-likelihood function (θ; D) = −log L(θ; D). Up to an additive constant, and Y = yy t . Since, (θ; D) is highly nonlinear and non-convex in θ, its minimization in a high dimensional parameter space is extremely challenging, time consuming and suboptimal as one should anticipate poor parameter estimates in such high dimensions (Rasmussen, 2004;Tajbakhsh et al., 2014;Mardia and Watkins, 1989).  (6), the number of parameters depends on Q which is the number of latent functions induced in the model. For instance, consider the case in Figure  2 where we share latent functions between each pair of outputs, i.e. Q = N (N − 1)/2 resulting in N (N − 1) kernels. This is comparable to two way interaction effects in ANOVA where (Wonnacott and Wonnacott, 1990) only pairwise interactions are considered through the shared latent GP. Now assume the kernels K iq (x) = α qi k qi (x) follow the most commonly used exponential kernel, k qi (x) = exp((x − x ) t Λ qi (x − x )) where Λ qi is a D × D positive definite diagonal matrix allowing different length scales for each dimension. As a result the total number of parameters in the model is N (N −1)(1+D)+N where the first part N (N −1) represents the number of kernels multiplied by the number of parameters (1 + D) in each kernel, while the second part N represents the number of parameters in σ. Note that this case is a bit conservative as we use the exponential kernel which is able to provide a large degree of flexibility with a small number of hyper parameters (Rasmussen, 2004;Colosimo et al., 2014). Even for a moderate scale case where N = 30 and D = 1, we are required to estimate 1770 parameters under a non-convex setting. In another case, considering the more restrictive approach in Li and Zhou (2016) and Li et al. (2018) where Q = 1 and all outputs possess the same noise, i.e. σ = σ 1 =, .., = σ N , the number of parameters still scales as N (N + 2D + 1)/2 + 1, therefore for N = 30 and D = 1 we are estimating 991 parameters. In conclusion, obtaining good estimates in such a high dimensional parameters space is an impractical task, for this reason the practical applications of the MGCP are limited. Note that once the parameters are learned, prediction complexity at any new test point x 0 is O(N p) for mean and O(N 2 p 2 ) for variance, which can be done rather efficiently.

Negative transfer
As previously mentioned, a major challenge in the MGP which to the best of our knowledge has not been tackled yet, is the negative transfer of knowledge that occurs when we integratively analyze outputs that share no commonalities. Similar to the second challenge, negative transfer is specifically important in nonseperabale approaches which are used when functions have unique features. We start with an example to illustrative the impact of negative transfer in MGCP models. Consider a simple case with two outputs and a one dimensional input x ∈ X ⊂ R. The outputs are generated according to y 1 (x) = 1 + sin(x) + 1 (x) and y 2 (x) = 4 + 0.5sin(1.5x) + 2 (x) and x ∈ [0, 10]. The number of observations per signal is p 1 = p 2 = 20 evenly spaced points and the measurement noise is set as σ 1 = σ 2 = 0.1. We analyze this data separately using a univariate GP applied to each output and integratively using an MGCP. In the univariate GP, we assume a Gaussian/squared exponential covariance where cov f ii (x, x ) = α 2 i exp(−d 2 /(4ν 2 )). In the MGCP we adopt the covariance function in Fricker et al. (2013) and Li et al. (2018), where . Results are shown in Figure 3. The results clearly indicate that separate modeling of each function is significantly better than there integrative analysis. This is specifically clear for y 2 , which interestingly, was predicted using a larger length scale than the truth due to sharing of information with a smoother function y 1 . Note that this problem occurred in our dense input example (p i = 20) which implies that the challenge of negative transfer becomes exceedingly significant with sparse data. For these reasons, and since negative transfer is a crucial issue in MGP models, the main goal of this article is to handle the negative transfer of knowledge while maintaining the scalability of the model.

Motivation: Pairwise estimation and the precision matrix
In this article we propose a pairwise distributed estimation scheme motivated by both the distributed estimation literature for univariate GP models (Tresp, 2000;Deisenroth and Ng, 2015) and pairwise modeling of longitudinal profiles (Fieuws and Verbeke, 2006;Li and Zhou, 2016). The proposed approach is based on distributing MGCP estimation through bivariate GP submodels which are individually estimated. Predictions are then made through combining predictions from the bivariate models within a Bayesian framework. Not only does this approach scale to arbitrarily large datasets by parallelization, also each bivariate model can be efficiently built with a limited of parameters and a small-scale covariance matrix.
While, pairwise modeling seems to address computational challenges, negative transfer remains an important issue. Interestingly, pairwise modeling turns out to possess unique characteristics which allows us to tackle the challenge of negative transfer. While few literature (Williams et al., 2009) have aimed to establish the number of latent functions to be shared, such approaches do not imply avoiding negative transfer. As we will show in this subsection, information sharing and independent predictions can only be avoided through sparsity on the precision matrix. However, since MGCP models are based on modeling the covariance through latent functions not the precision matrix then we can only control the precision matrix under specific structures. It is clear from (3) that prediction accuracy for the GP/MGP is dependent on the inverse covariance matrix, also known as precision matrix Ω = (C f ,f + Σ) −1 ∈ R P ×P . The precision matrix carries conditional independence information. This matrix consists of block matrices Ω ij ∈ R p i ×p j , where the (c, c ) th entry of each block is denoted as Ω c,c ij . One can directly show that cov(y c i , y c j |ÿ) = 0 if and only if Ω c,c ij = 0, whereÿ = {y}/{y c i , y c j } (y excluding y c i and y c j ). Thus, conditionally independent variables lead to zero entries in the precision matrix (Whittaker, 2009). As mentioned previously, GP/MGP models are characterized through a positive semidefinite covariance function (ex: cov f ij (x, x )) rather than a conditional covariance function to generate the precision matrix. However, the remarks below illustrate some useful properties in the case of a bivariate GP.
then Ω ij = 0 if and only if, the multivariate Gaussian random vectors y i and y j are inde-pendent, i.e. cov(y c i , y c j ) = 0 for every c ∈ {1, .., p i } and c ∈ {1, .., p j }.
A brief proof of both Lemmas and further references are provided in Appendix A. The key conclusion from the Lemmas is that through pairwise modeling, we are able to control the precision matrix through parameters in the covariance function used to construct the bivariate GP.
The utilization of pairwise modeling, distributed estimation and this direct mapping from the covariance function to the precision matrix is detailed in the following sections.

Model Development
The proposed framework presents a flexible alternative that can scale to a large number of outputs while avoiding the negative transfer of knowledge. The nature of our proposed pairwise approach circumvents any need to find or establish latent functions between pairs. While within each pair, model selection is automatically done through our regularization approach which is consistently able to infer whether information should be shared or not. In Section 4.1, we establish our pairwise model based on a CP construction. Our pairwise scheme is based on distributing the estimation of the high dimensional MGCP into bivariate GP's which are individually built. An MGCP with N outputs as a result decomposes into N (N − 1)/2 pairwise submodels to predict all outputs. However, for the sake of notational simplicity, and building on Figure 1, we focus on predicting one output through sharing information from the remaining N − 1 outputs as shown in Figure 4 below. In Section, 4.2 we provide some statistical guarantees for the proposed method. Section 4.3, provides a direct approach to applying our regularized pairwise approach to separable modeling. Finally, Section 4.4, provides the methodology to combine predictions from the bivariate submodels.

General Multivariate GP
Further, we assume that y ij = [y t i , y t j ] t represents the noisy observations corresponding to the latent function values In order to capture both the unique properties of each output and their interdependence we construct both f i (x) and f j (x) as the sum of a latent function unique to each output and a shared latent function which facilitates the sharing of information. The model structure is shown in Figure 5. As shown in Figure 5, independent features are encoded through dependence on a latent function that has no effect on the other output (X i , X j ) while dependent features are encoded through the common dependence on X 0 .

Methodology
Following the CP construction in (6) This model is quite flexible, as it provides both shared and unique latent functions for both outputs. Based on this modeling framework, the cross covariance function between the two outputs is given as where, τ ij is the Kronecker delta function, which is equal to one if i = j, and is zero otherwise. In a more general case to that of Section 2, we define the latent functions {X q : q = 0, i, j} as cov(X q (u), X q (u )) = ξ 2 q δ uu where ξ q ∈ R. Then given the fact that X q (u) and X q (u ) are independent latent functions which only covary if q = q and u = u and utilizing the commutativity of the convolution and the "sifting" property of the Dirac Note that the derivation in (9) is a wrapper function for cov f ii (x, x ), cov f jj (x, x ) and cov f ij (x, x ). Based on our model framework illustrated in Figure 5 and due to the shared dependence on only X 0 , the autocovariance and cross covariance functions simplify as This modeling framework is generic in terms of the choice of the kernel function. However, a general purpose kernel can be constructed through assuming the kernels follow a Gaussian form. As mentioned previously, the Gaussian kernel is the most common choice of kernels due to its flexibility and correspondence to Bayesian linear regression with an infinite number of basis functions (Rasmussen, 2004). Further, similar constructions using such kernels have been utilized in both a GP and MGP setting (Paciorek and Schervish, 2004;Álvarez and Lawrence, 2011).
, to be the identity for the product of two Gaussian distributions. Then, we have that where (11), is that the marginal process i.e. cov f ii (x, x ) = q={0,i} ξ 2 q α 2 qi exp(− 1 4 d t Λ qi d) has the most common Gaussian covariance function resulting from the convolution of two Gaussian kernels. Therefore, (11) can be viewed as the extension of the Gaussian covariance function to the multivariate case. Once again we note that (11) is a wrapper function where for i = j, we have that to be the set of all parameters in the bivariate submodel. Then, given our CP formulation and following (8), the marginal density of the bivariate submodel is expressed as pr where pr(f ij |θ f ij ) = N (0 P , C f ij ,f ij ), pr(y ij |f ij ) = N (0 P , Σ ij ) and p = p i + p j . Therefore such that C f ij ,f ij ∈ R p×p is the covariance matrix relating all input points of outputs i and j with cov f ij (x, x ) in (10) and (11). As previously mentioned, parameter estimates are obtained from minimizing the negative log-likelihood function (θ ij ; D i , D j ) = −log pr(y i , y j |X 1 , X 2 , θ ij ). Denoting Y ij = y ij y t ij , up to an additive constant, the bivariate likelihood and its derivatives are given as Further, through denoting the gradient and second derivatives with respect to any parameter θ (n) ij ∈ θ ij are then given as (in Appendix B we expand on ∂C ij /∂θ The computational complexity of learning from the bivariate likelihood in (14) is reduced to O((2p) 3 ). More importantly, the total number of parameters in the model is reduced to 4(1 + D) + 2 + 3 where 4 represents the number of kernels (K ii , K jj , K 0i , K 0j ) multiplied by the number of parameters (1 + D) in each kernel, while 2 represents σ i and σ j , and 3 represents the parameters (ξ 0 , ξ i , ξ j ) of the latent functions (X 0 , X i , X j ) respectively. Note that reductions in parameter number can be done through assuming K 0i = K 0j or σ i = σ j . As shown distributed estimation using bivariate submodels, is able to handle both the computational complexity and large parameter number to be estimated. Also, all pairwise models can be parallelized and thus our model can scale to an arbitrarily large number of outputs by parallelization. Now, recall Figure (5), X 0 , defined by ξ 0 , represents the latent function which facilitates the sharing of information between outputs i and j. Therefore, in order to handle negative transfer of knowledge we use the bivariate likelihood in (14) but with ξ 0 penalized. Following Lemma 2, we will show in the following section that shrinking ξ 0 decreases the cross correlation amongst the outputs and that ξ 0 = 0 ensures that each output is predicted independently. The penalized negative log-likelihood function is defined as where P λ (ξ 0 ) is a penalty function. Different types of penalty functions can be used, examples include: ridge penalty P λ (|ξ 0 |) = λξ 2 0 , 1 penalty P λ (|ξ 0 |) = λ|ξ 0 |, bridge penalty P λ (|ξ 0 |) = λ|ξ 0 | 0<·<1 , and scad penalty P λ (|ξ 0 |) = λ|ξ 0 | if |ξ 0 | ≤ λ, (ξ 2 0 −2γλ|ξ 0 |+λ 2 )/(2γ−2) if λ < |ξ 0 | ≤ γλ, λ 2 (γ + 1)/2 if |ξ 0 | > γλ. The tuning parameter λ (λ and γ in Scad) has an important effect on predictions. For instance in the Lasso, as λ increases, ξ 0 shrinks to zero. Typically, the optimal tuning parameter is found using a grid search such as generalized cross validation (GCV), specifically b-fold GCV (Friedman et al., 2001). The GCV is based on splitting the data into b groups. For a given λ, first, we exclude one group as the testing dataset, and the rest b − 1 groups are used as the training data, from which predictions at test locations are obtained. This is repeated b times for each group. The tuning parameter is then chosen as the minimizer of some accuracy measure (ex: absolute error, mean squared error, etc..) based on the GCV procedure.

Statistical properties
We now discuss some of the structural properties of our model. We first introduce the main theorem which guarantees that our regularized model in (15) is able to avoid the negative transfer of knowledge.
Theorem 1. Suppose that ξ 0 = 0, then the predictive distribution of the bivariate model at any new input x 0 ∈ X reduces to that of a univariate model where pr(y i (x 0 )|y ij ) = pr(y i (x 0 )|y i ) and pr(y j (x 0 )|y ij ) = pr(y j (x 0 )|y j ) such that The proof is detailed in Appendix C. The key feature of this theorem is that penalizing only one variable, in our initial parameter set ξ 0 ∈ θ f ij ⊂ θ ij , will lead to separating the bivariate model into two models equivalent to the univariate GP established through a CP in (4). Our regularization approach is flexible to any specified kernel function and not based on the Gaussian covariance derived in (11), where theorem 1 holds for any valid kernel function K iq . We note that the result of theorem 1, is based on the fact that for c ∈ {i, j} This is important to note since non-sparse penalties such as the ridge penalty can still minimize the negative transfer of knowledge through shrinking ξ 0 . This however comes at the expense of variable selection implied in sparse penalties.
Next, we discuss some asymptotic properties of our regularized bivariate model. In order to investigate the asymptotic properties of our regularized model, we first need to examine the properties of θ ij obtained from minimizing the unpenalized likelihood in (14). Here we note that one advantage of our model is that for both the penalized P (θ ij ) and unpenalized likelihood (θ ij ) we are minimizing over the same set of parameters θ ij since ξ 0 ∈ θ ij . Based on mild regularity conditions for dependent observations, it has been shown that the maximum likelihood estimator of θ ij is r p consistent. Refer to Appendix D, for more details (Basawa, 1980;Basawa et al., 1976;Shi and Choi, 2011). Now for the penalized model, let θ * t ij = {θ * t f ij , σ * t ij } t be the true parameter values corresponding to θ t ij = {θ t f ij , σ t ij } t , and letθ ij be the estimated parameters obtained from minimizing P (θ ij ). Hence ξ * 0 and ξ 0 respectively represent the true and estimated value of ξ 0 . For the penalty function P λ (|ξ 0 |), we assume that the penalty is non-negative; P λ (|ξ 0 |) ≥ 0 and P λ (0) = 0, and that larger coefficients are penalized no less than smaller ones; P λ (|ξ 0 |) ≥ P λ (|ξ 0 |) if |ξ 0 | ≥ |ξ 0 |. These are typical assumptions and are satisfied by the aforementioned penalties (Fan and Li, 2001). Further, we assume that the first and second derivatives of P λ (|ξ 0 |) are continuous at ξ * 0 = 0. We next provide two theorems that establish parameter estimation and selection consistency. The theroems provide similar results as in Fan and Li (2001), but defined within our model specifications and based on dependent observations. Note that the " " notation on a function implies a derivative.
The proof for theorems 2 and 3 is detailed in Appendix D. In the above theorem "lim inf " denoted the infimum of the limit points. It is clear from theorem 2 and 3, that if we choose a proper tuning parameter λ and penalty function P λ (|ξ 0 |) there exists an r p consistent estimator for the penalized likelihood P (θ ij ), which possesses the sparsity propertyξ 0 = 0, i.e. asymptotically performs as well as knowing that ξ 0 = 0 beforehand. This result is also known as an oracle property which provides consistency in variable selection (Donoho and Johnstone, 1995). Here variable selection implies selecting whether functions should be predicted independently or not, a stated in theorem 1.

Application to separable covariance
In this section, we provide a direct approach to applying our regularized pairwise approach to separable modeling. Fortunately, since separable models provide a simplified alternative to nonseparable models, avoiding negative transfer can be readily accomplished. As previously mentioned, the covariance function in a separable model is of the form cov f ij (x, x ) = T ij cov(x, x ), where T ij is the between-output covariance matrix and cov(., .) is a covariance function over inputs x ∈ R D , the same for all outputs. For cov f ij (x, x ) to be a valid covariance function, it is required that T ij = {t c,c } be a positive definite matrix with unit diagonal elements (PDUDE) (Zhou et al., 2011). Therefore, in a bivariate case with two outputs i and j, T ij = 1 t ij t ij 1 , such that t ij measures the correlation between output i and j. It is interesting to note that the covariance between outputs only varies through t ij , this is why outputs in separable functions are instantaneously mixed as they are directly derived by a scaling or a rotation to an output space of dimension D. Inspired by Theorem 1 and Appendix C, one can directly show that t ij = 0 ensures that each output is predicted independently. Therefore, in separable modeling, we only need to adjust the penalty function to P λ (|t ij |), and optimize the penalized likelihood while restricting T ij to be PDUDE. It is interesting to mention here the intrinsic coregionalization (IC) model (Helterbrand and Cressie, 1994) which is a simplified version of the LMC previously mentioned in the introduction. The IC is a separable construction that reduces to independent predictions over each output under an isotopic data case and if outputs are modeled as noise free. Unfortunately, despite its ability to avoid negative transfer, such a model cannot make use of commonalities across outputs.

Combining the predictions
Without loss of generality, we focus on predicting output N through sharing information from the remaining N − 1 outputs as shown in Figure 4. Based on (3), for each submodel in Figure 4, the predictive equation for any new input x 0 ∈ X for output N is expressed as where Our goal is to efficiently combine the predictions from the N − 1 bivariate submodels to form an overall result. To this end we utilize the product of GP experts (PoE) model, used in univariate GP's, however implemented within the specifications of our pairwise model (Ng and Deisenroth, 2014;Deisenroth and Ng, 2015). Here, we aim to combine predictions from N −1 "experts", where each expert is a regularized bivariate GP. The PoE model combines the predictions by the product of all expert predictions. In our pairwise model, the PoE implies that pr(y N (x 0 )|y) = N −1 c=1 pr(y N (x 0 )|y cN ). PoE models are straightforward and theoretically appealing as each expert is weighted by the inverse covariance, therefor experts which are uncertain about their predictions are automatically weighted less than experts that are certain about their predictions. However, a major shortcoming of PoE models is that as N increases, the combined prediction tends to be overconfident. For instance, assume that all functions in I −N are exactly equivalent, then we have that pr(y N (x 0 )|y cN ) = N (M, V) ∀c ∈ I −N for some mean M and variance V. Therefore,pr(y N (x 0 )|y) = N (M = M,V = V/(N − 1)) and as N → ∞ =⇒V → 0. Naturally, in such a case we would wantpr(y N (x 0 )|y) = N (M = M,V = V). To this end, we weight the contributions of each expert with a weight β c = 1/(N − 1) for c ∈ I −N . As a result, given that pr(y N (x 0 )|y cN ) = N (M c , V c ) ∀c ∈ I −N , and following the identity that the product of Gaussian distributions is Gaussian, we have thatpr(y N (x 0 )|y) = N (M,V), wherē This efficient closed form inference for combining the bivariate models is independent of the computational graph, and, consequently, facilitates the ability to scale to arbitrarily large datasets by parallelization, where each bivariate model is efficiently built with a limited of parameters and a small-scale covariance matrix. Note that (17) is similar to the log opinion pool model (Heskes, 1998) and the Generalized product of experts (Cao and Fleet, 2014;Deisenroth and Ng, 2015), therefore, the key feature of the PoE model is still retained as experts that are uncertain about their predictions are weighted less, also, since c∈I −N β c = 1, then it ensures a consistent model that falls back to the prior. Some slight modifications to the traditional PoE have been also proposed such as the Bayesian Committee Machine (BCM) or the robust BCM (Tresp, 2000;Schwaighofer and Tresp, 2003) . The BCM is based on adjusting the variance of the unweighted PoE, by a prior variance pr(y N (x 0 )), while the robust BCM adjusts the BCM with weights similar to those in our model. However, the BCM is based upon assuming a block diagonal covariance where all experts share the same parameters, which hinders its application in our model.

Numerical Case Studies
We conduct case studies to demonstrate the advantageous features of our regularized and distributed multivariate Gaussian convolution process denoted as MGCP-RD. In Section 5.1, we discuss benchmarked methods and the general setting for our numerical case studies. Then Section 5.2 uses simulated functions to demonstrate the performance of the proposed method under four different model settings. Further an illustrative example is provided in Section 5.2.

General settings
In this section, we discuss the settings used to assess the MGCP-RD performance using simulated data. To evaluate the performance of our proposed method, we randomly generate N signals from different model settings, in which the first N − 1 outputs are used as a training set, while the N th output is selected as the testing function. We repeat the study for W = 1000 times. For each replication, we use the mean absolute error (MAE) between the true signal value and its predicted value at p test = 50 points as the criterion to evaluate our prediction accuracy. We then report the distribution of the MAE across the W simulations using a group of boxplots, with respective means represented as black dots.
The MAE values are denoted as prediction errors in all boxplots. Kernel parameters are obtained through minimizing the negative log-likelihood using a scaled conjugate gradient algorithm (Li and Zhou (2016);Álvarez and Lawrence (2011); Rasmussen (2004)). For the MGCP-RD, we distribute the computations over only two systems, where each system was responsible for sequentially fitting (N − 1)/2 of the bivariate models. All computations are done on R-3.2.2 in a 64-bit Windows 7 setting. Further, in our simulation studies, we benchmark our method with four other reference methods for comparison: 1) The individual GP established using a CP, denoted as GCP, where the test function is fitted separately (Rasmussen, 2004); 2) The full MGCP model, denoted as MGCP, described in Section 3.2; 3) The inducing variable approximation, denoted as MGCP-I, which tackles the computational complexity challenge (Álvarez and Lawrence, 2011;Alvarez et al., 2012); 4) The pairwise model for longitudinal profiles, denoted as MGCP-P (Fieuws and Verbeke, 2006;Li and Zhou, 2016;Li et al., 2018). To provide consistent results we utilize the Gaussian kernel in Section 4.1 for each of the benchmarked methods. In the GCP, this kernel reduces to the Gaussian/squared exponential covariance function in the univariate case. For the MGCP-I and since no specific latent structure is proposed we use that of the full MGCP in Section 3.2. However, for the MGCP-P, we utilize the proposed latent structure in the paper which only involves one common latent function for each pairwise model. Finally, throughput the numerical study we use a 3-fold cross validation method to find the tuning parameters for our approach.

Results
We simulate functions from three different settings to demonstrate the benefits of the MGCP-RD. The model settings and results are shown below.

Setting I
In this setting, we aim to compare the performance of the MGCP-RD in a simple case with few number of outputs (N = 5) and no negative transfer. In order to establish a setting with no negative transfer, the multivariate output model for the N curves are generated according to the same functional form y i (x) = 1 + sin(x) + i (x) for x ∈ [0, 10] and i ∈ I. The number of observations per signal is p = p 1 =, .., = p N = 10 evenly spaced points, the test points are evenly spaced across [0, 10] and measurement noise standard deviation is set to σ = σ 1 =, .., = σ N = 0.1 for all outputs. For the MGCP-I, all p design points are used as inducing variables. Also, in this setting we implement the ridge penalization. The importance of this case is that the three challenges of nonseparable modeling can be readily handled in a low dimensional setting with no negative transfer. Therefore, this scenario is able to evaluate the performance of the MGCP-RD relative to the full MGCP. The results are shown in Figure 6. The results in Figure (6) indicate that there is an insignificant difference in the prediction error between the full model (MGCP) and our proposed method. First, amongst the models that considered multiple outputs the MGCP-P had the worst performance. The MGCP-P is based on averaging parameter estimates from paired submodels. Such an approach is extremely dangerous due to the multimodality of the GP/MGP likelihood (Mardia and Watkins, 1989). As mentioned in Rasmussen (2004), different parameters estimates correspond to different interpretations of the data, however, predictions will be moderately effected as the Bayesian correspondence in the GP implies that predictions should pass through or close to the design points. As will be shown in later settings, this approach of averaging parameter estimates becomes specifically dangerous with high noise levels, large parameter space and outputs with varying forms and characteristics where parameter estimates fluctuate widely between different submodels and iterations. Second, and following the intuition why we average predictions rather than parameter estimates, Figure (7) below shows the advantageous features of the weighted PoE model. In Figure (7), we compare the MGCP-RD results with MGCP-RD prediction errors before averaging the N − 1 pairwise submodels from each iteration (denoted as BIVARIATE). Note that since we consider the marginal errors from each submodel in the BIVARIATE then we have W (N −1) errors compared to MGCP-RD with G errors. As shown in the figure, the straightforward mechanism of PoE models, where experts that are uncertain about their predictions are automatically weighted less by the inverse covariance, provides both a simple yet efficient solution for distributed modeling.

Setting II
In this setting, we aim to compare the performance of the MGCP-RD in a case with a moderate number of outputs (N = 50) and no negative transfer. Also, we also aim to illustrate the importance of MGP models in extrapolation which has not been fully exploited in literature. We adopt the quadratic function example from Han et al. (2009) with some modifications. The multivariate output model for the N curves are generated according to y i (x) = 1 + e II x 2 + i (x) for i ∈ I and e II ∼ unif orm(0.8, 1.2). The number of observations for the N − 1 training output is p = 20 evenly spaced points for x ∈ [0, 10], while for the N th function to be predicted we generated p = 10 evenly spaced points for x ∈ [0, 7]. The test points are evenly spaced across [0, 10] and measurement noise standard deviation is set to σ = 1 for all outputs. Due to the long model building time, only W = 50 iterations for the MGCP and MGCP-I are conducted. For the MGCP-I, all p = 20 design points in [0, 10] are used as inducing variables. Also, in this setting we implement the 1 penalization. The importance of this case is that we are able to test all benchmarked models in a rather high dimensional parameter space with higher computational complexity. The results are shown in Figure 8, while an illustrative example of the MGCP-RD results in shown in 9. Based on the results we can obtain some important insights. First, from a computational perspective, and as shown in Figure 8, the model building time of the MGCP and MGCP-I is extremely prohibitive. It takes on average 13 hours to build one MGCP model and this is with a one dimensional input and a moderate number of outputs (N = 50). The mean building time for each iteration was half a second for the GCP, 6.8 minutes for the MGCP-RD and 104 minutes for the MGCP-P. Although the MGCP-P considers paired models however to predict the N th output we need to fit N (N − 1)/2 pairwise submodels, unlike our method where only N − 1 models need to be built. These results illustrate why non-separable MGP models are only used in low dimensional settings. As mentioned previously, for the MGCP-RD, we parallelized computational only over 2 systems (25 models sequentially fit in each system), however, with more computational power, the building time for the MGCP-RD can be significantly reduced.
Second, in addition to the severe computational drawback of the MGCP, its predictive accuracy greatly degrades in such a high dimensional space. This is intuitively understandable, as minimizing the negative likelihood in such a high dimensional (4950) search space is a prohibitive task for any search algorithm. The MLE will directly get trapped in a local minima and will not be able to move even if different starting points are tried, therefore leading to suboptimal parameter estimates with undesirable properties. Besides that, the computation is plagued by numerical issues associated with inverting the 990 × 990 covariance matrix at each iteration of the search algorithm. Similar results have also been shown in Li and Zhou (2016); Li et al. (2018); Shi and Choi (2011), where MGP models tend to loose accuracy in a high dimensional parameter space. This issue is also faced in the MGCP-I, which is not able to address the large parameter space challenge, despite tackling the computational complexity where it only requires the inversion of 20 × 20 covariance matrices. It is important to note that multivariate statistical modeling often encounters functional data with N >> 50 (McFarland et al., 2008;Ramsay, 2006), thus the aforemen-Finally, in Figure 8 we do not report the GCP prediction accuracy to maintain scale, since the mean prediction error is 21.78. The reason for this poor performance is related to the fact the GP in traditional settings cannot extrapolate as correlation goes to zero for data points that are far away from the observed data. However, the MGP is able to borrow strength from other observed outputs to predict the future evolution of a specific output. In other words, extrapolation in the MGP can be seen as interpolation across different output. To highlight this aspect, Figure 10 illustrates the MGCP-RD performance in extrapolation.
In Figure 10, for function 1 we perform interpolation, while for functions 2 and 3 we extrapolate following setting II specifications. However, in function 3, we generate p = 10 evenly spaced points for x ∈ [0, 5] instead of x ∈ [0, 7]. The figure clearly shows that one main advantage of MGP models is the ability to extrapolate an output when other correlated outputs are observed over a larger domain. This advantage of MGP models indeed has many practical applications in cases where extrapolation might be also of interest to the user.

Setting III
In this setting, we aim to compare the performance of the MGCP-RD when negative transfer of knowledge exists. We establish a simple setting with few number of outputs (N = 8). Motivated by M/M/1 queuing systems, the multivariate output for the N functions is generated according to: The number of observations per signal is p = 7 evenly spaced points for x ∈ [0, 0.8], the test points are evenly spaced across [0, 0.8] and measurement noise standard deviation is set to σ = 0.005 for all outputs. In this setting, if we assume x to be the system utilization, then y i (x) respectively define the steady state closed form equations of the expected Queue time and Queue length in an M/M/1 system, where the inter-arrival time is exponentially distributed with rate 2 (Kleinrock, 1976). All p = 7 design points in [0, 0.8] are used as inducing variables for the MGCP-I and we implement the 1 penalization for the MGCP-RD. Following our general settings, for each iteration we find the MAE of the N th function to be predicted which belongs to y i (x) and represents the expected queue length. Also, we benchmark with one other method denoted as MGCP-Sep. In MGCP-Sep, the MGCP is used to predict the N th output using only the training signals with the same functional form, i.e. we fit outputs i ∈ {7, 8} separately using the MGCP. Note that, the importance of our Model Setting III, is that the we are able to test all benchmarked models in a setting where outputs behave according to different functional forms. The results are shown in Figure 11. i (x) are highly correlated, therefore, MGCP-RD was able to learn this cross correlation of y i (x). Indeed, in this model setting we observe that ξ 0 = 0 for pairwise models including y i (x) which indicates that these outputs possess no common features and should be predicted independently. Finally, one important observation is that MGCP-P behaved worse than GCP. As previously mentioned, averaging parameter estimates is specifically dangerous in cases with different functional forms as parameter estimates from different submodels will greatly fluctuate.

Setting IV
In this setting, the goal is to predict the foreign exchange rate compared to the United States dollar currency. The data comes from the pacific exchange rate service (http:// fx.sauder.ubc.ca/data.html). Our analysis utilized the exchange rates of the top ten international currencies (Canadian Dollar CAD/USD, Euro EUR/USD, Japanese Yen JPY/USD, Great British Pound GBP/USD, Swiss Franc CHF/USD, Australian Dollar AUD/USD, Hong Kong Dollar HKD/USD, New Zealand Dollar NZD/USD, South Koreon Won KRW/ USD, Mexican Peso MXN/USD) during the 52 weeks of the 2017 calender year. The data is illustrated in Figure 12. Each output is adjusted to have zero mean and unit variance. We use a leave-one-out cross validation approach to evaluate the performance of the MGCP-RD and the benchmarked methods. We iteratively treat one exchange currency rate as the test output and the remaining 9 currencies as the training set. This procedure is repeated for the 10 currencies. For each test output we randomly remove 13 data points (25% of the data from a specific output) and test the model capability to recover the true underlying values at these test points. For the MGCP, all remaining input points are used as inducing variables for the MGCP-I and we implement the 1 penalization for the MGCP-RD. First we provide some illustrative results in Figure 13. As shown in the figure, the MGCP-RD was able to able to efficiently recover the underlying truth. For performance accuracy comparison we also use the standardized mean square error (SMSE) defined in (Rasmussen, 2004). The results in in Table 1 show that the MGCP-RD was significantly able to outperform the benchmarked methods. This result in intuitively understandable as based on Figure 12, the trends display clear heterogeneity and thus negative transfer is a key issue when integratively modeling the exchange rates.

Conclusion
MGP models established using a CP construction offer a general and flexible solution for multiple output regression. Despite that, the added flexibility arises serious computational challenges even with a moderate number of outputs due to the significant increase in computational demands and number of parameters to be estimated. Further, the integrative analysis of multiple outputs implicitly assumes that these outputs share some commonalities. However, if this does not hold, negative transfer of knowledge may occur. In this paper, we try to simultaneously address the computational (computational complexity, high dimensional parameter space) and negative transfer challenges. To do so, we propose a regularized pairwise modeling approach for MGCP models that has excellent scalability and minimizes the negative transfer of knowledge between uncorrelated outputs. The proposed approach is based on distributing MGCP estimation through bivariate GP submodels which are individually estimated. Predictions are then made through combining predictions from the bivariate models within a Bayesian framework. Interestingly, pairwise modeling turns out to possess unique characteristics which allows to tackle the challenge of negative transfer through penalizing shared latent functions. The modeling framework is generic in terms of the choice of the kernel function and can scale to arbitrarily large datasets by parallelization. We also provide statistical guarantees for the proposed method, extent our method to separable molding cases and demonstrate its advantageous features through numerical studies. The numerical studies illustrate that we can (1) achieve similar prediction performance as the full multivariate approach when the output dimension is low, (2) outperform the full multivariate approach, with only a fraction of its computational needs, when the output dimension is high, (3) outperform the full multivariate approach when some functions are uncorrelated even when the output dimension is low. One important extension of this model lies in the domain of functional graphical models. In such models, nodes are functions rather than random variables. In fact, and since a GP itself is an undirected graphical model, the MGP represents a fully connected undirected functional graphical model where each node represents an output. In our pairwise approach, we are encouraging independence between pairs of functions through our regularization framework. However, an interesting extension would be to extend this regularization framework, to build conditional independence amongst functions in a similar sense to Lemma 1. The main challenge however remains in providing a mapping between the covariance matrix and the precision matrix which control conditional independence between the outputs. We will work along this line and report the results in the future.
Using the inverse variance Lemma, we have that where H 2 = C t y i ,y j C −1 y i ,y i and cov(y j |y i ) = C y j ,y j − C t y i ,y j C −1 y i ,y i C y i ,y j . We now have that Ω ij = −H t 2 cov(y j |y i ) −1 = −C −1 y i ,y i C y i ,y j cov(y j |y i ) −1 . However, since cov(y i , y j ) −1 is positive definite then both C y i ,y i and cov(y j |y i ) are positive definite and thus Ω ij = 0 if and only if C y i ,y j = 0.
To prove Lemma 1, we consider the random vectors y i , y j andÿ whereÿ = {y}/{y c i , y c j }. Now replace y i byÿ and replace y j by the partitioned vector [y t i , y t j ] t . Following (18), we observe that Ω jj = cov(y i , y j |ÿ) −1 . By applying the inverse variance Lemma again, but this time to cov(y i , y j |ÿ) −1 instead of cov(y i , y j ) −1 we have that Ω jj = cov(y i , y j |ÿ) −1 = where H 3 = cov(y i , y i |ÿ) −1 +H t 1 cov(y j |y i ,ÿ) −1 H 1 and H 1 = cov(y i , y j |ÿ) t cov(y i , y i |ÿ) −1 . Then following Lemma 2, the off-diagonal block is zero in this case, if and only if the multivariate Gaussian random vectors y i and y j are conditionally independent where cov(y c i , y c j |ÿ) = 0 for every c ∈ {1, .., p i } and c ∈ {1, .., p j }.

Appendix B.
In this appendix we expand on ∂C ij /∂θ (n) ij from Section 4.1. Recall that θ .
We first expand on ∂C ii /∂θ (n) ij , where, as shown in Section 4.1, the covariance of the marginal process is cov We exclude due to similarity with their counterparts above.
We also exclude due to similarity with their counterparts above.

Appendix C.
In this appendix we prove the following theorem from Section 4.2: Theorem 1. Suppose that ξ 0 = 0, then the predictive distribution of the bivariate model at any new input x 0 ∈ X reduces to that of a univariate model where pr(y i (x 0 )|y ij ) = pr(y i (x 0 )|y i ) and pr(y j (x 0 )|y ij ) = pr(y j (x 0 )|y j ) such that where for c ∈ {i, j}, Ω cc = (C fc,fc + σ 2 c I pc ) −1 ∈ R pc×pc , C fc,f 0 Proof. Based on (9), ξ 0 = 0 implies that, for i = j, cov f ij (x, x ) = ξ 2 0 ∞ −∞ K 0i (u)K 0j (u − d)du = 0 for every d. Therefore, we have that Applying (3) under a bivariate setting we we have that Recall Similarly, we can obtain the proof pr(y j (x 0 )|y ij ).