{'fraction_non_alphanumeric': 0.03626775453346932, 'fraction_numerical': 0.010833063854132982, 'mean_word_length': 4.725486007263405, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Cluster-Specific Predictions with Multi-Task Gaussian Processes

A model involving Gaussian processes (GPs) is introduced to simultaneously handle multi-task learning, clustering, and prediction for multiple functional data. This procedure acts as a model-based clustering method for functional data as well as a learning step for subsequent predictions for new tasks. The model is instantiated as a mixture of multi-task GPs with common mean processes. A variational EM algorithm is derived for dealing with the optimisation of the hyper-parameters along with the hyper-posteriors' estimation of latent variables and processes. We establish explicit formulas for integrating the mean processes and the latent clustering variables within a predictive distribution, accounting for uncertainty on both aspects. This distribution is defined as a mixture of cluster-specific GP predictions, which enhances the performances when dealing with group-structured data. The model handles irregular grid of observations and offers different hypotheses on the covariance structure for sharing additional information across tasks. The performances on both clustering and prediction tasks are assessed through various simulated scenarios and real datasets. The overall algorithm, called MagmaClust, is publicly available as an R package.


Introduction
The classic context of regression aims at inferring the underlying mapping function associating input to output data. In a probabilistic framework, some strategies assume that this function is drawn from a prior Gaussian process (GP). According to Rasmussen and Williams (2006), a Gaussian process can be defined as a collection of random variables (indexed over a continuum), any finite number of which having a joint Gaussian distribution. From this definition, we may enforce some properties for the target function solely by characterising the mean and covariance parameters of the process. Since GPs are an example of kernel methods, a broad range of assumptions can be expressed through the definition of the covariance function. We refer to Duvenaud (2014) for a comprehensive review. Despite undeniable advantages, natural implementations for GPs scale cubically with the number of data points, which constitutes a major drawback in many applications. Thereby, the early literature focused on deriving tractable approximations to mitigate this problem (Schwaighofer et al., 2004;Snelson and Ghahramani, 2005;Titsias, 2009;Hensman et al., 2013a). Subsequent reviews (Quiñonero-Candela et al., 2007;Bauer et al., 2016) have also provided standardised formulations and comparisons on this topic. Besides, several approximations have been proposed (Neal, 1997) and implemented (Rasmussen and Nickisch, 2010;Vanhatalo et al., 2013) for tackling the issue of non-Gaussian errors and adapting GPs to a broad variety of likelihoods. Since a GP corresponds to a probability distribution over a functional space, alternate approaches for modelling functional data (Ramsay and Silverman, 2005) should also be mentioned, in particular for our clustering purpose.
Unsupervised learning of functional data, also called curve clustering, focuses on the definition of sub-groups of curves, making sense according to an appropriate measure of similarity. When dealing with functional data, the concept of basis functions expansion is of paramount importance for defining smooth and manageable representations of the data. Such a representation allows the adaptation of multivariate methods such as k-means, in combination with B-splines bases for instance (Abraham et al., 2003), to handle curve clustering problems. Different bases can be used, such as Fourier or wavelets (Giacofci et al., 2013), according to the context and the nature of the signal. Besides, model-based clustering methods aims at defining probabilistic techniques for this task, and many approaches (Jiang and Serban, 2012;Jacques and Preda, 2013) have been proposed in this sense for the past decade. In particular, the algorithms funHDDC (Bouveyron and Jacques, 2011) and funFEM (Bouveyron et al., 2015) establish a mixture model where representative coefficients of the curves are supposed to come from cluster-specific Gaussian distributions. Furthermore, the authors in Schmutz et al. (2020) introduced an extension to the case of multivariate functional data. A comprehensive review (Jacques and Preda, 2014) has been proposed to discuss and compare the major approaches of this active research area. We can also mention recent works leveraging generalised Bayesian predictors and PAC-Bayes for learning and clustering streams of data (Li et al., 2018;Li and Guedj, 2021), which later inspired a work on clustering (Cohen-Addad et al., 2021). In line with the previous methods that simultaneously analyse multiple curves, we also introduce a framework that takes advantage of similarities between resembling data.
The multi-task paradigm consists in using data from several tasks (also called batches or individuals) to improve the learning or predictive capacities compared to an isolated model. It was introduced by Caruana (1997) and then adapted in many fields of machine learning. An initial GP adaptation (Schwaighofer et al., 2004) came as a hierarchical Bayesian model using an expectation-maximisation (EM) algorithm for learning, and a similar approach can be found in Shi et al. (2005). Another hierarchical formulation using a GP to model the mean parameter of another GP was later proposed in Hensman et al. (2013b). Such modelling assumptions resemble those of the present paper, although the strategies used for learning and prediction largely differ. Meanwhile, Yu et al. (2005) offered an extensive study of the relationships between the linear model and GPs to develop a multi-task formulation. More recently, the expression multi-task GP has been coined by Bonilla et al. (2008) for referring to a covariance structure involving inputs and tasks in two separate matrices. Some further developments on this approach were proposed (Hayashi et al., 2012;Rakitsch et al., 2013;Chen et al., 2020) and we can also mention the work of Swersky et al. (2013) on Bayesian hyper-parameter optimisation in such models. Generally, as presented in the reviewÁlvarez et al. (2012) which favours the term multi-output GP, all these frameworks can be expressed as specific cases of the linear model of coregionalisation (LMC) introduced by Goovaerts (1997) in the field of geostatistics. Finally, let us emphasise the algorithm Magma from Leroy et al. (2022) that recently proposed a different multi-task paradigm for GPs, by transferring information through a latent mean process rather than the covariance structure, an intuition that partially appeared before in Shi et al. (2007). This approach offers enhanced performances in forecasting while scaling linearly in the number of tasks, which is noticeably lower than the previous multi-output methods which generally bear a cubic complexity. However, the assumption of a unique mean process might happen to be too restrictive and could benefit from a more general formulation. For instance, Shi and Wang (2008) proposed an idea close to our following model by introducing a curve clustering component to a hybrid splines-GPs multi-task framework, although their approach does not fully account for uncertainty and cannot handle irregular measurements. Moreover, no implementation has been released for their algorithm, and by deriving a unified multi-task GP framework that is more general, we aim at offering to practitioners a powerful tool for tackling the simultaneous clustering and prediction of multiple functional data.
Our contributions. The present paper contributes a significant extension of Magma (Leroy et al., 2022), by introducing a clustering component into the procedure. To this end, (i) we introduce a more general model involving multiple mean GPs, each one being associated with a particular cluster. These processes represent the prior mean trend, possibly different from one cluster to another, that is associated with an individual covariance structure for each functional data. Moreover, we propose 4 different modelling hypotheses regarding the kernels' hyper-parameters of the GPs. (ii) We derive a variational expectation-maximisation (VEM) algorithm called MagmaClust (available as an R package on the CRAN and at https://github.com/ArthurLeroy/MagmaClustR) for estimating the hyper-parameters along with the hyper-posterior distributions of the mean processes and latent clustering variables. A variational BIC criterion is proposed to estimate the number of clusters. (iii) We enrich this learning step with an additional EM algorithm and analytical formulas to determine both clusters probabilities and predictive distributions for any new, partially observed, individual. The final multi-task prediction can be expressed in terms of cluster-specific distributions or as an overall GPs mixture. The algorithmic complexity of learning and prediction steps are discussed as well. (iv) We illustrate the advantages of our approach on synthetic and three real-life data sets. The results exhibit that MagmaClust outperforms state-of-the-art alternatives on both curve clustering and prediction tasks, in particular for group-structured data sets.
Outline of the paper. We introduce the multi-task Gaussian processes mixture model in Sec. 2, along with notation. Sec. 3 is devoted to the inference procedure, with a Variational Expectation-Maximisation (VEM) algorithm to estimate hyper-parameters and approximation of hyper-posterior distributions along with mixture proportions. We leverage this strategy in Sec. 4 and derive both a mixture and cluster-specific GP prediction formulas, for which we provide an analysis along with computational costs in Sec. 5. The performances of our algorithm for clustering and prediction purposes are illustrated in Sec. 6 with a series of experiments on both synthetic and real-life data sets along with comparisons to competing state-of-the-art algorithms. We close with a summary of our work in Sec. 7. All proofs to original results are deferred to Sec. 8.

Motivation
Before diving into modelling considerations, let us provide a motivational example used throughout the article to illustrate the kind of problems we expect to tackle. Assume that we have observed results of swimming competitions for thousands of athletes from 10 to 20 years old. In order to improve talent detection, one might be interested in using those data to forecast future performances for new young swimmers (e.g observed only between 10 and 14 years old). Such a data set is composed of thousands of irregular age-performance time series, where each swimmer would have specific number and locations of their data points, as illustrated in Fig. 1. These examples come from a real data set, thoroughly described in Sec. 6, which was originally the applicative motivation for developing methods described in the present paper. The following multi-task GPs framework is tailored to simultaneously allocate individuals into clusters while learning model parameters, and then provide probabilistic predictions for future performances of any young swimmer, by sharing information between resembling individuals through common mean trends.

Notation
To remain consistent with this illustrative example, we refer throughout to the input variables as timestamps and use the term individual as a synonym of batch or task. However, although the temporal formulation helps intuition, the present framework still applies to the wide range of inputs one can usually think of in GP applications. As we suppose the data set to be composed of point-wise observations from multiple functions, the set of all indices is denoted by I ⊂ N, which in particular contains {1, . . . , M }, the indices of the observed individuals (i.e. the training set). Since the input values are defined over a continuum, let us name T this input space (we can assume T ⊂ R here for simplicity). Moreover, since the following model is defined for clustering purposes, the set of indices K = {1, . . . , K} refers to the K different groups of individuals. For the sake of concision, the notation is shortened as follows: for any object x, We assume data is collected from M different sources, such that a set of N i inputoutput values constitutes the observations for the i-th individual. Below follows additional convenient notation: 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 i , the pooled set of all timestamps among individuals, • N = card(t), the total number of observed timestamps.
Let us stress that the input values may vary both in number and location among individuals, and we refer as a common grid of timestamps to the case where t i = t, ∀i ∈ I. Otherwise, we call it an uncommon grid. Besides, in order to define a GP mixture model, a latent binary random vector Z i = (Z i1 , . . . , Z iK ) needs to be associated with each individual, indicating in which cluster it belongs. Namely, if the i-th individual comes from the k-th cluster, then Z ik = 1 and 0 otherwise. Moreover, we assume these latent variables to come from the same multinomial distribution: Z i ∼ M(1, π), ∀i ∈ I, with a vector of mixing proportions π = (π 1 , . . . , π K ) and K k=1 π k = 1.
Figure 2: Graphical model of dependencies between variables in the multi-task Gaussian Processes mixture model.

Model and Assumptions
Assuming that the i-th individual belongs to the k-th cluster, we can define its functional expression as the sum of a cluster-specific mean process and an individual-specific centred process: where: • µ k (·) ∼ GP(m k (·), c γ k (·, ·)) is the common mean process of the k-th cluster, is the specific process of the i-th individual, • ε i (·) ∼ GP(0, σ 2 i I) is the error term.
This general model depends upon several mean and covariance parameters, fixed as modelling choices, and hyper-parameters to be estimated: • ∀k ∈ K, m k (·) is the prior mean function of the k-th cluster, • ∀k ∈ K, c γ k (·, ·) is the covariance kernel with hyper-parameters γ k , • ∀i ∈ I, ξ θ i (·, ·) is the covariance kernel with hyper-parameters θ i , • ∀i ∈ I, σ 2 i ∈ R is the noise variance associated with the i-th individual, • ∀i ∈ I, we define the shorthand ψ θ i ,σ 2 i (·, ·) = ξ θ i (·, ·) + σ 2 i I, the set of all hyper-parameters of the model.
One can note that we assume here the error term to be individual-specific, although we could also assume it to be cluster-specific and thus indexed by k. Such a choice would result in a valid model since the upcoming developments remain tractable if we substitute ε k to ε i everywhere, and associate σ 2 k I with c γ k (·, ·) instead of ξ θ i (·, ·). A discussion about additionally available assumptions on the covariance structures follows in Sec. 2.4. In this paper, we seek an estimation for Θ among the above quantities, whereas the other objects are pre-specified in the model. For instance, the prior mean m k (·) is usually set to zero but could also integrate expert knowledge if available. Furthermore, we assume that: • ∀i ∈ I, ∀k ∈ K, µ k , f i , Z i , ε i and all pairwise combinations are independent.
We display a graphical model on Fig. 2 to enlighten the relationships between the different components. From these hypotheses, we can naturally integrate out f i and derive the conditional prior distribution of y i (·), providing a hierarchical formulation for the model: As a consequence, the output processes {y i (·) | {Z i } i , {µ k (·)} k } i are also independent (conditionally to the latent variables) from one another. Although this model is expressed in terms of infinite-dimensional GPs, we proceed to the inference using finite-dimensional sets of observations {t i , y i } i . Therefore, we can write the joint conditional likelihood of the model (conditioning on the inputs is omitted throughout the paper for clarity): The mean processes being common to all individuals in a cluster, we need to evaluate their prior distributions on the pooled grid of timestamps t: Finally, the joint distribution of the clustering latent variables also factorises over the individuals: From all these expressions, the complete-data likelihood of the model can be derived: This expression would usually serve to estimate the hyper-parameters Θ, although it depends here on latent variables that cannot be evaluated directly. Even if the prior distributions over {Z i } i and {µ k (t)} k are independent, the expressions of their respective posteriors would inevitably depend on each other. Nevertheless, it remains possible to derive variational approximations for these distributions that still factorise nicely over the terms Z i , ∀i ∈ I, and µ k (t), ∀k ∈ K. Consequently, the following inference procedure involves a variational EM algorithm that we shall detail after a quick discussion on the optional hypotheses for the model.

Assumptions on the Covariance Structure
Throughout this paper, we detail a common ground procedure that remains consistent regardless of the covariance structure of the considered GPs. Notice that we chose a parametric distinction of the covariance kernels through the definition of hyper-parameters, different from one individual to another. However, there are no theoretical restrictions on the underlying form of the considered kernels, and we indicate a differentiation on the sole hyper-parameters merely for convenience in writing. A usual kernel in the GP literature is known as the exponentiated quadratic kernel (also called sometimes squared exponential or radial basis function kernel). This kernel only depends upon two hyper-parameters θ = {v, } such as: Although beyond the scope of the present paper, let us mention the existence of a rich literature on kernel design, properties and combinations: see Rasmussen and Williams (2006, Chapter 4) or Duvenaud (2014) for comprehensive studies.
In the algorithm Magma (Leroy et al., 2022), the multi-task aspect is mainly supported by the mean process, although the model also allows information sharing among individuals through the covariance structure. These two aspects being constructed independently, we could think of the model as potentially double multi-task, both in mean and covariance. More precisely, if we assume {{θ i } i , σ 2 i i } = {θ 0 , σ 2 0 }, ∀i ∈ I, then all f i are assumed to be different realisations of the same GP, and thus all individuals contributes to the estimation of the common hyper-parameters. Hence, such an assumption that may appear restrictive at first glance actually offers a valuable way to share common patterns between tasks. Furthermore, the same kind of hypothesis can be proposed at the cluster level with {γ k } k = γ 0 , ∀k ∈ K. In this case, we would assume that all the clusters' mean processes {µ k } k share the same covariance structure. This property would indicate that the patterns, or the variations of the curves, are expected to be roughly identical from one cluster to another and that the differentiation would be mainly due to the mean values. Conversely, different covariance structures across kernels offer additional flexibility for the clusters to differ both in position and in trend, smoothness, or any property that could be Table 1 Summary of the 4 available assumptions on the hyper-parameters, with their respective shortening notation and the associated number of sets of hyper-parameters (HPs) to optimise.
coded in a kernel. Speaking rather loosely, we may think of these different settings as a trade-off between flexibility and information sharing, or as a choice between an individual or collective modelling of the covariance. Overall, our algorithm provides 4 different settings, offering a rather wide range of assumptions for an adequate adaptation to different applicative situations. Note that the computational considerations are also of paramount importance when it comes to optimising a likelihood over a potentially high number of parameters. Hence, we display on Table 1 a summary of the 4 different settings, providing a shortening notation along with the associated number of hyper-parameters (or sets of hyper-parameters in the case of θ i and γ k ) that are required to be learnt in practice.

Inference
Although a fully Bayesian point-of-view could be taken on the learning procedure by defining prior distributions of the hyper-parameters and directly use an MCMC algorithm (Rasmussen and Williams, 2006;Yang et al., 2016) for approximate inference on the posteriors, this approach remains computationally challenging in practice. Conversely, variational methods have proved to be highly efficient to conduct inference in difficult GP problems (Titsias, 2009;Hensman et al., 2013a) and may apply in our context as well. By introducing an adequate independence assumption, we are able to derive a variational formulation leading to analytical approximations for the true hyper-posterior distributions of the latent variables. Then, these hyper-posterior updates allow the computation of a lower bound of the true log-likelihood, thereby specifying the E step of the VEM algorithm (Attias, 2000) that conducts the overall inference. Alternatively, we can maximise this lower bound with respect to the hyper-parameters in the M step for optimisation purpose, to provide estimates. By iterating on these two steps until convergence (pseudo-code in Algorithm 1), the procedure is proved to reach local optima of the lower bound (Boyd and Vandenberghe, 2004), usually in a few iterations. Regarding the previous illustrative example, this learning step allocates swimmers presenting similar progression patterns into dedicated clusters while simultaneously estimating the underlying mean trend of those clusters (E step) and optimising all hyper-parameters controlling each specific covariance structure (M step). For the sake of clarity, the shorthand Z = {Z i } i and µ = {µ k (t)} k is used in this section when referring to the corresponding set of latent variables.

Variational EM Algorithm
We seek an appropriate and analytical approximation q(Z, µ) for the exact hyper-posterior distribution p(Z, µ | {y i } i , Θ). Let us first notice that for any distribution q(Z, µ), the following decomposition holds for the observed-data log-likelihood: with: Therefore, we expressed the intractable log-likelihood of the model by introducing the Kullback-Leibler (KL) divergence between the approximation q(Z, µ) and the corresponding true distribution p(Z, µ | {y i } i , Θ). The right-hand term L(q; Θ) in (2) defines a so-called lower bound for log p({y i } i | Θ) since a KL divergence is nonnegative by definition. This lower bound depends both upon the approximate distribution q(·) and the hyper-parameters Θ, while remaining tractable under adequate assumptions. By maximising L(q; Θ) alternatively with respect to both quantities, optima for the hyper-parameters shall be reached. To achieve such a procedure, the following factorisation is assumed for the approximated distribution: Colloquially, we could say that the independence property that lacks to compute explicit hyper-posterior distributions is imposed. Such a condition restricts the family of distributions from which we choose q(·), and we now seek approximations within this family that are as close as possible to the true hyper-posteriors.

E step
In the expectation step (E step) of the VEM algorithm, the lower bound of the marginal likelihood L(q; Θ) is maximised with respect to the distribution q(·), considering that initial or previously estimated values for Θ are available. Making use of the factorised form previously assumed, we can derive analytical expressions for the optimal distributions over q Z (Z) and q µ (µ). As the computing of each distribution involves taking an expectation with respect to the other one, this suggests an iterative procedure where whether the initialisation or a previous estimation serves in the current optimisation process. Therefore, two propositions are introduced below, respectively detailing the exact derivation of the optimal distributions q Z (Z) and q µ (µ) (all proofs are deferred to the corresponding Sec. 8).
Proposition 3.1. Assume that the hyper-parameters Θ and the variational distribution where: Proposition 3.2. Assume that the hyper-parameters Θ and the variational distribution The optimal variational approximation q µ (µ) of the true hyper-posterior p µ | {y i } i , Θ factorises as a product of multivariate Gaussian distributions: with: where the following shorthand notation is used: Notice that the forced factorisation we assumed between Z and µ for approximation purpose additionally offers an induced independence between individuals as indicated by the factorisation in (3), and between clusters in (5).

M step
At this point, we have fixed an estimation for q(·) in the lower bound that shall serve to handle the maximisation of L( q, Θ) with respect to the hyper-parameters. This maximisation step (M step) depends on the initial assumptions on the generative model (Table 1), resulting in four different versions for the VEM algorithm (the E step is common to all of them, the branching point is here).
the optimal values are given by: where E {Z,µ} indicates an expectation taken with respect to q µ (µ) and q Z (Z). In particular, optimal values for π can be computed explicitly with: The remaining hyper-parameters are estimated by numerically solving the following maximisation problems, according to the situation. Let us note: Then, for hypothesis H ki : For hypothesis H k0 : For hypothesis H 0i : For hypothesis H 00 : Let us stress that, for each sub-case, explicit gradients are available for the functions to maximise, facilitating the optimisation process with gradient-based methods (Hestenes and Stiefel, 1952;Bengio, 2000). The current version of our code implements those gradients and makes use of them within the L-BFGS-B algorithm (Nocedal, 1980;Morales and Nocedal, 2011) devoted to the numerical maximisation. As previously discussed, the hypothesis H ki necessitates to learn M + K sets of hyper-parameters. However, we notice in Proposition 3.3 that the factorised forms defined as the sum of a Gaussian log-likelihoods and trace terms offer a way to operate the maximisations in parallel on simple functions. Conversely, for the hypothesis H 00 , only 2 sets of hyper-parameters need to be optimised, namely γ 0 , and {θ 0 , σ 2 0 }. The small number of functions to maximise is explained by the fact that they are defined as larger sums over all individuals (respectively all clusters). Moreover, this context highlights a multi-task pattern in covariance structures, since each individual (respectively cluster) contributes to the learning of shared hyper-parameters. In practice, H 00 is far easier to manage, and we generally reach robust optima in a few iterations. On the contrary, the settings with many hyper-parameters to learn, using mechanically less data for each, may lead more often to computational burden or pathological results. The remaining hypotheses, H 0i and H k0 , are somehow middle ground situations between the two extremes and might be used as a compromise according to the problem being dealt with.

Initialisation
Below some modelling choices are discussed, in particular the initialisation of some quantities involved in the VEM algorithm: • {m k (·)} k ; the mean functions from the hyper-prior distributions of the associated mean processes {µ k (·)} k . As it may be difficult to pre-specify meaningful values in the absence of external or expert knowledge, these values are often assumed to be 0. However, it remains possible to integrate information in the model by this mean. However, as exhibited in Proposition 3.2, the influence of {m k (·)} k in hyper-posterior computations decreases rapidly when M grows in a multi-task framework.
• {γ k } k , {θ i } i and σ 2 i i ; the kernel hyper-parameters. We already discussed that the form itself of kernels has to be chosen as well, but once set, we would advise initiating {γ k } k and {θ i } i with close and reasonable values whenever possible. As usual in GP models, nearly singular covariance matrices and numerical instability may occur for pathological initialisations, in particular for the hypotheses, like H ki , with many hyper-parameters to learn. This behaviour frequently occurs in the GP framework, and one way to handle this issue is to add a so-called jitter term (Bernardo et al., 1998) on the diagonal of the ill-defined covariance matrices.
• {τ ik } ik ; the estimated individual membership probabilities (or π; the prior vector of clusters' proportions). Both quantities are valid initialisation depending on whether we start the VEM iterations by an E step or an M step. If we only want to set the initial proportions of each cluster in the absence of additional information, we may merely specify π and start with an E step. Otherwise, if we insert the results from a previous clustering algorithm as an initialisation, the probabilities τ ik for each individual and cluster can be fully specified before proceeding to an M step (or to the q µ (µ)'s computing and then the M step).
Let us finally stress that the convergence (to local maxima) of VEM algorithms partly depends on these initialisations. Different strategies have been proposed in the literature to manage this issue, among which simulated annealing (Ueda and Nakano, 1998) or repeated short runs (Biernacki et al., 2003).

Pseudocode
The overall algorithm is called MagmaClust (as an extension of the algorithm Magma to cluster-specific mean GPs) and we provide below the pseudo-code summarising the inference procedure. The corresponding R code is available at https://github.com/ ArthurLeroy/MAGMAclust.

Model Selection
The question of the adequate choice of K in clustering applications is a recurrent concern in practice. Many criteria have been introduced in the literature, among which those relying on penalisation of the likelihood like AIC (Akaike, 1974) or BIC (Schwarz, 1978) for instance. Whereas we seek a BIC-like formula, we can recall that the likelihood p({y i } i | Θ) cannot be computed directly in the present context. However, as for inference, we may still use the previously introduced lower bound L( q; Θ) to adapt a so-called variational-BIC (VBIC) quantity to maximise, as proposed in You et al. (2014). The expression of this criterion is provided below while we defer the full derivation of the lower bound to Sec. 8.4.
Proposition 3.4. After convergence of the VEM algorithm, a variational-BIC expression can be derived as: where: • α i is the number of hyper-parameters from the individual processes' kernels, • α k is the number of hyper-parameters from the mean processes' kernels, • K − 1 is the number of free parameters π k (because of the constraint K k=1 π k = 1 ).
Let us mention that the numbers α i and α k in the penalty term vary according to the considered modelling hypothesis (H 00 , H k0 , H 0i or H ki ), see Table 1 for details.

Prediction
At this point, we would consider that the inference on the model is completed, since the training data set of observed individuals {y i } i enabled to estimate the desired hyperparameters and the distributions of latent variables. For the sake of concision, we thus omit the writing of conditionings over Θ in the sequel. Recalling our illustrative example, we would have used competition's results over a long period from thousands of swimmers for training the model, and we now expect to make predictions of future performances Table 2 Summary of the different steps to perform in the prediction procedure, according to the model assumptions and the target grid of timestamps.
for any young athlete in the early stages of their career. Therefore, let us assume the partial observation of a new individual, denoted by the index * , for whom we collected a few data points y * (t * ) at timestamps t * . Defining a multi-task GPs mixture prediction consists in seeking an analytical distribution p(y * (·) | y * (t * ), {y i } i ), according to the information brought by: its own observations; the training data set; the cluster structure among individuals. As we aim at studying the output values y * (·) at arbitrarily chosen timestamps, say t p (the index p stands for prediction), we propose a new notation for the pooled vector of timestamps t p * = t p t * . This vector serves as a working grid on which the different distributions involved in the prediction procedure are evaluated. In the absence of external restrictions, we would strongly advise to include the observed timestamps of all training individuals, t, within t p * , since evaluating the processes at these locations allows for sharing information across tasks. Otherwise, any data points defined on timestamps outside of the working grid would be discarded from the multi-task aspect of the model. In particular, if t p * = t, we may even use directly the variational distribution q µ (µ) computed in the VEM algorithm, and thus skip one step of the prediction procedure that is described below. Throughout the section, we aim at defining a probabilistic prediction for this new individual, accounting for the information of all training data {y i } i . To this end, we manipulate several distributions of the type p(· | {y i } i ) and refer to them with the adjective multi-task. Additionally to highlighting the information-sharing aspect across individuals, this term allows us to distinguish the role of {y i } i from the one of the newly observed data y * (t * ), which are now the reference data for establishing if a distribution is called a prior or a posterior. Deriving a predictive distribution in our multi-task GP framework requires to complete the following steps.

Compute the hyper-posterior approximation of {µ
5. Deduce the multi-task GPs mixture prediction: We already discussed the influence of the initial modelling hypotheses on the overall procedure. Hence, let us display in Table 2 a quick reminder helping to keep track of which steps need to be performed in each context.

Posterior Inference on the Mean Processes
In order to integrate the information contained in the shared mean processes, we first need to re-compute the variational approximation of {µ k (·)} k 's hyper-posterior on the new N -dimensional working grid t p * . By using once more Proposition 3.2, it appears straightforward to derive this quantity that still factorises as a product of Gaussian distributions where we merely substitute the values of timestamps: with: where the following shorthand notation is used: We acknowledge that the subsequent analytical developments party rely on this variational approximate distribution q µ ({µ k (t p * )} k ), and may thus be considered, in a sense, as approximated as well. However, this quantity provides a valuable closed-form expression that we can substitute to the true hyper-posterior in Proposition 4.1 below.

Computation of the Multi-Task Prior Distributions
For a sake of completeness, let us recall the equivalence between two ways of writing conditional distributions that are used in the subsequent results: We may regularly substitute one to the other in the sequel depending on the handier in the context. Once the mean processes' distributions are re-computed on the working grid, their underlying influence shall be directly plugged into a marginalised multi-task prior over y * (t p * ) by integrating out the {µ k (t p * )} k . As the mean processes vanish, the new individual's outputs y * (t p * ) directly depends upon the training data set {y i } i , as highlighted in the proposition below.
Proposition 4.1. For a set of timestamps t p * , the multi-task prior distribution of y * knowing its clustering latent variable is given by: Proof.Let us recall that, conditionally to their mean process, the individuals are independent of one another. Then, for all k ∈ K, we have: The final line is obtained by remarking that such a convolution of Gaussian distributions remains Gaussian as well (Bishop, 2006, Chapter 2), and we refer to Leroy et al. (2022) for the detailed calculus in this exact context. Therefore, we finally get:

Optimisation of the New Hyper-Parameters and Computation of the Clusters' Probabilities
Now that the mean processes have been removed at the previous step, this section strongly resembles the classical learning procedure through an EM algorithm for a Gaussian mixture model. In our case, it allows us both to estimate the hyper-parameters of the new individual {θ * , σ * } and to compute the hyper-posterior distribution of its latent clustering variable Z * , which provides the associated clusters' membership probabilities τ * . As before, E steps and M steps are alternatively processed until convergence, but this time by working with exact formulations instead of variational ones.

E step
In the E step, hyper-parameters estimates are assumed to be known. Recalling that the latent clustering variable Z * is independent from the training data {y i } i , the multi-task hyper-posterior distribution maintains an explicit derivation: By inspection, we recognise the form of a multinomial distribution and thus retrieve the corresponding normalisation constant to deduce: with:

M step
Assuming to know the value of τ * , we may derive optimal values for the hyperparameters of the new individual through the following maximisation: Let us note L * (θ * , σ * ) = log p(y * (t * ), Z * | {y i } i , θ * , σ * , π). By remarking that π has already been estimated previously, we may easily derive the expression to maximise with respect to θ * and σ * in practice: where C 1 is a constant term. Thus, the optimisation in this case merely relies on the maximisation of a weighted sum of Gaussian log-likelihoods, for which gradients are wellknown.

3bis.
In the case where the hyper-parameters are supposed to be common across individuals (H 00 or H k0 ), there is no need to additional optimisation since we already have θ * = θ 0 and σ 2 * = σ 0 by definition. However, the probabilities of lying in each cluster τ * for the new individual still need to be computed, which shall be handled by directly using the expression (8) from the E step.

3ter.
Conversely, even if hyper-parameters for each individual are supposed to be different (H 0i or H ki ), it remains possible to avoid the implementation of an EM algorithm by stating τ * = π. Such an assumption intuitively expresses that we would guess the membership probabilities of each cluster from the previously estimated mixing proportions, without taking new individual's observations into account. Although we would not recommend this choice for getting optimal results, it still seems to be worth a mention for applications with a compelling need to avoid EM's extra computations during the prediction process.

Computation of the Multi-Task Posterior Distributions
Once the needed hyper-parameters have been estimated and the prior distribution established, the classical formula for GP predictions can be applied to the new individual, for each possible latent cluster. First, let us recall the prior distribution by separating observed from target timestamps, and introducing a shorthand notation for the covariance: and likewise for the other blocks of the matrices. Therefore, recalling that conditioning on the sub-vector of observed values y * (t * ) maintains a Gaussian distribution (Bishop, 2006;Rasmussen and Williams, 2006), we can derive the multi-task posterior distribution for each latent cluster: where:

Computation of the Multi-Task GPs Mixture Prediction
To conclude, by summing over all possible combinations for the latent clustering variable Z * , we can derive the final predictive distribution.
Proposition 4.2. The multi-task GPs mixture posterior distribution for y * (t p ) has the following form:

Proof.
Taking advantage of (9) and the multi-task hyper-posterior distribution of Z * as computed in (7), it is straightforward to integrate out the latent clustering variable: where we recall for the transition to the last line that Z * k = 1 if the * -th individual belongs to the k-th cluster and Z * k = 0 otherwise. Hence, summing a product with only one non-zero exponent over all possible combination for Z * is equivalent to merely sum over the values of k, and the variable Z * k simply vanishes.

Alternative predictions
Even though Proposition 4.2 provides an elegant probabilistic prediction in terms of GPs mixture, it remains important to notice that this quantity is no longer a Gaussian distribution. In particular, the distribution of an output value at any point-wise evaluation is expected to differ significantly from a classical Gaussian variable, by being multi-modal for instance. This property is especially true for individuals with high uncertainty about the clusters they probably belong to, whereas the distribution would be close to the a when τ * k ≈ 1 for one cluster and almost zero for the others. While we believe that such a GPs mixture distribution highlights the uncertainty resulting from a possible cluster structure in data and offers a rather original view on the matter of GP predictions, some applications may suffer from this non-Gaussian final distribution. Fortunately, it remains pretty straightforward to proceed to a simplification of the clustering inference by assuming that the * -individual only belongs to its more probable cluster, which is equivalent to postulate max{τ * k } k = 1 and the others to be zero. In this case, the final Gaussian mixture turns back into a Gaussian distribution, and we retrieve a uni-modal prediction, easily displayed by its mean along with credible intervals.

Complexity Analysis for Training and Prediction
It is customary to stress that computational complexity is of paramount importance in GP models as a consequence of their usual cubic (resp. quadratic) cost in the number of data points for learning (resp. prediction). In the case of MagmaClust, we use information from M individuals scattered into K clusters, each of them providing N i observations, and those quantities mainly specify the overall complexity of the algorithm. Moreover, N refers to the number of distinct timestamps (i.e. N ≤ M i=1 N i ) in the training data set and corresponds to the dimension of the objects involved in the kernel-specific mean processes computations. Typically, the learning complexity would be proportional to one iteration of the VEM algorithm, which requires O M × N 3 i + K × N 3 operations. Let us stress that this complexity is linear in the number of tasks, which is significantly lower than classical multi-output GP approaches in the literature. As detailed iń Alvarez et al. (2012), algorithms that can be formulated as linear models of coregionalisations typically present a O(M 3 N 3 ) training complexity when not resorting to sparse approximations. Several approximations have been developed to reduce this cost, for instance, by using pseudo-inputs to decrease the N 3 term, while lowering the M 3 term generally requires underlying independence assumptionsÁlvarez and Lawrence (2011). In the Experiments section, the computational advantage of MagmaClust is also empirically highlighted in Table 5.
From a practical point of view, the hypotheses formulated on the hyper-parameters may influence the constant of this complexity but generally not in more than an order of magnitude. For instance, the models under the assumption H 00 usually require less optimisation time in practice, although it does not change the number or the dimensions of the covariance matrices to inverse, which mainly control the overall computing time. The dominating terms in this expression depend on the context, regarding the relative values of M , N i , N and K. In contexts where the number of individuals M dominates, like with small common grids of timestamps for instance, the left-hand term would control the complexity, and clustering's additional cost would be negligible. Conversely, for a relatively low number of individuals or a large size N for the pooled grid of timestamps, the right-hand term becomes the primary burden, and the computing time increases proportionally to the number of clusters compared to the original Magma algorithm.
During the prediction step, the re-computation of {µ k (·)} k 's variational distributions implies K inversions of covariance matrices with dimensions depending on the size of the prediction grid t p * . In practice though, if we fix a fine grid of target timestamps in advance, this operation can be assimilated to the learning step. In this case, the prediction complexity remains at most in the same order as the usual learning for a single-task GP, that is O(K × N 3 * ) (this corresponds to the estimation of the new individual's hyperparameters, and would decrease to O(K × N 2 * ) for H k0 or H 00 ). We shall mention that the definition of a fine grid is generally desirable only for low-dimensional applications since we may quickly reach running time or memory limits as the input's dimension grows. In many contexts, most of the time-consuming learning steps can be performed in advance, and the immediate prediction cost for each new individual is negligible in comparison (generally comparable to a single-task GP prediction).

Experiments
The present section is dedicated to the evaluation of MagmaClust on both synthetic and real data sets. The performance of the algorithm is assessed in regards to its clustering and forecast abilities. To this purpose, we introduce below the simulation scheme generating the synthetic data along with the measures used to compare our method to alternatives quantitatively. Throughout, the exponentiated quadratic (EQ) kernel, as defined in Eq. (1), serves as covariance structure for both generating data and modelling. The manipulation of more sophisticated kernels remains a topic beyond the scope of the present paper, and the EQ proposes a fair common ground for comparison between methods. Thereby, each kernel introduced in the sequel is associated with two hyper-parameters. Namely, v ∈ R + represents a variance term whereas ∈ R + specifies the length-scale. The synthetic data sets are generated following the general procedure below, with minor modifications according to the model assumptions H 00 , H k0 , H 0i or H ki : 1. Define a random working grid t ⊂ [ 0, 10 ] of N = 200 timestamps to study M = 50 individuals, scattered into K clusters, 2. Draw the prior mean functions for {µ k (·)} k : m k (t) = at + b, ∀t ∈ t, ∀k ∈ K, where a ∈ [ −2, 2 ] and b ∈ [ 20, 30 ], 3. Draw uniformly hyper-parameters for {µ k (·)} k 's kernels : where v γ k ∈ 1, e 3 and γ k ∈ 1, e 1 , (or γ 0 = {v γ 0 , γ 0 }),

For all i ∈ I, draw uniformly the hyper-parameters for individual kernels
where v θ i ∈ 1, e 3 , θ i ∈ 1, e 1 , and σ 2 i ∈ [ 0, 0.1 ], (or θ 0 = {v θ 0 , θ 0 } and σ 2 0 ), 6. Define π = ( 1 K , . . . , 1 K ) and draw Z i ∼ M(1, π), ∀i ∈ I, 7. For all i ∈ I and Z ik = 1, draw uniformly a random subset t i ⊂ t of N i = 30 timestamps, and draw y This procedure offers data sets for both the individuals {t i , y i } i and the underlying mean processes {t, µ k (t)} k . In the context of prediction, a new individual is generated according to the same scheme, although its first 20 data points are assumed to be observed while the remaining 10 serve as testing values. While it may be argued that this repartition 20-10 is somehow arbitrary, a more detailed analysis with changing numbers of observed points in Leroy et al. (2022) revealed a low effect on the global evaluation. Unless otherwise stated, we fix the number of clusters to be K * = 3 and the model assumption to be H 00 for generating the data. Let us recall that we provided a variational-BIC formula in Proposition 3.4 to select an appropriate number of clusters K from data. Therefore, this measure is evaluated in following experiments and used for model selection purposes in the real-life application.
Besides, the adjusted rand index (ARI) (Hubert and Arabie, 1985) is used as a measure of adequacy for comparison between the groups obtained through the clustering procedure and the true clusters that generated the data. More specifically, the ARI is defined by counting the proportions of matching pairs between groups, and a value of 1 represents a perfect correspondence. One can note that the ARI still applies when it comes to evaluating clustering partitions with different numbers of clusters. On the matter of prediction, the mean square error (MSE) between predicted means and the true values offers a measure of the average forecast performance. Formally, we define the MSE in prediction on the 10 testing points for the new individual as: Moreover, an additional measure accounting for the validity of uncertainty quantification is defined in Leroy et al. (2022) as the percentage of true data effectively lying within the 95% credible interval (CI 95 ), which is constructed from the predictive distribution. We extend here this measure to the context of GPs mixture, where CI 95 is no longer available directly (as for any multi-modal distribution). Namely, the weighted CI 95 coverage (W CIC 95 ) is defined to be: where CI k 95 represents the CI 95 computed for the k-th cluster-specific Gaussian predictive distribution (9). In the case where K = 1, i.e. a simple Gaussian instead of a GPs mixture, the W CIC 95 reduces to the previously evoked CI 95 coverage. By averaging the weighted cluster-specific CI k 95 coverage, we still obtain an adequate and comparable quantification of the uncertainty relevance for our predictions. By definition, the value of this indicator should be as close as possible to 95%. Finally, the mean functions {m k (·)} k are set to be 0 in MagmaClust, as usual for GPs, whereas the membership probabilities τ ik are initialised thanks to a preliminary k-means algorithm. Fig. 3 provides a comparison on the same data set between a classical GP regression (top), the multi-task GP algorithm Magma (middle), and the multi-task GPs mixture approach MagmaClust (bottom). On each sub-graph, the plain blue line represents the mean parameter from the predictive distribution, and the grey shaded area covers the CI 95 . The dashed lines stand for the multi-task prior mean functions { m k (·)} k resulting from the estimation of the mean processes. The points in black are the observations for the new individual * , whereas the red points constitute the true target values to forecast. Moreover, the colourful background points depict the data of the training individuals, which we colour according to their true cluster in MagmaClust displays (bottom). For the sake of fairness, the prior mean for standard GP has been set as the average of training data points. As expected, a simple GP regression provides an adequate fit close to the data points before quickly diving to the prior value 0 when lacking information. Conversely, Magma takes advantage of its multi-task component to share knowledge across individuals by estimating a more relevant mean process. However, this unique mean process appears unable to account for the clear group structure, although adequately recovering the dispersion of the data. In the case of MagmaClust, we display the cluster-specific prediction (9) for the most probable group instead of the GPs mixture prediction, since max k (τ * ) ≈ 1 in this example. The model selection method based on maximum VBIC values correctly retrieved the true number of cluster K = 3. To illustrate the training procedure dynamics in this example, we provide on Fig. 4 a tracking, for each iteration of the VEM algorithm until convergence, of the evidence lower bound (ELBO) and the corresponding ARI between predicted and true clusters. As the ELBO increases, notice that the ARI also improves to finally reach 1 at convergence, which means we fully retrieved the true clusters at the end of training. Although this simple example only depicts a single-run instance, it provides an accurate intuition on the general behaviour of MagmaClust. In practice, the algorithm quickly improves the clustering structure and associated mean processes during the first two steps, and generally converges in a handful of iterations.

Illustration on Synthetic Examples
Overall, this illustrative example highlights the benefit we can get from considering group-structured similarities between individuals in GP predictions. Notice that our method offers both a significant improvement in mean prediction and a narrowed uncertainty around this value. Additionally, we display on Fig. 5 the specific predictions according to the two remaining clusters (although associated with nearly 0 probabilities). We remark that the predictions move towards the cluster specific mean processes as soon as the observations become too distant. In this idealised example, we displayed Gaussian predictive distributions for convenience though, in general, a Gaussian mixture might rarely be unimodal. Therefore, we propose in Fig. 6 another example with a higher variance and overlapping groups, where the VBIC still provides the correct number of clusters. While the ARI between predicted and true clusters was equal to 1 (perfect match) in the previous example, it now decreases to 0.78. Moreover, the vector of membership probabilities associated with the Fig. 6 for the predicted individual happens to be: τ * = (0.95, 0.05, 0). The left-hand graph provides an illustration of the predictive mean, acquired from the multi-task GPs mixture distribution described in Proposition 4.2. We may notice that this curve lies very close to one cluster's mean although not completely overlapping it, because of the τ * k = 0.05 probability for another cluster, which slightly pulls the prediction onto its own mean. Besides, the right-hand graph of Fig. 6 proposes a representation of the multi-task GPs mixture distribution as a heatmap of probabilities for the location of  Figure 5: Cluster-specific prediction curves (blue) with associated 95% credible intervals (grey) from MagmaClust, for two unlikely clusters. The dashed lines represent the mean parameters from the mean processes estimates. Observed data points are in black, testing data points are in red. Backward points are the observations from the training data set, coloured by clusters. Figure 6: Left: GPs mixture mean prediction curve (blue) from MagmaClust. Right: heatmap of probabilities for the GPs mixture predictive distribution from Mag-maClust. The dashed lines represent the mean parameters from the mean processes estimates. Observed data points are in black, testing data points are in red. Backward points are the observations from the training data set, coloured by clusters.
our predictions. This way, we can display, even in this multi-modal context, a thorough visual quantification for both the dispersion of the predicted values and the confidence we may grant to each of them.
Finally, let us propose on Fig. 7 an illustration of the capacity of MagmaClust to retrieve the shape of the underlying mean processes, by plotting their estimations { m k (·)} k (dotted lines) along with the true curves (plain coloured lines) generated by the simulation scheme. The ability to perform this task generally depends on the structure of the data as well as on the initialisation, although we may observe satisfactory results both on the previous fuzzy example (left) and on a well-separated case (right). It should be noted that the mean processes' estimations also come with uncertainty quantification, albeit not displayed on Fig. 7 for the sake of clarity.

Clustering Performance
Generally, many curve clustering methods struggle to handle irregularly observed data directly. Therefore, for the sake of fairness and to avoid introducing too many smoothing biases in alternative methods, the data sets used in the following are sampled on regular grids, although MagmaClust can deal with irregular measurements by construction (see also Leroy et al. (2022) for an empirical study). The competing algorithms are the Bsplines expansion associated with a kmeans algorithm proposed in Abraham et al. (2003), funHDDC (Bouveyron and Jacques, 2011), and funFEM (Bouveyron et al., 2015). The two latter methods were introduced to handle curve clustering problems, by taking advantage of a functional latent mixture modelling, and demonstrated their ability in several applications (Schmutz et al., 2020). A naive multivariate kmeans is used as initialisation for funHDDC, funFEM, and MagmaClust. We propose on Fig. 8 an evaluation of each algorithm in terms of ARI over 100 data sets, simulated from various schemes. First, the procedure detailed in Sec. 6 is applied for each of the 4 different hypotheses (H ki , H k0 , H 0i , H 00 ) to generate data in accordance with our modelling assumptions. Additionally, we propose an alternative simulation scheme, inspired by Bouveyron et al.
On another aspect, Fig. 9 provides some insights into the robustness of MagmaClust to an incorrect setting of K, the number of clusters. For 100 data sets with a true value K * = 3, the ARI has been computed between the true partitions and the ones estimated by MagmaClust initialised with different settings K = 2, . . . , 10. Except for K = 2 where the low number of clusters prevents from getting enough matching pairs by definition, we may notice relatively unaffected performances as K increases. Despite a nonnegligible variance in results, the partitions remain consistent overall, and the clustering performances of MagmaClust seem pretty robust to misspecification of K. Table 3 Percentage of data sets for which the true number of cluster is K * , and the number of cluster selected by the VBIC is K. A total of 50 data sets were simulated for different values K * = 1, . . . , 5, and MagmaClust was tested on each with varying settings K = 1, . . . , 6, to retain the configuration that reaches the highest VBIC value. To remain on the matter of clusters' number, the model selection abilities of the proposed VBIC (Proposition 3.4) maximisation procedure are assessed on simulated data. For different true numbers of clusters, 50 data sets were simulated according to the previous scheme, and MagmaClust was run successively on each with different settings of K. The setting that reaches the highest value of VBIC is selected as the best model. The percentages of selection for each true K * and the corresponding values K retained through VBIC are reported in Table 3. The procedure seems to adequately select the number of clusters in most context, with results that deteriorate as K grows and a tendency to overpenalise, which is a classical behaviour with BIC in practice (Weakliem, 2016). As the typical BIC formula relies on asymptotic approximations, we also ran the simulation for different numbers of individuals M = 50 and M = 100. It may be noticed that the VBIC performs better as M increases, as expected. Such a property appears reassuring since the following real-data applications involves data sets gathering around M = 10 3 individuals. Another piece of evidence for the robustness to a wrong selection of K is highlighted by Table 4 in the context of forecasting. The predictive aspect of MagmaClust remains the main purpose of the method and its performances of this task partly rely on the adequate clustering of the individuals. It may be noticed on Table 4 that both MSE and W CIC 95 regularly but slowly deteriorate as we move away from the true value of K. However, the performances remain of the same order, and we may still be confident about the predictions obtained through a misspecified running of MagmaClust. In particular, the values of MSE happen to be even better when setting K = 4, . . . , 6 (we recall that the same 100 data sets are used in all cases, which can thus be readily compared). Besides, the righthand part of the table provides indications on the relative time (in seconds) that is needed to train the model for one data set and to make predictions. As expected, both training and prediction times increase roughly linearly with the values of K, which seems consistent with the complexities exposed in Sec. 5. This property is a consequence of the extra mean processes and hyper-parameters that need to be estimated as K grows. Nonetheless, the influence of K is lesser on the prediction time, which yet remains negligible, even when computing many group-specific predictions.

Prediction Performance
To pursue the matter of prediction, we provide on Table 5 the comparison of forecasting performances between our method and several state-of-the art alternatives. We use the classical (i.e. single-task) GP regression as benchmark, as we would expect all the competing multi-tasks methods to perform better, considering the additional information that can be shared between individuals. The Magma algorithm (Leroy et al., 2022) (which is equivalent to MagmaClust in the specific case where K = 1) is also evaluated to measure the accuracy gain we can achieve thanks to the additional clustering on the simulated group-structured data sets. As detailed in Sec. 1, the existing multi-task approaches in the GP literature generally consider specific kernel structures accounting for explicit correlations both between inputs and outputs (Álvarez et al., 2012). The Multi-Output Gaussian Process Tool Kit (MOGPTK) (de Wolff et al., 2021) is a Python package that implements the main multi-output covariance kernels from literature in a unified framework. We relied on this package to run experiments, first applying the proposed combination (SM LCM) of the spectral mixture (SM) (Wilson and Adams, 2013) and the linear model of coregionalisation (LMC) (Goovaerts, 1997), which is the historical and very general formulation for multi-output kernels. Additionally, we made use of the more recent multi-Output spectral mixture (MOSM) algorithm (Parra and Tobar, 2017), which is defined as the default method in the package. For the sake of fair comparison, the L-BFGS-B algorithm (Nocedal, 1980;Morales and Nocedal, 2011) is the optimisation procedure used in all competing methods.
Regarding both mean prediction and uncertainty quantification, our approach outperforms the alternatives. In terms of MSE, MagmaClust takes advantage of its multiple mean processes to results in an order of magnitude enhancement compared to the best competitors (Magma and SM LMC). As expected, the simple GP regression performs rather poorly and surprisingly MOSM results are even worse, as in practice, its current implementation seems to reach pathological cases during training most of the time. Additionally, the empirical quantification of uncertainty of MagmaClust appears very convincing since there are on average exactly 95% of the observations lying within the weighted CI 95 , as expected.
In accordance with the theoretical complexity, the increase in training times displayed in Table 5 remains roughly proportional to the value of K (we recall that MagmaClust assumes K = 3 here, compared to Magma which is K = 1). In contrast, the multioutput GPs methods (SM LMC and MOSM) multiply tenfold both training and prediction times even in these reasonable experiments (50 individuals and 30 timestamps). This cost comes from the cubic complexity in the number of tasks (due to size-augmented covariance matrices) that results in a massive computational burden as M increases. As a consequence, in the following real data applications, where thousands of individuals are considered, these methods would merely be unusable in practice.

Application to real data sets
To evaluate the efficiency of our approach in real-life applications, we introduce 3 distinct data sets corresponding to various contexts: sports performances, weight follow-up during childhood, and missing data reconstruction in CO 2 emissions. The common aspect in all these problems lies in the presence of time series collected from multiple sources, possibly with irregular measurements, for which we expect to provide accurate forecasts by taking advantage of shared information and clustered structures in the data. For all experiments, the individuals (or countries for CO 2 ) are split into training and testing sets (in proportions 60% − 40%). In the absence of expert knowledge, the prior mean functions {m k (·)} k are set to be constant equal to 0. The hypothesis H 00 is specified along with random initialisations for the hyper-parameters. The hyper-parameters, the mean processes and the cluster's membership probabilities are learnt on the training data set. Then, the data points of each testing individual are split for evaluation purposes between observed (the first 60%) and testing values (the remaining 40%). Therefore, each new process y * (·) associated with a test individual is assumed to be partially observed, and its testing values are used to compute MSE and W CIC 95 for the predictions. As measuring clustering performances directly for real-life applications is a vain effort (Kleinberg, 2002), the results are provided for several values of K, from K = 1 (i.e. Magma) to K = 6, to evaluate the effect of clustered-data structures on predictions. In all the following experiments, predictive performances tend to reach a plateau as we increase the number of clusters and no substantial improvement is noticeable for K ≥ 7.

Talent identification in competitive swimming
As previously presented in Sec. 2.1, the 100m freestyle swimming data sets initially proposed in Leroy et al. (2018) and Leroy et al. (2022) is analysed below in the new light of MagmaClust. The data sets contain results from 1731 women and 7876 men, members of the French swimming federation, 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 fol- lowing, 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)). The analysis focuses on the youth period, from 10 to 20 years, where the progression is the most noticeable. Let us recall that we aim at modelling a curve of progression from competition results for each individual in order to forecast their future performances. We expect MagmaClust to provide relevant predictions by taking advantage of both its multi-task feature and the clustered structure of data, previously highlighted in Leroy et al. (2018).
As exhibited by Table 6, MagmaClust offers excellent performances on both data sets and slightly improves Magma predictions, as we increase the number of clusters. Values of both MSE and W CIC 95 appear satisfactory in all cases, and cluster-specific predictions provide additional accuracy though one may fairly argue that the gain remains moderate overall. One of the explaining reasons is highlighted in Fig. 10, where we displayed illustrative predictions for a random man and woman when K = 5. Although we can notice clear distinctions between the different mean curves at young ages, these differences tend to decrease afterwards, as adults' performances lie in narrow intervals, especially in regards to the overall signal-on-noise ratio. Nevertheless, MagmaClust provides several additional insights into this problem compared to Magma.
First, the clusters offer interesting results to assess the profile of young swimmers and to determine the individuals to whom they most resemble. In particular, it is also possible to differentiate future evolutions associated with each cluster, along with their probabilities to occur (we do not display all the cluster-specific predictions here for the sake of concision). On the other hand, our method leads to tighter predictive distributions in terms of uncertainty. Compared to Magma which uses all training data identically, MagmaClust somehow discards the superfluous information, through the weights τ * k , to only retain the most relevant cluster-specific mean processes.

Follow-up of children's weight in Singapore
In contrast with the previous application, we now study individual time series that are almost similar at young ages before diverging while growing older. This data set (collected through the GUSTO program, see https://www.gusto.sg/) corresponds to a weight follow-up of 342 children from Singapore at 11 timestamps between birth and 72 months. In this experiment, the goal is to predict the weight of a child at {24, 36, 48, 60, 72} months, using its observed weight at {0, 3, 6, 9, 12, 18} months and data from the training individ- Figure 10: Top: women data set. Bottom: men data set. Heatmap of probabilities obtained through the GPs mixture predictive distribution of MagmaClust with K = 5 for a random illustrative swimmer. The dashed lines represent the mean parameters from the mean processes estimates. Observed data points are in black, testing data points are in yellow. Backward points are a sample of observations from the training data set, coloured according to their most probable cluster. uals. Since the weight differences between toddlers are shallow until 18 months, providing accurate long-term forecasts, while clear morphology differences emerge, seems particularly challenging at first sight. However, MagmaClust still achieves impressive performances in this context as highlighted by Table 7. Once again, the accuracy of predictions seems to slightly improve as we increase the number of clusters. In this application, recovering an adequate cluster for a child is essential to anticipate which future weight pattern is the most likely, and assuming that more clusters exist appears to help in identifying specific patterns more precisely.
This being said, each new cluster increases the overall complexity of the model and whether it is worth adding one cluster regarding the relative gain in accuracy remains a practitioner's choice. For instance, the VBIC measure we proposed in Sec. 3.4 indicates K = 3 as the optimal number of clusters for this data set (keep in mind that this criterion maximises a penalised ELBO and tells us nothing regarding the predicting abilities). Using this value for K allows us to display on Fig. 11 the behaviour of MagmaClust predictions for a random testing child. Notice that the prediction remains nicely accurate all along even though the testing points (in yellow) are not close to any of the mean processes. This particularity comes from the learning of cluster weights τ * k at young ages, where the algorithm identifies that this child belonged in nearly equal proportions to the clusters K1 (in red) and K2 (in blue). Therefore, the multi-task GPs mixture posterior distribution defined in Proposition 4.2 defines a weighted trade-off between the two mean processes that remarkably predicts the true weight values of the child for almost 5 years. Although this nice behaviour on a single example does prevent pathological cases to occur in general, we know from Table 7 that those remain rare in practice, and that the computed credible intervals encompass true data with the correct degree of uncertainty.

Reconstructing missing data in CO 2 emissions
For this last application, we propose to use MagmaClust to tackle a different kind of problem, namely, missing data reconstruction. Contrarily to the previous data sets focused on forecasting, the present collection of time series consists in historical measurements of CO 2 emissions per capita for each country from 1750 to 2020 (freely available at https: //github.com/owid/co2-data). Naturally, most countries in the world have not collected such data regularly and many annual observations are missing, especially as we move back in the past (for instance, only Canada and the United Kingdom reported values before 1800). However, our method seems particularly well-suited to recover probable distributions of historical CO 2 emissions by exploiting similarities between countries and transferring knowledge from the densely observed time series to those that remain sparse. Figure 11: Heatmap of probabilities obtained through the GPs mixture predictive distribution of MagmaClust with K = 5 for a random illustrative swimmer. The dashed lines represent the mean parameters from the mean processes estimates. Observed data points are in black, testing data points are in yellow. Backward points are observations from the training data set, coloured according to their most probable cluster. To illustrate this ambition, we provide on Fig. 12 a visualisation of the raw observations for the first 16 countries appearing by alphabetical order in the data set (for the sake of clarity, we cannot display all countries at once). While studying per capita quantities allows us to compare countries with different populations, we can observe that CO 2 emissions can massively differ around the world, and so does the regularity of measurements. By simple inspection, some patterns seem to emerge for countries with similar topographical properties (oil and gas producing regions) or similar lifestyles (typically depending on the average wealth of populations). By applying MagmaClust on this data set, we have noticed its ability to automatically recover patterns that seem logical (although such a thing remains highly subjective). For instance, when setting K = 5, one cluster of 11 countries gathered United Kingdom, United States, Russia, France and Canada among others. Another cluster only counted Kuwait and Qatar as members, while the largest regrouped around 90 countries mainly from Africa, South Asia and South America, for which the CO2 emissions per capita have generally been really low during the past centuries. More generally, we reported the prediction performances on this dataset in Table 8. Notice that increasing the number of clusters dramatically improves accuracy in this context, as it seems clear on Fig. 12 that some countries present atypical patterns and should be treated in separated clusters to reach satisfying predictions of the missing values. Overall, these abilities to take advantage of group-structured data, even with irregular measurements, and to provide probabilistic statements, highlight once more the interest of MagmaClust to tackle various machine learning problems both in supervised and unsupervised contexts.

Discussion
We introduced a novel framework to handle clustering and regression purposes with a multi-task GPs mixture model. This approach, called MagmaClust, extends the pre-  (18) vious algorithm Magma (Leroy et al., 2022) to deal with group-structured data more efficiently. The method provides new insights on the matter of GP regression by introducing cluster-specific modelling and predictions while remaining efficiently tractable through the use of variational approximations for inference. Moreover, this nonparametric probabilistic framework accounts for uncertainty both in regards to the clustering and predictive aspects, which appears to be notable in the machine learning literature. We demonstrated the practical efficiency of MagmaClust on both synthetic and real data sets where it outperformed the alternatives, particularly in group-structured context. Even though the main concern of our method remains the predictive abilities, the clustering performances also deserve to be highlighted, compared to state-of-the-art functional clustering algorithms.
While we recall that computational cost is of paramount importance to ensure broad applicability of GP models, the present version of MagmaClust yet lacks a sparse approximation, which is not trivial to derive in this framework. However, one of the stateof-the-art sparse methods (Titsias, 2009;Bauer et al., 2016) also makes use of variational inference, both to select pseudo-inputs and learn hyper-parameters of GP models. Therefore, an interesting extension could come from simultaneously computing {µ k (·)} k 's hyperposteriors and pseudo-inputs in MagmaClust, allowing for a sparse approximation of the highest dimensional objects in our model. Besides, several additional features would be worth investigating in future work, such as the extension to non-Gaussian likelihoods or enabling online updates in the learning procedure.

Proofs
As a prerequisite, let us introduce an intermediate result that will be used several times in the proofs below.
Lemma 8.1. Let X ∈ R N be a random Gaussian vector X ∼ N (m, K), where E X denotes the expectation and V X the variance with respect to this distribution. Additionally, let b ∈ R N be an arbitrary vector and S a N × N covariance matrix. Then: Proof.
8.1 Proof of Proposition 3.1 Throughout, we note E µ the expectation with respect to the variational distribution q µ (µ). From Bishop (2006, Chapter 10), the optimal solution q Z (Z) to the variational formulation verifies: Applying Lemma 8.1 to the expectation in (10), we obtain: where (by inspection of both Gaussian and multinomial distributions): , ∀i ∈ I, ∀k ∈ K.

Proof of Proposition 3.2
Let E Z denote by the expectation with respect to the variational distribution q Z (Z). From Bishop (2006, Chapter 10), the optimal solution q µ (µ) to the variational formulation verifies: If we regroup the scalar coefficient τ ik with the covariance matrix Ψ t i θ i , σ 2 i −1 , we simply recognise two quadratic terms of Gaussian likelihoods on the variables µ k (·), although evaluated onto different sets of timestamps t and t i . By taking some writing cautions and expanding the vector-matrix products entirely, it has been proved in Leroy et al. (2022) that this expression factorises with respect to µ k (t) simply by expanding vectors y i and matrices Ψ t i θ i , σ 2 i with zeros, ∀t ∈ t, t / ∈ t i . Namely, we can notice that: • ∀i ∈ I,ỹ i = 1 [t∈t i ] × y i (t) t∈t , a N -dimensional vector, • ∀i ∈ I,Ψ i = 1 [t,t ∈t i ] × ψ θ i , σ 2 i (t, t ) t,t ∈t , a N × N matrix. Therefore: By inspection, we recognise a sum of a Gaussian log-likelihoods, which implies the underlying values of the constants. Finally: with: 8.3 Proof of Proposition 3.3 Let E Z,µ be the expectation with respect to the optimised variational distributions q Z (Z) and q µ (µ). From Bishop (2006, Chapter 10), we can figure out the optimal values for the hyper-parameters Θ by maximising the lower bound L( q; Θ) with respect to Θ: Moreover, we can develop the formulation of the lower bound by expressing the integrals as an expectation, namely E Z,µ . Recalling the complete-data likelihood analytical expression and focusing on quantities depending upon Θ, we can write: It is straightforward to see that some manipulations of linear algebra also allows the derivation of explicit gradients with respect to {γ k } k , {θ i } i and σ 2 i i . Hence, we may take advantage of efficient gradient-based methods to handle the optimisation process. We should stress that the quantity (14) is a sum on the sole values of k, whereas (15) also implies a sum on the values of i. Hence, each term of these sums involves only one hyper-parameter at a time, which thus may be optimised apart from the others. Conversely, if we assume all individuals (respectively all clusters) to share the same set of hyper-parameters, then the full sum has to be maximised upon at once. Therefore, recalling that we introduced 4 different settings according to whether we consider common or specific hyper-parameters for both clusters and individuals, we shall notice the desired maximisation problems that are induced by (14) and (15).

Proof of Proposition 3.4
Let us reconsider the expression of L( q; Θ) from the previous proof. As the model selection procedure takes place after convergence of the learning step, we can use the optimal variational approximation q Z,µ to compute the lower bound explicitly. Contrarily to the M step though, we now need to develop its full expression, and thus make use of Lemma 8.1 three times.
The result follows by considering the analogous expression L( q; Θ) in which the hyper-parameters are evaluated at their optimal value.