MAGMA: inference and prediction using multi-task Gaussian processes with common mean

A novel multi-task Gaussian process (GP) framework is proposed, by using a common mean process for sharing information across tasks. In particular, we investigate the problem of time series forecasting, with the objective to improve multiple-step-ahead predictions. The common mean process is defined as a GP for which the hyper-posterior distribution is tractable. Therefore an EM algorithm is derived for handling both hyper-parameters optimisation and hyper-posterior computation. Unlike previous approaches in the literature, the model fully accounts for uncertainty and can handle irregular grids of observations while maintaining explicit formulations, by modelling the mean process in a unified GP framework. Predictive analytical equations are provided, integrating information shared across tasks through a relevant prior mean. This approach greatly improves the predictive performances, even far from observations, and may reduce significantly the computational complexity compared to traditional multi-task GP models. Our overall algorithm is called Magma (standing for Multi tAsk GPs with common MeAn). The quality of the mean process estimation, predictive performances, and comparisons to alternatives are assessed in various simulated scenarios and on real datasets.

A novel multi-task Gaussian process (GP) framework is proposed, by using a common mean process for sharing information across tasks.In particular, we investigate the problem of time series forecasting, with the objective to improve multiple-step-ahead predictions.The common mean process is defined

Introduction
Gaussian processes (GPs) are a powerful tool, widely used in machine learning (Bishop, 2006;Rasmussen and Williams, 2006).The classic context of regression aims at inferring the underlying mapping function associating input to output data.In a probabilistic framework, a typical strategy is to assume that this function is drawn from a prior GP.Doing so, we may enforce some properties for the function solely by characterising the mean and covariance functions of the process, the latter often being associated with a specific kernel.This covariance function plays a central role and GPs are an example of kernel methods.We refer to Álvarez et al. (2012) for a comprehensive review.On the other hand, the mean function is generally set to 0 for all entries assuming that the covariance structure already integrates the desired relationship between observed data and prediction targets.In this paper, we consider a novel multi-task learning framework where a series of GPs share a common mean, expressed as a GP as well.We demonstrate that modelling the mean function as such can be key to obtain more relevant predictions.

Related work
The multi-task framework consists in using data from several tasks (or individuals) to improve learning or predictive capacities compared to an isolated model.It has been introduced by Caruana (1997) and then adapted in many fields of machine learning.GP versions of such models were introduced in Schwaighofer et al. (2004), which proposed an Expectation-Maximisation (EM) algorithm for learning.Similar techniques can be found in Shi et al. (2005).Meanwhile, Yu et al. (2005) offered an extensive study of the relationships between the linear model and GPs to develop a multi-task GP formulation.However, since the introduction in Bonilla et al. (2008) of the idea of two matrices, modelling covariance between inputs and tasks respectively, the term multi-task Gaussian process has mostly referred to the choice made regarding the covariance structure.Some further developments were discussed by Hayashi et al. (2012), Rakitsch et al. (2013) and Zhu and Sun (2014).In particular, an interesting approach in Nguyen and Bonilla (2014) proposed a sparse approximation for multi-task GP inference.More generally, these approaches are known as examples of linear models of coregionalization (LMC) in the geostatistics literature, and Álvarez and Lawrence (2011) provides a unified view on the topic as well as an efficient strategy for constructing computationally efficient approximations.Let us emphasise that the present paper is not based on the same assumptions and principles, and aims at defining a different multi-task paradigm for GPs, focusing on sharing information through the mean function rather than the covariance structure.Besides, the work of Swersky et al. (2013) on Bayesian hyper-parameter optimisation in such LMC models is also worth a mention.Real applications were tackled by similar models in Williams et al. (2009) and Alaa and van der Schaar (2017), while Clingerman and Eaton (2017) and Moreno-Muñoz et al. (2019) developed continual learning methods for multi-task GP.
As we focus on multi-task time series forecasting, a connection can be drawn to the study of multiple curves, or functional data analysis (FDA).As initially proposed in Rice and Silverman (1991), it is possible to model and learn mean and covariance structures simultaneously in this context.We refer to the monographs Ramsay and Silverman (2005) and Ferraty and Vieu (2006) for a comprehensive introduction to FDA.In particular, these books introduced several usual ways for modelling a set of functional objects in frequentist frameworks, for example by using a decomposition in a basis of functions (such as Bsplines, wavelets, Fourier).This kind of B-splines decomposition was used in Shi et al. (2007) for modelling the mean function in a generative model that somehow resembles ours.Subsequently, some Bayesian alternatives were developed in Thompson andRosen (2008), andCrainiceanu andGoldsmith (2010).

Our contributions
A multi-task GP framework with a common mean process is introduced, allowing reliable probabilistic forecasts even in multiple-step-ahead problems, or for sparsely observed individuals.For this purpose, (i) we introduce a GP model where the specific covariance structure of each task is defined through a separate kernel and its associated set of hyperparameters, whereas the common mean function µ 0 allows sharing information across tasks and overcomes the weaknesses of classic GPs in making predictions far from observed data.To account for uncertainty, we propose a hierarchical formulation to define the common mean process µ 0 as a GP as well.(ii) We derive an algorithm called Magma (available as an R package at https://github.com/ArthurLeroy/MagmaClustR) to compute µ 0 's hyper-posterior distribution together with the estimation of hyper-parameters in an EM fashion, and discuss its computational complexity.(iii) We enrich Magma with explicit formulas to make predictions for any new, partially observed, task.The hyper-posterior distribution of µ 0 provides a prior belief on what we would expect to observe before seeing any new data, acting as an already-informed mean process, integrating both trend and uncertainty coming from other tasks.(iv) We illustrate the performance of our method on synthetic and two real-life datasets and obtain state-of-the-art results compared to alternative approaches.

Outline
The paper is organised as follows.We introduce our multi-task Gaussian process model in Sec. 2, along with notation.Sec. 3 is devoted to the inference procedure, with an Expectation-Maximisation (EM) algorithm to estimate the Gaussian process hyperparameters and µ 0 's hyper-posterior.We leverage this strategy in Sec. 4 and derive a prediction algorithm.In Sec. 5, we analyse and discuss the computational complexity of both the inference and prediction procedures.Our methodology is illustrated in Sec.6, with a series of experiments on both synthetic and real-life datasets, and a comparison to competing state-of-the-art algorithms.On those tasks, we provide empirical evidence that our algorithm outperforms other approaches.Sec.7 draws perspectives for future work, and we defer some proofs to original results claimed in the paper to Sec. 8.

Notation
While GPs can handle many types of data, their continuous nature makes them particularly well suited to study temporal phenomena.Throughout, the term individual is used as a synonym of task or batch, and we adopt notation and vocabulary of time series to remain consistent with the application on real dataset provided in Sec.6.5, which addresses young swimmers performances' forecast.
We are provided with functional data coming from M ∈ I different individuals, where I ⊂ N.For each individual i, we observe a set of inputs {t 1 i , . . ., t N i i } and associated outputs {y i (t 1 i ), . . ., y i (t N i i )}, where N i is the number of data points for the i-th individual.Since many objects are defined for all individuals, we shorten our notation as follows: for any object x existing for all i, we denote {x i } i = {x 1 , . . ., x M }.Moreover, as we work in a temporal context, the inputs are referred to as timestamps.In the specific case where all individuals are observed at the same timestamps, we call the grid of observations common.On the contrary, a grid of observations is uncommon if the timestamps are different in number and/or location among the individuals.Some convenient notation follows: , the set of timestamps for the i-th individual, • y i = y i (t i ), the vector of outputs for the i-th individual, • t = M i=1 t i , the pooled set of timestamps among individuals, • N = card(t), the total number of observed timestamps.

Model and hypotheses
Suppose that functional data are coming from the sum of a mean process, common to all individuals, and an individual-specific centred process.To clarify relationships in the generative model, we illustrate our graphical model in Fig. 1.Let T be the input space, our model is are respectively the common mean and individual specific processes.Moreover, the error term is supposed to be i (•) ∼ N (0, σ 2 i I).The following notation is used for parameters: • m 0 (•), an arbitrary prior mean function, • k θ 0 (•, •), a covariance kernel of hyper-parameters θ 0 , • ∀i ∈ I, c θ i (•, •), a covariance kernel with hyper-parameters θ i , • σ 2 i ∈ R + , the noise variance associated with the i-th individual, • ∀i ∈ I, we define the shorthand Figure 1: Graphical model of dependencies between variables in the Multi-task Gaussian Process model.
• Θ = {θ 0 , {θ i } i , σ 2 i i }, the set of all hyper-parameters to learn in the model.
We also assume that: • ∀i ∈ I, µ 0 , f i and i are independent.
It follows that {y i | µ 0 } i=1,...,M are independent from one another, and for all i ∈ I: Let us emphasise that this property only holds conditionally to µ 0 .Otherwise, once µ 0 is integrated out, the y i are no longer independent.Here, we do not assume any specific covariance structure between individuals contrarily to standard LMC approaches.As we shall see in the next sections, the process µ 0 will be key to handle the dependencies and share information across the individuals.
Although this model is based on infinite-dimensional GPs, the inference will be conducted on a finite grid of observations.According to the aforementioned notation, we observe {(t i , y i )} i , and the corresponding likelihoods are Gaussian: might be different among individuals, we also need to evaluate µ 0 on the pooled grid of timestamps t: An alternative hypothesis consists in considering hyper-parameters {θ i } i and σ 2 i i equal for all individuals.We call this hypothesis Common HP (where HP stands for hyper-parameters) in the Sec. 6.This particular case represents a context where individuals correspond to different trajectories of the same process, whereas different hyper-parameters indicate different covariance structures and thus a more flexible model.For the sake of generality, the remainder of the paper is written with θ i and σ 2 i notation, when there are no differences in the procedure.Moreover, the model above and the subsequent algorithm may use any form of covariance function, often parametrised by a finite set (usually small) of hyper-parameters.For example, a common kernel in the GP literature is known as the Exponentiated Quadratic kernel (also called sometimes Squared Exponential or Radial Basis Function kernel).It solely depends on two hyper-parameters θ = {v, } and is defined as: The Exponentiated Quadratic kernel is simple and enjoys useful smoothness properties.This is the kernel used in the current version of our implementation (see Sec. 6 for details).Note that there is a rich literature on kernel choice, their construction and properties, which is beyond the scope of the present work: we refer to Rasmussen and Williams (2006) or Duvenaud (2014) for comprehensive studies.

Learning
Several approaches to learn hyper-parameters for Gaussian processes have been proposed in the literature, we refer to Rasmussen and Williams (2006) for a comprehensive study.One classical approach, called empirical Bayes (Casella, 1985), is based on the maximisation of an explicit likelihood to estimate hyper-parameters.This procedure avoids sampling from intractable distributions, usually resulting in additional computational cost and complicating practical use in moderate to large sample sizes.As previously stated, once µ 0 is marginalised out, the log-likelihood cannot be written as a sum of Gaussian log-likelihoods any more.Therefore, we propose an EM algorithm (see the pseudocode in Algorithm 1) to learn the hyper-parameters Θ in this context.The procedure alternatively computes the hyper-posterior distribution p(µ 0 | (y i ) i , Θ) with current hyper-parameters, and then optimises Θ according to this hyper-posterior distribution.This EM algorithm converges to local maxima (Dempster et al., 1977), typically in a handful of iterations.

E step
For the sake of simplicity, we assume in that section that ∀i, j ∈ I, t i = t j = t, i.e. the individuals are observed on a common grid of timestamps.We provide a generalisation of the following proposition in Sec. 4 (Proposition 4.1), where the result holds for uncommon grids.The E step then consists in computing the hyper-posterior distribution of µ 0 (t).
Proposition 3.1.Assume the hyper-parameters Θ known from initialisation or estimated from a previous M step.The hyper-posterior distribution of µ 0 remains Gaussian: with Proof.
We omit specifying timestamps in what follows since each process is evaluated on t.Therefore, we can write: The term L 1 = −(1/2) log p(µ 0 | {y i } i , Θ) may then be written as where the constant terms are gathered into C 1 , C 2 ∈ R. Identifying terms in the quadratic form with the Gaussian likelihood, we get the desired result.
The maximisation step depends on the assumptions on the generative model, resulting in two versions for the EM algorithm (the E step is common to both, the branching point is here).

M step: different hyper-parameters
Assuming each individual has its own set of hyper-parameters {θ i , σ 2 i }, the M step is given by the following procedure.
inducing M + 1 independent maximisation problems: Proof.One simply has to distribute the conditional expectation in order to get the right likelihood to maximise, and then notice that the function can be written as a sum of M +1 independent (with respect to the hyper-parameters) terms.Moreover, by rearranging, one can observe that each independent term is the sum of a Gaussian likelihood and a correction trace term.See Sec.8.2 for details.

M step: common hyper-parameters
Alternatively, assuming all individuals share the same set of hyper-parameters {θ, σ 2 }, the M step is given by the following procedure.
For a set of hyper-parameters Θ = {θ 0 , θ, σ 2 }, optimal values are given by inducing two independent maximisation problems: where Proof.We use the same strategy as for Proposition 3.2, see Sec. 8.2 for details.
In both cases, explicit gradients associated with the likelihoods to maximise are available, facilitating the optimisation with gradient-based methods.

Initialisation
To implement the EM algorithm described above, several constants must be (appropriately) initialised: • m 0 (•), the mean parameter from the hyper-prior distribution of the process µ 0 (•).A somewhat classical choice in GP is to set its value to a constant function, typically 0 in the absence of external knowledge.Notice that, in our multi-task framework, the influence of m 0 (•) in hyper-posterior computation decreases as M grows anyway (see Proposition 3.1).
• Initial values for kernel parameters θ 0 and {θ i } i .Those strongly depend on the chosen kernel and its properties.We advise initiating θ 0 and {θ i } i with close values, as a too large difference might induce nearly singular covariance matrices and result in numerical instability (typical in GPs applications).In such pathological regime, the influence of a specific individual tends to overtake others in the calculus of µ 0 's hyper-posterior distribution.
• Initial values for the variance of the error terms σ 2 i i .This choice mostly depends on the context and properties of the dataset.We suggest avoiding initial values with more than an order of magnitude different from the variability of data.In particular, a too high value might result in a model mostly capturing noise.
As a final note, let us stress that the EM algorithm depends on the initialisation and is only guaranteed to converge to local maxima of the likelihood function (McLachlan and Krishnan, 2007).Several strategies have been considered in the literature to tackle this issue such as simulated annealing (Ueda and Nakano, 1998) or repeated short runs (Biernacki et al., 2003).In this work, we chose the latter option.

Pseudocode
We wrap up this section with the pseudocode of the EM component of our complete algorithm, which we call Magma (standing for Multi tAsk Gaussian processes with common MeAn).The corresponding code is available at https://github.com/ArthurLeroy/MAGMA.

Discussion of EM algorithms and alternatives
Let us stress that even though we focus on prediction purpose in this paper, the output of the EM algorithm already provides results on related FDA problems.The generative model in Yang et al. (2016) describes a Bayesian framework that resembles ours to smooth multiple curves simultaneously.However, modelling variance structure with an Inverse-Wishart process forces the use of an MCMC algorithm for inference or the introduction of a more tractable approximation in Yang et al. (2017).One can think of the learning through Magma and applying a single task GP regression on each individual as an empirical Bayes counterpart to their approach.Meanwhile, µ 0 's hyper-posterior distribution also provides the probabilistic estimation of a mean curve from a set of functional data.The closest method to our approach can be found in Shi et al. (2007) and the following book Shi and Choi (2011).The authors also work in the context of a multi-task GP model, and one can retrieve the idea of defining a mean function µ 0 to overcome the weaknesses of classic GPs in making predictions far from observed data.However, since their model uses B-splines to estimate this mean function, the method only works if all individuals share the same grid of observations, and does not account for uncertainty over µ 0 .

Prediction
Once the hyper-parameters of the model have been learned, we can focus on our main goal: prediction for new individuals at unobserved timestamps.Since Θ is known and for the sake of concision, we omit conditioning on Θ in the sequel.Note there are two cases for prediction (referred to as Type I and Type II in Shi and Cheng, 2014, Section 3.2.1),depending on whether we observe some data or not for any new individual we wish to predict on.We denote by the index * a new individual for whom we want to make a prediction, say at timestamps t p .If there are no available data for this individual, we have no * -specific information, and the prediction is merely given by p(µ 0 (t p ) | {y i } i ).This quantity may be considered as the 'generic' (or Type II ) prediction according to the trained model, and only informs us through the mean process.Computing p(µ 0 (t p ) | {y i } i ) is also one of the steps leading to the prediction for a partially observed new individual (Type I ).The latter being the most compelling case, we consider Type II prediction as a particular case of the full Type I procedure, described below.
If we observe {t * , y * (t * )} for the new individual, the multi-task GP prediction is obtained in our model by computing the posterior distribution p(y * (t p ) | y * (t * ), {y i } i ).Note that the conditioning is taken over y * (t * ), as for any GP regression, but also on {y i } i , which is specific to our multi-task setting.The procedure for computing this distribution requires to successively complete the following steps: 1. choose a grid of prediction t p and define the pooled vector of timestamps t

Posterior inference on the mean process
As mentioned above, we observed a new individual at timestamps t * .The GP regression consists in arbitrarily choosing a vector t p of timestamps for which we aim at making predictions.Then, we define new notation for the pooled vector of timestamps t p * = t p t * , which will serve as a working grid to define the prior and posterior distributions involved in the prediction process.One can note that, although not mandatory in theory, it is often a good idea to include the observed timestamps of training individuals, t, within t p * since they match locations that contain information for the mean process to 'help' the prediction.In particular, if t p * = t, the computation of µ 0 's hyper-posterior distribution is not necessary since p(µ 0 (t) | {y i } i ) has previously been obtained from the EM algorithm.However, in general, it is necessary to compute the hyper-posterior p(µ 0 (t p * ) | {y i } i ) at the new timestamps.The idea remains similar to the E step aforementioned, and we obtain the following result.
Proposition 4.1.Let t p * be a vector of timestamps of size Ñ .The hyper-posterior distribution of µ 0 remains Gaussian: with: where we used the shortening notation: Proof.The sketch of the proof is similar to Proposition 3.1 in the E step.The only technicality consists in dealing carefully with the dimensions of vectors and matrices involved, and whenever relevant, to define augmented versions of y i and Ψ θ i , σ 2 i with 0 elements at unobserved timestamps' position for the i-th individual.Note that if we pick a vector t p * including only some of the timestamps from t i , information coming from y i at the remaining timestamps is ignored.We defer details to Sec. 8.1.

Computing the multi-task prior distribution
According to our generative model, given the mean process, any new individual * is modelled as: . Therefore, for any finite-dimensional vector of timestamps, and in particular for t p * , p(y * (t p . (3) Proof.To compute this prior, we need to integrate out the mean process µ 0 in p(y * | µ 0 , {y i } i ), whereas the multi-task aspect remains through the conditioning over {y i } i .We omit the writing of timestamps, by using the simplified notation µ 0 and y * instead of µ 0 (t p * ) and y * (t p * ), respectively.We first use the assumption that {y i | µ 0 } i∈{1,...,M } ⊥ ⊥ y * | µ 0 , i.e., the individuals are independent conditionally to µ 0 .Then, one can notice that the two distributions involved within the integral are Gaussian, which leads to the explicit Gaussian target distribution after integration.
This convolution of two Gaussians remains Gaussian (Bishop, 2006, Chapter 2.3.3).The mean parameter is then given by Following the same idea, the second-order moment is given by Note that the process y * (•) | {y i } i is not strictly a GP, although its finite-dimensional evaluation (3) remains Gaussian.The covariance structure cannot be expressed as a kernel that could be directly evaluated at any timestamps: the process is known as a degenerated GP.In practice however, this does not bear much consequence as any arbitrary vector of timestamps τ can be chosen at first, and computing hyper-posterior p(µ 0 (τ ) | {y i } i ) still yields to the Gaussian distribution p(y * (τ ) | {y i } i ) as above.For the sake of simplicity, we now rename the covariance matrix of the multi-task prior distribution: where the indices in the blocks of the matrix correspond to the associated timestamps t p and t * .

Learning the new hyper-parameters
When we collect data points for a new individual, as in the single-task GPs setting, we would need to learn the hyper-parameters of its covariance kernel before making predictions.A salient fact in our multi-task approach is that we consider this step being part of the prediction process, for two main reasons.First, the model is already trained for individuals i = 1, . . ., M , and this training is independent of the future individual * or the choice of prediction timestamps.Since learning these new hyper-parameters requires knowledge of µ(t p * ) and thus of the prediction timestamps, we cannot compute them beforehand.Second, learning these hyper-parameters with the empirical Bayes approach only requires maximisation of a Gaussian likelihood which is negligible in computing time compared to the previous EM algorithm.As for single-task GP, we have the following estimates for hyper-parameters: Note that this step is optional depending on the modelling assumption: in the common hyper-parameters model (i.e.(θ, σ 2 ) = (θ i , σ 2 i ), ∀i ∈ I), any new individual will also share the same hyper-parameters and we already have Θ * = ( θ * , σ 2 * ) = ( θ, σ 2 ) from the EM algorithm.

Prediction
We can rewrite the multi-task prior distribution, by separating observed and prediction timestamps, as: As usual, the conditional distribution remains Gaussian, and the multi-task posterior distribution is given by: where: Although this predictive distribution presents a formulation nicely analogous to standard GPs, let us emphasise on the terms m 0 (t p * ) and Γ p * , which embed crucial information from training individuals for the mean prediction to be more relevant even in far from the observed points y * (t * ).

Complexity analysis for training and prediction
Computational complexity is of paramount importance in GPs as it quickly scales with large datasets.The classical cost to train a GP is O(N 3 ), and O(N 2 ) for prediction (Rasmussen and Williams, 2006) where N is the number of data points (although there exist various sparse approximations, see Sec. 7 for references).Moreover, multi-task GP models lying on LMC approaches typically present a complexity of O(M 3 N 3 ) in training, which can be diminished when using sparse approximations ( Álvarez and Lawrence, 2011).As detailed below, our model reaches a reduction to O((M +1)N 3 ) for the training complexity in a similar context (common grid of timestamps for all individuals), without using any sparse approximation.
More specifically, since Magma uses information from M individuals, each of them providing N i observations, these quantities determine the overall complexity of the algorithm.
If we recall that N is the number of distinct timestamps (i.e.N ≤ e. the complexity of each EM iteration).As usual with GPs, the cubic costs come from the inversion of the corresponding matrices, and here, the constant is proportional to the number of iterations of the EM algorithm.The dominating term in this expression depends on the values of M , relatively to N .For a large number of individuals with many common timestamps (M N i N ), the first term dominates.For diverse timestamps among individuals (M N i N ), the second term becomes the primary burden, as in any GP problem.During the prediction step, the re-computation of µ 0 's hyper-posterior implies the inversion of a Ñ × Ñ (dimension of t p * ) which has a O( Ñ 3 ) complexity while the new hyper-parameters estimation's cost is O(N 3 * ).In practice, the most computationally-expensive steps can be performed in advance to allow for quick onthe-fly prediction when collecting new data.If we observe the training dataset once and pre-compute the hyper-posterior of µ 0 on a fine grid on which to predict later, the immediate computational cost for each new individual is identical to the one of the single-task GP regression.

Experimental results
We evaluate our Magma algorithm on synthetic data and two real datasets.The classical GP regression on single tasks separately is used as the baseline alternative for predictions.While it is not expected to perform well on the dataset used, the comparison highlights the interest of multi-task approaches.To our knowledge, the only alternative to Magma is the GPFDA algorithm from Shi et al. (2007); Shi and Choi (2011), described in Sec.3.4, and the associated R package GPFDA, which is applied during the experiments.Throughout the section, the standard Exponentiated Quadratic kernel (see Eq. ( 1)) is used both for simulating the data and for modelling the covariance structures in the three algorithms.Hence, each kernel is associated with θ = {v, }, v, ∈ R + , a set of variance and lengthscale hyper-parameters, respectively.Each simulated dataset has been drawn from the sampling scheme below:  6. ∀i ∈ I, draw a subset t i ⊂ t of N i = 30 timestamps uniformly, and draw This procedure provides a synthetic dataset {t i , y i } i , and its associated mean process µ 0 (t).Those quantities are used to train the model, make predictions with each algorithm, and then compute errors in µ 0 estimation and forecasts.We recall that the Magma algorithm enables two different settings depending on the model's assumption over hyper-parameters (HP), and we refer to them as Common HP and Different HP in the following.In order to test these two contexts, differentiated datasets have been generated, by drawing Common HP data or Different HP data for each individual at step 5.We previously presented the idea of the model used in GPFDA, and, although the algorithm has many features (in particular about the type and number of input variables), it is not yet usable when timestamps are different among individuals.Therefore, two frameworks are considered, Common grid and Uncommon grid, to take this specification into account.Thus, the comparison between the different methods can only be performed on data generated under the settings Common HP and Common grid, and the effect of those different settings on Magma is analysed separately.Moreover, the initialisation for the prior mean function, m 0 (•), is set to be constant, equal to 0 for each algorithm.Except in some experiments, where the influence of the number of individuals is analysed, the generic value is M = 20.In the case of prediction on unobserved timestamps for a new individual, the first 20 data points are used as observations, and the remaining 10 are taken as test values.Optimisation of the hyper-parameters is performed by likelihood maximisation, using the L-BFGS-B algorithm (Nocedal, 1980;Morales and Nocedal, 2011) in all methods.The convergence criterion for all algorithms is reached if the difference of log-likelihood between two iterations is lower than 10 −2 .In general, the EM algorithm in Magma converges in a few iterations, typically fewer than 5 with the Common HP setting, and rarely more than 15 even with the Different HP setting.

Illustration on a simple example
To illustrate the multi-task approach of Magma, Fig. 2 displays a comparison between standard GP regression and Magma on a simple example, from a dataset simulated according to the scheme above and using the Uncommon grid /Common HP setting.Given the observed data (in black), values on a thin grid of unobserved timestamps are predicted and compared, in particular, with the true test values (in red).As expected, the GP regression provides a good fit close to the data points and then dives rapidly to the prior 0 with increasing uncertainty.Conversely, although the initialisation for the prior mean is 0 in Magma as well, the hyper-posterior distribution of µ 0 (dashed line) is estimated thanks to all individuals in the training dataset.This process acts as an informed prior helping GP prediction for the new individual, even far from its own observations.More precisely, 3 phases can be distinguished according to the level of information coming from the data: in the first one, close to the observed data (t ∈ [ 1, 7 ]), the two processes behave similarly, except for a slight increase in the variance for Magma, which is logical since the prediction also takes uncertainty over µ 0 into account (see Eq. ( 3)); in the second one, on intervals of unobserved timestamps containing data points from the training dataset (t ∈ [ 0, 1 ]∪[ 7, 10 ]), the prediction is guided by the information coming from other individuals through µ 0 .In this context, the mean trajectory remains coherent and the uncertainty increases only slightly.In the third phase, where no observations are available, neither from the new individual nor from the training dataset (t ∈ [ 10, 12 ]), the prediction behaves as expected, with a slow drifting to the prior mean 0, with highly increasing variance.Overall, the multi-task framework provides reliable probabilistic predictions on a wider range of timestamps, potentially outside of the usual scope of GPs.

Performance comparison on simulated datasets
We confront the performance of Magma to alternatives in several situations and for different datasets.In the first place, the classical GP regression (GP), GPFDA and Magma are compared through their performance in prediction and estimation of the true mean process µ 0 .In the prediction context, the performances are evaluated according to the following indicators: • the mean squared error (MSE) which compares the predicted values to the true test values of the 10 last timestamps: 1 10 Figure 3: MSE with respect to the number M of training individuals (boxplots are displayed from 100 runs in each case).Left: prediction error on 10 testing points.
Right: estimation error of the true mean process µ 0 .The CIC 95 provides insights on the reliability of the predictive variance and should be as close to the value 95% as possible.Other values would indicate a tendency to underestimate or overestimate the uncertainty.Let us recall that GPFDA uses B-splines to estimate the mean process and does not account for uncertainty, contrarily to a probabilistic framework as Magma.However, a measure of uncertainty based on an empirical variance estimated from training curves is proposed (see Shi and Cheng, 2014, Section 3.2.1).In practice, this measure constantly overestimates the true variance, and their 95% empirical interval coverage is generally equal or close to 100%.
In the estimation context, the performances are evaluated thanks to another MSE, which compares the estimations to the true values of µ 0 at all timestamps: Table 1 presents the results obtained over 100 datasets, where the models are trained on M = 20 individuals, each of them observed on N = 30 common timestamps.As expected, both multi-task methods lead to better results than GP.However, Magma outperforms GPFDA, both in the estimation of µ 0 and in predictive performance.In terms of error as well as in uncertainty quantification, Magma provides more accurate results, in particular with a CI 95 coverage close to the 95% expected value.Each method presents a quite high standard deviation for MSE in prediction, which is due to some datasets with particularly difficult values to predict, although most of the cases lead to small errors.This behaviour is reasonably expected since such 10-timestamps-ahead forecasts might sometimes be tricky.It can also be noticed on Fig. 3 that Magma consistently provides lower errors as well as less pathological behaviour, as it may sometimes occur with the B-splines modelling used in GPFDA.
To highlight the effect of the number of individuals M on the performance, Fig. 3 provides the same 100 runs trial as previously, for different values of M .The boxplots exhibit, for each method, the behaviour of the prediction and estimation MSE as information is added in the training dataset.Let us mention the absence of discernible changes as soon as M > 200.As expected, we notice on the right panel that adding information from new individuals improves the estimation of µ 0 , leading to shallow errors for high values of M , in particular for Magma.Meanwhile, the left panel exhibits reasonably unchanged prediction performance with respect to the values of M , excepted for some random fluctuations.This property is expected for GP regression since no external information is used from the training dataset in this context.For both multi-tasks algorithms though, the estimation of µ 0 improves the prediction by one order of magnitude below the typical errors, even with only a few training individuals.Furthermore, since a new individual behaves independently through f * , it is natural for a 10-points-ahead forecast to present intrinsic variations, despite an adequate estimation of the shared mean process.
To illustrate the advantage of multi-task methods, even for M = 20, we display on Fig. 4 the evolution of MSE according to the number of timestamps N that are assumed to be observed for the new individual on which we make predictions.These predictions remain computed on the last 10 timestamps, although in this experiment, we only observe the first 5, 10, 15, or 20 timestamps, in order to change the volume of information and the distance from training observations to targets.We observe on Fig. 4 that, as expected in a GP framework, the closer observations are to targets, the better the results.However, for multi-tasks approaches and in particular for Magma, the prediction remains consistently adequate even with few observations.Once more, sharing information across individuals significantly helps the prediction, even for small values of M or few observed data.

Magma's specific settings
As we previously discussed, different settings are available for Magma according to the nature of data and the model hypotheses.First, the Common grid setting corresponds to cases where all individuals share the same timestamps, whereas Uncommon grid is used otherwise.Moreover, Magma enables to consider identical hyper-parameters for all individuals or specific ones, as previously discussed in Sec.2.2.To evaluate the effect of the different settings, performances in prediction and µ 0 's estimation are evaluated in the following cases in Table 2: • Common HP, when data are simulated with a common set of hyper-parameters for all individuals, and Proposition 3.3 is used for inference in Magma, • Different HP, when data are simulated with its own set of hyper-parameters for each individual, and Proposition 3.2 is used for inference in Magma, • Common HP on different HP data, when data are simulated with its own set of hyper-parameters for each individual, and Proposition 3.3 is used for inference in Magma.
Note that the first line of the table (Common grid / Common HP ) of Table 2 is identical to the corresponding results in Table 1, providing reference values, significantly better than for other methods.The results obtained in Table 2 indicate that the Magma performance is not significantly altered by the settings used or the nature of the simulated data.To confirm the robustness of the method, the setting Common HP was applied to data generated by drawing different values of hyper-parameters for each individual (Different HP data).In this case, performances in prediction and estimation of µ 0 are slightly deteriorated, although Magma still provides quite reliable forecasts.This experience also highlights a particularity of the Different HP setting: looking at the estimation of µ 0 performance, we observe a significant decrease in the CI 95 coverage, due to numerical instability in some pathological cases.Numerical issues, in particular during matrix inversions, are classical problems in the GP literature and, because of the potentially large number of different hyper-parameters to train, the probability for at least one of them to lead to a nearly singular matrix increases.In this case, one individual might overwhelm others in the calculus of µ 0 's hyper-posterior (see Proposition 4.1), and thus lead to an underestimated posterior variance.This problem does not occur in the Common HP settings, since sharing the same hyper-parameters prevents the associated covariance matrices from running over each other.Thus, except if one specifically wants to smooth multiple curves presenting really different behaviours, keeping Common HP as a default setting appears as a reasonable choice.Let us notice that the estimation of µ 0 is slightly 92.5 (15.9) 3.2 (4.5) 93.4 (9.8) better for common than for uncommon grid since the estimation problem on the union of different timestamps is generally more difficult.However, this feature only depends on the nature of data.

Running times comparisons
The counterpart of the more accurate and general results provided by Magma is a natural increase in running time.Table 3 exhibits the raw and relative training times for GPFDA and Magma (prediction times are negligible and comparable in both cases), on data coming from the simulation scheme with varying values of M on a Common grid of N = 30 timestamps.The algorithms were run under the 3.6.1 R version, on a laptop with a dualcore processor cadenced at 2.90GHz and an 8GB RAM.The reported computing times are in seconds, and for small to moderate datasets (N 10 3 , M 10 4 ) the procedures ran in few minutes to few hours.The difference between the two algorithms is due to GPFDA modelling µ 0 as a deterministic function through B-splines smoothing, whereas Magma accounts for uncertainty.The ratio of computing times between the two methods tends to decrease as M increases, and stabilises around 2 for higher numbers of training individuals.This behaviour comes from the E step in Magma, which is incompressible and quite insensitive to the value of M .Roughly speaking, one needs to pay twice the computing price of GPFDA for Magma to provide (significantly) more accurate predictions and uncertainty over µ 0 .Table 4 provides running times of Magma according to its different settings, with M = 20.Because the complexity is linear in M in each case, the ratio in running times would remain roughly similar no matter the value of M .Prediction time appears negligible compared to training time, and generally takes less than one second to run.Besides, the Different HP setting increases the running time since in this context M maximisations (instead of one for Common HP ) are required at each EM iteration.In this case, the prediction also takes slightly longer because of the necessity to optimise hyper-parameters for the new individual.Although the nature of the grid of timestamps does not matter in itself, a key limitation lies in the dimension N of the pooled set of timestamps, which tends to get bigger when individuals have different timestamps from one another.

Data and problematic
We consider the problem of performance prediction in competition for french swimmers.The French Swimming Federation provided us with an anonymised dataset, compiling the age and results of its members between 2000 and 2016.For each competitor, the race times are registered for competitions of 100m freestyle (50m swimming-pool).The database contains results from 1731 women and 7876 men, each of them compiling an average of 22.2 data points (min = 15, max = 61) and 12 data points (min = 5, max = 57), respectively.In the following, age of the i-th swimmer is considered as the input variable (timestamp t) and the performance (in seconds) on a 100m freestyle as the output (y i (t)).For reasons of confidentiality and property, the raw dataset cannot be published.The analysis focuses on the youth period, from 10 to 20 years, where the progression is the most noticeable.In order to get relevant time series, we retained only individuals having a sufficient number of data points (N i ≥ 5) on the considered time period.For a young swimmer, observed during its first years of competition, we aim at modelling its progression curve and make predictions on its future performance in the subsequent years.Since we consider a decision-making problem involving irregular time series, the GP probabilistic framework is a natural choice to work on.Thereby, assuming that each swimmer in the database is a realisation y i defined as previously, we expect Magma to provide multi-task predictions for a new young swimmer, that will benefit from information of other swimmers already observed at older ages.To study such modelling, and validate its efficiency in practice, we split the individuals into training and testing datasets with respective sizes: • M F train = 1039, for the female training set, • M F test = 692, for the female testing set, • M M train = 4726, for the male training set, • M M test = 3150, for the male testing set.
Inference on the hyper-parameters is performed thanks to the training dataset in both cases.Considering the different timestamps and the relative monotony of the progression curves, the settings Uncommon grid /Common HP has been used for Magma.The overall training lasted around 2 hours with the same hardware configuration as for simulations.To compute MSE and the CI 95 coverage, the data points of each individual in the testing set has been split into observed and testing timestamps.Since each individual has a different number of data points, the first 80% of timestamps are taken as observed, while the remaining 20% are considered as testing timestamps.Magma's predictions are compared with the true values of y i at testing timestamps.As previously, both GP and Magma have been initialised with a constant 0 mean function.Initial values for hyper-parameters are also similar for all i, θ ini 0 = θ ini i = (exp(1), exp(1)) and σ ini i = 0.4.Those values are the default in Magma and remain adequate in the context of these datasets.

Results and interpretation
The overall performance and comparison are summarised in Table 5.We observe that Magma still provides excellent results in this context, and naturally outperforms predictions provided by a standard GP regression.As the progression curves present relatively monotonic variations and thus avoid pathological behaviours that could occur with synthetic data, the MSE in prediction remains very low.The CI 95 coverage sticks close to the 95% expected value for Magma, indicating an adequate quantification of uncertainty.To illustrate these results, an example is displayed on Fig. 5 for both men and women.For a randomly chosen testing individual, we plot its predicted progression curve (in blue), where its first 15 data points are used as observations (in black), while the remaining true data points (in red) are displayed for comparison purpose.As previously observed in the simulation study, the simple GP quickly drifts to the prior 0 mean, as soon as data lack.However, for both men and women, the Magma predictions remain close to the true data, which also lie within the 95% credible interval.Even for long term forecast, where the mean prediction curve tends to overlap the mean process (dashed line), the true data remain in our range of uncertainty, as the credible interval widens far from observations.For clarity, we displayed only a few individuals from the training dataset (colourful points) in the background.The mean process (dashed line) seems to represent the main trend of progression among swimmers correctly, even though we cannot numerically compare µ 0 to any real-life analogous quantity.From a more sportrelated perspective, we can note that both genders present similar patterns of progression.However, while performances are roughly similar in mean trend before the age of 14, they start to differentiate afterwards and then converge to average times with approximatively a 5 seconds gap.Interestingly, the difference between men and women in terms of world records in swimming competitions for the 100m freestyle is currently 4.8 seconds (46.91 versus 51.71).These results, obtained under reasonable hypotheses on several hundreds of swimmers, seem to indicate that Magma would give quite reliable predictions for a new young swimmer.Furthermore, the uncertainty provided through the predictive posterior distribution offers an adequate degree of caution in a decision-making process.

Discussion
We have introduced a unified multi-task framework integrating a mean Gaussian process prior in the context of GP regression.While we believe that this process is an interesting object in itself, it also allows individuals to borrow information from each other and provides more accurate predictions, even far from data points.Furthermore, our method accounts for uncertainty in the mean process and remains applicable no matter eventual irregular timestamps in the grid of observations.The proposed algorithm, Magma, also presents a reduced computational complexity compared to previous multi-task GPs frameworks.Both on simulated and real-life datasets, we exhibited the efficiency of such an approach and studied some of its properties and possible settings.Magma outperforms the alternatives in estimation of the mean process as well as in prediction, and leads to a reliable quantification of uncertainty.We also displayed evidence of its predictive efficiency for real-life problems and provided some insights on practical interpretations about the mean process.
Despite the extensive literature on these aspects of GPs, our model does not yet include sparse approximations.While these aspects remain beyond the scope of the present paper, we might aim at adapting existing approaches (Snelson and Ghahramani, 2006;Quiñonero-Candela et al., 2007;Titsias, 2009) in our model to widen its applicability.Another possible avenue is an adaptation to the classification context, which is presented in Rasmussen and Williams (2006, Chapter 3).Besides, this work leaves the door open to improvement as we tackled here the problem of unidimensional regression: enabling either multidimensional or mixed type of inputs as in Shi and Choi (2011) would be of interest.To conclude, the hypothesis of a unique underlying mean process might be considered too restrictive for some datasets, and enabling cluster-specific mean processes would be a relevant extension.

Proofs
Note that the proof of Proposition 3.1 is a particular case of the proof below, where τ = t exactly (where τ is the set of timestamps the hyper-posterior is to be computed on).Moreover, in order to keep an analytical expression for µ 0 's hyper-posterior distribution, we discard the superfluous information contained in {y i } i at timestamps on which the hyper-posterior is not to be computed.Hence, the proof below states that the remaining data points are observed on subsets {τ i } i of τ .

Proof of Proposition 4.1
Let τ be a finite vector of timestamps, and {τ i } i such as ∀i ∈ I, τ i ⊂ τ .We define convenient notation: • m τ 0 = m 0 (τ ), • µ τ i 0 = µ 0 (τ i ), ∀i ∈ I, • Moreover, for a covariance matrix C, and u, v ∈ τ , we note [ C ] −1 uv the element of the inverse matrix at row associated with timestamp u, and column associated with timestamp v.We also ignore the conditionings over Θ, τ i and τ to maintain simple expressions.By construction of the models, we have: The term L 1 = −(1/2) log p(µ τ 0 | {y τ i i } i ) associated with the hyper-posterior remains quadratic and we may find the corresponding Gaussian parameters by identification: where we entirely decomposed the vector-matrix products.We factorise the expression according to the common timestamps between τ i and τ .Since for all i, τ i ⊂ τ , let us introduce a dummy indicator function 1 τ i = 1 {u,v∈τ i } to write: where the y i and Ψ i are completed by zeros: By identification of the quadratic form, we reach: p(µ τ 0 | {y τ i i } i ) = N µ τ 0 ; m 0 (τ ), K , with, 8.2 Proof of Proposition 3.2 and Proposition 3.3 Since the central part of the proofs is similar for both propositions, we detail the calculus by denoting Θ = {θ 0 , {θ i } i , σ 2 i i } for generality, and dissociating the two cases only when necessary.Before considering the maximisation, we notice that the joint density can be developed as: As we note that X and b play symmetrical roles in the calculus of the conditional expectation, we can apply the lemma regardless of the position of µ 0 in the M +1 equalities involved.Applying Lemma 8.1 to our previous expression of f (Θ), we obtain: We recall that, at the M step, m 0 (t) is a known constant, computed at the previous E step.Thus, we identify here the characteristic expression of several Gaussian loglikelihoods and associated correction trace terms.Moreover, each set of hyper-parameters is merely involved in independent terms of the whole function to maximise.Hence, the global maximisation problem can be separated into several maximisations of sub-functions according to the hyper-parameters getting optimised.Regardless to additional assumptions, the hyper-parameters θ 0 , controlling the covariance matrix of the mean process, appears in a function which is exactly a Gaussian log-likelihood, log N m 0 (t), m 0 (t), K t θ 0 , added to a corresponding trace term, − 1 2 Tr K t K t θ 0 −1 .This function can be maximised independently from the other parameters, giving the first part of the results in Proposition 3.2 and Proposition 3.3.
Although the idea is analogous for the remaining hyper-parameters, we have to discriminate here regarding the assumption on the model.If each individual is supposed to have its own set {θ i , σ i }, which thus can be optimised independently from the observations and hyper-parameters of other individuals, we identify a sum of M Gaussian log-likelihoods, log N y i , m 0 (t i ), Ψ t i θ i ,σ 2 i , and the corresponding trace terms, − 1 2 Tr( K t Ψ t i θ i ,σ 2 i −1 ).This property results in M independent maximisation problems on the corresponding functions, proving Proposition 3.2.Conversely, if we assume that all individuals in the model shares their hyper-parameters (i.e.θ, σ 2 = θ i , σ 2 i , ∀i ∈ I), we can no longer divide the problem into M sub-maximisations, and the whole sum on all individual should be optimised thanks to observations from all individuals.This case corresponds to the second part of Proposition 3.3.

Figure 2 :
Figure 2: Prediction curves (blue) of a new individual with associated 95% credible intervals (grey) for GP regression (left) and Magma (right).The dashed line represents the mean function m 0 , from the hyper-posterior p(µ 0 | {y i } i ).Observed data points are in black, and testing data points are in red.The colourful backward points are the observations from the training dataset, each colour corresponding to a different individual.

Figure 4 :1
Figure 4: MSE prediction error on the 10 last testing points with respect to the increasing number N of observed timestamps, among the first 20 points (boxplots are displayed from 100 runs in each case).

Figure 5 :
Figure 5: Prediction curves (blue) for a testing individual with associated 95% credible intervals (grey) for GP regression (left) and Magma (right), for both women (top) and men (bottom).The dashed lines represent the mean functions m 0 , from the hyper-posterior p(µ 0 | {y i } i ).Observed data points are in black, and testing data points are in red.The colourful backward points are observations from the training datasets, each colour corresponding to a different individual.

Table 1
Average MSE (standard deviation) and average CIC 95 (standard deviation) on 100 runs for GP, GPFDA and Magma.( : 99.6 (2.8), the measure of incertitude from the GPFDA package is not a genuine credible interval)

Table 2
Average MSE (standard deviation) and average CIC 95 (standard deviation) on 100 runs for the different settings of Magma.

Table 3
Average (standard deviation) training time (in seconds) for Magma and GPFDA on 100 runs for different numbers M of individuals in the training dataset.The relative running time between Magma and GPFDA is provided on the line Ratio.

Table 5
Average MSE (standard deviation) and average CIC 95 (standard deviation) for prediction on french swimmer testing datasets.