Regularized Multi-output Gaussian Convolution Process with Domain Adaptation
Wang Xinming Wang Chao Song Xuan Kirby Levi Wu Jianguo






IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

Shell et al.: Regularized Multi-output Gaussian Convolution Process with Domain Adaptation



§ ABSTRACT 
Multi-output Gaussian process (MGP) has been attracting increasing attention as a transfer learning method to model
multiple outputs. Despite its high flexibility and generality, MGP still faces two critical challenges when applied to transfer
learning. The first one is negative transfer, which occurs when there exists no shared information among the outputs.
The second challenge is the input domain inconsistency, which is commonly studied in transfer learning yet not explored in MGP. In this paper, we propose a regularized MGP modeling framework with domain adaptation to overcome these challenges. More
specifically, a sparse covariance matrix of MGP is proposed by using convolution process, where penalization terms are added to
adaptively select the most informative outputs for knowledge transfer. To deal with the domain inconsistency, a domain adaptation
method is proposed by marginalizing inconsistent features and expanding missing features to align the input domains among different
outputs. Statistical properties of the proposed method are provided to guarantee the performance practically and asymptotically. The
proposed framework outperforms state-of-the-art benchmarks in comprehensive simulation studies and one real case study of a
ceramic manufacturing process. The results demonstrate the effectiveness of our method in dealing with both the negative transfer and the domain inconsistency.




Gaussian process, transfer learning, convolution process, domain adaptation.




Regularized Multi-output Gaussian Convolution Process with Domain Adaptation
    Xinming Wang,  Chao Wang,    Xuan Song,    Levi Kirby,    Jianguo Wu  
X. Wang and J. Wu are with the Department of Industrial Engineering and Management, Peking University, Beijing 100089, China.

E-mail: wang-xm20@stu.pku.edu.cn, j.wu@pku.edu.cn.
C. Wang, X. Song and L. Kirby are with the Department of Industrial and Systems Engineering, the University of Iowa, Iowa City 52242, America.

E-mail: {chao-wang-2,  xuan-song,  levi-kirby}@uiowa.edu.
(Corresponding authors: Chao Wang and Jianguo Wu)

    September 9, 2024
=========================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================









§ INTRODUCTION



§ INTRODUCTION




Gaussian process regression (GPR) model has been gaining widespread applications in many fields,  e.g., computer experiments, geostatistics, and robot inverse dynamics <cit.>.
As a powerful nonparametric method, it possesses many desirable and important properties, including excellent fitting capability for various functional relationships under some regularity conditions, providing not only predictions but also uncertainty quantification, and more importantly having closed-form expressions for both tasks.


The conventional GPR models are designed for single-output cases, i.e. the output is a scalar, which have been extensively studied in various applications <cit.>. 
Recently, there has been a growing interest in extending GPR models to multiple outputs, which are ubiquitous nowadays. 
A straightforward way to deal with multiple outputs is known as multi-kriging, which constructs models for each output independently <cit.>. 
It is clear that the multi-kriging impairs the modeling of covariance among outputs, especially when there is strong evidence for the existence of such relationship resulting from physics or constraints. 
Hence, multi-output Gaussian process, which can model correlation among outputs, has been attracting more attention as a joint prediction model. 
The study of MGP begins in geostatistics community known as Co-Kriging in the past few decades <cit.>. 
In today's machine learning society, it is usually known as a multi-task learning method <cit.>, which aims to learn all tasks/outputs simultaneously to achieve better model generalization. 
For example, in the anomaly detection for a manufacturing process, the joint modeling of multiple closely related sensor signals using MGP can help detect the anomaly in each signal more efficiently <cit.>.


However, despite its wide applications, MGP-based multi-task learning requires that the data among these outputs are balanced, which might not be the case in practice. 
For example, in <cit.>,  jointly modeling two correlated pressure signals may not improve the anomaly detection efficiency for both signals if the samples from one signal are much less than the other. 
Combining MGP and transfer learning is an effective way of handling such problems. 
When the observed data for one output is rare or expensive to collect, it is needed to exploit useful information from other outputs whose data are abundant. 
In this work, we focus on making predictions for one output which is denoted as target, by leveraging data in some related outputs which are denoted as sources. 
For instance, in the robot inverse dynamics problem, the target is the torque at a joint when the robot is working with a new load, and the sources are the torques at the same joint when the robot works with other loads <cit.>. 
The key to the MGP-based transfer learning is to extract and represent the underlying similarity among outputs and leverage information from source outputs to target output so as to improve the prediction accuracy <cit.>. 
Specifically, this information transfer in MGP is achieved by constructing a positive semi-definite covariance matrix describing the correlation of data within and across the outputs <cit.>. 
There are two categories of models for the covariance structure: separable models and non-separable models. 
Separable models are most widely used approaches, including intrinsic coregionalization model (ICM) <cit.>, linear model of coregionalization (LMC) <cit.>, and their extensions. 
These models use the Kronecker products of a coregionalization matrix and a covariance matrix of single GP to represent the covariance matrix of MGP. 
It is clear that the separable models are not suitable for transfer learning since it restricts the same covariance structure (from single GP) for both the sources and the target. 
On the other hand, the non-separable models overcome this limitation by using convolution process (CP) to construct the MGP and its covariance structure.
They build non-separable covariance function through a convolution operation and allow modeling each output with individual hyperparameters <cit.>. 
This property makes them more flexible and superior to the separable models <cit.>. 


Nevertheless, there are two critical issues to be considered when apply non-separable MGP to transfer learning. 
The first issue is  negative transfer, which occurs when the assumption of 'existence of shared information' breaks, i.e.,  learning source outputs will have negative impacts on the learning of target output. <cit.>. 
The root cause of this issue is the excessive inclusion of data into the learning process, which is an increasingly severe issue in the big data environment. In such conditions, only a portion of the source data is correlated with the target data  and it is desired to select the sources that yield the best transfer <cit.>. 
A recent work <cit.> proposed a two-stage strategy to alleviate the negative transfer in MGP. 
In this method, the first step is to train a two-output Gaussian process model between each source and the target. 
In the second step, the inverse of predictive standard deviation of each two-output model is adopted as the index to evaluate the negative transfer of each source and integrate the results of transfer learning.
However, the two-stage approach raises significant concerns of losing global information as it only measures the pairwise transferability. 
<cit.> establishes a mixed-effect MGP model which has the ability to infer the behavior of the target output when it is highly similar with some sources. 
However, this method cannot guarantee the optimal selection of related sources, which will be demonstrated in our case study.

The second critical issue is that the input domains of the source processes might be inconsistent with that of the target process. 
For example, in multilingual text categorization, data in different languages have different features and we can't directly combine them to train a classifier for the target data <cit.>. 
Another example is shown in our case study, where the goal is to conduct transfer learning for predicting product density between dry pressing process and additive manufacturing process. 
It is clear that two different manufacturing processes will have different process parameters (inputs) that contribute to the product density, e.g., the dry pressing process is dominated by temperature and pressure while the additive manufacturing process is influenced by solids loading  percent and temperature <cit.>. 
The two processes share one common process parameter (temperature), yet they also have distinct process parameters, which makes the transfer learning of product density (output) a non-trivial task. 
Indeed, the input domain inconsistency is a common issue in transfer learning, and domain adaptation is usually used to overcome this issue. The basic idea of domain adaptation is to align the domains between source and target by transforming data into certain feature domain, and it mainly applies to classification methods such as logistic regression and support vector machine (SVM)  <cit.>. 
These domain adaptation methods aim to find feature mapping by minimizing sum of the training loss of learner and the difference among inconsistent domains, through solving a convex optimization problem. 
However, the training loss of MGP is the negative log-likelihood function, which is a strongly non-convex function. 
Therefore, a unique estimation of the parameters in feature mapping cannot be guaranteed.
More importantly, applying the existing domain adaptation methods directly to MGP might fail the transfer learning due to the existence of negative transfer, i.e., minimizing the difference of features between a negative source and the target will aggravate the severity of negative transfer. To the best of our knowledge, there is no research simultaneously handling issues of domain inconsistency and negative transfer in the context of MGP.


To overcome the above challenges, we propose a comprehensive regularized multi-output Gaussian convolution process (MGCP) modeling framework. In some literature <cit.>, MGP with CP-based covariance is also named as multi-output convolution Gaussian process (MCGP).
Our method focuses on mitigating negative transfer of knowledge while at the same time adapting inconsistent input domain. 
In this work, we assume that there is at least one shared input feature between the sources and the target. 
This assumption is also necessary to facilitate the transferability, i.e., there is nothing to transfer if all the inputs in the sources and the target are different. 
Instead of learning all outputs equivalently, the proposed framework is based on a special CP structure that emphasizes the knowledge transfer from all source outputs to the target output,  which features the unique characteristic of transfer learning and differentiates it from many existing multi-task learning methods. 
The computation complexity is also significantly  reduced from O((qn+n_t)^3) to O(qn^3+n_t^3) when modeling q sources with n data points in each and one target with n_t data points due to this special CP structure. 
The major contributions of this work include:


  * Building upon this special CP structure, a global regularization framework is proposed, which can penalize un-correlated source outputs so that the selection of informative source outputs and transfer learning can be conducted simultaneously.

  * We provide some theoretical guarantees for our method, including the connection between penalizing parameters and selecting source outputs, and the asymptotic properties of the proposed framework.

  * We propose to marginalize extra input features and expand missing input features in the source to align with the input domain of the target, so that the domain inconsistency can be solved. 

Both the simulation studies and real case study demonstrate the effectiveness of our framework in selecting informative sources and transferring positive information even when the target is not quite similar to all the sources.

The remainder of this article is organized as follows. 
The general multi-output Gaussian process and convolution process modeling framework are stated in Section <ref>. 
In Section <ref>, a detailed description of our regularized MGCP modeling framework for transfer learning is presented, including some statistical properties and domain adaptation technique. 
Section <ref> presents numerical studies to show the superiority of the proposed method using both simulated data and real manufacturing data. 
The conclusion is given in Section <ref>. Technical proofs are relegated to the appendix. 




§ PRELIMINARIES






 §.§ Related works on MGCP

As mentioned above, several multi-output Gaussian convolutional process has been investigated for multi-task learning recently.  To handle different kinds of outputs, e.g., continuous output and categorical output, <cit.> proposes a heterogeneous multi-output Gaussian process and conduct variational inference in training and forecasting. Considering that each output may have its unique feature which is not shared with other outputs, <cit.> constructs a MGCP model, where each output consists of two parts: one part is correlated with other output, while the remaining part is independent of others. Compared with these works, our method tackles the problem of inconsistent input domain rather than heterogeneous outputs, and focuses on selecting informative sources in one output (target) prediction. 

Besides for multi-task learning, there are two works using MGCP for information transfer to one output <cit.>. Both of them  pay no attention to the problem of  inconsistent input domain, which limits the available source data for them. In addition, negative transfer is not explored in <cit.>, and the two-strategy method in <cit.> only realizes sub-optimal performance in reducing negative transfer. More detailed comparison can be found in Section 4.5.

Computational load is a severe limitation for multi-output Gaussian process when dealing with large amount of data. In addition to the popular sparse approximation method using inducing variables <cit.>, <cit.> assumes that all q outputs lie in a low-dimensional linear subspace, which can be represented by q̃≪ q orthogonal basis process, and the computational complexity can be reduced to O(q̃ n^3). <cit.> proposes an approach based on local GP experts, which partitions the input and output space into segments to train local experts, then combines them to form a model on full space. These techniques can be applied to our method when extending it to a big data environment. In this paper, we focus more on the effectiveness of our method in reducing negative transfer and handling inconsistent inputs.

Furthermore, convolutional-kernel-based Gaussian process has been applied to high-dimensional and structural data, e.g., image, graph and point cloud data <cit.>. 
In these works, discrete convolution operation is applied on patch of pixels to construct covariance between two data samples. This type of Gaussian process can be applied to image or 3D mesh classification. However, in most of MGCP methods, including our proposed, the convolution operation  is continuous and applied on latent processes.




 §.§ Multi-output Gaussian Process


In this subsection, we will review some basic theories of Multi-output Gaussian process. 
Consider a set of q source outputs f_i : 𝒳↦ℝ, i=1,...,q and one target output f_t : 𝒳↦ℝ, where 𝒳 is an input domain applied to all outputs. 
The q+1 outputs jointly follow some multi-output Gaussian process as

    (f_1,f_2,...,f_q,f_t)^T ∼𝒢𝒫(0, 𝒦(x,x^')),

where the covariance matrix 𝒦(x,x^') is defined as 

    {𝒦(x,x^') }_ij= cov_ij^f (x,x^')= cov( f_i(x), f_j(x^') ),

i,j ∈ℐ= {1,2,...,q,t} and x,x^'∈𝒳. Let ℐ^S={1,2,...,q} denote the index set of source outputs. The element {𝒦(x,x^') }_ij corresponds to the dependency between f_i(x) and f_j(x^'). 

Assume that the observation at point x is

    y_i(x)=f_i(x)+ϵ_i, i ∈ℐ,

where ϵ_i ∼𝒩(0,σ_i^2) is independent and identically distributed (i.i.d) Gaussian noise assigned to the ith output.
Denote the observed data for the ith output as 𝒟_i={X_i, y_i}, where X_i=(x_i,1,...,x_i,n_i), y_i=(y_i,1,...,y_i,n_i)^T are the collections of input points and associated observations,  and n_i is the number of observations for the ith output.  
Suppose that N=∑_i ∈ℐn_i.  
Let 𝒟^S={𝒟_i | i∈ℐ^S} denote the observed data of q source outputs and 𝒟={𝒟^S,𝒟_t} denote all data.
 Define the matrix X and vector y for all input points and observations as
X=(X_1,X_2,...,X_q,X_t ), 
y=(y_1^T,y_2^T,...,y_q^T,y_t^T )^T.

Since GP is a stochastic process wherein any finite number of random variables have a joint Gaussian distribution, for any new input point x_* associated with the target output f_t, the joint distribution of all observations y and the target function value f_t^*=f_t(x_*) is

    [     y; f_t^* ]∼𝒩[ [ 0; 0 ],
    	[           K(X,X)+Σ         K (X, x_*);       K (X, x_*)^T cov_tt^f(x_*, x_*) ] ],

where K(X,X) ∈ℝ^N × N is a block partitioned covariance matrix whose i,jth block, K_i,j∈ℝ^n_i × n_j, represents the covariance matrix between the output i and output j; Σ is a block diagonal noise covariance matrix with Σ_i,j=σ_i^2I_n_i if i=j and 0 otherwise; K(X, x_*)=( K_1,*^T, K_2,*^T,..., K_q,*^T,K_t,*^T )^T and K_i,*=(  cov_i,t^f (x_i,1,x_*), cov_i,t^f (x_i,2,x_*),..., cov_i,t^f (x_i,n_i,x_*) )^T. 
To simplify the notations, we introduce a compact form that K=K(X,X), K_*=K (X, x_*) and C=K+Σ.

Based on the multivariate normal theory, the posterior distribution of f_t(x_*) given data {X,y} can be derived as

    f_t(x_*)| X,y∼𝒩( μ(x_*), V_f(x_*) ),

where the predictive mean μ(x_*) and variance V_f(x_*) can be expressed as

    μ(x_*)   =K_*^T C^-1y, 
    
    V_f(x_*)   = cov_tt^f(x_*, x_*)-K _*^T C^-1K_*.

It can be seen that the mean prediction eq:mean prediction is a linear combination of the observations y , while the variance prediction eq:variance prediction does not depend on y. 
The first term in variance, cov_tt^f(x_*, x_*), is the prior covariance while the second term is the variance reduction due to the mean prediction. 
For the predictive variance of target observation at x_*, we can simply add the noise variance σ_t^2 to that of f_t(x_*).

Equation (<ref>) implies that the key feature of multi-output Gaussian process is borrowing strength from a sample of  q source outputs {f_1,f_2,...f_q} to predict the target output f_t more precisely. 
This effect is achieved by combining the observed source outputs and target output in a linear form wherein the weight is characterized by covariance matrix C and K_*. 
We would like to mention again that the key assumption for the desired function of multi-output Gaussian process is that the source outputs and the target output are correlated, and this correlation can be represented by C and K_*.



 §.§ Convolution Process

From previous studies <cit.>, it is known that the convolution of a Gaussian process and a smoothing kernel is also a Gaussian process. 
Based on this property, we can construct a non-separable generative model which builds valid covariance function for MGP by convolving a base process Z(x) with a kernel g(x). 
More precisely, as shown in fig:convolution process, for output i ∈ℐ, f_i(x) can be expressed as

    f_i(x)=g_i(x)∗ Z (x)=∫_-∞^∞g_i (x-u) Z (u) du,

where ∗ denotes a convolution operation, g_i(x) is the output-dependent kernel function and Z (x) is the shared process across all outputs f_i(x), i ∈ℐ. 

We assume that  Z(x) is a commonly used white Gaussian noise process, i.e., cov(Z(x), Z(x^'))=δ(x-x^') and 𝔼(Z(x))=0, where δ(·) is the Dirac delta function. Note that f_i(x) is also zero-mean GP, thus the cross covariance can be derived as 

    cov_ij^f (x, x^')    = cov{ g_i(x)∗ Z (x), g_j(x^')∗ Z (x^')}
       =∫_-∞^∞ g_i(u)g_j(u-v)d u,

where v=x-x^'. The calculation detail is in Appendix <ref> 
Equation (<ref>) implies that the correlation between f_i(x) and f_j(x^') is dependent on the difference  x-x^' and the hyperparameters in kernels g_i and g_j when they are constructed by a common process. 


 
Specially, if we use the Dirac delta function δ(x) as the smoothing kernel, i.e., g_i(x)=a_i δ(x), the convolution process will degenerate to the LMC model with single shared latent process, i.e. f_i(x)=a_i Z(x) where a_i ∈ℝ is specific to each output i <cit.>. 
So the convolution process can be considered as a dynamic version of LMC because of the smoothing kernel, which also illustrates the superiority of the non-separable MGP model. 

More generally, we can combine the influence of multiple latent processes and extend eq:single convolution process to a more flexible version as

    f_i(x)=∑_e=1^lg_ie(x)∗ Z_e (x),

where l is the number of different latent processes. 
This expression can capture the shared and output-specific information by using a mixture of common and specific latent processes <cit.>. 




§ MODEL DEVELOPMENT


The proposed framework presents a flexible alternative which can simultaneously reduce negative source information transfer and handle inconsistent input domain. 
In Section <ref>, our regularized multi-output Gaussian process model is established using a convolution process under the assumption of consistent input domain. 
The structure of our model enables the separate information sharing between the target and each source. More importantly, our regularized model can realize the selection of informative sources globally. 
Section <ref> provides some statistical properties for the proposed model, including the consistency and sparsity of estimators.
Section <ref> presents the domain adaptation method to deal with the inconsistent input domain of sources. 
In Section <ref>, we discuss the implementation of our model using Gaussian kernel and L_1 norm regularization.



 §.§ Regularized MGCP modeling framework


In this and the following subsection, we focus on the circumstance that the source input domain is consistent with the target input domain. Note that we will relax this assumption in Section <ref>.

As described in Section <ref>, we are provided with q source outputs {f_i | i∈ℐ^S}, one target output f_t, and the observed data 𝒟={𝒟^S,𝒟_t}. 
Under the framework of MGP, we use CP to construct the covariance functions as shown in eq:cov in convolution process. The structure of our model is illustrated in fig:structure. 
With the aim of borrowing information from the source outputs to predict the target output more accurately, the latent process Z_i, i ∈ℐ^S and kernels g_ii,g_it serve as the information-sharing channel between the outputs f_t and f_i. 
On the other hand, Z_i's are set independent of each other so that no information is shared among source outputs, which significantly reduces the computation complexity that will be analyzed later. 
Considering the existence of target-specific behavior, for simplicity yet without loss of generality, a single latent process Z_t(x) is added to the construction of f_t. 



Based on the structure illustrated in fig:structure, the observation of outputs can be expressed as

    y_i(x)   =f_i(x)+ϵ_i(x)=g_ii(x)∗ Z_i(x)+ϵ_i(x), i ∈ℐ^S 
    
    y_t(x)   =f_t(x)+ϵ_t(x)=∑_j ∈ℐg_jt(x) ∗ Z_j(x)+ϵ_t(x),

where g_ii is the kernel connecting latent process Z_i and the output f_i, and g_it is the kernel connecting the latent process Z_i and the target output f_t. For the q source outputs, individual kernel for each source enables an accurate approximation for their feature. 
For the target output f_t, its shared features with source outputs are encoded in Z_i and g_it, i∈ℐ^S, while its specific feature is encoded in Z_t and g_tt.
Based on the assumption that f(x) is independent with ϵ(x), the covariance between any two observations of the outputs i,j ∈ℐ can be decomposed as: 
cov(y_i(x),y_j(x^')) 
= cov(f_i(x),f_j(x^')) +  cov(ϵ_i(x),ϵ_j(x^') ). 
To keep the notational consistency, denote cov(y_i(x),y_j(x^')) as cov_ij^y(x,x^'). As ϵ_i(x), i∈ℐ are i.i.d Gaussian noises, the covariance of two observations y_i(x),y_j(x^') can be expressed as

    cov_ij^y(x,x^')= cov_ij^f(x,x^')+σ_i^2 τ_ij(x-x^'), ∀ i,j ∈ℐ

where τ_ij(x-x^') is equal to 1 if i=j and x=x^', and 0 otherwise. 
Note that every output is a zero-mean GP and {Z_i(x)|,i∈ℐ} are independent white Gaussian noise processes, so cov_ij^f (x,x^')=0 for i,j ∈ℐ^S and i ≠ j, i.e. the covariance across sources are set as zero. And the source-target covariance cov_it^f(x,x^') can be calculated as

    cov_it^f(x,x^')
    =∫_-∞^∞ g_ii(u)g_it(u-v)du,   i ∈ℐ^S

where the last equality is based on eq:cov in convolution process, and v=x-x^'. The detailed calculation can be found in Appendix <ref>. In the same way, we can derive the auto-covariance as

    cov_ii^f(x,x^')   =∫_-∞^∞ g_ii(u)g_ii(u-v)du, i∈ℐ^S 
     cov_tt^f(x,x^')   =∑_j ∈ℐ∫_-∞^∞ g_jj(u)g_jt(u-v)du.

Finally based on the above results, we can obtain the explicit expression of covariance matrix C=K(X,X)+Σ in eq:joint distribution as


    C=
    [ [   C_1,1       0       ⋯       0   C_1,t;       0   C_2,2       ⋯       0   C_2,t;       ⋮       ⋮       ⋱       ⋮       ⋮;       0       0       ⋯   C_q,q   C_q,t; C_1,t^T C_2,t^T       ⋯ C_q,t^T   C_t,t ] ]:=
    [ [   Ω_s,s   Ω_s,t; Ω_s,t^T Ω_t,t , ] ]


where C_i,j=K_i,j+Σ_i,j∈ℝ^n_i × n_j consists of elements {C_i,j}_a,b= cov_ij^y(x_i,a, x_j,b), 1≤ a ≤ n_i, 1 ≤ b ≤ n_j. 
Re-partition the covariance matrix into four blocks: Ω_s,s, a block diagonal matrix, representing the covariance of source outputs' data; Ω_s,t representing the cross covariance between the source outputs and the target output, which realizes the information transfer from sources to the target; Ω_t,t representing the covariance within the target output.

Regarding the structure shown in eq:covariance matrix, there are two interesting points worth of further discussion. 
The first point is setting the covariance across the source outputs to zero, which is the result of the independency among {Z_i}_i=1^q. Ignoring the interactions among sources may cause some  loss of prediction accuracy for the sources, especially when the amount of observed data n_i is small. 
However, we aim to improve the prediction accuracy only for the target and assume the sample data for each source are sufficient, which guarantees the prediction performance of our method. 
Another one is about the covariance between the target and each source, which reveals the advantage of our proposed framework in dealing with negative transfer. 
The source-target covariance function in eq:cov_it illustrates that f_t can share information with each source  through the kernels  g_it and g_ii with different hyperparameters.
It can be intuitively understood that if g_it(x) is equal to 0, the covariance between f_i and f_t is also zero.  
As a result, the prediction of f_t will not be influenced by f_i, i.e., no information transfer between them. 
Further, we derive the following theorem which presents the ability of our model in reducing negative transfer.
theoremTheorem


Suppose that g_it(x)=0, ∀ i ∈𝒰⊆ℐ^S for all x∈𝒳. 
For notational convenience, suppose 𝒰={1,2,...,h|h≤ q}, then the predictive distribution of the model at any new input x_* is unrelated with {f_1,f_2,...,f_h} and is reduced to:

    p(y_t(x_*) | y)=𝒩(   k_+^T C_+^-1y_+, 
        cov_tt^f(x_*,x_*)+σ_t^2-k_+^T C_+^-1k_+),

where k_+=(K_h+1,*^T,...,K_q,*^T,K_t,*^T)^T, y_+=(y_h+1^T,...,y_q^T,y_t^T)^T, and 


    C_+=
    [ C_h+1,h+1         ⋯         0   C_h+1,t;         ⋮         ⋱         ⋮         ⋮;         0         ⋯     C_q,q     C_q,t; C_h+1,t^T         ⋯   C_q,t^T     C_t,t ].



The proof is detailed in Appendix <ref>. 
This theorem demonstrates one key property of our framework. If we penalize the smoothing kernels {g_it(x)}_i ∈𝒰 to zero, the MGCP is actually reduced to a marginalized version, which only contains source outputs {g_it(x)}_i ∈ℐ^S / 𝒰 and the target output. 
The possible negative transfer between {f_i}_i ∈𝒰 and f_t can thus be avoided completely. This result is based on the fact that if g_it(x)=0, 

    cov_it^f (x)=∫_-∞^∞ g_ii(u)g_it(u-v)du=0,

and C_it=0. 

To apply the idea in <ref> to model regularization, we denote that g_it(x)=θ_i0g̃_it(x), where θ_i0 satisfies the condition that g_it(x)=0, ∀x if and only if θ_i0=0. 
Let θ be the collection of all parameters in the model and θ_0={θ_i0| i ∈ℐ^S}⊂θ . 
Then, based on the results of <ref>, our regularized model can be derived as:

    θ max L_ℙ(θ| y)
    =   L(θ| y)-ℙ_γ(θ_0)
    
    =   -1/2y^TC^-1y-1/2 log|C| 
       -N/2log(2π)-ℙ_γ(θ_0),

where L_ℙ(θ| y) denotes the regularized log-likelihood, L(θ| y) denotes the normal log-likelihood for Gaussian distribution, and ℙ_γ(θ_0) is a non-negative penalty function. 
To make the smoothing kernel connecting target and uncorrelated source to 0, common choices of the regularization function include: L_1 norm ℙ_γ(θ_0)=γ |θ_0| and smoothly clipped absolute deviation (SCAD) function <cit.>. 
The validity of our method is ensured by two claims. 
Firstly, based on the theory of multivariate Gaussian distribution, if the source f_i is uncorrelated with the target f_t, then the corresponding covariance matrix block C_it should be zero. 
Secondly, <ref> guarantees that by shrinking some elements of θ_0, {θ_i0}_i∈𝒰, to zero, {C_it}_i∈𝒰=0 and the target output can be predicted without the influence of the source outputs {f_i}_i ∈𝒰. 
Another unique advantage of the proposed method is that it is a global regularized model, since the shrinkage of parameters are applied to all the sources simultaneously, which is different from the local regularization over a subset of data in <cit.>.

Besides the property of global regularization over all the sources, the computational complexity of our method in parameter optimization  is greatly reduced because of the sparse covariance matrix. 
Based on the partitioned covariance matrix
C=
[   Ω_s,s   Ω_s,t; Ω_s,t^T   Ω_t,t ] and using the inversion lemma of a partitioned matrix, the log-likelihood function can be decomposed as:

    L(θ| y)=   -1/2[ỹ^T Ω_s,s^-1ỹ+(Aỹ-y_t)^T B^-1(Aỹ-y_t) ] 
       -1/2[ log |Ω_s,s|+ log |B| ]-N/2log(2π),

where ỹ={y_1^T,...,y_q^T}^T, A=Ω_s,t^T Ω_s,s^-1, B= Ω_t,t-AΩ_s,t is the Schur complement. 
The computational load of MLE is mainly on calculating the inverse of covariance matrix Ω_s,s and B. As Ω_s,s is a diagonal blocked matrix with q square matrix C_i,i∈ℝ^n × n, the complexity for Ω_s,s^-1 is O(qn^3). 
As B∈ℝ^n_t× n_t, the complexity for B^-1 is O(n_t^3). 
As a result, the complexity of our method is O(qn^3+n_t^3). 

However, in the ordinary MGP methods <cit.>, C is a full matrix without zero blocks, so the complexity increases to O((qn)^3). Therefore, the whole computational complexity is O((qn)^3+n_t^3). The above complexity calculation still holds when some sources have different input domains with the target, which will be analyzed in Section <ref>.



 §.§ Statistical properties for regularized MGCP


In Section <ref>, we have discussed that if the source output f_i is uncorrelated with the target output f_t, the cross-covariance between them should be zero, i.e. C_i,t=0.
On the other hand, if the kernel g_it(x)=0, then C_i,t=0 and thus the predictive distribution of f_t(x_*) is uncorrelated with the observations from the source output f_i according to <ref>. 
Therefore, to avoid negative transfer, the estimated parameter θ̂_i0 should be zero, which can be realized through the regularized estimation.

In this subsection, we provide some asymptotic properties of the regularized maximum likelihood estimator θ̂. 
Same as the  last subsection, suppose there are q elements in θ_0, denoted by {θ_10,θ_20,...,θ_q0}, which correspond to the q smoothing kernels {g_1t,g_2t,...,g_qt} respectively. 
Denote the true parameter values of θ_0 and θ in eq:regularized log-likelihood as θ_0^* and θ^*. Suppose there are h zero elements in θ_0^*. 
 Regarding to the penalty function, we assume that ℙ_γ(θ_0) ≥ 0, ∀θ_0; ℙ_γ(0) = 0;  ℙ_λ(θ_0^') ≥ℙ_λ(θ_0) if |θ_0^'| ≥ |θ_0|. 
 These typical assumptions are easily satisfied by the previously mentioned penalty functions.

Before discussing the statistical properties of the regularized model eq:regularized log-likelihood, we first need to introduce the consistency of the maximum log-likelihood estimator (MLE), θ̂_#, for the unpenalized L(θ|y). 
Note that the observations of Gaussian process are dependent. 
So based on some regularity conditions for stochastic process, it has been proved that θ̂_# asymptotically converges to θ^* with rate r_N s.t. r_N →∞ as N→∞, i.e.,

    θ̂_#-θ^* = O_P(r_N^-1).
 
For more details of the regular conditions and consistency proof, please refer to Appendix <ref> and the chapter 7 in <cit.>. 

We first discuss the consistency of the MLE for the regularized log-likelihood L_ℙ(θ| y). 


Suppose that the MLE for L(θ|y), θ̂_#, is r_N consistent, i.e., satisfying non-penalized maximum log-likelihood estimator consistency. If max{ |ℙ^''_γ(θ_i0^*)|: θ_i0^* ≠ 0}→ 0, then there exists a local maximizer θ̂ of  L_ℙ(θ|y) s.t. θ̂-θ^*=O_P(r_N^-1+r_0), where r_0=max{ |ℙ^'_γ(θ_i0^*)|: θ_i0^* ≠ 0}.


The proof is detailed in Appendix <ref>. 
This theorem states that if the derivative of penalty function satisfies some conditions, the estimator of the regularized log-likelihood is also consistent. 
If we take a proper sequence of γ for the penalty, for example choose γ to make r_0 satisfies r_0=o_P(r_N^-1), θ̂ is also r_N consistent as θ̂_#. 
The condition in this theorem, max{ |ℙ^''_γ(θ_i0)|: θ_i0≠ 0}→ 0, is easily satisfied for common regularization functions. 
For example, if ℙ_γ(θ_i0)=γ|θ_i0|, then |ℙ^''_γ(θ_i0)|=0 satisfies. 

Besides the consistency, another key property of θ̂ is that it possess the sparsity, which is provided in <ref> as follows.


Let θ_10^* and θ_20^* contain the zero and non-zero components in  θ_0^* respectively.
Assume the conditions in <ref> also hold, and θ̂ is r_N consistent by choosing proper γ in ℙ_γ(θ_0). If 
N →∞liminf θ→ 0^+liminf γ^-1ℙ^'_γ(θ) >0 and (r_N γ)^-1→ 0, 
then 
    N →∞limP ( θ̂_10=0)=1.


The proof is detailed in Appendix <ref>. 
This theorem implies that by choosing proper penalty functions and tuning parameters, the regularized MGCP model can realize variable selection, i.e. the estimator θ̂ can perform as well as if θ_10=0 is known in advance. More importantly, in our model, the variable selection of θ means the  selection of informative source based on <ref>. 
The conditions in this theorem can also be satisfied easily. Again taking the example ℙ_γ(θ_i0)=γ|θ_i0|, if we let γ=r_N^-1/2, then N →∞liminf θ→ 0^+liminf γ^-1ℙ^'_γ(θ)=1 and (r_Nγ)^-1=r_N^-1/2→ 0, which satisfies the conditions. 



 §.§ Domain adaptation through marginalization and expansion (DAME)


The discussions in Section <ref> and <ref> are based on the assumption that the target and source data share the same input domain. 
However, as we emphasized in the introduction, domain inconsistency is a common issue in transfer learning. 
In this subsection, we propose an effective domain adaptation method for dealing with the domain inconsistency in our MGCP model. 
The general assumption for the proposed domain adaptation method is that there is at least one commonly shared input feature between each  source and the target. 
But we do not require that all sources share the same input, i.e., different sources can share different dimensions with the target.

The basic idea of our domain adaptation method is to first marginalize extra features in the sources, then expand missing features to align with the target input domain. 
More specifically, our method aims to find the marginal distribution of the source data in the shared input domain  with the target,  then create a pseudo dataset in the target input domain. Thus, we name the method as DAME.
This newly created pseudo dataset will be in the same input domain as the target data and have the same marginal distribution as the original source data, which can be used as the new source data to plug in the proposed MGCP model. 
Figure <ref> shows the adaptation procedure using the normalized density data of ceramic product, where the source input domain contains two features x^(c) and x^(s), and the target input domain contains features x^(c) and x^(t). 
In  fig:domain adaptation 1, the source data are marginalized to the domain which only has feature x^(c), and a marginal distribution is obtained based on the marginalized data. 
In  fig:domain adaptation 2, several data are induced according to the marginal distribution, then we expand them to get the pseudo data which have the same features with the target data. 

To generalize the example in fig:domain adaptation, we slightly abuse the notation and focus on one source 𝒟_i={X_i,y_i}, i∈ℐ^S. 
Note that the proposed method will be applied to every source that does not have consistent domain with the target. Let x^(c)∈ℝ^d_c denote the shared features in both the target and source input domain, x^(s)∈ℝ^d_i denote the unique features in the input domain of the ith source, and x^(t)∈ℝ^d_t denote the unique features in the target input domain. 
Then, any source and target data can be expressed as

    x_i,·=[ x_i,·^(c); x_i,·^(s) ], 
    x_t,·=[ x_t,·^(c); x_t,·^(t) ].


The first step is to marginalize the extra features. Define the shared input domain as 𝒳^P which is represented by x^(c), and a projection matrix 
P=[       I_d_c 0_d_c × d_i ]. 
Then, we can get marginalized source data 𝒟_i^P={X_i^P, y_i }, where X_i^P=PX_i =(x^(c)_i,1,...,x^(c)_i,n_i). 
The projected data usually get too dispersed in the shared domain, e.g., the blue triangle in fig:domain adaptation, and the dispersion will be recognized as large measurement noise of the data. 
As a result, a smoothing method is needed to extract the overall trend of the marginalized data and generate induced data with smaller dispersion. 
Many non-parametric methods are available for this purpose, such as kernel regression, B-spline, and GP model, etc. 
In our work, kernel regression is chosen to model the marginalized data, and n_i^' samples {x_i^', a^(c), y_i^', a}_a=1^n_i^' are induced based on the trained model

    y_i^', a=∑_b=1^n_iK_λ(x^(c)_i,b,x_i^', a^(c))y_i,b/∑_b=1^n_iK_λ(x^(c)_i,b,x_i^', a^(c)),

where K_λ is the kernel function and λ is the estimated hyperparameter through cross-validation. 
Note we use i^' to denote a new (marginalized) source resulting from the original source i. 
For example, the mean of marginal distribution is represented by the orange curve in fig:domain adaptation, and the induced data are represented by the orange triangle in fig:domain adaptation 2.



The second step of DAME is to expand the {x_i^', a^(c), y_i^', a}_a=1^n_i^' to include the unique features in the target domain, i.e. x^(t). To realize this idea, we expand the marginalized data along the dimension of x^(t) by adding i.i.d noise which simulates the measurement error. 
For example, if x^(t) is one-dimensional and the observed target data have lower bound x_ low^(t) (e.g., x_ low^(t)=2 in fig:domain adaptation 2) and upper bound x_ up^(t) (e.g., x_ low^(t)=8 in fig:domain adaptation 2), choose n_i^'' values, {x_ i^'',b^(t)}_b=1^n_i^'', spaced in [x_ low^(t), x_ up^(t)]. 
The pseudo data of the source i can be expressed as:

    𝒟_i^ new={[  x_i^', a^(c); x_i^'', b^(t) ], y_i^'', a, b| a∈ [1, n_i^'], b ∈ [1, n_i^''] },

where y_i^'',a, b=y_i^',a+ϵ_a,b and ϵ_a,b is an i.i.d Gaussian noise. For the standard deviation of ϵ_a,b, we can choose the estimated noise deviation of target data using a single-output GP.
In fig:domain adaptation 2, we take n_i^''=4 and the pseudo data is represented by the orange dots.
This expanding approach is intuitively practicable as no prior information about the missing features for source data is given. 
If domain knowledge is available for the expansion step, e.g., there is a monotonically increasing trend along x^(t) in fig:domain adaptation 1, it is also straightforward to incorporate such information.

For the pseudo dataset, our DAME method preserves the marginal information of sources in the shared domain and does not introduce other information in the target-specific domain. This is the unique property of our DAME method. This method will benefit the target if the pseudo dataset is informative. On the other hand, even if the pseudo dataset provides negative information on the target prediction after the marginalization and expansion, our regularized model can identify it and exclude in the learning process. It is worth noting that the existing domain adaptation methods by minimizing the difference between the transformed source and target data might not be able to mitigate the potential negative transfer. Besides, finding an optimal feature mapping (like the existing methods) for the source and target within MGCP training would be extremely difficult, if not impossible. 


In addition, it is worth analyzing the complexity of our method under the circumstance of inconsistent input domains. For the domain adaptation process, the main computational load is on applying the smooth method to obtain the marginal distribution. It will not exceed O(n^3) for each source whether we use kernel regression or GP model. Thus, the complexity of domain adaptation process is not larger than that of constructing the MGCP model in our method. Therefore, the computational complexity of our method is still O(qn^3+n_t^3) when domain inconsistency occurs.




 §.§ Implemetation using Gaussian kernel and L_1 norm


In this subsection, we use Gaussian kernel to implement the modeling framework introduced in Section <ref>-<ref>. 
Gaussian kernel is a very popular choice which is flexible for various spatial characteristics with a small number of hyperparameters. In order to obtain a neat closed form of the convolved covariance function, we take the smoothing kernel as:

    g_ij(x)=α_ijπ^-d/4|Λ|^-1/4exp(-1/2x^T Λ^-1x)  i,j ∈ℐ,

where α_ij is the scaling parameter and Λ is the diagonal matrix representing the length-scale for each input feature.
Using the domain adaptation introduced in Section <ref>, we can assume that every source i ∈ℐ^S is transformed to have the same input domain as the target.
By plugging the kernel eq:smooth kernel in eq:cov_it-eq:cov_tt we obtain


    cov_it^f(x,x^')
    = 2^d/2α_iiα_it|Λ_ii|^1/4 |Λ_it|^1/4/|Λ_ii+Λ_it|^1/2×
                   exp[-1/2(x-x^')^T (Λ_ii+Λ_it)^-1(x-x^') ], 
        cov_ii^f(x,x^')
    = α_ii^2 exp[-1/4(x-x^')^T Λ_ii^-1(x-x^') ], 
        cov_tt^f(x,x^')
    =∑_j ∈ℐα_jt^2 exp[-1/4(x-x^')^T Λ_jt^-1(x-x^')],


where i ∈ℐ^S. eq:cov collection shows that by using the kernel in eq:smooth kernel, the covariance functions are similar to the traditional Gaussian kernel, especially for the auto-covariance within each source.

Then, based on these results,  in regularized log-likelihood of eq:regularized log-likelihood, the collection of all parameters is θ={α_ii, α_it, Λ_ii, Λ_it, σ_i | i∈ℐ}, and the sparsity parameters is θ_0={α_it| i∈ℐ^S}. 
Moreover, if we take L_1 norm as the regularization function and consider the transformed source data, eq:regularized log-likelihood will become


    θ max L_ℙ(θ| y^ new)=   L(θ| y^ new)-γ∑_i=1^q|α_it|,


where y^ new={y^ nda,y^ da}, y^ nda includes the data from sources which do not conduct domain adaptation, y^ da includes the data from sources conducting domain adaptation; γ is the tuning parameter and has critical effect on the optima. 
Typically, the choice of tuning parameter is made through a grid search with cross validation(CV) or generalized CV, such as leave-one-out  (generalized) CV and 5-fold (generalized) CV.

The estimated parameter θ̂ are obtained through solving the problem in eq:regularized log-likelihood with L1.  
However, there are two noticeable details in practice. 
Firstly, as this optimization problem is not a convex problem and multiple local optima exist with high probability, we usually need to set random initial values for several times. 
Secondly, as the commonly used gradient methods, such as L-BFGS method and conjugate gradient method, require the objective function to be smooth, they cannot be applied to solving eq:regularized log-likelihood with L1 because L_1 norm function is not smooth at zero point. 
To solve this issue, we take a Huber smooth approximation as


    γ∑_i=1^q|α_it| ≈γ∑_i=1^q{1/2ηα_it^2, |α_it|≤η
    
    |α_it|-η/2, |α_it|> η.


where η is a small constant, e.g. 10^-4. 
As the maximum bias between the approximation and original function is η, it brings little influence to the optima and makes the common gradient method applicable. 
Finally for prediction, calculate C, K_* and cov_tt^f(x_*,x_*) in eq:mean prediction-eq:variance prediction with θ̂ at point x_*. 
Then, the predictive distribution of f_t(x_*) is in the form of eq:predictive distribution.

The implementation of the regularized MGCP modeling in this work is summarized in Algorithm <ref>.





 §.§ Unique Methodology Contribution

As mentioned in the introduction, the works of  <cit.>,<cit.> also focus on predicting one output through multi-output Gaussian process. Although the covariance structure in <cit.> is similar to us, the work in <cit.> mainly focuses on realizing the transfer from multiple sources to one target. Moreover, works in <cit.> does not consider negative transfer and the corresponding theoretical guarantee on the regularization, which is the major focus in our work.

For the two-stage strategy in <cit.>, it conducts regularization pair-wisely between the target and each source, and combines each sub-model's prediction linearly with different weights. This way has two main drawbacks. 
First, the correlation in one pair might be influenced by other sources. In other words, the strong correlation in one pair might be the results of other sources, and such strong correlation might disappear when considering all sources together. 
Second, the integration of all pairs is conducted by the predictive variance, which is a sub-optimal way for both performance and interpretability.

Finally, compared with existing MGP models, this work is the first one considering inconsistent input domain problem. This technique can increase available source data for transfer learning and multi-task learning with MGP. 





§ NUMERICAL STUDIES


We apply the proposed regularized multi-output Gaussian convolution process model, referred as MGCP-R, to two simulation cases and one real case. In Section <ref>, we introduce the general settings and benchmark methods for our numerical studies. Section <ref> demonstrates the advantages of our method in reducing negative transfer when the sources have the consistent input domain with the target. Section <ref> presents the effectiveness of our framework in dealing with the inconsistent source input domain. And in <ref>, we test and verify the performance with a moderate number of sources and input dimensions. Finally in Section <ref>, we apply the proposed modeling framework to the density prediction of ceramic product.



 §.§ General settings


In this section, we discuss the general settings for assessing the benefits of MGCP-R using simulated data. To evaluate the performance in selecting informative sources and mitigating the negative transfer of knowledge, we randomly generate observations from q source outputs and 1 target output, in which only q_1 source outputs share information with the target output. For simplicity, the q source outputs have equal  number of observations n_1=...=n_q=n, and the target output have less observations, i.e., n_t < n.  These observations form the training set and n_test samples from the target output form the test set. 

For comparison, we take four other reference methods as benchmarks: 


  * The non-regularized MGCP model, whose covariance structure is the same as the proposed model but without regularization, denoted as MGCP; 

  * A regularized MGCP model with a full covariance structure, i.e., constructing the covariance among sources, denoted as MGCP-RF;

  * The two-stage method <cit.> denoted as BGCP-R, which first trains two-output GP models with regularization for each source and the target, then integrate the results of each sub-model in an empirical way;

  * The single GP model constructed by a convolution process, in which only observations from the target output are used for training, denoted as GCP.

In MGCP-RF, the sources and target are modeled as follows:

    y_i(x)   =g_ii(x)∗ Z_i(x)+g_0i(x)∗ Z_0(x)+ϵ_i(x), i ∈ℐ^S 
    
    y_t(x)   =∑_j ∈ℐg_jt(x) ∗ Z_j(x)+g_0t(x) ∗ Z_0(x)+ϵ_t(x),

where Z_0(x) is used for capturing shared information among sources, and Z_i(x) is for unique information in each source/target. This structure refers to <cit.>, but is tailored for transfer learning. To realize similar source-selection effect as MGCP-R, we penalize scale parameters both in g_0i(x) and g_it(x) as a group, i.e., ℙ_γ(θ_0)=γ∑_i=1^q √(α_0i^2+α_it^2). More details can be found in Appendix <ref>.
In BGCP-R, q regularized two-output GP models are trained using the data from each source and the target. The predictive distribution of each sub-model can be expressed as

    f_t(x_*)| X_it, y_it∼𝒩(μ_i (x_*), V_if(x_*) ),

where X_it=(X_i,X_t), y_it=(y_i^T,y_t^T)^T, μ_i (x_*)=K_*^T(X_it,x_*)C(X_it,X_it)^-1y_it, V_if(x_*)= cov_tt^f(x_*,x_*)-K_*^T(X_it,x_*)C(X_it,X_it)^-1K_*(X_it,x_*). Then, the integrated results for BGCP-R is derived as


    f_t(x_*)| X, y∼𝒩( ∑_i=1^qμ_i  (x_*) V_if^-1(x_*)/∑_i=1^q V_if^-1(x_*) , q/∑_i=1^q V_if^-1(x_*) ),


which is an empirical combination of the predictions of each sub-model.
Gaussian kernel in eq:smooth kernel is used for all methods and L_1 norm is used as the regularization function in MGCP-R and BGCP-R.

Regarding the model parameter settings, scaling parameters {α_ii,α_it|i=1,...,q,t} and noise parameter σ for all outputs are initialized with random values in [0,1]. The length-scale diagonal matrix {Λ_ii,Λ_it|i=1,...,q,t} are also initialized with random values in [0,1]. Regarding the hyperparameter learning, we use L-BFGS method in GPflow <cit.>, which is a Python library based on TensorFlow, to maximize the log-likelihood. For the smoothing of L_1 norm regularization function, the value of parameter η in eq:smooth approximation is set to 10^-5. 

Finally, to assess the prediction accuracy, we adopt the mean absolute error (MAE) criterion,

    MAE=1/n_test∑_i=1^n_test| f_t(x_*,i)-f̂_t(x_*,i)|,

where f̂_t(x_*,i) is the predicted mean at x_*,i. We repeat each case for G=100 times and present the distribution of the four methods' MAE in a group of boxplots.



 §.§ Simulation case 1


In order to assess the performance of different methods when the negative transfer of knowledge exists, i.e., learning some sources will bring negative influence on the learning of the target, we adopt an example with one-dimensional input.


The 1D example has q=4 source outputs defined in 𝒳_1=[0,5]:

    f_1(x)   =0.3(x-3)^3,     f_2(x)   =0.3x^2+2sin(2x), 
    
    f_3(x)   =(x-2)^2,     f_4(x)   =(x-1)(x-2)(x-4) ,

and one target output:

    f_t(x)=0.2(x-3)^3+0.15x^2+sin(2x).

The standard deviation of the measurement noise is set as σ=0.2. It can be found that the target output is a linear combination of the outputs f_1 and f_2. The other source outputs, which have different order (f_3) or zero points (f_4), are set as less-correlated sources.  The n=30 observations for each source are evenly spaced in 𝒳_1, and n_t=10 observations for the target are evenly spaced in the left domain, x ∈ [0,3]. The n_test=60 test points are sampled uniformly in 𝒳_1. Note that under such settings, the MAE at these test points contains both the interpolation error and the extrapolation error. 












Considering that the target is a combination of two sources, we benchmark with another method denoted as MGCP-T, which only uses the source outputs f_1 and f_2 to construct a non-regularized MGCP model. MGCP-T is set as the underlying true model and possesses the true covariance structure, wherein negative transfer will not happen. It presents the optimal predictive performance in all introduced methods. 
Figure <ref> shows the boxplots of MAE in these two examples, and fig:case1 1D results shows the data of each source and the predicted trends of the target in one repetition of the 1D example.










Firstly, we focus on three methods, MGCP-R, MGCP-T and MGCP.
The results shown in fig:case1 boxplot illustrate the superior performance of our method. 

MGCP-R performs similarly with MGCP-T and provides much more accurate and stable prediction than MGCP. Note that MGCP-T is the true model with the smallest median and variance value of MAE. This result exactly verifies the conclusion claimed in Section <ref> that our regularized model possesses the ability of selecting informative sources. The negative transfer of information caused by f_3 and f_4 is greatly reduced in the proposed method. 
To state the above conclusion more clearly, we compare part of the estimated parameters of MGCP-T, MGCP-R and MGCP in one repetition of the 1D example. <ref> shows the parameters belonging to the smooth kernels g_it,i∈ℐ^S, which connect the target and each source. As shown in the table, three methods provide similar estimators except the scaling parameters in g_3t and g_4t, which are shrunk to nearly zero in MGCP-R but not in MGCP. We can also directly observe from fig:case1 1D results, a visualization of <ref>, that the predictive mean of the target by MGCP presents an obvious linear shift-up in the right domain, and the sources f_3 and f_4 also shift up linearly at the same area. Moreover, we would like to mention that our method is robust towards both the linear correlation and the non-linear correlation. As the correlation analysis shown in <ref>, in the 1D example, the Pearson correlation between f_4 and f_t is very high while the Kendall correlation between them is low. MGCP-R is not misled by the high linear correlation between f_4 and f_t since it can comprehensively consider both the linear and non-linear relationships and their combinations.
In addition, all source outputs are correlated with each other, which is not listed in <ref>.









An interesting observation is that the MGCP-RF performs only comparable with the MGCP-R in fig:case1 boxplot, although it predicts a little more accurately and also selects the two informative sources in one repetition presented in <ref> and fig:case1 1D results. 

We believe this is due to its larger parameter space. Under same circumstances, MGCP-RF needs 50% more parameters to construct, which poses great challenge in parameter estimation. The following analysis on higher dimension inputs and/or outputs in Sections <ref> and <ref> also confirms our findings. 

The median of prediction error of BGCP-R is the largest. This is because BGCP-R focuses on the  information transfer from each individual source pair-wisely and cannot incorporate the information of all sources globally. As a result, the negative transfer happens to BGCP-R when the target contains combination of sources, i.e., which leads to larger prediction error than only using the target data (GCP). This is one of major shortcomings of BGCP-R since there are very rare cases in practice that the target and source share the same functional form. From the results in fig:case1 1D results, we can observe this influence more clearly: the predictive mean of BGCP-R has a similar valley shape with f_1 in the right domain. 


Finally, the influence of the tuning parameter γ is worth attention, which serves as a similar role of the tuning parameter in LASSO, i.e., there will be a continual selection path as we increase the value of γ. The larger γ is, the less sources will be selected, which means too-large γ may bring negative influence due to the exclusion of some relatively-weak-informative sources. To demonstrate this, we apply MGCP-R to model f_1, f_2 and f_t with varying values of γ. More details and experiment results can be found in Appendix <ref>.



 §.§ Simulation case 2


In this subsection, we apply the proposed framework to transfer information from the source with  inconsistent input domain to the target. We adopt p=3 source outputs:

    f_1(x_1)   =3sin(x_1), 
     
    f_2(x_1,x_2)   =4cos(2x_1)+x_2^2+x_2, 
     
    f_3(x_1,x_2)   =2sin(2x_1)+x_2^2,

and one target output:

    f_t(x_1,x_2)=2sin (x_1)+x_2^2+x_2.

where x_1 ∈𝒳_1 = [-2,2] and x_2 ∈𝒳_2 = [-2,2]. The standard deviation of measurement noise is also set as σ=0.2. In this case, the source f_1 has the inconsistent input domain with the target. Besides, f_1 is set as the obtained mean of marginal distribution after our domain adaptation method. 

According to the notations in Section <ref>, for the source f_1, the common input feature is x^(c)=x_1 and the unique input feature is x^(t)=x_2. Thus, following the procedure of DAME, we firstly generate n_1^'=8  induced data for f_1, {x^(c)_1^',a, y_1^',a}_a=1^8 evenly spaced in 𝒳_1. Then, choose another eight points {x^(t)_1^'',b}_b=1^8 evenly spaced in 𝒳_2. 
The 64 pseudo data of the source f_1 can be obtained through eq:pseudo data, where ϵ_a,b∼𝒩(0, 0.2^2) is the same as the measurement noise of target.  
For the other two sources, n=64 sample points are generated at the same location in 𝒳_1 ×𝒳_2. 
For the target oputput, n_t=24 sample points are located at the nodes of a 3 × 8 grid in [0,2]×[-2,2], and n_test=100 test points are uniformly spaced in 𝒳_1 ×𝒳_2. 
In order to identify the effect of our domain adaptation method, the first source's data are not used in MGCP and BGCP-R.





The results shown in fig:case2 boxplot and fig:case2 results demonstrate the effectiveness of our modeling framework, especially the DAME approach. As the observations of the target are located in the half domain [0,2]×[-2,2], the information of the target's behavior along x_1 can only be borrowed from the source f_1, which leads to the superior performance of the proposed method in the boxplot of MAE. 
Predictive results of one repetition in fig:case2 results also verify the above conclusion, where we can clearly see that MGCP-R is the only method providing accurate fitting in both x_1 and x_2 directions. 
Besides, as shown in the boxplot, the predictive accuracy of MGCP and BGCP-R are better than GCP, because f_2 and f_3 can also provide some beneficial information to the target prediction. However, as more informative information along x_1 is contained in the first source, our method can effectively leverage this knowledge and predict the target more accurately. 




 §.§ Simulation case III


The above simulation cases have demonstrated the effectiveness of our method with small number of sources and input dimensions. In this section, we aim to verify the performance of MGCP-R with more sources and higher input dimensions. 

Setting 1: To test the performance with more sources, we adopt similar setting in the 1D example of Section <ref> and define the following four kinds of functions:

    f_k(x)   =0.3(x-2.5-e_1^k)^3, 
    
    f_n_e+k(x)   =0.3x^2+2sin(2x+e_2^k), 
    
    f_2n_e+k(x)   =(x-1.5-e_3^k)^2, 
    
    f_3n_e+k(x)   =(x-1)(x-2)(x-3.5-e_4^k) ,

where n_e is the number for each kind of sources, e_i^k are uniformly sampled from [0,1], and k={1,...,n_e}. 
We define the target output as:

    f_t(x)=0.2(x-2.5-e_1^1)^3+0.15x^2+sin(2x+e_2^1).

In this setting, we take n_e = 2,4,10, thus the maximum  number of outputs are 41 (including the target). We keep the other settings same as simulation case I and repeat the experiments 50 times.


Setting 2: Regarding the ability of our method with higher input dimensions, we define the following three kinds of sources:

    f_k(x)   =3∑_j=1^2sin(x_j+e_1j^k) , 
    
    f_n_e+k(x)   =4∑_j=1^2cos(2x_j+e_2j^k)+x_3^2+x_3+2x_4-x_5, 
    
    f_2n_e+k(x)   =2∑_j=1^2sin[2(x_1+e_3j^k)]+x_3^2-x_4+2x_5,

where e_ij^k are uniformly sampled from [-0.25,0.25]. 
The target output is:

    f_t(x)=2∑_j=1^2sin(x_j+e_1j^1)+x_3^2+x_3+2x_4-x_5.

In this setting, we take n_e = 1,2, thus the maximum number of outputs are 7. The number of input dimension is 5 and the inconsistent dimensions are 3. 
Similarly to the simulation case II, {f_i(x)}_i=1^n_e is set as the mean of two-dimensional marginal distribution. To apply the  domain adaptation method to these sources, we first generate n_1^'=10 induced data in 𝒳_1={x_1,x_2} from x∼𝒩(0,I_2). Then, we randomly choose 10 points in 𝒳_2={x_3,x_4,x_5} (x∼𝒩(0,I_3)) for each induced data to generate 100 pseudo data. For the other sources, n=100 observations for each source are sampled randomly from x∼𝒩(0,I_5). 
For the target, 150 samples are generated in the same way and 50 of them, which satisfies x_1>0, are picked as training data. Then, another 150 samples are randomly generated as test data. 

The average prediction error shown in <ref> reveals that MGCP-R still performs the best under more sources and higher input dimensions. For a moderate number of sources (40 when n_e=10), MGCP has severe negative transfer effect compared with GCP. As n_e increases,  the difference between MGCP-RF and MGCP-R gets larger, which is expected due to the higher parameter space of MGCP-RF, and thus a large number of hyper-parameters to be optimized.


We also provide the average optimization and prediction time in <ref> for one random start (five random starts in one repetition). The computational load of MGCP-RF is much heavier than other methods, which is a severe drawback. Comparing MGCP-R and MGCP, their prediction time is close but the former's optimization time is less, which is another advantage of regularization. For BGCP-R, it needs less optimization time than MGCP-R in setting 1 due to its smaller parameter space, which makes local optima more easy to reach. However, BGCP-R's optimization time in setting 2 is more than MGCP-R. This is because the difference between their parameter dimensions is smaller than setting 1, and inversing the covariance matrix (O(qn^3+n_t^3) for MGCP-R, O(q(n+n_t)^3) for BGCP) in optimization dominates the computational complexity.







 §.§ Real case of ceramic manufacturing


In this real case study, the goal is to predict the response surface of ceramic product's density. 



  §.§.§ Data description

The data are collected through two groups of experiments differing in manufacturing methods and process parameters. 
The first group contains 28 (4 × 7) experiments using dry pressing manufacturing technique under 4 pressures and 7 temperatures. 
The second group contains 16 (4 × 4)  experiments using stereolithography-based additive manufacturing under 4 solids loadings and 4 temperatures. 
Table <ref> summarizes the values of controlled parameters and all other process parameters are kept fixed in each group. 
As only the temperature is the shared input parameter, the domain adaptation is needed for leveraging information from the data of one manufacturing method to the other.

Two methods, mass-volume method and Archimedes method, are used to measure the density, so there are two sets of measurement for each group. 
The overall 4 datasets are shown in fig:real case data, where the first index of dataset represents the manufacturing method and the second index represents the measurement method. 
Density data of each dataset are standardized to have zero mean and unit variance. 
Note that for the same group of experiments, two measurement methods give different response surfaces because the size of ceramic product is small. 
In this case, the measurement error of Archimedes method might be higher, resulting in negative transfer if we incorporate it in the transfer learning.  
We treat the density data of 1-1 as the target output and the remaining 3 datasets as source outputs. 
For the target, only 8 data points are randomly chosen as observations and the rest 20 points are used for testing. 



For MGCP and BGCP-R, only 1-2, which have the same input domain as the target, is used as the source data in the model. In such condition, MGCP degenerates to a two-output Gaussian Process model, so the main difference between it and BGCP-R is that the regularization in the latter model provides the ability to reduce the negative transfer of knowledge. 
For our method, we apply DAME to the sources 2-1 and 2-2 as follows. Firstly marginalize  the original data to the input domain only with `temperature' feature. Then conduct kernel regression to obtain the mean of marginal distribution.  In this case, we use 7 induced points to keep up with the number of target data. Then expand them to the target input domain and obtain 28 pseudo data. Note that the adaptation process of 2-1 has been shown in fig:domain adaptation before. Thus, we have equal number of data for the original dataset 1-2 and the pseudo dataset of 2-1, 2-2, i.e., n_1-2=n_2-1=n_2-2=28. Finally, they are taken as 3 source outputs to establish a regularized MGCP model, where L_1 norm regularization is implemented. 







  §.§.§ Performance Evaluation

Firstly we provide some intuitive understanding of the advantage of our method. 
From the experimental data, it can be found that temperature is the key factor affecting the density of ceramic products. 
The trend of density in 1-1 is similar to that in 2-1, which makes it feasible to transfer information from 2-1 to 1-1, i.e., from one manufacturing method to another. 
This application is highly desirable in real world as the cost of lab experiment is highly expensive. 
For example, each data point in fig:real case data takes 20 hours to produce. 
Borrowing knowledge from previous research or experiments can greatly reduce the generation of new samples, and acquire a more accurate response surface efficiently and cheaply.



We repeat the case 50 times and the results are shown in Table <ref>. 
The mean and variance of MAE illustrates that MGCP-R outperforms the other benchmarks, with the help of regularization and data from other manufacturing method. 

The results of full-covariance method MGCP-RF are close to MGCP-R, as the optimization problem in large parameter space is not severe for MGCP-RF with a small number of data. Nevertheless, the larger variance of MGCP-RF comparing to MGCP-R and the lower computational complexity of MGCP-R still demonstrate the superiority of our proposed method.
The performance of BGCP-R is almost the same as GCP, while MGCP performs worst in all methods. 
This suggests that the information transferred from the source 1-2 misleads the prediction for the target in MGCP, but BGCP-R reduces its influence to nearly zero through regularization. 
From the predictive results in fig:caseReal results, we can clearly see that MGCP-R is capable of recovering the response surface more accurately with only a few experimental samples, when some historical sources can offer some informative information. 








§ CONCLUSION


We propose a regularized MGCP modeling framework that can select informative source outputs globally and transfer information from sources with both consistent input domain and inconsistent input domains. 
Our work firstly conducts convolution process to establish a special covariance structure that models the similarity within and across outputs. Then, a regularized maximum log-likelihood estimation is performed based on the structure. 
Some statistical properties are also derived to guarantee the effectiveness of our method. 
A domain adaptation approach based on marginalization and expansion deals with the inconsistent input domain of sources successfully. 
Both simulation cases and the real case of ceramic manufacturing demonstrate the superiority of our method.

There are several open topics worthy of investigation in the future based on our work. 
The first one is to apply our method to classification problems, where the posterior distribution needs to be approximated as it doesn't have an explicit form. 
One important issue of classification problems is that the data usually contain considerable amount of features, e.g., gene expression, which increases the complexity and computational burden of GP model.
Therefore, the selection of informative sources and critical features should be combined together, and computationally-efficient algorithms are needed to train the model with high-dimensional data. 
The second one is considering the correlated noise.
For example, in time series analysis, the auto-correlated noise should be considered, which may greatly increase the prediction accuracy and improve the flexibility of the MGCP model. 
Thirdly, in the proposed approach, MGCP modeling and domain adaptation are treated as two separate tasks. Jointly optimizing these two tasks in a unified framework will be studied in future.




§ ACKNOWLEDGMENTS



§ ACKNOWLEDGMENT





10

url@samestyle








Williams2006
C. K. Williams and C. E. Rasmussen, Gaussian processes for machine learning.1em plus 0.5em minus 0.4emMIT press Cambridge, MA, 2006, vol. 2, no. 3.

Shi2011
J. Q. Shi and T. Choi, Gaussian process regression analysis for functional data.1em plus 0.5em minus 0.4emCRC Press, 2011.

Williams1996
C. K. Williams and C. E. Rasmussen, “Gaussian processes for regression,” 1996.

Stegle2008
O. Stegle, S. V. Fallert, D. J. MacKay, and S. Brage, “Gaussian process robust regression for noisy heart rate data,” IEEE Transactions on Biomedical Engineering, vol. 55, no. 9, pp. 2143–2151, 2008.

Boyle2005
P. Boyle and M. Frean, “Dependent gaussian processes,” Advances in Neural Information Processing Systems, vol. 17, pp. 217–224, 2005.

Haas1996
T. C. Haas, “Multivariate geostatistics: an introduction with applications,” Journal of the American Statistical Association, vol. 91, no. 435, pp. 1375–1377, 1996.

Pan2009
S. J. Pan and Q. Yang, “A survey on transfer learning,” IEEE Transactions on Knowledge and Data Engineering, vol. 22, no. 10, pp. 1345–1359, 2009.

Cho2019
W. Cho, Y. Kim, and J. Park, “Hierarchical anomaly detection using a multioutput gaussian process,” IEEE Transactions on Automation Science and Engineering, vol. 17, no. 1, pp. 261–272, 2019.

Chai2008
K. M. Chai, S. Klanke, C. Williams, and S. Vijayakumar, “Multi-task gaussian process learning of robot inverse dynamics,” 2008.

Liu2018
H. Liu, J. Cai, and Y.-S. Ong, “Remarks on multi-output gaussian process regression,” Knowledge-Based Systems, vol. 144, pp. 102–121, 2018.

Conti2010
S. Conti and A. O’Hagan, “Bayesian emulation of complex multi-output and dynamic computer models,” Journal of Statistical Planning and Inference, vol. 140, no. 3, pp. 640–651, 2010.

Goovaerts1997
P. Goovaerts et al., Geostatistics for natural resources evaluation.1em plus 0.5em minus 0.4emOxford University Press on Demand, 1997.

Goulard1992
M. Goulard and M. Voltz, “Linear coregionalization model: tools for estimation and choice of cross-variogram matrix,” Mathematical Geology, vol. 24, no. 3, pp. 269–286, 1992.

Majumdar2007
A. Majumdar and A. E. Gelfand, “Multivariate spatial modeling for geostatistical data using convolved covariance functions,” Mathematical Geology, vol. 39, no. 2, pp. 225–245, 2007.

Alvarez2011
M. A. Alvarez and N. D. Lawrence, “Computationally efficient convolved multiple output gaussian processes,” The Journal of Machine Learning Research, vol. 12, pp. 1459–1500, 2011.

Rosenstein2005
M. T. Rosenstein, Z. Marx, L. P. Kaelbling, and T. G. Dietterich, “To transfer or not to transfer,” in NIPS 2005 workshop on transfer learning, vol. 898, 2005, pp. 1–4.

Weiss2016
K. Weiss, T. M. Khoshgoftaar, and D. Wang, “A survey of transfer learning,” Journal of Big Data, vol. 3, no. 1, pp. 1–40, 2016.

Kontar2020
R. Kontar, G. Raskutti, and S. Zhou, “Minimizing negative transfer of knowledge in multivariate gaussian processes: A scalable and regularized approach,” IEEE Transactions on Pattern Analysis and Machine Intelligence, 2020.

Kontar2018
R. Kontar, S. Zhou, C. Sankavaram, X. Du, and Y. Zhang, “Nonparametric modeling and prognosis of condition monitoring signals using multivariate gaussian convolution processes,” Technometrics, vol. 60, no. 4, pp. 484–496, 2018.

Li2013
W. Li, L. Duan, D. Xu, and I. W. Tsang, “Learning with augmented features for supervised and semi-supervised heterogeneous domain adaptation,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 36, no. 6, pp. 1134–1148, 2013.

He2019
L. He, F. Fei, W. Wang, and X. Song, “Support-free ceramic stereolithography of complex overhanging structures based on an elasto-viscoplastic suspension feedstock,” ACS Applied Materials & Interfaces, vol. 11, no. 20, pp. 18 849–18 857, 2019.

Daume2009
H. Daumé III, “Frustratingly easy domain adaptation,” arXiv preprint arXiv:0907.1815, 2009.

Shi2010
X. Shi, Q. Liu, W. Fan, S. Y. Philip, and R. Zhu, “Transfer learning on heterogenous feature spaces via spectral transformation,” in 2010 IEEE International Conference on Data Mining.1em plus 0.5em minus 0.4emIEEE, 2010, pp. 1049–1054.

Xiao2015
M. Xiao and Y. Guo, “Feature space independent semi-supervised domain adaptation via kernel matching,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 37, no. 1, pp. 54–66, 2015.

Kasarla2021
P. Kasarla, C. Wang, T. L. Brown, and D. McGehee, “Modeling and prediction of driving performance measures based on multi-output convolutional gaussian process,” Accident Analysis & Prevention, vol. 161, p. 106360, 2021.

Moreno2018
P. Moreno-Muñoz, A. Artés, and M. Alvarez, “Heterogeneous multi-output gaussian process prediction,” Advances in neural information processing systems, vol. 31, 2018.

Bruinsma2020
W. Bruinsma, E. Perim, W. Tebbutt, S. Hosking, A. Solin, and R. Turner, “Scalable exact inference in multi-output gaussian processes,” in International Conference on Machine Learning.1em plus 0.5em minus 0.4emPMLR, 2020, pp. 1190–1201.

Yu2021
Z. Yu, M. Zhu, M. Trapp, A. Skryagin, and K. Kersting, “Leveraging probabilistic circuits for nonparametric multi-output regression,” in Uncertainty in Artificial Intelligence.1em plus 0.5em minus 0.4emPMLR, 2021, pp. 2008–2018.

van2017
M. Van der Wilk, C. E. Rasmussen, and J. Hensman, “Convolutional gaussian processes,” Advances in Neural Information Processing Systems, vol. 30, 2017.

Walker2019
I. Walker and B. Glocker, “Graph convolutional gaussian processes,” in International Conference on Machine Learning.1em plus 0.5em minus 0.4emPMLR, 2019, pp. 6495–6504.

Barry1996
R. P. Barry, M. Jay, and V. Hoef, “Blackbox kriging: spatial prediction without specifying variogram models,” Journal of Agricultural, Biological, and Environmental Statistics, pp. 297–322, 1996.

Fricker2013
T. E. Fricker, J. E. Oakley, and N. M. Urban, “Multivariate gaussian process emulators with nonseparable covariance structures,” Technometrics, vol. 55, no. 1, pp. 47–56, 2013.

Fan2001
J. Fan and R. Li, “Variable selection via nonconcave penalized likelihood and its oracle properties,” Journal of the American Statistical Association, vol. 96, no. 456, pp. 1348–1360, 2001.

Basawa2014
I. V. Basawa, Statistical Inferences for Stochasic Processes: Theory and Methods.1em plus 0.5em minus 0.4emElsevier, 2014.

GPflow2017
A. G. d. G. Matthews, M. van der Wilk, T. Nickson, K. Fujii, A. Boukouvalas, P. León-Villagrá, Z. Ghahramani, and J. Hensman, “GPflow: A Gaussian process library using TensorFlow,” Journal of Machine Learning Research, vol. 18, no. 40, pp. 1–6, apr 2017.




[
    < g r a p h i c s >
]Xinming Wang
received the B.S. degree in Mechanical Engineering from Tsinghua University, Beijing, China, in 2020. He is currently working towards the Ph.D. degree in Industrial and System Engineering with Peking University, Beijing, China. His current research interests include data science,  transfer learning, and intelligent manufacturing.

[
    < g r a p h i c s >
]Chao Wang
is an Assistant Professor in the Department of Industrial and Systems Engineering at the University of Iowa. He received his B.S. from the Hefei University of Technology in 2012, and M.S. from the University of Science and Technology of China in 2015, both in Mechanical Engineering, and his M.S. in Statistics and Ph.D. in Industrial and Systems Engineering from the University of Wisconsin-Madison in 2018 and 2019, respectively. His research interests include statistical modeling, analysis, monitoring and control for complex systems. He is member of INFORMS, IISE, and SME.

[
    < g r a p h i c s >
]Xuan Song
is an assistant professor at the department of industrial and systems engineering at the University of Iowa. His research interest is additive manufacturing process development and optimization as well as novel applications of AM technologies in various areas, such as biomedical imaging, tissue engineering, energy harvest, robotics, etc. At UIowa, Dr. Song’s research focuses on the development of next- generation additive manufacturing processes with multi-material, multi-scale or multi-directional capabilities. He obtained his Ph.D. degree in industrial and systems engineering from the University of Southern California in 2016.

[
    < g r a p h i c s >
]Levi Kirby
is a PhD student at the University of Iowa. He obtained his Bachelor’s and Master’s Degree from Western Illinois University in engineering technology. At Iowa, his research focuses on various forms of additive manufacturing, including printing of energetic composites and highly dense ceramics. Throughout his collegiate career, he has been awarded the E. Wayne Kay Scholarship, the Departmental Scholar, Magna Cum Laude, and 3MT finalist.

[
    < g r a p h i c s >
]Jianguo Wu
received the B.S. degree in Mechanical Engineering from Tsinghua University, China in 2009, the M.S. degree in Mechanical Engineering from Purdue University in 2011, and M.S. degree in Statistics in 2014 and Ph.D. degree in Industrial and Systems Engineering in 2015, both from University of Wisconsin-Madison.
Currently, he is an Assistant Professor in the Dept. of Industrial Engineering and Management at Peking University, Beijing, China. He was an Assistant Professor at the Dept. of IMSE at UTEP, TX, USA from 2015 to 2017.

His research interests are mainly in quality control and reliability engineering of intelligent manufacturing and complex systems through engineering-informed machine learning and advanced data analytics. He is a recipient of the STARS
Award from the University of Texas Systems, Overseas Distinguished Young Scholars from China, P&G Faculty Fellowship, BOSS Award from MSEC, and several Best Paper Award/Finalists from INFORMS/IISE Annual Meeting.  He is an Associate Editor of the Journal of Intelligent Manufacturing, and a member of IEEE, INFORMS, IISE, and SME.




§ DERIVATION OF COVARIANCE FUNCTION IN CONVOLUTION PROCESS



For the convolution process:

    f_i(x)=g_i(x)∗ Z (x)=∫_-∞^∞g_i (x-u) Z (u) du,

If Z(x) is a commonly used white Gaussian noise process, i.e., cov(Z(x), Z(x^'))=δ(x-x^') and 𝔼(Z(x))=0, then the cross covariance is derived as:

    cov_ij^f (x, x^') = cov{ g_i(x)∗ Z (x), g_j(x^')∗ Z (x^')}
       =𝔼{∫_-∞^∞g_i (x-u) Z (u) du∫_-∞^∞g_j (x^'-u^') Z (u^') du^'}
       =∫_-∞^∞∫_-∞^∞g_i (u)g_j (u^') 𝔼{ Z (x-u)  Z (x^'-u^') } dudu^'
       =∫_-∞^∞∫_-∞^∞ g_i(u)g_j(u^') δ(x-u-x^'+u^') d udu^'
       =∫_-∞^∞ g_i(u)g_j(u-v)d u,

where v=x-x^' and the last equality is based on the property of Dirac function that ∫ g(u^')δ(u^'-x)du^'=g(x).

For our MGCP structure:

    y_i(x)   =f_i(x)+ϵ_i(x)=g_ii(x)∗ Z_i(x)+ϵ_i(x), i ∈ℐ^S 
    
    y_t(x)   =f_t(x)+ϵ_t(x)=∑_j ∈ℐg_jt(x) ∗ Z_j(x)+ϵ_t(x),

the source-target covariance function can be calculated as:

    cov_it^f(x,x^') 
       = cov (f_i(x),f_t(x^')) 
       = cov{g_ii(x)∗ Z_i(x), ∑_j ∈ℐg_jt(x^') ∗ Z_j(x^') }
       = ∑_j ∈ℐ cov{ g_ii(x)∗ Z_i(x), g_jt(x^') ∗ Z_j(x^') }
       =∫_-∞^∞ g_ii(u)g_it(u-v)du,   i ∈ℐ^S

where the last equality is based on eq:cov in convolution process, and v=x-x^'. In the same way, we can derive the auto-covariance as

    cov_ii^f(x,x^')   =∫_-∞^∞ g_ii(u)g_ii(u-v)du, i∈ℐ^S 
     cov_tt^f(x,x^')   =∑_j ∈ℐ∫_-∞^∞ g_jj(u)g_jt(u-v)du.





§ PROOF OF THEOREM 1


Suppose that g_it(x)=0, ∀ i ∈𝒰⊆ℐ^S for all x∈𝒳. 
For notational convenience, suppose 𝒰={1,2,...,h|h≤ q}, then the predictive distribution of the model at any new input x_* is unrelated with {f_1,f_2,...,f_h} and is reduced to:

    p(y_t(x_*) | y)=𝒩(   k_+^T C_+^-1y_+, 
        cov_tt^f(x_*,x_*)+σ_t^2-k_+^T C_+^-1k_+),

where k_+=(K_h+1,*^T,...,K_q,*^T,K_t,*^T)^T, y_+=(y_h+1^T,...,y_q^T,y_t^T)^T, and 

    C_+=
    [ C_h+1,h+1         ⋯         0   C_h+1,t;         ⋮         ⋱         ⋮         ⋮;         0         ⋯     C_q,q     C_q,t; C_h+1,t^T         ⋯   C_q,t^T     C_t,t ].


Proof. Recall that 

    cov_jt^y(x,x^')   = cov_jt^f(x,x^') 
       =∫_-∞^∞ g_jj(u)g_jt(u-v)du, 
     cov_tt^y(x,x^')   = cov_tt^f(x,x^')+σ_t^2 δ(x-x^') 
       =∑_h ∈ℐ∫_-∞^∞ g_hh(u)g_ht(u-v)du+σ_t^2 δ(x-x^'),

for all j ∈{1,2,...,q}, so g_it(x)=0,i∈{1,2,...,h| h ≤ q} implies that cov_it^y (x,x^')=0 for all i ∈{1,2,...,h} and 

    cov_tt^y(x,x^')=∑_i=h+1^t ∫_-∞^∞ g_ii(u)g_it(u-v)du+σ_t^2 δ(x-x^').

Therefore, we have that C_i,t=0, i ∈{1,2,...,h} and partition covariance matrix 
C=
[ C_-   0;   0 C_+ ], 
where
C_-=
[ C_1,1     0     ⋯     0;     0 C_2,2     ⋯     0;     ⋮     ⋮     ⋱     ⋮;     0     0     ⋯ C_h,h;       ].

The predictive distribution at point x_* is 

    y_t(x_*)∼ N(K_*^T C^-1y,  cov_tt^f (x_*,x_*)+σ_t^2-K_*^T C^-1K_*).

Also, based on that cov_it^y(x,x^')=0 for all i ∈{1,2,...,h}, we have that K_*=(0, k_+^T)^T. Let y_-=(y_1^T,...,y_h^T)^T, then y=(y_-^T,y_+^T)^T. Therefore, 

    K_*^T C^-1y   =(0, k_+^T)
    [ C_-   0;   0 C_+ ]^-1
    (y_-^T,y_+^T)^T  
       =(0, k_+^T)
    [ C_-^-1      0;      0 C_+^-1 ]
    (y_-^T,y_+^T)^T 
       =k_+^T C_+^-1y_+, 
    K_*^T C^-1K_*    =(0, k_+^T)
    [ C_-   0;   0 C_+ ]^-1
    (0, k_+^T)^T  
       =k_+^T C_+^-1k_+.


Note that the auto-covariance matrix of target output f_t, C_tt, is also unrelated with observed data {X_i | i=1,2,...,h} which from source output {f_i| i=1,2,...,h}. As a result, the predictive distribution is totally independent on these outputs. Proof completes.



§ REGULARITY CONDITIONS


In this part, we state the regularity conditions for the consistency theorem of the MLE θ̂_#, which are formulated in <cit.>.

Denote y with total N observations as y^N, and let 

    p_k(θ) = p(y^k| θ) /p(y^k-1| θ)

for each k. Assume p_k(θ) is twice differentiable with respect to θ in a neighborhood of θ^*. Also assume that the support of p(y^N| θ) is independent of θ in the neighborhood. Define ϕ_k(θ) = log p_k(θ), and its first derivative ϕ_k^'(θ), second derivative ϕ_k^''(θ). 

For simplicity and without loss of generality, we only consider the conditions for one-dimensional case. Define ϕ_k^*'=ϕ_k^'(θ^*) and ϕ_k^*''=ϕ_k^''(θ^*). Let ℱ_N be the σ-field generated by y_j, 1 ≤ j ≤ N, and ℱ_0 be the trivial σ-field. Define the random variable i_k^* = var(ϕ_k^*' | ℱ_k-1)=𝔼[(ϕ_k^*')^2 | ℱ_k-1] and I_N^* = ∑_k=1^N i_k^*. Define S_N = ∑_k=1^N ϕ_k^*' and S_N^* = ∑_k=1^N ϕ_k^*''+I_N^*. If the following conditions hold:


  (c1) ϕ_k(θ)  is thrice differentiable in the neighborhood of θ^*. Let ϕ_k^*'''=ϕ_k^'''(θ^*) be the third derivative,

  (c2) Twice differentiation of ∫ p(y^N | θ) d μ^N(y^N) with respect to θ of exists in the neighborhood of θ^*,

  (c3) 𝔼|ϕ_k^*''| < ∞ and 𝔼|ϕ_k^*'' +  (ϕ_k^*')^2| < ∞.

  (c4) There exists a sequence of constants K(N) →∞ as N →∞ such that:
			
			
  (i) K(N)^-1 S_N p→ 0,
			
  (ii) K(N)^-1 S_N^* p→ 0,
			
  (iii) there exists a(θ^*) > 0 such that ∀ϵ > 0,  P[K(N)^-1I_N^* ≥ 2a(θ^*)] ≥ 1-ϵ for all N ≥ N(ϵ),
			
  (iv) K(N)^-1∑_k=1^N 𝔼|ϕ_k^*'''| < M < ∞ for all N,
			

then the MLE θ̂_# is consistent for θ^*. There exists a sequence r_N such that r_N →∞ as N→∞, i.e.,

    θ̂_#-θ^* = O_P(r_N^-1).
 











§ PROOF OF THEOREM 2


Suppose that the MLE for L(θ|y), θ̂_#, is r_N consistent, i.e., satisfying non-penalized maximum log-likelihood estimator consistency. If max{ |ℙ^''_γ(θ_i0^*)|: θ_i0^* ≠ 0}→ 0, then there exists a local maximizer θ̂ of  L_ℙ(θ|y) s.t. θ̂-θ^*=O_P(r_N^-1+r_0), where r_0=max{ |ℙ^'_γ(θ_i0^*)|: θ_i0^* ≠ 0}.

Proof. Recall the assumptions in Section <ref>. For the unpenalized log-likelihood L(θ), the MLE θ̂_# is r_N consistent where r_N is a sequence such that r_N →∞ as N →∞. And we have that L^'(θ^*)=O_P(r_N) and I_N(θ^*)=O_P(r_N^2), which are the standard argument based on the consistency of estimator. Based on that, we aim to study the asymptotic properties of the penalized likelihood L_ℙ(θ)=L(θ)-r_N^2ℙ_γ(θ_0). Here we multiply the penalty function by r_N^2 to avoid that penalty term degenerates as N →∞. The following proof is similar to that of Fan and Li <cit.> but based on dependent observations.

To prove theorem 2, we need to show that for any given ϵ>0, there exists a large constant  U such that:

    P{u=Usup L_ℙ(θ^*+r_N^+u)<L_ℙ(θ^*) }≥ 1-ϵ,

where r_N^+=r_N^-1+r_0. This implies that with probability at least 1-ϵ there exists a local maximum in the ball {θ^*+r_N^+u: u≤ U}. So the local maximizer θ̂ satisfies that θ̂-θ^*=O_P(r_N^+).

By ℙ_γ(0)=0, we have

    L_ℙ(θ^*+r_N^+u)-L_ℙ(θ^*) 
       ≤ L(θ^*+r_N^+u)-L(θ^*) 
          -r_N^2∑_i=h+1^q [ ℙ_γ(|θ_i0^*+r_N^+u_i0|)-ℙ_γ(|θ_i0^*|) ],

where h and q are the number of zero components and all components in θ_i0^*, and u_i0 is the element corresponding to θ_i0 in u. Let I_N(θ^*) be the finite and positive definite information matrix at θ^* with N observations. Applying a Taylor expansion on the likelihood function, we have that 

    L_ℙ (θ^*+r_N^+u)-L_ℙ(θ^*) 
       ≤ r_N^+ L^'(θ^*)^Tu-1/2(r_N^+)^2u^T I_N(θ^*)u[1+o_P(1)] 
          -r_N^2∑_i=h+1^q { r_N^+ ℙ_γ^'(|θ_i0^*|)  sign (θ_i0^*) u_i0
                    +1/2(r_N^+)^2ℙ_γ^''(|θ_i0^*|)u_i0^2[1+o_P(1)] },

Note that L^'(θ^*)=O_P(r_N) and I_N(θ^*)=O_P(r_N^2). so the first term on the right-hand side of t2-proof-Taylor is on the order O_P(r_N^+r_N), while the second term is O_P( (r_N^+r_N)^2 ). By choosing a sufficient large U, the first term can be dominated by the second term uniformly in u=U. Besides, the absolute value of the third term is bounded by

    √(q-h)r_N^2r_N^+r_0 u+(r_Nr_N^+)^2 max{ |ℙ^''_γ(θ_i0)|: θ_i0≠ 0}u^2,

which is also dominated by second term as it is on the order of o_P((r_Nr_N^+)^2 ). Thus, t2-proof-target holds and the proof completes.



§ PROOF OF THEOREM 3


Let θ_10^* and θ_20^* contain the zero and non-zero components in  θ_0^* respectively.
Assume the conditions in <ref> also hold, and θ̂ is r_N consistent by choosing proper γ in ℙ_γ(θ_0). If 
N →∞liminf θ→ 0^+liminf γ^-1ℙ^'_γ(θ) >0 and (r_N γ)^-1→ 0, 
then 
    N →∞limP ( θ̂_10=0)=1.


Proof. To prove this theorem, we only need to prove that for a small ϵ_N=Ur_N, where U is a given constant and i=1,...,s,

    ∂ L_ℙ (θ)/∂θ_i0θ_i0<0, 0<|θ_i0|<ϵ_N.

By Taylor's expansion,

    ∂ L_ℙ (θ)/∂θ_i0 
    =∂ L (θ)/∂θ_i0-r_N^2 ℙ_γ^'(|θ_i0|)  sign (θ_i0) 
       =∂ L (θ^*)/∂θ_i0+[ ∂( ∂ L (θ^*)/∂θ_i0)/ ∂θ]^T(θ-θ^*)[1+o_P(1)]
          -r_N^2ℙ_γ^'(|θ_i0|)  sign (θ_i0).

As ∂ L (θ)/∂θ_i0=O_P(r_N), ∂( ∂ L (θ^*)/∂θ_i0)/ ∂θ_j =O_P(r_N^2) by the standard argument for r_N consistent estimator, thus

    ∂ L_ℙ (θ)/∂θ_i0   =O_P(r_N)-r_N^2ℙ_γ^'(|θ_i0|)  sign (θ_i0) 
       =r_N^2γ( O_P(1/r_Nγ)-γ^-1ℙ_γ^'(|θ_i0|)  sign (θ_i0) ).

Because that N →∞liminf θ→ 0^+liminf γ^-1ℙ^'_γ(θ) >0 and (r_N γ)^-1→ 0, ∂ L_ℙ (θ)/∂θ_i0 will be positive while θ_i0 is negative and vise versa. As a result, t3-proof-target follows. Proof completes.



§ INTERPRETATION OF THE BENCHMARK: MGCP-RF



The illustration of MGCP-RF is shown in fig:structure-MGCP-RF.


In this structure, target f_t is generated by three kinds of latent process: Z_0(x), {Z_i (x)}_i=1^q and Z_t(x). As Z_0(x) is the common process shared by sources, the covariance matrix blocks between source f_i and the other outputs are zero only when the scale parameters in g_0i(x) and g_it(x)  are zero simultaneously. Thus, the marginalized covariance matrix C_+ in Theorem 1 will be:

    C_+=
    [ C_h+1,h+1         ⋯   C_h+1,q   C_h+1,t;         ⋮         ⋱         ⋮         ⋮; C_h+1,q^T         ⋯     C_q,q     C_q,t; C_h+1,t^T         ⋯   C_q,t^T     C_t,t ].

The difference to MGCP-R is that covariance among the remaining sources {f_i}_i=h+1^q can be modeled. This structure is indeed more comprehensive but with the cost of a half more parameters than MGCP-R. The cost will increase if we use more latent process to model the correlation among sources.

To realize the effect of shrinking g_0i(x) and g_it(x) at the same time, group-L1 penalty is used and the penalized log-likelihood function is:

    θ max L_ℙ(θ| y)=   L(θ| y)-γ∑_i=1^q √(α_0i^2+α_it^2),





§ INFLUENCE OF TUNING-PARAMETER



To test the influence of the tuning-parameter γ in our model, we conduct the following experiment. Based on the same dataset in the 1D example of simulation case I, we construct MGCP-R model only with sources f_1 and f_2, and let γ vary from 0 to 10 at a step of 1. Note that MGCP-T is equal to the model with γ=0. The boxplot of MAE with respect to different values of γ is shown in Fig. <ref>. The estimated value of α_1t, α_2t in one repetition is presented in Fig. <ref>. It can be seen that as γ increases, source f_2 will be excluded from the prediction of target, leading to an increased prediction error. In practice, cross-validation can be used to select an optimal tuning-parameter.






