Faïcel Chamroukhi

Nhat Thien Pham

Van Hà Hoang

Geoffrey J Mclachlan

Functional Mixtures-of-Experts

Keywords: Mixtures-of-Experts, Functional Data Analysis, EM algorithms, Maximum likelihood Estimation, Lasso regularization, Clustering 1

We consider the statistical analysis of heterogeneous data for prediction in situations where the observations include functions, typically time series. We extend the modeling with Mixtures-of-Experts (ME), as a framework of choice in modeling heterogeneity in data for prediction with vectorial observations, to this functional data analysis context. We first present a new family of ME models, named functional ME (FME) in which the predictors are potentially noisy observations, from entire functions. Furthermore, the data generating process of the predictor and the real response, is governed by a hidden discrete variable representing an unknown partition. Second, by imposing sparsity on derivatives of the underlying functional parameters via Lasso-like regularizations, we provide sparse and interpretable functional representations of the FME models called iFME. We develop dedicated expectation-maximization algorithms for Lasso-like (EM-Lasso) regularized maximum-likelihood parameter estimation strategies to fit the models. The proposed models and algorithms are studied in simulated scenarios and in applications to two real data sets, and the obtained results demonstrate their performance in accurately capturing complex nonlinear relationships and in clustering the heterogeneous regression data.

Introduction

Mixture-of-experts (ME), introduced by [START_REF] Jacobs | Adaptive mixtures of local experts[END_REF], is a successful and flexible supervised learning architecture that allows one to efficiently represent complex non-linear relationships in observed pairs of heterogeneous data (X, Y ). The ME model relies on the divide and conquer principle, so that the response Y is gathered from soft-association of several expert responses, each targeted to a homogeneous sub-population of the heterogeneous population, given the input covariates (predictors or features) X. From the statistical modeling point of view, a ME model is an extension of the finite mixture model [START_REF] Mclachlan | Finite Mixture Models[END_REF] which explores the unconditional (mixture) distribution of a given set of the features X. It is thus more tailored to unsupervised learning than to supervised learning and it has the fully conditional mixture model of the form

ME(y|x) = K k=1 G k (x)E k (y|x)•
(1)

In model ( 1) the ME distribution of the response y given the predictors x is a conditional mixture distribution with predictor-dependent mixing weights, referred to as gating functions, G k (x), and conditional mixture components, referred to as experts E k (y|x), K being the number of experts. Mixture of experts (ME) models thus allow one to better capture more complex relationships between y and x in heterogeneous situations in non-linear regression y ∈ R, classification y ∈ {1, . . . , G}, and in clustering the data by associating each expert component to a cluster. The richness of the class of ME models in terms of conditional density approximation capabilities has been recently demonstrated by proving denseness results (Nguyen et al., 2021a[START_REF] Nguyen | Approximation results regarding the multiple-output Gaussian gated mixture of linear experts model[END_REF].

They have been investigated in their simple form, as well as in their hierarchical form [START_REF] Jordan | Hierarchical mixtures of experts and the EM algorithm[END_REF], for non-linear regression and model-based cluster and discriminant analyses and in different application domains. The inference in this case can be performed by maximum likelihood estimation (MLE) via the expectation-maximization (EM) algorithm [START_REF] Jordan | Hierarchical mixtures of experts and the EM algorithm[END_REF][START_REF] Mclachlan | The EM algorithm and extensions[END_REF][START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] or, when p is possibly larger than the sample size n, by regularized MLE via dedicated EM-Lasso algorithms as in [START_REF] Khalili | New estimation and feature selection methods in mixture-of-experts models[END_REF]; [START_REF] Montuelle | Mixture of gaussian regressions model with logistic weights, a penalized maximum likelihood approach[END_REF]; Chamroukhi et al. (2019); [START_REF] Chamroukhi | Regularized maximum likelihood estimation and feature selection in mixtures-of-experts models[END_REF]; [START_REF] Huynh | Estimation and feature selection in mixtures of generalized linear experts models[END_REF]; [START_REF] Nguyen | An l 1oracle inequality for the lasso in mixture-of-experts regression models[END_REF] which include Lasso-like penalties [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. For a more complete review of ME models, the reader is referred to [START_REF] Yuksel | Twenty years of mixture of experts[END_REF] and [START_REF] Nguyen | Practical and theoretical aspects of mixtureof-experts modeling: An overview[END_REF]. More recent theoretical results about the ME estimation and model selection for different families of ME models, can be found in [START_REF] Nguyen | A nonasymptotic model selection in block-diagonal mixture of polynomial experts models[END_REF]Nguyen et al. ( ,c, 2020)).

To the best of our knowledge, ME models have been exclusively studied in multivariate analysis when the inputs are vectors, i.e. X ∈ X = R p . However, in many problems, the predictors and/or the responses are observed from smooth functions. Indeed, in many situations, unlike in predictive and cluster analyses of multivariate and potentially highdimensional heterogeneous data, which have been studied with the ME modeling in (1), the observed data may arise from continuously observed processes, e.g. time series. Thus, a multivariate (vectorial) analysis does not allow one to enough capture the inherent functional structure of the data. In such situations, classical multivariate models are not adapted as they ignore the underlying intrinsic nature and structure of the data. Functional Data Analysis (FDA) [START_REF] Ramsay | Functional Data Analysis[END_REF][START_REF] Ferraty | Nonparametric Functional Data Analysis: Theory and Practice[END_REF] in which the individual data units are assumed to be functions, rather than vectors, offers an adapted framework to deal with continuously observed data, including in regression, classification and clustering. FDA considers the observed data as (discretized) values of smooth functions, rather than multivariate observations represented in the form of "simple" vectors.

The study of functional data has been considered in most of the statistical modeling and inference problems including regression, classification, clustering, functional graphical models [START_REF] Qiao | Functional graphical models[END_REF], among others. In regression, functional linear models have been introduced including penalized functional regression ( Élodie [START_REF] Brunel | Non-asymptotic adaptive prediction in functional linear models[END_REF][START_REF] Goldsmith | Penalized functional regression[END_REF] and in particular the FLiRTI approach, a functional linear regression constructed upon interpretable regularization [START_REF] James | Functional linear regression that's interpretable[END_REF], and more generally generalized linear models with functional predictors [START_REF] Müller | Generalized functional linear models[END_REF][START_REF] James | Generalized linear models with functional predictor variables[END_REF], which cover functional logistic regression for classification. In classification, we can also cite functional linear discriminant analysis [START_REF] James | Functional linear discriminant analysis for irregularly sampled curves[END_REF], and, as a penalized model, Lasso-regularized functional logistic regression [START_REF] Mousavi | Multinomial functional regression with wavelets and lasso penalization[END_REF]. To deal with heterogeneous functional data, the construction of mixture models with functional data analytic aspects have been introduced for model-based clustering [START_REF] Liu | Simultaneous curve registration and clustering for functional data[END_REF], Jacques and Preda (2014a)) including Lasso-regularized mixtures for functional data [START_REF] Devijver | Model-based clustering for high-dimensional data. Application to functional data[END_REF][START_REF] James | Clustering for sparsely sampled functional data[END_REF][START_REF] Jacques | Model-based clustering for multivariate functional data[END_REF]Chamroukhi and Nguyen, 2019). The resulting functional mixture models are better able to handle functional data structures compared to standard multivariate mixtures.

The problem of clustering and prediction in presence of functional observations from heterogeneous populations, leading to complex distributions, is still however less investigated. In this paper, we investigate the framework of Mixtures-of-Experts (ME) models, as models of choice in modeling heterogeneity in data for prediction and clustering with vectorial observations, and extend it to the functional data framework. The main novelty of our paper is to get interpretable results for functional ME (FME). First, the ME framework is extended to a functional data setting to learn from functional predictors, and the statistical inference in the resulting setting (ie. of the FME model) looks for estimating a sparse and interpretable FME model parameters. The key technical challenges brought by this setting are addressed by considering sparsity on the derivatives of the underlying functional parameters of the FME model, thanks to Lasso-like regularizations, solved by a dedicated EM-Lasso algorithm. To the best of our knowledge, this is we first time the ME model is constructed upon functional predictors, and provides interpretable sparse estimates.

Firstly, we introduce in Section 2 a new family of ME models, referred to as FME, to relate a functional predictor to a scalar response, and develop a dedicated EM algorithm for the maximum-likelihood parameter estimation. Secondly, to deal with potential high-dimensional setting of the introduced FME model, we develop in Section 3 a Lassoregularized approach, which consists of a penalized MLE estimation via an hybrid EM-Lasso algorithm, which integrates an optimized coordinate ascent procedure to efficiently implement the M-Step. Thirdly, we in particular present in Section 3.2 and extended FME model, which is constructed upon a sparse and highly-interpretable regularization of the functional expert and gating parameters. The resulting model, abbreviated as iFME, is fitted by regularized MLE via an dedicated EM algorithm. The developed algorithms for the two introduced ME models are applied and evaluated in Section 4 on several simulated scenarios and on real data sets, in both clustering and non-linear regression.

2 Functional Mixtures-of-Experts (FME)

We wish to derive and fit new mixture-of-experts (ME) models in presence of functional predictors, and potentially, functional responses. In this paper, we first consider ME models with a functional predictor X(•) and a real response Y where the pair arises from heterogeneous population composed of unknown K homogeneous sub-populations. To the best of our knowledge, this is the first time ME models are considered for functional data.

ME with functional predictor and scalar response

Let {X i (t), Y i } n i=1 , be a sample of n independently and identically distributed (i.i.d.) data pairs where Y i is a real-valued response and X i (t) is a functional predictor with t ∈ T ⊂ R, for example the time in time series. First, to model the conditional relationships between the continuous response Y and the functional predictor X(•), given an expert z, we formulate each expert component E z (y|x) in (1) as a functional regression model (cf. [START_REF] Müller | Generalized functional linear models[END_REF], [START_REF] James | Functional linear regression that's interpretable[END_REF]). The resulting functional expert regression model for the ith observation takes the following stochastic representation

Y i = β z i ,0 + T X i (t)β z i (t)dt + ε i , i ∈ [n], (2)
where β z i ,0 is an unknown constant intercept, β z i (t), t ∈ T is the function of unknown coefficients of functional expert z i , and ε i ∼ N (0, σ 2 z i ) are independent Gaussian errors, z i ∈ [K] being the unknown label of the expert generating the ith observation. The notation z i ∈ [K] means z i = 1, . . . , K which is used throughout this paper. In this context, the response Y is related to the entire trajectory of X(•). Let β = {β z,0 , β z (t), t ∈ T } K z=1 represents the set of unknown functional parameters for the experts network. Now consider the modeling of the gating network in the proposed FME model. As in the context of ME for vectorial data, different choices are possible to model the gating network function, typically softmax-gated or Gaussian-gated ME (e.g. see [START_REF] Nguyen | Practical and theoretical aspects of mixtureof-experts modeling: An overview[END_REF], [START_REF] Xu | An alternative model for mixtures of experts[END_REF], Chamroukhi et al. (2019)). A standard choice as in [START_REF] Jacobs | Adaptive mixtures of local experts[END_REF] to model the gating network G z (x) in ( 1) is to use the multinomial logistic (softmax) function as a distribution of the latent variable Z. In this functional data modeling context with K ≥ 2 experts, we use a multinomial logistic function as an extension of the functional logistic regression presented in [START_REF] Mousavi | Functional logistic regression: a comparison of three methods[END_REF] for linear classification. The resulting functional softmax gating network then takes the following form

π z (X(t), t ∈ T ; α) = P(Z = z|X(t), t ∈ T ; α) = exp{α z,0 + T X(t)α z (t)dt} 1 + K-1 z ′ =1 exp{α z′,0 + T X(t)α z′ (t)dt} , (3)
where α = {α z,0 , α z (t), t ∈ T } K z=1 is the set of unknown constant intercept coefficients α z,0 and functional parameters α z (t), t ∈ T for each expert z ∈ [K]. Note that model ( 3) is equivalent to assuming that each expert z is related to the entire trajectory X(•) via the following functional linear predictor for the gating network

h z (X(t), t ∈ T ; α) = log π z (X(t), t ∈ T ; α) π K (X(t), t ∈ T ; α) = α z,0 + T X(t)α z (t)dt. ( 4
)
The objective is to estimate the functional parameters α and β of the FME model defined by ( 2)-(3), from an observed sample. In this setting with functional predictors, this requires estimating a possibly infinite number of coefficients (as many as the number of temporal observations for the predictor). In order to reduce the complexity of the problem, the observed functional predictor can be projected onto a fixed number of basis functions so that we sufficiently capture enough the functional structure of the data, and sufficiently reduce enough the number of coefficients to estimate.

Smoothing representation of the functional experts

Here we consider the case of fixed design, that is, the covariates X i (t) are non-random functions. We suppose that the X i (•)'s are measured with error at any given time t. Hence, instead of observing directly X i (t), one has a noisy version of it U i (t), defined as

U i (t) = X i (t) + δ i (t), i ∈ [n],
where δ i (•) ∼ N (0, σ 2 δ ) are measurement errors assumed to be independent of the X i (•)'s and the Y i 's. Since the functional predictors X i (t) are not directly observed, we first construct an approximation of X i (t) from the noisy predictors U i (t) by projecting the latter onto a set of continuous basis functions. Let b r (t) = [b 1 (t), . . . , b r (t)] ⊤ be a r-dimensional (B-spline, Fourier, Wavelet) basis, then with r sufficiently large X i (t) can be represented as

X i (t) = r j=1 x ij b j (t) = x ⊤ i b r (t), (5)
where

x ij = T X i (t)b j (t)dt for j ∈ [r] and x i = (x i1 , . . . , x ir ) ⊤ . Since X i (t)
is not observed, the representation coefficients x ij 's are unknown. Hence we propose an unbiased estimator of x ij defined as

x ij := T U i (t)b j (t)dt.
Thus, an estimate X i (t) of X i (t) is given by

X i (t) = x ⊤ i b r (t), i ∈ [n], (6)
with

x i = ( x i1 , . . . , x ir ) ⊤ .
Similarly, to represent the regression coefficient functions

β z (•), consider a p-dimensional basis b p (t) = [b 1 (t), b 2 (t), . . . , b p (t)] ⊤ .
Then the function β z (t) can be represented as

β z (t) = η ⊤ z b p (t) (7)
where η z = (η z,1 , η z,2 , . . . , η z,p ) ⊤ is the vector of unknown coefficients and the choice of p should ensure the tradeoff between smoothness of the functional predictor and complexity of the estimation problem. We select r ≥ p to satisfy the identifiability constraint (see for instance [START_REF] Goldsmith | Penalized functional regression[END_REF], [START_REF] Ramsay | Applied Functional Data Analysis: Methods and Case Studies[END_REF]). Furthermore, rather than assuming a perfect fit of β z (t) by b p (t) as in ( 7), we use for each Gaussian expert regressor z, the following error model as proposed by [START_REF] James | Functional linear regression that's interpretable[END_REF] for functional linear regression

β z (t) = η ⊤ z b p (t) + e(t)
where e(t) represents the approximation error of β z (t) by the linear projection (7). As we choose p ≫ n, |e(t)| can be assumed to be small.

Smoothing representation of the functional gating network

Since here we are examining functional predictors, an appropriate representation has also to be given for the gating network (3) with functional parameters {α z (t), t ∈ T } K z=1 . Due to the infinite number of these parameters, we also represent the gating network by a finite set of basis functions similarly as for the experts network. For the representation of the functional predictors X i (t), i ∈ [n], we use X i (t) established in (6). The coefficients function α z (t) is represented similarly as for the β coefficients function of the experts network, by using a q-dimensional basis b q (t) = [b 1 (t), b 2 (t), . . . , b q (t)] ⊤ , q ≤ r, via the projection

α z (t) = ζ ⊤ z b q (t), (8)
where ζ z = (ζ z1 , ζ z2 , . . . , ζ zq ) ⊤ is the vector of softmax coefficients function. Note that here we use the same type of basis functions for both representations of β z and α z , but one can use different types of bases if needed. Then, by using the representations ( 6) and (8) of X(t) and α z (t), respectively, in the linear predictor h z (•) defined in (4) for i ∈ [n], the latter is thus approximated as

h z (U i (t), t ∈ T ; α) = α z i ,0 + ζ ⊤ z i r i , (9)
where

r i = T b r (t)b q (t) ⊤ dt ⊤ x i .
Thus, following its definition in (3), the functional softmax gating network is approximated as -1) is the unknown parameter vector of the functional gating network to be estimated.

π k (r i ; ξ) = exp {α k,0 + ζ ⊤ k r i } 1 + K-1 k ′ =1 exp {α k ′ ,0 + ζ ⊤ k ′ r i } , ( 10
)
where ξ = (α 1,0 , ζ ⊤ 1 ), . . . , (α K-1,0 , ζ ⊤ K-1 ) ⊤ ∈ R (q+1)(K

The FME model conditional density

We now have appropriate representations for the functional predictors, as well as for both the functional gating network and the functional experts network, involved in the construction of the functional ME (FME) model ( 2)-(3). Gathering ( 6) and ( 7), the stochastic representation (2) of the FME model can thus be defined as follows,

Y i |u i (•) = β z i ,0 + η ⊤ z i x i + ε ⋆ i , i ∈ [n], (11)
where

x i = T b r (t)b p (t) ⊤ dt ⊤ x i and ε ⋆ i = ε i + x ⊤ i T b r (t)e(t)dt.
We can see that the problem is now reduced at the standard version of the ME model. From this stochastic representation under the Gaussian assumption for the error variable ε i , the conditional density of each approximated functional expert z i = k is thus given by

f y i |u i (•), z i = k; θ k = ϕ(y i ; β k,0 + η ⊤ k x i , σ 2 k ), (12)
where ϕ( • ; µ, v) is the Gaussian probability density function with mean µ and variance v, β k,0 + η ⊤ k x i is the mean of the approximated functional regression expert, σ 2 k its variance, and

θ k = (β k,0 , η ⊤ k , σ 2 k ) ⊤ ∈ R p+2
the unknown parameter vector of expert density k, k ∈ [K] to be estimated. Finally, combining ( 12) and ( 10) in the ME model (1), leads to the the following conditional density defining the FME model,

f (y i |u i (•); Ψ ) = K k=1 π k (r i ; ξ)ϕ(y i ; β k,0 + η ⊤ k x i , σ 2 k ), (13)
where Ψ = (ξ ⊤ , θ ⊤ 1 , . . . , θ ⊤ K ) ⊤ is the parameter vector of the model to be estimated.

Maximum likelihood estimation via the EM algorithm

The FME model ( 13) is now defined upon an adapted finite representation of the functional predictors, and its parameter estimation can then be performed given an observed data sample. We first consider the maximum likelihood estimation framework via the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF][START_REF] Jacobs | Adaptive mixtures of local experts[END_REF] which has many desirable properties including stability and convergence guarantees (e.g. see [START_REF] Mclachlan | The EM algorithm and extensions[END_REF] for more details). Note that here we use the term maximum likelihood estimation to not unduly clutter the clarity of the text, while as it will be specified later, we refer to the conditional maximum likelihood estimator.

In practice, the data are available in the form of discretized values of functions. The noisy functional predictors U i (t) are usually observed at discrete sampling points t i1 < . . . < t im i with t ij ∈ T for j ∈ [m i ]. We suppose that U i (t) is scaled such that 0 ≤ t ≤ 1 and divide the time period [0, 1] up into a fine grid of m i points t i1 , . . . , t im i . Thus, in (11) we have

x i = m i j=1 b r (t ij )b p (t ij ) ⊤ ⊤ x i , r i = m i j=1 b r (t ij )b q (t ij ) ⊤ ⊤ x i , where x ij = m i j=1 U i (t ij )b j (t ij ).
Note that if we choose p = q = r, then x i = r i = x i . Let D = {(u 1 , y 1 ), . . . , (u n , y n )} be an i.i.d. sample of n observed data pairs where u i = (u i,1 , . . . , u i,m i ) is the observed functional predictor for the ith response y i .

We use D to estimate the parameter vector Ψ by iteratively maximizing the observed data log-likelihood,

log L(Ψ ) = n i=1 log K k=1 π k (r i ; ξ)ϕ(y i ; β k,0 + η ⊤ k x i , σ 2 k ), (14)
via the EM algorithm. As detailed in Appendix, the EM algorithm for the FME model is implemented as follows. After starting with an initial solution Ψ (0) , it alternates, at each iteration s, between the two following steps, until convergence (when there is no longer a significant change in the values of the log-likelihood ( 14)).

E-step. Calculate the following conditional probability memberships τ (s)

ik (for all i ∈ [n]), that the observed pair (u i , y i ) originates from the kth expert, given the observed data and the current parameter estimate Ψ (s) , τ (s)

ik = P(Z i = k|y i , u i (•); Ψ (s) ) = π k (r i ; ξ (s) )ϕ(y i ; β (s) k,0 + x ⊤ i η (s) k , σ 2 k (s) ) f (y i |u i (•); Ψ (s) ) • (15)
M-step. Update the value of the parameter vector Ψ by maximizing the Q-function (43) with respect to Ψ . The maximization is performed by separate maximizations with respect to (w.r.t.) the gating network parameters ξ and, for each expert k, w.r.t. the expert network parameters θ k , for each of the K experts.

Updating the gating network's parameters ξ consists of maximizing w.r.t. ξ the part of (43) that is a function of ξ. Since we use a softmax-gated expert network in (10), this maximization problem consists of a weighted multinomial logistic problem for which there is no a closed-form solution. We then use a Newton-Raphson (NR) procedure, which iteratively maximizes (44) after starting from an initial parameter vector ξ (0) , by updating, at each NR iteration t, the values of the parameter vector ξ according to the following updating formula:

ξ (t+1) = ξ (t) -H(ξ; Ψ (s) ) -1 ξ=ξ (t) g(ξ; Ψ (s) ) ξ=ξ (t) (16)
where H(ξ; Ψ (s) ) and g(ξ; Ψ (s) ) are, respectively, the Hessian matrix and the gradient vector of Q(ξ; Ψ (s) ), and are provided in Appendix A. At each NR iteration, the Hessian matrix and gradient vector are evaluated at the current value of ξ. We keep updating the gating network parameter ξ according to (16) until there is no significant change in Q(ξ; Ψ ). The maximization then provides ξ (s+1) for the next EM iteration. Alternatively, instead of using the NR procedure, one can employ the Minorize-Maximization (MM) algorithm to update the parameter ξ of the gating network. In this approach, the Hessian matrix in ( 16) is approximated by a square matrix that is independent of ξ. This can offer an improved computational efficiency, although it comes at the cost of requiring a greater number of iterations. For further details, refer to [START_REF] Gormley | A mixture of experts model for rank data with applications in election studies[END_REF].

Updating the experts network parameters θ k consists of solving K independent weighted regression problems where the weights are the conditional expert memberships τ (s) ik given by (15). The updating formulas for the regression parameters (β k,0 , η k ) and the noise variances σ 2 k for each expert k are straightforward and correspond to weighted versions of those of standard Gaussian linear regression, i.e., weighted ordinary least squares. The updating rules for the experts network parameters are given by the following formulas:

β (s+1) k,0 = 1 n (s) k n i=1 τ (s) ik (y i -x ⊤ i η (s) k ), η (s+1) k = 1 n i=1 τ (s) ik x i x ⊤ i n i=1 τ (s) ik (y i -β (s+1) k,0 )x i , σ 2 k (s+1) = 1 n (s) k n i=1 τ (s) ik y i -(β (s+1) k,0 + x ⊤ i η (s+1) k ) 2 , (17)
where n

(s) k = n i=1 τ (s)
ik represents the expected cardinal number of component k. This EM algorithm, designed here for the FME that is constructed upon smoothing of the functional data, can be seen as a direct extension of the vectorized version the ME model. While it can hence be expected to provide accurate estimations as in the vector predictors setting, the number of parameters to estimate here in the case of the FME can still be high, for example when a big number of basis functions is used to have more accurate approximation of the functional predictors. In that case, it is better to regularize the maximum likelihood estimator in order to establish a compromise between the quality of fit and complexity.

Regularized maximum likelihood estimation

We rely on the LASSO [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] as a successful procedure to encourage sparse models in high-dimensional linear regression based on an ℓ 1 -penalty, and include it in this mixture of experts modeling framework for functional data. The ℓ 1 -regularized ME models have demonstrated their performance from a computational point of view [START_REF] Chamroukhi | Regularized maximum likelihood estimation and feature selection in mixtures-of-experts models[END_REF][START_REF] Huynh | Estimation and feature selection in mixtures of generalized linear experts models[END_REF] and enjoy good theoretical properties of universal approximation and consistent model selection in the high-dimensional setting. [START_REF] Nguyen | An l 1oracle inequality for the lasso in mixture-of-experts regression models[END_REF][START_REF] Nguyen | A non-asymptotic penalization criterion for model selection in mixture of experts models[END_REF].

ℓ 1 -regularization and the EM-Lasso algorithm

We propose an ℓ 1 -regularization of the observed-data log-likelihood (18) to be maximized, along with coordinate ascent algorithms to deal with the high-dimensional setting when updating the parameters within the resulting EM-Lasso algorithm. The objective function in this case is given by the following ℓ 1 -regularized observed-data log-likelihood,

L(Ψ ) = log L(Ψ ) -Pen λ,χ (Ψ ), ( 18
)
where log L(Ψ ) is the observed-data log-likelihood of Ψ defined by ( 14), and Pen λ,χ (Ψ ) is a LASSO regularization term encouraging sparsity for the expert and the gating network parameters, defined by

Pen λ,χ (Ψ ) = λ K k=1 p j=1 |η k,j | + χ K-1 k=1 q j=1 |ζ k,j |, (19)
where λ and χ are positive real values representing tuning parameters. The maximization of (18) cannot be performed in a closed form but again the EM algorithm can be adapted to iteratively solve the maximization problem. The resulting algorithm for the FME model, called EM-Lasso, takes the following form, as detailed in Appendix B. After starting with an initial solution Ψ (0) , it alternates between the two following steps, until convergence, i.e., when there is no longer a significant change in the values of the ℓ 1 -penalized log-likelihood (18).

E-step. The E-Step in this EM-Lasso algorithm is unchanged compared to the previously presented EM algorithm, and only requires the computation of the conditional expert memberships τ

(s)

ik according to (15).

M-step.

In this regularized MLE context, the parameter vector Ψ is now updated by maximizing the regularized Q-function (46), i.e.,

Ψ (s+1) = arg max Ψ Q(Ψ ; Ψ (s) ) -Pen λ,χ (Ψ ) .
This is performed by separate maximizations w.r.t. the gating network parameters ξ and, for each expert k, w.r.t. the expert network parameters

θ k , k ∈ [K].
Updating the gating network parameters at iteration s of the EM-Lasso algorithm consists of maximizing the following regularized Q-function w.r.t. ξ,

Q χ (ξ; Ψ (s) ) = Q(ξ; Ψ (s) ) -χ K-1 k=1 ∥ζ k ∥ 1 , (20)
where

Q(ξ; Ψ (s) ) = n i=1 K-1 k=1 τ (s) ik α k,0 + ζ ⊤ k r i -log 1 + K-1 k ′ =1 exp{α k ′ ,0 + ζ ⊤ k ′ r i } .
One can see this is equivalent to solving a weighted regularized multinomial logistic regression problem for which Q χ (ξ; Ψ (s) ) is its penalized log-likelihood, the weights being the conditional probabilities τ

(s)

ik . There is no closed-form solution for this kind of problem. We then use an iterative optimization algorithm to seek for a maximizer of Q χ (ξ; Ψ (s) ), i.e., an update for the parameters of the gating network. To be effective when the number of parameters to estimate is high, we propose a coordinate ascent algorithm to update the softmax gating network coefficients in this regularized context.

Coordinate ascent for updating the gating network. The idea of the coordinate ascent algorithm (e.g. see [START_REF] Hastie | Statistical learning with sparsity: the lasso and generalizations[END_REF], [START_REF] Huynh | Estimation and feature selection in mixtures of generalized linear experts models[END_REF]), implemented in our context at the sth EM-Lasso iteration to maximize Q χ (ξ; Ψ (s) ) at the M-Step, is as follows. The gating function parameter vectors ξ k = (α k,0 , ζ ⊤ k ) ⊤ as components of the whole gating network parameters ξ = (ξ ⊤ 1 , . . . , ξ ⊤ K-1 ) ⊤ , are updated one at a time, while fixing the other gate's parameters to their previous estimates. Furthermore, to update a single gating parameter vector ξ k , we only update its coefficients ξ kj one at a time, while fixing the other coefficients to their previous values. More precisely, for each single gating function k, we partially approximate the smooth part of Q χ (ξ; Ψ (s) ) with respect to ξ k at the current EM-Lasso estimate, say ξ (t) , then optimize the resulting objective function (with respect to ξ k ). This corresponds to solving penalized weighted least squares problems using coordinate ascent. Thus, this results into an inner loop, indexed by m, that cycles over the components of ξ k and updates them one by one, while the others are kept fixed to their previous values, i.e., ξ

(m+1) kh = ξ (m)
kh for all h ̸ = j, until the objective function ( 48) is not significantly increased. The obtained closed form updates for each coefficient ζ kj , j ∈ [q], and for the intercept α k,0 , are as follows

ζ (m+1) kj = S χ n i=1 w ik r ij (c ik - c(m) ikj ) n i=1 w ik r 2 ij for j ∈ [q], α (m+1) k,0 = n i=1 w ik (c ik -r ⊤ i ζ (m+1) k ) n i=1 w ik , where c(m) ikj = α (m) k0 + r ⊤ i ζ (m) k -ζ (m)
kj r ij is the fitted value excluding the contribution from the jth component of the ith vector r ij in the design matrix of the gating network and S χ (•) is a soft-thresholding operator defined by S χ (u) = sign(u)(|u| -χ) + with (v) + is a shorthand for max{v, 0}. The values (α

(m+1) k,0 , ζ (m+1) k
) obtained at convergence of the coordinate ascent inner loop for the kth gating function are taken as the fixed values of that gating function, in the procedure of updating the next parameter vector ξ k+1 . Finally, when all the gating functions have their values updated, i.e., after K -1 inner coordinate ascent loops, to avoid numerical instability, we perform a backtracking line search, before actually updating the gating network's parameters for the next EM-Lasso iteration. More precisely, the update is ξ (t+1) = (1 -ν)ξ (t) + ν ξ(t) , where ξ(t) is the output after K -1 inner loops and ν is backtrackingly determined to ensure

Q χ (ξ (t+1) ; Ψ (s) ) ≥ Q χ (ξ (t) ; Ψ (s) ).
We keep cycling these updated iterates for the parameter vectors ξ k , until convergence of the whole coordinate ascent (CA) procedure inside the M-Step, i.e., when the relative increase in the Lasso-regularized objective Q χ (ξ; Ψ (s) ) related to the softmax gating network is not significant, e.g., less than a small tolerance. Then, the next EM-Lasso iteration is calculated with the updated gating network's parameters ξ

(s+1) = ( α 1,0 , ζ ⊤ 1 , . . . , α K-1,0 , ζ ⊤ K-1 )
⊤ where the values α k,0 and ζ kj for all k ∈ [K -1], j ∈ [q] are those obtained for the α k,0 's and the ζ kj 's at convergence of the CA algorithm.

Updating the experts network parameters The maximization step for updating the expert parameters θ k consists of maximizing the function Q λ (θ k ; Ψ (s) ) given by

Q λ (θ k ; Ψ (s) ) = Q(θ k ; Ψ (s) ) -λ∥η k ∥ 1 , (21)
where

Q(θ k ; Ψ (s) ) = -1 2σ 2 k n i=1 τ (s) ik y i -(β k,0 + η ⊤ k x i ) 2 -n 2 log(2πσ 2 k )
• This corresponds to solving a weighted LASSO problem where the weights are the conditional experts memberships τ (s) ik . We then solve it by an iterative optimization algorithm similarly to the previous case of updating the gating network parameters. As it can be seen in Appendix B.2, updating (β k,0 , η k ) according to ( 21) is obtained by coordinate ascent as follows. For each j ∈ [p], the closed-form update for η kj at the mth iteration of the coordinate ascent algorithm within the M-Step of EM-Lasso, is given by

η (m+1) kj = S λσ 2 k (s) n i=1 τ (s) ik x ij (y i - ỹ(m) ij ) n i=1 τ (s) ik x 2 ij , β (m+1) k,0 = n i=1 τ (s) ik (y i -x ⊤ i η (m+1) k ) n i=1 τ (s) ik , in which ỹ(m) ij = β (m) k,0 + x ⊤ i η (m) k -η (m)
kj x ij is the fitted value excluding the contribution from x ij and S λσ 2 k (s) (•) is the soft-thresholding operator. We keep updating the components of (β k,0 , η k ) cyclically until no enough increase in objective function (21). Then, once (β k,0 , η k ) are updated while fixing the variance σ 2 k , the latter is then updated straightforwardly as in the case of standard weighted Gaussian regression, and it is update is given by

σ 2 k (s+1) = n i=1 τ (s) ik y i -β (s+1) k,0 -x ⊤ i η (s+1) k 2 n i=1 τ (s) ik
, where (β

(s+1) k,0 , η (s+1) k ) = ( β k,0 , η k )
is the solution obtained at convergence of the CA algorithm, which is taken as the update in the next EM-Lasso iteration.

This completes the parameter vector update

Ψ (s+1) = ξ (s+1) , θ (s+1) 1 , . . . , θ (s+1) K
of the regularized FME model, where ξ (s+1) and θ

(s+1) k , k ∈ [K]
, are, respectively, the updates of the gating network parameters and the experts network parameters, calculated by the coordinate ascent algorithms.

The EM-Lasso algorithm provides an estimate of the FME parameters with sparsity constraints on the values of the parameter vectors ξ and θ k , k ∈ [K]. Actually, since here these parameter vectors do not relate directly the original functional inputs, to the output, assuming some of their values is zero is not necessarily justified, as there is no indeed any reason that these values are zero, nor easily interpretable, compared to the sparsity in vectorial generalized linear models, mixture of regressions and ME models.

From now on, we refer to FME and FME-Lasso, respectively, the FME model fitted by EM algorithm in Section 2.5 and the regularized FME model fitted by EM-Lasso algorithm, in Section 3.1. In the following Section, we introduce a regularization that is interpretable and encourages sparsity in the FME model.

Interpretable Functional Mixture of Experts (iFME)

In this section, we provide a sparse and, especially highly-interpretable fit, for the coefficient functions {β k (t), t ∈ T } and {α k (t), t ∈ T } representing each of the K functional experts and gating network. We call our approach Interpretable Functional Mixture of Experts (iFME). The presented iFME allows us to control the expert and gating parameter functions while still providing performance as with the standard FME model presented previously.

Interpretable sparse regularization

We rely on the methodology of Functional Linear Regression That's Interpretable (FLiRTI) presented in [START_REF] James | Functional linear regression that's interpretable[END_REF]. The idea of the FLiRTI methodology is as follows. We use variable selection with sparsity assumption on appropriate chosen derivatives of the coefficient function, say β z i (t) here, in the case of the functional expert network, to produce a highly interpretable estimate for the coefficient functions β z i (t). For instance, β (0) z i (t) = 0 implies that the predictor X i (t) has no effect on the response Y i at t, β

(1)

z i (t) = 0 means that β z i (t) is constant in t, β (2)
z i (t) = 0 shows that β z i (t) is linear in t, etc. Assuming sparsity in higher-order derivatives of β z i (t), for instance when d = 3 or d = 4, will however give us a less easily interpretable fit. Hence, for example, if one believes that the expert parameter function β z i (t) is exactly zero over a certain region and exactly linear over other region of t, then it is necessary to estimate β z i (t) such that β (0) z i (t) = 0 and β

(2) z i (t) = 0 over those regions, respectively. In this situation, we need to model β z i (t) assuming that its zeroth and second derivatives are sparse. However, with the EM-Lasso method derived above via the Lasso regularization, there is no actually any reason that we could obtain those desired properties, which may result in an estimate for β z i (t) that is rarely exactly zeros (and/or linear), and making the sparsity and coefficient curves hard to interpret. The same situation may occur with the gating parameter functions. To handle this, we describe in what follows the construction of our iFME model that produces flexible-shape and highly-interpretable estimates for the expert and gating coefficient functions, by simultaneously constraining any two of their derivatives to be sparse.

We start by selecting a p-dimensional basis b p (t) and a q-dimensional basis b q (t) for approximating the experts and gating networks, respectively. For the expert network, if we divide the time domain into a grid of p evenly spaced points, and let D d be the dth finite difference operator defined recursively as

Db p (t j ) = p [b p (t j ) -b p (t j-1 )] , D 2 b p (t j ) = D [Db p (t j )] = p 2 [b p (t j ) -2b p (t j-1 ) + b p (t j-2 )] ,
. . .

D d b p (t j ) = D D d-1 b p (t j ) , then D d b p (t j ) provides an approximation for b (d) p (t j ) = [b (d) 1 (t j ), . . . , b (d)
p (t j )] ⊤ , j ∈ [p]. Let A p = D d 1 b p (t 1 ), D d 1 b p (t 2 ), . . . , D d 1 b p (t p ) A [d 1 ] p , D d 2 b p (t 1 ), D d 2 b p (t 2 ), . . . ,
γ z i = A p η z i (22) with γ z i = γ [d 1 ] z i ⊤ , γ [d 2 ] z i ⊤ ⊤ , where γ [d 1 ] z i = γ [d 1 ] 1,z i , . . . , γ [d 1 ] p,z i ⊤ , γ [d 2 ] z i = γ [d 2 ] 1,z i , . . . , γ [d 2 ]
p,z i ⊤ , and η z i defined as in relation to β z i (t) as in ( 7), then one can see that γ

[d 1 ]
z i provides an approximation to β

(d 1 ) z i (t), i.e. the d 1 th derivative of β z i (t). Respectively, γ [d 2 ] z i provides an approximation to β (d 2 ) z i (t), the d 2 th derivative of β z i (t). Therefore, enforcing sparsity in γ z i will constrain β (d 1 ) z i (t) and β (d 2 ) z i (t) to be zero at most time points. Similarly, let A q = A [d 1 ] q , A [d 2 ] q ⊤ ∈
R 2q×q be the corresponding matrix defined for the gating network. If we consider the following representation for the gating network coefficient function

ω z i = A q ζ z i (23)
with

ω z i = ω [d 1 ] z i ⊤ , ω [d 2 ] z i ⊤ ⊤
, where ω

[d 1 ] z i = ω [d 1 ] 1,z i , . . . , ω [d 1 ] q,z i ⊤ , ω [d 2 ] z i = ω [d 2 ]
1,z i , . . . , ω

[d 2 ] q,z i ⊤
, and ζ z i defined as in relation to α z i (t) as in ( 8), then we can derive the same interpretation as above for the coefficient vector ω z i and the gating parameter function α z i (t).

Using simple calculations, the representations ( 22) and ( 23) imply the following relations that subsequently used in the iFME model:

η z i = A [d 1 ] p -1 γ [d 1 ] z i , γ [d 2 ] z i = A [d 2 ] p A [d 1 ] p -1 γ [d 1 ] z i , (24)
and

ζ z i = A [d 1 ] q -1 ω [d 1 ] z i , ω [d 2 ] z i = A [d 2 ] q A [d 1 ] q -1 ω [d 1 ] z i , (25)
respectively.

In fact, one can construct A p and A q with only one derivative. Then the constraints involved to the d 2 th derivative will be eliminated making the estimation easier, but also limiting the flexibility in the shapes of the functions. That is why in this construction and in our experimental studies, A p and A q are constructed with multiple derivatives in order to produce curves of β z i (•) and α z i (•) with such many desired properties.

The iFME model

Combining the stochastic representation of the FME model in ( 11) for the experts model, the linear predictor definition in (9), and the relations ( 24)-( 25), we obtain the following iFME model construction

Y i |u i (•) = β z i ,0 + γ [d 1 ] z i ⊤ v i + ε ⋆ i , (26)
h z i (u i (•); α) = α z i ,0 + ω [d 1 ] z i ⊤ s i , (27)
subject to

γ [d 2 ] z i = A [d 2 ] p A [d 1 ] p -1 γ [d 1 ] z i and ω [d 2 ] z i = A [d 2 ] q A [d 1 ] q -1 ω [d 1 ] z i , (28)
where

v i = A [d 1 ] p -1 ⊤
x i is the new design vector for the experts and

s i = A [d 1 ] q
-1 ⊤ r i the new one for the gating network. Hence, from ( 26) and ( 27), the conditional density of each expert and the gating network are now written as

f (y i |u i (•); ψ k ) = ϕ(y i ; β k,0 + γ [d 1 ] k ⊤ v i , σ 2 k ) (29)
and

π k (s i ; w) = exp {α k,0 + ω [d 1 ] k ⊤ s i } 1 + K-1 k ′ =1 exp {α k ′ ,0 + ω [d 1 ] k ′ ⊤ s i } , ( 30
)
where

ψ k = β k,0 , γ [d 1 ] k ⊤ , σ 2 k ⊤ is the unknown parameter vector of expert component density k and w = α 1,0 , ω [d 1 ] 1 ⊤ , . . . , α K-1,0 , ω [d 1 ] K-1 ⊤ ⊤ , with α K,0 , ω [d 1 ] K
⊤ ⊤ a null vector, is the unknown parameter vector of the gating network. Finally, gathering ( 29) and (30) as for ( 13), the iFME model density is now given by

f (y i |u i (•); Ψ ) = K k=1 π k (s i ; w)ϕ(y i ; β k,0 + γ [d 1 ] k ⊤ v i , σ 2 k ), (31)
where Ψ = (w ⊤ , ψ ⊤ 1 , . . . , ψ ⊤ K ) ⊤ is the parameter vector of the model to be estimated. Thus, the iFME model constructed upon the parameter vectors γ k 's and ω k 's, for which the sparsity is assumed to obtain interpretable estimates.

Regularized MLE via the EM-iFME algorithm

In order to fit the iFME model and to maintain the sparsity in γ k and ω k , the following EM-iFME algorithm is then developed to maximize the penalized log-likelihood function

L(Ψ ) = n i=1 log f (y i |u i (•); Ψ ) + Pen λ,χ (Ψ ) (32)
with the conditional iFME density f (y i |u i (•); Ψ ) is defined in (31) and the new sparse and interpretable regularization term is given by

Pen λ,χ (Ψ ) = λ K k=1 ∥γ [d 1 ] k ∥ 1 + ρ∥γ [d 2 ] k ∥ 1 + χ K-1 k=1 ∥ω [d 1 ] k ∥ 1 + ϱ∥ω [d 2 ] k ∥ 1 , (33)
where ρ and ϱ are, respectively, the weights associated to the d 2 th derivative of the expert and the gating parameter function. The appearance of the weighting parameters ρ and ϱ, besides the usual regularization parameters λ and χ, is motivated by the fact that one may wish to place a greater emphasis on sparsity in the d 2 th derivative than in the d 1 th derivative of the parameter functions, or vice versa (we will see the usage of them in the subsequent section of experimental study). In practice, the selection of ρ and ϱ is more about whether they are greater than or less than one (i.e., the emphasis on d 2 th) rather than select an exact value.

Note that, firstly, unlike the previous FME-Lasso, in iFME model the regularization operates on the functional derivative γ k 's rather than the functional coefficients η k 's for the experts, and on the functional derivatives ω k 's rather than the functional coefficients ζ k 's for the gating network. Secondly, maximizing the penalized log-likelihood function (32) with penalization in (33) in iFME model must be coupled with the constrains (28). Follows are the two steps of the proposed EM-iFME algorithm.

E-

Step. The E-Step for the new iFME model calculates for each observation the conditional probability memberships of each expert k τ (s)

ik = P(Z i = k|y i , u i (•); Ψ (s) ) = π k (s i ; w (s) ) ϕ(y i ; β (s) k,0 + v ⊤ i γ [d 1 ] k (s) , σ 2 k (s) ) f (y i |u i (•); Ψ (s) ) , (34)
where f (y i |u i (•); Ψ (s) ) is now calculated according to the iFME density given by (31).

M-

Step. The maximization is performed by separate maximizations w.r.t. the gating network parameters w and the experts network parameters ψ k 's.

Updating the gating network parameters. The maximization step for updating the gating network parameters w = α 1,0 , ω

[d 1 ] 1 ⊤ , . . . , α K-1,0 , ω [d 1 ] K-1 ⊤ ⊤
consists of maximizing the function Q χ (w; Ψ (s) ) given by

Q χ (w; Ψ (s) ) = Q(w; Ψ (s) ) -χ K-1 k=1 ∥ω [d 1 ] k ∥ 1 + ϱ∥ω [d 2 ] k ∥ 1 , (35)
subject to

ω [d 2 ] k = A [d 2 ] q A [d 1 ] q -1 ω [d 1 ] k , ∀k ∈ [K -1], (36)
where

Q(w; Ψ (s) ) = n i=1 K-1 k=1 τ (s) ik α k,0 + ω [d 1 ] k ⊤ s i -log 1 + K-1 k ′ =1 exp{α k ′ ,0 + ω [d 1 ] k ′ ⊤ s i } •
This is a constrained version of the weighted regularized multinomial logistic regression problem, where the weights are the conditional probabilities τ (s)

ik . To solve it, in the same spirit as when updating the gating network in the previous EM-Lasso algorithm, we use an outer loop that cycles over the gating function parameters to update them one by one. However, to update a single gating function parameter

w k = (α k,0 , ω [d 1 ] k ⊤ ) ⊤ , k ∈ [K -1],
since the maximization problem ( 35) is now coupled with an additional constraint (36), rather than using a coordinate ascent algorithm as in EM-Lasso, we use the following alternative approach. For each single gating network k, using a Taylor series expansion, we partially approximate the smooth part of Q χ (w; Ψ (s) ) defined in (35) w.r.t. w k at the current estimate w (t) , then maximize the resulting objective function (w.r.t. w k ), subject to the corresponding constraint (w.r.t. k) in (36). It corresponds to solving the following penalized weighted least squares problem with constraints, max

(α k,0 , ω [d 1 ] k ,ω [d 2 ] k ) - 1 2 n i=1 w ik c ik -α k,0 -s ⊤ i ω [d 1 ] k 2 -χ ∥ω [d 1 ] k ∥ 1 + ϱ∥ω [d 2 ] k ∥ 1 subject to ω [d 2 ] k = A [d 2 ] q A [d 1 ] q -1 ω [d 1 ] k , (37)
where w ik = π k (w (t) ; s i ) 1 -π k (w (t) ; s i ) are the weights and c ik = α

(t) k,0 + s ⊤ i ω [d 1 ] k (t) + τ (s) ik -π k (w (t) ;s i ) w ik
are the working responses computed given the current estimate w (t) . This problem can be solved by quadratic programming (see, e.g., [START_REF] Gaines | Algorithms for fitting the constrained lasso[END_REF]) or by using the Dantzig selector [START_REF] Candes | The dantzig selector: Statistical estimation when p is much larger than n[END_REF], which we opt to use in our experimental studies. The details of using Dantzig selector to solve problem (37) are given in Appendix C.1.

Therefore, if ( α k,0 , ω [d 1 ] k , ω [d 2 ]
k ) is an optimal solution to problem (37), then w k = ( α k,0 , ω

[d 1 ] k ⊤ ) ⊤
is taken as an update for the gating parameter vector w k . We keep cycling over k ∈ [K -1] until there is no significant increase in the regularized Q-function (35). Updating the experts network parameters. The maximization step for updating the expert parameter vector

ψ k = (β k,0 , γ [d 1 ] k ⊤ , σ 2 k ) ⊤ consists of solving the following problem: max (β k,0 , γ [d 1 ] k ,γ [d 2 ] k ,σ 2 k ) n i=1 τ (s) ik log ϕ y i ; β k,0 + v ⊤ i γ [d 1 ] k , σ 2 k -λ ∥γ [d 1 ] k ∥ 1 + ρ∥γ [d 2 ] k ∥ 1 subject to γ [d 2 ] k = A [d 2 ] p A [d 1 ] p -1 γ [d 1 ] k . ( 38
)
As in the previous EM-Lasso algorithm, we first fix σ 2 k to its previous estimate and perform the update for (β k,0 , γ

[d 1 ]
k ), which corresponds to solving a penalized weighted least squares problem with constraints. This is be performed by using the Dantzig selector, in the same manner as previously for solving problem (37). The corresponding technical details can be found in Appendix C.2.

Once the (β k,0 , γ

[d 1 ]
k ) are updated, the straightforward update for the variance σ 2 k is given by the standard estimate of a weighted Gaussian regression. More specifically, let ( β k,0 , γ

[d 1 ] k , γ [d 2 ]
k ) be the solution to the problem (38) (with σ 2 k fixed to σ 2 k (s) ), the updates for expert parameter vector ψ k are given by (β

(s+1) k,0 , γ [d 1 ] k (s+1) ) = ( β k,0 , γ [d 1 ] k ), σ 2 k (s+1) = n i=1 τ (s) ik y i -β (s+1) k,0 -v ⊤ i γ [d 1 ] k (s+1) 2 n i=1 τ (s) ik
• Thus, at the end of the M-Step, we obtain a parameter vector update

Ψ (s+1) = (w (s+1) , ψ (s+1) 1 , . . . , ψ (s+1) K
) for the next EM iteration, where w (s+1) and ψ

(s+1) k , k ∈ [K]
, are, respectively, the updates of the gating network parameters and the experts network parameters, calculated by the two procedures described above. Alternating the E-Step and M-Step until convergence, i.e., when there is no longer a significant change in the values of the penalized log-likelihood (32), leads to a penalized maximum likelihood estimate Ψ for Ψ . Estimating the coefficient functions: Finally, once the parameter vector of iFME model has been estimated, the coefficient functions of the gating network α k (t), k ∈ [K -1] and the ones of the experts network β k (t), k ∈ [K], can be reconstructed by the following formulas

α k (t) = b q (t) ⊤ A -1 q ω [d 1 ] k , β k (t) = b p (t) ⊤ A -1 p γ [d 1 ] k , (39)
where ω

[d 1 ] k and γ [d 1 ] k
are respectively the regularized maximum likelihood estimates for ω

[d 1 ] k and γ [d 1 ]
k .

Non-linear regression and clustering with FME models

Once the model parameters have been estimated, one can then construct an estimate of the unknown functional non-linear regression function, as well as a clustering of the data into groups of similar pairs of observations. For the aim of functional non-linear regression, the unknown non-linear regression function with functional predictors is given by the following conditional expectation y|u(

•) = E[Y |U (•); Ψ ], which is defined by y i |u i (•) = K k=1 π k (r i ; ξ)( β k,0 + η ⊤ k x i
) for the FME model (13), and by

y i |u i (•) = K k=1 π k (s i ; w)( β k,0 + γ [d 1 ] k ⊤ v i
) for the iFME model (31). For clustering, a soft partition of the data into K clusters, represented by the estimated posterior probabilities τ ik = P(Z i = k|u i , y i ; Ψ ), is obtained. A hard partition can also be computed according to the Bayes' optimal allocation rule. That is, by assigning each curve to the component having the highest estimated posterior probability τ ik , defined by (15) for FME or by (34) for the iFME model, using the MLE Ψ of Ψ :

z i = arg max 1≤k≤K τ ik ( Ψ ), i ∈ [n],
where z i denotes the estimated cluster label for the ith curve.

Tuning parameters and model selection

In practice, appropriate values of the tuning parameters should be chosen. In using FME, this cover the selection of K, the number of experts, and r, p, and q, the dimensions of B-spline bases used to approximate, respectively, the predictors, the experts, and the gating network functions, although they can be chosen to be equal. For the FME-Lasso, additionally the ℓ 1 penalty constants χ and λ in (19) should be chosen. For the iFME model, the tuning parameters include also d 1 , d 2 , the two derivatives related to the sparsity constraints, and ρ, ϱ, the weights associated to the d 2 th derivative of the expert and gating functions (see ( 33)).

The selection of the tuning parameters can be performed by a cross-validation procedure with a grid search scheme to select the best combination. An alternative is to use the wellknown BIC criterion [START_REF] Schwarz | Estimating the dimension of a model[END_REF] or, in our paper, its extension, called modified BIC [START_REF] Städler | ℓ1-penalization for mixture regression models[END_REF] defined as

mBIC = L( Ψ ) -df( Ψ ) log n 2 , ( 40
)
where Ψ is the obtained log-likelihood estimator (for the FME model) or penalized loglikelihood estimator (for the FME-Lasso and iFME models), and the number of degrees of freedom df( Ψ ) is the effective number of parameters of the model, given by df( Ψ ) = df(ζ) + (K -1) + df(η) + K + K for the FME and FME-Lasso models, df(ω) + (K -1) + df(γ) + K + K for the iFME model, in which the quantities df(ζ), df(η), df(ω) and df(γ) are, respectively, the counts for nonzero free coefficients in the vectors ζ, η, ω, and γ. Note that, because of the constraints (28) for the iFME model, free coefficients in ω and γ consist of only the part related to the

d 1 derivative. That is, df(ζ) = K-1 k=1 q j=1 1 {ζ kj ̸ =0} , df(η) = K k=1 p j=1 1 {η kj ̸ =0} , df(ω) = K-1 k=1 q j=1 1 {ω [d 1 ] kj ̸ =0} and df(γ) = K k=1 p j=1 1 {γ [d 1 ] kj ̸ =0} .
We apply both the BIC and the modified BIC in our experimental study.

Experimental study

We study the performances of the FME, FME-Lasso, and iFME models in regression and clustering problems by considering simulated scenarios and real data with functional predictors and scalar responses. The interests of this study consist of the prediction performance as well as the estimation of the functional parameters, i.e., the expert and gating functions in the ME model, along with the clustering partition of the considered heterogeneous data.

Evaluation criteria

We will use the following criteria, for where applicable, to assess and compare the performances of the models and the related algorithms. For regression evaluation, we use the relative prediction error (RPE) and the correlation (Corr) index as standard criteria to evaluate the prediction performance, i.e., to quantify the relationship between the true and the predicted values of the scalar outputs. The RPE is defined by RPE = n i=1 (y i -y i ) 2 / n i=1 y 2 i , where y i and y i are, respectively, the true and the predicted response of the ith observation in the testing set. For clustering evaluation, we use the Rand index (RI), the adjusted Rand index (ARI), and the clustering error (ClusErr), as standard criteria to evaluate the clustering performance, i.e., to quantify how similar the testing observations are presented in the true partition compared to the predicted partition. To evaluate the parameters estimation performance, we compute the mean squared error (MSE) between the true and the estimated functional parameters. The MSE between a true function g(•) and its estimate g(•) is defined by

MSE g(•) = 1 m m j=1 g(t j ) -g(t j ) 2 , ( 41
)
with m being the number of time points taken into account. The function g(•) here corresponds to an expert function β(•), or a gating function α(•). The parameter functions are reconstructed from the model parameters using ( 7) and ( 8) for both FME and FME-Lasso models, and using (39) for iFME model. The values of these criteria are averaged over N Monte Carlo runs (N = 100 for simulation, for the real data, see the corresponding section). Note that, the average over N sample replicates of MSE g(•) is equivalent to the usual Mean Integrated Squared Error (MISE)

defined by MISE g(•) = E T g(t) -g(t)
2 dt , where the integral here is calculated numerically by a Riemann sum over the grid t 1 , . . . , t m .

Simulation studies

Simulation parameters and experimental protocol

We generate data from a K-component functional mixture of Gaussian experts model that relates a scalar response y ∈ R to a univariate functional predictor X(t), t ∈ T defined on a domain T ⊂ R. The data generation protocol is detailed in Appendix D.1. The parameters that were used in this data generating process (53) are described and provided in the simulation parameters and experimental protocol section in Appendix D.2. We study different simulation scenarios (sample sizes n, number of observations per time-series input m, noise levels σ 2 δ ,..) summarised in Section D.2 and Table. 9. Figure 1 displays, for each scenario, 10 randomly taken predictors colored according to their corresponding clusters. For each run, the concerned dataset is randomly split into a training set and testing set of equal size, the model parameters are estimated using training set, with the tuning parameters selected by maximizing the modified BIC (40). The evaluation criteria are computed on testing set and reported for each model accordingly.

Some implementation details

For all scenarios, for all datasets, we implemented the three proposed models with 10 EM runs and with a tolerance of 10 -6 . For the iFME model, in principle, for each parameter function the two derivatives d 1 and d 2 to be penalized, and the weights for the latter (i.e., ρ and ϱ) can

S4 predictors

Figure 1: 10 randomly taken predictors in scenarios S1 (large m and small σ δ ), S2 (small m and small σ δ ), S3 (large m and large σ δ ), and S4 (small m and large σ δ ). Here, the noisy predictors U i (•) are colored (blue, red, yellow) according to their true cluster labels Z i ∈ {1, 2, 3}.

be seen as tuning parameters. However, such an implementation could be computationally expensive in this simulation study with 400 datasets in total and 10 EM runs for each dataset. Therefore, we opted to fix d 1 and d 2 for all implementations (d 1 = 0 and d 2 = 3 for both expert and gating networks.) and left ρ and ϱ to be selected in some sets of targeted values. The choices of the targeted values were made by the following straightforward arguments. Since β 1 (•) and β 2 (•) have zero-valued regions, the weight for their zeroth derivative in penalization term should be large, equivalently, the weight for the third derivative should be small, so ρ is selected in a set of small values: ρ ∈ {10 -2 , 10 -3 , 10 -4 }. On the other hand, for α 1 (t) and α 2 (t), we select ϱ ∈ {10, 10 2 , 10 3 } as we should emphasize sparsity in their third derivative.

Simulation results

Clustering and prediction performances. We report in Table 1 the results of regression and clustering tasks on simulated datasets in the four considered scenarios. The mean and standard error of the relative predictions error (RPE) and correlation (Corr) summarize the regression performance, while the mean and standard error of the Rand index (RI), adjusted Rand index (ARI) and clustering error (ClusErr) summarize the clustering performance.

As we can see from Table 1, all the models have very good performances on both regression and clustering tasks, with high correlation, RI, ARI, and small RPE and clustering error. The iFME appears to slightly have a better performance than the others in all scenarios. The low standard errors confirm the stability of the algorithms.

Figure 2 shows the clustering results obtained by the models with highest values of the modified BIC criterion, on a dataset selected in scenario S1. against the predictors at two specific time points: t 1 = 0 and t 50 = 0.5. The highly accurate predictions (in both regression and clustering) can be seen visually through Figure 2. This figure also shows that it is difficult to cluster these data according to a few number of time observations, for example in R 2 , according to {(U i (t 0 ), y i )} n i=1 or {(U i (t 50 ), y i )} n i=1 , which suggests using functional approaches.

Comparison with functional regression mixtures (FMR): Finally, we compare our proposed models with the functional mixture regression (FMR) model proposed in [START_REF] Yao | Functional mixture regression[END_REF]. In their approach, the functional predictors are first projected onto its eigenspace, then the obtained new coordinates are fed to the standard mixture regression model to estimate the weights and the coefficients of the β(•) in that eigenspace. They performed functional principal component analysis (FPCA) to obtain estimates for the eigenfunctions and the principal component scores (which serve as predictors). The number of relevant FPCA components are chosen automatically for each dataset by selecting the minimum number of components that explain 90% of the total variation of the predictors. It is noticed that, in their simulation studies, the authors computed the scalar responses by using conditional prediction, i.e., the true y i were used to determine which cluster the observation belongs to. Then the predicted y i is calculated as the conditional mean of the density of the corresponding cluster. For comparison with that approach, we also used this strategy to make predictions in our models. We further employed the FMR model with the B-spline functional representation, instead of the functional PCA, the number of B-spline functions is set to be the same as the number of basis functions used in our models. Table 2 shows the evaluation criteria corresponding to the considered models evaluated on 100 datasets in scenario S3. Here, FMR-PC is the original model of [START_REF] Yao | Functional mixture regression[END_REF] and FMR-B is the modified one with B-spline bases. As expected, FME, FME-Lasso, and iFME, which are more flexible compared to the the FMR model, allows to capture more complexity in the Parameter estimation performance. To evaluate the parameter estimation performance, we consider the functional parameter functions estimation error as defined by ( 41). This error between the true function and the estimated one, provides an evaluation of how well the proposed models reconstruct the hidden functional gating and expert networks. In this evaluation, we considered scenario S1, i.e., m = 300, σ 2 δ = 1. Moreover, in order to provide an idea of the impact of training size on parameter estimation, we run the models with different training sizes (share the same scenario S1) and report the MSE for each function, for each model in Table 3. It shows that there are significant improvements, even with small training size, when using iFME model in estimating the gating network, compared to the FME and FME-Lasso model. Now, in order to evaluate the designed sparsity of zeroth and third derivatives of the reconstructed functions, we compute the MSEs versus their true values, i.e., zeros, on each designed intervals. In particular, we divide the domain T = [0, 1] into three parts: T 1 = [0, 0.3), T 2 = [0.3, 0.7), and T 3 = [0.7, 1], then for each model the MSEs on the corresponding intervals are reported in Table 4. 1.55e 08

β 1 (•) β 2 (•) β 3 (•) α 1 (•) α 2 (•)
on T 2 on T 1 ∪ T 3 on T β (0) 1 (•) β (0) 2 (•) β (3) 1 (•) β (3) 2 (•) β (3) 3 (•) α (3) 1 (•) α (3) 2 (•)
(1.34e 08)

1.03e 08

(1.58e 08)

8.34e 09

(1.20e 08)

Table 4: MSE of the derivatives of reconstructed functions on the corresponding interested intervals, in which,

T 1 = [0, 0.3), T 2 = [0.3, 0.7), T 3 = [0.7, 1] and T = [0, 1].
The reported values are the averaged of 100 Monte Carlo samples with standard errors in parentheses.

The reported values show that, as expected, iFME model is better than both FME and FME-Lasso in providing sparse solutions with respect to the derivatives of the parameter functions.

Table 5 shows the means and the standard errors of the estimated intercepts and variances. Note that these coefficients are not considered in the penalization. For the intercepts β k,0 , all the models estimated them very well, while for the intercepts α k,0 , iFME is slightly better than the others. For the variances, FME-Lasso gave the estimated values closest with the true values on average. The initialization is crucial for the EM algorithm. In all of our experimental studies, the model parameter vector Ψ of the FME model was initialized as follows (similar for the FME-Lasso and iFME models). Firstly, we perform a K-means algorithm on the predictors {x i } n i=1 , then for each estimated group k, (β

β 1,0 β 2,0 β 3,0 α 1,0 α 2,0 σ 1 σ 2 σ 3 True value -5 0
(0) k,0 , η (0)
k ) is initialized as the solution to the linear regression problem 4.17

y i = β k,0 + x ⊤ i η k ,

3.91

(2.08) initialized as the estimated variance within group kth. For the gating parameter, we simply drawn each component of ξ (0) randomly from the uniform distribution U(0, 1). However, to see the impact of different initialization strategies, we performed a side experiment with data taken from S1 scenario. In particular, the initialization for expert parameters is fixed as described before, while that for gating parameter varies in "rand.", "zeros", and "LR". Here, "rand." is the random strategy described above, "zeros" is the strategy where all coefficients of ξ are initialized as zeros, i.e., we are putting equal weights for the experts; and "LR" is the strategy where we perform a logistic regression with the predictors are r i (for FME, FME-Lasso) or s i (for iFME), and the responses are the labels resulted by K-means on them. Table 6 shows the time to convergence, and the log-likelihood, RPE on testing sets corresponding with the three strategies. We can see that, in terms of log-likelihood and RPE, simple strategies as "rand." and "zeros" perform quite well, while in terms of time to convergence, LR is a good choice for the iFME model.

Finally, to illustrate the selection of the number of expert components using BIC and/or modified BIC in this simulation study, we provide, in Figure 3, the plots of these criteria against the number of experts for each model. Here, we implemented the models with all fixed tuning parameters, except K which varies in the set {1, . . . , 6}. We can observe that BIC selects the correct number K = 3 for both FME-Lasso and iFME, while it selects K = 4 for FME. However, the modified BIC selects K = 4 for both FME and FME-Lasso, and selects the true number of components K = 3 for iFME.

Application to real data

In this section, we apply the FME, FME-Lasso and iFME models to two well-known real datasets, Canadian weather and Diffusion Tensor Imaging (DTI). For each dataset, we perform clustering and investigate the prediction performance, estimate the functional mixture of experts models with different number of experts K and perform the selection of K using modified BIC, and discuss the obtained results.

Canadian weather data

Canadian weather data has been introduced in [START_REF] Ramsay | Functional Data Analysis[END_REF] and is also available in the R package fda. This dataset consists of m = 365 daily temperature measure- curves (functional predictors) to predict the precipitations (scalar responses) at each station. Moreover, in addition to predicting the precipitation values, we are interested in clustering the temperature curves (therefore the stations), as well as identifying the periods of time of the year that have effect on prediction for each group of curves. Firstly, in order to assess the prediction performance of the FME, FME-Lasso and iFME models on this dataset, we implement the models by selecting the tuning parameters, including the number of expert components K in the set {1, 2, 3, 4, 5, 6}, by maximizing the modified BIC criterion, given its performance as shown in the simulation study. We report in Table 7 the results in terms of correlations, sum of squared errors (SSE) and relative prediction errors. According to the obtained results, iFME provides the best results w.r.t. all the criteria. The cross-validated RPE provided by iFME is only of 1.2%, the next is FME-Lasso with 4.5%, while the FME model has the worst RPE value. Note that in [START_REF] James | Functional linear regression that's interpretable[END_REF], the authors applied their proposed model to Canadian weather data and obtained a 10-fold cross-validated SSE of 4.77. Clearly, with the smaller cross-validated SSEs, the FME-Lasso and iFME models significantly improve the prediction.

Figure 5 shows the obtained results with the FME, FME-Lasso and iFME models, with K = 4, and with the derivatives d 1 = 0, d 2 = 3 for the iFME model. The estimated experts functions and gating functions are presented in the two top panels of the curve, while the clustering for the temperature curves and the stations are shown in the two bottom panels. As we can see, all models provide reasonable clustering for the curves which may be corresponding to different complicated underlying meteorological forecasting mechanisms. Particularly, although not using any spatial information, merely temperature information, the obtained clustering for the stations is also comparable with the original labels of the stations. For example, the FME and FME-Lasso models identify exactly the Arctic stations, while iFME identifies exactly the Pacific stations, and all of the models provide reasonable spatially organized clusters. However, what is interesting here is the shape of the expert and gating functions α(•)'s and β(•)'s obtained by the models. While FME and FME-Lasso gave less interpretable estimations, iFME appears to give, as it can be seen in the two top-right panels, piece-wise zero-valued and possibly quadratic estimated functions, which have a wide range of flat relationship from January to February and from June to September.

Motivated by the above results, on direction of identifying the periods of time of the year that truly have an effect on prediction, we implement the iFME model with K = 2, and the derivative levels d 1 and d 2 are set to be the zeroth and the third derivatives. The reason for the choices of d 1 and d 2 is that the penalization on the zeroth derivative would take into account zero ranges in the expert and gating functions, while the penalization on the third derivative, would take into account the smoothness for the changes between the periods of times in the functions. The obtained results are shown in Figure 6. As we can see, there are differences in the prediction mechanisms of the models between the northern stations and the southern stations. At southern stations, the obtained β 2 (t) shows that there is a negative relationship in the spring and a positive relationship in the late fall, but no relationship in the remaining period of the year. This phenomenon is concordant with the result obtained in [START_REF] James | Functional linear regression that's interpretable[END_REF], where the authors obtained the same relationships in the same periods of time. However, our iFME model additionally suggests that, at the northern stations, the relationship between temperature and precipitation may differ from that of southern stations. This may be explained by the differences in mean temperatures and climatic characteristics between the two regions. Finally, Figure 7 displays the values of modified BIC for varying number of expert component for the proposed models on Canadian weather data. According to these values, FME-Lasso and iFME select K = 2, while FME selects K = 3.

Diffusion tensor imaging data for multiple sclerosis subjects

We now apply our proposed models to the diffusion tensor imaging (DTI) data for subjects with multiple sclerosis (MS), discussed in [START_REF] Goldsmith | Longitudinal penalized functional regression for cognitive outcomes on neuronal tract measurements[END_REF]. The data come from a longitudinal study investigating the cerebral white matter tracts of subjects with multiple sclerosis, recruited from an outpatient neurology clinic and healthy controls. We are interested in the underlying relationship between the fractional anisotropy profile (FAP) from the corpus callosum and the paced auditory serial addition test (PASAT) score, which is a commonly used examination of cognitive function affected by MS. The FAP curves are derived from DTI data, which are obtained by a Magnetic Resonance Imaging (MRI) scanner. Each curve is recorded at 93 locations along the corpus callosum. The PASAT score is the In [START_REF] Ciarleglio | Wavelet-based scalar-on-function finite mixture regression models[END_REF], the authors applied their Wavelet-based functional mixture regression (WBFMR) model with two components to this dataset, and observed that there is one group in which there is no association between the FAP and the PASAT score for those subjects belonging to it. We accordingly fix K = 2 in our models. Figure 8 displays the obtained results for each of the three models, and Figure 9a shows the functional predictors FAP curves clustered with the iFME model. In this implementation, we tried iFME model with two different combinations of d 1 and d 2 : (d 1 , d 2 ) = (0, 2) and (d 1 , d 2 ) = (0, 3). As expected, when d 2 is the second derivative, the reconstructed parameter functions are piecewise zero and linear, while when d 2 is the third derivative, the reconstructed functions have smooth changes along the tract location.

In Figure 8, we have the three following observations. First, as it can be observed in Figure 8 right-panel, all models give a threshold of 50 that clusters the PASAT scores, this is the same as the threshold observed in [START_REF] Ciarleglio | Wavelet-based scalar-on-function finite mixture regression models[END_REF]. Second, the absolute values of the coefficient functions β 2 (t)'s are significantly smaller than those of β 1 (t)'s, this is again the same with the result obtained by the WBFMR model. Third, when β 2 (t) is estimated as zero in the FME-Lasso and iFME models, the shape of β 1 (t) is almost the same as the shape obtained in [START_REF] Ciarleglio | Wavelet-based scalar-on-function finite mixture regression models[END_REF], particularly, the peak at around the tract location of 0.42. These confirm the underlying relationship between the fractional anisotropy and the cognitive function: higher fractional anisotropy values between the locations about 0.2 to 0.7 results in higher PASAT scores for subjects in Group 1. The clustering of the FAP curves, resulted by the iFME model with d 1 = 0, d 2 = 2, is shown in Figure 9 (b)-(d).

Next, to compare with [START_REF] Ciarleglio | Wavelet-based scalar-on-function finite mixture regression models[END_REF], we investigate the prediction performance of the proposed models with respect to the leave-one-out cross validated relative prediction errors defined by CVRPE = n i=1 (y i -y

(-i) i
) 2 / n i=1 y 2 i , where y i is the true score for subject i and y

(-i) i
is the score predicted by the model fit on data without subject i. In this implementation, we keep fixing K = 2 and select the other tuning parameters by maximizing the modified BIC. The CVRPEs corresponding to the models are provided in

Estimated groups

Group 1 Group 2

Figure 8: The estimated expert and gating coefficient functions, the estimated groups of the PASAT scores, resulted by FME, FME-Lasso and iFME models with K = 2 for the DTI dataset. For iFME model, the upper is implemented with penalization on the zeroth and third derivatives, while the lower is with penalized zeroth and second derivatives.

Conclusion and discussion

The first algorithm for mixtures-of-experts constructed upon functional predictors is presented in this paper. Beside the classic maximum likelihood parameter estimation, we proposed two other regularized versions that allow sparse and notably interpretable solutions, by regularizing in particular the derivatives of the underlying functional parameters of the model, after projecting onto a set of continuous basis functions. The performances of the proposed approaches are evaluated in data prediction and clustering via experiments involving simulated and two real-world functional data.

The presented FME models can be extended in different ways. First direct extensions of the modeling framework presented here can be considered with categorical response, to perform supervised classification with functional predictors, or with vector response, to perform multivariate functional regression. Then, it may be interesting to consider the extension of the FME model to setting involving vector (or scalar) predictors and functional responses [START_REF] Chiou | Functional response models[END_REF].

Another extension, which we intend also to investigate in the future, concerns the case when we observe pairs of functional data, i.e., a sample of n functional data pairs {X i (t), Y i (u)} n i=1 , t ∈ T ⊂ R, u ∈ U ⊂ R, where Y i (•) is a functional response, explained by a functional predictor X i (•). The modeling with such FME extension then takes the form Y i (u) = β z i ,0 (u) + t X i (t)β z i (t, u)dt + ε i (u), to explain the functional response Y by the functional predictor X via the unknown discrete variable z. The particularity with this model is that, for the clustering, as well as for the prediction, we model the relation between Y at any time u and the entire curve of X, or the entire curve of each variable X ij in the case of multivariate functional predictor X i .

Since we look for sparse estimates in the space of the derivatives of the functional parameters, and provided the resulting estimates with highly structured shapes, we believe that using the fused LASSO [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF] rather than the LASSO could be expected to take accommodate more the targeted sparsity. The choice of the regularization constants for the derivatives is more about the targeted penalization magnitude (i.e., very large or very small) rather than the exact values. Typically, in the case one has no idea about the shape of the functional experts, but still want them to be interpretable, then it is practical to consider the zeroth and third derivatives (since they produce smooth changes), with the grid for ρ containing both very small and very large values. The same for the functional gating network. Furthermore, if we believe that the relation of the response on the predictor is very sparse, except over some small regions where they are linearly related, then we will penalize the zeroth and second derivatives of the functional experts, and we will put very small regularization on the second derivative. This could produce curves such as those in bottom panels of Figure 8. For a more general tuning of the regularization constant values, we recommend performing a complete cross-validation study to select the best values of the regularization constants.

with

Q(θ k ; Ψ (s) ) = - 1 2σ 2 k n i=1 τ (s) ik y i -(β k,0 + η ⊤ k x i ) 2 - n 2 log(2πσ 2 k ),
where θ k = (β k,0 , η ⊤ k , σ 2 k ) ⊤ ∈ R p+2 is the unknown vector and Ψ (s) is the current estimation of model's parameters. There is no closed-form solution for this optimization problem, we then solve it by an iterative optimization algorithm similarly to updating the gating network parameters. We first perform the update for (β k,0 , η k ) while fixing σ 2 k . This corresponds to solving a weighted LASSO problem where the weights are the the posterior experts memberships τ

•

We keep updating the components of (β k,0 , η k ) cyclically until the change in objective function of ( 49) is small enough. So, the update for (β k,0 , η k ) in this EM iteration is then (β

(s+1) k,0 , η (s+1) k ) = (β * k,0 , η * k )
where the latter is the optimal solution of (49). Finally, the update for σ 2 k is given by

σ 2 (s+1) k = n i=1 τ (s) ik (y i -β (s+1) k,0 -x ⊤ i η (s+1) k ) 2 n i=1 τ (s) ik •
Hence, the update for the value of parameters vector Ψ at M-step, i.e., the solution to problem (47), is Ψ (s+1) = (ξ (s+1) , θ ) where ξ (s+1) and θ (s+1) k , k ∈ [K], are solved by maximizing Q χ (ξ; Ψ (s) ) and Q λ (θ k ; Ψ (s) ), respectively, using the algorithms described above. The EM algorithm monotonically increases (18). Furthermore, the sequence of parameter estimates generated by the EM algorithm converges toward a local maximum of the log-likelihood function [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF][START_REF] Mclachlan | The EM algorithm and extensions[END_REF][START_REF] Wu | On the convergence properties of the em algorithm[END_REF]. Similarly to the problem in the gating network update, by decomposing γ into its positive and negative parts, γ = γ + -γ -, the above problem becomes max ( γ + , γ -)∈R 4p+2 -[0, 1 ⊤ p , ρ1 ⊤ p , 0, 1 ⊤ p , ϱ1 ⊤ p ]

γ + γ - subject to            X ⊤ στ X στ -X ⊤ στ X στ -X ⊤ στ X στ X ⊤ στ X στ γ + γ - ≤ λ + X ⊤ στ y στ λ -X ⊤ στ y στ , A -A γ + γ - = 0 4p+2 , γ + ≥ 0 2p+1 , γ -≥ 0 2p+1 , (52)
which is a standard linear program with 4q + 2 variables and can be solved similarly to the gating network case. The solution γ * = (β * 0 , γ [d 1 ] * ⊤ , γ [d 2 ] * ⊤ ) ⊤ to the original problem is retrieved by the relation γ * = γ * + -γ * -where ( γ * ⊤ + , γ * ⊤ -) ⊤ is the solution of (52). Finally, the update for σ 2 k is

σ 2 k (s+1) = n i=1 τ (s) ik (y i -β (s+1) k,0 -v ⊤ i γ [d 1 ] k (s+1) ) 2 n i=1 τ (s) ik , in which β (s+1) k,0 , γ [d 1 ] k (s+1)
are the new updates for β k,0 and γ

[d 1 ]
k .

Figure 2 :

 2 Figure 2: Scatter plots of y i against U i (t 1 ) (top panels) and U i (t 50 ) (bottom panels) on a randomly selected dataset. Here, the clustering errors are 5.5%, 4.75% and 5.0% for FME, FME-Lasso and iFME models, respectively.

Figure 3 :Figure 4 :

 34 Figure 3: Values of BIC (top) and modified BIC (bottom) for (a) FME, (b) FME-Lasso and (c) iFME model versus the number of experts K, fitted on a randomly taken dataset in scenario S1. Here, the square points correspond to the highest values.

Figure 5 :Figure 6 :

 56 Figure 5: Results obtained by FME (left panels), FME-Lasso (middle panels) and iFME (right panels) on Canadian weather data with K = 4. For each column, the panels are respectively the estimated functional experts network, estimated functional gating network, estimated clusters of the temperature curves and estimated clusters of the stations.

Figure 7 :

 7 Figure 7: Values of modified BIC for (a) FME, (b) FME-Lasso and (c) iFME model, versus the number of experts K, fitted on Canadian weather data. The square points correspond to highest values. Here, iFME is implemented with d 1 = 0, d 2 = 3 and ρ = ϱ = 100.

Figure 9 :Figure 10 :

 910 Figure 9: DTI data with (a) FAP curves of all subjects, (b) Cluster 1 and (c) Cluster 2, obtained by iFME model with d 1 = 0, d 2 = 2, and (d) the point-wise average of the curves in each of the two clusters.

  (β k,0 , η k ) has new value, the variance σ 2 k is updated straightforwardly by the standard update of a weighted Gaussian regression.More specifically, when σ 2 k , the variance of expert k, is fixed toσ 2 k (s) , updating (β k,0 , η k )consists of solving the following weighted LASSO problem:max (β k,0 ,η k ) ik y i -(β k,0 + η ⊤ k x i ) kj |, λ > 0, (49)which can be solved by coordinate ascent algorithm. For each j ∈ [p], the closed-form update for η kj is given by η kj x ij is the fitted value excluding the contribution from x ij and S χ (•) is the soft-thresholding operator (see[START_REF] Hastie | Statistical learning with sparsity: the lasso and generalizations[END_REF], sec. 5.4). Here m denotes the mth iteration of the coordinate ascent algorithm. The update for β k,0 is β

  problem with constraints. From now on, let us simplify the subscript k in the notations and rewrite the problem under matrix form as follows max to A γ = 0 p ,where γ = (β 0 , γ [d 1 ] ⊤ , γ [d 2 ] ⊤ ) ⊤ is the unknown coefficients vector, y στ = (σ √ τ 1 y 1 , . . . , σ √ τ n y n ) ⊤ ∈ R n is the weighted response vector, X στ = σ ⊙ [ √ τ |V τ |0 n×p ] ∈ R n×(2p+1) is the weighted design matrix, Λ = diag(0, 1 ⊤ p , ρ1 ⊤ p ) is the diagonal weighting matrix, n k = n i=1 τ ik and A = [p ] ∈ R p×(2p+1) is the constraints matrix. Here,√ τ = ( √ τ 1 , . . . , √ τ n ) ⊤ ∈ R n , V τ = [ √ τ 1 v 1 , . . . , √ τ n v n ] ⊤ ∈ R n×p, with v i are the design vectors (see (28)), 0 n×p ∈ R n×p contains 0's, 0 p ∈ R p contains 0's, 1 p ∈ R p contains 1's and I p is the identity matrix in R p×p . As the last term in the objective function is independent of γ, this problem is similar to the problem (50) and then can be solved by Dantzig selector estimate as the solution to the following problem max γ∈R 2p+1-∥Λ γ∥ 1 subject to |X ⊤ στ (y στ -X στ γ)| ≤ λ1 2p+1 , A γ = 0 p .

  D d 2 b p (t p )

	A [d 2 ] p

⊤

be the matrix of approximate d 1 th and d 2 th derivative of the basis b p (t). We denote A [d 1 ] p the first p rows of A p and A [d 2 ] p the remainder, i.e., A p = A [d 1 ] p , A [d 2 ] p ⊤ . One can see such matrix A p is in R 2p×p and A [d 1 ] p is a square invertible matrix. Now, if we consider the following representation for the expert network coefficient function

Table 1 :

 1 Evaluation criteria of FME, FME-Lasso and iFME models for test data in scenarios S1, . . . , S4. The reported values are the averages of 100 Monte Carlo runs with standard errors in parentheses. The bold values correspond to the best solution.

	Here, we plotted the responses

  y i against U i (t 1 ) (top panels) and U i (t 50 ) (bottom panels) on a randomly selected dataset. Here, the clustering errors are 5.5%, 4.75% and 5.0% for FME, FME-Lasso and iFME models, respectively. data, thanks to the predictor-depending mixture weights, and provides clearly better results than the FMR alternatives.

	RPE	Corr	ARI	ClusErr
	FME .0614 (.0896)	.9681 (.0469)	.8523 (.0934)	.0618 (.0490)
	FME-Lasso .0145 (.0254)	.9924 (.0133)	.9582 (.0641)	.0174 (.0325)
	iFME .0139 (.0086) .9928 (.0046) .9594 (.0430) .0151 (.0166)
	FMR-B .0460 (.0399)	.9768 (.0201)	.7064 (.1044)	.1124 (.0434)
	FMR-PC .0191 (.0331)	.9902 (.0164)	.8345 (.0794)	.0612 (.0361)

Table 2 :

 2 Performance comparison of the models for datasets in scenarios S3. The reported values are the averages of 100 Monte Carlo samples with standard errors in parentheses.

Table 3 :

 3 Average of 100 Monte Carlo runs of MSE between the estimated functions resulted by FME, FME-Lasso and iFME models in S1 scenario.

Table 5 :

 5 Intercepts and variances obtained by FME, FME-Lasso and iFME models in scenario S1. The reported values are the averages of 100 Monte Carlo samples with standard errors in parentheses.

Table 6 :

 6 Comparison of different initialization strategies. The reported values are the mean and standard error (in parentheses) over 100 Monte Carlo runs.

Table 8 .

 8 Note that, for comparison, in Ciarleglio and Ogden (2016), the CVRPE of their WBFMR model is 0.0315 and of the wavelet based functional linear model (FLM) is 0.0723.Finally, we present in Figure10the selection of the number of experts K with modified BIC. In this case, FME and FME-Lasso select K = 2, and iFME selects K = 4.

				All FAPs in dataset				53 FAPs in group 1
		0.8						0.8
	Fractional anisotropy	0.3 0.4 0.5 0.6 0.7					Fractional anisotropy	0.3 0.4 0.5 0.6 0.7
		0.2						0.2	Average profile
		0	0.2	0.4	0.6	0.8	1	0	0.2	0.4	0.6	0.8	1
				Tract location					Tract location
				(a)						(b)
				46 FAPs in group 2				Average profiles of two groups
		0.8						0.8
	Fractional anisotropy	0.4 0.5 0.6 0.7					Fractional anisotropy	0.3 0.4 0.5 0.6 0.7	Group 1
		0.3	Average profile				0.2	Group 2
		0	0.2	0.4	0.6	0.8	1	0	0.2	0.4	0.6	0.8	1
				Tract location					Tract location
				(c)						(d)

Table 8 :

 8 CVRPEs of the models on the DTI data.

Acknowledgements

This research was partly supported by an Ethel Raybould Fellowship, University of Queensland (FC), the Australian Research Council grant number DP180101192 (GJM), the French National Research Agency ANR grant SMILES ANR-18-CE40-0014 (FC and NTP), Région Normandie grant RIN AStERiCs (FC and VHH), Vietnam National University -Ho Chi Minh City (VNU-HCM) grant number C2022-18-01 (VHH), and the exploratory research facility "EXPLO" at SystemX (FC).

A EM for the FME model

The FME model can be fitted by iteratively maximizing the observed-data log-likelihood (14) iteratively via the EM algorithm. For FME, the EM takes the following form. The complete-data log-likelihood upon which the EM principle is constructed is defined by

Z ik being an indicator binary-valued variable such that Z ik = 1 if Z i = k (i.e., if the ith pair (x i , y i ) is generated from the kth expert component) and Z ik = 0, otherwise.

E-step

This step computes at each EM iteration s the expectation of the complete-data log-likelihood (42), given the observed data D, and the current parameter vector Ψ (s) :

where τ (s) ik = P(Z i = k|y i , u i (•); Ψ (s) ) is the conditional probability that the observed pair (u i (•), y i ) is generated by the kth expert. This step therefore only requires the computation of the conditional probabilities τ (s) ik as defined by (15). M-step This step updates the value of the parameter vector Ψ by maximizing the Qfunction (43) with respect to Ψ , that is Ψ (s+1) = arg max Ψ Q(Ψ ; Ψ (s) ), via separate maximizations w.r.t. the gating network parameters, and the experts network parameters as follows.

Updating the the gating network parameters Updating the the gating network's parameters ξ consists of maximizing w.r.t. ξ the following function

This consists of a weighted multinomial logistic problem for which there is no closed-form solution. This can be performed by the Newton-Raphson (NR) algorithm which iteratively maximizes (44) according to the procedure ( 16).

Let us denote by ξ 1 , . . . , ξ K-1 the parameter vectors (α 1,0 , ζ ⊤ 1 ) ⊤ , . . . , (α K-1,0 , ζ ⊤ K-1 ) ⊤ . Since there are K -1 parameter vectors to be estimated, the Hessian matrix H(ξ; Ψ (s) ) is a block-matrix, consists of (K -1) × (K -1) blocks, in which each block H kℓ (ξ; Ψ (s) ), for k, ℓ ∈ [K -1], is given by:

where δ kℓ is the Kronecker symbol (δ kℓ = 1 if k = ℓ, 0 otherwise). The gradient vector consists of K -1 gradients corresponding to the vectors ξ k , for k ∈ [K -1], and is given by

where, for k

Updating the experts network parameters Updating the experts network's parameters θ k consists of maximizing the function Q(θ k ; Ψ (s) ) given by

, consists of a weighted Gaussian regression problem where the weights are the conditional memberships τ (s) ik , and the updates are given by (17).

B EM-Lasso for ℓ 1 -regularized MLE of the FME model

The EM-Lasso algorithm for the maximization of (18) firstly requires the construction of the penalized complete-data log-likelihood

where log L c (Ψ ) is the non-regularized complete-data log-likelihood log-likelihood defined by (42). Thus, the EM algorithm for the FME model is implemented as follows. After starting with an initial solution Ψ (0) , it alternates between the two following steps, until convergence (when there is no longer a significant change in the values of the penalized log-likelihood (18)).

E-step. This step computes the expectation of the complete-data log-likelihood (45), given the observed data D, using the current parameter vector Ψ (s) :

which only requires the computation of the posterior probabilities of component membership τ

), for each of the K experts as defined by (15). M-step. This step updates the value of the parameter vector Ψ by maximizing the Qfunction ( 46) with respect to Ψ , that is, by computing the parameter vector update

The maximization is performed by separate maximizations w.r.t. the gating network parameters and the experts network parameters.

B.1 Updating the gating network parameters

Updating the gating network parameters at the sth EM iteration consists of maximizing the following function

is the gating network parameter vector and Ψ (s) is the current estimation of model's parameters. One can see this is equivalent to solving a weighted regularized multinomial logistic problem for which Q χ (ξ; Ψ (s) ) is its penalized log-likelihood. There is no closed-form solution for this kind of problem. We then use an iterative optimization algorithm to seek for a maximizer of Q χ (ξ; Ψ (s) ), i.e., an update for the parameters of the gating network. The idea is to update only a single gate at a time, while fixing other gate's parameters to their previous estimates. Again, to update that single gate, we only update one component at a time, while fixing the other components to their previous values. This procedure for updating the gating network parameters is supported by the methodology of coordinate ascent algorithm: if the objective function consists of a concave, differentiable function and a sum of concave functions then the maximizer can be achieved by iteratively maximizing with respect to each coordinate direction at a time.

Coordinate ascent for updating the gating network. Suppose at the sth EM iteration, we wish to update the gates one by one such that it maximizes Q χ (ξ; Ψ (s) ). To do that, we create an outer loop, indexed by t, which cycles over the gates. For each single gate, say gate k, we partially approximate the smooth part of Q χ (ξ; Ψ (s) ) with respect to (α k,0 , ζ k ) at ξ (t) , then optimize the obtained objective function (with respect to (α k,0 , ζ k )) by solving a penalized weighted least square problem using coordinate ascent algorithm. Note that ξ (t) denotes the current value of ξ at the iteration tth of the outer loop, while ξ (s) is the value of ξ before entering the outer loop.

In particular, using Taylor expansion, one has a quadratic approximation for smooth part of Q χ (ξ; Ψ (s) ) with respect to (α k,0 , ζ k ) at ξ (t) given by

where

and C(ξ (t) ) is a function of ξ (t) . After calculating the partial quadratic approximation l k (α k,0 , ζ k ) about the current estimator ξ (t) , we then solve the following penalized weighted least square problem max

to obtain an update for the parameters of gate k.

As mentioned above, this problem could be solved by coordinate ascent algorithm. This means we will create an inner loop, indexed by m, cycles over the components of (α k,0 , ζ k ) and update them one by one until the objective function of (48) does not gain any significant increase. For each j ∈ [q], using the soft-thresholding operator (see [START_REF] Hastie | Statistical learning with sparsity: the lasso and generalizations[END_REF], sec. 5.4), one can obtain the closed form update for ζ kj as follows

kj r ij is the fitted value excluding the contribution from r ij , S χ (•) is a soft-thresholding operator defined by S χ (u) = sign(u)(|u| -χ) + and (v) + a shorthand for max{v, 0}. Note that at each iteration of the inner loop, only one component is updated while the others are kept to their previous values, that means ζ

kh for all h ̸ = j. For α k,0 , the closed-form update is given by

Once the inner loop converges, the new values of (α k,0 , ζ k ) are used for the updating procedure of the next gate. When all the gates have their new values, i.e., after K -1 inner loops, we perform a backtracking line search before actually updating the gating network's parameters for the next t-indexed iteration. More precisely, the update is

where ξ(t) is the output after K -1 inner loops and ν is backtrackingly determined to ensure

We keep running the t-indexed loop until convergence, i.e., when there is no significant relative variation in Q χ (ξ; Ψ (s) ). Once α k,0 and ζ kj reach their optimal values α k,0 and ζ kj for all k ∈ [K -1], j ∈ [q], the update for gating network's parameters is then ξ

B.2 Updating the experts network parameters

The maximization step for updating the expert parameters θ k consists of maximizing the function Q λ (θ k ; Ψ (s) ) given by

C EM-iFME for updating iFME model parameters C.1 Updating the gating network parameters

This section presents the using of Dantzig selector to solve problem (37). Let us simplify the subscript k in the notations and rewrite the problem under matrix form as follows max

where

is the weighted working response vector,

is the diagonal weighting matrix and

s i are the design vectors (see ( 28)), 0 n×q ∈ R n×q contains 0's, 0 q ∈ R q contains 0's, 1 q ∈ R q contains 1's and I q is the identity matrix in R q×q . This problem can be viewed as the problem of finding a sparse solution via Lasso estimate for a linear regression model with constraints. Therefore, we can solve it alternatively by Dantzig selector estimate as the solution to the following problem max

where the absolute value operator is understood componentwise. By decomposing ω into its positive and negative parts, ω = ω + -ω -, the above problem becomes

which is a standard linear program with 4q + 2 variables, therefore can be easily solved by available toolboxes for linear programming, e.g., Matlab's linprog function. Finally, the solution

to the original problem could be retrieved by the relation

⊤ is the solution of (51).

C.2 Updating the expert network parameters

This section presents how to solve problem (38). Firstly, we fix σ 2 k to its previous estimate and perform an update for (β k,0 , γ

k ), it corresponds to solving a penalized weighted least square D Detailed Data generating protocol, simulation parameters and experimental protocol

D.1 Data generating protocol

In the simulated data, the data generation protocol is as follows. We consider a K-component functional mixture of Gaussian experts model that relates a scalar response y ∈ R to a univariate functional predictor X(t), t ∈ T defined on a domain T ⊂ R. Given the model parameters (defined in the next paragraph)

, t ∈ T , we first construct the functional predictors X i (•) for i ∈ [n] using the representation defined in (5), with detailed parameterization in (54). Then, for each i ∈ [n], conditional on the functional predictor X i (•), a hidden categorical random variable

is given by (3). Finally, conditional on Z i = z i and X i (•), the scalar response Y i is obtained by sampling from the Gaussian distribution with mean β z i ,0 + T X i (t)β z i (t)dt and variance σ 2 z i . The value z i is then the true cluster label of the predictor X i (•). This hierarchical generative process can be summarized as follows:

The true generated data, i.e., X i (•), Y i , Z i n i=1 , are used for evaluating the prediction and clustering performance.

Finally, to mimic real data, in our simulations, the generated predictors X i (•) are contaminated with measurement errors, that means we will not use X i (•) for analysis, but

with δ i (t) is an independent Gaussian noise with mean zero and constant variance σ 2 δ for all t ∈ T . We considered different noise levels σ 2 δ (see Table 9).

D.2 Simulation parameters and experimental protocol

The parameters that were used in the data generating process (53) are as follows. We consider K = 3 components and a time domain T = [0, 1].

Functional experts and gating parameters: The functional experts parameters are given by

3 ) ⊤ = (5, 5, 5) ⊤ , and the functional gating network parameters are given by

Note that to satisfy the identifiability condition (see [START_REF] Jiang | On the identifiability of mixtures-of-experts[END_REF]), the experts are ordered, for instance (β 1,0 , β 1 , σ 2 1 ) ≺ . . . ≺ (β K,0 , β K , σ 2 K ), and the last gating network parameters, α K,0 , α K (t), are initialized, i.e. fixed to zeros. Here, "≺" is the lexicographical order on R p+2 .

As it can be seen in Figure 11, the expert parameter functions β 1 (t) and β 2 (t) have a flat region in the interval 0.3 ≤ t < 0.7, outside of which they are quadratic, while β 3 (t) and the gating parameter functions α 1 (t), α 2 (t) are all quadratic on the whole domain. By choosing these functional experts and gating parameters, we can later compare the sparsity in the zeroth and third derivatives of the solutions obtained by the proposed models. Functional predictors parameters: In this simulation, the functional predictors X i (•) are constructed using the following formula,

in which b(t) = [b 1 (t), . . . , b 10 (t)] ⊤ is a 10-dimensional B-spline basis x i ∈ R 10 is a coefficient vector defined as x i = W v i , where W ∈ R 10×10 is a matrix of i.i.d. random values from the uniform distribution U(0, 1) and v i ∈ R 10 is a vector of i.i.d. random values from the normal distribution N (0, 10). Here, the matrix W acts as a factor to facilitate the fluctuation of the generated X i (•). Noisy functional predictors: Since in real practical situations we do not usually directly observe the true functional predictors X(•)'s, but only a noisy and discretized version of them, we thus consider several scenarios with different noise and sampling levels. To this end, in the data generating protocol we consider noisy versions U i (t j ) = X i (t j ) + δ i (t j ) of the functional predictors values X i (t j ), where δ i (t j ) ∼ N (0, σ 2 δ ) is a centred Gaussian noise with variance σ 2 δ , for all i ∈ [n], j ∈ [m]. We investigate simulated scenarios S1, . . . , S4 with curve length m and the noise level σ 2 δ of the functional predictors, including two levels of sampling m ∈ {300, 100}, and two levels of noise σ 2 δ ∈ {1, 4}. The resulting four considered scenarios of simulated data are presented in Table 9. For each considered scenario, we generate N = 100 datasets, each dataset contains n pairs of U i (t), Y i , with n ∈ {200, 400, 600, 800, 1000}.