Non-stationary and Sparsely-correlated Multi-output Gaussian Process with Spike-and-Slab Prior
Wang Xinming Li Yongxiang Yue Xiaowei Wu Jianguo






















Wang et al. 






Non-stationary and Sparsely-correlated MGP




Non-stationary and Sparsely-correlated Multi-output Gaussian Process with Spike-and-Slab Prior[The code capsule has been submitted to Code Ocean with provisional DOI: 10.24433/CO.4010696.v1. ]






Xinming Wang
Department of Industrial Engineering and Management, Peking University, China. wang-xm20@stu.pku.edu.cn 
Yongxiang Li
Department of Industrial Engineering and Management, Shanghai Jiaotong University, China
Xiaowei Yue
Department of Industrial Engineering, Tsinghua University, China.
Jianguo Wu
Department of Industrial Engineering and Management, Peking University, China. j.wu@pku.edu.cn 

 


Multi-output Gaussian process (MGP) is commonly used as a transfer learning method to leverage information among multiple outputs. A key advantage of MGP is providing uncertainty quantification for prediction, which is highly important for subsequent decision-making tasks. However, traditional MGP may not be sufficiently flexible to handle multivariate data with dynamic characteristics, particularly when dealing with complex temporal correlations. Additionally, since some outputs may lack correlation, transferring information among them may lead to negative transfer. To address these issues, this study proposes a non-stationary MGP model that can capture both the dynamic and sparse correlation among outputs. Specifically, the covariance functions of MGP are constructed using convolutions of time-varying kernel functions. Then a dynamic spike-and-slab prior is placed on correlation parameters to automatically decide which sources are informative to the target output in the training process. An expectation-maximization (EM) algorithm is proposed for efficient model fitting. Both numerical studies and a real case demonstrate its efficacy in capturing dynamic and sparse correlation structure and mitigating negative transfer for high-dimensional time-series data.  Finally, a mountain-car reinforcement learning case highlights its potential application in decision making problems.







Transfer learning, Gaussian process, non-stationary correlation, negative transfer.

Search for singly charmed dibaryons in baryon-baryon scattering
    Jialun Ping^1
    
School of Physics, Nankai University, Tianjin, 300071, China 

===================================================================




















§ INTRODUCTION


Gaussian process (GP) provides an elegant and flexible Bayesian non-parametric framework for modeling nonlinear mappings <cit.>. Characterized solely by mean and covariance functions, it is capable of capturing complex input-output relationships, as well as measuring prediction uncertainty which is critical for decision-making. As a result, GP has been widely applied in various fields, such as 
Bayesian optimization <cit.>, experiment design  <cit.>, and product quality monitoring <cit.>. 
However, standard GP is designed for only one single output, which limits its use in multi-output or multi-task scenarios arising in various fields, such as Bayesian optimization <cit.>, traffic network <cit.>, and computer simulation emulator <cit.>. Consequently, multi-output Gaussian process (MGP) has been gaining increasing attention from researchers and has emerged as an important member in the vast family of transfer learning <cit.> and multi-task learning <cit.> methods. 

Stationary GP and MGP models are commonly used with covariance function depending only on the distance among data points. However, this invariance to input space translation makes them unsuitable for non-stationary environments, where the data characteristics vary across the input domain <cit.>. 
This phenomenon is quite common in time-series data. For instance, in the energy field, the mean of power consumption in a household is different in every season <cit.>. 
In clinical studies, sepsis is very likely to cause changes in the cross-correlation among vital signs in the early onset <cit.>. 
In kinesiology, the cooperation patterns of human joints vary across different gestures, e.g., both hands move jointly in a `shoot' action but separately in a `throw' action <cit.>. 
In such cases, non-stationary models, that allow all or a subset of parameters to vary are generally more appropriate. Modeling and capturing such a structural change are important in subsequent decision making tasks, such as identifying the risk of disease and taking a medical care <cit.>.

Mainly two kinds of methods have been proposed to capture the dynamic characteristics of the non-stationary data. 
The first category assumes that the parameters are the same within local regions but different across regions. For example, a Bayesian tree-based GP <cit.> uses a tree structure to partition the input space of computer simulation emulator. Another method, called jump GP, cuts the input space into several segments to model piece-wise continuous functions <cit.>. This model is optimized using the expectation-minimization (EM) algorithm or variational inference. Besides, a clustering-based GP <cit.> partitions the spatial data into groups by calculating a cluster dissimilarity and constructs stationary GPs for each group of data. A space-partitioning based GP is further extended to the active learning area to accelerate the design of experiments of heterogeneous systems <cit.>. However, these methods are not suitable for data with gradually-changing characteristics. 
To address this issue, the methods in the second category abandon the locally-stationary assumption. They allow all or some parameters to be input/time-dependent, and model those parameters by additional GPs.  For instance, non-stationary GP introduced in <cit.> and <cit.> applies GP priors on the amplitude and length-scale parameters of a square-exponential covariance function. 
In addition, based on these single-output non-stationary GPs, researchers have explored MGP models for multivariate data with dynamic characteristics. For example, a non-stationary MGP is established to model the varying correlation between vital signals, where a GP prior is imposed on the time-dependent correlation matrix <cit.>. 
However, in this state-of-the-art MGP model, the GP prior has no shrinkage effect, which does not encourage a sparse estimation of cross-correlation among multiple outputs.

Pursuing sparse estimation of cross-correlation is rooted in the negative information transfer, which is another critical challenge when using MGP. Transfer learning is very promising for leveraging information to data-insufficient domain, called target domain, from other correlated domains, called source domain. However, not all the data from the source domain are necessarily correlated with the target domain. If knowledge is transferred from uncorrelated domains, it may reduce the performance of target learning, known as negative transfer <cit.>. For example, the motion signal of a specific human joint may only be correlated with a subset of the other joints. In order to recover the joint's motion information by borrowing information from others, it is necessary to detect which joints share similar moving trajectories with the target joint. Therefore, it is crucial that researchers or engineers can make the best choice on using which sources to transfer information. 

Negative transfer exists widely in transfer learning, often stemming from the excessive inclusion of source data. To handle this issue, one straightforward approach is to measure the relatedness of each source to the target and choose the most related one for information transfer. For example, the method proposed in <cit.> takes Jensen-Shannon (JS) divergence as a criterion and selects the source with the least divergence for knowledge transfer. However, such a method only takes the pairwise transferability into account and ignores the global structure between the target and the sources. And this choice is made independently on the specific model before training, which is far from an optimal decision. An alternative approach is the regularized MGP <cit.>, which jointly models all outputs and selects informative sources during the training process. However, all these approaches assume that the source-target cross-correlation is fixed in the time space, and thus cannot model the dynamic and sparse structure among multiple outputs.


To this end, we propose a non-stationary MGP model to capture the varying characteristics of data and mitigate the negative transfer simultaneously. Specifically, we focus on modeling the dynamic and sparse correlations between the sources and the target. 
In the proposed framework, we first construct a convolution-process-based MGP for transfer learning, whose covariance function parameters are allowed to vary in time space. We then apply a spike-and-slab prior to the parameters that are related to the sparse correlation between the sources and the target. The slab part mainly accounts for smoothly-changing or constant correlation parameters, while the spike part is responsible for shrinking some parameters to zero, thereby removing the corresponding uninformative sources.
To the best of our knowledge, this is the first research on MGP that simultaneously handles dynamic relationship and negative transfer.
Our contributions can be summarized as follows:


  * A novel non-stationary MGP model is established using the convolution of latent processes and time-dependent kernel functions, which is suitable for modeling multiple outputs with varying characteristics.

  * A dynamic spike-and-slab prior is applied to capture the temporal and sparse correlations among outputs, deciding from which sources to transfer information to the target.

  * The mixture weight of the spike and slab priors is automatically adjusted during the training process using an EM-based optimization algorithm, which can effectively prevent placing shrinkage effects on non-zero elements. 


The rest of this paper is organized as follows. In Section <ref>, we revisit the related literature and the static MGP model. Section <ref> presents the proposed non-stationary MGP and an efficient EM algorithm for model training.  In Section <ref>, we evaluate the effectiveness of our model on simulated data. In Section <ref>, we perform one time-series analysis case on human gesture data <cit.> and one control policy optimization case on the mountain-car problem <cit.>. In Section <ref>, we conclude the paper with a discussion.



§ PRELIMINARIES


In this section,  we first review researches that are related to our work. We then introduce the static MGP based on the convolution process, which has been widely applied in various areas due to its flexibility <cit.>. 



 §.§ Related work

To deal with non-stationary data, a natural extension of GPs is to release the restriction that the parameters of the covariance functions are invariant throughout the input space. Most of the existing approaches either encourage the parameters to be constant in a local area and construct a piece-wise model (the first category, e.g., <cit.>), 
or allow them to vary at each point and model them using other GPs (the second category, e.g.,  <cit.>). The methods in the second category have a similar structure to that of a two-layer Deep Gaussian Process (Deep GP), where the input is first transformed by the first GP layer into a latent input, and then fed into the second GP layer to obtain the output <cit.>. However, the parameters of Deep GP are stationary, which differs from the second category where the covariance parameters are dynamic.


Non-stationary GPs mainly focus on modeling the dynamic mean, smoothness, and amplitude parameters. With regards to MGP, dynamic correlation is another key characteristic that needs to be considered.  In classical Linear model of Coregionalization (LMC), each output is a linear combination of several latent Gaussian processes, and the covariance matrix is modeled by a Kronecker product of a correlation matrix and a single GP’s covariance matrix <cit.>.The existing non-stationary MGPs are mainly extensions of the classical LMC model. 
For example, the approach in <cit.> allows the correlation matrix to vary with inputs to model the dynamic relationship among outputs. 
In <cit.>, a non-stationary MGP combines the time-varying correlation (across outputs) and smoothness (within each output) together. However, as extensions of the traditional LMC, these methods also suffer from the limitation that all outputs possess the same covariance structure. More flexible MGPs are proposed by constructing each output through the convolution process and modeling them with individual parameters <cit.>. However, these approaches are for stationary data. Furthermore, all existing approaches fail to capture a sparse correlation structure in a non-stationary environment.



 Spatial-temporal modeling of non-stationary data is closely related to our work. In comparison with normal time-series modeling, spatial-temporal analysis requires to model the spatial correlation to enhance the prediction accuracy. A large number of spatial-temporal models have been investigated, such as spatial-temporal auto-regressive integrated moving average method (ST-ARIMA) <cit.>, spatial-temporal k-nearest neighbors (ST-KNN) <cit.>, spatial-temporal random fields <cit.>, and spatial-temporal deep neural networks <cit.>. Based on the aforementioned methods, a number of recent works try to extend them to handle non-stationary spatial-temporal data. One popular and efficient solution is utilizing some change detection algorithm to partition the time-series into several stationary periods, and then applying the stationary model for each period, e.g., a ST-KNN with a wrapped K-means partition algorithm <cit.>, an auto-regressive model coupled with a block-fused-Lasso change-point detection algorithm <cit.>. Besides partitioning the time-series into stationary parts, the method proposed by <cit.> maps the non-stationary space-time process into a high-dimensional stationary process through augmenting new dimensions. Another type of non-stationary spatial-temporal model is Bayesian random fields with non-stationary space-time kernels <cit.>, whose hyper-parameters change over time or location. Deep learning methods are also explored on non-stationary data recently, such as non-stationary recurrent neural networks <cit.>, long short-term memory networks <cit.>, and transformer-based networks <cit.>. In contrast to the spatial-temporal model, MGP does not impose a restriction that the source outputs must be sampled during the same period as the target outputs. Furthermore, it does not depend on spatial distance to establish correlations among outputs. As a result, the MGP model is capable of accommodating a wider range of scenarios. 


It is important to mention that we use a dynamic spike-and-slab prior in our model. The classical spike-and-slab prior is a Bayesian selection approach <cit.> that has been used for feature selection in (generalized) linear models, additive models, and Gaussian processes <cit.>. With this prior, smaller parameters tend to be more influenced by the spike prior to reach zero, while the larger ones are mainly dominated by the slab part and bear little shrinkage. However, the classical spike-and-slab prior cannot account for the modeling of dynamic and sparse correlation parameters in our model. Therefore, we propose to extend this prior to a dynamic version. Although current works have explored the dynamic variable selection for the varying-coefficient linear models <cit.>, no work utilizes a dynamic spike-and-slab prior to model the dynamic and sparse correlation among outputs in a non-stationary MGP.







 §.§ Static MGP based on convolution process


Consider a set of m outputs f_i : 𝒳↦ℝ, i=1,...,m, where 𝒳 is a d-dimensional input domain applied to all outputs. Suppose that the observation y_i is accompanied with independent and identically distributed (i.i.d.) noise ϵ_i ∼𝒩(0, σ_i^2), i.e.,

    y_i (x) = f_i(x) + ϵ_i.

where x∈ℝ^d is the input.
Denote the n_i observed data for the ith output as 𝒟_i={X_i, y_i}, where X_i=(x_i,1,...,x_i,n_i)^T and y_i=(y_i,1,...,y_i,n_i)^T are the collections of input points and associated observations respectively. Let the total number of observations be represented by N = ∑_i n_i for m outputs. Denote the data of all outputs as 𝒟={X, y}, where X=(X_1^T,...,X_m^T)^T ∈ℝ^N × d and y=(y_1^T,...,y_m^T )^T ∈ℝ^N. 
In an MGP model, the observation vector y follows a joint Gaussian distribution:

    y| X∼𝒩[ 0, K ],

where K=K(X,X) ∈ℝ^N × N is a block-partitioned covariance matrix. 














The (i,i^')-th block of K, K_i,i^'= cov_i,i^'^f(X_i,X_i^') + τ_i, i^'σ_i^2 I∈ℝ^n_i × n_i^' represents the covariance matrix between the output i and output i^' (τ_i, i^' equals to 1 if i = i^', and 0 otherwise). The function cov_i,i^'^f(x, x^') measures the covariance between f_i(x) and f_i^'(x^'). In the covariance matrix K, the cross-covariance block K_i,i^'(i ≠ i^') is the most important part to realize information transfer, as it models the correlation between different outputs.

As the convolution of a GP and a smoothing kernel is still a GP, we can construct each output f_i through convolving a group of shared latent processes {z_j(x)}_j=1^h and kernel functions {g_ji(x)}_j=1^h  in the following way <cit.>:

    f_i(x) = ∑_j = 1^h α_jig_ji(x)∗ z_j(x) =  ∑_j = 1^h α_ji∫_-∞^∞ g_ji(x-u) z_j(u) d u

where ∗ represents convolution operation, α_ji is the amplitude parameter, and h is the number of shared latent processes. Usually, {z_j(x)}_j=1^h are independent white Gaussian noise processes with cov(z(x), z(x^')) = δ(x-x^'), where δ(·) is the Dirac delta function. Thus, the covariance function can be derived as:

    cov_i,i^'^f (x, x^')  =  cov[f_i(x),  f_i^'(x^')]  
       = ∑_j=1^h cov{α_ji g_ji(x)∗ z_j (x), α_ji^' g_ji^'(x^')∗ z_j (x^')}
       =∑_j=1^h  α_jiα_ji^'∫_-∞^∞ g_ji(u)g_ji^'(u-v)d u,

where v=x-x^'. 
It such a way, the covariance between f_i(x) and f_i^'(x^') is dependent on their difference x-x^', the amplitude parameters, and the hyperparameters in kernels g_ji and g_ji^'. Compared with the classical LMC model f_i(x) = ∑_j = 1^h α_ji q_j(x),
where q_j(x) is a latent GP with covariance k_j(x, x^'). The convolution-process-based MGP is more flexible than LMC, as it does not restrict all outputs to having the same auto-covariance pattern.

At a new point x_*, the posterior distribution of y_i(x_*) given data {X,y} is:

    y_i(x_*)| X,y∼𝒩( μ(x_*), Σ(x_*) ),

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

    μ(x_*)   =K_*^T K^-1y, 
    Σ(x_*)   = cov_ii^f(x_*, x_*)  + σ_i^2 - K _*^T K^-1K_*,

where K_*^T = [ cov_i1^f(x_*, X_1)^T, ...,  cov_im^f(x_*, X_m)^T] is the covariance between the new point x_* and all observed data. From the posterior distribution, we can find that the covariance function plays a crucial role in prediction. For instance, the predicted mean is the linear combination of output data, where the weight is decided by the covariance matrix. However, in the static MGP, the covariance between two data points depends solely on their distance and does not change dynamically. Additionally, some outputs may be uncorrelated with others, therefore the estimated covariance matrix should possess a sparse structure to avoid negative transfer between the uncorrelated outputs. In the following section, we will propose a novel non-stationary MGP to simultaneously address both problems.
















































§ MODEL DEVELOPMENT


We propose a non-stationary MGP model for transfer learning that can capture sparse source-target correlation in a dynamic environment. Specifically, we assume that the correlations between the target and each source vary over time. Besides, some sources may not be related to the target during certain time periods. 

Under such a circumstance, a spike-and-slab prior is utilized to model the varying and sparse correlation structure. 



 §.§ The proposed model.

The structure of our hierarchical model is illustrated in fig: graphical structure. The first layer constructs outputs through the convolution of time-dependent kernel functions and latent white Gaussian noise processes, and the second layer consists of priors on function parameters designed to encourage desired properties, such as smoothness and sparsity.





  §.§.§ Dynamic MGP

In this subsection, we introduce the major part of the proposed non-stationary MGP, which corresponds to the first layer in fig: graphical structure.
To simplify the notation, we slightly abuse the notation used in the previous section. Specifically, we take the first m-1 outputs f_i : 𝒳↦ℝ, i=1,...,m-1 as the sources, and the last one f_m : 𝒳↦ℝ as the target. We still assume that the observation y_i is accompanied with the i.i.d. measurement noise ϵ_i ∼ N(0, σ_i^2). Let ℐ={1,2,...,m} be the index set of all outputs, and  ℐ^S=ℐ/{m} contain the indices of all sources. For simplicity yet without loss of generality, we assume the source and target data are sampled at the time t ∈{1, 2, ... ,n}, i.e., X_i = (x_i, 1, ..., x_i, t, ..., x_i, n)^T and y_i = (y_i,1,...,y_i,t,..,y_i,n)^T. In Appendix A, we will show a more general case where each output is observed only at a subset of time stamps. 


Based on the above assumption, our dynamic MGP model is formulated as:

    y_i (x_t)     = f_i(x_t) + ϵ_i =α_ii, t g_ii, t(x_t)∗ z_i(x_t) + ϵ_i , i ∈ℐ^S 
    
    y_m (x_t)     = f_m(x_t) + ϵ_m =∑_j ∈ℐα_jm, tg_jm, t(x_t) ∗ z_j(x_t) + ϵ_m






where α_ii, t and α_jm,t are time-varying amplitude parameters, g_ii, t(x) and g_jm, t(x) are time-varying kernel functions, and {z_i(x)}_i=1^m are latent white Gaussian noise processes independent of each other. 
This model is highly flexible as various types of kernel functions can be utilized. We choose to employ a Gaussian kernel which is widely used due to its flexibility <cit.>. The kernel is given by 

    g_ij, t(x)= (2π)^-d/4|θ_ij, t|^-1/4exp(-1/2x^Tθ_ij, t^-1x),




where θ_ij, t is a diagonal matrix representing the length-scale for each input feature.  More importantly, such a Gaussian kernel can yield closed-from covariance functions through the convolution operation in eq:cov in convolution process <cit.>. 

This flexible model allows each source to have its own kernel g_ii, thereby allowing for heterogeneity among the sources. In order to transfer knowledge from the sources to the target, the target is connected to {z_i}_i=1^m-1 though the kernel function g_im, t.  Regarding the parameters, {α_ii, t, θ_ii, t}_i=1^m are responsible for the non-stationary behavior within each output, while {α_im, t, θ_im, t}_i=1^m-1 capture the dynamic correlation between target and sources. More specifically,  the amplitude parameter α_im, t controls the knowledge transfer. For example, if α_im, t=0, then  f_i will not transfer information to f_m at time t.

As the latent processes are independent of each other, the covariance matrix among sources is zero-valued. Therefore, the covariance matrix can be re-partitioned as:

    K=
    [ [     K_1,1         ⋯         0     K_1,m;         ⋮         ⋱         ⋮         ⋮;         0         ⋯ K_m-1,m-1   K_m-1,m;   K_1,m^T         ⋯ K_m-1,m^T     K_m,m ] ]
    =
    [ [   K_(ss)   K_(sm); K_(sm)^T     K_mm ] ],

where the (i,i^')-th block K_i,i^'= cov_i,i^'^f(X_i,X_i^') + τ_i, i^'σ_i^2 I, the block-diagonal matrix K_(ss) represents the covariance of source outputs, and K_(sm) represents the cross-covariance between the sources and the target. 
Based on eq:cov in convolution process,





we can obtain covariance functions for the proposed non-stationary MGP model, as shown below:














    cov_ii^f(x_t, x_t^') 
    
    =   α_ii, tα_ii, t^'|θ_ii, t|^1/4 |θ_ii, t^'|^1/4/|θ_ii, t + θ_ii, t^'|^1/2exp[ -1/2(x_t - x_t^')^T(θ_ii, t + θ_ii, t^')^-1(x_t - x_t^')], 
    
        cov_im^f(x_t, x_t^') 
    
    =   α_ii, tα_im, t^'|θ_ii, t|^1/4 |θ_im, t^'|^1/4/|θ_ii, t + θ_im, t^'|^1/2exp[ -1/2(x_t - x_t^')^T(θ_ii, t + θ_im, t^')^-1(x_t - x_t^')], 
    
        cov_mm^f(x_t, x_t^') 
    
    =   ∑_j ∈ℐα_jm, tα_jm, t^'|θ_jm, t|^1/4 |θ_jm, t^'|^1/4/|θ_jm, t + θ_jm, t^'|^1/2exp[ -1/2(x_t - x_t^')^T(θ_jm, t + θ_jm, t^')^-1(x_t - x_t^')],
























where i ∈ℐ^S.
Equations (<ref>- <ref>) represent the covariance within the sources, between the sources and the target, and within the target, respectively. To ensure the positivity of those hyper-parameters, we utilize a soft-plus transformation for them <cit.>:
α_ij,t = log[1 + exp(α̃_ij,t)],  θ_ij,t = log[1 + exp(θ̃_ij,t)],
where α̃_ij,t, θ̃_ij,t are underlying parameters to estimate, whose range is [-∞, ∞]. 
 The proposed covariance functions can be viewed as an extension of the non-stationary kernels <cit.> from the single-output to the multi-output case. Specifically, the auto-covariance of each source in eq: dynamic source auto-cov is the same as the covariance for a single-output non-stationary GP. From the cross-covariance between each source and the target, we can clearly see that the amplitude parameter α_im, t controls whether the cross-correlation is zero or not. The validity of the proposed covariance functions is outlined in Proposition <ref>:

The proposed non-stationary MGP covariance matrix in eq: covariance matrix is positive-definite, i.e., ∀y≠0,

    y^T Ky > 0.



The proof is provided in Appendix B. 

Based on eq: covariance matrix, the  joint distribution of all sources and the target is expressed as:

    [ [ y_(s);   y_m ] ]| X∼𝒩(
    [ [ 0; 0 ] ],
    [ [   K_(ss)   K_(sm); K_(sm)^T    K_ mm ] ]),

where y_(s) is the collection of all source data. For notational convenience, we partition the parameters in a similar way, 







α_(s) = {α_(s), t}_t=1^n,
θ_(s) = {θ_(s), t}_t=1^n,
σ_(s) = {σ_i}_i=1^m-1,
α_m = {α_m, t}_t=1^n ,
θ_m = {θ_m, t}_t=1^n,
where α_(s), t={α_ii, t}_i=1^m-1, θ_(s), t={θ_ii, t}_i=1^m-1, α_m, t={α_im, t}_i=1^m,  and θ_m, t = {θ_im, t}_i=1^m. Furthermore, we denote the collection of all parameters as Φ = {Φ_(s), Φ_m}, where Φ_(s) = {α_(s), θ_(s), σ_(s)} and Φ_m = {α_m, θ_m, σ_m}. 

















In the proposed model, the most important and challenging task is to estimate those time-varying parameters. If no restriction is applied, model training may suffer from a serious over-fitting problem. To address this issue, Gaussian processes are typically employed to model the kernel parameters, e.g., 
log(α_t) ∼𝒢𝒫(0, k_α(t, t^')), log(θ_t) ∼𝒢𝒫(0, k_θ(t, t^')), where k_α,  k_θ are covariance functions for the amplitude and length-scale parameters respectively. Although this technique can force the parameters to vary smoothly and reduce over-fitting, it cannot model a sparse correlation between the sources and the target. Consequently, this approach cannot avoid the negative transfer caused by unrelated sources. 













  §.§.§ Dynamic spike-and-slab prior.

The classical spike-and-slab prior <cit.> only handles the shrinkage of one single parameter and cannot model the smooth-varying parameters.
To take both the dynamic and sparse property of correlation into account, we propose a dynamic spike-and-slab prior placing on α_m: 

    α_im, t|γ_i, t, α_im, t-1   ∼(1-γ_i, t) p_spike(α_im, t)+γ_i, tp_slab(α_im, t| α_im, t-1) 
    γ_i, t|η   ∼  Bern(η),

where γ_i, t∈{0,1} is a binary sparse indicator for α_im, t following a Bernoulli distribution, p_spike(α_im, t) is a zero-mean spike prior pushing parameter to zero, p_slab(α_t, im| α_t-1, im) is a slab prior connecting α_im, t-1 and α_im, t, and η is a prior weight between the spike and slab. If there is no prior information regarding the weight, we can set η to 0.5. The spike-and-slab prior is shown in the second layer of the graphical structure in fig: graphical structure. As for all the other parameters, we do not force them to possess sparsity, so only the slab prior is placed on them to control the smoothness, i.e.,

    α_ii, t| α_ii, t-1   ∼p_slab(α_ii, t| α_ii, t-1),  
    θ_ii, t| θ_ii, t-1   ∼p_slab(θ_ii, t| θ_ii, t-1),  
    θ_im, t| θ_im, t-1   ∼p_slab(θ_im, t| θ_im, t-1).

Note that θ_ij, t is a diagonal matrix, so that the slab prior is placed on its d diagonal elements independently, i.e., p_slab(θ_ii, t| θ_ii, t-1) = ∏_l=1^d p_slab({θ_ii, t}_l| {θ_ii, t-1}_l).
By using the conditional distributions as priors, we can control the change of amplitude or smoothness from the previous time step to the current one, e.g., from α_ii, t-1 to α_ii, t.
Compared with the dynamic spike-and-slab prior used in linear models <cit.>, our method does not constrain the slab prior to be a stable autoregressive process. Besides, we use a simpler but more flexible prior for γ_i,t, while the work in  <cit.> uses a prior conditional on the coefficients of the linear model.

The spike prior is responsible for shrinking the parameters to zero and cutting down the information transfer channel to the target. Common choices for this prior include point mass distribution, Laplace distribution, and Gaussian distribution with a small variance. Considering the shrinkage performance and optimization convenience, we choose the Laplace distribution as the spike prior, i.e.,

    p_spike(α_im, t) = 1/2ν_0exp(-| α_im, t |/ν_0),

where ν_0 is the length-scale for Laplace distribution. If we maximize the log-likelihood function to optimize the model, this prior will become a L_1 norm penalty and have the ability to shrink parameters.
The slab prior encourages the smoothness of parameter change. In this work, we consider two types of slab priors. The first one is a hard slab prior, 

    p_slab^hard(α_im, t| α_im, t-1) = 1/2ν_1exp(-| α_im, t -  α_im, t-1|/ν_1),

which encourages α_im, t to remain constant in a continuous period, approximating a piecewise model. In the second one, the parameters are allowed to change smoothly,

    p_slab^soft(α_im, t| α_im, t-1) =  1/√(2 πν_1)exp(-( α_im, t - ρα_im, t-1)^2/2ν_1),

where ν_1 is variance of Gaussian distribution, and ρ < 1 is an autoregressive coefficient.  A similar smoothing approach can also be found in <cit.>. These two slab priors make the current parameter value exactly or roughly concentrated around the previous value. Typically, we set ν_0 to be much smaller than √(ν_1) in the soft slab prior (or ν_1 in the hard slab prior) to make the two priors more separable and to put more penalty on sparsity. Besides, the values of α_im at multiple time steps before t can be included in the soft slab prior, e.g., p_slab^soft(α_im, t| α_im, t-1, α_im, t-2, ...). We choose the simplest form p_slab^soft(α_im, t| α_im, t-1) due to its wide application and robust performance in practice. 




At time t, η can be interpreted as the prior probability that α_im, t belongs to the slab process. It influences the strength of shrinkage effect that α_im, t bears. Ideally, for non-zero α_im, t, the posterior mean of γ_i,t should be close to one so that α_im, t is barely impacted by the spike prior. In the optimization algorithm developed in the next subsection, we will show that the estimated mean of γ_i,t is automatically adjusted based on the estimated α_im, t to avoid shrinking non-zero elements. This makes our method superior to traditional Lasso methods where the sparse penalty weights are identical for zero and non-zero parameters.













Finally, based on the above discussion, the whole hierarchical model of the proposed non-stationary MGP can be expressed as follows:

    y_(s), y_m | Φ_(s), Φ_m ∼𝒩 (0, K) 
       α_ii, t| α_ii, t-1∼p_slab(α_ii, t| α_ii, t-1),  i∈ℐ^S,
       θ_ii, t| θ_ii, t-1∼p_slab(θ_ii, t| θ_ii, t-1),  i∈ℐ^S, 
       α_im, t|γ_i, t, α_im, t-1∼ (1-γ_i, t) p_spike(α_im, t)  +γ_i, tp_slab(α_im, t| α_im, t-1), i∈ℐ, 
       γ_i, t|η∼  Bern(η), i∈ℐ, 
       θ_im, t| θ_im, t-1∼p_slab(θ_im, t| θ_im, t-1), i∈ℐ.






 §.§ Expectation-maximization-based optimization algorithm


The widely-used algorithm for a Bayesian model is Markov Chain Monte Carlo (MCMC) sampling, but it is computationally-inefficient for the proposed non-stationary model with considerable time-varying parameters. Therefore, we develop an efficient EM algorithm. Instead of directly maximizing the posterior p(Φ|y)=p(Φ_(s), Φ_m | y_(s), y_m), we proceed iteratively in terms of the complete log posterior log p(Φ, γ | y), where the binary parameters γ are treated as “missing data”. Since this function is not observable, in the Expectation-step (E-step), we calculate its conditional expectation given the observed data and the current estimated parameters. Then in the Maximization-step (M-step), we maximize the expected complete log-posterior with respect to Φ.
More precisely, the E-step and M-step at the (k+1)th iteration can be expressed as:

    E-step:  Q(Φ | Φ^(k)) = E_γ | Φ^(k), y{log p(Φ, γ | y)}, 
        M-step:  Φ^(k+1) = Φ argmax{ Q(Φ | Φ^(k)) }

where E_γ | Φ^(k), y (·) is the conditional expectation on posterior of γ, and Φ^(k) is the optimized parameters at the kth iteration. For simplicity, we use E_γ (·) to denote E_γ | Φ^(k), y (·). 

Based on Bayes' Theorem and the property of multivariate normal distribution, the expectation of log p(Φ, γ | y) can be as (derivation details can be found in Appendix C):

    E_γ{log p(Φ, γ | y) } 
    
    
    
    
    =    -1/2{y_(s)^T K_(ss)^-1y_(s) + log|K_(ss)| + (y_m-μ)^T Σ^-1 (y_m-μ) + log|Σ| }
        + ∑_i=1^m-1∑_t=2^n [  logp_slab(θ_ii, t| θ_ii, t-1) + logp_slab(α_ii, t| α_ii, t-1)   ] 
        + ∑_i=1^m ∑_t=2^n [ logp_slab(θ_im, t| θ_im, t-1) + (1-E_γγ_i, t) logp_spike(α_im, t)  
                  + E_γγ_i, tlogp_slab(α_im, t| α_im, t-1) ] + const.,

where μ=K_(sm)^T K_(ss)^-1y_(s) is the conditional mean of target given the sources and Σ=K_mm-K_(sm)^T K_(ss)^-1K_(sm) is the conditional covariance. 


In the E-step, since γ is only dependent on Φ_m, the posterior of γ_i, t is calculated as:

    p(γ_i, t| Φ^(k)_m) = p(α_im, t^(k) | γ_i, t) p(γ_i, t)/p(α_im, t^(k))∝ p(α_im, t^(k) | γ_i, t) p(γ_i, t) 
       =[ (1-γ_i, t) p_spike(α_im, t^(k)) +γ_i, tp_slab(α_im, t^(k)| α_im, t-1^(k)) ]  ·η^γ_i, t(1-η)^(1-γ_i, t).

Then the conditional expectation of γ_i, t can be updated as:

    E_γγ_i, t = ηp_slab(α_im, t^(k)| α_im, t-1^(k))/(1-η) p_spike(α_im, t^(k)) +ηp_slab(α_im, t^(k)| α_im, t-1^(k)),

















The posterior mean E_γγ_i, t can be interpreted as the posterior probability of classifying the current parameter α_im, t into a slab process as opposed to a spike process, based on the past value α_im, t-1. For example, we set η to 0.5 as a non-informative prior, and take a small ν_0 (e.g., 0.01) for the spike prior and a large ν_1 (e.g., 0.1) for the soft slab prior. Based on eq: gamma update, supposing that (α_im, t-α_im, t-1)^2 is small and |α_im, t| is large, the expectation E_γγ_i, t will tend towards one, indicating that α_im, t is more likely from the slab prior. 
On the other hand, if |α_im, t| is small, E_γγ_i, t is close to zero, enforcing strong shrinkage on α_im, t. 

In the M-step, we can optimize the objective function in eq: EM objective with various gradient ascent methods, such as (stochastic) ADAM, L-BFGS, etc. <cit.>. 
This objective function is actually a standard Gaussian process log-likelihood with additional regularization terms. The regularization terms penalize the difference between parameters at successive time points and shrink the amplitude parameters to facilitate source selection. The weights of the regularization terms are modulated by the expectation of γ. Ideally, for non-zero α_im, t, E_γγ_i,t will equal to one, so no shrinkage effect will be placed on it. In other words, the strength of sparse penalty is automatically adjusted through eq: gamma update, and this adjustment has explicit statistical interpretability. 

In our case studies, we find the algorithm converges rapidly, e.g., achieving convergence after only five iterations. 
The whole algorithm is summarized in Algorithm <ref>, where an ADAM method <cit.> is utilized in the M-step.
Note that the large parameter space of Φ poses a challenge in identifying the modes of the full posterior. To speed up the convergence, we propose to initialize the source parameter Φ_(s) by maximizing the sum of sources’ marginal log-likelihood and source parameter prior:

    Φ_(s) max    log p(y_(s)|  Φ_(s)) + log p ( Φ_(s) )


For target parameters, we find simple random initialization can achieve satisfactory performance in experiments.



 §.§ Computational Challenge



There are three main computational challenges that we need to address. The first one is the calculation of integration for a convolution kernel. To avoid an intractable integration for the covariance, we utilize the Gaussian kernel in eq: Gaussian kernel and derive closed-form covariance functions in Equations (<ref>- <ref>).

The second challenge is calculating the inverse of covariance matrix, which is a critical issue for all GPs. The computational complexity of Algorithm <ref> is approximately O(∑_i=1^m n_i^3), where n_i is the number of data points for each output. In comparison, the computational cost for traditional MGP is O((∑_i=1^m n_i)^3) with the same size of data and non-separable covariance.  The main computational load of our model is in calculating the inverse of K_(ss) and Σ in the source marginal distribution  log p(y_(s)|  Φ_(s) ) and the target conditional distribution log  p(y_m | Φ_m, y_(s), Φ_(s)) respectively in eq: EM objective. Since all the latent processes are independent of each other, K_(ss) is a block-diagonal matrix and the complexity is reduced from O(∑_i=1^m-1 n_i^3) to O(∑_i=1^m-1 n_i^3). The calculation of the inverse of Σ is O(n_m^3). Therefore, the overall computational complexity is reduced to O(∑_i=1^m n_i^3).

Finally, it is a hard task to estimate a considerable number of time-varying parameters. Therefore, we develop the EM-based algorithm to fit the model rather than using a sampling method. Based on the results of a non-stationary linear model <cit.>, the MCMC and EM algorithm lead to very close prediction errors, but the running time of MCMC is about ten times longer than that of the EM.





 §.§ Tuning Parameter Selection

 The tuning parameters for our model is ν_0, ν_1 (for the hard slab prior), and ν_2 (for the soft slab prior). Here, since the key of our method is selecting the most informative sources, we propose to maximize the following criterion:

    B(ν) = N log p(y|X) - log(N) c_ν(α_m),

where N is the number of data, log p(y|X) is the log-likelihood for both the source and the target outputs, and c_ν(α_m) is the number of nonzero elements in α_m given ν. Similar criterion is proposed in <cit.>. Note that ν = {ν_0, ν_1} for the hard slab prior and  ν = {ν_0, ν_2} for the soft slab prior. This criterion is similar to the Bayesian Information Criterion (BIC). The first term tends to select a more complex model with a larger likelihood for all outputs, while the second term prefers simpler models where less sources chosen to transfer information to the source. To reduce the computation, we first determine the ratios r_1 = ν_1/ν_0 and r_2 = ν_2/ν_0 to make the spike and slab priors separable <cit.>. Then we design a two-dimension search grid for (ν_1, ν_1/r_1) with the hard slab or (ν_2, ν_2/r_1) with the soft slab. The optimal value of ν is searched over the two-dimensional grid. For example, we set r_1∈{5, 10} and (ν_1, ν_0) ∈{(1/5, 1/25), (1/5, 1/50), (1/10, 1/50), (1/10, 1/100), (1/15, 1/75), (1/15, 1/150)} for a hard slab prior.



 §.§ Model Prediction

Since the EM algorithm only estimates the value of parameters at the observed time points, given a new input x_t^* of interest, we first need to estimate the target parameter α_m, t^* and  θ_m, t^* at the new time point t^*, then derive the predictive distribution of y_m (x_t^*). 



  §.§.§ Forecasting.

In the forecasting task, t^* > n. The estimated value of θ_im, t^* is,

    θ_im, t^* = {θ_im, n,      for  hard  slab  prior, 
    ρ^t^*-nθ_im, n,     for soft slab prior,
    .

which is actually the mode of p_slab(θ_im, t^* | θ_im, n). As for α_im, t^*, if E_γγ_i, n≥ 0.5, we consider it from the slab process and estimate it using the same method as θ_im, t^*. Otherwise, it is classified to a spike process and shrunk to zero <cit.>. 



  §.§.§ Recovery.

In the recovery task, the target data are unobserved at some specific time points and we aim to recover the missing data, i.e., 1<t^*<n. Define t_be and t_af to be the nearest observed time points before and after t^* respectively. Denote the nearest observation time to t^* as t_near, i.e., 

    t_near = t∈{t_be, t_af} argmin |t-t^*| .

As the parameters before and after t_* are already optimized by the EM algorithm, the estimation of θ_im, t^* becomes:

    θ_im, t^* = {θ_im, t_near,      for  hard  slab  prior, 
     LSE (θ_im, t_be, θ_im, t_af),     for soft slab prior,
    .

where LSE( · ) represents a least-square estimation for auto-regressive process introduced in <cit.>.
We also let α_im, t^* = 0, if E_γγ_i, t_near<0.5. Otherwise, we estimate its value in the same way as θ_im, t^*. 

Then, given the parameter α_m, t^* and θ_m, t^*, and the new input point x_t^*, the joint distribution of y_m( x_t^* ) and observations y_m can be expressed as:

    [          y_m; y_m( x_t^* ) ]∼𝒩[ [     μ; μ_t^* ],
    	[          Σ      Σ_t^*;    Σ_t^*^T Σ_t^*, t^* ] ],

where μ_t^*=K_(s*)^TK_(ss)^-1y_(s), Σ_t^* = K_m*-K_(sm)^TK_(ss)^-1K_(s*), and Σ_t^*, t^* = K_t^*, t^*-K_(s*)^TK_(ss)^-1K_(s*).





In the above equations, K_(s*) is the cross-covariance matrix of the sources and the new input x_t^*, K_(m*) is the covariance of the target observation and the new point, and K_t^*, t^* is the variance at x_t^*.
Then, the posterior distribution of y_m( x_t^* ) can be derived as:

    y_m( x_t^* ) ∼𝒩(μ_t^*+Σ_t^*^T Σ^-1 (y_m - μ),
    Σ_t^*, t^* - Σ_t^*^T Σ^-1Σ_t^*).




§ NUMERICAL STUDY


In this section, we verify the effectiveness of the proposed non-stationary MGP with the spike-and-slab prior (denoted as DMGP-SS) using synthetic data. In Section <ref>, we briefly describe the general settings for the numerical study and benchmark methods. In Section <ref>, we introduce the design of simulation cases, where the cross-correlation of the sources and the target are dynamic and sparse. In Section <ref>, we demonstrate our model's capability in detecting the underlying correlation pattern as well as improving target prediction performance on the synthetic data.



 §.§ General settings


Similar to the assumption made in the model development section, we generate m sequences consisting of m-1 sources and 1 target, sampled at the same timestamps. The input space is simply time. To investigate the source selection capability of DMGP-SS, we assign only m_t < m-1 sources to be correlated with the target at time t. Besides, a source will remain either correlated or uncorrelated continuously for a certain period of time.

For comparison, we consider three benchmarks:


  * GP. The target is modeled using a single-output GP, with a squared-exponential covariance function.

  * MGP-L1. It is a state-of-art static method introduced in <cit.>. MGP-L1 models the target and sources in one MGP model, with the same covariance structure as in eq: covariance matrix. The scaling parameters {α_im}_i=1^m-1 are penalized by L_1 term to achieve source selection. The regularized log-likelihood of this model is:

    logℱ = -1/2y^T K^-1y - 1/2log |K|  - λ∑_i=1^m-1|α_im| - const.

where K is calculated using the static covariance functions in <cit.> and λ is a tunning parameter.  

  * DMGP-GP. This is a state-of-art non-stationary MGP model, which constructs an LMC model for all outputs and assumes the hyper-parameters follow other GPs <cit.>. More details can be found in the Appendix D. In this model, the covariance for the m outputs is

    cov[y(x_t),  y(x_t^')] = A_t A_t'^T k(x_t, x_t^') +  diag{σ_i},

where A_t A_t^'^T ∈ℝ^m × m is the correlation matrix of m outputs, and diag{σ_i} is the diagonal matrix with  {σ_i}_i=1^m. In this study, we focus on the correlation between the sources and target, which corresponds to the last column of A_t A_t^T (except the last element (A_t A_t^T)_mm), i.e., (A_t A_t^T)_0:m-1, m.  


All methods are implemented in Python 3.8 and executed on a computer with an Inter(R) Core(TM) i5-7400 CPU with 3.00GHz and 16GB RAM. Both GP and MGP-L1 are implemented using gpflow <cit.>.  DMGP-GP is implemented using TensorFlow and optimized with ADAM <cit.>, with a maximum iteration limit of 500. The EM algorithm for DMGP-SS is also based on Tensorflow, and we use stochastic ADAM with four batches in the M-step. We set the k_out and k_in in Algorithm <ref> to 5 and 400 respectively. 

For MGP-L1, the weight of L_1 penalty λ is a tuning parameter. For DMGP-GP, we use square exponential functions for k_α and k_θ. For simplicity, we apply the same tuning parameters for both kernels, i.e., the amplitude α^# and length-scale θ^#. Those parameters are tuned by cross-validation. In the case of DMGP-SS, the prior sparsity parameter η is set to 0.5. We repeat each case 50 times and report the prediction performance through averaging the results. 



 §.§ Simulation cases


The main objective of this section is to demonstrate the effectiveness of our method in capturing the non-stationary and sparse cross-correlation between the sources and the target. For simplicity, we hold the other characteristics constant over time, e.g., the smoothness of each output. 
Specifically, we design two simulation cases with different cross-correlation patterns. The first case involves a piece-wise constant cross-correlation, while the second case has a smoothly-changing correlation. In each case, the input data {x_t}_t=1^130 are evenly spaced in [1, 130]. The observed data are generated from sine functions with measurement noise ϵ_t ∼𝒩(0,0.3^2).  

Case 1. 
In this case, we define four kinds of source functions:

    Y_1(x_t)=3sin(π x_t /20 + e_1) + ϵ_t,  
       Y_2(x_t)=2sin(2 π x_t /20 + e_2) exp[0.5(x_t%40-1)] + ϵ_t, 
       Y_3(x_t)=3sin(4 π x_t /20 + e_3) + ϵ_t, 
       Y_4(x_t)=2sin(5 π x_t /20 + e_4) + ϵ_t,

where e_i∼𝒩(0,0.2^2) is a random phase to make the sampled outputs different from each other, and “%” represents the reminder operation. The term exp[0.5(x_t%40-1)] is used to deviate the shape of Y_2 from the standard sine function. In each experiment, we generate 4k sources through sampling each kind of function k times, i.e., m=4k+1. Specifically, the sources {y_i+4j |  0 ≤ j ≤ k-1 } are sampled from Y_i. 
Then, we define the dynamic target output as: 

    f_m(x_t)=a_1, tsin(π x_t/20)+a_2, tsin(2 π x_t/20) + a_3, tsin(4π x_t/20) + a_4,tsin(5 π x_t /20)

In this case, we simulate a piece-wise constant cross-correlation by setting: 

    a_1, t = (2 + 2a_1) I_t<40, 
       a_2, t = (2+ 2a_2) I_40 ≤ t < 80 + (1 + a_2) I_80 ≤ t ≤ 130,
       a_3, t = (1+ a_3) I_80 ≤ t ≤ 130, 
       a_4,t = 0.

Therefore, there are three segments, [0,40), [40, 80), and [80,130]. Only the 1st, the 2nd, and the 2nd and 3rd sources are correlated to the target in the three periods, respectively. The other sources remain uncorrelated to the target at all times.

Case 2.
Compared with Case 1, we only change the coefficients {a_i, t}_i=1^3 into smoothly-changing ones in this case. Specifically, we let them vary along sine-cosine trajectories, 

    a_1, t =  [(2+a_1)cos(π t / 120) + 0.5] I_t<40,
       a_2, t = [(2+a_2)sin(π t / 120 - π/6) +0.5] I_40 ≤ t < 130,
       a_3, t = [(2+a_3) sin(π t / 120 - π/2) +0.5] I_80 ≤ t < 130.










In all cases, we set k=1, 4 to generate four and sixteen source outputs for each case. In order to compare the prediction performance of different methods, we randomly remove three short sequences of length 10 from [10, 30], [50, 70] and [90, 110] respectively. These 30 data points are treated as missing data, while the others are used as training data. 



 §.§ Simulation results


To begin with, we demonstrate that the proposed DMGP-SS is capable of capturing the dynamic and sparse correlation between the sources and the target. 
fig: simulation scale illustrates the estimated {α_im}_i=1^4 for MGP-L1, the {α_im, t}_i=1^4 for DMGP-SS, and the estimated (A_t A_t^T)_1:4, m for DMGP-GP in Case 1 and 2 with four sources. And fig: case1 prediction visualizes the sources and target prediction results in Case 1 with four sources. DMGP-SS is implemented with the hard and soft slab priors for Case 1 and 2 respectively.  Note that the value of a_i,t and α_im, t are not identical, since a_i,t is a linear combination weight rather than the correlation parameter in MGP.  

Overall, DMGP-SS successfully recovers the true dynamic structural sparsity of correlation shown in the first column of fig: simulation scale. Firstly, DMGP-SS tracks closely the periods of correlation between each source and the target, achieving a dynamic selection of sources. Since the target's characteristics do not change abruptly (as shown in fig: case1 prediction), it is reasonable that the estimated correlation change points are about ten time-steps before or after the designed change times. Second, the hard slab prior encourages nearly piece-wise correlation, while the soft slab prior allows smoothly changing correlation. Due to the appropriate selection of sources at different times in fig: case1 prediction, the proposed DMGP-SS achieves precise prediction with the lowest prediction uncertainty. This highly improves the confidence of decision making when using the recovered series. The difference on confidence interval of three MGP models is due to that the uncorrelated sources 'poison' the correlation structure and decrease the influence of the truly-correlated sources, resulting in a lower value of variance reduction term Σ_t,*^T Σ^-1Σ_t,* in posterior prediction eq:posterior distribution.

In contrast, another non-stationary method DMGP-GP fails to estimate a sparse structure in the source-target correlation since the GP prior on parameters lacks the shrinkage effect. The proposed DMGP-SS addresses this problem through combining the smooth slab prior and the shrinking spike prior.
Regarding MGP-L1, although it can estimate a sparse structure of α_m, the estimated values of non-zero parameters are constant over time and cannot reflect the change of correlation. As a result, it even performs worse than GP in recovering the target output in  fig: case1 prediction.



Another advantage of DMGP-SS is the adaptive adjustment of the spike-and-slab combination weight, E_γγ_i, t. In our settings, ν_0^-1 in the spike prior is much larger than ν_1^-1 in the hard slab prior (or √(ν_1^-1) in the soft slab prior), to put more penalty on the correlation sparsity. For example, we set ν_0 = 0.02 and ν_1 = 0.1 in Case 1. However, this sparse penalty does not cause significant bias on non-zero α_im, t because of the automatically adjusted E_γγ_i, t in the EM algorithm. Specifically, we initialize it with 0.99 to barely shrink parameters at the beginning. Then E_γγ_i, t is updated in the E-step of the EM algorithm. fig: simulation gamma shows its estimated value after five iterations. For the correlated sources (e.g., the first source during t∈[0,50]), their corresponding E_γγ_i, t is very close to 1, so they bear negligible shrinkage effect from the spike prior. In contrast, for the uncorrelated sources (e.g., the first source after t = 50), E_γγ_i, t is approximately 0.2, implying a substantial shrinkage effect. The value 0.2 can be derived based on eq: gamma update. For consecutive zero elements, p_slab=η (2ν_1)^-1, and p_spike=(1-η) (2ν_0)^-1, resulting in E_γγ_i, t≈ν_1^-1/ν_0^-1.





tab: simulation MAE summarizes the results of 40 repetitions for each case.  
DMGP-SS outperforms all the other methods in both cases. And the increasing of source number does not affect the prediction accuracy,   demonstrating its remarkable robustness in dealing with moderate data. 
MGP-L1 exhibits slightly higher prediction accuracy than GP. Because the static covariance structure limits its ability to transfer accurate information at all times. Under some circumstances, this limitation will cause a negative transfer effect on the learning of target (for example, the result shown in fig: case1 prediction).
DMGP-GP has a better performance on target prediction than GP and MGP-L1, due to the ability to model dynamic correlation. However, it does not achieve the same prediction accuracy as DMGP-SS in these cases. On the one hand, it lacks the ability to exclude the impact of uncorrelated sources. On the other hand, DMGP-GP is an extension of LMC and uses the same function k(x_t, x_t^') to model the auto-covariance of every output. This feature makes it unsuitable for our cases where the source covariance functions have four kinds of length-scales. Nevertheless, the proposed DMGP-SS models each source with separate kernels and latent functions, highly increasing its flexibility.

 tab: computational time lists the computational time of the four methods in Case 1. Between the two non-stationary methods, DMGP-SS requires much less computation time than DMGP-GP. This exactly verifies the analysis in <ref> that our model can save a large mount of time than the traditional non-stationary MGP method, due to a block-sparse covariance matrix. In fact, our model is scalable for larger size of data, which is described in Appendix E.




§ CASE STUDY



In this section, we apply DMGP-SS to two cases: human movement signal modeling and control policy optimization. In the first case, these signals are recorded by sensors attached to different joints, such as hands and feet. As the movement of joints are dynamically correlated across different gestures <cit.>, it is reasonable to utilize a non-stationary method to capture the varying correlation and leverage information among the signals of joints. Regarding the control policy iteration, we study on a classical reinforcement learning problem, mountain-car. We evaluate the performance of DMGP-SS on leveraging knowledge between difference systems when the environment is non-stationary. 



 §.§ Movement Signal Modeling



  §.§.§ Data Description

We use the MSRC-12 gesture dataset <cit.> consisting of sequences of human skeletal body movement. Twenty sensors are distributed across the body, and each sensor can record three-dimensional coordinates of joint positions. The motion signal has a sampling frequency of 30 Hz and an accuracy of ± 2cm. The dataset comprises 12 different gestures performed by 30 participants. The experiment involves each participant starting from a standing position, performing one gesture, and returning to the standing state. This process repeats several times continuously within each sample. 

To demonstrate the effectiveness of DMGP-SS, we connect the instances of three gestures (“Goggles”, “Shoot”, and “Throw”) performed by the same individual.  fig: gesture skl shows the snapshots of the standing position and the selected gestures.
In the first two gestures, the participant stretches both arms in front of him to perform searching or shooting motions. In the third gesture, the participant only uses his left arm to make an overarm throwing movement. In these gestures, the main acting joints are hands, wrists, and elbows, where the trajectories of hands and wrists are almost identical. Therefore, we select the movement signals of four joints (left wrist, left elbow, right wrist and right elbow) as twelve outputs. We choose the z-coordinate movement of the left elbow as a target output, while using the remaining eleven movement signals as sources.







We select two 120-frame-long instances for each gesture and down-sample each instance to 30 frames. Therefore, there are 180 points for each output. To eliminate the difference in initial coordinate values, we reset the initial three-dimensional coordinate value to (0, 0, 2) across different recordings. Additionally, we re-scale all outputs to [-2, 2].  fig: real data displays the 12 outputs and the change of joints' positions.












  §.§.§ Results

Intuitively, the cross-correlation between the source and target signals should remain constant for a single gesture, so a hard slab prior is used in DMGP-SS. All other settings for this case are identical to those used in the simulation studies. We still simulate consecutive data missing for the target. From the 60-points-long time-series of each gesture, we randomly remove a sequence of length ten as missing data.

First of all, fig: real correlation shows the estimated correlation between the sources and the target. MGP-L1 selects six sources (the third 'R-W-z', the forth 'R-E-x', the sixth 'R-E-z', the seventh 'L-W-x', the ninth 'L-W-z', and the eleventh 'L-E-y')  as the correlated signals for the whole time period, in which the ninth source has the strongest correlation. DMGP-SS selects the sixth source (`R-E-z') and the ninth source (`L-W-z') as correlated sources when 0 ≤ t ≤ 120 and 120 ≤ t ≤ 180, respectively. DMGP-GP does not provide a sparse estimation of cross-correlation. Among them, the proposed DMGP-SS accurately captures the underlying physical movements of each gesture. In the “Google” and “Shoot” gestures, both elbows have almost the same trajectory. In the “throw” gesture, the left elbow's movement is highly correlated with that of the left wrist. These findings align well with the signals shown in fig: real data. On the contrary, limited by the static correlation structure, MGP-L1 can only select some sources and force them to be correlated with the target all the time, but such signals do not exist in the dynamic environment. For DMGP-GP, the results cannot provide us an intuitive understanding on the joints' relationship.

















fig: real prediction displays the recovered target signal in one experiment. Notably, DMGP-SS accurately recovers the missing data with high precision. Besides, it has the minimal uncertainty because it can selects the most correlated sources at each time and such a high correlation improves the confidence on prediction results.  Conversely, the predictions of MGP-L1 and DMGP-GP display an obvious bias in the first gap and higher uncertainty for all gaps. Besides, similar to the results of numerical studies, the proposed DMGP-SS also gives predictions with the lowest uncertainty, which significantly improves the confidence of decision making.

 We further repeat the experiments 36 times. <ref> compares the prediction accuracy of four methods in terms of both the MAE and the continuous-ranked-probability-score (CRPS). The CRPS measure is a widely used metric to evaluate the probabilistic forecast for a real-valued outcome <cit.>: 

    CRPS = n_test^-1∑_i=1^n_test∫[Φ( ŷ_i ) - 1_ŷ_i ≥ y_i)]^2 dŷ_i

where ŷ_i is the predicted output and Φ is the cumulative density function (CDF) of the predicted distribution. A low CRPS value suggests that the predicted posterior are concentrated around the true value, indicating a high probabilistic forecast performance. As expected, DMGP-SS outperforms the other methods due to its ability to capture the dynamic and sparse correlation accurately. MGP-L1 performs better than GP, benefiting from the information borrowed from the other sources. However, it cannot model a dynamic correlation, resulting in lower prediction accuracy than DMGP-SS. Regarding DMGP-GP, although it captures the change of correlations between the target and some sources, its prediction accuracy is even lower than GP. This result may be attributed to that non-sparse correlations lead to potential negative transfer effects. Besides, since the sources' smoothness is heterogeneous, it is inappropriate to use the same auto-covariance function to model all outputs. 





 §.§ Control Policy Optimization


In reinforcement learning problems, there are five important elements for the agent: state s, action a, reward r, value V and policy W. Starting from one state s, the agent takes an action a, transitions into a new state u and gets an immediate reward r(u) from the environment. In general, a reinforcement learning framework includes two primary procedures: 

    
  * Model estimation: estimating the state transition function U(s, a) based on the observed transition samples [(s, a), u].
    
  * Policy iteration: iteratively estimating the state value V(s) (the long-term expected reward) for each state and improving the control policy W(s).

More details could be found in the works on GP-based reinforcement learning <cit.>. 

In this study, we employ the well-known reinforcement learning problem, mountain-car, to demonstrate the application of our model in decision making. fig: mountain car illustrates such a problem. A car begins in a valley and aims to reach a goal position of the right-most hill. Due to the steep slope of the hill, the car cannot directly accelerate to the goal position. Instead, it needs to drive up the opposite side, turn back, and accelerate to reach the goal position. In system, the agent state is described by the position s^pos∈[-1.2, 1.0] and the velocity s^vel∈[-0.07, 0.07] of the car, i.e., s=(s^pos, s^vel). The agent action is a horizontal force a ∈[-1, 1]. The car starts at the state s_init=(s^vel_init, s^vel_init) with s^pos_init∈ [-0.6, -0.5] and s^vel_init=0, aiming to reach the goal state s_goal=(0.45, 0). The reward function is the probability density function of N(s_goal, diag{0.05^2, 0.0035^2}). The dynamic equation of this system is approximated by:

    s^vel_t = s^vel_t-1 + P · a_t-1 - G ·cos(3 · s^pos_t-1) 
        s^pos_t = s^pos_t-1 + s^vel_t

where P is the horizon power unit and G is the vertical force unit.


We consider the control problem involving a car in a non-stationary environment with limited transition samples. During t∈[0,20), the target car runs in an environment with (P,G)=(0.009, 0.0015) and we get 20 random samples. However, at t=20, an unknown change occurs, altering  the environment factors to (P, G)=(0.0011, 0.0026). We sampled another 20 random samples during t∈ [20,40). After t=40, we start to control the car to reach the goal position and the environment does not change any more. Given that those samples are too limited to build an accurate a transition model, we transfer information from two historical source datasets, each consisting of 200 samples from stationary environments with (P,G)=(0.01, 0.0015) and (P, G)=(0.001, 0.0025), respectively. We can see that the target environment is close to the first source environment before t=20, and is close to the second one after that.



Algorithm <ref> summarizes the workflow of reinforcement learning, where we employ the DMGP-SS as a transition model. For simplicity, we model the source transition data with stationary GPs in DMGP-SS. The sources and the target are expressed as:

    u_i (s, a)     =α_ii g_ii(s, a)∗ z_i(s, a) + ϵ_i , i ∈ℐ^S, 
    
    u_m (s, a)     = ∑_j ∈ℐα_jm, tg_jm, t(s, a) ∗ z_j(s, a) + ϵ_m.

Here, u is the position or velocity state to where the agent transitioned from a state s after taking an action a. After fitting this model both the target and source samples, we train a GP for the value model V(s) with even-spaced support points S_supp and rewards r(S_supp). Based on the transition model U(s, a), we iteratively improve the policy W(s) and value model V(s) until the prediction of V(s) converges. Details on the policy improvement procedure can be found in <cit.>. While we focus on the offline setting in this case, it is worth noting that this framework can be readily extended to accommodate an online setting, wherein the training sample consists of the visited states and the transition models are updated every few steps.





















We choose two reinforcement learning benchmark methods <cit.> with the stationary GP and MGP as the state transition model, respectively. The maximum execution steps are 600 for each method. 
In tab: RL MAE and Distance, we report the the mean of absolute distances to the goal state for three methods. The DMGP-SS-based control policy has the shortest average distance to the goal state in 600 steps. Specifically, fig: rl case path compares the position of the car controlled by the three policies. With DMGP-SS as the transition model, the car reaches the goal position and stays there after about 250 time steps. However, the other methods cannot find a good policy to reach the goal position within 600 moves, since their stationary transition models cannot account for the environment change and make the policy iteration hard to converge.





In tab: RL MAE and Distance, we further compare the predictive MAE on state transition for three methods. The proposed method has a significantly lower prediction error on velocity state transition than the the other methods, since it can capture the change of environment and transfer information from the related sources. fig: rl correlation illustrates the estimated α_m from DMGP-SS trained on the velocity transition data. We can find that the proposed method successfully finds that the correlation between the sources and the target changes at t=20. Therefore, during the policy improvement stage, it can leverage information from the similar source (the second one) and avoid the negative transfer from the uncorrelated source (the first one). Regarding the prediction error on position transition, although RL-DMGP-SS has a higher MAE than the other benchmarks, the difference is minor considering the position range [-1.2, 0.6]. Therefore, our method can provide the best control policy.







§ CONCLUSION


This paper proposes a flexible non-stationary multi-output Gaussian process for modeling multivariate data in transfer learning. The novelty of our approach lies in its ability to capture dynamic and sparse cross-correlations between sources and targets. We achieve this by allowing correlation parameters to follow a spike-and-slab prior, where the slab prior ensures correlation variation over time, and the spike prior encourages parameters to shrink to zero, eliminating negative transfer effects from uncorrelated sources. The ratio of these two priors is automatically adjusted in the proposed EM algorithm, preventing the shrinkage effect on non-zero correlation parameters. Through the experiments on simulation and human gesture dataset, we demonstrate that our method is well-suited for both capturing non-stationary characteristics and mitigating negative transfer. 

The proposed data-driven method provides a powerful tool for researchers and engineers to select the most informative sources to transfer knowledge. Except high-dimensional time-series modeling, our approach could also find applications in both sequential decision making and change point detection. For instance, transfer learning has arisen to handle the critical challenge of sparse feedbacks in reinforcement learning (RL), a popular framework for solving sequential decision making problems. However, negative transfer is still a notable challenge for multi-task reinforcement learning and it is risky to naively share information across all tasks <cit.>. Therefore, we embedded  our model with the offline RL transfer framework to automatically select correlated sources to share knowledge in a non-stationary environment <cit.>. In the future, we can further develop our method to account for an online RL task. Besides that, the estimated dynamic correlation among outputs can help us to understand the structure of time-series even with some missing data, as well as to detect some structural change points for subsequent decision making.

Many extensions are possible for our model. Although the proposed methodology is flexible enough to capture complex and dynamic correlation, scaling it up to large datasets is computationally challenging. Possible solutions include utilizing a sparse approximation for the covariance matrix or developing a more efficient optimizing algorithm. 

In addition, the proposed EM optimization algorithm solely provides estimated values of the model parameters without incorporating any uncertainty measurement. To address this issue, we can use variational inference to approximate the true posterior distribution of parameters, thereby capturing the inherent uncertainty associated with these parameters. The approximated uncertainty would be propagated into the prediction at new points.





[4]



§ DMGP-SS WITH MISSING DATA

In general, observation may be missing at some time points within the range of 1 ≤ t ≤ n. Under such a circumstance, the observation number for each output is n_i ≤ n. To account for the missing data, we re-denote the data for the ith output as X_i = (x_i, t_1, ..., x_i, t_i, ..., x_i, t_n_i)^T and y_i = (y_i, t_1, ..., y_i, t_i, ..,y_i, t_n_i)^T, where t_i represents the observation time index.

Compared with the model presented in the main text, the general DMGP-SS requires only adjustments on the parameter priors. The general model can be expressed as follows,  

    y_(s), y_m | Φ_(s),  Φ_m∼𝒩 (0, K) 
       α_ii, t_i| α_ii, t_i-1∼p_slab(α_ii, t_i| α_ii, t_i-1),  
       θ_ii, t_i| θ_ii, t_i-1∼p_slab(θ_ii, t_i| θ_ii, t_i-1),  
       α_im, t_i|γ_i, t_i, α_im, t_i-1∼ (1-γ_i, t_i) p_spike(α_im, t_i)  
                        +γ_i, t_ip_slab(α_im, t_i| α_im, t_i-1), 
       γ_i, t_i|η∼  Bern(η), 
       θ_im, t_i| θ_im, t_i-1∼p_slab(θ_im, t_i| θ_im, t_i-1).


The spike prior is the same as that in the main text. The slab priors are re-defined as,

    p_slab^hard(α_im, t_i| α_im, t_i-1) = 1/2ν_1exp(-| α_im, t_i -  α_im, t_i-1|/ν_1),


    p_slab^soft(α_im, t_i| α_im, t_i-1) =  1-ρ/(1-ρ^t_i - t_i-1)√(2 πν_1)exp(-( α_im, t_i - ρ^t_i - t_i-1α_im, t_i-1)^2 (1-ρ^2)/2ν_1 (1-ρ^2(t_i - t_i-1))),

where the soft prior is derived based on the property of an auto-regressive process. 




§ PROOF OF PROPOSITION 1

The proposed non-stationary MGP covariance matrix Eq. (9) is positive-definite, i.e., ∀y≠0,

    y^T Ky > 0.


Proof of Proposition 1
Recall that the covariance functions are generated by the convolution of kernel functions:

    cov_ii^f(x_t, x_t^')     = α_ii,tα_ii,t'∫ g_ii,t(x_t - u) g_ii, t'(x_t' - u) du
     cov_im^f(x_t, x_t^')     = α_ii,tα_im,t'∫ g_ii,t(x_t - u) g_im, t'(x_t' - u) du
     cov_mm^f(x_t, x_t^')     = ∑_j=1^m α_jm,tα_jm,t'∫ g_jm,t(x_t - u) g_jm, t'(x_t' - u) du


Decompose this quadratic form as follows,

    y^T Ky    = ∑_1≤ i ≤ m∑_1≤ j ≤ m∑_t ∑_t' y_i,ty_j,t' [ cov_i,j^f (x_t, x_t^') + σ_i^2 𝕀_i=j, t = t'] 
       =  ∑_1≤ i ≤ m-1∑_t ∑_t' y_i,ty_i,t' cov_ii^f (x_t, x_t^')  + 2 ∑_1≤ i ≤ m-1∑_t ∑_t' y_i,t y_m, t' cov_im^f (x_t, x_t^') 
          + ∑_t∑_t' y_m,t y_m, t' cov_mm^f (x_t, x_t^') + ∑_i ∑_t y_i,t^2 σ_i^2 
       =  ∑_1 ≤ i ≤ m-1{∑_t ∑_t' y_i,ty_i,t'α_ii,tα_ii,t'∫ g_ii,t(x_t - u) g_ii, t'(x_t^' - u) du. 
                    + ∑_t ∑_t' y_m,ty_m,t'α_im,tα_im,t'∫ g_im,t(x_t - u) g_im, t'(x_t^' - u) du
                   . + 2∑_t ∑_t' y_i,ty_m,t'α_ii,tα_im,t'∫ g_ii,t(x_t - u) g_im, t'(x_t^' - u) du}
          + ∑_t ∑_t' y_m,ty_m,t'α_mm,tα_mm,t'∫ g_mm,t(x_t - u) g_mm, t'(x_t^' - u) du + ∑_i ∑_t y_i,t^2 σ_i^2 
       =  ∑_1 ≤ i ≤ m-1{∫[ ∑_t  y_i,tα_ii,t g_ii,t(x_t - u) ∑_t' y_i,t'α_ii,t' g_ii, t'(x_t^' - u)  . . 
                        +  ∑_t y_m,tα_im,t g_im,t(x_t - u) ∑_t' y_m,t'α_im,t'  g_im, t'(x_t^' - u)  
                       . . + 2 ∑_t y_i,tα_ii,t g_ii,t(x_t - u) ∑_t' y_m,t'α_im,t' g_im, t'(x_t^' - u) ] du}
          + ∫∑_t  y_m,tα_mm,t g_mm,t(x_t - u) ∑_t'y_m,t'α_mm,t' g_mm, t'(x_t^' - u) du + ∑_i ∑_t y_i,t^2 σ_i^2 
       =  ∑_1 ≤ i ≤ m-1{∫[ ∑_t  y_i,tα_ii,t g_ii,t(x_t - u) + ∑_t' y_m,t'α_im,t' g_im, t'(x_t^' - u) ]^2 du}
          + ∫[ ∑_t  y_m,tα_mm,t g_mm,t(x_t - u) ]^2 du + ∑_i ∑_t y_i,t^2 σ_i^2  > 0

Proof completes.




§ DERIVATION OF THE OBJECTIVE FUNCTION IN M-STEP

Based on Bayes theorem, the parameter posterior can be expressed as:

    p(Φ, γ | y) ∝ p(y|Φ)p(Φ|γ)p(γ).

And based on the theory of multivariate Gaussian distribution, we have:

    p(y|Φ) = p(y_(s), y_m | Φ_(s), Φ_m) 
       =𝒩(
    [ [ y_(s);   y_m ] ] |
    [ [ 0; 0 ] ],
    [ [   K_(ss)   K_(sm); K_(sm)^T    K_ mm ] ]) 
       = 𝒩( y_(s) | 0, K_(ss)) 𝒩( y_m | μ, Σ),

where μ=K_(sm)^T K_(ss)^-1y_(s) is the conditional mean of target given the sources and Σ=K_mm-K_(sm)^T K_(ss)^-1K_(sm) is the conditional covariance. 

Therefore, the objection function can be derived as:

    E_γ{log p(Φ, γ | y) }
    
    =    E_γ{log p(y | Φ, γ) p(Φ, γ) }  + const. 
    
    =   log p(y_(s)|  Φ_(s) )  +  log  p(y_m | Φ_m, y_(s), Φ_(s))   
       + log p(Φ_(s)) + E_γ{log p(Φ_m| γ) + log p(γ) }  + const. 
    
    
    =   log p(y_(s)|  Φ_(s) )  +  log  p(y_m | Φ_m, y_(s), Φ_(s))   
       + log p(θ_(s))  + log p(α_(s))  + log p(θ_m) + E_γ{log p(α_m| γ) }  + const. 
    
    =    -1/2{y_(s)^T K_(ss)^-1y_(s) + log|K_(ss)| + (y_m-μ)^T Σ^-1 (y_m-μ) + log|Σ| }
        + ∑_i=1^m-1∑_t=2^n [  logp_slab(θ_ii, t| θ_ii, t-1) + logp_slab(α_ii, t| α_ii, t-1) + logp_slab(θ_im, t| θ_im, t-1)  ] 
        + ∑_i=1^m ∑_t=2^n [  (1-E_γγ_i, t) logp_spike(α_im, t) + E_γγ_i, tlogp_slab(α_im, t| α_im, t-1) ] + const..




§ DETAILS OF DMGP-GP

DMGP-GP is a state-of-art non-stationary MGP model, which constructs a LMC model for all outputs and assumes the hyper-parameters follow other GPs <cit.>:

    y(x_t)    = A_t q(x_t) + ϵ
    log(A_ii, t)    ∼𝒢𝒫(0, k_α(t, t^')) 
    
    A_ij, t   ∼𝒢𝒫(0, k_α(t, t^')), i ≠ j 
    
    q_i (x_t)    ∼𝒢𝒫(0, k(x_t, x_t^')),  
    
    
    k(x_t, x_t^')    = √(2 θ_tθ_t^'/θ_t^2 + θ_t'^2)exp[ (x_t - x_t^')^2/2(θ_t^2 + θ_t^'^2)] 
    log(θ_t)    ∼𝒢𝒫(0, k_θ(t, t^'))

where y(x_t)  = [y_1(x_t), ..., y_m(x_t)]^T are m outputs, A_t ∈ R^m × m is the time-varying coefficient matrix, q(x_t) = [q_1(x_t), ..., q_m(x_t)]^T are m i.i.d. latent Gaussian processes with zero mean and the same covariance function k(x_t, x_t^'), and ϵ=(ϵ_1,...ϵ_m) is measurement noise with ϵ_i ∼ N(0, σ_i^2). The covariance for the m outputs is

    cov[y(x_t),  y(x_t^')] = A_t A_t'^T k(x_t, x_t^') +  diag{σ_i},

where A_t A_t^'^T ∈ℝ^m × m is the correlation matrix of m outputs, and diag{σ_i} is the diagonal matrix with elements {σ_i}_i=1^m.



§ SCALABILITY OF DMGP-SS

Since our model supposes that the latent processes z_i(x) are independent on each other, the computational complexity is O(mn^3) when all the m outputs have an equal length of n. In comparison, the computational complexity of the classical MGP is O(m^3n^3), much larger than that of ours. Therefore, the proposed model can handle a number of outputs much easier. 
For example, we test our method on one numerical case with up to mn=8580 points. tab: scalable experiments show the prediction error and model fitting time. The prediction accuracy of the proposed method is better than that of GP in all the three experiments. Besides, although the third experiments have two times as many points as the first one, the fitting time of the third one (m=33, n=260) is only two times that of the first one (m=17, n=260), which identifies that the computational complexity of our method is O(mn^3). Besides, the fitting time of the second experiment (m=17, n=520) is only four times that of the first one (m=17, n=260), which means the second experiment takes less gradient descent steps to converge than the first one does. 





This work was supported by NSFC under Grants NSFC-72171003, NSFC-71932006.














69




[AlBahar et al.(2022)AlBahar, Kim, Wang,  Yue]BOdeepGP2022
AlBahar A, Kim I, Wang X, Yue X (2022) Physics-constrained bayesian optimization for optimal actuators placement in composite structures assembly. IEEE Transactions on Automation Science and Engineering .

[Alvarez  Lawrence(2011)]EfficientConvolvedMGP2011
Alvarez MA, Lawrence ND (2011) Computationally efficient convolved multiple output gaussian processes. The Journal of Machine Learning Research 12:1459–1500.

[Bai et al.(2022)Bai, Safikhani,  Michailidis]NonStaAR2022
Bai Y, Safikhani A, Michailidis G (2022) Hybrid modeling of regional covid-19 transmission dynamics in the u.s. IEEE Journal of Selected Topics in Signal Processing 16(2):261–275, <http://dx.doi.org/10.1109/JSTSP.2022.3140703>.

[Boyle  Frean(2004)]DependentGP2004
Boyle P, Frean M (2004) Dependent gaussian processes. Advances in neural information processing systems 17.

[Cheng et al.(2021)Cheng, Lu,  Peng]NonStaSTKNN2021
Cheng S, Lu F, Peng P (2021) Short-term traffic forecasting by mining the non-stationarity of spatiotemporal patterns. IEEE Transactions on Intelligent Transportation Systems 22(10):6365–6383, <http://dx.doi.org/10.1109/TITS.2020.2991781>.

[Choong et al.(2009)Choong, Charbit,  Yan]AutoregressiveImpute2009
Choong MK, Charbit M, Yan H (2009) Autoregressive-model-based missing value estimation for dna microarray time series data. IEEE Transactions on information technology in biomedicine 13(1):131–137.

[Christakos(2012)]STRandomFields2012
Christakos G (2012) Random field models in earth sciences (Courier Corporation).

[Damianou  Lawrence(2013)]DeepGP2013
Damianou A, Lawrence ND (2013) Deep gaussian processes. Artificial intelligence and statistics, 207–215 (PMLR).

[Dance  Paige(2022)]SSGP2022
Dance H, Paige B (2022) Fast and scalable spike and slab variable selection in high-dimensional gaussian processes. International Conference on Artificial Intelligence and Statistics, 7976–8002 (PMLR).

[Fairchild et al.(2017)Fairchild, Lake, Kattwinkel, Moorman, Bateman, Grieve, Isler,  Sahni]VitalSign2017
Fairchild KD, Lake DE, Kattwinkel J, Moorman JR, Bateman DA, Grieve PG, Isler JR, Sahni R (2017) Vital signs and their cross-correlation in sepsis and nec: a study of 1,065 very-low-birth-weight infants in two nicus. Pediatric research 81(2):315–321.

[Fothergill et al.(2012)Fothergill, Mentis, Kohli,  Nowozin]GestureData2012
Fothergill S, Mentis H, Kohli P, Nowozin S (2012) Instructing people for training gestural interactive systems. Proceedings of the SIGCHI conference on human factors in computing systems, 1737–1746.

[Frazier(2018)]BayesianOptimization2018
Frazier PI (2018) Bayesian optimization. Recent advances in optimization and modeling of contemporary problems, 255–278 (Informs).

[Fricker et al.(2013)Fricker, Oakley,  Urban]MGPemulator2013
Fricker TE, Oakley JE, Urban NM (2013) Multivariate gaussian process emulators with nonseparable covariance structures. Technometrics 55(1):47–56.

[Garg et al.(2012)Garg, Singh,  Ramos]NonStaSTGP2012
Garg S, Singh A, Ramos F (2012) Learning non-stationary space-time models for environmental monitoring. Proceedings of the AAAI Conference on Artificial Intelligence, volume 26, 288–294.

[Gelfand et al.(2004)Gelfand, Schmidt, Banerjee,  Sirmans]VaryingCoregionalization2004
Gelfand AE, Schmidt AM, Banerjee S, Sirmans C (2004) Nonstationary multivariate process modeling through spatially varying coregionalization. Test 13(2):263–312.

[George  McCulloch(1993)]SS1993
George EI, McCulloch RE (1993) Variable selection via gibbs sampling. Journal of the American Statistical Association 88(423):881–889.

[Goulard  Voltz(1992)]LMC1992
Goulard M, Voltz M (1992) Linear coregionalization model: tools for estimation and choice of cross-variogram matrix. Mathematical Geology 24(3):269–286.

[Gramacy(2020)]Surrogates2020
Gramacy RB (2020) Surrogates: Gaussian process modeling, design, and optimization for the applied sciences (CRC press).

[Gramacy  Lee(2008)]TreeGP2008
Gramacy RB, Lee HKH (2008) Bayesian treed gaussian process models with an application to computer modeling. Journal of the American Statistical Association 103(483):1119–1130.

[Heaton et al.(2017)Heaton, Christensen,  Terres]ClusterGP2017
Heaton MJ, Christensen WF, Terres MA (2017) Nonstationary gaussian process models using spatial hierarchical clustering from finite differences. Technometrics 59(1):93–101.

[Heinonen et al.(2016)Heinonen, Mannerström, Rousu, Kaski,  Lähdesmäki]Non-stationaryGP2016
Heinonen M, Mannerström H, Rousu J, Kaski S, Lähdesmäki H (2016) Non-stationary gaussian process regression with hamiltonian monte carlo. Artificial Intelligence and Statistics, 732–740 (PMLR).

[Hersbach(2000)]CRPS2000
Hersbach H (2000) Decomposition of the continuous ranked probability score for ensemble prediction systems. Weather and Forecasting 15(5):559–570.

[Hu  Wang(2021)]OnlineMGP
Hu Z, Wang C (2021) Nonlinear online multioutput gaussian process for multistream data informatics. IEEE transactions on industrial informatics 18(6):3885–3893.

[Huber et al.(2021)Huber, Koop,  Onorante]SparseTVP2021
Huber F, Koop G, Onorante L (2021) Inducing sparsity and shrinkage in time-varying parameter models. Journal of Business & Economic Statistics 39(3):669–683.

[Ishwaran  Rao(2005)]SS2005
Ishwaran H, Rao JS (2005) Spike and slab variable selection: frequentist and bayesian strategies. The Annals of Statistics 33(2):730–773.

[Kalli  Griffin(2014)]SparseTVP2014
Kalli M, Griffin JE (2014) Time-varying sparsity in dynamic regression models. Journal of Econometrics 178(2):779–793.

[Kingma  Ba(2015)]ADAM2014
Kingma DP, Ba JL (2015) Adam: A method for stochastic optimization. Proceedings of 3rd International Conference on Learning Representations.

[Ko  Kim(2022)]DeepGP2022
Ko J, Kim H (2022) Deep gaussian process models for integrating multifidelity experiments with nonstationary relationships. IISE Transactions 54(7):686–698.

[Kontar et al.(2018)Kontar, Zhou, Sankavaram, Du,  Zhang]RULMGP2018
Kontar R, Zhou S, Sankavaram C, Du X, Zhang Y (2018) Nonparametric modeling and prognosis of condition monitoring signals using multivariate gaussian convolution processes. Technometrics 60(4):484–496.

[Kuss  Rasmussen(2003)]GPRL
Kuss M, Rasmussen C (2003) Gaussian processes in reinforcement learning. Advances in neural information processing systems 16.

[Lee et al.(2023)Lee, Wang, Wu, Cai,  Yue]ActivePartitionGP2023
Lee C, Wang K, Wu J, Cai W, Yue X (2023) Partitioned active learning for heterogeneous systems. Journal of Computing and Information Science in Engineering 23(4):041009.

[Liu et al.(2020)Liu, Gong, Yang,  Chen]NonStaRNN2020
Liu Y, Gong C, Yang L, Chen Y (2020) Dstp-rnn: A dual-stage two-phase attention-based recurrent neural network for long-term and multivariate time series prediction. Expert Systems with Applications 143:113082.

[Liu et al.(2022)Liu, Wu, Wang,  Long]NonStaTransformers2022
Liu Y, Wu H, Wang J, Long M (2022) Non-stationary transformers: Exploring the stationarity in time series forecasting. Advances in Neural Information Processing Systems 35:9881–9893.

[Matthews et al.(2017)Matthews, van der Wilk, Nickson, Fujii, Boukouvalas, León-Villagrá, Ghahramani,  Hensman]GPflow2017
Matthews AGdG, van der Wilk M, Nickson T, Fujii K, Boukouvalas A, León-Villagrá P, Ghahramani Z, Hensman J (2017) GPflow: A Gaussian process library using TensorFlow. Journal of Machine Learning Research 18(40):1–6.

[Meng et al.(2021)Meng, Soper, Lee, Liu, Greene,  Ray]Non-stationaryMGP2021
Meng R, Soper B, Lee HK, Liu VX, Greene JD, Ray P (2021) Nonstationary multivariate gaussian processes for electronic health records. Journal of Biomedical Informatics 117:103698.

[Moore(1990)]MountainCar1990
Moore AW (1990) Efficient memory-based learning for robot control. Technical report, University of Cambridge, Computer Laboratory.

[Paciorek  Schervish(2003)]Non-stationaryGP2003
Paciorek C, Schervish M (2003) Nonstationary covariance functions for gaussian process regression. Advances in neural information processing systems 16.

[Padakandla et al.(2020)Padakandla, KJ,  Bhatnagar]NonStationaryRL2020
Padakandla S, KJ P, Bhatnagar S (2020) Reinforcement learning algorithm for non-stationary environments. Applied Intelligence 50:3590–3606.

[Pan  Yang(2009)]TransferSurvey2009
Pan SJ, Yang Q (2009) A survey on transfer learning. IEEE Transactions on knowledge and data engineering 22(10):1345–1359.

[Park(2022)]JumpGP2022
Park C (2022) Jump gaussian process model for estimating piecewise continuous regression functions. Journal of Machine Learning Research 23(278):1–37.

[Paun et al.(2023)Paun, Husmeier,  Torney]Non-stationaryGP2023
Paun I, Husmeier D, Torney CJ (2023) Stochastic variational inference for scalable non-stationary gaussian process regression. Statistics and Computing 33(2):44.

[Rangapuram et al.(2018)Rangapuram, Seeger, Gasthaus, Stella, Wang,  Januschowski]NonStaRNN2018
Rangapuram SS, Seeger MW, Gasthaus J, Stella L, Wang Y, Januschowski T (2018) Deep state space models for time series forecasting. Advances in neural information processing systems 31.

[Ročková  George(2014)]EMVS2014
Ročková V, George EI (2014) Emvs: The em approach to bayesian variable selection. Journal of the American Statistical Association 109(506):828–846.

[Rockova  McAlinn(2021)]DynamicSS2021
Rockova V, McAlinn K (2021) Dynamic variable selection with spike-and-slab process priors. Bayesian Analysis 16(1):233–269.

[Rodrigues et al.(2019)Rodrigues, Henrickson,  Pereira]MGPtraffic2019
Rodrigues F, Henrickson K, Pereira FC (2019) Multi-output gaussian processes for crowdsourced traffic data imputation. IEEE Transactions on Intelligent Transportation Systems 20(2):594–603, <http://dx.doi.org/10.1109/TITS.2018.2817879>.

[Scheipl et al.(2012)Scheipl, Fahrmeir,  Kneib]FunctionSelection2012
Scheipl F, Fahrmeir L, Kneib T (2012) Spike-and-slab priors for function selection in structured additive regression models. Journal of the American Statistical Association 107(500):1518–1532.

[Shand  Li(2017)]NonStaAugment2017
Shand L, Li B (2017) Modeling nonstationarity in space and time. Biometrics 73(3):759–768.

[Shen et al.(2023)Shen, Gnanasambandam, Wang,  Kong]MGPBO2023
Shen B, Gnanasambandam R, Wang R, Kong ZJ (2023) Multi-task gaussian process upper confidence bound for hyperparameter tuning and its application for simulation studies of additive manufacturing. IISE Transactions 55(5):496–508.

[Stathopoulos  Karlaftis(2003)]ARIMA2003
Stathopoulos A, Karlaftis MG (2003) A multivariate state space approach for urban traffic flow modeling and prediction. Transportation Research Part C: Emerging Technologies 11(2):121–135.

[Ton et al.(2018)Ton, Flaxman, Sejdinovic,  Bhatt]NonStaFourier2018
Ton JF, Flaxman S, Sejdinovic D, Bhatt S (2018) Spatial mapping with gaussian processes and nonstationary fourier features. Spatial statistics 28:59–78.

[Verstraeten et al.(2020)Verstraeten, Libin,  Nowé]LMCRL
Verstraeten T, Libin PJ, Nowé A (2020) Fleet control using coregionalized gaussian process policy iteration. ECAI, 1571–1578.

[Wang et al.(2020a)Wang, Hamelijnck, Damoulas,  Steel]NonStaRandomFields2020
Wang K, Hamelijnck O, Damoulas T, Steel M (2020a) Non-separable non-stationary random fields. International Conference on Machine Learning, 9887–9897 (PMLR).

[Wang et al.(2020b)Wang, Cao,  Philip]STDeep2020
Wang S, Cao J, Philip SY (2020b) Deep learning for spatio-temporal data mining: A survey. IEEE transactions on knowledge and data engineering 34(8):3681–3700.

[Wang et al.(2022)Wang, Wang, Song, Kirby,  Wu]RegularizedMGP2022
Wang X, Wang C, Song X, Kirby L, Wu J (2022) Regularized multi-output gaussian convolution process with domain adaptation. IEEE Transactions on Pattern Analysis and Machine Intelligence 45(5):6142–6156.

[Wang et al.(2019)Wang, Zhang, Zhu, Long, Wang,  Yu]NonStaLSTM2019
Wang Y, Zhang J, Zhu H, Long M, Wang J, Yu PS (2019) Memory in memory: A predictive neural network for learning higher-order non-stationarity from spatiotemporal dynamics. Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, 9154–9162.

[Wen et al.(2023)Wen, Zhou, Zhang, Chen, Ma, Yan,  Sun]TransformersSurvey2023
Wen Q, Zhou T, Zhang C, Chen W, Ma Z, Yan J, Sun L (2023) Transformers in time series: a survey. Proceedings of the Thirty-Second International Joint Conference on Artificial Intelligence, 6778–6786.

[Williams  Rasmussen(2006)]GP2006
Williams CK, Rasmussen CE (2006) Gaussian processes for machine learning, volume 2 (MIT press Cambridge, MA).

[Xia et al.(2016)Xia, Wang, Li, Li,  Zhang]STKNN2016
Xia D, Wang B, Li H, Li Y, Zhang Z (2016) A distributed spatial–temporal weighted model on mapreduce for short-term traffic flow forecasting. Neurocomputing 179:246–263.

[Xu et al.(2022)Xu, Wu, Yue,  Li]DynamicSubspace2022
Xu R, Wu J, Yue X, Li Y (2022) Online structural change-point detection of high-dimensional streaming data via dynamic sparse subspace learning. Technometrics 1–14.

[Yang et al.(2020)Yang, Xu, Wu,  Wang]MultiTaskRL2020
Yang R, Xu H, Wu Y, Wang X (2020) Multi-task reinforcement learning with soft modularization. Advances in Neural Information Processing Systems 33:4767–4777.

[Ye  Dai(2021)]TransferForecasting2021
Ye R, Dai Q (2021) Implementing transfer learning across different datasets for time series forecasting. Pattern Recognition 109:107617.

[Yoon  Li(2018)]PositiveTransfer2018
Yoon H, Li J (2018) A novel positive transfer learning approach for telemonitoring of parkinson’s disease. IEEE Transactions on Automation Science and Engineering 16(1):180–191.

[Yun et al.(2022)Yun, Zhang,  Li]STRandomFields2022
Yun S, Zhang X, Li B (2022) Detection of local differences in spatial characteristics between two spatiotemporal random fields. Journal of the American Statistical Association 117(537):291–306.

[Zhang et al.(2021)Zhang, Yan, Lee,  Shi]DynamicSubspace2021
Zhang C, Yan H, Lee S, Shi J (2021) Dynamic multivariate functional data modeling via sparse subspace learning. Technometrics 63(3):370–383.

[Zhang et al.(2016)Zhang, Wang,  Chen]QualityMonitoring2016
Zhang L, Wang K, Chen N (2016) Monitoring wafers’ geometric quality using an additive gaussian process model. IIE Transactions 48(1):1–15.

[Zhang  Yang(2021)]Multi-taskSurvey2021
Zhang Y, Yang Q (2021) A survey on multi-task learning. IEEE Transactions on Knowledge and Data Engineering 34(12):5586–5609.

[Zhu et al.(1997)Zhu, Byrd, Lu,  Nocedal]BFGS1997
Zhu C, Byrd RH, Lu P, Nocedal J (1997) Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale bound-constrained optimization. ACM Transactions on mathematical software (TOMS) 23(4):550–560.

[Zhu et al.(2023)Zhu, Lin, Jain,  Zhou]TransferRL2023
Zhu Z, Lin K, Jain AK, Zhou J (2023) Transfer learning in deep reinforcement learning: A survey. IEEE Transactions on Pattern Analysis and Machine Intelligence .

[Zou et al.(2023)Zou, Careem, Dutta,  Thawdar]NonStaRandomFields2023
Zou Z, Careem M, Dutta A, Thawdar N (2023) Joint spatio-temporal precoding for practical non-stationary wireless channels. IEEE Transactions on Communications 71(4):2396–2409.





