Nonstationary and Sparsely-Correlated Multioutput Gaussian Process with Spike-and-Slab Prior

1. Introduction

Gaussian process (GP) provides an elegant and flexible Bayesian nonparametric framework for modeling nonlinear mappings (Williams and Rasmussen 2006). 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 (Frazier 2018), experiment design (Gramacy 2020), and product quality monitoring (Zhang et al. 2016). However, standard GP is designed for only a single output, which limits its use in multioutput or multitask scenarios arising in various fields, such as Bayesian optimization (Shen et al. 2023), a traffic network (Rodrigues et al. 2019), and a computer simulation emulator (Fricker et al. 2013). Consequently, multioutput Gaussian process (MGP) has been gaining increasing attention from researchers and has emerged as an important member in the vast family of transfer learning (Pan and Yang 2009) and multitask learning (Zhang and Yang 2021) 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 nonstationary environments, where the data characteristics vary across the input domain (Williams and Rasmussen 2006). 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 (Paun et al. 2023). In clinical studies, sepsis is very likely to cause changes in the cross-correlation among vital signs in the early onset (Fairchild et al. 2017). In kinesiology, the cooperation patterns of human joints vary across different gestures: for example, both hands move jointly in a “shoot” action but separately in a “throw” action (Xu et al. 2022). In such cases, nonstationary models, which 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 providing medical care (Paun et al. 2023).

Mainly two kinds of methods have been proposed to capture the dynamic characteristics of the nonstationary 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 (Gramacy and Lee 2008) uses a tree structure to partition the input space of a computer simulation emulator. Another method, called jump GP, cuts the input space into several segments to model piece-wise continuous functions (Park 2022). This model is optimized using the expectation-minimization (EM) algorithm or variational inference. Further, a clustering-based GP (Heaton et al. 2017) 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 with heterogeneous systems (Lee et al. 2023). 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, nonstationary GP introduced in Heinonen et al. (2016) and Paun et al. (2023) applies GP priors on the amplitude and length-scale parameters of a square-exponential covariance function. In addition, based on these single-output nonstationary GPs, researchers have explored MGP models for multivariate data with dynamic characteristics. For example, a nonstationary MGP is established to model the varying correlation between vital signals, where a GP prior is imposed on the time-dependent correlation matrix (Meng et al. 2021). 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 negative information transfer, which is another critical challenge when using MGP. Transfer learning is very promising for leveraging information to a data-insufficient domain, called a target domain, from other correlated domains, called the 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 (Pan and Yang 2009, Yoon and Li 2018). 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 closely related one for information transfer. For example, the method proposed in Ye and Dai (2021) 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 (Wang et al. 2022), 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 nonstationary 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 the 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, whereas 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 relationships and negative transfer. Our contributions can be summarized as follows:

1. A novel nonstationary 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.

2. 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.

3. 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 nonzero elements.

The rest of this paper is organized as follows. In Section 2, we revisit the related literature and the static MGP model. Section 3 presents the proposed nonstationary MGP and an efficient EM algorithm for model training. In Section 4, we evaluate the effectiveness of our model on simulated data. In Section 5, we perform one time-series analysis case on human gesture data (Fothergill et al. 2012) and one control policy optimization case on the mountain-car problem (Moore 1990). In Section 6, we conclude the paper with a discussion.

2. Preliminaries

In this section, we first review research that is related to our work. We then introduce the static MGP based on the convolution process, which has been widely applied in various areas because of its flexibility (Boyle and Frean 2004, Kontar et al. 2018, Hu and Wang 2021).

2.1. Related Work

To deal with nonstationary 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., Gramacy and Lee 2008, Heaton et al. 2017, Park 2022, Lee et al. 2023), or allow them to vary at each point and model them using other GPs (the second category, e.g., Paciorek and Schervish 2003, Heinonen et al. 2016, Paun et al. 2023). 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 (Damianou and Lawrence 2013, AlBahar et al. 2022, Ko and Kim 2022). However, the parameters of Deep GP are stationary, which differs from the second category where the covariance parameters are dynamic.

Nonstationary 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 the 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 (Goulard and Voltz 1992). The existing nonstationary MGPs are mainly extensions of the classical LMC model. For example, the approach in Gelfand et al. (2004) allows the correlation matrix to vary with inputs to model the dynamic relationship among outputs. In Meng et al. (2021), a nonstationary 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 (Boyle and Frean 2004). However, these approaches are for stationary data. Furthermore, all existing approaches fail to capture a sparse correlation structure in a nonstationary environment.

Spatial-temporal modeling of nonstationary data is closely related to our work. In comparison with normal time-series modeling, spatial-temporal analysis requires modeling the spatial correlation to enhance the prediction accuracy. A large number of spatial-temporal models have been investigated, such as the spatial-temporal autoregressive integrated moving average method (ST-ARIMA) (Stathopoulos and Karlaftis 2003), spatial-temporal k-nearest neighbors (ST-KNN) (Xia et al. 2016), spatial-temporal random fields (Christakos 2012, Yun et al. 2022), and spatial-temporal deep neural networks (Wang et al. 2020b, Wen et al. 2023). Based on the aforementioned methods, a number of recent works try to extend them to handle nonstationary 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, for example, an ST-KNN with a wrapped K-means partition algorithm (Cheng et al. 2021) or an autoregressive model coupled with a block-fused-Lasso change-point detection algorithm (Bai et al. 2022). Besides partitioning the time-series into stationary parts, the method proposed by Shand and Li (2017) maps the nonstationary space-time process into a high-dimensional stationary process through augmenting new dimensions. Another type of nonstationary spatial-temporal model is Bayesian random fields with nonstationary space-time kernels (Garg et al. 2012, Ton et al. 2018, Wang et al. 2020a, Zou et al. 2023), whose hyper-parameters change over time or location. Deep learning methods are also explored on nonstationary data recently, such as nonstationary recurrent neural networks (Rangapuram et al. 2018, Liu et al. 2020), long short-term memory networks (Wang et al. 2019), and transformer-based networks (Liu et al. 2022, Wen et al. 2023). 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.

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 (George and McCulloch 1993) that has been used for feature selection in (generalized) linear models, additive models, and Gaussian processes (Ishwaran and Rao 2005, Scheipl et al. 2012, Dance and Paige 2022). With this prior, smaller parameters tend to be more influenced by the spike prior to reach zero, whereas 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 (Kalli and Griffin 2014, Huber et al. 2021, Rockova and McAlinn 2021), no work utilizes a dynamic spike-and-slab prior to model the dynamic and sparse correlation among outputs in a nonstationary MGP.

2.2. Static MGP Based on Convolution Process

Consider a set of m outputs $f_i:\mathcal{X}\mapsto \mathbb{R},\hspace{0.17em}i=1,\dots ,m$, where $\mathcal{X}$ 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\sim \mathcal{N}(0,{\sigma }_i^2)$, that is,

where $\mathit{x}\in {\mathbb{R}}^d$ is the input. Denote the n_i observed data for the ith output as ${\mathcal{D}}_i=\{{\mathit{X}}_i,{\mathit{y}}_i\}$, where ${\mathit{X}}_i={({\mathit{x}}_{i,1},\dots ,{\mathit{x}}_{i,n_i})}^T$ and ${\mathit{y}}_i={(y_{i,1},\dots ,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={\sum }_i{n}_i$ for m outputs. Denote the data of all outputs as $\mathcal{D}=\{\mathit{X},\mathit{y}\}$, where $\mathit{X}={({\mathit{X}}_1^T,\dots ,{\mathit{X}}_m^T)}^T\in {\mathbb{R}}^{N\times d}$ and $\mathit{y}={({\mathit{y}}_1^T,\dots ,{\mathit{y}}_m^T)}^T\in {\mathbb{R}}^N$. In an MGP model, the observation vector y follows a joint Gaussian distribution:where $\mathit{K}=\mathit{K}(\mathit{X},\mathit{X})\in {\mathbb{R}}^{N\times N}$ is a block-partitioned covariance matrix. The $(i,i^{\prime })$-th block of K, ${\mathit{K}}_{i,i^{\prime }}={\text{cov}}_{i,i^{\prime }}^f({\mathit{X}}_i,{\mathit{X}}_{i^{\prime }})+{\tau }_{i,i^{\prime }}{\sigma }_i^2\mathit{I}\in {\mathbb{R}}^{n_i\times n_{i^{\prime }}}$ represents the covariance matrix between the output i and output $i^{\prime }$ (${\tau }_{i,i^{\prime }}$ equals to one if $i=i^{\prime }$, and zero otherwise). The function ${\text{cov}}_{i,i^{\prime }}^f(\mathit{x},{\mathit{x}}^{\prime })$ measures the covariance between $f_i(\mathit{x})$ and $f_{i^{\prime }}({\mathit{x}}^{\prime })$. In the covariance matrix K, the cross-covariance block ${\mathit{K}}_{i,i^{\prime }}(i\ne i^{\prime })$ 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(\mathit{x})\}}_{j=1}^h$ and kernel functions ${\{g_{ji}(\mathit{x})\}}_{j=1}^h$ in the following way (Boyle and Frean 2004, Alvarez and Lawrence 2011):

where * represents the convolution operation, α_ji is the amplitude parameter, and h is the number of shared latent processes. Usually, ${\{z_j(\mathit{x})\}}_{j=1}^h$ are independent white Gaussian noise processes with $\text{cov}(z(\mathit{x}),z({\mathit{x}}^{\prime }))=\delta (\mathit{x}-{\mathit{x}}^{\prime })$, where $\delta (·)$ is the Dirac delta function. Thus, the covariance function can be derived aswhere $\mathit{v}=\mathit{x}-{\mathit{x}}^{\prime }$. In such a way, the covariance between $f_i(\mathit{x})$ and $f_{i^{\prime }}({\mathit{x}}^{\prime })$ is dependent on their difference $\mathit{x}-{\mathit{x}}^{\prime }$, the amplitude parameters, and the hyper-parameters in kernels g_ji and $g_{j{i}^{\prime }}$. Compared with the classical LMC model $f_i(\mathit{x})={\sum }_{j=1}^h{\alpha }_{ji}{q}_j(\mathit{x})$, where $q_j(\mathit{x})$ is a latent GP with covariance $k_j(\mathit{x},{\mathit{x}}^{\prime })$. The convolution-process-based MGP is more flexible than LMC, as it does not restrict all outputs to having the same autocovariance pattern.

At a new point ${\mathit{x}}_{*}$, the posterior distribution of $y_i({\mathit{x}}_{*})$ given data $\{\mathit{X},\mathit{y}\}$ is

where the predictive mean $\mu ({\mathit{x}}_{*})$ and variance $\mathrm{\Sigma }({\mathit{x}}_{*})$ can be expressed aswhere ${\mathit{K}}_{*}^T=[{\text{cov}}_{i1}^f{({\mathit{x}}_{*},{\mathit{X}}_1)}^T,\dots ,{\text{cov}}_{im}^f{({\mathit{x}}_{*},{\mathit{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 nonstationary MGP to simultaneously address both problems.

3. Model Development

We propose a nonstationary 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.

3.1. The Proposed Model

The structure of our hierarchical model is illustrated in Figure 1. 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.

3.1.1. Dynamic MGP.

In this subsection, we introduce the major part of the proposed nonstationary MGP, which corresponds to the first layer in Figure 1. To simplify the notation, we slightly abuse the notation used in the previous section. Specifically, we take the first m − 1 outputs $f_i:\mathcal{X}\mapsto \mathbb{R},\hspace{0.17em}i=1,\dots ,m-1$ as the sources, and the last one $f_m:\mathcal{X}\mapsto \mathbb{R}$ as the target. We still assume that the observation y_i is accompanied with the i.i.d. measurement noise ${ϵ}_i\sim N(0,{\sigma }_i^2)$. Let $\mathcal{I}=\{1,2,\dots ,m\}$ be the index set of all outputs, and ${\mathcal{I}}^S=\mathcal{I}/\{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\in \{1,2,\dots ,n\}$, that is, ${\mathit{X}}_i={({\mathit{x}}_{i,1},\dots ,{\mathit{x}}_{i,t},\dots ,{\mathit{x}}_{i,n})}^T$ and ${\mathit{y}}_i={(y_{i,1},\dots ,y_{i,t},\mathrm{.}.,y_{i,n})}^T$. In Online 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

where ${\alpha }_{ii,t}$ and ${\alpha }_{jm,t}$ are time-varying amplitude parameters, $g_{ii,t}(\mathit{x})$ and $g_{jm,t}(\mathit{x})$ are time-varying kernel functions, and ${\{z_i(\mathit{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 because of its flexibility (Boyle and Frean 2004, Alvarez and Lawrence 2011). The kernel is given bywhere ${\mathit{\theta }}_{ij,t}$ is a diagonal matrix representing the length-scale for each input feature. More importantly, such a Gaussian kernel can yield closed-form covariance functions through the convolution operation in Equation (3) (Paciorek and Schervish 2003, Wang et al. 2020a).

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, ${\{{\alpha }_{ii,t},{\mathit{\theta }}_{ii,t}\}}_{i=1}^m$ are responsible for the nonstationary behavior within each output, whereas ${\{{\alpha }_{im,t},{\mathit{\theta }}_{im,t}\}}_{i=1}^{m-1}$ capture the dynamic correlation between target and sources. More specifically, the amplitude parameter ${\alpha }_{im,t}$ controls the knowledge transfer. For example, if ${\alpha }_{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 repartitioned as

where the $(i,i^{\prime })$-th block ${\mathit{K}}_{i,i^{\prime }}={\text{cov}}_{i,i^{\prime }}^f({\mathit{X}}_i,{\mathit{X}}_{i^{\prime }})+{\tau }_{i,i^{\prime }}{\sigma }_i^2\mathit{I}$, the block-diagonal matrix ${\mathit{K}}_{(ss)}$ represents the covariance of source outputs, and ${\mathit{K}}_{(sm)}$ represents the cross-covariance between the sources and the target. Based on Equation (3), we can obtain covariance functions for the proposed nonstationary MGP model, as shown below:where $i\in {\mathcal{I}}^S$. Equations (10a)–(10c) 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 (Heinonen et al. 2016): ${\alpha }_{ij,t}=\log [1+\exp ({\tilde{\alpha }}_{ij,t})],{\mathit{\theta }}_{ij,t}=\log [1+\exp ({\tilde{\mathit{\theta }}}_{ij,t})],$ where ${\tilde{\alpha }}_{ij,t},{\tilde{\mathit{\theta }}}_{ij,t}$ are underlying parameters to estimate, whose range is $[-\infty ,\infty ]$. The proposed covariance functions can be viewed as an extension of the nonstationary kernels (Paciorek and Schervish 2003) from the single-output to the multioutput case. Specifically, the autocovariance of each source in Equation (10a) is the same as the covariance for a single-output nonstationary GP. From the cross-covariance between each source and the target, we can clearly see that the amplitude parameter ${\alpha }_{im,t}$ controls whether the cross-correlation is zero or not. The validity of the proposed covariance functions is outlined in Proposition 1:

Proposition 1.

The proposed nonstationary MGP covariance matrix in Equation (9) is positive-definite, that is, $\forall \mathit{y}\ne \mathbf{0}$,

$${\mathit{y}}^T\mathit{\text{Ky}}>0.$$

The proof is provided in Online Appendix B.

Based on Equation (9), the joint distribution of all sources and the target is expressed as

where ${\mathit{y}}_{(s)}$ is the collection of all source data. For notational convenience, we partition the parameters in a similar way, ${\mathit{\alpha }}_{(s)}={\{{\mathit{\alpha }}_{(s),t}\}}_{t=1}^n,{\mathit{\theta }}_{(s)}={\{{\mathit{\theta }}_{(s),t}\}}_{t=1}^n,{\mathit{\sigma }}_{(s)}$ $={\{{\sigma }_i\}}_{i=1}^{m-1},{\mathit{\alpha }}_m={\{{\mathit{\alpha }}_{m,t}\}}_{t=1}^n,{\mathit{\theta }}_m={\{{\mathit{\theta }}_{m,t}\}}_{t=1}^n,$ where ${\mathit{\alpha }}_{(s),t}={\{{\alpha }_{ii,t}\}}_{i=1}^{m-1},\hspace{0.17em}{\mathit{\theta }}_{(s),t}={\{{\theta }_{ii,t}\}}_{i=1}^{m-1},\hspace{0.17em}{\mathit{\alpha }}_{m,t}={\{{\alpha }_{im,t}\}}_{i=1}^m$, and ${\mathit{\theta }}_{m,t}={\{{\theta }_{im,t}\}}_{i=1}^m$. Furthermore, we denote the collection of all parameters as $\mathbf{\Phi }=\{{\mathbf{\Phi }}_{(s)},{\mathbf{\Phi }}_m\}$, where ${\mathbf{\Phi }}_{(s)}=\{{\mathit{\alpha }}_{(s)},{\mathit{\theta }}_{(s)},{\mathit{\sigma }}_{(s)}\}$ and ${\mathbf{\Phi }}_m=\{{\mathit{\alpha }}_m,{\mathit{\theta }}_m,{\sigma }_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 overfitting problem. To address this issue, Gaussian processes are typically employed to model the kernel parameters, for example, $\log ({\alpha }_t)\sim \mathcal{G}\mathcal{P}(0,k_{\alpha }(t,t^{\prime })),\log ({\mathit{\theta }}_t)\sim \mathcal{G}\mathcal{P}(0,k_{\mathit{\theta }}(t,t^{\prime }))$, where $k_{\alpha },k_{\mathit{\theta }}$ are covariance functions for the amplitude and length-scale parameters, respectively. Although this technique can force the parameters to vary smoothly and reduce overfitting, it cannot model a sparse correlation between the sources and the target. Consequently, this approach cannot avoid the negative transfer caused by unrelated sources.

3.1.2. Dynamic Spike-and-Slab Prior.

The classical spike-and-slab prior (George and McCulloch 1993, Dance and Paige 2022) only handles the shrinkage of one single parameter and cannot model the smooth-varying parameters. To take both the dynamic and sparse properties of correlation into account, we propose a dynamic spike-and-slab prior placing on ${\mathit{\alpha }}_m$:

where ${\gamma }_{i,t}\in \{0,1\}$ is a binary sparse indicator for ${\alpha }_{im,t}$ following a Bernoulli distribution, $p_{\mathit{\text{spike}}}({\alpha }_{im,t})$ is a zero-mean spike prior pushing the parameter to zero, $p_{\mathit{\text{slab}}}({\alpha }_{t,im}|{\alpha }_{t-1,im})$ is a slab prior connecting ${\alpha }_{im,t-1}$ and ${\alpha }_{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 Figure 1. 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, that is,

Note that ${\mathit{\theta }}_{ij,t}$ is a diagonal matrix, so that the slab prior is placed on its d diagonal elements independently, that is, $p_{\mathit{\text{slab}}}({\mathit{\theta }}_{ii,t}|{\mathit{\theta }}_{ii,t-1})={\prod }_{l=1}^d{p}_{\mathit{\text{slab}}}({\{{\mathit{\theta }}_{ii,t}\}}_l|{\{{\mathit{\theta }}_{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, for example, from ${\alpha }_{ii,t-1}$ to ${\alpha }_{ii,t}$. Compared with the dynamic spike-and-slab prior used in linear models (Rockova and McAlinn 2021), our method does not constrain the slab prior to be a stable autoregressive process. Besides, we use a simpler but more flexible prior for ${\gamma }_{i,t}$, whereas the work in Rockova and McAlinn (2021) 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, that is,

where ν₀ is the length-scale for Laplace distribution. If we maximize the log-likelihood function to optimize the model, this prior will become an L₁ 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,which encourages ${\alpha }_{im,t}$ to remain constant in a continuous period, approximating a piece-wise model. In the second one, the parameters are allowed to change smoothly,where ν₁ is variance of Gaussian distribution, and $\rho <1$ is an autoregressive coefficient. A similar smoothing approach can also be found in Hu and Wang (2021). These two slab priors make the current parameter value exactly or roughly concentrated around the previous value. Typically, we set ν₀ to be much smaller than $\sqrt{{\nu }_1}$ in the soft slab prior (or ν₁ in the hard slab prior) to make the two priors more separable and to put more penalty on sparsity. Besides, the values of ${\mathit{\alpha }}_{im}$ at multiple time steps before t can be included in the soft slab prior, for example, $p_{\mathit{\text{slab}}}^{\mathit{\text{soft}}}({\alpha }_{im,t}|{\alpha }_{im,t-1},{\alpha }_{im,t-2},\dots )$. We choose the simplest form $p_{\mathit{\text{slab}}}^{\mathit{\text{soft}}}({\alpha }_{im,t}|{\alpha }_{im,t-1})$ because of its wide application and robust performance in practice.

At time t, η can be interpreted as the prior probability that ${\alpha }_{im,t}$ belongs to the slab process. It influences the strength of the shrinkage effect that ${\alpha }_{im,t}$ bears. Ideally, for nonzero ${\alpha }_{im,t}$, the posterior mean of ${\gamma }_{i,t}$ should be close to one so that ${\alpha }_{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 ${\gamma }_{i,t}$ is automatically adjusted based on the estimated ${\alpha }_{im,t}$ to avoid shrinking nonzero elements. This makes our method superior to traditional Lasso methods where the sparse penalty weights are identical for zero and nonzero parameters.

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

3.2. 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 nonstationary model with considerable time-varying parameters. Therefore, we develop an efficient EM algorithm. Instead of directly maximizing the posterior $p(\mathbf{\Phi }|\mathit{y})=p({\mathbf{\Phi }}_{(s)},{\mathbf{\Phi }}_m|{\mathit{y}}_{(s)},{\mathit{y}}_m)$, we proceed iteratively in terms of the complete log posterior $\log p(\mathbf{\Phi },\mathit{\gamma }|\mathit{y})$, where the binary parameters $\mathit{\gamma }$ are treated as “missing data.” Because 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 $\mathbf{\Phi }$. More precisely, the E-step and M-step at the $(k+1)$th iteration can be expressed as

where $E_{\mathit{\gamma }|{\mathbf{\Phi }}^{(k)},\mathit{y}}(·)$ is the conditional expectation on posterior of $\mathit{\gamma }$, and ${\mathbf{\Phi }}^{(k)}$ is the optimized parameters at the kth iteration. For simplicity, we use $E_{\mathit{\gamma }}(·)$ to denote $E_{\mathit{\gamma }|{\mathbf{\Phi }}^{(k)},\mathit{y}}(·)$.

Based on Bayes’ theorem and the property of multivariate normal distribution, the expectation of $\log p(\mathbf{\Phi },\mathit{\gamma }|\mathit{y})$ can be as follows (derivation details can be found in Online Appendix C):

where $\mathit{\mu }={\mathit{K}}_{(sm)}^T{\mathit{K}}_{(ss)}^{-1}{\mathit{y}}_{(s)}$ is the conditional mean of target given the sources and $\mathbf{\Sigma }={\mathit{K}}_{mm}-{\mathit{K}}_{(sm)}^T{\mathit{K}}_{(ss)}^{-1}{\mathit{K}}_{(sm)}$ is the conditional covariance.

In the E-step, because $\mathit{\gamma }$ is only dependent on ${\mathbf{\Phi }}_m$, the posterior of ${\gamma }_{i,t}$ is calculated as

Then the conditional expectation of ${\gamma }_{i,t}$ can be updated as

The posterior mean $E_{\mathit{\gamma }}{\gamma }_{i,t}$ can be interpreted as the posterior probability of classifying the current parameter ${\alpha }_{im,t}$ into a slab process as opposed to a spike process, based on the past value ${\alpha }_{im,t-1}$. For example, we set η to 0.5 as a noninformative prior, and take a small ν₀ (e.g., 0.01) for the spike prior and a large ν₁ (e.g., 0.1) for the soft slab prior. Based on Equation (21), supposing that ${({\alpha }_{im,t}-{\alpha }_{im,t-1})}^2$ is small and $|{\alpha }_{im,t}|$ is large, the expectation $E_{\mathit{\gamma }}{\gamma }_{i,t}$ will tend toward one, indicating that ${\alpha }_{im,t}$ is more likely from the slab prior. On the other hand, if $|{\alpha }_{im,t}|$ is small, $E_{\mathit{\gamma }}{\gamma }_{i,t}$ is close to zero, enforcing strong shrinkage on ${\alpha }_{im,t}$.

In the M-step, we can optimize the objective function in Equation (19) with various gradient ascent methods, such as (stochastic) ADAM, L-BFGS, etc. (Zhu et al. 1997, Kingma and Ba 2015). 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 $\mathit{\gamma }$. Ideally, for nonzero ${\alpha }_{im,t},\hspace{0.17em}{E}_{\gamma }{\gamma }_{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 Equation (21), and this adjustment has explicit statistical interpretability.

In our case studies, we find the algorithm converges rapidly, for example, achieving convergence after only five iterations. The whole algorithm is summarized in Algorithm 1, where an ADAM method (Kingma and Ba 2015) is utilized in the M-step. Note that the large parameter space of $\mathbf{\Phi }$ poses a challenge in identifying the modes of the full posterior. To speed up the convergence, we propose to initialize the source parameter ${\mathbf{\Phi }}_{(s)}$ by maximizing the sum of sources’ marginal log-likelihood and source parameter prior:

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

3.3. 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 Equation (8) and derive closed-form covariance functions in Equations (10a)–(10c).

The second challenge is calculating the inverse of the covariance matrix, which is a critical issue for all GPs. The computational complexity of Algorithm 1 is approximately $O({\sum }_{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({({\sum }_{i=1}^m{n}_i)}^3)$ with the same size of data and nonseparable covariance. The main computational load of our model is in calculating the inverse of ${\mathit{K}}_{(ss)}$ and $\mathbf{\Sigma }$ in the source marginal distribution $\log p({\mathit{y}}_{(s)}|{\mathbf{\Phi }}_{(s)})$ and the target conditional distribution $\log p({\mathit{y}}_m|{\mathbf{\Phi }}_m,{\mathit{y}}_{(s)},{\mathbf{\Phi }}_{(s)})$, respectively, in Equation (19). Because all the latent processes are independent of each other, ${\mathit{K}}_{(ss)}$ is a block-diagonal matrix and the complexity is reduced from $O({\sum }_{i=1}^{m-1}{n}_i^3)$ to $O({\sum }_{i=1}^{m-1}{n}_i^3)$. The calculation of the inverse of $\mathbf{\Sigma }$ is $O(n_m^3)$. Therefore, the overall computational complexity is reduced to $O({\sum }_{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 nonstationary linear model (Rockova and McAlinn 2021), the MCMC and EM algorithm lead to very close prediction errors, but the running time of MCMC is about 10 times longer than that of the EM.

Algorithm 1

(The Optimization Algorithm for the Nonstationary MGP Model)

Input: Data ${\{{\mathit{X}}_i,{\mathit{y}}_i\}}_{i=1}^m$, ν₀, ν₁, ρ, η.

• 1: Set starting value: ${\mathbf{\Phi }}_{(s)}^{(-1)},{\mathbf{\Phi }}_m^{(0)}$.

• 2: Initialization: obtain ${\mathbf{\Phi }}_{(s)}^{(0)}$ through Equation (22).

• 3: for k_out iterations do

• 4:  E-step given ${\mathbf{\Phi }}^{(k)}$: Update ${\{E_{\mathit{\gamma }}{\gamma }_{i,t}\}}_{i=1,t=2}^{m,n}$ through Equation (21).

• 5:  M-step given $E_{\mathit{\gamma }}{\gamma }_{i,t}$:

• 6:  for k_in epochs do

• 7:   Calculate ${\mathit{K}}_{(ss)},\hspace{0.17em}{\mathit{K}}_{(sm)}$, and ${\mathit{K}}_{mm}$ according to Equation (10).

• 8:   Calculate the value and gradient of objective function in Equation (19).

• 9:   Obtain ${\mathbf{\Phi }}^{(k+1)}$ using the ADAM method.

• 10:  end for

• 11: end for

3.4. Tuning Parameter Selection

The tuning parameters for our model are ν₀, ν₁ (for the hard slab prior), and ν₂ (for the soft slab prior). Here, because the key of our method is selecting the most informative sources, we propose to maximize the following criterion:

where N is the number of data, $\log p(\mathit{y}|\mathit{X})$ is the log-likelihood for both the source and the target outputs, and $c_{\mathit{\nu }}({\mathit{\alpha }}_m)$ is the number of nonzero elements in ${\mathit{\alpha }}_m$ given $\mathit{\nu }$. A similar criterion is proposed in Zhang et al. (2021). Note that $\mathit{\nu }=\{{\nu }_0,{\nu }_1\}$ for the hard slab prior and $\mathit{\nu }=\{{\nu }_0,{\nu }_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, whereas the second term prefers simpler models where fewer sources are chosen to transfer information to the source. To reduce the computation, we first determine the ratios $r_1={\nu }_1/{\nu }_0$ and $r_2={\nu }_2/{\nu }_0$ to make the spike-and-slab priors separable (Rockova and McAlinn 2021). Then we design a two-dimensional search grid for $({\nu }_1,{\nu }_1/r_1)$ with the hard slab or $({\nu }_2,{\nu }_2/r_1)$ with the soft slab. The optimal value of $\mathit{\nu }$ is searched over the two-dimensional grid. For example, we set $r_1\in \{5,10\}$ and $({\nu }_1,{\nu }_0)\in \{(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.

3.5. Model Prediction

Because the EM algorithm only estimates the value of parameters at the observed time points, given a new input ${\mathit{x}}_{t^{*}}$ of interest, we first need to estimate the target parameter ${\mathit{\alpha }}_{m,t^{*}}$ and ${\mathit{\theta }}_{m,t^{*}}$ at the new time point $t^{*}$, then derive the predictive distribution of $y_m({\mathit{x}}_{t^{*}})$.

3.5.1. Forecasting.

In the forecasting task, $t^{*}>n$. The estimated value of ${\mathit{\theta }}_{im,t^{*}}$ is

which is actually the mode of $p_{\mathit{\text{slab}}}({\mathit{\theta }}_{im,t^{*}}|{\mathit{\theta }}_{im,n})$. As for ${\alpha }_{im,t^{*}}$, if $E_{\mathit{\gamma }}{\gamma }_{i,n}\ge 0.5$, we consider it from the slab process and estimate it using the same method as ${\mathit{\theta }}_{im,t^{*}}$. Otherwise, it is classified to a spike process and shrunk to zero (Ročková and George 2014).

3.5.2. Recovery.

In the recovery task, the target data are unobserved at some specific time points and we aim to recover the missing data, that is, $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, that is,

As the parameters before and after $t_{*}$ are already optimized by the EM algorithm, the estimation of ${\mathit{\theta }}_{im,t^{*}}$ becomes

where $\mathit{\text{LSE}}(·)$ represents a least-square estimation for autoregressive process introduced in Choong et al. (2009). We also let ${\alpha }_{im,t^{*}}=0$, if $E_{\mathit{\gamma }}{\gamma }_{i,t_{\mathit{\text{near}}}}<0.5$. Otherwise, we estimate its value in the same way as ${\mathit{\theta }}_{im,t^{*}}$.

Then, given the parameter ${\mathit{\alpha }}_{m,t^{*}}$ and ${\mathit{\theta }}_{m,t^{*}}$, and the new input point ${\mathit{x}}_{t^{*}}$, the joint distribution of $y_m({\mathit{x}}_{t^{*}})$ and observations ${\mathit{y}}_m$ can be expressed as

where ${\mu }_{t^{*}}={\mathit{K}}_{(s*)}^T{\mathit{K}}_{(ss)}^{-1}{\mathit{y}}_{(s)},\hspace{0.17em}{\mathbf{\Sigma }}_{t^{*}}={\mathit{K}}_{m*}-{\mathit{K}}_{(sm)}^T{\mathit{K}}_{(ss)}^{-1}{\mathit{K}}_{(s*)}$, and ${\mathbf{\Sigma }}_{t^{*},t^{*}}={\mathit{K}}_{t^{*},t^{*}}-{\mathit{K}}_{(s*)}^T{\mathit{K}}_{(ss)}^{-1}{\mathit{K}}_{(s*)}$. In the above equations, ${\mathit{K}}_{(s*)}$ is the cross-covariance matrix of the sources and the new input ${\mathit{x}}_{t^{*}},\hspace{0.17em}{\mathit{K}}_{(m*)}$ is the covariance of the target observation and the new point, and ${\mathit{K}}_{t^{*},t^{*}}$ is the variance at ${\mathit{x}}_{t^{*}}$. Then, the posterior distribution of $y_m({\mathit{x}}_{t^{*}})$ can be derived as

4. Numerical Study

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

4.1. General Settings

Similar to the assumption made in the model development section, we generate m sequences consisting of m − 1 sources and one target, sampled at the same time stamps. 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:

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

2. MGP-L1. It is a state-of-art static method introduced in Wang et al. (2022). MGP-L1 models the target and sources in one MGP model, with the same covariance structure as in Equation (9). The scaling parameters ${\{{\alpha }_{im}\}}_{i=1}^{m-1}$ are penalized by L₁ term to achieve source selection. The regularized log-likelihood of this model is

   $$\log \mathcal{F}=-\frac{1}{2}{\mathit{y}}^T{\mathit{K}}^{-1}\mathit{y}-\frac{1}{2}\log |\mathit{K}|-\lambda {\displaystyle \sum _{i=1}^{m-1}|}{\alpha }_{im}|-\mathit{\text{const}}.,$$
   where K is calculated using the static covariance functions in Wang et al. (2022) and λ is a tuning parameter.

3. DMGP-GP. This is a state-of-the-art nonstationary MGP model, which constructs an LMC model for all outputs and assumes the hyper-parameters follow other GPs (Meng et al. 2021). More details can be found in Online Appendix D. In this model, the covariance for the m outputs is

   $$\text{cov}[\mathit{y}(x_t),\mathit{y}(x_{t^{\prime }})]={\mathit{A}}_t{\mathit{A}}_{t\prime }^{T}k(x_t,x_{t^{\prime }})+\text{diag}\left\{{\sigma }_i\right\},$$
   where ${\mathit{A}}_t{\mathit{A}}_{t^{\prime }}^T\in {\mathbb{R}}^{m\times m}$ is the correlation matrix of m outputs, and $\text{diag}\{{\sigma }_i\}$ is the diagonal matrix with ${\{{\sigma }_i\}}_{i=1}^m$. In this study, we focus on the correlation between the sources and target, which corresponds to the last column of ${\mathit{A}}_t{\mathit{A}}_t^T$ (except the last element ${({\mathit{A}}_t{\mathit{A}}_t^T)}_{mm}$), that is, ${({\mathit{A}}_t{\mathit{A}}_t^T)}_{0:m-1,m}$.

All methods are implemented in Python 3.8 and executed on a computer with an Intel® Core™ i5-7400 CPU with 3.00 GHz and 16 GB RAM. Both GP and MGP-L1 are implemented using GPflow (Matthews et al. 2017). DMGP-GP is implemented using TensorFlow and optimized with ADAM (Kingma and Ba 2015), 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 1 to 5 and 400, respectively.

For MGP-L1, the weight of L₁ penalty λ is a tuning parameter. For DMGP-GP, we use square-exponential functions for $k_{\alpha }$ and $k_{\theta }$. For simplicity, we apply the same tuning parameters for both kernels, that is, the amplitude ${\alpha }^{\#}$ and length-scale ${\theta }^{\#}$. 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.

4.2. Simulation Cases

The main objective of this section is to demonstrate the effectiveness of our method in capturing the nonstationary and sparse cross-correlation between the sources and the target. For simplicity, we hold the other characteristics constant over time, for example, 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, whereas 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\sim \mathcal{N}(0,{0.3}^2)$.

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

  $$\begin{array}{l}{Y}_1(x_t)=3\sin (\pi x_t/20+e_1)+{ϵ}_t,\hspace{1em}{Y}_2(x_t)=2\sin (2\pi x_t/20+e_2)\exp [0.5(x_t\%40-1)]+{ϵ}_t,\\ Y_3(x_t)=3\sin (4\pi x_t/20+e_3)+{ϵ}_t,\hspace{1em}{Y}_4(x_t)=2\sin (5\pi x_t/20+e_4)+{ϵ}_t,\end{array}$$
  where $e_i\sim \mathcal{N}(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₂ from the standard sine function. In each experiment, we generate 4k sources through sampling each kind of function k times, that is, $m=4k+1$. Specifically, the sources $\{{\mathit{y}}_{i+4j}|0\le j\le k-1\}$ are sampled from Y_i. Then, we define the dynamic target output as
  $$f_m(x_t)=a_{1,t}\sin (\pi x_t/20)+a_{2,t}\sin (2\pi x_t/20)+a_{3,t}\sin (4\pi x_t/20)+a_{4,t}\sin (5\pi x_t/20).$$

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

  $$\begin{array}{l}{a}_{1,t}=(2+2a_1)I_{t<40},\hspace{1em}{a}_{2,t}=(2+2a_2)I_{40\le t<80}+(1+a_2)I_{80\le t\le 130},\\ a_{3,t}=(1+a_3)I_{80\le t\le 130},\hspace{1em}\\ a_{4,t}=0.\end{array}$$

  Therefore, there are three segments, $[0,40),\hspace{0.17em}[40,80)$, and $[80,130]$. Only the first, the second, and the second and third sources are correlated with the target in the three periods, respectively. The other sources remain uncorrelated with 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,

  $$\begin{array}{l}{a}_{1,t}=[(2+a_1)\cos (\pi t/120)+0.5]I_{t<40},\\ a_{2,t}=[(2+a_2)\sin (\pi t/120-\pi /6)+0.5]I_{40\le t<130},\\ a_{3,t}=[(2+a_3)\sin (\pi t/120-\pi /2)+0.5]I_{80\le t<130}.\end{array}$$

In all cases, we set k = 1, 4 to generate 4 and 16 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],\hspace{0.17em}[50,70]$, and $[90,110]$ respectively. These 30 data points are treated as missing data, whereas the others are used as training data.

4.3. 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. Figure 2 illustrates the estimated ${\{{\alpha }_{im}\}}_{i=1}^4$ for MGP-L1, the ${\{{\alpha }_{im,t}\}}_{i=1}^4$ for DMGP-SS, and the estimated ${({\mathit{A}}_t{\mathit{A}}_t^T)}_{1:4,m}$ for DMGP-GP in Cases 1 and 2 with four sources. Figure 3 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 Cases 1 and 2, respectively. Note that the values of $a_{i,t}$ and ${\alpha }_{im,t}$ are not identical, because $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 Figure 2. Firstly, DMGP-SS tracks closely the periods of correlation between each source and the target, achieving a dynamic selection of sources. Because the target’s characteristics do not change abruptly (as shown in Figure 3), it is reasonable that the estimated correlation change points are about 10 time steps before or after the designed change times. Second, the hard slab prior encourages nearly piece-wise correlation, whereas the soft slab prior allows smoothly changing correlation. Because of the appropriate selection of sources at different times in Figure 3, 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 because 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 ${\mathbf{\Sigma }}_{t,*}^T{\mathbf{\Sigma }}^{-1}{\mathbf{\Sigma }}_{t,*}$ in posterior prediction Equation (25).

In contrast, another nonstationary method, DMGP-GP, fails to estimate a sparse structure in the source-target correlation because 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 ${\mathit{\alpha }}_m$, the estimated values of nonzero parameters are constant over time and cannot reflect the change of correlation. As a result, it performs even worse than GP in recovering the target output in Figure 3.

Another advantage of DMGP-SS is the adaptive adjustment of the spike-and-slab combination weight, $E_{\mathit{\gamma }}{\gamma }_{i,t}$. In our settings, ${\nu }_0^{-1}$ in the spike prior is much larger than ${\nu }_1^{-1}$ in the hard slab prior (or $\sqrt{{\nu }_1^{-1}}$ in the soft slab prior), to put more penalty on the correlation sparsity. For example, we set ${\nu }_0=0.02$ and ${\nu }_1=0.1$ in Case 1. However, this sparse penalty does not cause significant bias on nonzero ${\alpha }_{im,t}$ because of the automatically adjusted $E_{\mathit{\gamma }}{\gamma }_{i,t}$ in the EM algorithm. Specifically, we initialize it with 0.99 to barely shrink parameters at the beginning. Then $E_{\mathit{\gamma }}{\gamma }_{i,t}$ is updated in the E-step of the EM algorithm. Figure 4 shows its estimated value after five iterations. For the correlated sources (e.g., the first source during $t\in [0,50]$), their corresponding $E_{\mathit{\gamma }}{\gamma }_{i,t}$ is very close to one, so they bear a negligible shrinkage effect from the spike prior. In contrast, for the uncorrelated sources (e.g., the first source after t = 50), $E_{\mathit{\gamma }}{\gamma }_{i,t}$ is approximately 0.2, implying a substantial shrinkage effect. The value 0.2 can be derived based on Equation (21). For consecutive zero elements, $p_{\mathit{\text{slab}}}=\eta {(2{\nu }_1)}^{-1}$, and $p_{\mathit{\text{spike}}}=(1-\eta ){(2{\nu }_0)}^{-1}$, resulting in $E_{\mathit{\gamma }}{\gamma }_{i,t}\approx {\nu }_1^{-1}/{\nu }_0^{-1}$.

Table 1 summarizes the results of 40 repetitions for each case. DMGP-SS outperforms all the other methods in both cases. Further, the increase 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. 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 Figure 3). DMGP-GP has a better performance on target prediction than GP and MGP-L1, because of 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^{\prime }})$ to model the autocovariance 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.

Table 1. Prediction Error for Cases 1 and 2

Table 2 lists the computational time of the four methods in Case 1. Between the two nonstationary methods, DMGP-SS requires much less computation time than DMGP-GP. This exactly verifies the analysis in Section 3.3 that our model can save more time than the traditional nonstationary MGP method, because of a block-sparse covariance matrix. In fact, our model is scalable for larger sizes of data, which is described in Online Appendix E.

Table 2. Computational Time for the Different Methods in Case 1

5. 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 is dynamically correlated across different gestures (Xu et al. 2022), it is reasonable to utilize a nonstationary method to capture the varying correlation and leverage information among the signals of joints. Regarding the control policy iteration, we study a classical reinforcement learning (RL) problem, the mountain-car problem. We evaluate the performance of DMGP-SS on leveraging knowledge between different systems when the environment is nonstationary.

5.1. Movement Signal Modeling

5.1.1. Data Description.

We use the MSRC-12 gesture data set (Fothergill et al. 2012) 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 $\pm 2\hspace{0.17em}\text{cm}$. The data set 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. Figure 5(a) 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 12 outputs. We choose the z-coordinate movement of the left elbow as a target output, while using the remaining 11 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 rescale all outputs to $[-2,2]$. Figure 5(b) displays the 12 outputs and the change of joints’ positions.

5.1.2. 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-point-long time-series of each gesture, we randomly remove a sequence of length 10 as missing data.

First of all, Figure 6(a) shows the estimated correlation between the sources and the target. MGP-L1 selects six sources (the third “R-W-z,” the fourth “R-E-x,” the sixth “R-E-z,” the seventh “L-W-x,” the ninth “L-W-z,” and the 11th “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\le t\le 120$ and $120\le t\le 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 “Goggles” 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 Figure 5(b). 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 of the joints’ relationship.

Figure 6(b) displays the recovered target signal in one experiment. Notably, DMGP-SS accurately recovers the missing data with high precision. Besides, it has the least uncertainty because it can select 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. Table 3 compares the prediction accuracy of four methods in terms of both the mean absolute error (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 (Hersbach 2000):

Table 3. Prediction Error for Real Case

where ${\widehat{y}}_i$ is the predicted output and $\mathrm{\Phi }$ is the cumulative density function (CDF) of the predicted distribution. A low CRPS value suggests that the predicted posteriors are concentrated around the true value, indicating a high probabilistic forecast performance. As expected, DMGP-SS outperforms the other methods because of 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 the fact that nonsparse correlations lead to potential negative transfer effects. Besides, because the sources’ smoothness is heterogeneous, it is inappropriate to use the same autocovariance function to model all outputs.

5.2. 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:

1. Model estimation: estimating the state transition function U(s, a) based on the observed transition samples $[(s,a),u]$.

2. 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 can be found in the works on GP-based reinforcement learning (Kuss and Rasmussen 2003, Verstraeten et al. 2020).

In this study, we employ the well-known reinforcement learning problem, the mountain-car problem, to demonstrate the application of our model in decision making. Figure 7 illustrates such a problem. A car begins in a valley and aims to reach a goal position of the rightmost hill. Because of 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 the system, the agent state is described by the position $s^{\text{pos}}\in [-1.2,1.0]$ and the velocity $s^{\text{vel}}\in [-0.07,0.07]$ of the car, that is, $\mathit{s}=(s^{\text{pos}},s^{\text{vel}})$. The agent action is a horizontal force $a\in [-1,1]$. The car starts at the state ${\mathit{s}}_{\text{init}}=(s_{\text{init}}^{\text{vel}},s_{\text{init}}^{\text{vel}})$ with $s_{\text{init}}^{\text{pos}}\in [-0.6,-0.5]$ and $s_{\text{init}}^{\text{vel}}=0$, aiming to reach the goal state ${\mathit{s}}_{\text{goal}}=(0.45,0)$. The reward function is the probability density function of $N({\mathit{s}}_{\text{goal}},\text{diag}\{{0.05}^2,{0.0035}^2\})$. The dynamic equation of this system is approximated by

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

We consider the control problem involving a car in a nonstationary environment with limited transition samples. During $t\in [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\in [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 transition model, we transfer information from two historical source data sets, 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 2 summarizes the workflow of reinforcement learning, where we employ the DMGP-SS as a transition model.

Algorithm 2

(Control Policy Optimization for Mountain-Car Case)

Input: Source transition samples ${\{({\mathit{s}}_{i,n},a_{i,n}),{\mathit{u}}_{i,n}\}}_{i,n=1}^{2,200}$, target transition samples ${\{({\mathit{s}}_{3,t},a_{3,t}),{\mathit{u}}_{3,t}\}}_{t=1}^{40}$, reward function $\{r(\mathit{s})\}$, 256 support points $S_{\text{supp}}$.

• 1: Initialize: policy $W(\mathit{s})\leftarrow$ random policy, value $V_{\text{supp}}=r(S_{\text{supp}})$.

• 2: Fit two DMGP-SS models for the state transition using both the source and the target transition samples, one for the position state and the other for the velocity state.

• 3: Fit a GP for the state value model $V(\mathit{s})$ using $\{S_{\text{supp}},V_{\text{supp}}\}$.

• 4: Improve policy $W(\mathit{s})$ and value model $V(\mathit{s})$ iteratively based on the state-transition model until $V(S_{\text{supp}})$ converges.

• 5: Execute the optimized policy starting at the state ${\mathit{s}}_{\text{init}}$:

For simplicity, we model the source transition data with stationary GPs in DMGP-SS. The sources and the target are expressed as

Here, u is the position or velocity state to which the agent transitioned from a state s after taking an action a. After fitting this model to both the target and source samples, we train a GP for the value model $V(\mathit{s})$ with even-spaced support points $S_{\text{supp}}$ and rewards $r(S_{\text{supp}})$. Based on the transition model $U(\mathit{s},\mathit{a})$, we iteratively improve the policy $W(\mathit{s})$ and value model $V(\mathit{s})$ until the prediction of $V(\mathit{s})$ converges. Details on the policy improvement procedure can be found in Kuss and Rasmussen (2003) and Verstraeten et al. (2020). Although 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 (Kuss and Rasmussen 2003, Verstraeten et al. 2020) with the stationary GP and MGP as the state transition model, respectively. The maximum execution steps are 600 for each method. In Table 4, we report 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, Figure 8 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, because their stationary transition models cannot account for the environment change and make the policy iteration hard to converge.

Table 4. Predictive Accuracy on State Transition and the Mean of Absolute Distances to the Goal State in the Mountain-Car Case

In Table 4, 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 other methods, because it can capture the change of environment and transfer information from the related sources. Figure 9 illustrates the estimated ${\mathit{\alpha }}_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.

6. Conclusion

This paper proposes a flexible nonstationary multioutput 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 nonzero correlation parameters. Through the experiments on simulation and a human gesture data set, we demonstrate that our method is well suited for both capturing nonstationary 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 for 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 RL, a popular framework for solving sequential decision-making problems. However, negative transfer is still a notable challenge for multitask reinforcement learning and it is risky to naively share information across all tasks (Yang et al. 2020, Zhu et al. 2023). Therefore, we embedded our model with the offline RL transfer framework to automatically select correlated sources to share knowledge in a nonstationary environment (Padakandla et al. 2020). 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 correlations, scaling it up to large data sets 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.