Náoh Êl-Huyen

Anh Tuyen

Nhat, Ban, chi Thanh

Huyen, Thu Hong

Ngoc

Dat, Hung, Nam, Tin, Thanh Thao Vit

Fiancée Thanh

Pham Nhat Thien

K} Operators

⊗ Kronecker

Bayesian Information Criterion ClusErr Clustering Error Corr Correlation EM Expectation-Maximization FDA Functional Data Analysis FME Functional Mixture of Experts FPCA Functional Principle Component Analysis FF-FME Function-on-Function Functional Mixture of Experts FS-FME Function-on-Scalar Functional Mixture of Experts HME Hierarchical mixture-of-experts iFME Interpretable Functional Mixture of Experts i.i.d. Independent and Identically Distributed KL Kullback-Leibler LASSO Least Absolute Shrinkage and Selection Operator LS Least Square ME Mixture-of-experts MISE Mean Integrated Squared Error MLE Maximum Likelihood Estimation MSE Mean Squared Error MM Minorization-Maximization PCA Principle Component Analysis PDF Probability Density Function RI Rand Index RPE Relative Prediction Error SSE Sum of Squared Srrors

Keywords: Mixture of Experts, Functional Data Analysis, Regularized Estimation, LASSO, EM algorithm, MM Algorithm, Distributed Learning, Optimal Transport, Prediction, Clustering Mixture of Experts, Functional Data Analysis, Regularized Estimation, LASSO, EM algorithm, MM Algorithm, Dantzig Selector, Distributed Learning, Optimal Transport Divergence, Clustering, Classification, Regression

Mixture of experts (ME) models are popular in statistics and machine learning and have been studied for high-dimensional vectorial data and that are centralized. However, in many problems, we observe time series data and applying the existing ME models directly to these raw data may limit the performance because they ignore the correlation between variables, an intrinsic nature of functional data, and therefore do not adequately capture the inherent functional structure of the data. In many applications the data may also be not available at a centralized mode, and there is therefore a need to develop adapted distributed strategies allowing for efficient parallel computations. From the modeling point of view, there is indeed a lack of ME models for functional data. From the estimation point of view, there is a lack of approaches to estimate ME models in a distributed fashion. This work addresses both of these problems. This thesis studies mixture of experts (ME) models in functional data and large-scale problems, in a heterogeneous scenario. The main objectives are i) to deal with situations in which we are given functional predictors (e.g., time series) to predict a potentially functional response, and ii) to learn from distributed data. We first propose in this thesis a new family of ME models, called functional mixture of experts (FME), which includes scalar-on-function FME, function-on-scalar FME, and function-on-function FME. We introduce regularized model estimation approaches via appropriate regularizations that encourage sparse and smooth estimates, while being interpretable. We develop efficient EM algorithms to maximize the corresponding observed-data log-likelihood functions and conduct extensive experimental studies to highlight the performance of the proposed models and algorithms. Then, to scale-up the ME estimation to a potentially distributed data, we develop a distributed learning strategy of ME models. It performs standard inference on local machines to obtain local estimators, then transmits them to a central machine where they are aggregated. Based on minimizing a proposed expected transportation divergence, the local estimators are aggregated to obtain a reduced estimator that is consistent with the global one, i.e., the estimator that could be possibly constructed upon the full dataset. Experimental studies demonstrate the performance of our approach. Mixture of experts (ME) models are popular in statistics and machine learning and have been studied for high-dimensional vectorial data and that are centralized. However, in many problems, we observe time series data and applying the existing ME models directly to these raw data may limit the performance because they ignore the correlation between variables, an intrinsic nature of functional data, and therefore do not adequately capture the inherent functional structure of the data. In many applications the data may also be not available at a centralized mode, and there is therefore a need to develop adapted distributed strategies allowing for efficient parallel computations. This thesis studies mixture of experts (ME) models in functional data and largescale problems, in a heterogeneous scenario. The main objectives are i) to deal with situations i in which we are given functional predictors (e.g., time series) to predict a potentially functional response, and ii) to learn from distributed data. We first propose in this thesis a new family of ME models, called functional mixture of experts (FME), which includes scalar-on-function FME, function-on-scalar FME, and function-on-function FME. We introduce regularized model estimation approaches via appropriate regularizations that encourage sparse and smooth estimates, while being interpretable. We develop efficient EM algorithms to maximize the corresponding observed-data log-likelihood functions and conduct extensive experimental studies to highlight the performance of the proposed models and algorithms. Then, to scale-up the ME estimation to a potentially distributed data, we develop a distributed learning strategy of ME models. It performs standard inference on local machines to obtain local estimators, then transmits them to a central machine where they are aggregated. Based on minimizing a proposed expected transportation divergence, the local estimators are aggregated to obtain a reduced estimator that is consistent with the global one, i.e., the estimator that could be possibly constructed upon the full dataset. Experimental studies demonstrate the performance of our approach.

Introduction

Context of study

Nowadays, learning from data is becoming more and more popular. For example, in medicine statistical models have been used to classify patients and detect tumors based on magnetic resonance imaging. In business, sophisticated machine learning algorithms have been used to make business decisions in seconds. In image processing, denoising is no longer done by filters as in classical techniques, but can be automated by machine learning algorithms. And many other applications in various fields such as speech recognition, natural language processing, product recommendation, etc. The above applications all have one thing in common: they are all enabled by learning from data. Moreover, in this big data age, data is cheaper and easier to obtain, which requires more efficient methods to deal with data in large-scale. Therefore, statistical models that are efficient and can be learned in parallel fashion are of great interest. One of the common challenges of modern statistical models is the increasing occurrence of highdimensional data and, more recently, functional data. High-dimensional data have been considered for many models. The common approach is to impose a penalization, usually the ℓ 1 -norm of the parameter vector, when performing model estimation. However, for functional data, classical multivariate models are not appropriate because 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 for dealing with continuously observed data. FDA has been shown to be efficient and flexible in many applications, including regression, classification, and clustering.

For the analysis of heterogeneous data, mixture models [START_REF] Titterington | Statistical Analysis of Finite Mixture Distributions[END_REF][START_REF] Mclachlan | Finite mixture models[END_REF] have become a standard and have shown numerous applications in various fields 1.2 CONTRIBUTIONS AND OUTLINE OF THE THESIS

Contributions and outline of the thesis

The manuscript is organized as follows. Chapter 2 is dedicated to the state of the art. Chapter 3 presents our first contribution to the modeling of mixtures-of-experts for data with functional predictor and continuous scalar response. Chapter 4 is the second contribution which extends the proposed models to the case of classification and to the cases of functional response. There are three new models proposed in this chapter. Chapter 5 presents our last contribution, in which a distributed learning approach for mixture-of-experts models is proposed. Finally, Chapter 6 is for conclusions and discussions of future work. Technical details related to the mathematical developments of our contributions are provided in Appendices A, B and C.

More particularly, in Chapter 2, we attempt to give a comprehensive overview of the current state of the art of models and estimation methods on the research topics of this thesis. First, we present an overview of the mixture-of-experts models, focusing on the modeling of the expert and gating networks, as well as the inference of the models using maximum likelihood estimation. We then introduce some important concepts of functional data analysis, including the mathematical framework, basis expansions, functional principal component analysis, etc. We also provide a review for the models for regression, clustering, as well as classification of functional data. We opt to present in details the estimation methods for the function-on-function and function-on-scalar regression models since they will be used later in our proposed model as extensions.

In Chapter 3, we propose a new family of mixtures-of-experts models, called functional mixtures of experts (FME), to relate a functional predictor to a scalar response. A dedicated EM algorithm for the maximum likelihood parameter estimation is developed for FME model. To deal with potential high-dimensional setting of the proposed FME model, we also consider a Lasso-regularized approach, which consists of a penalized MLE estimation via a hybrid EM-Lasso algorithm. Finally, to obtain a sparse and highly-interpretable regularization of the functional expert and gating parameters, interpretable FME (iFME) model is proposed, which is constructed upon the coefficients of the derivatives of the functional parameters. The resulting model is fitted by regularized MLE via an dedicated EM-iFME algorithm. Extensive experimental studies are then constructed to compare the proposed approaches to the main competitive state of art methods for the subject. The model selection via modified BIC and the implementation details are also discussed. On simulated data, our models outperforms its competitive in both prediction and clustering tasks. The models are then applied to two real-world datasets, on which they also outperforms its competitive on prediction. Specially, the iFME model provides very promising results as it not only produces highly interpretable fits but also high prediction performance. This chapter has leads to the a submission to Journal of Statistics and Computing [START_REF] Chamroukhi | Functional mixturesof-experts[END_REF] (under revision).

Then, in Chapter 4, we extend the FME model proposed in Chapter 3 to classification problem, and to the case of functional responses. In particular, in the FME model for classification, the experts will be modeled by the functional multinomial regression model. In the same spirit with Chapter 3, an interpretable model will also be developed.

To deal with functional responses, we propose in Chapter 4 two ME models: function-on-scalar

Chapter 2

State of the art Contents

Introduction

The goal of this chapter is to provide an overview of the topics of mixture of experts (ME) models and functional data analysis (FDA). The following sections are selective and carefully arranged to introduce the main related concepts of the research topics in this thesis. First, in Section 2.2, the general notions, the modelings and approximation capacity of ME models are presented to provide a first introduction to the research topic. Then, in Section 2.3, we present the estimation of ME models via Expectation-Maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], and review some applications of ME models. In Section 2.4, we present the most basis and important tools of FDA, including the mathematical framework, the basis expansion, data pre-processing methods, and notably the regression models concerning functional data which we will use oftenly through the thesis. In this section, we also review the clustering and classification models for functional data with some illustrative examples. Finally, in preparation for the study in Chapter 5, we present in Section 2.9 some preliminaries on divergence and distance on the space of probability measures.

Mixture of experts models

In this section, we give a general introduction to a class of probabilistic models known as Mixturesof-experts (ME). By the end of this section, the reader will have an overview of this research topic, and will be familiarized with the notations regarding ME models that we will use throughout this thesis.

General notions of ME models

Introduced by [START_REF] Jacobs | Adaptive mixtures of local experts[END_REF], ME is a successful and flexible supervised learning architecture for efficiently representing complex nonlinear relationships in observed pairs of heterogeneous data. The ME model is based on the principle of divide and conquer, so that the response y is obtained from the soft-association of multiple expert responses, each targeting a homogeneous sub-population of the heterogeneous population, given the input covariates (predictors or features). 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 features.

Let x and y be the generic notations for the input and output, respectively, we wish to relate together via ME models. In the most general notion, a ME model can be written as

ME(y|x) = K k=1 Gate k (x) Expert k (y|x),
(2.1)

MIXTURE OF EXPERTS MODELS

in which ME(y|x), the distribution of output y given the input x, is modeled as a mixture distribution with input-dependent mixing proportions Gate k (x) and conditional mixture components Expert k (y|x). ME is therefore a fully conditional mixture model that allows the mixing proportions to be functions of the input. In ME model terminology, Gate k (x) is referred to as gating function, and Expert k (y|x) is referred to as expert function, while K is the number of experts. The general notion (2.1) will be referred many times throughout this thesis as a starting point for developing any specific ME model. Note that the generic variables x and y in (2.1), also commonly referred to as covariate (or predictor) and response, respectively, can be of any data type (such as continuous scalars, vectors, binary, categorical, functions, multivariate functions, etc) and of arbitrary dimension, leading to many different settings and issues for studies and applications regarding ME models. At the end of this chapter, we will try to provide a list of ME models that have been studied so far, as well as point out the open important questions in this research area.

In the notion of probability theory, a ME model can be explained as follows. Let (X, Y ) ∈ X ×Y be a pair random variables which follows some probability model, where X ∈ X is the input, and Y ∈ Y is the output. The X × Y is a subset of some probability space, e.g., R p × R q for some p, q ∈ N * . Then, a ME model can be written as the following decomposition of the output variable Y and the input variable X

P(Y |X) = Z P(Z|X)P(Y |X, Z), (2.2)
in which, Z is a hidden variable that represents the missing data, takes values in the finite set {1, . . . , K}, and indicates which expert was responsible for generating the data point.

The decomposition (2.2) suggests that the data generating process of a ME model is assumed to be as follows: there is a hidden process that specifies the state of Z given X, then the output Y is generated conditional on the values of Z and X. From this decomposition point of view, the functions Gate k (x) and Expert k (y|x) in (2.1) are thus the ways we model the probabilities P(Z|X) and P(Y |X, Z), i.e., how the mixing proportions vary and how the outputs were generated. Because of that, any modeling for the gating function must satisfy the constraint K k=1 Gate k (x) = 1. Thus, depending on the specific parametric forms of the gating function Gate k (x) and the expert functions Expert k (y|x), we have different specifications for ME models. From the application point of view, the main task for us is to select suitable models for the functions Gate k (x) and Expert k (y|x), which depends on the problem and the nature of the data, and to develop efficient algorithms for parameter estimations. From theoretical point of view, researchers focus on studying the approximation capacities, convergence rates of the models, as well as addressing many issues regarding model complexity.

Note that while it is possible to have K different forms for the K expert functions Expert k (y|x), they are usually modeled by a same parametric family, e.g., Gaussian, t-Student, multinomial distribution, etc. Thus, if K is assumed to be known, estimating a ME model consists of estimating the parameters of the gating and expert functions. Figure 2.1 describes the schematic diagram of a K-component ME model.

ME models thus allow one to better capture more complex relationships between y and x in CHAPTER 2 STATE OF THE ART

Expert 1 Expert k Expert K Gate x x x x y g 1 g k g K • • • • • • Figure 2
.1: Visualization of the architecture of a ME model. Here, the gating network and each of K expert receive x as input, then outputs of the experts are weighted by outputs of the gating network, i.e., g 1 , . . . , g K , to produce the final output y.

heterogeneous situations in nonlinear regression, classification, and clustering the data by associating each expert component to a cluster. In situations where the standard mixture of regression models fails to approximate well the conditional probability distribution of data, ME model may be a good substitute because it can capture a wide variety of data distributions thanks to their universal approximation property well established for almost of popular gating and expert functions. The approximation capacity of ME will be briefly reviewed in Subsection 2.2.4.

Remark 2.2.1. We also use the terms "gating network" and "expert network" to refer to the set of all gating functions and the set of all expert functions, respectively. In addition, throughout the thesis, we will frequently say "gate k" and "expert k" instead of gating function k and expert function k, respectively.

So far, x and y represent the generic variables. However, in order to recall the concepts, in the following subsections we will present the modelings of gating and expert functions in the usual vectorial setting, which has been considered in the original proposed papers [START_REF] Jacobs | Adaptive mixtures of local experts[END_REF][START_REF] Xu | An alternative model for mixtures of experts[END_REF][START_REF] Jordan | Hierarchical mixtures of experts and the EM algorithm[END_REF] and many others, although they have been extended to a functional setting that adopts other parametric formulas that we will consider later.

Modeling the gating functions

Let {x i , y i } n i=1 be a sample of n independently and identically distributed (i.i.d.) data, where x i ∈ R p is the input and y i ∈ R q is the output. The most popular models for the function Gate k (x) in (2.1) are softmax gating and Gaussian gating function.

MIXTURE OF EXPERTS MODELS ■ Softmax gating function

With softmax gating modeling, each function Gate k (x) is parameterized by

Gate k (x i ; α) = exp(x ⊤ i α k ) 1 + K-1 k ′ =1 exp(x ⊤ i α k ′ ) , (2.3)
where α = (α ⊤ 1 , . . . , α ⊤ K-1 ) ⊤ , α k ∈ R p for k ∈ [K -1], is the unknown parameter vector to be estimated, here and throughout the thesis the superscript ⊤ denotes the transposition operator. The notation k ∈ [K -1] means k = 1, . . . , K -1, we will use this convention throughout the thesis. Here, from the constraint K k=1 Gate k (x) = 1 in the general notion of ME model, one can see that in softmax gating modeling, the Kth gate is constrained to have the form Gate K (x; α) = 1

1 + K-1 k ′ =1 exp(x ⊤ α k ′ )
, equivalently, the vector α K is constrained to be a null vector. The reason for this parameterization is that, without constraining one of the K vectors α 1 , . . . , α K to a fixed vector, the mixing proportion Gate k (x; α) can admit an arbitrary number of ways of parameterizing, e.g., α = αγ for any γ ∈ R p.K . So, constraining the parameter vector of the Kth gate to be null ensures the identification condition for the estimation of the gating network.

■ Gaussian gating function

An frequently used alternative to the softmax gating function is the Gaussian gating function which has the following parametric form

Gate k (x i ; ζ) = π k ϕ(x i ; ν k , Σ k ) K k ′ π k ′ ϕ(x i ; ν k ′ , Σ k ′ ) , (2.4)
where

ϕ(x; ν, Σ) = det(2πΣ) -1/2 exp - 1 2 (x -ν) ⊤ Σ -1 (x -ν)
is the multivariate Gaussian distribution in R p with mean ν ∈ R p and symmetric positive definite covariance matrix Σ ∈ R p×p . Here, π k > 0 for all k ∈ [K] and K k=1 π k = 1, and ζ = (π 1 , ν ⊤ 1 , vech ⊤ Σ 1 , . . . , π K , ν ⊤ K , vech ⊤ Σ K ) ⊤ is the unknown parameter vector of the gating network to be estimated, where vech(•) is the operator that extracts the unique elements of a symmetric matrix [START_REF] Henderson | Vec and vech operators for matrices, with some uses in jacobians and multivariate statistics[END_REF]). Here, one can easily verify the constraint K k=1 Gate k (x i ; ζ) = 1 for all x i .

The modeling (2.4) was first considered by [START_REF] Xu | An alternative model for mixtures of experts[END_REF], since then there have been many using examples in the literature, e.g., [START_REF] Mak | Estimation of elliptical basis function parameters by the em algorithm with application to speaker verification[END_REF] with application to speaker verification, aki [START_REF] Aki Sato | On-line em algorithm for the normalized gaussian network[END_REF] with application to robot dynamics problems using their proposed online EM algorithm, [START_REF] Lima | Hybridizing mixtures of experts with support vector machines: Investigation into nonlinear dynamic systems identification[END_REF] where Gaussian gating combined with support vector machine (SVM) experts was applied to problems of nonlinear dynamic system identification, etc.

The two softmax and Gaussian gating modelings are shown to be equivalent under some restrictions [START_REF] Nguyen | Practical and theoretical aspects of mixture-of-experts modeling: An overview[END_REF]. In principle, any modeling of the gating network that satisfies the conditions Gate k (x) ⩾ 0 for all x ∈ X , k ∈ [K], and K k=1 Gate k (x) = 1 for all x ∈ X can be used as an alternative to (2.3) and (2.4). However, only a few were considered, including Xu et al. (1994a) where the gating network was modeled using the exponential family (the Gaussian gating network is an instance), and [START_REF] Perthame | Non linear robust regression in high dimension[END_REF] where Student-t gating function was used in combination with Student-t expert functions to model the air quality on the subways in Paris.

Modeling the expert functions

The modeling of expert functions depends on the problem at hand. For regression problems with continuous responses, the Gaussian expert function is often used, while for classification problems the multinomial expert function is the appropriate choice.

■ Gaussian expert function

When the response y is continuous, Expert k (y|x) is usually modeled by the conditional density function of the Gaussian distribution

Expert k (y i |x i ; θ k ) = ϕ(y i ; x ⊤ i β k , σ 2 k ) = (2πσ 2 k ) -1/2 exp - (y i -x ⊤ i β k ) 2 2σ 2 k , (2.5)
where θ k = (β ⊤ k , σ 2 k ) ⊤ is the unknown parameter vector of the expert k to be estimated. This modeling is based on the assumption that the error terms ε i in the linear relationship y i = x ⊤ i β k +ε i have normal distribution. Alternatively, for the data containing a group or groups of observations with asymmetric behavior, heavy tails or atypical observations, other assumptions can be put on the error terms, resulting in skew-normal experts, robust t experts, and skew t experts [START_REF] Chamroukhi | Non-normal mixtures of experts[END_REF][START_REF] Chamroukhi | Skew t mixture of experts[END_REF].

ME models with experts modeled by (2.5) have been widely investigated in many contexts and applied to many regression problems and clustering analysis. We will discuss some notable studies later in this section. For now, we proceed to introduce another important modeling for the experts. When the response y is categorical, i.e., in classification problems, the appropriate choice for the Expert k (y|x) in (2.1) is the conditional density function of the categorical distribution1 given by Expert(y i |x i ;

β k ) = P (y i |x i ; β k ) = G g=1 exp(x ⊤ i β kg ) 1 + G-1 g ′ =1 exp(x ⊤ i β kg ′ ) y ig , (2.6)
where β k = (β ⊤ k1 , . . . , β ⊤ k,G-1 ) ⊤ is the unknown parameter to be estimated, and y ig = I (y i =g) , g ∈ [G]. Here, I denotes the indicator function. Note that, similar to the parameterization of the softmax gating function, the parameter vector responsible for the class G, i.e., β kG , is constrained to be a null vector, this allows the estimation of β k to be identified.

Remark 2.2.2. The name "multinomial expert" comes from the fact that in machine learning literature, it is common to speak of a "multinomial distribution", although a "categorical distribution" would be more precise. In addition, estimating parameter of the conditional density in (2.6) usually appears in regression analysis literature under the name "multinomial logistic regression", rather than "categorical logistic regression". Moreover, it is also known by a variety of other names such as softmax regression, multiclass logistic regression, maximum entropy classifier, etc.

In addition to the above mentions, other interesting models for the experts have been considered in the literature include SVM experts [START_REF] Cao | Support vector machines experts for time series forecasting 3 lijuan cao[END_REF][START_REF] Lima | Hybridizing mixtures of experts with support vector machines: Investigation into nonlinear dynamic systems identification[END_REF], Gaussian processes experts [START_REF] Tresp | Mixtures of gaussian processes[END_REF][START_REF] Jeon | Bayesian mixture of gaussian processes for data association problem[END_REF][START_REF] Yang | An efficient em approach to parameter learning of the mixture of gaussian processes[END_REF]. Notably, [START_REF] Jordan | Hierarchical mixtures of experts and the EM algorithm[END_REF] models each expert by another ME model, leading to the so-called hierarchical mixture of experts (HME) models, an important architecture with many applications, as described below.

■ Hierarchical mixture of experts

In HME models, each expert Expert k (y|x) in (2.1) broadly studied and has found many interesting applications in various fields. In [START_REF] Jordan | Hierarchical mixtures of experts and the EM algorithm[END_REF], the authors illustrated the advantages of HME compared to other approaches such as back propagation network, classification and regression tree (CART), multivariate adaptive regression splines (MARS). In particular, a four-level HME model was used to solve a nonlinear system identification problem, it yielded a relative error better than CART, MARS, and with a running time much shorter than that of the back propagation network. In [START_REF] Ramamurti | Advances in using hierarchical mixture of experts for signal classification[END_REF], the authors presented the use of HME models with multinomial experts for signal classification. In particular, both softmax and Gaussian gating functions were investigated, gathering with multinomial experts, to identify the types of glass in the Glass dataset.

Remark 2.2.3. The softmax-gated mixture of Gaussian regressions and the softmaxgated mixture of multinomial regressions are the most popular choices for many real-world applications. They and their extensions to functional data will be considered in the thesis.

Approximation capacity and applications of ME models

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]. In this section, we review the theoretical results on the approximation capacity of the ME models. Firstly,

ESTIMATION OF ME MODELS

we have an obvious question that: given a continuous function over a compact subset of some Euclidean space, whether we can approximate that function with a ME models? The answer is yes. Moreover, the experts in this case are quiet simple, i.e., linear experts. In particular, consider the ME model in (2.1) in which the gating function Gate k (x) is the soft-max function and the experts Expert k (y|x) are the linear experts. Then the mean function of this ME model can be written as a function of x as

E[y|x] = K k=1 exp(α k0 + x ⊤ α k ) K k ′ exp(α k ′ 0 + x ⊤ α k ′ ) (β k0 + x ⊤ β k ) := h(x).
(2.8)

Let C(X ) be the class of continuous functions and M(X ) = {h(x) : h has form 2.8} be the class of mean functions of obtained from the mixture of linear experts described above. Then provided X ⊂ R p is compact, the class M(X ) is dense within the class C(X ) [START_REF] Nguyen | A universal approximation theorem for mixture-of-experts models[END_REF].

Some successful applications of ME architecture may be mentioned here include ME for time series prediction [START_REF] Zeevi | Time series prediction using mixtures of experts[END_REF][START_REF] Yümlü | Financial time series prediction using mixture of experts[END_REF], segmentation [START_REF] Chamroukhi | Joint segmentation of multivariate time series with hidden process regression for human activity recognition[END_REF][START_REF] Chamroukhi | Time series modeling by a regression approach based on a latent process[END_REF], ME for social network data [START_REF] Gormley | A mixture of experts latent position cluster model for social network data[END_REF], for classification of gender and pose of human faces [START_REF] Gutta | Mixture of experts for classification of gender, ethnic origin, and pose of human faces[END_REF], among many others. For an overview of practical and theoretical aspects of ME modeling, reader is referred to [START_REF] Nguyen | Practical and theoretical aspects of mixture-of-experts modeling: An overview[END_REF]; [START_REF] Yuksel | Twenty years of mixture of experts[END_REF].

Estimation of ME models

So far, we have only introduced the ME models by describing their concepts, the modelings for gating and expert functions, and occasionally some applications, but we have not mentioned the approaches for model inference. This section is devoted to presenting a general overview of the estimation for ME models. At the end of this section we provide a list of notable applications of ME models that have been published in the literature.

In this thesis, the ME models will be considered frequently are the softmax gating Gaussian experts, and the softmax gating multinomial experts (and their variants). By softmax gating Gaussian experts we refer the ME model with the gating function given by (2.3) and the expert function given by (2.5). Whereas, by softmax gating multinomial experts, we refer the ME model with the gating function given by (2.3) and the expert function given by (2.6). For simplicity, from now on we will refer to them respectively as ME model with Gaussian regression experts and ME model with multinomial regression experts. In the subsequent parts, we will individually present the maximum likelihood estimation (MLE) of these two models using EM algorithm.

CHAPTER 2 STATE OF THE ART

MLE for ME model with Gaussian regression experts

Let D = {(x i , y i )} n i=1 be an observed sample of n i.i.d. data pairs (x i , y i ), where x i ∈ R p is the input and y i ∈ R is the output. Suppose that the relationship of x i and y i follows the ME model

f (y i |x i ; θ) = K k=1 π k (x i ; α)ϕ(y i ; x ⊤ i β k , σ 2 k ),
where π k (x i ; α) and ϕ(y i ; x ⊤ i β k , σ 2 k ) are given as in (2.3) and (2.5), respectively. The objective here is to estimate the unknown parameter θ = (α, β 1 , . . . , β k ) via MLE. This can be done by maximizing the observed-data log-likelihood function defined as

log L(θ) = n i=1 log K k=1 π k (x i ; α)φ(y i ; x ⊤ i β k , σ 2 k ) (2.9)
using the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. The EM algorithm starts with an initial solution θ (0) , then alternates the two following steps until convergence, i.e., there is no significant change in the values of the function log L(θ).

■ E-step. For i = 1, . . . , n, calculate the conditional probability memberships τ (s) ik that the observed pair (x i , y i ) originates from the kth expert, given the observed data and the current estimate θ (s) , defined by

τ (s) ik = π k (x i ; α (s) )φ(y i ; x ⊤ i β (s) k , σ 2 k (s) ) K k=1 π k (x i ; α (s) )φ(y i ; x ⊤ i β (s) k , σ 2 k (s) ) • ■ M-step. Given the memberships τ (s)
ik , update the value of the parameter θ by maximizing w.r.t. θ the so-called Q-function defined by

Q(θ; θ (s) ) = n n=1 K k=1 τ (s) ik log π k (x i ; α)φ(y i ; x ⊤ i β k , σ 2 k ) .
(2.10)

The maximization can be performed separately w.r.t. the gating network parameter α and the expert network parameters β k , σ k as follows.

Updating the gating network parameter α consists of maximizing w.r.t. α the function Q(α; θ (s) ) that is the part of (2.10) involved α. This is a weighted multinomial logistic regression problem with the weights are the conditional probability memberships τ (s)

ik . Since there is no closedform solution for such problem, one common approach is to use the iteratively reweighted least squares (IRLS) algorithm [START_REF] Jordan | Hierarchical mixtures of experts and the EM algorithm[END_REF][START_REF] Ng | Using the em algorithm to train neural networks: misconceptions and a new algorithm for multiclass classification[END_REF]) that starts from an 2.3 ESTIMATION OF ME MODELS initial estimate α (0) , then updates α iteratively via the following updating formula (2.11) at each iteration t, where H(α; θ (s) ) and g(α; θ (s) ) are respectively the Hessian matrix and the gradient vector of Q(α; θ (s) ), given in detailed, for example, in the Appendix B.1. We keep running this inner-loop until there is no significant change in Q(α; θ (s) ), the obtained solution is then an estimate α (s) for the next EM iteration.

α (t+1) = α (t) -H(α; θ (s) ) -1 α=α (t) g(α; θ (s) )| α=α (t)
Updating the expert network parameters β k and σ k consists of solving K independent weighted regression problems with the weights are the conditional probability memberships τ (s)

ik . This can be done easily with the following closed-form updating formulas:

β (s+1) k = n i=1 τ (s) ik x i x ⊤ i -1 n i=1 τ (s) ik y i x i σ 2 k (s+1) = 1 n i=1 τ (s) ik n i=1 τ (s) ik (y i -x ⊤ i β (s+1) k ) 2 .
The EM algorithm increases the value of the log likelihood function log L(θ) at each iteration, therefore it is guaranteed to converge to at least a local minimum after a finite number of iterations.

MLE for mixture of multinomial regression experts model

ME model with multinomial regression experts is suitable for multiclass classification problems. Let D = {(x i , y i )} n i=1 be an observed sample of n i.i.d. data pairs (x i , y i ), where x i ∈ R p is the input, and y i = (y i1 , . . . , y iG ) ∈ R G is the categorical output. Suppose that the relationship of x i and y i follows the ME model

f (y i |x i ; θ) = K k=1 π k (x i ; α)P (y i |x i ; β k ),
where π k (x i ; α) and P (y i |x i ; β k ) are given as in (2.3) and (2.6), respectively. The objective here is to estimate the unknown parameter θ = (α, β 1 , . . . , β k ) via MLE. Similarly to the ME model with Gaussian regression experts, this can be done by maximizing the following observed-data log-likelihood function

log L(θ) = n i=1 log K k=1 π k (x i ; α)P (y i |x i ; β k ) (2.
= π k (x i ; α (s) )P (y i |x i ; β (s) k ) K k=1 π k (x i ; α (s) )P (y i |x i ; β (s) k ) •
In the M-step of the EM algorithm, the update for the gating network parameter α is identified with the update for the gating network of the ME model with Gaussian regression experts. This is equivalent to solving a weighted multinomial logistic regression problem, where the weights are the conditional probability memberships τ (s) ik . However, for each multinomial expert, its parameter vector can no longer be updated in a closed-form as that of the Gaussian experts. It is instead updated using Newton-Raphson method [START_REF] Fletcher | Practical methods of optimization[END_REF], which performs an update formula analogous to the update formula of IRLS algorithm in (2.11), but the Hessian matrix and gradient vector are now weighted by the posterior probabilities τ (s) ik of the current EM iteration. The calculations for these Hessian matrix and gradient vector are provided, for example, in the Appendix B.1. In particular, during the Newton-Raphson loop, the weights τ (s) ik are fixed, and we keep running this loop until its convergence. The point obtained at convergence is then the update for β k for the next EM iteration.

Regularized MLE for ME model

Alternatively to the MLE, to encourage sparsity in the parameter vectors α and β k 's in highdimensional settings, one may estimate the regularized MLE via maximizing the ℓ 1 -regularized observed-data log-likelihood L(θ) = log L(θ) -Pen λ,χ (θ), (2.13) where log L(θ) is the observed-data log-likelihood of θ defined in (2.9) or in (2.12), and Pen λ,χ (θ) is a LASSO regularization term defined by

Pen λ,χ (θ) = λ K k=1 ∥β k ∥ 1 + χ K-1 k=1 ∥α k ∥ 1 ,
where λ, χ ∈ R are the tuning parameters and ∥ • ∥ 1 is the ℓ 1 -norm operator. The solution to the maximization of (2.13) cannot be obtained in closed form. However, the EM algorithm can be adapted to iteratively maximize (2.13), e.g., see [START_REF] Chamroukhi | Regularized Maximum Likelihood Estimation and Feature Selection in Mixtures-of-Experts Models[END_REF].

The MLE and regularized MLE presented above will be revisited frequently throughout this thesis. In the next chapters, our proposed ME models will oftenly refer to these frameworks for model estimation, the corresponding EM algorithms will also be presented in detailed particularly for each model.

FUNCTIONAL DATA ANALYSIS

Functional data analysis

Functional data analysis (FDA) is a branch of statistics that analyzes and provides information about data of types curves, surfaces or anything else varying over a continuum. FDA has roots going back to work by Grenander and Karhunen in the 1940s and 1950s (Grenander, 1950;[START_REF] Müller | Functional data analysis and random objects[END_REF]. Since then, FDA has been rapidly developed to be a rich subfield of statistics, especially in the last two decades. The term Functional Data Analysis was first coined by [START_REF] Ramsay | Functional Data Analysis[END_REF].

One of the distinctive features of FDA is its ability to incorporate the correlation within the functional data as well as their derivatives. The key techniques of FDA include data smoothing, dimensionality reduction, functional data clustering, classification, functional linear modeling and forecasting methods. In the subsequent sections, we will review the regression models and the methods of clustering as well as classification for functional data. An excellent review paper on the applications of FDA can be found in [START_REF] Ullah | Applications of functional data analysis: A systematic review[END_REF].

This section aims first to outline some very basic and important concepts of FDA, e.g., the mathematical framework, the basis expansions, and some tools for pre-processing functional data (such as smoothing via penalizing derivatives, curve alignment, etc). Secondly, this section aims to review the current state of the art works on FDA, including regression, clustering, and classification for functional data, as well as some related theoretical results. The popular linear models concern functional data will be presented in details in order to prepare for the ME models proposed in the next chapters. For convenience, in this thesis, we use the notation, e.g., X(t) (or sometimes X(•) if there is no ambiguity) to refer to a functional object (more clearly, a function of t), not to a particular value of X at any t. If we want to refer to a particular value, we will make that clear.

Mathematical framework

By mathematical framework, we aim at presenting here the general view of the elements (i.e., the functions) that we will work with in FDA, and the basic properties of the space they live in. Firstly, in FDA, each sample element of functional data is considered to be a random function. The random functions can be viewed as random elements taking values in some Hilbert space. The common one is L 2 (T ), the space of all square integrable functions defined on some interval T ⊂ R, and we will mainly working on it.

A function X(t) is said to be square integrable on

T , i.e., X ∈ L 2 (T ), if T [X(t)] 2 dt < ∞.
As a Hilbert space, L 2 (T ) enjoys all properties and theoretical results those have been developed for Hilbert spaces. Let X, Y be two functions in L 2 (T ), their inner product is defined by

⟨X, Y ⟩ = T X(t)Y (t)dt.
CHAPTER 2 STATE OF THE ART Analogous to the Euclidean space, each function X in L 2 (T ) has an associated norm, define by ∥X∥ = ⟨X, X⟩, and the distance between two functions, e.g., X and Y , is also defined as the norm of their difference, i.e., .14) In addition to the Hilbertian point of view, it is common to view random functions from the applied perspective as stochastic processes, which is more suitable in our context of study in this thesis. In particular, each random function X(t) can be viewed as a collection of random variables {X(t)} t∈T , with T ⊂ R. From this perspective, we can also see that functional data are of infinite dimensional intrinsically.

∥X -Y ∥ = T (X(t) -Y (t)) 2 dt 1/2 . ( 2
Mean and covariance functions. Let X be a random function. Similar to the multivariate case, we have the notions of mean and covariance for random functions. The mean of X is also a function of t, usually denoted by µ, namely, µ(t)

:= E[X(t)]. The covariance of X is defined by c(t, u) = E [(X(t) -µ(t))(X(u) -µ(u))] , for t, u ∈ T .
Let X 1 , . . . , X n are realizations of X. The sample mean and sample covariance functions, respectively, are given by

μ(t) = 1 n N i=1 X i (t), ĉ(t, u) = 1 n -1 n i=1 (X i (t) -μ(t)) (X i (u) -μ(u)) .
The sample mean and sample covariance functions can be viewed as the estimators of µ(t) and c(t, u) defined above. It has been shown that μ and ĉ converge to µ and c in L 2 -norm with convergences rate of O(n -1/2 ) [START_REF] Jacques | Functional data clustering: A survey[END_REF][START_REF] Deville | Méthodes statistiques et numériques de l'analyse harmonique[END_REF]. These functions play an important role in summary statistics for functional data since they provide a view of the population, similar to the sample mean and covariance in multivariate data analysis. We will revisit these functions oftenly throughout the thesis.

Basis expansions

The very first step in working with functional data is to express them via a basis expansion as

X i (t) ≈ K k=1 x ik ω k (t), i ∈ [n],
( 2.15) in which, x ik are the coefficients, and ω k (t) are some collection of basis functions. There are some popular collections that are oftenly used in FDA, including B-splines, wavelet, Fourier, Gaussian 2.4 FUNCTIONAL DATA ANALYSIS and sigmoidal. Figure 2.3 shows some examples of the basis functions collections. Since we will meet basis expansions many times later, we would mention here some useful conventions on using basis expansions. Firstly, throughout the thesis the expansions as in (2.15) are usually written with an equal sign, i.e., "=". In practice, functional data is often available in form of discretized values. However, through the representation with basis functions, the functional data will become available on the whole domain. Therefore, when the time points differ between the functions, using the basis expansions we can put the functions into a common domain. Analogous to the principal component analysis (PCA) in the classical multivariate analysis, FPCA is an important tool in FDA that helps explaining the variance in the population, via the principal components (PCs), extracting the features, visualizing the data, and is a widely used technique for dimensionality reduction. The concept of FPCA can be described from the application point of view as follows. Let X be a random function, and {X i } n i=1 be an i.i.d. sample of X. We wish to find a decomposition for CHAPTER 2 STATE OF THE ART the X i as (2.16) such that the coefficients ξ ik are i.i.d., uncorrelated, and have mean zero; the functions µ(t) and υ k (t) are independent of i, υ k (t) are orthonormal in L 2 -space, and importantly X i can be well approximated by only a small number of the functions υ k (t). One can see, there are many benefits can be obtained from this decomposition, especially, the dimensionality reduction aspect. Nowadays, such decomposition in (2.16) is well-known as the Karhunen-Loève expansion of X, and was in fact conceived from 1946 in [START_REF] Karhunen | Zur Spektraltheorie stochastischer Prozesse[END_REF]; [START_REF] Loève | Fonctions aléatoires à décomposition orthogonale exponentielle[END_REF]. Since then, a framework, termed as FPCA, was developed to compute such coefficients and functions as in (2.16) and use them for further data explorations.

X i (t) = µ(t) + ∞ k=1 ξ ik υ k (t),
Let c(t, u) be the covariance function of X. Consider the following eigenvalue equation

T c(t, u)υ(u)du = λυ(t).
Denote by λ 1 , λ 2 , . . . the non-increasing ordered eigenvalues of the above equation, and υ 1 , υ 2 , . . . their corresponding eigenfunctions. Let µ(t) be the mean function of X and ξ ik be computed by

ξ ik = T (X i (t) -µ(t))υ k (t)dt.
(2.17)

Then µ(t), ξ ik and υ k (t) fulfill the expansion (2.16). In other words, each X i can be approximated using the mean function, and a linear combination of the eigenfunctions of c(t, u) with coefficients are computed as in (2.17). The eigenfunctions υ k (t) are called the functional principal components (FPCs). The coefficients ξ ik are called the scores of X i w.r.t. to the FPCs υ k (t).

The important properties of the scores ξ ik are that they are independent across i, uncorrelated across k, i.e.,

E[ξ ik ] = 0, Cov(ξ ik , ξ iℓ ) = 0 if k ̸ = ℓ, moreover, we have Var[ξ ik ] = λ k , and E ∥X -µ∥ 2 = ∞ k=1 λ k .
The latter equation shows that the total variance of X can be decomposed into the sum of the eigenvalues, i.e., variance of its scores. Because λ k 's are arranged in non-increasing order, we can quantify the explanation of the corresponding principal components. We will use FPCA in Chapter 3 as an alternative technique to expand the functions in our models. Applications of FPCA. Several works have used FPCA as the starting point to represent the functional data in their models. By using the scores of the functions corresponding to first few principal components, the construction of the parametric models becomes more parsimonious. In Yao et al. (2010), the authors uses FPCA and proposed a mixture regression model for scalar-onfunction regression.

Scalar-on-function regression

Functional regression models come to deal with situations where one or both of the input and output variables are functions. The coefficient parameters must therefore be appropriately defined for such situations. This section and the next two ones are dedicated to present the three most considered models involving functional data, including: scalar-on- is a functional predictor with t ∈ T ⊂ R. The scalar-on-function regression model [START_REF] Müller | Generalized functional linear models[END_REF][START_REF] Ramsay | Functional Data Analysis[END_REF] is given by

y i = β 0 + T X i (t)β(t)dt + ε i , i ∈ [n], (2.18)
where β 0 is the unknown intercept and β(t) is the unknown functional parameter to be estimated,

ε i ∼ N (0, σ 2 ) are independent Gaussian errors.
The popular approach is to expand the functional parameter β(t) using a collection of basis functions ω 1 (t), . . . , ω K (t), i.e.,

β(t) = K k=1 b k ω(t), (2.19)
whee K is generally chosen to be smaller the number of time points at which we observe X i (t), but sufficiently large to ensure to capture the vary of the functions. Note that, the equation is only approximation, but we use the equality sign "=" to be able to manipulate such expansion without 2.6 FUNCTION-ON-SCALAR REGRESSION adding additional error terms. Using the expansion (2.19), the model can be rewritten as

y i = β o + K k=1 b k T X i (t)ω k (t) + ε i = β 0 + K k=1 b k x ik + ε i = β 0 + x i b + ε i , i ∈ [n],
where

x i = (x [ i1], . . . , x iK ) ⊤ with x ik = cT X i (t)ω k (t)
dt is now the new design vector, and b = (b 1 , . . . , b K ) ⊤ is the unknown parameter to be estimated. The intercept β 0 and the parameter vector b are then can be estimated by (

β 0 , b ⊤ ) ⊤ = (X ⊤ X) -1 X ⊤ Y , where X and Y denote X =       1 x 11 • • • x 1K 1 x 21 • • • x 2K . . . . . . . . . . . . 1 x n1 • • • x nK       , Y =       y 1 y 2 . . . y n       .
Above is the easy-to-use approach for fitting parameter of scalar-on-function regression model. But it is not the only one. Other approaches include functional principal component regression method [START_REF] Reiss | Functional principal component regression and functional partial least squares[END_REF], hybrid method that consists of a FPCA to pre-process predictor functions and penalized splines to model the coefficient function (Goldsmith et al., 2011a). In [START_REF] Yao | Functional quadratic regression[END_REF], the authors proposed a functional quadratic regression model whose model parameter is estimated using the FPCA approach. Scalar-on-function regression has shown various applications in many fields, including medicine, bioinformatics, etc. In the next section we will present the two less-common but still important models in the family of regression models for functional data, they are function-on-scalar, and function-on-function regression models.

Function-on-scalar regression

Many problems today involve modeling functional responses via scalar explanatory variables, especially in the biological sciences and engineering design. In this section, we present the formulation and estimation methods for the so-called function-on-scalar regression model, which links a functional response to scalar covariates via a linear relationship. This regression model will serve as one of the extensions for the expert regression model that we will consider in Chapter 4.

Model formulation

Y i (t) = β 0 (t) + x i1 β 1 (t) + . . . + x ip β p (t) + ε i (t), i ∈ [n],
( 2.20) where β 0 (t) is the unknown functional intercept and β 1 (t), . . . , β p (t), t ∈ T , are the unknown functional parameters, often referred to as effect functions, to be estimated. Here and throughout the thesis, the notations x ij 's are implicitly understood as the components of the vector x i , i.e., x i = (x i1 , . . . , x ip ). The error function ε i (t) is commonly assumed to be independent and has a normal distribution for all t ∈ T and for all i ∈ [n]. However, this assumption on the error function could also be replaced by the assumption that ε i (t) is a Gaussian process with mean zero and some covariance function Σ ε (t, s), i.e., in this case we allow the errors to be correlated within individuals, but not between individuals.

The model (2.20) was considered in [START_REF] Faraway | Regression analysis for a functional response[END_REF] and was applied to an ergonomics design problem. In [START_REF] Faraway | Regression analysis for a functional response[END_REF], the parameters were estimated pointwisely by the least squares approach. However, it is a pointwise estimation and is only appropriate when the response functions are smooth. When the responses are noisy, a penalty method, as presented in the next subsection, is preferable.

For convenience, let us stack the variables into matrix form as follows

Y (t) =       Y 1 (t) Y 2 (t) . . . Y n (t)       , X =       x 11 x 12 • • • x 1p x 21 x 22 • • • x 2p . . . . . . . . . . . . x n1 x n2 • • • x np       , β(t) =       β 1 (t) β 2 (t) . . . β p (t)       , ε(t) =       ε 1 (t) ε 2 (t) . . . ε n (t)       • (2.21)
Then, the equations (2.20) can be written in the familiar form as

Y (t) = Xβ(t) + ε(t). (2.22)
Note that, in (2.21) and throughout the thesis, in terms of data organization, each curve is aligned as a row of observed values. For example, in practice if the responses are observed in a same grid t 1 , . . . , t J then such Y (t) will be available as a n × J matrix. Like that, the matrix multiplication in the equation (2.22) is ensured to be correct in practice.

2.6 FUNCTION-ON-SCALAR REGRESSION

Penalized least squares estimation

Provided that X has full rank as in the usual regression situation, the least squares estimator of β(t) can be obtained pointwisely by

β LS (t) = (X ⊤ X) -1 X ⊤ Y (t), ∀t ∈ T .
However, as mentioned earlier, such pointwise least squares estimation may be not suitable for functional data, especially, for noisy data. In this case, the penalized least squares estimation of β(t) (e.g., (Kokoszka and Reimherr, 2017, Chapter 5)) defined as follow is more preferred:

β(t) = min β(t)   T ∥Y (t) -Xβ(t)∥ 2 dt + p j=1 λ j T [(Dβ j )(t)] 2 dt   , (2.23)
where λ j are the penalty constants and D is some differential operator, for example, the second derivative. Here, ∥f(t

)∥ 2 = n i=1 f 2 (t) for all f(t) = [f 1 (t), . . . , f n (t)] ⊤ .
The following part presents an approach to solve the penalized least squares estimator β(t) in (2.23).

Let ω(t) = [ω 1 (t), . . . , ω K (t)] ⊤ be a vector of basis functions. Then for each i ∈ [n] and j ∈ [p], one can approximate the functions Y i (•) and β j (•) respectively by

Y i (t) = K k=1 r ik ω k (t) =: r ⊤ i ω(t), and β j (t) = K k=1 b jk ω k (t) =: b ⊤ j ω(t), (2.24)
where r i , b j ∈ R K are the coefficient vectors. Here, K should ensure the tradeoff between smoothness of the functions and complexity of the estimation problem. By using these basis expansions and denoting R = [r 1 , . . . ,

r n ] ⊤ ∈ R n×K , B = [b 1 , . . . , b p ] ⊤ ∈ R p×K , we can rewrite the equation (2.22) as R ω(t) = XB ω(t) + ε(t).
Similarly, the function (Dβ j )(t) can also be expressed in terms of B and ω(t) as

(Dβ j )(t) = K k=1 b jk (Dω)(t) =: b ⊤ j (Dω)(t),
here (Dω)(t) obviously denotes the vector obtained by operating D on the components of ω(t), i.e., (Dω)(t) = [(Dω 1 )(t), . . . , (Dω K )(t)] ⊤ . Then, the integral in the penalty term can be approximated by

T (Dβ j )(t) 2 dt = T b ⊤ j (Dω)(t).[(Dω)(t)] ⊤ b j dt =: b ⊤ j Db j ,
where D denotes the remaining cross-product matrix after factorizing b j , that is

D = T (Dω k )(t)(Dω ℓ )(t)dt 1⩽k,ℓ⩽K ∈ R K×K . (2.25)
It is worth noting that D is a non-negative definite matrix. Hence, the objective function in (2.23) can now be rewritten by

T ∥R ω(t) -XB ω(t)∥ 2 dt + p j=1 λ j b ⊤ j Db j .
(2.26)

The coefficient matrix B (hence b j 's) minimizing (2.26) can be obtained by solving a first order differential equation as for ridge regression in multivariate analysis. Particularly, the penalized least squares estimator for the coefficient matrix B can be calculated in a condensed form using the vectorization operator and the Kronecker product operator as follow:

vec( B ⊤ ) = U ⊤ U + Λ ⊗ D -1 U ⊤ vec(I 1/2 ω R ⊤ ), U = X ⊗ I 1/2 ω , (2.27)
in which, Λ = diag(λ 1 , . . . , λ p ) and

I 1/2 ω
is the square root of the non-negative definite matrix I ω defined by

I ω = T ω k (t)ω ℓ (t)dt 1⩽k,ℓ⩽K ∈ R K×K .
Note that if the chosen basis is orthonormal then I ω is reduced to an identity matrix. The details of the derivation of (2.27) can be found in Appendix A.1.1.

Remark 2.6.1. The solution for B given in (2.27) will serve in Section 4.4 as the update formula for the experts in our proposed function-on-scalar ME model.

Function-on-function regression

This type of regression problem appears in the situation where both the responses and predictors are functions. In this case, the relation between functional response Y i (t), t ∈ T , and functional 2.7 FUNCTION-ON-FUNCTION REGRESSION predictor X i (u), u ∈ U, can be modeled by the following fully functional linear model

Y i (t) = β 0 (t) + U X i (u)β(t, u)du + ε i (t), i ∈ [n],
where β 0 (t), defined on T , is the unknown functional intercept and β(t, u), defined on T × U, is the unknown bivariate regression coefficient function, often referred to as the kernel, to be estimated. Similarly to the function-on-scalar regression case, the error function ε i (t) is usually assumed to be independent for all i ∈ [n], t ∈ T and has Gaussian distribution with mean zero and covariance matrix Σ ε , i.e., ε i (t) ∼ N (0, Σ ε ). From now on we will refer to this model by function-on-function regression model (FFRM).

FFRM model was first considered by [START_REF] Ramsay | Functional Data Analysis[END_REF], in which the functional parameter β(t, u) is estimated by minimizing the integrated sum of squares. [START_REF] Shimokawa | An expansion of functional regression analysis (in japanese)[END_REF] extended it to the case of there are more than one covariate, i.e., functional multiple regression model. [START_REF] Yamanishi | 8. functional data analysis geographically weighted functional multiple regression analysis: A numerical investigation[END_REF] suggested a geographically weighted regression model to explore the functional relationship between the variables. [START_REF] Müller | Functional additive models[END_REF] proposed an adaptive approach to FFRM where the functional parameters are estimated with regularization. In [START_REF] Matsui | Regularized functional regression modeling for functional response and predictors[END_REF], the authors proposed a maximum penalized likelihood approach to estimate β(t, u) since the maximum likelihood approach and least squares approach have been shown to potentially produce unstable estimates for the functional parameter β(t, u). Recently, [START_REF] Beyaztas | On function-on-function regression: Partial least squares approach[END_REF] proposed a partial least squares approach to estimate the parameter of the FFRM model, overcoming the singular matrix problem that occurs when a large number of functional predictors are included in the FFRM.

Model formulation and least squares estimation

Since one can replace the functions X i and Y i by their centered versions, the functional intercept β 0 (t) can be simplified from the model. Therefore, in the remainder we consider the following model .28) with assumption that EX = 0 and consequently EY = 0.

Y i (t) = U X i (u)β(t, u)du + ε i (t), i ∈ [n], ( 2
Let ω ⋆ (t, u) = {ω ⋆ 1 (t, u), . . . , ω ⋆ K (t, u)} be a collection of bivariate basis functions defined on T × U. Then with K sufficiently large, the bivariate functional parameter β(t, u) can be expressed by

β(t, u) = K k=1 b k ω ⋆ k (t, u), (2.29)
CHAPTER 2 STATE OF THE ART with b k are the coefficients. Therefore, the model (2.28) can be rewritten as 2.30) where X ⋆ ik (t) := U X i (u)ω ⋆ k (t, u)du are functions of t which serve as new predictors of the model, and b k are now the unknown coefficients to be estimated. Let us denote

Y i (t) = U X i (u) K k=1 b k ω ⋆ k (t, u) du + ε i (t) = K k=1 b k U X i (u)ω ⋆ k (t, u)du + ε i (t) =: K k=1 b k X ⋆ ik (t) + ε i (t), i ∈ [n], (
X ⋆ (t) =       X ⋆ 11 (t) X ⋆ 12 (t) • • • X ⋆ 1K (t) X ⋆ 21 (t) X ⋆ 22 (t) • • • X ⋆ 2K (t) . . . . . . . . . . . . X ⋆ n1 (t) X ⋆ n2 (t) • • • X ⋆ nK (t)       , b =       b 1 b 2 . . . b K       •
Then the equations (2.30) can be rewritten simply under matrix form as

Y (t) = X ⋆ (t) b + ε(t), (2.31)
where Y (t) and ε(t) are as in (2.21), and the matrix multiplication is understood formally, i.e., each X ⋆ ik (t) is treated as an element of the matrix X ⋆ (t). Such equations (2.22) and (2.31) are the yet simplest and condensed representations for the considered models.

Suppose that the Y i (t) are observed on a finite grid which does not depend on i, and X ⋆ i (t) can also be calculated on the same grid. If we organize all the curves in Y (t) and X * (t) vertically (to make the equation (2.31) holds in practice) then we can estimate b directly with raw data by

b = X ⋆⊤ (t)X ⋆ (t) -1 X ⋆⊤ (t)Y (t).
However, analogously as in functional-on-scalar model, such raw approach is sensitive with noise and potentially produce rough estimate for β(t, u). Therefore, here again, we express Y i (t) and X ⋆ ik (t) using basis expansions and work on their coefficient vectors as done previously for the functional-on-scalar model. Moreover, a penalization on the derivatives of β(t, u) is necessary to produce a smooth estimation for β(t, u).

Remark 2.7.1. The collection of bivariate basis functions ω ⋆ (t, u) are generally constructed using two collections of univariate basis functions. In particular, if {ω k (t), k ⩾ 1} and {ψ ℓ (u), ℓ ⩾ 1} are bases in L 2 (T ) and L 2 (U), then {ω k (t)ψ ℓ (u), k, ℓ ⩾ 1} is a basis in L 2 (T × U). The orthonormal property can also be inherited.

FUNCTION-ON-FUNCTION REGRESSION

Penalized least squares estimation

Similarly to (2.23), the penalized least squares estimation of the parameter β(t, u) is defined as follow 2.32) where Pen λ (β(t, u)) is the penalty term that responsible for the smoothness of β(t, u). It it common to penalize the t and u directions separately with a penalty term taking the following form

β(t) = min β(t) T ∥Y (t) -X ⋆ (t) b∥ 2 dt + Pen λ (β(t, u)) , (
Pen λt,λu (β(t, u)) = λ t [(D t β)(t, u)] 2 dtdu + λ u [(D u β)(t, u)]
Y i (t) = L ℓ=1 r iℓ ψ ℓ (t) =: r ⊤ i ψ(t), X ⋆ ik (t) = L ℓ=1 x ⋆ ikℓ ψ ℓ (t) =: x ⋆⊤ ik ψ(t),
where r i ∈ R L and x ⋆ ik ∈ R L are the vectors of coefficients. In practice, the coefficients x ⋆ ikℓ can be calculated directly using X i (t) and the basis responsible for covariate u that construct ω ⋆ (t, u), see Remark 2.7.1, particularly,

x ⋆ ik = U X i (u)ω 1 (u)du, . . . , U X i (u)ω L (u)du ⊤ ∈ R L . Let us denote R = [r ⊤ 1 , . . . , r ⊤ n ] ⊤ ∈ R n×L , and
G =       x ⋆⊤ 11 x ⋆⊤ 12 • • • x ⋆⊤ 1K x ⋆⊤ 21 x ⋆⊤ 22 • • • x ⋆⊤ 2K . . . . . . . . . . . . x ⋆⊤ n1 x ⋆⊤ n2 • • • x ⋆⊤ nK       ∈ R n×KL .
Then, the model (2.31) can be rewritten in the basis expansion fashion as 2.34) in which, we conventionally treat each vector x ⋆⊤ ik as an element when performing the matrix multiplication. Thus, the function-on-function regression problem (2.28) now has been approximated by the tractable finite dimensional vector-on-vector regression problem (2.34).

R ψ(t) = Gb ψ(t) + ε(t), (
Next, analogously to the penalized least squares estimation for function-on-scalar model, the penalty term in (2.32) can also be expressed in terms of the coefficient vector b and derivatives of the bivariate basis functions. In particular, we have

[(D t β)(t, u)] 2 dtdu = K k,ℓ=1 b k b ℓ (D t ω k )(t, u).(D t ω ℓ )(t, u)dtdu =: b ⊤ D t b, (2.35)
where D t denotes the K × K matrix whose entries calculated by

(D t ω k )(t, u).(D t ω ℓ )(t, u)dtdu, with k, ℓ ∈ [K].
Similarly, the smoothness penalty term corresponding to the u direction can also be expressed by

[(D u β)(t, u)] 2 dtdu = K k,ℓ=1 b k b ℓ (D u ω k )(t, u).(D u ω ℓ )(t, u)dtdu =: b ⊤ D u b, (2.36)
with D u is defined correspondingly. Finally, given the above notations, the penalized least squares estimator of b takes the following formula (2.37) in which V is defined by

b = V ⊤ V + λ t D t + λ u D u -1 V ⊤ vec(I 1/2 ψ R ⊤ ),
V =       I 1/2 ψ x ⋆ 11 I 1/2 ψ x ⋆ 12 • • • I 1/2 ψ x ⋆ 1K I 1/2 ψ x ⋆ 21 I 1/2 ψ x ⋆ 22 • • • I 1/2 ψ x ⋆ 2K . . . . . . . . . . . . I 1/2 ψ x ⋆ n1 I 1/2 ψ x ⋆ n2 • • • I 1/2 ψ x ⋆ nK       ∈ R nL×K ,
where

I 1/2 ψ
is the square root of the non-negative definite matrix I ψ defined by

I ψ = T ψ k (t)ψ ℓ (t)dt 1⩽k,ℓ⩽L ∈ R L×L .
The detailed calculation to establish (2.37) can be found in Appendix A.1.2.

2.7 FUNCTION-ON-FUNCTION REGRESSION

Maximum likelihood estimation

Another way to approach the model estimation is via maximum likelihood method. Firstly, we assume that the error functions ε i (t) can also be represented by linear combination of basis functions in ψ(t), i.e.,

ε i (t) = L ℓ=1 ε iℓ ψ ℓ (t) =: ε ⊤ i ψ(t),
with ε i ∼ N L (0, Σ), for some unknown symmetric matrix Σ ∈ R L×L . Then, the equation (2.34) can be rewritten by

R ψ(t) = Gb ψ(t) + Eψ(t)
,

where E = (ε ⊤ 1 , . . . , ε ⊤ n ) ⊤ ∈ R n×L
is the matrix of error coefficients. Observe that the matrix ψ(t)ψ ⊤ (t) is non-singular, therefore the above equation implies

R = Gb + E, (2.38)
which takes a familiar form as in ordinary multivariate regression if the matrix multiplication is conventionally understood as in (2.34). Alternatively, if we reframe the matrices R and E into column vectors, equation (2.38) can also be rewritten as vec(R ⊤ ) = Hb + vec(E ⊤ ), (2.39) where the matrix H, differs from G, is formed by the vectors x ⋆ ik ∈ R L as follow

H =       x ⋆ 11 x ⋆ 12 • • • x ⋆ 1K x ⋆ 21 x ⋆ 22 • • • x ⋆ 2K . . . . . . . . . . . . x ⋆ n1 x ⋆ n2 • • • x ⋆ nK       ∈ R nL×K . (2.40)
Note, however, that (2.39) is not the usual linear equation met in linear regression with continuous response, because the noise here may not be i.i.d., and more importantly, the observations are not independent of each other since the matrices are correlated within each L-rows block due to the correlation between the components of r i and x ⋆ ik . From the representation (2.38), we see that the model estimation can be viewed as maximum likelihood estimation for the following probability density function

f (r i |x i ; θ) = det(2πΣ) 1/2 exp 1 2 (r i -Z i b) ⊤ Σ -1 (r i -Z i b) , (2.41)
where .43) where the matrix H is given in (2.40).

Z i := [x ⋆ i1 , . . . , x ⋆ iK ] ∈ R L×K ,
= H ⊤ H -1 H ⊤ vec(R ⊤ ), (2.42) Σ = 1 n R -G b ⊤ R -G b . ( 2
Remark 2.7.2. The solution for b and Σ presented here will serve in Section 4.4 as the update formula for the experts in our proposed function-on-scalar ME model.

Clustering and classification of functional data

Clustering and classification are important topics of FDA. The purpose of this section is to review the most important concepts and methods related to clustering and classification of functional data, via illustrative examples, to help the reader have a picture about the fields.

Clustering of functional data

Similarly to the classical multivariate data analysis, clustering functional data concerns an upsupervised learning process that group a set of functional data into homogenous subgroups, i.e., clusters, such that whose members are more similar than across clusters with respect to some metric. The concepts and algorithms developed for clustering multivariate data are therefore can be extended to functional data, certainly with additional considerations regarding to, e.g., the infinite dimension, the approximations of distance measures, etc.

There are various algorithms for clustering functional data. They can be categorized into four groups based on their approaches: raw-data clustering, two-stage clustering, non-parametric clustering, and model-based clustering.

Raw-data clustering. This is the most simple approach that uses directly the observed values of the functions as vectorial data for clustering. However, due to the so-called "curse of dimensionality", i.e., these raw-data are typically of high-dimension, the techniques in this direction mostly consist of a dimension reduction step followed by an algorithm for clustering vectorial data. The most often used dimension reduction methods includes principal component analysis (PCA) and factor analysis (FA). For a complete review of PCA and FA methods, we would refer reader to [START_REF] Bishop | Pattern recognition and machine learning[END_REF], Chapter 12.

2.8 CLUSTERING AND CLASSIFICATION OF FUNCTIONAL DATA Two-stage clustering. This approach consists of a filtering step followed by a classical clustering algorithm for finite dimensional data. Although it seems to be similar to the raw-data clustering that, their difference is that: rather than using PCA or FA, in the filtering step, we use the tools of FDA to represent the functional data. In particular, the representations via Spline basis (e.g., in [START_REF] Abraham | Unsupervised curve clustering using b-splines[END_REF]; [START_REF] Rossi | Clustering functional data with the som algorithm[END_REF]) and FPCA (e.g., in Jacques and Preda (2013a); [START_REF] Peng | Distance-based clustering of sparsely observed stochastic processes, with applications to online auctions[END_REF]) are the most common choices in this direction. Once the coefficients are obtained, they are fitted to the usual clustering algorithms such as k-means (see, e.g., [START_REF] Hartigan | Clustering Algorithms[END_REF]; [START_REF] Hartigan | A K-means clustering algorithm[END_REF]) or Self-Organised Map (SOM) [START_REF] Kohonen | Self-Organizing Maps[END_REF].

Non-parametric clustering. This terminology covers all clustering techniques that involve nonparametric algorithms. Analogously for multivariate data, the techniques in this approach require a dissimilar measure for two objects. The common choice is the measure d ℓ defined based on the L 2 -distance in (2.14), given by

d ℓ (X i , X j ) = T (X (ℓ) i (t) -X (ℓ) j (t)) 2 dt 1/2
, where X (ℓ) i denotes the ℓ-derivative of X. Note that, with ℓ = 0 the measure d ℓ becomes exactly the distance defined in (2.14), which is equivalent to using the raw-data approach.

In this direction, there are several interesting works, for example, in [START_REF] Ieva | Multivariate functional clustering for the morphological analysis of electrocardiograph curves[END_REF] the authors used K-means combing with d 0 , d 1 , and (d 2 0 + d 2 1 ) 1/2 to cluster ECG 2 traces, aim at early detect the heart failure, replacing a traditional clinical observation. In [START_REF] Tarpey | Clustering functional data[END_REF], K-means and d 0 were used for clustering Gaussian processes. Moreover, cluster centres are shown to be a linear combinations of the FPCA eigenfunctions. Figure 2.6 shows the result of doing clustering for Tecator data (see Appendix A.2) using d 0 , d 1 and d 2 measures combined with K-means algorithm.

Model-based clustering. Let us briefly recall the ideas of the method of model-based clustering for multivariate data. In the view of model-based clustering, the observations {x i } n i=1 are assumed to be generated according to a K-component mixture distribution. Let z i = (z i1 , . . . , z iK ) be the unknown one-hot vector that encodes the membership of x i , and the aim here is to estimated z i 's. This can be done by treating z i 's as missing data and maximizing the so-called complete-data log-likelihood function given by where Ψ = (θ 1 , . . . , θ K , π 1 , . . . , π K ) is the model parameter to be estimated, f k (x i |θ k ) is the density corresponding to the kth cluster, parameterized by θ k , and π k is the probability that an observation belongs to the kth cluster. The common choice of f k (x i |θ k ) is the multivariate normal distribution with mean µ k and variance Σ k (i.e., θ k = (µ k , Σ k )), for example in [START_REF] Banfield | Model-based gaussian and non-gaussian clustering[END_REF], [START_REF] Celeux | Gaussian parsimonious clustering models[END_REF]. The inference is usually carried by an iterative procedure such as EM algorithm. Finally, the cluster labels are estimated using the Maximum A Posteriori (MAP) rule

L(Ψ; x 1 , . . . , x n , z 1 , . . . , z n ) = n i=1 log K k=1 P(z ik = 1)f k (x i ; θ k ) z ik = n i=1 K k=1 z ik log(π k f k (x i ; θ k )),
ẑi = argmax k∈[K] πk f k (x i ; θk ) K k=1 πk f k (x i ; θk )
.

In this model-based direction for clustering functional data, the techniques are generally consists of two tasks: one performs function expansion (using basis functions or FPCA) and one performs model-based clustering. However, as highlighted in [START_REF] Jacques | Functional data clustering: A survey[END_REF], it differs from the two-stage direction in that it performs the two mentioned tasks simultaneously, whereas in twostage techniques, clustering is applied once the functions were already "converted" to vectors of coefficients.

The first model-based clustering algorithm for functional data was proposed by [START_REF] James | Clustering for sparsely sampled functional data[END_REF], in which the spline expansion coefficients of the curves are assumed to be distributed according to a mixture Gaussian distribution with cluster-specific means µ k and common variance Σ. Later, [START_REF] Giacofci | Wavelet-based clustering for mixed-effects functional models in high dimension[END_REF] proposed using wavelet expansion as an alternative, since the spline expansion is only convenient when the curves are regular but not appropriate for irregular (e.g., peak-like) data. The authors assumes a Gaussian model on the wavelet decomposition of the curves whose parameters are estimated via MLE, and applies the approach on mass spectrometry data, as well as proposes an original application on microarray CGH3 data. In Chapter 3, via the FME model, we also propose a model-based clustering technique and apply for Canadian weather and DTI4 data.

Classification of functional data

Similar to regression and clustering for functional data, due to the practical needs, classification for functional data also has a long history of study with many notable applications. Many approaches have been proposed for functional data classification, they are mostly based on functional regression models that treat class labels as responses and the observed functional data as predictors. For example, in [START_REF] Müller | Generalized functional linear models[END_REF], the authors used the PCA in combination with logistic regression model for classification problems, and applied to simulated and medfly data. The generalized linear model has also been extended to functional data [START_REF] James | Generalized linear models with functional predictor variables[END_REF][START_REF] Müller | Functional modelling and classification of longitudinal data*[END_REF][START_REF] Müller | Generalized functional linear models[END_REF]. Recently, in [START_REF] Mousavi | Functional logistic regression: a comparison of three methods[END_REF], the authors applied a wavelet basis expansion in combination with multinomial functional regression model to classify acceleration signals in the problem of lameness diagnosis for horses, and periodogram signals in phoneme data. In their proposed model, the parameter vector is estimated via regularized MLE which produces a certain level of sparsity. In Chapter 4, our proposed FME model for classification also use the phoneme data as a benchmark and compare the result with that of [START_REF] Mousavi | Functional logistic regression: a comparison of three methods[END_REF].

Let X(t), t ∈ T , be a functional predictor and y be its categorical response, y ∈ [G]. In multinomial functional regression model, the probability that y belongs to category g is modeled by .44) where β g,0 are the unknown constant intercepts, and β g (t), t ∈ T , are the unknown coefficient functions to be estimated. Similar to the softmax modeling for the gating function in ME models, model (2.44) is over-parameterized since it can get a same probability by adding a constant to all of the intercepts and coefficient functions. Therefore, one can also constraint, e.g., the parameter responsible for gth class (β g,0 , β g (t)) to be null, which yields the model

P(y = g|X(t)) = exp{β g,0 + T β g (t)X(t)dt} G g ′ exp{β g ′ ,0 + T β g ′ (t)X(t)dt} , ( 2
P(y = g|X(t)) = exp{β g,0 + T β g (t)X(t)dt} 1 + G-1 g ′ exp{β g ′ ,0 + T β g ′ (t)X(t)dt} •
We will use this parameterization for the functional gating function of the FME model proposed in Chapter 3. Alternatively, over-parameterization in (2.44) can be dealt by the LASSO regularization as can be seen in [START_REF] Mousavi | Functional logistic regression: a comparison of three methods[END_REF]. Classification of functional data can also be approached using discriminant analysis. The early works include [START_REF] Hall | A functional data-analytic approach to signal discrimination[END_REF] and James and Hastie (2001a), where linear discrimination analysis (LDA) was applied on the scores obtained from PCA technique, and on the coefficients from spline expansion, respectively. In [START_REF] Ferraty | Curves discrimination: a nonparametric functional approach[END_REF], the authors proposed a non-parametric approach for functional data. Later, [START_REF] Chang | Functional data classification: a wavelet approach[END_REF] developed a kernel-based non-parametric approach, applied to classify the major depressive disorders based on positron emission tomography images. [START_REF] Berlinet | Functional supervised classification with wavelets[END_REF] proposed a supervised wavelet-based functional data classification method. Along this line of discriminant analysis, theoretical results have been also developed for linear and quadratic decision boundaries (Delaigle andHall, 2012, 2013).

Preliminary on divergence

In order to prepare for Chapter 5 where we propose a distributed learning approach for ME models, in this section we briefly present some preliminaries on general notations of divergence and distance on the space of probability measures.

Let us recall the general notions of divergence and distance defined on the space of probability 2.10 SUMMARY measures.

Definition 2.9.1 (Divergence). Let Θ be a space and ρ(•, •) be a bivariate function defined on

Θ. Then ρ(•, •) is called a divergence if ρ(θ 1 , θ 2 ) ⩾ 0 for all θ 1 , θ 2 ∈ Θ, with equality holding if and only if θ 1 = θ 2 .
Kullback-Leibler (KL) is one of the most used divergences in machine learning to measures the similarity between two probability distributions, or can be interpreted as the information gain from using p instead of q [START_REF] Cover | Elements of Information Theory[END_REF]. Definition 2.9.2 (KL-divergence). Let p and q are two probability distributions, KLdivergence between p and q is defined by

KL (p∥q) = p(u) log p(u) q(u) du.
In our case, u being (x, y), the KL-divergence KL (p∥q) can be factorized as

KL (p(x, y)∥q(x, y)) = p(x, y) log p(x, y) q(x, y) dxdy = p(x) log p(x) q(x) dx + p(x) p(y|x) log p(y|x) q(y|x) dxdy = KL (p(x)∥q(x)) + E p(x) KL (p(y|x)∥q(y|x)) , (2.45)
where KL (p(y|x)∥q(y|x)) := p(y|x) log p(y|x) q(y|x) dy.

(2.46)

The last term in (2.45) is known as conditional KL-divergence, which is the expected value of the KL-divergence between conditional distributions p(y|x) and q(y|x) where the expectation is taken with respect to p(x).

Summary

This chapter has introduced the main framework we will working on in the next chapters, i.e., ME models. Notably, the general notion, the modelings, as well as the estimation procedure for gating and expert networks have been presented, these will make the following chapters more consistent and easier to follow. This chapter has also present the necessary tools and concepts for FDA such as data pre-processing, FPCA, functional linear models. Interested reader is referred CHAPTER 2 STATE OF THE ART to other comprehensive books on functional regression models [START_REF] Kokoszka | Introduction to functional data analysis[END_REF], [START_REF] Ramsay | Functional Data Analysis[END_REF], [START_REF] Ramsay | Spatial functional data analysis[END_REF]. While the ME models with vectorial setting has been widely studied in literature, their functional versions are still lack of investigations and applications. This thesis therefore proposes, in the subsequent chapters, three new ME models to deal with situations where data are available under functional forms.

Chapter 3

Functional mixtures of experts

Introduction

Recall that the general notion of a ME is given by

ME(y|x) = K k=1 Gate k (x) Expert k (y|x). (3.1)
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 high-dimensional heterogeneous data, which have been studied with the ME modeling in (3.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 | 2.1 Visualization of the architecture of a ME model. Here, the gating network and each of K expert receive x as input, then[END_REF]Goldsmith et al., 2011b) 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] including Lasso-regularized mixtures for functional data [START_REF] Bibliography 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 chapter, 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. Firstly, we introduce in Section 3.2 a new family of functional ME models 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 highdimensional setting of the introduced FME model, we develop in Section 3.2.6 a Lasso-regularized approach, which consists of a penalized MLE estimation via a 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.3 and extended FME model, which is constructed upon a sparse and highly-interpretable regularization of the functional expert and gating parameters.

FUNCTIONAL MIXTURES-OF-EXPERTS (FME)

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 3.4 on several simulated scenarios and on real-world data sets, in both clustering and non-linear regression. This chapter is a detailed version of the paper [START_REF] Chamroukhi | Functional mixturesof-experts[END_REF].

Functional Mixtures-of-Experts (FME)

We wish to derive and fit new ME models in presence of functional predictors and functional responses. In this chapter, 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.

ME with functional predictor and scalar response

Let {X i (t), Y i } n
i=1 , be a sample of n 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 Expert z (y|x) in (3.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], (3.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. In this context, the response Y is related to the entire trajectory of X(•). Let β = {β z,0 , β z (t), t ∈ T } K z=1 represent the set of unknown functional parameters for the experts network. Now consider the modeling of the gating network in the proposed functional ME 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 mixture-of-experts modeling: An overview[END_REF], Xu et al. (1994a), Chamroukhi et al. (2019)). A standard choice as in [START_REF] Jacobs | Adaptive mixtures of local experts[END_REF] to model the gating network Gate z (x) in (3.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.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.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. (3.4)
The objective is to estimate the functional parameters α and β of the FME model defined by (3.2)- (3.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), (3.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

FUNCTIONAL MIXTURES-OF-EXPERTS (FME)

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], (3.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) (3.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 Goldsmith et al. (2011b), [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 (3.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 (3.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.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 (3.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

X i (t), i ∈ [n], we use X i (t) established in
q (t) = [b 1 (t), b 2 (t), . . . , b q (t)] ⊤ , (q ⩽ r) via the projection α z (t) = ζ ⊤ z b q (t), (3.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 (3.6) and (3.8) of X(t) and α z (t), respectively, in the linear predictor h z (.) defined in (3.4) for i ∈ [n], the latter is thus approximated

CHAPTER 3 FUNCTIONAL MIXTURES OF EXPERTS as h z (U i (t), t ∈ T ; α) = α z i ,0 + ζ ⊤ z i r i , (3.9)
where

r i = T b r (t)b q (t) ⊤ dt ⊤ x i .
Thus, following its definition in (3.3), the functional softmax gating network is approximated as

π k (r i ; ξ) = exp {α k,0 + ζ ⊤ k r i } 1 + K-1 k ′ =1 exp {α k ′ ,0 + ζ ⊤ k ′ r i } (3.10) where ξ = (α 1,0 , ζ ⊤ 1 ), . . . , (α K-1,0 , ζ ⊤ K-1 ) ⊤ ∈ R (q+1)(K-1
) is the unknown parameter vector of the functional gating network to be estimated.

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 ( 3.2)- (3.3). Gathering (3.6) and (3.7), the stochastic representation (3.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], (3.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.
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 ), (3.12)
where ϕ(.; µ, v) is the Gaussian probability density function (pdf) 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 (3.12) and (3.10) in the ME model (3.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 ), (3.13)
where Ψ = (ξ ⊤ , θ ⊤ 1 , . . . , θ ⊤ K ) ⊤ is the parameter vector of the model to be estimated.

Maximum likelihood estimation via the EM algorithm

The FME model (3.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

FUNCTIONAL MIXTURES-OF-EXPERTS (FME)

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 (eg. 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 (3.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 ), (3.14)
via the EM algorithm. As detailed in Appendix B.1, 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 (3.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) ) • (3.15) ■ M-step.
Update the value of the parameter vector Ψ by maximizing the Q-function (B.2) with respect to Ψ . The maximization 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 , for each of the K experts.

Updating the gating network's parameters ξ consists of maximizing w.r.t ξ the part of (B.2) that is a function of ξ. Since we use a softmax-gated expert network in (3.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 (B.3) 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) , (3.16)
where H(ξ; Ψ (s) ) and g(ξ; Ψ (s) ) are, respectively, the Hessian matrix and the gradient vector of Q(ξ; Ψ (s) ), and are provided in Appendix B.1. 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 (3.16) until there is no significant change in Q(ξ; Ψ ). The maximization then provides ξ (s+1) for the next EM iteration.

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 (3.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 = (n (s) k ) -1 n i=1 τ (s) ik (y i -x ⊤ i η (s) k ); η (s+1) k = n i=1 τ (s) ik x i x ⊤ i -1 n i=1 τ (s) ik (y i -β (s+1) k,0 )x i , σ 2 k (s+1) = (n (s) k ) -1 n i=1 τ (s) ik y i -(β (s+1) k,0 + x ⊤ i η (s+1) k ) 2 ,
(3.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 mixture-of-experts 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 functional ME 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]Huynh and Chamroukhi, 2019a) and enjoy good theoretical properties [START_REF] Nguyen | An l 1-oracle 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].

FUNCTIONAL MIXTURES-OF-EXPERTS (FME)

ℓ 1 -regularization and the EM-Lasso algorithm

We propose an ℓ 1 -regularization of the observed-data log-likelihood (3.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, (3.18) where log L(Ψ ) is the observed-data log-likelihood of Ψ defined by (3.14), and Pen λ,χ (Ψ ) is a LASSO regularization term encouraging sparsity for the expert and the gating network parameters, defined by

L(Ψ ) = log L(Ψ ) -Pen λ,χ (Ψ ),
Pen λ,χ (Ψ ) = λ K k=1 p j=1 |η k,j | + χ K-1 k=1 q j=1 |ζ k,j |, (3.19)
where λ and χ are positive real values representing tuning parameters. The maximization of (3.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.2. 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 ℓ 1 -penalized log-likelihood (3.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 (3.15).

■ M-step. In this regularized MLE context, the parameter vector Ψ is now updated by maximizing the regularized Q-function (B.5), 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 (k ∈ [K]), w.r.t the expert network parameters θ 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 , (3.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 (eg. see [START_REF] Hastie | Statistical learning with sparsity: the lasso and generalizations[END_REF], Huynh and Chamroukhi (2019a)), 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 (B.7) 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 softthresholding operator defined by S χ (u) = sign(u)(|u|χ) + and (v) + 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 10 -6 . 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.

FUNCTIONAL MIXTURES-OF-EXPERTS (FME)

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 , (3.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.2, updating (β k,0 , η k ) according to (3.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 (3.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 rely directly the original functional inputs, to the output, assuming some of their values is zero is not 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 3.2.5 and the regularized FME model fitted by EM-Lasso algorithm, in Section 3.2.7. 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.

Motivation

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.

Interpretable sparse regularization

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

INTERPRETABLE FUNCTIONAL MIXTURE OF EXPERTS (IFME)

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 ), D d 2 b p (t 1 ), D d 2 b p (t 2 ), . . . , D d 2 b p (t p ) ⊤ (3.22)
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. Similarly, let A q = [A [d 1 ] q A [d 2 ]
q ] ⊤ ∈ R 2q×q be the corresponding matrix defined for the gating network. Now, if we consider the following representation for the expert network coefficient function

γ z i = A p η z i (3.23) 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 ) ⊤ and γ [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) in (3.7), then γ

[d 1 ] z i (resp. γ [d 2 ]
z i ) provides an approximation to β

(d 1 ) z i (t), the d 1 th derivative of β z i (t) (resp. β (d 2 )
z i (t), the d 2 th derivative of β z i (t)). Hence, enforcing sparsity in γ z i constrains β

(d 1 ) z i (t) and β (d 2 )
z i (t) to be zero at most time points. Similarly, if we consider the following representation for the gating network coefficient function

ω z i = A q ζ z i (3.24)
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) in (3.8), then we can derive the same interpretation for the coefficient vector ω z i and the gating parameter function α z i (t). The constraints (3.23) and (3.24) imply

η 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 , 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 , 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 (3.11) for the experts model and the linear predictor definition in (3.9), we obtain the following iFME model construction,

Y i |u i (•) = β z i ,0 + γ [d 1 ] z i ⊤ v i + ε ⋆ i , (3.25)
h z i (u i (•); α) = α z i ,0 + ω [d 1 ] z i ⊤ s i , (3.26) 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 , (3.27)
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 (3.25) and (3.26), 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 ) (3.28) 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 } , (3.29)
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 (3.28) and (3.29) as for (3.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 ), (3.30)
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 λ,χ (Ψ ) (3.31)
with the conditional iFME density f (y i |u i (•); Ψ ) is defined in (3.30) 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 ), (3.32) 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. 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. It is worth noting two points: first, 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; and second, maximizing the penalized log-likelihood function (3.31) with penalization in (3.32) in iFME model must be coupled with the constrains (3.27). 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) ) , (3.33)
where f (y i |u i (•); Ψ (s) ) is now calculated according to the iFME density given by (3.30).

■ 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 ), (3.34) subject to ω

[d 2 ] k = A [d 2 ] q A [d 1 ] q -1 ω [d 1 ] k , ∀k ∈ [K -1], (3.35)
where Q(w;

Ψ (s) ) = n i=1 K-1 k=1 τ (s) ik α k,0 + ω [d1] k ⊤ s i -log 1 + K-1 k ′ =1 exp{α k ′ ,0 + ω [d1] 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 (3.34) is now coupled with an additional constraint (3.35), 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 (3.34) 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 (3.35). 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 , (3.36)
where t) ; s i ) are the weights and

w ik = π k (w (t) ; s i ) 1 -π k (w (
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 (3.36)

are given in Appendix B.3.1. Therefore, if ( α k,0 , ω [d 1 ] k , ω [d 2 ]
k ) is an optimal solution to problem (3.36), 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 (3.34).

INTERPRETABLE FUNCTIONAL MIXTURE OF EXPERTS (IFME)

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 .

(3.37)

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 (3.36). The corresponding technical details can be found in Appendix B. 3.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 (3.37) (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 (3.31), 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 , (3.38)
where ω

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

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

Clustering and non-linear regression with FME models

Once the model parameters have been estimated, 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 (3.15) for FME or by (3.33) for the iFME model, uing 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.

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 , σ 2 k )
for the FME model (3.13), and by

y i |u i (•) = K k=1 π k (s i ; w)( β k,0 + γ [d 1 ] k ⊤ v i , σ 2 k )
for the iFME model (3.30).

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 (3.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 (3.32)).

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 well-known BIC criterion [START_REF] Schwarz | Estimating the dimension of a model[END_REF] or, in our context, its extension, called modified BIC [START_REF] Bibliography Städler | ℓ1-penalization for mixture regression models[END_REF] defined as

mBIC = L( Ψ ) -df( Ψ ) log n 2 , (3.39)
where Ψ is the obtained log-likelihood estimator (for the FME model) or penalized log-likelihood estimator (for the FME-Lasso and iFME models), and the number of degrees of freedom df( Ψ ) is

(ζ) = 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 studies

We study the performances of the FME, FME-Lasso, and iFME models in regression and clustering problems by considering simulated scenarios and real-world 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 mixture-ofexperts 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 predictions error (RPE) and the correlation (Corr) index 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 iy 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 adjusted Rand index (ARI), and the clustering error (ClusErr), 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 the model parameters using (3.7, 3.8) for both FME and FME-Lasso models, and using (3.38) for iFME model. The values of these criteria are averaged over N Monte Carlo runs (N = 100 for simulation, for the real-world 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) 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 .

MSE( g(•)) = 1 m m j=1 (g(t j ) -g(t j )) 2 , ( 3

Simulation studies

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)

β = {β k,0 , β k (t), σ 2 k } K k=1 and α = {α k,0 , α k (t)} K k=1
, t ∈ T , we first construct the functional predictors X i (•) for i ∈ [n] using the representation defined in (3.5), with detailed parameterization in (3.42). Then, for each i ∈ [n], conditional on the functional predictor X i (•), a hidden categorical random variable (3.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:

Z i ∈ [K] is generated following the multinomial distribution M 1, π 1 (X i (t); α), . . . , π K (X i (t); α) , where π k (X i (t); α) for k ∈ [K] is given by
Y i |Z i = z i , X i (t) ∼ N β z i ,0 + T X i (t)β z i (t)dt; σ 2 z i , Z i |X i (t) ∼ M 1, π 1 (X i (t); α), . . . , π K (X i (t); α) . (3.41)
The true generated data (X i (•), Y i , Z i ) are used for evaluating the prediction and clustering performance. Finally, to mimic real-world 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

U i (t) = X i (t) + δ i (t),
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 3.1).

Simulation parameters and experimental protocol

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

EXPERIMENTAL STUDIES

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

β 1 (t) =      -50(t -0.5) 2 + 4 if 0 ⩽ t < 0.3, 0 if 0.3 ⩽ t < 0.7, 50(t -0.5) 2 -4 if 0.7 ⩽ t ⩽ 1, β 2 (t) = -β 1 (t), β 3 (t) = 100(t -0.5) 2 -10, 0 ⩽ t ⩽ 1, (β 1,0 , β 2,0 , β 3,0 ) ⊤ = (-5, 0, 5) ⊤ , (σ 2 1 , σ 2 2 , σ 2
3 ) ⊤ = (5, 5, 5) ⊤ , and the functional gating network parameters are given by

α 1 (t) = 80(t -0.5) 2 -8, α 2 (t) = -α 1 (t), α 3 (t) = 0, 0 ⩽ t ⩽ 1, (α 1,0 , α 2,0 , α 3,0 ) ⊤ = (-10, -10, 0) ⊤ .
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 zero. Here, "≺" is the lexicographical order on R p+2 .

As it can be seen in Figure 3.1, 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 considering this network, 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,

X i (t) = x ⊤ i b(t), t ∈ T , (3.42) 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 ∈ {100, 300}, and two levels of noise σ 2 δ ∈ {1, 4}. The resulting four considered scenarios of simulated data are presented in Table 3.1. For each considered scenario, we generate N = 100 datasets, each dataset contains n = 200 pairs of (U

i (t), Y i ), i ∈ [n].

EXPERIMENTAL STUDIES

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 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 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 3.2, 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 3.3 shows the clustering results obtained by the models with highest values of the modified BIC criterion, on a dataset selected in scenario S1. Here, we plotted the responses 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 3.3. 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. 3: Scatter plots of y i against U i (t 1 ) (top panels) and U i (t 150 ) (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.

Comparison with functional regression mixtures (FMR). Finally, we compare our proposed models with the functional mixture regression (FMR) model proposed in Yao et al. (2010). 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. Parameter estimation performance. To evaluate the parameter estimation performance, we consider the functional parameter functions estimation error as defined by (3.40). 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.4. 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 3.5.

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 3.6 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. Table 3.7 shows a comparison between three different initialization strategies for the EM algorithm. In particular, the initialization strategies regards only the gating network parameter vector, since for the expert network a K-means approach is more suitable. With the "rand." 4.17 strategy, all coefficients of the gating network parameter vectors are drawn randomly in N (0, 1). With the "zeros" strategy, they are initialized as zeros. And with the "LR" we basically perform a logistic regression where the predictors are the design vectors associated with the gating network (i.e., r i for FME, FME-Lasso, and s i for iFME model), and the responses are the labels resulted by K-means on the gating design matrix. 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.4, 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. Here, the square points correspond to the highest values.

β 1 (•) β 2 (•) β 3 (•) α 1 (•) α 2 (•)
T 2 on T 1 ∪ T 3 on T β (0) 1 (•) β (0) 2 (•) β (3) 1 (•) β (3) 2 (•) β (3) 3 (•) α (3) 1 (•) α (3) 2 (•)
β 1,0 β 2,0 β 3,0 α 1,0 α 2,0 σ 1 σ 2 σ 3 True value -5 0

Application to real-world data

In this section, we apply the FME, FME-Lasso and iFME models to two well-known real-world 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. 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 3.8 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. Finally, in this cross-validation study, the obtained number of components was K = 4 for FME and iFME, while for the FME-Lasso model, the selected one is K = 5.

Figure 3.6 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 3.7. 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 3.8 displays the values of modified BIC for varying number of expert component for the proposed models on Canadian weather data. According to these plots, 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 number of correct answers out of 60 questions, and thus ranges from 0 to 60. In our context, the FAP curves serve as the noisy functional predictors U i (•) and the PASAT scores serve as the scalar responses Y i . So this dataset consists of n = 99 pairs

(U i (•), Y i ), i ∈ [n],
with each U i (•) contains m = 93 fractional anisotropy values. A visualization of the predictors (i.e., the FAP curves) and the responses (i.e., the PASAT Scores) can be seen in Figure 3.9 (top-left panel) and Figure 3 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 3.4 EXPERIMENTAL STUDIES belonging to it. We accordingly fix K = 2 in our models. Figure 3.10 displays the obtained results for each of the three models, and Figure 3.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.

Estimated groups

Group 1 Group 2

Figure 3.10: 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.

In Figure 3.10, we have the three following observations. First, as it can be observed in Figure 3.10 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 3.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 (CVRPE) defined by CVRPE = n i=1 (y iy

(-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 Table 3.9. Note that, for comparison, in [START_REF] Ciarleglio | Wavelet-based scalar-on-function finite mixture regression models[END_REF], the CVRPE of their WBFMR model is 0.0315 and of the wavelet based functional linear model (FLM) is 0.0723.

CVRPE FME

0.0273 FME-Lasso 0.0280 iFME (with 0th and 3nd derivatives) 0.0271 iFME (with 0th and 2nd derivatives) 0.0267 Table 3.9: CVRPEs of the models on the DTI data.

Finally, we present in Figure 3.11 the 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. 3.5 SUMMARY

Summary

The first algorithm for mixtures-of-experts constructed upon functional predictors is presented in this chapter. Beside the classic maximum likelihood parameter estimation, we proposed two other regularized versions that allow sparse and 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 (u), Y i (t)} 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 (t) = β z i ,0 (t) + t X i (u)β z i (t, u)du + ε i (t), 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 .

Introduction

In this chapter, we have two goals. First, we extend the FME and iFME models proposed in previous chapter to the case of categorical responses, i.e., for multiclass classification problems. Second, we develop ME models to deal with the situations where the responses are functions.

To achieve the first goal, in Section 4.2 we model the experts in FME models by the functional multinomial logistic regression model, and develop the corresponding conditional density functions for the FME and iFME models. The numerical experiments are also presented to illustrate the performance of the extended models.

FME FOR CLASSIFICATION

In Section 4.3, a ME model for function-on-function regression is proposed. It can be viewed as an extension of the FME model to the case of functional responses. We refer to it by function-onfunction FME model (FF-FME). The model estimation is performed via maximizing a penalized log-likelihood function, with the penalized term is constructed such that estimated kernels are smoothened. An efficient EM-like algorithm is developed to estimate the model parameters. The proposed model is examined on both simulated and real-world data.

Finally, in Section 4.4, to complete the family of FME models, we present a ME model for function-on-scalar regression. We refer to this model by function-on-scalar FME model (FS-FME).

A penalized MLE is considered, in which the penalization is given to the roughness of the functional parameters. We also develop the corresponding EM-like algorithm for model estimation. The model performance is evaluated on simulated data, which show promising results for other applications.

FME for classification

In this section, we extend the FME framework for multiclass classification, derive adapted EM-like algorithms to obtain sparse and interpretable fit of the gating and experts network coefficients functions. Let {X i (t), t ∈ T ; Y i } n i=1 , be a sample of n i.i.d. data pairs where Y i ∈ {1, . . . , G} is the class label of a functional predictor X i (•), G being the number of classes. In this case of functional inputs, a natural choice to model the conditional distribution Expert k (y|x) = P(Y = y|X i (•)) in (2.1) is to use the functional multinomial logistic regression modeling, e.g., see [START_REF] Müller | Generalized functional linear models[END_REF]; James (2002), that is Similarly, a typical choice for the functional gating network Gate k (x) = P(Z = k, X(•)) in ( 2.1), where Z ∈ [K] is a hidden within-class clustering label, acting as weights for potential clusters {k} in the heterogeneous functional inputs X(•) and which we denote as π k (X(•)), is to use a functional softmax function defined by

P (y i |X i (•); β k ) = G g=1 exp β kg,0 + T X i (t)β kg (t)dt 1 + G-1 g ′ =1 exp β kg ′ ,0 + T X i (t)β kg ′ (t)dt y ig , ( 4
π k (X i (•); α) = exp{α k,0 + T X i (t)α k (t)dt} 1 + K-1 k ′ =1 exp{α k ′ ,0 + T X i (t)α k ′ (t)dt} , ( 4
P (y i |X i (•); ψ) = K k=1 π k (X i (•); α)P (y i |X i (•); β k ), ψ = (α, β 1 , . . . , β K ). (4.3)
where ψ = (α, β 1 , . . . , β K ) is the unknown parameter of the model. Here the coefficient functions α k (•), β kg (•) have the roles as the coefficient vectors as in the classic ME model for classification.

Smooth functional representation

In practice, X i (•) is observed at a finite but large number of points on T ⊂ R. In the perspective of parameter estimation, this results in estimating a very large number of coefficients β and α.

In order to handle this high-dimensional problem, we consider a usual approach that projects the predictors and coefficient functions onto a family of reduced number of basis functions. Let b r (t) = [b 1 (t), . . . , b r (t)] ⊤ be a r-dimensional basis (B-spline, Wavelet, ...). Then, with r, p, q ∈ N sufficiently large, one can approximate X i (•), α k (•) and β kg (•) respectively by

X i (t) = x ⊤ i b r (t), α k (t) = ζ ⊤ k b p (t), β kg (t) = η ⊤ kg b q (t). (4.4)
Here,

x i = (x i1 , . . . , x ir ) ⊤ , with x ij = T X i (t)b j (t)dt for j ∈ [r], is the vector of coefficients of X i (•) in the basis b r (t), ζ k = (ζ k,1 , . . . , ζ k,p
) ⊤ , and η kg = (η kg,1 , . . . , η kg,q ) ⊤ are the unknown coefficient vectors associated with the gating coefficient function α k (•) and the expert coefficient function β kg (•) in the corresponding basis. In our case, we used B-spline bases. Using the approximation of X i (•) and α k (•) in (4.4), the functional softmax gating network (4.2) can be represented by

π k (r i ; ξ) = exp{α k,0 + r ⊤ i ζ k } 1 + K-1 k ′ =1 exp{α k ′ ,0 + r ⊤ i ζ k ′ } , (4.5)
where r i = T b r (t)b p (t) ⊤ dt ⊤ x i is the design vector associated with the gating network and -1) is the unknown parameter vector of the gating network, to be estimated.

ξ = ((α 1,0 , ζ ⊤ 1 ), . . . , (α K-1,0 , ζ ⊤ K-1 )) ∈ R (p+1)(K
In the same manner, using the approximations of X i (•) and β kg (•) in (4.4), the expert conditional distribution (4.1) can be represented by

P (y i |x i ; θ k ) = G g=1 exp β kg,0 + x ⊤ i η kg 1 + G-1 g ′ =1 exp β kg ′ ,0 + x ⊤ i η kg ′ y ig , (4.6)
where x i = T b r (t)b q (t) ⊤ dt ⊤ x i is the design vector associated with the expert network, and

θ k = (θ ⊤ k1 , . . . , θ ⊤ k,G-1 ) ⊤ , with θ kg = (β kg,0 , η ⊤ kg ) ⊤ ∈ R q+1 for g ∈ [G -1],
is the unknown parameter vector to be estimated of the expert distribution k.

Finally, combining (4.5) and (4.6), the conditional distribution P (y i |X i (•); ψ) in ( 4.3) can be 4.2 FME FOR CLASSIFICATION rewritten as

P (y i |X i (•); Ψ ) = K k=1 π k (r i ; ξ)P (y i |x i ; θ k ),
where Ψ = (ξ ⊤ , θ ⊤ 1 , . . . , θ ⊤ K ) ⊤ is the unknown parameter vector of the model.

Parameter estimation via EM algorithm

A maximum likelihood estimate (MLE) Ψ of Ψ can be obtained by using the EM algorithm for ME model for classification with vector data as in [START_REF] Chen | Improved learning algorithms for mixture of experts in multiclass classification[END_REF]. We will refer to this approach as FME-EM.

Similarly to the FME model, to encourage sparsity in the model parameters Ψ , one can perform penalized MLE by using the EM-Lasso algorithm as in [START_REF] Huynh | Estimation and feature selection in mixtures of generalized linear experts models[END_REF]. We refer to this approach as FME-EM-Lasso.

Interpretable FME for classification

Although fitting the FME model via EM-Lasso can accommodate sparsity in the parameters, it unfortunately does not ensure the reconstructed coefficient functions α k (•) and β kg (•) are sparse and enjoy easy interpretable sparsity.

To obtain interpretable and sparse fits for the coefficient functions, we simultaneously estimate the model parameters while constraining some targeted derivatives of the coefficient functions to be zero [START_REF] Chamroukhi | Functional mixturesof-experts[END_REF]. The construction of the interpretable FME model which we will fit with an adapted EM algorithm, is as follows.

First, in order to calculate the derivative of the gating coefficient functions α k (•), let A p be the matrix of approximate d 1 th and d 2 th derivative of b p (t), defined as in (3.22) by

A p = [A [d 1 ] p A [d 2 ] p ] ⊤ = D d 1 b p (t 1 ), . . . , D d 1 b p (t p ), D d 2 b p (t 1 ), . . . , D d 2 b p (t p ) ⊤ ,
where D d is the dth finite difference operator. Here A

[d j ] p
is a square invertible matrix and A p ∈ R 2p×p . Similarly, to calculate the derivatives of the expert coefficient functions

β kg (•), let A q = [A [d 1 ] q A [d 2 ]
q ] ⊤ ∈ R 2q×q be the corresponding matrix defined for the b q (t)'s. Now, if we define ω k = A p ζ k and denote ω k = (ω

[d 1 ] k ⊤ , ω [d 2 ] k ⊤ ) ⊤ , then ω [d 1 ] k
and ω k (•) to be zero at most of time points. Similarly, if we define γ kg = A q ζ k and denote by γ kg = (γ

[d 1 ] kg ⊤ , γ [d 2 ] kg ⊤
) ⊤ , then we can derive the same regularization for the coefficient functions

β kg (•).
From the definitions of ω k and γ kg we can easily get the following relations:

   ζ k = A [d 1 ] p -1 ω [d 1 ] k and ω [d 2 ] k = A [d 2 ] p A [d 1 ] p -1 ω [d 1 ] k η kg = A [d 1 ] q -1 γ [d 1 ]
kg and γ

[d 2 ] kg = A [d 2 ] q A [d 1 ] q -1 γ [d 1 ]
kg .

(4.7a) (4.7b)

Plugging the relation (4.7a) into (4.5) one gets the following new representation for π k (r i ; ξ)

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

s i = (A [d 1 ] p -1
) ⊤ r i is now the new design vector 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. Similarly, plugging (4.7b) into (4.6) one obtains the new representation for P (y i |x i ; θ k ):

P (y i |v i ; Γ k ) = G g=1   exp β kg,0 + v ⊤ i γ [d 1 ] kg 1 + G-1 g ′ =1 exp β kg ′ ,0 + v ⊤ i γ [d 1 ] k ′ g   y ig , (4.9) in which, v i = (A [d 1 ] q -1
) ⊤ x i is now the new design vector and Γ k = (β kg,0 , γ

[d 1 ] k ′ g ⊤ ) ⊤
is the unknown parameter vector of the expert network. Finally, gathering the gating network (4.8) and the expert network (4.9), the iFME model for classification is given by (4.10) where Υ = (w ⊤ , Γ ⊤ 1 , . . . , Γ ⊤ K ) ⊤ is the unknown parameter vector to be estimated. We perform penalized MLE by penalizing the ML via a Lasso penalization on the derivative coefficients ω k 's and γ kg 's of the form Pen

P (y i |X i (•); Υ) = K k=1 π k (s i ; w)P (y i |v i ; Γ k ),
χ,λ (Υ) = χ K-1 k=1 ∥ω k ∥ 1 + λ K k=1 G-1
g=1 ∥γ kg ∥ 1 , with χ and λ regularization constants. The estimation is performed by using an adaptation to this classification context of the EM algorithm developed in [START_REF] Chamroukhi | Functional mixturesof-experts[END_REF]. The only difference resides in the maximization with respect to the expert network parameters Γ k .

Numerical experiments ■ Experiments on simulated data

We conducted experiments by considering a 3-class classification problem (i.e., G = 3) with a K = 2-component FME model. In particular, we generate n = 1000 functions X i (•) using

X i (t) = x ⊤ i b r (t), with x i = W v i , W ∈ R 8×8 is a matrix of i.i.d. random values from the uniform distribution U(0, 1), v i ∈ R 8 is

FME for function-on-function regression

In this section, we extend the FME model proposed in Chapter 3 to the case where the responses are functions, i.e., to relate input and output when they are both functions. In particular, we assume that in the population, the relationship between the functional output and the functional input is governed by some kernels that vary among the subpopulations of the input space. We refer to this proposed model as function-on-function FME model (FF-FME). This model can also be viewed as an extension of the mixture model for functional regression proposed in [START_REF] Devijver | High-dimensional mixture regression models, application to functional data[END_REF], where the author considered a mixture of functional regression models to deal with function-on-function regression problems. In [START_REF] Devijver | High-dimensional mixture regression models, application to functional data[END_REF], the model parameters are estimated using MLE and penalized MLE approaches, and the performance was evaluated on both simulated and real functional datasets. In our proposed FF-FME model, similar to the way the mixture of experts regression model extends the mixture model, we allow the mixing proportions to vary across predictors via a functional gating network, so that the model can potentially capture a more complex structure in the population.

This section begins with the construction of the FF-FME model. Then the model estimation via EM algorithm is presented. Finally, the numerical experiments on both simulated and real functional data are presented to illustrate the use of the proposed model.

The FF-FME model

Let {X i (u), Y i (t)} n
i=1 be a sample of n i.i.d. data pairs where Y i (t), t ∈ T , is a functional response, and X i (u), u ∈ U, is a functional predictor, where T and U are some intervals on R. An example of such context is: Y i (t) is the demand of electricity, and X i (u) is the temperature during the day at a certain house or city. In this case, both the response and the predictor are time series and U and T represent the time, e.g., from 0h to 24h.

Firstly, we formulate each expert component Expert k (y|x) in ( 2.1) by a fully functional linear model as in (2.28). In particular, the regression model for the ith observation takes the following stochastic representation 4.11) where β z i ,0 (t), defined on T , is the unknown functional intercept, β z i (t, u), defined on T × U, is the unknown kernel (i.e., unknown bivariate function) associated with expert z i to be estimated, and ε i (t) is an independent random error function follows a normal distribution with mean zero and covariance function Σ z i (t, u). Here z i ∈ [K] denotes the unknown label of the expert responsible for the ith observation. However, the appearance of β z i ,0 (t) in the model can be simplified by assuming that the functional predictors and responses are all centered (e.g., as mentioned in Section 2.7). Therefore, for simplicity, from now on we take this assumption and denote by β = {β 1 (•, * ), . . . , β K (•, * )} the set of K unknown kernels of the expert network to be estimated.

Y i (t) = β z i ,0 (t) + U X i (u)β z i (t, u)du + ε i (t), i ∈ [n], (

Modeling the functional gating network

The modeling for the functional gating network is unchanged compared to the FME model. That is, we also model the gating network by the functional logistic regression model for linear classification [START_REF] Mousavi | Multinomial functional regression with wavelets and lasso penalization[END_REF], with the note that the intercept α z,0 is now simplified. The resulting functional softmax gating network is then given by

π k (X i (u), u ∈ U; α) = exp{ U X i (u)α k (u)du} 1 + K-1 k ′ =1 exp{ U X i (u)α k′ (u)du} , (4.12)
where α k (u) is the unknown functional parameter of gating network, and

α = {α k (u), u ∈ U} K-1 k=1
is the set of unknown functional parameters α k (u), u ∈ U to be estimated. Again, the function α K (u) is necessarily zero to satisfy the identifiability condition. Analogous to the construction for the finite representation of the functional gating network as in (3.10) or in (4.5), here we can also obtain such representation for (4.12). In particular, let ω(u) = [ω 1 (u), . . . , ω q (u)] ⊤ be a vector of basis functions defined on U. Then with q sufficiently large, we can approximate α k (u) by

α k (u) = ζ ⊤ k ω(u), (4.13)
where ζ k ∈ R q is the unknown vector of coefficients of the functional gating parameter α k (u). Substituting (4.13) into (4.12) we obtain a new representation for the gating network:

π k (r i ; ξ) = exp(ζ ⊤ k r i ) 1 + K-1 k ′ =1 exp(ζ ⊤ k r i ) , (4.14)
where r i = U X i (u)ω 1 (u)du, . . . , U X i (u)ω q (u)du ⊤ ∈ R q is the new design vector, and ξ =

(ζ ⊤ 1 , . . . , ζ ⊤ K-1 ) ⊤ ∈ R q(K-1
) is now the unknown parameter vector to be estimated.

Modeling the functional expert network

It is known that in general there is no notion of probability density for functional data [START_REF] Delaigle | Defining probability density for a distribution of random functions[END_REF], so in particular in our case it is not obvious to write down a conditional density function for the functional responses Y i (t) as in (3.13). However, as can be seen in Section 2.7, using projections onto appropriate bases, and assuming a normal distribution for the coefficients of the error functions, we can also introduce the conditional density function for the FF-FME model analogously to (3.13) as follows.

Let ω ⋆ (t, u) = [ω ⋆ 1 (t, u), . . . , ω ⋆ M (t, u)] ⊤ be a vector of bivariate M basis functions defined on T × U. Then with M sufficiently large, the bivariate functional parameter β k (t, u) can be expressed 4.15) where b k = (b k1 , . . . , b kM ) ⊤ ∈ R M is the vector of unknown coefficients associated with the expert kth, to be estimated. Recall that, analogous to the equations related to basis expansions stated throughout the thesis, the equation ( 4.15) is only approximation, but we use the equality sign "=" to be able to manipulate such expansion without adding additional error terms. The stochastic representation (4.11), with the intercept has been simplified as discussed earlier and with z i = k, can be rewritten as 4.16) where X ⋆ im (t) := U X i (u)ω ⋆ m (t, u)du are functions of t which now serve as new predictors. Let ω(t) = [ω 1 (t), . . . , ω q (t)] ⊤ be a vector of basis functions defined on T . Then with q sufficiently large, the functions Y i (t), X ⋆ im (t) and ε i (t) can also be represented by

4.3 FME FOR FUNCTION-ON-FUNCTION REGRESSION by β k (t, u) = M m=1 b km ω ⋆ m (t, u) =: b ⊤ k ω ⋆ (t, u), (
Y i (t) = U X i (u) M m=1 b km ω ⋆ m (t, u) du + ε i (t) = M m=1 b km U X i (u)ω ⋆ m (t, u)du + ε i (t) = M m=1 b km X ⋆ im (t) + ε i (t), i ∈ [n], (
Y i (t) = q ℓ=1 y iℓ ω ℓ (t) =: y ⊤ i ω(t), (4.17a) X ⋆ im (t) = q ℓ=1
x ⋆ imℓ ω ℓ (t) =: x ⋆⊤ im ω(t), (4.17b)

ε i (t) = q ℓ=1 ε iℓ ω ℓ (t) =: ϵ ⊤ i ω(t), (4.17c)
respectively, where y i , x ⋆ im and ϵ i belong to R q are the coefficient vectors, and ε i ∼ N q (0, Σ) for some unknown symmetric matrix Σ ∈ R q×q .

Note that the matrix ω(t)ω ⊤ (t) is non-singular, therefore by substituting equations (4.17) into (4.16) and simplifying ω(t), we obtain a new stochastic representation for the model:

y i = M m=1 b km x ⋆ im + ϵ i = X ⋆⊤ i b k + ϵ i , i ∈ [n],
where

X ⋆ i denotes the matrix [x ⋆ i1 , . . . , x ⋆ iM ] ⊤ ∈ R M ×q .
From this representation and the Gaussian assumption for the error ε i , the conditional density of each approximated functional expert z i = k is then given by (4.18) where ϕ(•; µ, ν) is the pdf of the q-dimensional Gaussian distribution with mean µ and covariance matrix ν, and θ k = (b ⊤ k , vech ⊤ Σ k ) ⊤ is the unknown parameter vector of the expert k, k ∈ [K], to be estimated. Here, vech(•) is the operator that extracts the unique elements of a symmetric matrix [START_REF] Henderson | Vec and vech operators for matrices, with some uses in jacobians and multivariate statistics[END_REF].

f (y i |X i (•), z i = k; θ k ) = ϕ(y i ; X ⋆⊤ i b k , Σ k ),

Conditional density function of the FF-FME model

Finally, combining the gating network (4.14) and the expert density (4.18) leads to the following conditional density function for our FF-FME model 4.19) where Ψ = (ξ ⊤ , θ ⊤ 1 , . . . , θ ⊤ K ) ⊤ is the parameter vector of the model to be estimated. There are some remarks on the constructions of the design vectors in practice. First, although the formulas are written using the functions X i (u) and Y i (t), in practice they are only available in form of discretized values. However, for simplicity we use the same notations to indicate the noisy version. Second, as seen in the previous models, all the involving integrals here can be calculated numerically by a Riemann sum over the grid where data was sampled. Finally, in practice the design vectors x ⋆ im can be calculated directly from the predictors X i (t) and the basis responsible for covariate u that construct ω ⋆ (t, u), as we have seen in Section 2.7. For example, if the bivariate basis ω ⋆ (t, u) is constructed by two B-spline bases ω(t) and ω(u), where ω(u) = [ω 1 (u), . . . , ω q (u)] ⊤ , then x ⋆ im is simply

f (y i |X i (u); Ψ) = K k=1 π k (r i ; ξ)ϕ(y i ; X ⋆⊤ i b k , Σ k ), (
x ⋆ im = U X i (u)ω 1 (u)du, . . . , U X i (u)ω q (u)du ⊤ ∈ R q .

Penalized maximum likelihood estimation via EM algorithm

The FF-FME model (4.19) is now defined upon a tractable finite representation of the functional predictors and responses. In this subsection, we consider the MLE of the model via the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF].

With the maximum likelihood estimation method, the observed-data log-likelihood function for FF-FME model to be maximized is given by log

L(Ψ) = n i=1 log K k=1 π k (r i ; ξ)ϕ(y i ; X ⋆⊤ i b k , Σ k ),
where

Ψ = (ξ ⊤ , θ ⊤ 1 , . . . , θ ⊤ K ) ⊤ .
L(Ψ ) = log L(Ψ ) -Pen χ (Ψ ) -Pen λ (Ψ ).
In the above definition, Pen χ (Ψ ) is the penalty term with respect to the gating network, defined by

Pen χ (Ψ ) = χ K-1 k=1 ∥ξ k ∥ 1 ,
and Pen λ (Ψ ) is a roughness penalty term for the expert network, defined by (4.20) where D t and D u are the second-order partial derivatives with respect to t and u directions, and χ, λ t , λ u are the tuning parameters. If χ = λ t = λ u = 0, the penalty term does not contribute to L(Ψ ), we end up with a maximum likelihood estimation. Note that, the dependence of Pen λ (Ψ ) on the model parameter Ψ is understood implicitly, since the bivariate functions (D t β k )(t, u) and (D u β k )(t, u) can be expressed in terms of the coefficient vector b k as we have seen in ( 2.35) and (2.36). The two steps of the EM algorithm for FF-FME model are performed as follows.

Pen λ (Ψ ) = K k=1 λ t [(D t β k )(t, u)] 2 dtdu + λ u [(D u β k )(t, u)] 2 dtdu ,

E-step. The conditional probability memberships τ (s)

ik that the observed pair (X i (u), Y i (t)) originates from the kth expert at EM iteration sth is calculated by

τ (s) ik = P(z i = k|X i (u), Y i (t); Ψ (s) ) = π k (r i ; ξ (s) )ϕ(y i ; X ⋆⊤ i b (s) k , Σ (s) k ) f (y i |X i ( * ); Ψ (s) ) , (4.21)
where

Ψ (s) = ( ξ (s)⊤ , θ (s)⊤ 1 , . . . , θ (s)⊤ K )
⊤ is the estimate of Ψ at iteration sth of the EM algorithm.

M-step.

The maximization is performed by separately with respect to the gating network parameter ξ and the expert network parameters θ k for each of K experts. First, the maximization with respect to ξ for such a functional gating function has been done earlier in this thesis, for example in (3.16) via the MLE, or in (3.20) via the regularized MLE. Therefore, we refer to the updating gating network parameters part given in Appendix B.1 for the maximization with respect to ξ, with a minor difference that in FF-FME model we have simplified the intercept α k,0 .

The maximization with respect to θ k , on the other hand, consists of solving a weighted functionon-function regression problem, where the weights are the conditional expert memberships τ (s) ik given in (4.21). More specifically, if penalized MLE is considered, then by expressing the penalized log-likelihood function L(Ψ) in terms of the parameters of the expert kth, i.e., b k and Σ k , the maximization with respect to b k can be recognized as a penalized least square problem. Therefore, we can take advantage from the solution formula for the function-on-function regression model presented in Subsection 2.7.2.

In particular, let us denote Y = [y 1 , . . . ,

y n ] ⊤ , τ (s) k = (τ (s) k1 , . . . , τ (s)
kn ) ⊤ . The update formulas for the parameters of expert k for the next EM iteration are then given by b

(s+1) k = V ⊤ τ V τ + λ t D t + λ u D u -1 V ⊤ τ vec(I 1/2 ω Y ⊤ τ ), (4.22) Σ (s+1) k = 1 n i=1 τ (s) ik Y τ -G τ b (s+1) k ⊤ Y τ -G τ b (s+1) k . (4.23)
Follows is the definitions of the matrices appear in the above formulas. First, D t and D t are the matrices related to the penalty term, defined similarly as in (2.35) and (2.36). Particularly, if we consider D t and D u the second-order partial derivatives, which is the case in our experimental studies later, D t and D u take the following forms

D t = ∂ 2 ω g ∂t 2 (t, u). ∂ 2 ω h ∂t 2 (t, u)dtdu 1⩽g,h⩽M ∈ R M ×M , D u = ∂ 2 ω g ∂u 2 (t, u). ∂ 2 ω h ∂u 2 (t, u)dtdu 1⩽g,h⩽M ∈ R M ×M . Next, I 1/2 ω
is the principal square root of the matrix I ω given by k ⊙ Y , where the square root is componentwise and ⊙ is the Hadamard product operator. The matrix V τ is the weighted design matrix given by

I ω = T ω g (t)ω h (t)dt
V τ = τ (s) k ⊗ 1 q ⊙       I 1/2 ω x ⋆ 11 I 1/2 ω x ⋆ 12 • • • I 1/2 ω x ⋆ 1M I 1/2 ω x ⋆ 21 I 1/2 ω x ⋆ 22 • • • I 1/2 ω x ⋆ 2M . . . . . . . . . . . . I 1/2 ω x ⋆ n1 I 1/2 ω x ⋆ n2 • • • I 1/2 ω x ⋆ nM       ∈ R nq×M ,
where ⊗ is the Kronecker product operator. Finally, G τ is the weighted design matrix given by

G τ = τ (s) k ⊙       x ⋆⊤ 11 x ⋆⊤ 12 • • • x ⋆⊤ 1M x ⋆⊤ 21 x ⋆⊤ 22 • • • x ⋆⊤ 2M . . . . . . . . . . . . x ⋆⊤ n1 x ⋆⊤ n2 • • • x ⋆⊤ nM       ∈ R n×M q ,
and the matrix multiplication

G τ b (s+1) k in equation (4.23) is performed by treating each x ⋆⊤ ij , i ∈ [n], j ∈ [M ]
, as an element, the same way as discussed in (2.34).

FME FOR FUNCTION-ON-FUNCTION REGRESSION

We can see, the right hand sides of (4.22) and (4.23) depend on the current iteration of the EM algorithm implicitly via the posterior probabilities τ (s) k . Finally, if MLE is considered, then (4.22) can also be used to update the parameter b k by simply setting λ t = λ u = 0.

Numerical experiments

In this subsection, we illustrate the use of the proposed FF-FME model on both simulated and real-world data to show its practical utility. For model evaluation, we compute the MISE between the true and the predicted functional responses, which quantifies how well our model was used for predicting. To evaluate the approximation performance, we also compute the MISE between the true and the estimated functional gating networks, as well as between the true and the estimated functional expert networks. Here, the MISE for bivariate functions is simply defined as follow

MISE( β(t, u)) = E T U ( β(t, u) -β ⋆ (t, u)) 2 dtdu ,
where β ⋆ (t, u) and β(t, u) respectively represents the true and the estimated functions. Moreover, we use ARI to evaluate the clustering performance of the model.

■ Experiment on simulated data

The data generating process is as follows. Firstly, we fix two kernel functions β 1 (t, u), β 2 (t, u) and a gating function α 1 (t) those to be estimated. The domains for both t and u are fixed to be the interval [0, 1]. Next, we generate the predictors X i (u) using the formula

X i (u) = x ⊤ i ω(u), t ∈ [0, 1],
in which ω(u) is a 5-dimensional B-spline basis, and x i ∈ R 5 is the coefficient vector follows the normal distributed N 5 (0, 1). Then, the functional responses Y i (t) are constructed using the usual data generating protocol for ME models. In particular, for each i ∈

[n], sample Z i from on the multinomial distribution M 1, π 1 (X i (t); α), . . . , π K (X i (t); α) , then conditional on Z i = z i , the response Y i (t) is generated by Y i (t j ) = 1 0 X i (u)β z i (t j , u)du + ε i (t j ), i ∈ [n],
where ε i (t j ) ∼ N (0, 1) at all time points t j . Furthermore, to mimic real-world data, we contaminate the predictors X i (u) with measurement noises before performing model estimation. That means we will not use the X i (u), but

X i (u) = X i (u) + δ i (u),
for analysis. Here, δ i (u) is an independent Gaussian noise, and for all u j , δ i (u j ) follows the distribution N (0, 0.04). We generated 200 pairs (U i (u), Y i (t)) for training and another sample of 200 pairs for evaluation.

For illustration, Figure 4.2 displays some randomly taken predictors X i (u) and their responses Y i (t) that we have just generated. Results. Figure 4.3 shows the true and the estimated kernels obtained by FF-FME model. We can see, given the generated dataset, the estimated kernels approximate very well the true ones.

The results estimated without roughness penalization fit almost exactly the true kernels. MISEs for each of β 1 (t, u) and β 2 (t, u) are 0.02 and 0.01, respectively. When estimating with roughness penalization, the resulting kernels clearly less wiggle than the the ones obtained without roughness penalization. Here the simulated true kernels are smooth and the noise level in data is quite small, so the estimated ones in both cases have more chance to be smooth. However, in real-world applications, when the data are more noisy, and the underlying true kernels may not be smooth, the proposed roughness penalization is promising. The true and estimated functional gating networks are also shown in Figure 4.4 for illustration. The MISE of the predicted functional responses Y i (t) on testing set is 0.3928. Finally, the ARI between the true labels and the estimated labels using the maximum a posteriori (MAP) rule is 0.7822.

■ Application to Berkeley growth data

In this part, we apply the proposed FF-FME approach to the Berkeley growth data, which is a well-known data originally published in [START_REF] Tuddenham | Physical growth of california boys and girls from birth to eighteen years. Publications in child development[END_REF] and have broadly been analyzed in the literature since then. The data consists of the height records for 39 boys and 54 girls from age 1 to 18. In original data, the measurements were taken quarterly from ages 1 to 2, annually from 2 to 8, and semiannually from 8 till 18. Therefore, as in many related studies, we performed a simple interpolation such that the trajectories are all available quarterly from age 1 to 18. The true gating function α 1 (t) and its estimated α 1 (t).

In [START_REF] Ramsay | Comparison of height acceleration curves in the fels, zurich, and berkeley growth data[END_REF], the authors investigate the height acceleration to reveal the dynamics of human growth. In their study, they predicted the height at age 18 based on the dynamic pattern till age 9, since it is known that the growth patterns of boys and girls during their pubertal spurts differ significantly in terms of magnitude and timing. The height trajectories of the two groups and their heights at 18 can be found in the appendix (Figure A.2).

In this experiment for FF-FME model, we study from a different perspective by examining the dependence of the heights at age 9 to 18 (functional response) on the heights at age 1 to 9 (functional predictor), i.e., before the pubertal spurts. set are then fed into our algorithm. Here, we used K = 2, a tolerance of 10 -8 and five random initializations for the EM algorithm.

For simplicity, the penalty constant χ for the functional gating function is fixed to a constant. Then, we implement the FF-FME model with four different values for the roughness penalty constant to see the effect of them on the shapes of the estimated kernels. In particular, λ t = λ u = γ, with γ ∈ {0, 10 -7 , 10 -5 , 10 -2 } were used as illustrations. Figure 4.6 shows the corresponding results. As we can see, the larger roughness penalty constant, the smoother kernels were obtained. The shapes of β 1 (t, u) and β 2 (t, u) can give us some interpretations about the relationship between the heights before and after age 9. The kernels estimated without roughness penalization suggest that, for example, for boys, the height at age 18 is highly related to the height at age 6, whereas for girls it is at age 3. Further deeper interpretations could be made about the growth patterns, but that is beyond the scope of this experiment. The estimated gating function, for example when λ t = λ u = 0, is shown in Figure 4.7a.

Finally, we evaluate the prediction performance by compute the MISE between the true and the predicted responses on the testing set. Given X i (u), the predicted response Y i (t) is computed using the estimated kernels β 1 (t, u), β 2 (t, u) and the estimated gating parameter α as

Y i (t) = E Y i (t)|X i (u) = µ Y (t) + K k=1 π k (X i (u), α) U X i (u) β k (t, u)du, t ∈ T . Note that, in the notation E [Y i (t)|X i (u)]
, the independent variables t and u were written to clarify the defined domains of Y and X, although it is more correctly to write E [Y i (t)|X i ( * )].

Figure 4.7b show the predicted responses obtained by the model estimated with λ t = λ u = 0. If we compare Figure 4.7b and Figure 4.5b, we can see that the predicted responses agree with the true responses about the average heights of two groups, i.e., higher for boys and lower for girls. Moreover, the roots of the MISEs for training and testing data are 1.64 and 1.54, respectively. Here, the root of MISE is reported because it is directly proportional to the unit of the height (cm). This means that, on testing set, given the heights at ages 1 to 9, the average error of predictions for heights at ages 9 to 18 is 1.54 cm.

FME for function-on-scalar regression

To complete the family of FME models, this sections extends the FME model proposed in Chapter 3 to the case where the responses are functions and the predictors are scalars. In particular, we consider a function-on-scalar FME (FS-FME) in which the responses are univariate functions, e.g., curves or time series, and we wish to relate them with the explanatory variables via a mixture of experts model. The dedicated algorithm is developed, via maximizing log-likelihood function, to estimate the expert and gating networks. Numerical experiments are also presented to illustrate the usage of our proposed model. 2.1) by the function-on-scalar linear model as in (2.20). In particular, the regression model for the ith observation takes the following stochastic representation

The FS-FME model

Let {x i , Y i (t)} n i=1 be a sample of n i.i.d. data pairs where x i ∈ R p is the input and Y i (t) is a functional response with t ∈ T ∈ R. For each k ∈ [K], we formulate the expert component Expert k (y|x) in (
Y i (t) = β z i ,0 (t) + p j=1 x ip β z i ,j (t) + ε i (t), i ∈ [n], (4.24)
where ε i (t) is an independent random error function follows a normal distribution with mean zero and covariance function Σ z i (t, s), and z i ∈ [K] denotes the unknown label of the expert responsible for the ith observation.

For simplicity, we can simplify the functional intercept β z i ,0 by assuming that the number one was already incorporated into the predictor vector x i . Then conditional on z i = k the above CHAPTER 4 EXTENSIONS OF FUNCTIONAL MIXTURES OF EXPERTS equations can be rewritten under matrix form simply as

Y i (t) = x ⊤ i β k (t) + ε i (t), i ∈ [n], (4.25)
where β k (t) = [β k1 , . . . , β kp ] ⊤ is the vector of functional parameters to be estimated.

Modeling the gating network

The gating function in this case can take the same form as one of those in the vectorial setting, e.g., softmax gating function as in (2.3) or Gaussian gating function as in (2.4). Here, we opt to model the gating network for FS-FME model using softmax modeling

π k (x i ; α) = exp(x ⊤ i α k ) 1 + K-1 k ′ =1 exp(x ⊤ i α k ′ ) , (4.26)
where α = (α 1 , . . . , α K-1 ) ⊤ ∈ R p(K-1) is the gating parameter vector to be estimated.

Modeling the expert network

Analogously to the modeling of expert network in the FF-FME model, we assume that the responses Y i (t), the parameters β kj (t), and the error functions ε i (t) can be represented by a collection of basis functions ψ(t) = [ψ 1 (t), . . . , ψ q (t)] ⊤ as

Y i (t) = y ⊤ i ψ(t) ⊤ , β kj (t) = b ⊤ kj ψ(t), and ε i (t) = ϵ ⊤ i ψ(t),
where y i , b kj , ϵ i ∈ R q are the corresponding coefficient vectors, and ϵ i ∼ N (0, Σ z i ). Then, the equations in (4.25) implies

y i = x ⊤ i B k + ϵ ⊤ i , i ∈ [n], (4.27)
where

B k = [b k1 , . . . , b kp ] ⊤ ∈ R
p×q is now the unknown coefficient matrix associated with the expert kth to be estimated. From (4.27), the model for a coefficient vector y i of a functional response Y i (t) given x i can be expressed as a probability density function as follows

f (y i |x i , z i = k; θ k ) = ϕ(y i ; x ⊤ i B k , Σ k ), (4.28)
where ϕ(•; µ, ν) is the pdf of the q-dimensional Gaussian distribution with mean µ and covariance matrix ν, and θ k = (vec(B ⊤ k ) ⊤ , vech ⊤ Σ k ) ⊤ is the unknown parameter vector of the expert k, k ∈ [K], to be estimated.

FME FOR FUNCTION-ON-SCALAR REGRESSION

Conditional density function of the FS-FME model

Gathering the gating function (4.26) and the expert density (4.28) we obtain the following conditional density function for the FS-FME model:

f (y i |x i ; Ψ) = K k=1 π k (x i ; α)ϕ(y i ; x ⊤ i B k , Σ k ), (4.29)
where Ψ = (α ⊤ , θ ⊤ 1 , . . . , θ ⊤ K ) ⊤ is the parameter vector of the model.

Penalized maximum likelihood estimation via EM algorithm

The maximum likelihood estimation for FS-FME model can be obtained via the EM algorithm in the same manner as in the FF-FME model considered in previous section. In addition, to obtain smooth fits for the functional parameters β kj (t) we impose a roughness penalty term to the log-likelihood function and consider the following penalized log-likelihood function

L(Ψ ) = log L(Ψ ) -Pen λ (Ψ ), (4.30)
where log L(Ψ ) is the observed-data log-likelihood of Ψ and Pen λ (Ψ ) is a roughness penalty term that encouraging smoothness for the functional experts defined by

Pen λ (Ψ ) = K k=1 p j=1 λ kj T (Dβ kj )(t) 2 dt, (4.31)
where D is some differential operator, for example, the second derivative. Here, similarly to the FF-FME model, the dependence of the penalty term of Ψ is implicitly. The two steps of the EM algorithm for FS-FME model are performed as follows.

E-step. We calculate the conditional probability memberships τ (s) ik given by τ (s)

ik = P(z i = k|x i , y i ; Ψ (s) ) = π k (x i ; α (s) )ϕ(y i ; x ⊤ i B (s) k , Σ (s) k ) f (y i |x i ; Ψ (s) ) , (4.32)
where

Ψ (s) = ( ξ (s)⊤ , θ (s)⊤ 1 , . . . , θ (s)
⊤ K ) ⊤ is the estimate of Ψ at the iteration sth of the EM algorithm.

M-step. In the M-step, the gating parameter and the expert parameters will be estimated separately as we have seen many times before. Firstly, for the gating network, similarly to the FF-FME model, the maximization with respect to α can be done as in (3.16) via the MLE, or in (3.20) via the regularized MLE. Therefore, we again refer to the updating gating network parameters part given in Appendix B.1 for the maximization with respect to α.

For each functional expert k, because of the stochastic representation (4.24), we can recognize that the maximization with respect to θ k consists of solving a weighted function-on-scalar regression problem, where the weights are the current conditional expert memberships τ (s) ik given in (4.32). Therefore, we can take advantage from the solution formula for the function-on-scalar regression model presented in Subsection 2.6.2. In particular, let us denote

Y = [y 1 , . . . , y n ] ⊤ , X = [x 1 , . . . , x n ] ⊤ , τ (s) k = (τ (s) k1 , . . . , τ (s)
kn ) ⊤ , the update formula for the parameter vector of expert kth for the next EM iteration is given by vec

( B ⊤ k ) (s+1) = U ⊤ τ U τ + Λ k ⊗ D -1 U ⊤ τ vec(I 1/2 ω Y ⊤ τ ), U τ = X τ ⊗ I 1/2 ω , (4.33)
Σ (s+1) k = 1 n i=1 τ (s) ik Y τ -X τ B (s+1) k ⊤ Y τ -X τ B (s+1) k , (4.34)
where

Λ k = diag(λ k1 , . . . , λ kp ) ∈ R p×p , D = T (Dω g )(t)(Dω h )(t)dt g,h ∈ R q×q , I 1/2 ω I 1/2 ω = T ω g (t)ω h (t)dt g,h ∈ R q×q , X τ = τ (s) k ⊙ X ∈ R n×p , Y τ = τ (s) k ⊙ Y ∈ R n×q
where the square root is understood componentwise, ⊗ is the Kronecker product operator and ⊙ is the Hadamard product operator. Similarly to the update formulas in FF-FME model, the dependence of the right hand side of (4.33) on the sth iteration of the EM algorithm is implicit via τ (s) k .

Numerical experiments

Data generating process

In this part, we illustrate the usage of the proposed FS-FME model on simulated data, to relate functional responses to scalar predictors. The data generating process is follows. First, we fix a true functional expert network that consists of the functions β kj (t) as can be seen in Figure 4.8.

Here K = 3, p = 3, and we intentionally created rough curves so that we can see the effect of the smooth penalization on the estimation. Next, we generated 200 predictors x i with whose components follow N (0, 1), N (1, 1), N (2, 1) with equal proportions. The functional responses Y i (t) are then generated by

Y i (t)|Z i = z i , x i = x ⊤ i β z i (t) + ε i (t), Z i |x i ∼ M 1, π 1 (x i ; α), . . . , π K (x i ; α) .
We can see, the wiggles in Y i come from both the expert functions and the error functions.

Implementation and results

We split the dataset into training and testing sets with equal proportions, estimate the functional experts using the training set, and report the evaluation metrics on the testing set. Specially, we implement the FS-FME model with different roughness penalty constants to illustrate the effect of the smooth penalization on the estimation. In this study, we set the operator D to the second derivative.

To evaluate the estimation of the functional experts, we compute the following quantity

MISE( β λ ) = 1 Kp K k=1 p j=1 T ( β λ kj (t) -β kj (t)) 2 dt,
where β kj (t) is the true parameter function, and β λ kj (t) is the parameter function estimated with penalty constant λ. To evaluate the prediction, we compute the MISE between the true and the predicted responses in the testing set.

Figure 4.9 shows the estimated expert functions β kj (t) with three different penalty constants: λ 1 = 0, λ 2 = 0.001, λ 3 = 1, and the corresponding predicted responses Y (t) for the predictors in the testing set. As we can see, the larger λ results in the smoother for the estimated functions. With λ = 1, the estimated experts tends to straight lines. The corresponding values for MISE( β λ ) are 0.0151, 0.0107, and 0.1607. The prediction on testing set can be seen visually in Figure 4.9. The MISEs in prediction on the testing set are 24.6612, 24.4827, and 25.4218, respectively for the three estimated models.

Summary

In this chapter, we extended the FME, FME-Lasso and iFME models proposed in Chapter 3 to the case of categorical responses, i.e., for multiclass classification problems. The corresponding EM algorithm were also developed for the extended models. Simulation study was constructed to illustrate the performance of the models. We also applied the models to the phoneme data and obtained promising results in terms of both classification correct rate and accuracy of parameter fitting.

The first time ME models are developed to deal with regression problems with functional responses. In particular, we proposed FS-ME and FF-ME models to model the relationship between functional responses and vetorial predictors (FS-ME model)/ functional predictors (FF-ME model) in heterogeneous populations. The dedicated algorithms were presented for parameter estimations. These models have many potential applications as illustrated in the experimental studies subsections, especially the effect of roughness penalization on producing smooth fits for the underlying true kernels.

There are many other extensions can be considered. For example, we can develop interpretable variants for the FF-FME and FS-FME models. In particular, the functional gating network in the FF-FME model can be defined upon the finite representation of the gating function's derivatives. Like that, the estimated kernels can be truly sparse and get interpretable shapes as seen in the iFME models.

INTRODUCTION

stations, or by privacy restriction such as with data in finance and medicine. Even if a single full dataset is available, it often requires a huge processing resource for handling and running statistical inference methods, this may be evenly more expensive.

Such restrictions may prevent one from obtaining global estimators for the models assumed on the data. Parallelizing or sequentializing the existing learning algorithms and statistical methods is therefore an active research area. Most of works are based on a two-step divide and conquer procedure: • Local inference (or local computation) and • Aggregation. In the first step, standard inference or computation is performed on subsets of data available at local machines, then in the second step, the local results are collected to a central machine where they are aggregated.

Although relying on the same divide and conquer principle, these attempts can be categorized into two groups: one parallelizes the model computation and one parallelizes the model inference. Parallelizing the model computation means the estimator for the model parameter is maintained and updated on a central machine, while the computations for that update procedure are distributed and carried out on local machines. However, not all the computations of the existing algorithms can be parallelized effectively, by the nature of the formulas, or by the high cost of communication across the machines during the computations. Therefore, parallelizing the model inference is a preferred approach. It performs standard inference on local machines to obtain local estimators, then transmits them to a central machine where they are aggregated together, based on some efficient strategies, to produce an overall estimator for the model parameter.

There are many successful attempts in this direction of parallelizing the existing learning algorithms and statistical methods. Those may be mentioned here include parallelizing stochastic gradient descent [START_REF] Zinkevich | Parallelized stochastic gradient descent[END_REF], parallelizing multiple linear regression [START_REF] Mingxian | Parallelizing multiple linear regression for speed and redundancy: an empirical study[END_REF], parallel K-Means clustering based on MapReduce [START_REF] Zhao | Parallel k-means clustering based on mapreduce[END_REF], distributed learning for heterogeneous data via model integration [START_REF] Merugu | A distributed learning framework for heterogeneous data sources[END_REF], split-and-conquer approach for penalized regressions (Chen and ge Xie, 2014), for logistic regression [START_REF] Shofiyah | Split and conquer method in penalized logistic regression with lasso (application on credit scoring data)[END_REF], distributed learning for finite Gaussian mixtures [START_REF] Zhang | Distributed learning of finite gaussian mixtures[END_REF], among others.

To make this chapter self-complete, let us briefly recall the notion of the ME framework. The main idea of ME models themselves is also the divide-and-conquer principle, in which a complex problem is divided into smaller simple sub-problems, then each one can be solved by a special expert. A ME model can be written in general notation as

f (y|x) = K k=1 Gate k (x) Expert k (y|x),
in which f (y|x), the distribution of response y given the covariate x, is modeled as a mixture distribution with covariate-dependent mixing proportions Gate k (x) and conditional mixture components Expert k (y|x). In ME model terminology, Gate k (x) is referred to as gating function, and Expert k (y|x) is referred to as expert function, while K ∈ N is the number of experts. ME is therefore a fully conditional mixture model that allows the mixing proportions to be functions of the covariates.

The most popular models for the function Gate k (x) are the softmax gating function and the Gaussian gating function, while for the conditional mixture components Expert k (y|x), there are many choices depend on the considered problem. When the response y is continuous and assumed to be Gaussian given the input x, Expert k (y|x) can be modeled by the conditional density function of Gaussian distribution. When the response y is binary, resp. categorical, the suitable choice for Expert k (y|x) is the conditional density function of the binary, resp. multinomial, logistic regression distribution. In addition to these popular ones, many other models for the conditional mixture components have been studied [START_REF] Chamroukhi | Non-normal mixtures of experts[END_REF]. In situations where standard mixture of regression models fails to approximate the conditional probability distribution of data, ME model may be a good replacement as it can capture a wide variety of data distributions thanks to their universal approximation property well established for almost of popular gating and experts functions [START_REF] Nguyen | Approximation results regarding the multiple-output Gaussian gated mixture of linear experts model[END_REF](Nguyen et al., , 2021a)). Some successful applications of ME architecture may be mentioned here include ME for time series prediction [START_REF] Zeevi | Time series prediction using mixtures of experts[END_REF][START_REF] Yümlü | Financial time series prediction using mixture of experts[END_REF], segmentation [START_REF] Chamroukhi | Joint segmentation of multivariate time series with hidden process regression for human activity recognition[END_REF][START_REF] Chamroukhi | Time series modeling by a regression approach based on a latent process[END_REF], ME for social network data [START_REF] Gormley | A mixture of experts latent position cluster model for social network data[END_REF], for classification of gender and pose of human faces [START_REF] Gutta | Mixture of experts for classification of gender, ethnic origin, and pose of human faces[END_REF], etc. For an overview of practical and theoretical aspects of ME modeling, reader is referred to [START_REF] Nguyen | Practical and theoretical aspects of mixture-of-experts modeling: An overview[END_REF]; [START_REF] Yuksel | Twenty years of mixture of experts[END_REF].

In this chapter, we particularly introduce a distributed learning approach for the softmax-gated mixtures of Gaussian experts model because of its popularity in applications. Let (x, y) ∈ X × Y be a random observation from some probability model, where x ∈ X ⊂ R d (d ∈ N) is the input and y ∈ Y ⊂ R is the output. A softmax-gated mixture of Gaussian experts model links y and x via the following conditional density function

f (y|x; θ) = K k=1 π k (x; α)φ(y; x ⊤ β k , σ 2 k ), (5.1)
where θ = (α ⊤ , β ⊤ 1 , . . . , β ⊤ K , σ 2 1 , . . . , σ 2 K ) ⊤ is the parameter vector of the model. Here, π k (x; α) is the softmax gating function given by

π k (x; α) = exp(x ⊤ α k ) 1 + K-1 k ′ =1 exp(x ⊤ α k ′ )
,

where α = (α ⊤ 1 , . . . , α ⊤ K-1 ) ⊤ , α k ∈ R d for k ∈ [K -1],
is the parameter vector of the gating network; φ(•; µ, σ 2 ) is the conditional Gaussian density function with mean µ and variance σ 2 , and

(β ⊤ k , σ 2 k ) ⊤ ∈ R d+1
is the parameter vector of the expert k for k ∈ [K]. Such setting has been exclusively studied in multivariate analysis. However, in practice when the dataset (based on it we estimate the model parameter) is large, the computation for parameter estimation may be costly or evenly infeasible. Therefore, a distributed learning approach for this ME model is particularly interesting and will be developed in this chapter.

The rest of this chapter is organized as follows. In Section 5.2 we formulate the problem of aggregating the local estimates of a ME model and we will propose a numerical strategy for solving it in Section 5.3. Section 5.5 is dedicated for numerical study and Section 5.6 is for summary and 5.2 AGGREGATING DISTRIBUTED ME MODELS discussion.

Aggregating distributed ME models

Problem setting

Let D = {(x i , y i )} N i=1 be a random sample of N (with N ∈ N supposed to be large) independently identically distributed (i.i.d.) data pairs from a ME model f (•|x, θ * ) as in (5.1), where (x i , y i ) ∈ X × Y, X ⊂ R d , Y ⊂ R q , and θ * is the true parameter vector. For simpler notation, we will denote by f * the true conditional density function f (•|x, θ * ). Note that x i , y i denote observed data, whereas x, y denote the generic variables.

With the assumption that N is so large that it is difficult to estimate θ * on a single computer within a reasonable amount of time, a natural approach is to use the divide-and-conquer principle to estimate θ * . Suppose that D is randomly divided into M disjoint subsets D 1 , . . . , D M and stored on M local machines. Let f (•|x, θ m ), abbreviated simply by f m , be the conditional density function (with parameter vector θ m ) of the ME model estimated based on sample D m , namely,

f m := f (•|x, θ m ) = K k=1 π k (x; α (m) ) φ(•|x; β (m) k ), (5.2)
where θ m = ( α (m) , β

(m) 1 , . . . , β (m)
K ) obtained by maximizing the observed-data log-likelihood function (2.9) with the data in D m .

The question is how to approximate the true global density f * from the local densities f m , and more importantly, how to aggregate the local estimators θ m to produce a single aggregated estimator for the true parameter vector θ * . Here and after, we use "local density" and "local estimator" to refer to the conditional density function f m and the parameter vector θ m , respectively, of the ME model estimated at local machine m based on D m , m ∈ [M ]. Besides, for the reader's convenience, from now on the local models and local parameters will be denoted with hat notations.

Aggregation strategies

Let N m be the sample size of the subset D m , i.e., M m=1 N m = N . A natural strategy to approximate the true density f * is to use the weighted average density f W defined by

f W = M m=1 λ m f m , (5.3)
where λ m denote the sample proportions Nm N those sum to one, and f m are the local densities defined in (5.2). This average density f W fits our intuition about an average model, and in fact, approximates very well the true density f * . However, there are two issues with this aggregating strategy. First, it only gives us an expression to approximate the conditional probability values, i.e., for each x ∈ X we have f W (y|x) ≈ f * (y|x), it does not give us the approximation expression for the true parameter θ * . In other words, it is not clear how to obtain a "weighted average" approximation for θ * from the local estimators θ m 's. Second, since each f m is a mixture of K components, we can express the density f W in (5.3) as

f W = M m=1 λ m K k=1 π k (x; α (m) ) φ(•|x; β (m) k ) = M m=1 K k=1 λ m π k (x; α (m) ) φ(•|x; β (m) k ).
(5.4)

One can see that M

m π k (x; α (m) ), for k ∈ [K], m ∈ [M ]
. Therefore, although this M K-component mixture approximates well the true density f * , it has a wrong number of components, i.e., M K instead of K, this will make the clustering results useless or difficult to interpret. Hence, this aggregating approach is not what we are aiming for.

Let M K denotes the space of all K-component ME models as in (5.1). Two other common strategies to aggregate the local densities f 1 , . . . , f M are via

f B := f (y|x; θB ) = arg inf g∈M K M m=1 λ m ρ( f m , g), (5.5
)

and f R := f (y|x; θR ) = arg inf g∈M K ρ( f W , g), (5.6)
where ρ(•, •) is some divergence defined on the space of finite mixture distributions. The solution f B and f R are often known as barycenter solution and reduction solution, respectively, with their interpretations are given as follows.

In the case of f B , we are finding a K-component ME model g that can be viewed as a barycenter the local models f 1 , . . . , f M with respect to the weights λ 1 , . . . , λ M and the divergence ρ(•, •). Whereas, in the case of f R , we are finding a K-component ME model g that is closest with the weighted average density f W defined in (5.3), with respect to the divergence ρ(•, •). Because we are searching for solutions in M K , both f B and f R solve the problem of wrong number of components in f W .

Both solutions are desirable, and in fact, can be shown to be connected under specific choices of ρ(•, •). Indeed, for any bivariate function ρ(•, •) that is linear in the first argument, f B and f R are identical. However, in this chapter we prefer the reduction solution f R for many reasons. First, 5.2 AGGREGATING DISTRIBUTED ME MODELS as we have already remarked, the M K-component mixture f W is a good approximation to the true density f * , so it makes more sense to find a K-component ME model that approximates f W , which we already know is good. Second, the barycenter approach may lead to a counter-intuitive solution under some specific ρ(•, •), e.g., as already shown in [START_REF] Zhang | Distributed learning of finite gaussian mixtures[END_REF] for the case of univariate Gaussian mixture models and the 2-Wasserstein divergence with Euclidean ground distance. Finally, as we can see, given the same divergence ρ(•, •), the computation in (5.5) will be more expensive than that of (5.6), especially in our large data context. Thus, our objective in the next sections is to develop an efficient algorithm to find the reduction estimator θR . We shall denote the components of the reduction estimator by θR = ( ᾱR , βR 1 , . . . , βR K ). A description of our considered problem can be seen in Figure 5.1, in which, the gating and expert parameters are explicitly indicated for each component of the local and the reduced models.

λ 1 , π 1 (x; α (1) ), ϕ(y|x; β (1) 1 ) . . . λ 1 , π K (x; α (1) ), ϕ(y|x; β (1) K ) . . . λ M , π 1 (x; α (M ) ), ϕ(y|x; β (M ) 1 ) . . . λ M , π K (x; α (M ) ), ϕ(y|x; β (M ) K ) machine 1, D 1 machine M , D M π 1 (x; ᾱR ), ϕ(y|x; βR k ) . . . π K (x; ᾱR ), ϕ(y|x; βR K ) central machine . . . M K components K components

Some remarks and notation conventions

Recall that, as we can see in equation (5.4), the weighted average mixture f W has totally M K terms of the form λ m π k (x;

α (m) )φ(•|x; β (m) k ), for m ∈ [M ], k ∈ [K]
. However, due to the presence of the weights λ m and the fact that the local gating parameters α (1) , . . . , α (M ) are not yet integrated into a single gating parameter vector, f W do not has the standard form of a ME model as in (5.1). For mathematical convenience, we will incorporate the weights λ m into the local gating functions π k (x; α (m) ) and write them simply by { π ℓ (x)} ℓ∈[M K] , which will be also referred to as gating functions. The corresponding experts in f W will be abbreviated by { φ ℓ (y|x)} ℓ∈[M K] . Similarly, for g ∈ M K , we will abbreviate its K gating functions by {π k (x)} k∈[K] , and its K experts by

{φ k (y|x)} k∈[K] . We then have f W = M K ℓ=1 π ℓ (x) φ ℓ (•|x) and g = K k=1 π k (x)φ k (•|x),
(5.7)

for g ∈ M K . Here, the dependence of the functions on the parameters is understood implicitly.

U The interval [0, 1] ⊂ R f W , f B , f R
Weighted average, barycenter and reduction densities

{π k (x)} k∈[K] Gating functions π k (x; α), k ∈ [K], with unknown α {φ k (•|x)} k∈[K] Expert functions φ(•|x; β k ), k ∈ [K], with unknown β k { π ℓ (x)} ℓ∈[M K] Gating functions π k (x; α (m) ), m ∈ [M ], k ∈ [K] { φ ℓ (•|x)} ℓ∈[M K] Expert functions φ(•|x; β (m) k ), m ∈ [M ], k ∈ [K] T c (h, g)
Transportation divergence between mixtures h and g

T c (g) Abbreviation of T c ( f W , g) R c (g)
Relaxation of T c (g) from the constraints involving π S c (g, g (t) )

Majorant function of R c (g) at g (t)

Table 5.1: Summary of notation

Optimal transport approach for the reduction strategy

Let us substitute (5.3) into (5.6) and write the problem more explicitly as

f R = arg inf g∈M K ρ M m=1 λ m f m , g , (5.8)
that is, we are seeking for a K-component ME model g of form as in (5.1) that is closest to the M K-component mixture f W = M m=1 λ m f m with respect to some divergence ρ(•, •). In general, an analytical solution is difficult to obtain for such reduction problem, especially in the context of ME models where the conditional density function is made up of several different gating and expert functions. Therefore, the aim of this section is to propose a framework to find a numerical solution

OPTIMAL TRANSPORT APPROACH FOR THE REDUCTION STRATEGY

for the the reduction model f R .

We observe that, firstly, the solution depends heavily on the choice of the divergence ρ(•, •), which measures the goodness of a candidate model g ∈ M K . A good choice of ρ should take into account the dissimilarity between the experts, as well as the weights incorporated in the gating functions between the mixtures. Then, the best candidate mixture g should minimize the divergence from the large mixture f W . Secondly, compute the divergence between two mixtures is difficult, while it would be easier to compute the divergence between two component densities, which may has explicit formula in some specific cases. Therefore, borrowing the idea in Zhang and Chen (2021) that finds a reduction estimator for finite Gaussian mixtures based on minimizing transportation divergence between the large starting mixture and the desired mixture, we establish a framework for solving a reduction estimator for ME models, which is generally formalized as in (5.8). Particularly, in the case φ(y|x; β) being Gaussian experts, if conditionally on one single observation of x, problem (5.8) can be viewed as the problem of finding reduction estimator for finite mixtures of 1D Gaussian distributions.

Expected transportation divergence

Let h = L ℓ=1 π ℓ (x) φ ℓ (y|x), L ∈ N, and g = K k=1 π k (x)φ k (y|x) be two mixture models of L and K components, respectively. Here, h will play the role as the mixture f W in problem (5.6). We wish to define a divergence that measures the dissimilarity between the two mixtures h and g. Let Φ denotes the family of the component conditional densities φ(•|x; β) =: φ(•|x). Let π and π be, respectively, the column vectors of gating functions π ℓ (•)'s and π k (•)'s. For x ∈ X , we denote

Π x ( π, π) := Π ( π(x), π(x)) = P ∈ U L×K : P 1 K = π(x), P ⊤ 1 L = π(x) ,
where U denotes the interval [0, 1], 1 K and 1 K are vectors of all ones. In other words, for x ∈ X , Π x ( π, π) denotes the set of all matrices P of size L × K, with entries P ℓk ∈ [0, 1], satisfying the marginal constraints

K k=1 P ℓk = π ℓ (x), for all ℓ ∈ [L],
(5.9a)

L ℓ=1 P ℓk = π k (x), for all k ∈ [K].
(5.9b) Given x ∈ X , one can think of such matrix P ∈ Π x ( π, π) as a plan that transports an amount of material distributed according to π(x) to distributed according to π(x), for short, we will say "transports π(x) to π(x)". Suppose that the transportation is costing for each unit of material. We are then interested in the optimal way, i.e., with a minimum cost, to transports π(x) to π(x). The optimal plan is clearly depends on x. We will write P (x) to indicate such optimal plan.

P 1 =        0.2 0 0 0 0.1 0 0 0.1 0 0 0.1 0.1 0 0 0.3 0 0 0.1        P 2 =        0.1 0.1 0 0 0 0.1 0 0.1 0 0 0.1 0.1 0.1 0 0.2 0 0 0.1        Figure 5
.2: Example of transportation plans with L = 6 and K = 3. Given x ∈ X , the vectors π(x) and π(x) can be viewed as the distributions of the weights (sum to one). Here, for example, the weights are distributed as in the left and right panels. P 1 and P 2 are two valid plans those transport π(x) in two different ways to π(x).

Definition 5.3.1 (Expected transportation divergence). Given the above notations, the expected transportation divergence between two mixtures h and g is defined by 5.10) where c(•, •) is a real-valued bivariate function defined on Φ × Φ satisfies c(φ 1 , φ 2 ) ⩾ 0 for all φ 1 , φ 2 ∈ Φ, and the equality holds if and only if

T c (h, g) = E inf P ∈Πx( π,π) L ℓ=1 K k=1 P ℓk c ( φ ℓ (•|x), φ k (•|x)) =: E T c (h, g, x) , (
φ 1 = φ 2 .
Here, the expectation is taken with respect to x and the function c is often referred to as the cost function. Assume that, for a given x ∈ X , P * (x) is an optimal solution to the optimization problem inside the expectation operator in (5.10). Then the value of T c (h, g, x) can be interpreted as the optimal total cost of transporting π(x) to π(x) with the unit cost of transportation is proportional to the values c ( φ ℓ (•|x), φ k (•|x))'s. This interpretation is common in optimal transportation literature under the name Kantorovich formulation, e.g., [START_REF] Oberman | An efficient linear programming method for optimal transportation[END_REF]. The value T c (h, g) is therefore defined as the expected value of these optimal transportation costs.

The cost function c must be chosen such that it not only reflects the dissimilarity between the conditional densities φ ℓ (•|x) and φ k (•|x), but also has to be easy to calculate. Later we will see that c being KL-divergence is a suitable choice for both the cases when the experts are Gaussian regression models (i.e., for regression problems), and when the experts are logistic regression models (i.e., for classification problems). 5.11) From now on, we will use T c (g) to refer to T c ( f W , g), since our objective now is to minimize the transportation divergence T c ( f W , g) as a function of g ∈ M K . Problem ( 5.11) appears to involve two optimizations: one over Π x ( π, π) and one over M K . However, we will show that the constraint P ∈ Π x ( π, π) in fact can be relaxed to P ∈ Π x ( π, •). This means we only need P to satisfy (5.9a), the constraints (5.9b) are redundant. In other word, instead of searching, for each x, a plan P that satisfies π(x), we can move π, the gating functions make up g, to match P .

f R = arg inf g∈M K T c ( f W , g) = arg inf g∈M K E inf P ∈Πx( π,π) M K ℓ=1 K k=1 P ℓk c ( φ ℓ (•|x), φ k (•|x)) • (
Indeed, for f W and g as in (5.7), let us define R c ( f W , g) as a function of g as follow

R c ( f W , g) = E inf P ∈Πx( π,•) M K ℓ=1 K k=1 P ℓk c ( φ ℓ (•|x), φ k (•|x)) =: E R c ( f W , g, x) .
(5.12)

Here, the M K × K matrix P is now free of the constraints involving the gating network π of g, namely, for all x ∈ X , the equality P ⊤ 1 M K = π(x) need not to be hold. Then, the following proposition will claim that solving f R = arg inf Before presenting the proposition, we define P( f W , g, x) a function of g ∈ M K and x ∈ X as follow

P( f W , g, x) = arg inf P ∈Πx( π,•) M K ℓ=1 K k=1 [P ℓk c ( φ ℓ (•|x), φ k (•|x))] .
( 5.13) This function returns the optimal plan for the optimization problem inside the expectation operator in the definition of R c ( f W , g). For simplicity of notation, similarly to T c (g), from now on we let the dependence of the functions on f W in the background, i.e., R c ( f W , g), R c ( f W , g, x) and P( f W , g, x) will be written by R c (g), R c (g, x) and P(g, x), respectively.

Proposition 5.3.1. Let T c (g), R c (g) and P(g, x) defined as above. Then

inf g∈M K T c (g) = inf g∈M K R c (g).
The reduction solution is hence given by f R = arg inf g∈M K R c (g), (5.14) and the gating functions of f R can be calculated by (5.15) where P ℓk ( f R , x) denotes the entry (ℓ, k) of P( f R , x).

π k (x) = M K ℓ=1 P ℓk ( f R , x), ∀x ∈ X ,
The proof of Proposition 5.3.1 is given in Appendix C.1.1. It is important to note here that the equation ( 5.15) only gives us an expression for computing the values of the gating network for each x ∈ X , it does not yet have the parametric form as in the standard model (5.1). Thus, to obtain the desired gating parameter ᾱR as described in Figure 5.1, we need to perform an additional estimation step as in equation ( 5.21), which will become clearer later in Subsection 5.4.2. Thus, thanks to Proposition 5.3.1, our objective now is to minimize R c (g) with respect to g. Replacing T c (g) by R c (g), problem (5.11) 5.16) This optimization can be done with the help of majorization-minimization (MM) algorithm (e.g., [START_REF] Lange | MM optimization algorithms[END_REF]). We defer the numerical approach for solving the problem (5.16) in Section 5.4, now we will show that the problem is well-posed and state the conditions for the consistency of the reduction estimator.

becomes f R = arg inf g∈M K R c (g) = arg inf g∈M K E inf P ∈Πx( π,•) M K ℓ=1 K k=1 P ℓk c ( φ ℓ (•|x), φ k (•|x)) • (

Well-posedness and consistency of the reduction estimator

We will make the following standard assumptions. ■ Well-posedness. First, we observe that for any x ∈ X , the minimization problem inside the expectation operator in (5.16) is in fact a linear programming (LP) problem. Moreover, for x ∈ X , we see that Π x ( π, •) is a nonempty set, and the sum M K ℓ=1 K k=1 P ℓk c ( φℓ (•|x), φ k (•|x)) is bounded below by zero for all P ∈ Π x ( π, •). Therefore, this LP problem has a global minimizer (e.g., see Boyd and Vandenberghe (2004)). Hence, for all g ∈ M K , there exists non-negative finite expectation for the random quantity R c (g, x).

Now, for a mixture g ∈ M K as in (5.7), we let C(P , φ) be a function of the transport plan with respect to the component densities of g, given by

C(P , φ) = M K ℓ=1 K k=1 P ℓk c ( φ ℓ (•|x), φ k (•|x)) , where P ∈ U M K×K , U = [0, 1] ⊂ R, and φ = (φ 1 (•|x), . . . , φ K (•|x)) ∈ Φ K .
Here, we can think of the family Φ as the space of parameter vectors that make up the density φ(•|x) (typically R κ , for some κ ∈ N * ). Then, the assumptions A2 and A3 on the continuity and convexity of the cost function c immediately lead to the continuity and convexity in φ of the function C(P , φ). Moreover, C(P , φ) as an affine function is obviously continuous and convex in P .

Next, for x ∈ X , let I(g, x) be a function of g and x defined by

I(g, x) = inf P ∈Πx( π,•)
C(P , φ),

i.e., the optimal cost given g and x inside the expectation operator in (5.16). We see that, the value of I(g, x) depends on the mixture g only through the component densities φ, but the gating functions. Then, the objective function of the problem (5.16) can be written as R c (g) = E [I(g, x)].

We have the following proposition.

Proposition 5.3.2. Given the notations as above. Assume in addition that there exists ∆ ∈ R + such that I(g, x) ⩽ ∆ for all g ∈ M K , x ∈ X . Then R c (g) is continuous and convex as a function of g ∈ M K . It follows that the problem (5.16) has a global solution.

The proof of the Proposition 5.3.2 is given in Appendix C.1.2. Note that the condition on the boundedness of the optimal transportation costs is common, and the solution to the problem (5.16) is not necessarily unique.

■ Consistency. The reduction estimator θR has a desired property that it is a consistent estimator of the true parameter θ * as soon as the local estimators are consistent estimators of θ * . We have the following proposition. t) be the model obtained at MM iteration t-th. Let 5.17) where P ℓk (g (t) , x) is given by

S c (g, g (t) ) = E M K ℓ=1 K k=1 P ℓk (g (t) , x) c( φ ℓ (•|x), φ k (•|x)) , (
P ℓk (g (t) , x) =    π ℓ (x) if k = arg inf k ′ ∈[K] c( φ ℓ (•|x), φ (t)
k ′ (•|x)) 0 otherwise.

(5.18)

Then S c (g, g (t) ) is a majorant function of R c (g) at g (t) .

The proof of Proposition 5.4.1 is in Appendix C.1.4. Here we see that, the plan with entries P ℓk (g (t) , x) given in (5.18) is an optimal plan for the optimization problem (5.13) when g = g (t) (this will be clear in the proof of the proposition). Note that, similarly to R c (g), the function S c (g, g (t) ) varies only through the parameters of the conditional densities φ k (•|x)'s.

Then, the next step is to minimize S c (g, g (t) ), with respect to g, to obtain a next point g (t+1) for the MM algorithm. As we can see, the optimization for S c (g, g (t) ) can be performed separately. At each iteration tth, the parameter vector of the expert k can be optimized via .19) According to each specification of the experts and the cost function c(•, •), problem (5.19) has a certain form that is supposed to be easy to solve, as we will see later for the cases of the Gaussian and logistic regression experts, with c is the KL-divergence. Hence, our algorithm will alternate between two following steps until convergence:

φ (t+1) k (•|x) = arg inf φ∈Φ E M K ℓ=1 P ℓk (g (t) , x) c( φ ℓ (•|x), φ(•|x)) . ( 5
• calculating majorant function S c (g, g (t) ) as in (5.17), which in fact only involves updating P ℓk (g (t) , x) according to (5.18);

• updating the experts parameters by solving (5.19).

At the convergence, the parameters of all the experts are optimized, so it remains to update the parameter of the gating functions. To do that, we rely on the equation (5.15), in which the theoretical solution f R is replaced by the model g * obtained at the convergence of the MM algorithm. This procedure of updating the gating parameter will be presented in Subsection 5.4.2.

As we can see, up to this point all the optimizations carried on the central machine are involving the expectations with respect to x. In practice, the distribution of x is unknown. Therefore, our approach is to create a small sample D S = {(x s , y s )} S s=1 of size S by randomly drawing from the full dataset D. This D S acts as a supporting sample for the calculations on the central machine. Then, the expectation E µ(x) [ • ] will be understood as the empirical expectation calculated based on the empirical distribution µ(•) = S -1 S s=1 δ xs (•). The size of D S is also an issue for further consideration. However, let us assume that D S has a sufficiently large size to characterize the distribution of the variable x. We will denote by X S the set of all observations x s in the supporting sample D S , and where there is no ambiguity we also use X S to denote the matrix with each row corresponds to one row vector x ⊤ s . The Algorithm 5.1 gives a pseudo-code that summarizes our main algorithm of aggregating ME models.

Algorithm 5.1 Aggregating ME local estimators

Input:

-local gating parameters α (m) , m ∈ [M ]; -local expert parameters β (m) k , m ∈ [M ], k ∈ [K]; -sample proportions λ m , m ∈ [M ].
Output: Reduction estimator θR = ( ᾱR , βR 1 , . . . , βR K ).

1: Initialize a ME model g (0) , i.e., initialize θ (0) = (α (0) , β

1 , . . . , β

K ). Assign t ←-0 2: Sampling a supporting sample D S from D 3: repeat 4: for k ∈ [K] do 5: for ℓ ∈ [M K] do 6: (0)
For all x ∈ D S , compute P ℓk (g (t) , x) according to (5.18) 7:

end for 8:

Update expert parameters by solving (5.19). Depending on weather the expert is Gaussian or logistic regression model, update formula takes the form (5.23) or (5.23), respectively 9: end for 10:

Assign t ←t + 1 11: until the objective function (5.17) converges. 12: Estimate the gating parameter according to (5.21).

Updating the gating network parameter

As mentioned earlier, although the equation (5.15) gives a formula to compute the mixing weights for every covariate x, to obtain the final desired ME model as described in Figure 5.1, it is necessary to find the parameter ᾱR for the gating functions π k (•, ᾱR ). Therefore, we propose to estimate ᾱR via solving a softmax regression problem in which the predictors are X S and the responses are computed according to (5.15). In particular, let A = [a sk ] be a matrix in R S×K defined by a sk = M K ℓ=1 P ℓk (g * , x s ), (5.20) where g * is the model obtained at convergence of the MM algorithm, x s ∈ D S , and P ℓk (g * , x s ) is computed according to (5.18). Then ᾱR can be estimated via the maximum likelihood approach ᾱR = arg max .21) This problem can be efficiently solved by the Newton-Raphson procedure, for example as used for updating the gating network in the M-Step of the EM algorithm in Section 2.3.1. Here, we note that the parameter vector is also initialized, i.e., the Kth gate's parameter is constrained to be zero, this is necessary for proving the consistency of θR . Now, we are left with the question of how to obtain the update for the expert parameter β (t+1) k

α S s=1 K k=1 a sk log π k (x s ; α) = arg max α S s=1 K-1 k=1 a sk α ⊤ k x s -log 1 + K-1 k ′ =1 exp{α ⊤ k ′ x s } . ( 5
as the solution of (5.19). In the next subsections, we will present the update formulas for the two cases of the Gaussian regression expert and logistic regression expert when the cost function c(•, •) is KL-divergence.

Updating the Gaussian regression experts parameter

Now, we consider the case of φ(•|x) being Gaussian regression experts, i.e., φ(•|x) = φ(•|x ⊤ β; σ 2 ), where β ∈ R d and σ 2 ∈ R + are the parameters of the expert. The KL-divergence between two Gaussian regression experts takes the following form

KL φ 1 (•|x)∥φ 2 (•|x) = 1 2 log σ 2 2 σ 2 1 + σ 2 1 σ 2 2 + (β 2 -β 1 ) ⊤ xx ⊤ (β 2 -β 1 ) σ 2 2 -1 .
Therefore, with the experts being Gaussian and the cost function c(•, •) being KL(•∥•), the objective function of the problem (5.19) can be written as a function of β and σ 2 as .22) where P (t) ℓk (x) denotes P ℓk (g (t) , x), which is calculated by (5.18) (5.19) can be simply rewritten by (β

K(β, σ 2 ) := E M K ℓ=1 P (t) ℓk (x) 1 2 log σ 2 σ 2 ℓ + σ 2 ℓ σ 2 + (β -β ℓ ) ⊤ xx ⊤ (β -β ℓ ) σ 2 -1 , ( 5
with c( φ ℓ (•|x), φ (t) k ′ (•|x)) is now replaced by KL( φ ℓ (•|x)∥φ (t) k ′ (•|x)). Then,
(t+1) k , σ 2 k (t+1) ) = arg inf (β,σ 2 )∈R d ×R + K(β, σ 2 ).
The following proposition gives a solution to this problem, i.e., the updating formula for the parameters of the Gaussian regression experts. The detailed calculations can be found in Appendix C.1.6.

β (t+1) k = (X ⊤ S X S ) -1 X ⊤ S (log(V (t) k ) -log(1 -V (t) k )), V k = D (t) k -1 W (t)⊤ k U ⊤ , ( 5

Experimental study

Models in comparison

In this section, we illustrate the effectiveness of our proposed approach on simulated and real-world datasets of various sizes. In particular, we compare the performances of our reduction estimator θR and the following estimators:

• Global estimator, denoted by θ G . It is the MLE of the ME model estimated based on the full dataset on the central machine.

• Middle estimator, denoted by θMid . It is the parameter of the local density (among M local densities) that gives the smallest sum of transportation divergences, weighted by λ m , with the other local densities, i.e.,

f M id = f (y|x; θMid ) = arg inf g∈{ f 1 ,..., f M } M m=1 λ m T c ( f m , g).
• Weighted average estimator, denoted by θW . It is an ad-hoc estimator defined as the pointwise weighted average of the local estimators, i.e.,

θW = M m=1 λ m θ m .
One can see this estimator is sensitive with label switching and non-i.i.d. data partitions, i.e., when there are differences between the partitions in terms of sample sizes of in terms of the equilibrium between clusters. However, when the partitions are i.i.d., have equal sizes, In this simulation study, we opted to consider the regression problem. To have more fairly data generating process, these true parameter vectors and the cluster centers are also selected randomly from the set of integers in the interval [-5, 5]. For each of the N 's considered above, we also generated 100 corresponding testing sets such that the ratio of training-testing split is 80% -20%, more precisely, the testing sets will have 50k, 125k and 500k observations, respectively. It is worth mentioning that each of the datasets with one million observations consumes roundly 1.6 GB of memory. The estimators are computed using training sets and the metrics (except the learning time) are evaluated on testing sets.

Implementation details and results

For the distributed approaches, i.e., the reduction, middle and weighted average estimators, we consider four settings of M , respectively, 4, 16, 64 and 128 machines. The data are distributed equally to the local machines, and the size of the supporting set is taken to be equal to the size of local data, i.e., S = N m = N/M for all m ∈ [M ].

The Global estimator is obtained by running the EM algorithm, on a single computer. The learning time for each of them is therefore the total time of running the corresponding algorithm. On the other hand, for the reduction estimator, since the local inferences can be performed parallelly on multiple machines, the learning time is recorded as the sum of the maximum local time and the aggregating time, i.e., the time consumed by the Algorithm 5.1. Similarly, the learning time of Middle and Weighted average estimator is the sum of the maximum local time and the time of the corresponding calculations to obtain them. All EM algorithms involving the global and local inferences are run with five random initializations. Finally, for each of the obtained estimators, we compute the transportation distance form the true model, the log-likelihood, MSE, RPE and ARI on the corresponding testing sets. model evaluation metrics described in Subsection 5.5.2.

Figure 5.3 shows the box plots of 100 Monte Carlo runs of the evaluation metrics when N = 100000. For each panel, e.g., the transport distance, the x-axis has two rows. The lower one represents the number of machines used for training, while the upper one represents the leaning approaches in comparison. As we can see, the prediction and clustering performance of our reduction estimator is as good as that of the global estimator when M is 4 or 16 machines. When M is 64 or 128 machines, the errors of the reduction estimator are slightly higher than the global estimator, but still better than the middle and weighted average estimators. In terms of transportation distance, we can see that with the reduction estimator, the obtained ME model is as close as the ME model formed by the global estimator, even when M = 128 machines. Finally, in terms of learning time, the distributed approach requires much less time than the global approach, for example 25 times on average when 4 machines are used.

Figure 5.4 and Figure 5.5 show the results when N is 500k and one million, respectively. Here, the evaluation metrics behave similarly to those in Figure 5.3, except in the MSE panels where the reduction estimator appears to be much worse than the middle and weighted average estimators when M = 128. However, given the magnitudes of the MSE in these cases, the difference is actually not too large. This behavior can also be explained by the fact that for N of 500k and one million, when using the distributed approaches, especially for datasets with one million observations. The prediction and clustering performance metrics such as MSE, RPE and ARI behave as expected across estimators and across sample sizes.

EXPERIMENTAL STUDY

Model evaluation, implementation details and results

Beside the RPE metric as used before in many regression problems throughout the thesis, we use the root mean square error (RMSE) value to evaluate the HR prediction since it is directly proportional to the unit of the predicted values. In addition, to compare with the related work in [START_REF] Oyeleye | A predictive analysis of heart rates BIBLIOGRAPHY using machine learning techniques[END_REF], the scatter index (SI), which is simply the RMSE divided by the average value of the observed value, was also computed. It is commonly understood that SI < 10% is a good model, and SI < 5% is a very good model. Now, we describe how we integrate our distributed leaning approach into the FME model proposed in Chapter 3, it is quite simple as follows. First, the functional predictors, the functional expert and the gating networks are represented using some appropriate basses (here we used Bspline). Then the conditional density function of the model can be written in a vectorial fashion as in (3.13). Note that, if we use a same basis for the functional expert and gating networks (which is our case here), the design vectors r i and x i are identical. Therefore, the model density (3.13) is now taking the same formula as the density (5.1), which is our starting point for developing the distributed learning approach. This means we can employ the Algorithm 5.1 to estimate the parameter vector Ψ in (3.13) of the FME model in a distributed manner.

The RPE, RMSE, SI and the learning time of each approach are given in Table 5.2. We can see, the distributed approach significantly reduced the learning time while other prediction metrics are not too different. The estimated functional expert and gating functions are shown in Figure 5

Summary

In this chapter, a distributed learning approach for ME models was proposed. In particular, the learning approach is based on minimizing the expected transportation divergence between two ME models: one is the large ME model that consists of all local mixtures obtained from the local machines, and one is the targeted K-component ME model. An algorithm based on MM principle was also proposed to solve the resulting minimization problem. The proposed distributed learning approach was successfully applied to simulated and real-world datasets. This distributed approach was shown to have a leaning time much less than the global approach, but still give results as good as those based on training with full dataset. Specially, we applied the proposed distributed learning approach to the FME model for functional data proposed in Chapter 3 with promising remarks: the integration is easy and the prediction performance is as good as the global approach, while the learning time has been much improved.

Conclusion

This thesis was written regarding to the ME models for functional data and a distributed learning approach for ME models. The goal of developing efficient ME models for functional data was achieved through the introduction of a new class of ME models, referred to by functional mixturesof-experts (FME), whose members including i. FME for regression, i.e., scalar-on-function FME, ii. FME for classification, also a scalar-on-function FME, iii. FF-FME for regression, i.e., function-on-function FME, iv. FS-FME for regression, i.e., function-on-scalar FME.

In addition, for scalar-on-function FME models, we have developed their interpretable variants: iFME for regression and iFME for classification, which were able to produce highly-interpretable fits for the functional expert and gating networks. The goal of developing a distributed learning approach for ME models was also achieved through the proposed reduction strategy and an efficient developed MM algorithm whose performance was demonstrated via application to simulated and real-world data.

More specifically, in Chapter 3, we proposed a ME for regression problems that takes functional data as inputs and continuous scalars as outputs. We refer to it by FME. Model estimation was performed using MLE and regularized MLE, which results in the FME-Lasso model that imposes a certain degree of sparsity on the estimated parameters. Furthermore, to obtain truly sparse fits for the functional expert and gating networks, we proposed an interpretable version of the FME CHAPTER 6 CONCLUSION AND FUTURE DIRECTIONS model, referred to by iFME model, to produce fits that are not only sparse itself, but also sparse with respect to their derivatives. The proposed FME, FME-Lasso, and iFME models for regression were applied to simulated and real-world datasets. The implementation of the algorithms, model selection, as well as the interpretation of the results were also discussed in detail. Matlab codes for the proposed FME model (as well as FME-Lasso and iFME models) will be available at the following link: https://github.com/fchamroukhi/FunME.

In Chapter 4, we extended the FME, FME-Lasso, and iFME models to the classification problems. In this case, we modeled the experts by the functional multinomial regression models. The resulting models were also outperform the current state of the art methods as the case of regression in the previous chapter. The corresponding EM algorithms were also developed for model estimations in the classification case. Furthermore, the FME model was extended to the case of functional responses via the introduction of two new models: FF-FME and FS-FME. Specially, the proposed models were estimated using penalized MLE that penalizes the roughness (via appropriate derivative operators) and results in smooth fits for the functional parameters. The corresponding EM-like algorithms were also developed for model estimations. Numerical experiments were shown to illustrate the application of the models.

Finally, in Chapter 5, a distributed learning approach was proposed for ME models to deal with the situation where the data is so large that it is difficult to estimate the model parameter on a single computer, within a reasonable amount of time. The approach is based on the divideand-conquer principle. In particular, the local estimators (of ME models on local machines) are aggregated by minimizing an expected transportation divergence between two ME models. The MM algorithm was also developed particularly for Gaussian and logistic regression experts. The experiments on simulated and real data verified the efficiency of our approach.

Future directions

Given the promising results of our proposals and the motivation from a methodological point of view, there are many further interesting extensions that can be considered in future work. First, with the FME models for regression and classification, it is natural to extend them to the multiple case, i.e., for each expert/gate k there are more than one functional parameter. This can also be applied to extend the iFME models.

Second, an interpretable approach to the FF-FME and FS-FME models can be explored. In particular, as we can see in Chapter 4, the models were able to generate smooth fits for the functional parameters via roughness penalties with second-order (partial) derivatives, they do not ensure highly sparse fits as, for example, in Figure 3.7. Therefore, according to our idea, a representation upon the coefficients of the derivatives of the functional parameter can be established for the model. In this way, constraining sparsity on the resulting parameter vectors can lead to truly sparse fits 6.2 FUTURE DIRECTIONS for the experts and gating networks.

Regarding the distributed learning approach for ME models, some interesting considerations can also be made. For example, in high-dimensional setting, the local parameter vectors are oftenly estimated with regularization, so the update formula such as the one in (5.23) must be adapted to preserve (at some level) the sparsity in the local estimators. This can be approached using a majority voting method as in Chen and ge Xie (2014). Moreover, the proposed distributed learning approach can be applied to the iFME, FF-FME, and FS-FME models.

Although numerous extensive experiments were constructed throughout the thesis to demonstrate the performance of the proposed models, real-world data experiment for the FS-FME model unfortunately was not included. However, we would like to mention that, given the performance on the simulated data, we strongly believe that FS-FME model can also produce good results when applied to real-world data. Numerical experiment for the distributed learning approach in classification problems was also not provided, it therefore is one of our future works. Finally, although the proposed models and the proposed learning approach work quite well, the theoretical results, however, need further investigation. is the square root of the non-negative definite matrix I ψ defined by

T ∥R ω(t) -XB ω(t)∥ 2 dt = T n i=1 (r ⊤ i ω(t) -x ⊤ i Bω(t)) 2 dt = n i=1 T (r ⊤ i ω(t) -x ⊤ i Bω(t)) 2 dt = n i=1 T (r ⊤ i -x ⊤ i B)ω(t) [ω(t)] ⊤ (r ⊤ i -x ⊤ i B) ⊤ dt = n i=1 (r ⊤ i -x ⊤ i B) T ω(t) [ω(t)] ⊤ dt (r ⊤ i -x ⊤ i B) ⊤ . (A.
I ψ = T ψ k (t)ψ ℓ (t)dt 1⩽k,ℓ⩽L ∈ R L×L .
Then (A.8) can be reformulated by

∥RI 1/2 ω -GbI 1/2 ω ∥ 2 = ∥ vec(I 1/2 ω R ⊤ ) -V b∥ 2 , (A.9)
which can be verified simply by the definitions of the matrices. Thus, the problem (A.7) can thus be rewritten by

minimize b∈R K ∥ vec(I 1/2 ω R ⊤ ) -V b∥ 2 + λ t b ⊤ D t b + λ u b ⊤ D u b . (A.10)
By taking the first derivative of the objective function and setting it to zero, the solution is obtained as given in (2.37).

A.2 Descriptions of data

Tecator data

Tecator is a well-known data in FDA. This dataset consists of 215 spectrometric curves and the corresponding Water, Fat and Protein content of meat samples. Each spectrometric curve corresponds to the absorbance measured at 100 wavelengths (100 discretization points from 850mm to 1050mm). The curves are split according to [START_REF] Ferraty | Nonparametric functional data analysis : theory and practice[END_REF] into two classes: with small (less than 20) and large fat content obtained by an analytical chemical processing.

Berkeley growth data

Berkeley growth data was originally published in [START_REF] Tuddenham | Physical growth of california boys and girls from birth to eighteen years. Publications in child development[END_REF] and was analyzed by many studies of functional data since then. The data consists of the height records for 39 boys and 54 girls. In original data, the measurements were taken quarterly from ages 1 to 2, annually from 2 to 8, and semiannually from 8 till 18. Therefore, it is common to perform a simple interpolation such that the trajectories are available quarterly from age 1 to 18.

ik (α k,0 + ζ ⊤ k r i ) -log 1 + K-1 k ′ =1 exp{α k ′ ,0 + ζ ⊤ k ′ r i } • (B.3)
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 (B.3) according to the procedure (3.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:

H kℓ (ξ; Ψ (s) ) = ∂ 2 Q(ξ; Ψ (s) ) ∂ξ k ∂ξ ⊤ ℓ = - n i=1 π k (r i ; ξ (t) ) δ kℓ -π l (r i ; ξ (t) ) r i r ⊤ i ,
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 g(ξ; Ψ (s) ) = ∂Q(ξ; Ψ (s) ) ∂ξ = g 1 (ξ (t) ), . . . , g K-1 (ξ (t) ) ⊤ , where, for k ∈ [K -1], g k (ξ (t) ) = ∂Q(ξ;Ψ (s) )

∂ξ k = n i=1 τ (s)
ikπ k (ξ (t) ; r i ) r ⊤ i .

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

Q(θ k ; Ψ (s) ) = n i=1 τ (s) ik log ϕ(y i ; β k,0 + η ⊤ k x i , σ 2 k ) = - 1 2σ 2 k n i=1 τ (s) ik y i -(β k,0 + η ⊤ k x i ) 2 - n 2 log(2πσ 2 k )•
Thus, updating θ k = (β k,0 , η ⊤ k , σ 2 k ) ⊤ , consists of a weighted Gaussian regression problem where the weights are the conditional memberships τ (s) ik , and the updates are given by (3.17).

B.2 EM-LASSO FOR ℓ 1 -REGULARIZED MLE OF THE FME MODEL a weighted regularized multinomial logistic problem for which Q χ (ξ; Ψ (s) ) is its penalized loglikelihood. 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

l k (α k,0 , ζ k ) = - 1 2 n i=1 w ik (c ik -α k,0 -r ⊤ i ζ k ) 2 + C(ξ (t) ),
where w ik = π k (ξ (t) ; r i ) 1π k (ξ (t) ; r i ) , (weights)

c ik = α (t) k,0 + r ⊤ i ζ (t) k + τ (s)
ikπ k (ξ (t) ; r i ) w ik , (working response)

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

(α k,0 ,ζ k ) l k (α k,0 , ζ k ) -χ∥ζ k ∥ 1 , χ > 0, (B.7)
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 (B.7) 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. is the weighted working response vector, X w = [ √ w|S w |0 n×q ] ∈ R 2q+1 is the weighted design matrix, Ω = diag(0, 1 ⊤ q , ϱ1 ⊤ q ) is the diagonal weighting matrix and

A = [0 q |A [d 2 ] q A [d 1 ] q -1
| -I q ] is the constraints matrix. Here, √ w = ( √ w 1 , . . . , √ w n ) ⊤ , S w = [ √ w 1 s 1 , . . . , √ w n s n ] ⊤ , with s i are the design vectors (see (3.27)), 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 ω∈R 2q+1

-∥Ω ω∥ 1 subject to |X ⊤ w (c w -X w ω)| ⩽ χ1 2q+1 , A ω = 0 q , where the absolute value operator is understood componentwise. By decomposing ω into its positive and negative parts, ω = ω +ω -, the above problem becomes max This section presents how to solve problem (3.37). Firstly, we fix σ 2 k to its previous estimate and perform an update for (β k,0 , γ Similarly to the problem in the gating network update, by decomposing γ into its positive and negative parts, γ = γ +γ -, the above problem becomes max Finally, the update for σ 2 k is

( ω + , ω -)∈R 4q+2 -[0, 1 ⊤ q , ϱ1 ⊤ q , 0, 1 ⊤ q , ϱ1 ⊤ q ] ω + ω - subject to                X ⊤ w X w -X ⊤ w X w -X ⊤ w X w X ⊤ w X w ω + ω - ⩽ χ + X ⊤ w c w χ -X ⊤ w c w , A -A ω + ω - = 0
( γ + , γ -)∈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
σ 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 .
C.1 PROOFS hold will lead to a contradiction to (C.2). Indeed, suppose that there exists k ∈ [K] such that the expert parameter βR k does not tend to β * k at some event. Then due to the property of the cost function c(•, •), the value T ( f R , f * ) must be bounded from below by some positive constant, which obviously contradicts the fact (C.2).

We are left to show that ᾱR also converge to α * in probability. The estimation of ᾱR via MLE in (5.21) makes itself a consistent estimator for the parameter of the logistic regression model with predictors x s ∈ X S and responses a sk given in (5.20). Since α (m) p → α * and the fact that the gating function is continuous with respect to the parameter, for all x s ∈ X S we have π k (x s ; α (m) ) Moreover, because at the convergence of the MM algorithm, the obtained experts parameters are supposed to be the theoretical solution βR k , the matrix P ℓk (g * , x) in ( 5.20) can be written by

P ℓk (g * , x s ) =    π ℓ (x s ) if k = arg inf k ′ ∈[K] c( φ ℓ (•|x s ), φ R k ′ (•|x s )) 0 otherwise, (C.4)
where φ R k ′ (•|x s ) denotes φ k ′ (•|x s ; βR k ), the kth expert density of the reduced ME model. Therefore, by noting the convention about the connection between λ m π k (x s ; α (m) ) and π ℓ (x s ) (discussed at the beginning of the Subsection 5.2.3), we can see that the responses a sk are in fact the approximations of π k (x s ; α * ). This means ᾱR is also a consistent estimation of α * .

C.1.4 Proof of Proposition 5.4.1

We will use the definition to prove that S c (g, g (t) ) given in (5.17) is a majorant function of R c (g) given in (5.12) at g (t) . That is, we will show that S c (g, g (t) ) ⩾ R c (g) and the equality holds when g = g (t) .

First, for all g ∈ M K and x ∈ X , from the definition of P ℓk (g (t) , x) we see that This means the transportation plan P = [P ℓk (g (t) , x)] ℓk is a member of Π x ( π, •). Therefore, by the definition of S c (g, g (t) ) and R c (g) we have S c (g, g (t) ) ⩾ R c (g).

C APPENDIX C

Finally, at g = g (t) , by definition the functions S c (g, g (t) ) and R c (g) are respectively given by R c (g (t) ) = E inf Obviously, the transportation P = [P ℓk (g (t) , x)] ℓk with P ℓk (g (t) , x) given in (5.18) is a minimizer of the minimization in (C.5). This means that the two quantities above are equal.

C.1.5 Proof of Proposition 5.4.2

To see this, it is necessary to minimize the function K(β, σ 2 ) given in (5.22). Firstly, we observe that K(β, σ 2 ) is continuously differentiable and convex with respect to β, so any stationary point is a global minimizer of K(β, σ 2 ) as a function of β. On the other hand, K(β, σ 2 ) is coercive with respect to σ 2 , then it has a global minimizer that satisfies the first-order optimality condition. The partial derivatives of K(β, σ 2 ) with respect to β and σ 2 are given by ∂K ∂β = E where U = [0, 1] ⊂ R. Then the above equation can be written under matrix form as

X ⊤ S D (t) k X S β = X ⊤ S X S ⊙ (W ⊤ k B ⊤ )1 p ,
where ⊙ denotes the Hadamard product. This follows the root for β is given by

β = (X ⊤ S D (t) k X S ) -1 X ⊤ S X S ⊙ (W (t)⊤ k B ⊤ )1 p .
Similarly, setting ∂K ∂σ 2 to 0 one gets

k )σ 2 = Σ ⊤ W (t) k 1 S + 1 ⊤ M K W (t) k ⊙ (β -B) ⊤ X ⊤ S X S (β -B)1 M K ,
where Σ = ( σ 2 1 , . . . , σ 2 M K ) ⊤ ∈ R M K + . The conclusion is followed immediately.

C.1.6 Proof of Proposition 5.4.3 To see the update formula for the parameter vector of the logistic regression expert, it is necessary to minimize the function K Bi (β) given in (5.25). The first order derivative of K Bi (β) is given by

∂K ∂β = E M K ℓ=1 P (t) ℓk (x)
x exp(x ⊤ β) 1 + exp(x ⊤ β) -x exp(x ⊤ β ℓ ) 1 + exp(x ⊤ β ℓ ) .

Replacing the expectation operator by the empirical expectation and setting ∂K ∂β to 0 we obtain the following equation

U =       exp(x ⊤ 1 β 1 ) 1+exp(x ⊤ 1 β 1 ) exp(x ⊤ 1 β 2 ) 1+exp(x ⊤ 1 β 2 )
. . . . . .

exp(x ⊤ S β M K ) 1+exp(x ⊤ S β M K )       ∈ U S×M K . (C.11)
Then the equation (C.9) can be rewritten under matrix form as

X ⊤ S D (t) k E1 S = X ⊤ S W (t)⊤ k U ⊤ 1 S .
We see that, the above equation holds if

E = D (t) k -1 W (t)⊤ k U ⊤ , or equivalently, X S β = log(V (t) k ) - log(1 -V (t) k ) with V k := D (t) k -1 W (t)⊤ k U ⊤
. Therefore, we obtain β = (X ⊤ S X S ) -1 X ⊤ S (log(V

(t) k ) - log(1 -V (t) k )).
Modeling and Learning with Mixtures of Experts for Functional Data and Distributed Data Abstract: Mixture of experts (ME) models are popular in statistics and machine learning and have been studied for high-dimensional vectorial data and that are centralized. However, in many problems, we observe time series data and applying the existing ME models directly to these raw data may limit the performance because they ignore the correlation between variables, an intrinsic nature of functional data, and therefore do not adequately capture the inherent functional structure of the data. In many applications the data may also be not available at a centralized mode, and there is therefore a need to develop adapted distributed strategies allowing for efficient parallel computations. This thesis studies mixture of experts (ME) models in functional data and large-scale problems, in a heterogeneous scenario. The main objectives are i) to deal with situations in which we are given functional predictors (e.g., time series) to predict a potentially functional response, and ii) to learn from distributed data. We first propose in this thesis a new family of ME models, called functional mixture of experts (FME), which includes scalar-on-function FME, function-on-scalar FME, and functionon-function FME. We introduce regularized model estimation approaches via appropriate regularizations that encourage sparse and smooth estimates, while being interpretable. We develop efficient EM algorithms to maximize the corresponding observed-data log-likelihood functions and conduct extensive experimental studies to highlight the performance of the proposed models and algorithms. Then, to scale-up the ME estimation to a potentially distributed data, we develop a distributed learning strategy of ME models. It performs standard inference on local machines to obtain local estimators, then transmits them to a central machine where they are aggregated. Based on minimizing a proposed expected transportation divergence, the local estimators are aggregated to obtain a reduced estimator that is consistent with the global one, i.e., the estimator that could be possibly constructed upon the full dataset. Experimental studies demonstrate the performance of our approach.

Keywords

  Contents 1.1 Context of study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Contributions and outline of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

Figure 2 . 3 :

 23 Figure 2.3: Examples of (a) B-splines, (b) Fourier, (c) Gaussian and (d) sigmoidal basis functions

Figure 2 . 4 :

 24 Figure 2.4: The first four estimated eigenfunctions and their associated eigenvalues of the Tecator data.

Figure 2 . 5 :

 25 Figure 2.5: The first four estimated eigenfunctions and their associated eigenvalues of the Berkeley growth data.

  and θ = {b, Σ} is the unknown parameter of the model to be CHAPTER 2 STATE OF THE ART estimated. The maximum likelihood estimators of b and Σ are given by b

Figure 2 . 6 :

 26 Figure 2.6: Examples of clustering Tecator data using non-parametric technique. Here, the curves are clustered into two groups by K-means using (a) d 0 measure, (b) d 1 measure, and (c) d 2 measure. Bottom panels: (d) Tecator data, (e) the first derivative, and (f) the second derivative curves.

3 . 4

 34 Experimental studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Evaluation criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2 Simulation studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.3 Application to real-world data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52 3 . 3

 33 INTERPRETABLE FUNCTIONAL MIXTURE OF EXPERTS (IFME)

3. 4

 4 EXPERIMENTAL STUDIES 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 non-zero free coefficients in the vectors ζ, η, ω, and γ. Note that, because of the constraints (3.27) for the iFME model, free coefficients in ω and γ consist of only the part related to the d 1 derivative. That is, df

  .40) 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

Figure 3 . 1 :

 31 Figure 3.1: The true expert and gating functions used in simulations.

Figure 3 . 2 :

 32 Figure 3.2: 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}. 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(3.39). The evaluation

  Figure 3.3: Scatter plots of y i against U i (t 1 ) (top panels) and U i (t 150 ) (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 . 4 :

 34 Figure 3.4: 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.

3. 4 . 3 . 1 Figure 3 . 5 :

 43135 Figure 3.5: (a) 35 daily mean temperature measurement curves; (b) Geographical visualization of the stations, in which the sizes of the bubbles corresponds to their log of precipitation values and the colors correspond to their climate regions.

CHAPTER 3 FUNCTIONALFigure 3 . 6 :Figure 3 . 7 :

 33637 Figure 3.6: 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 3 . 8 :

 38 Figure 3.8: 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 3 .

 3 Figure 3.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.

Figure 3 .

 3 Figure 3.11: Values of modified BIC for (a) FME, (b) FME-Lasso and (c) iFME, versus the number of experts K, fitted on DTI data. The square points correspond to highest values.

  .1) where β k represents the set of coefficient functions {β kg (t), t ∈ T } and intercepts {β k,0 } for k ∈ [K] = {1, . . . , K} and g ∈ [G], and y ig = I {y i =g} .

[d 2

 2 ] k provide approximations to the d 1 and the d 2 derivatives of the coefficient function α k (•), respectively, which we denote as α (d 1 ) k (•) and α (d 2 ) k (•). Therefore, enforcing sparsity in ω k will constrain α (d 1 ) k (•) and α (d 2 )

  a vector of i.i.d. random values from the normal distribution N (1, 10), and b r (t) is the 8-dimensional B-spline basis defined on [0, 1]. Then, Y i is generated conditional on X i (•) and the true coefficient functions. The time domain is divided into a grid of functions for phonemes (log-odd to "ao")

Figure 4 . 1 :

 41 Figure 4.1: Results of (top) FME-EM-Lasso and (bottom) iFME-EM on phoneme data.

∈

  R q×q . The matrix Y τ denotes the weighted version of Y , i.e., Y τ = τ (s)

ResponsesFigure 4 . 2 :

 42 Figure 4.2: 50 randomly taken predictors and their responses.

4. 3

 3 FME FOR FUNCTION-ON-FUNCTION REGRESSION

Figure 4 . 3 :

 43 Figure 4.3: Panels (a) and (b) are the true kernels β 1 (t, u) and β 2 (t, u). Panels (c) and (d) are the estimated kernels β 1 (t, u) and β 2 (t, u), estimated without roughness penalization. Panels (e) and (f) are the estimated kernels β 1 (t, u) and β 2 (t, u), estimated with roughness penalty constant λ t = λ u = 0.001.

Figure 4 .

 4 5 displays our functional predictors X i (u) and functional responses Y i (t), where u ∈ [1, 9] and t ∈ [9, 18].

Figure 4 . 5 :

 45 Figure 4.5: The functional predictors and the true functional responses from Berkeley data.

Figure 4

 4 Figure 4.7: (a) Estimated functional gating function resulted by FF-FME model and (b) the functional responses predicted by FF-FME model, colored in their predicted cluster.

Figure 4 .8 4 . 4 Figure 4 . 8 :

 44448 Figure 4.8: The true functional experts and the noisy functional responses in the simulated data.

Figure 4 . 9 :

 49 Figure 4.9: The estimated functional experts resulted by different roughness penalty constants, and the predicted responses for the predictors in the testing set.

  π k (x; α (m) ) = 1 for all x ∈ X , so f W can be viewed as a mixture model with M K components in which the component densities are φ(•|x; β (m) k ) and the gating functions are λ

Figure 5 . 1 :

 51 Figure 5.1: Description of the considered aggregation problem. Left panels are the components of the local ME models with their corresponding estimated parameters. Right panels are the components of the aggregated K-component ME model, i.e., the desired solution. The unknowns are therefore ᾱR , βR 1 , . . . , βR K .

Figure 5 .

 5 2 gives an example of transportation plans. CHAPTER 5 DISTRIBUTED LEARNING FOR MIXTURES OF EXPERTS π

5. 3

 3 OPTIMAL TRANSPORT APPROACH FOR THE REDUCTION STRATEGY Thus, by considering ρ(•, •) being the expected transportation divergence T c (•, •), problem (5.6) becomes

  g∈M K T c ( f W , g) can actually be reduced to solving f R = arg inf g∈M K R c ( f W , g), which is much easier.

A1

  The dataset D = {(x i , y i )} N i=1 is an i.i.d. sample from the K-component ME model f (y|x, θ * ), in which the parameters are ordered and initialized. A2 The cost function c(•, •) is continuous in both arguments, and c(φ 1 , φ 2 ) → 0 if and only if φ 1 → φ 2 in distribution. 5.3 OPTIMAL TRANSPORT APPROACH FOR THE REDUCTION STRATEGY A3 The cost function c(•, •) is convex in the second arguments. This property holds for most of popular divergences such as Kullback Liebler, Hellinger distance, ..., see e.g., Dragomir (2013).

5. 4

 4 AN MM ALGORITHM FOR CONSTRUCTING THE REDUCTION ESTIMATOR Proposition 5.4.1. Let g (

5. 5

 5 EXPERIMENTAL STUDY Proposition 5.4.3. Let the notations be as defined above. If c(•, •) is chosen to be KL (•∥•), ) in (5.19) are logistic regression experts, then the update formulas for their parameters are given by

  those in Proposition 5.4.2, and U is given in (C.11).

5. 5

 5 EXPERIMENTAL STUDY problem, or θ * := (α * ⊤ , β * 1 ⊤ , . . . , β * K ⊤ ) ⊤ , if classification problem.

Figure 5 . 4 :

 54 Figure 5.4: Same interpretation as of Figure 5.3 but with N = 500000.

Figure 5 . 5 :Figure 5 . 6 :

 5556 Figure 5.5: Same interpretation as of Figure 5.3 but with N = 1000000.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 6.2 Future directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128

T

  ω k (t)ω ℓ (t)dt 1⩽k,ℓ⩽K ∈ R K×K ,then I ω is a non-negative define matrix, hence there exits I 1/2 ω ∈ R K×K such that I ω = I

A

  Figure A.1: Spectrometric curves in the Tecator data.

Figure A. 2 :

 2 Figure A.2: Left panel: height trajectories from age 1 to 18 of 39 boys and 54 girls. Right panel: height at age 18 of the two groups.

X w ω∥ 2 2 -

 2 5.4), one can obtainB.3 EM-iFME for updating iFME model parameters B.3.1 Updating the gating network parametersThis section presents the using of Dantzig selector to solve problem(3.36). Let us simplify the subscript k in the notations and rewrite the problem under matrix form as follows max χ∥Ω ω∥ 1 subject to A ω = 0 q , (B.9)where ω = (α 0 , ω [d 1 ] ⊤ , ω [d 2 ] ⊤ ) ⊤ is the unknown coefficients vector, c w = ( √ w 1 c 1 , . . . , √ w n c n ) ⊤

  4q+2 , ω + ⩾ 0 2q+1 , ω -⩾ 0 2q+1 , (B.10)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ω * = (α * 0 , ω [d 1 ] * ⊤ , ω [d 2 ] * ⊤ ) ⊤ tothe original problem could be retrieved by the relation ω * = ω * +ω * where ( ω * ⊤ + , ω * ⊤ -) ⊤ is the solution of (B.10). 140 B.3 EM-IFME FOR UPDATING IFME MODEL PARAMETERS B.3.2 Updating the expert network parameters

[d 1

 1 ] k ), it corresponds to solving a penalized weighted least square problem with constraints. From now on, let us simplify the subscript k in the notations and rewrite the problem under matrix form as followsmax 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 weighteddesign 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 (3.27)), 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 (B.9) 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 .

  4p+2 , γ + ⩾ 0 2p+1 , γ -⩾ 0 2p+1 , (B.11) 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 (B.11).

p→

  π k (x s ; α * ) by the Slusky theorem. It followsM m=1 λ m π k (x s ; α (m) ) p → π k (x s ; α * ), ∀x s ∈ X S .(C.3)

P

  ℓk (g (t) , x) = π ℓ (x), ∀ℓ ∈ [M K].

P

  ℓk (g (t) , x) c( φ ℓ (•|x), φ (t) k (•|x)) .

σ 2 4 .

 24 σ 2 ℓ -(ββ ℓ ) ⊤ xx ⊤ (ββ ℓ ) 2σReplacing the expectation operator by the empirical expectation and setting ∂K ∂β to s ) β ℓ .Let us denoteB = β 1 , . . . , β M K ∈ R p×M K , k (x 1 ) P (t) M K,k (x 2 ) . . . P (t) M K,k (x S )

  s ) σ 2 ℓ + (ββ ℓ ) ⊤ xx ⊤ (ββ ℓ ) ,which can be written under matrix form as trace(D (t)

  is modeled by another ME model, i.e., Expert k (y|x; θ k ) = is the number of experts connected to the kth lower-level gating network. Here, the function Expert kj (y|x) itself can be further modeled by another ME model and so on, resulting in higher level HME models. An illustration of a two-level HME architecture for regression is shown in Figure2.2.

				y
				g 2
				Gate
				g
				x
	g 1|1			g 1|2
	Gate			Gate
	g 2|1			g 2|2
	x			x
	where J k Expert 11 x	Expert 12 x	J k j=1	Gate j|k (x) Expert kj (y|x), Expert 21 Expert 22 x x	(2.7)

HME architecture was first introduced in Jordan and Jacobs (1994), since then it has been CHAPTER 2 STATE OF THE ART 1 Figure 2.2: Example of a two-level HME model. The two ME (rounded by the two dashed rectangles) are combined with the gate at the top to produce the final response. Each ME model in the dashed rectangles is itself composed of two experts and a gating network.

  Suppose that we observe n i.i.d. realizations {Y i (t)} n i=1 of a real-valued smooth function Y (t), t ∈ T , corresponding to n i.i.d. realizations {x i } n i=1 of a p-vector valued covariate x. The functionon-scalar regression model is given by

  (resp. D u ) is the derivative operator w.r.t. t (resp. u), and λ t , λ u are the smoothness penalty constants. For example, D t and D u are usually chosen to be the second-order partial derivative operators ∂ 2 ∂t 2 and ∂ 2 ∂u 2 . Let ψ(t) = [ψ 1 (t), . . . , ψ L (t)] ⊤ be a vector of basis functions defined on T . Firstly, we can express Y i (t) and X ⋆ ik (t) as follows

	2 dtdu,	(2.33)
	where D t

Table 3

 3

.1: Simulation scenarios with different noise level σ 2 δ and curve sampling level m.

Table 3

 3 .2 the results of regression and clustering tasks on simulated datasets in the four considered scenarios.

	RPE	Corr	ARI	ClusErr
	S1 (m = 300, σ 2 δ = 1)
	FME .1552 (.1282)	.9188 (.0602)	.7852 (.0934)	.0899 (.0434)
	FME-Lasso .1390 (.1074)	.9224 (.0603)	.7670 (.1436)	.1112 (.0891)
	iFME .1334 (.0977) .9287 (.0504) .7997 (.0796) .0792 (.0340)
	S2 (m = 100, σ 2 δ = 1)
	FME .1600 (.1698)	.9164 (.0712)	.7966 (.0881) .0852 (.0427)
	FME-Lasso .1656 (.1316)	.9071 (.0733)	.7455 (.1520)	.1236 (.0961)
	iFME .1540 (.0950) .9192 (.0476) .7955 (.0822)	.0820 (.0341)
	S3 (m = 300, σ 2 δ = 4)
	FME .1724 (.1807)	.9110 (.0778)	.7798 (.0838)	.0918 (.0415)
	FME-Lasso .1457 (.1504)	.9197 (.0785)	.7629 (.1524)	.1115 (.0928)
	iFME .1432 (.0860) .9228 (.0448) .7987 (.0899) .0783 (.0361)
	S4 (m = 100, σ 2 δ = 4)
	FME .2251 (.4462)	.9048 (.0822)	.7816 (.1040) .0927 (.0497)
	FME-Lasso .1432 (.1229) .9188 (.0683) .7526 (.1654)	.1215 (.1055)
	iFME .1639 (.0978)	.9125 (.0521)	.7798 (.0848)	.0864 (.0325)

Table 3 .

 3

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

Table 3 .

 3 3 shows the evaluation criteria corresponding to the considered models evaluated on 100 datasets in scenario S3. Here, FMR-PC is the original model ofYao et al. (2010) 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 data, thanks to the predictor-depending mixture weights, and provides clearly better results than the FMR alternatives. iFME .0139 (.0086) .9928 (.0046) .9594 (.0430) .0151 (.0166) FMR-B .0460 (.0399)

	RPE	Corr	ARI	ClusErr
	FME .0614 (.0896)	.9681 (.0469)	.8523 (.0934)	.0618 (.0490)
	FME-Lasso .0145 (.0254)	.9924 (.0133)	.9582 (.0641)	.0174 (.0325)
		.9768 (.0201)	.7064 (.1044)	.1124 (.0434)
	FMR-PC .0191 (.0331)	.9902 (.0164)	.8345 (.0794)	.0612 (.0361)

Table 3 .

 3 3: 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 4: Average of 100 Monte Carlo runs of MSE between the estimated functions resulted by FME, FME-Lasso and iFME models in S1 scenario.

	3.4 EXPERIMENTAL STUDIES

Table 3 .

 3 5: 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.

Table 3 .

 3 6: 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 3 .

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

  .10 (for example, top-right panel).

				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)

  However, to encourage the smoothness of the bivariate functions 4.3 FME FOR FUNCTION-ON-FUNCTION REGRESSION β k (t, u), we propose to consider the following penalized log-likelihood function

Table 5 .

 5 1 summarizes our notations for this chapter.

	Symbol	Explanation
	θ	Parameter vector of the model
	x i , y i	Training data. Without subscript they refer to generic variables
	K, M	Number of components of the true model, number of local machines

Table 5 .

 5 2: Comparison of Global estimator and Reduction estimator on MMASH data.

	.8.

125 CHAPTER 5 DISTRIBUTED LEARNING FOR MIXTURES OF EXPERTS

:

  Mixture of Experts; Functional Data Analysis; Regularized Estimation; LASSO; EM algorithm; MM Algorithm; Distributed Learning; Optimal Transport; Prediction; Clustering. Modélisation et Apprentissage avec des Mélanges d'Experts pour des Données Fonctionnelles et des Données Distribuées Résumé: Les modèles de mélange d'experts (ME) sont populaires en statistique et en apprentissage automatique et ont été étudiés pour des données vectorielles de haute dimension et qui sont centralisées. Cependant, dans de nombreux problèmes, nous observons des données de séries temporelles et l'application des modèles ME existants directement à ces données brutes peut limiter les performances car ils ignorent la corrélation entre les variables, une prropriété intrinsèque à ces données fonctionnelles, et ne capturent donc pas adéquatement la structure fonctionnelle inhérente à ces données. Dans de nombreuses applications, les données peuvent également ne pas être disponibles en mode centralisé, et il est donc nécessaire de développer des stratégies distribuées adaptées permettant des calculs parallèles efficaces. Cette thèse étudie les modèles de mélange d'experts (ME) pour les données fonctionnelles et les problèmes à grande échelle, dans un scénario hétérogène. Les objectifs principaux sont i) d'adresser des situations dans lesquelles on cherche à prédire une réponse potentiellement fonctionnelle à partir d'un prédicteur fonctionnel (par exemple, des séries temporelles), et ii) d'apprendre à partir de données distribuées. Dans cette thèse, nous proposons premièrement une nouvelle famille de modèles ME, appelée mélanges d'experts fonctionnels (FME), qui inclut le FME scalaire-sur-fonction, le FME fonction-sur-scalaire, et le FME fonction-sur-fonction. Nous introduisons des approches d'estimation régularisée via des régularisations appropriées qui encouragent des solutions parsimonieuses et des modèles lisses, tout en étant interprétables. Nous développons des algorithmes EM efficaces pour maximiser les fonctions de log-vraisemblance pour des données observées et menons des études expérimentales approfondies pour mettre en évidence les performances des modèles et algorithmes proposés. Ensuite, pour étendre l'estimation des mélanges d'experts ME au cas des données potentiellement distribuées, nous développons une stratégie d'apprentissage distribué des modèles ME. Elle effectue une inférence standard sur des machines locales pour obtenir des estimateurs locaux qui sont agrégés au niveau central. En se basant sur la minimisation d'une divergence de transport espérée proposée, les estimateurs locaux sont agrégés pour obtenir un estimateur réduit qui est consistant avec l'estimateur global, i.e., l'estimateur qui pourrait éventuellement être construit à partir de l'ensemble complet de données dans la version non-distribuée. Des études expérimentales démontrent la performance de notre approche.

*** Mots-clés: Mélanges d'Experts; Analyse de Données Fonctionnelles; Techniques de régularisation; LASSO; Algorithmes EM; Algorithmes MM; Apprentissage Distribué; Transport Optimal; Prédiction; Clustering.

occasionally called discrete distribution

Comparative Genomic Hybridization. CGH profiles are oftenly used to identify molecular subtypes of cancer

Diffusion Tensor Imaging

Acknowledgements

Chapter 4

Extensions of functional mixtures of experts Contents 4.2 FME FOR CLASSIFICATION m = 100 evenly spaced points. To mimic real-world data, for each i ∈ [n], we add to X i (t j ), j ∈ [m], measurement error δ i (t j ) ∼ N (0, σ 2 δ ), with σ 2 δ the noise level. We consider two situations: σ 2 δ = 1 and σ 2 δ = 5. We compare our proposed models, including FME fitted by EM, FME fitted by EM-Lasso, and iFME-EM model, with the classic functional multinomial logistic regression (FMLR) model that receives the coefficient vectors x i 's as inputs.

Each dataset is split into training and testing set of equal sizes. Parameters are fits on training set and the correct classification rates are evaluated on testing set. This procedure is repeated 100 times.

The classification results obtained with the described algorithms FME-EM, FME-EM-Lasso and iFME-EM, as well as with functional multinomial logistic regression (FMLR), are given in Table 4.1 and show higher classification performance of the iFME-EM approach.

Model Correct classification rate

Noise level: σ 2 δ = 1 Noise level: σ 2 δ = 5 FME-EM .8560 (.0199) .8474 (.0196) FME-EM-Lasso .9332 (.0104) .9178 (.0142) iFME-EM .9346 (.0108) .9219 (.0127) FMLR .7951 (.0249) .7922 (.0270) Table 4.1: Correct classification rates obtained on testing data. The reported values are averages on 100 samples with standard errors in parentheses.

■ Application to phonemes data

We then applied the two algorithms allowing for sparsity (FME-EM-Lasso and iFME-EM) to the well-known phoneme data [START_REF] Hastie | Penalized Discriminant Analysis[END_REF]. The data consists of n = 1000 log-periodogram recordings of length 256 each, used here as the univariate functional predictors, of five phonemes transcribed (i.e., the corresponding class labels): "sh" as in she, "dcl " as in dark, "iy" as the vowel in she, "aa" as the vowel in dark, and "ao" as the first vowel in water.

The obtained averaged correct classification rate for the two approaches are more than 0.94 in mean, which are competitive with the correct rate of 0.929 resulted by the wavelet-based multinomial functional regression (MFR) model [START_REF] Mousavi | Multinomial functional regression with wavelets and lasso penalization[END_REF]. Figure 4.1 shows the estimated coefficient functions for the gating network α k (t) (left) and those for the expert network β kg (t) as functions of sampling time t, obtained by FME-EM-Lasso (top) and the iFME-EM (bottom). Here the iFME-EM is fitted with constraints on the zero and the second derivatives of the coefficients functions. The results show clearly sparse and piece-wise-linear gating and experts functions when using the iFME-EM approach.

Chapter 5

Distributed learning for mixtures of experts Then θR is a consistent estimator of θ * .

The assumption that the parameters are ordered and initialized in Assumption A1 is necessary because under it the ME models are identifiable [START_REF] Jiang | On the identifiability of mixtures-of-experts[END_REF]. The proof of Proposition 5.3.3 is given in Appendix C.1.3.

An MM algorithm for constructing the reduction estimator

The main algorithm

We see that, the optimization problem (5.16) is accompanied by the objective function R c (g) defined upon another optimization, with respect to the transportation plan P , that makes the approach such as gradient descent hard to apply directly. In this section, we present the approach of using the MM algorithm, which requires the definition of a so-called majorant function, to solve the problem (5.16). We will apply it to the two common specifications of ME models: softmax-gated mixtures Gaussian regression and softmax-gated mixtures logistic regression models.

MM algorithm [START_REF] Lange | The mm algorithm[END_REF]) is an iterative procedure that consists of alternating between two steps: Majorize the objective function at the current iterate with a majorant function, and Minimize the majorant function to define the next iterate. The sequence of minimizers of the majorant functions is guaranteed to converge to a stationary point of the objective function. Employing MM algorithm to our problem means that we will starts from an initial model g (0) ∈ M K , then at each iteration tth, find a majorant function for R c (g) at g (t) , and minimize it to obtain the next model g (t+1) . The generated sequence (g (t) ) t⩾1 satisfies R c (g (0) ) ⩾ R c (g (1) ) ⩾ . . . We keep running until convergence is reached, i.e., there is no significant change in the value of R c (g). The model obtained at convergence is then supposed to be the desired model f R . Note that, for a general objective function, MM algorithm only guarantees to converge to a stationary point, not a global or local minimum, so multiple initializations are often required in order to obtain high quality solution [START_REF] Nguyen | An introduction to mm algorithms for machine learning and statistical estimation[END_REF].

The following proposition gives us a majorant function for R c (g). 5.23b) in which, the matrices

k and B are given in (C.6)-(C.8).

The detailed calculations are given in Appendix C.1.5.

Updating the logistic regression experts parameter

When the problem under consideration is classification, it is common to use logistic regression model to model the component densities. For simplicity, we consider here the case of binary responses, i.e., y ∈ {0, 1}. The conditional density in this case takes the following form

, where β ∈ R d is the parameter vector. Let φ 1 (y|x; β 1 ) and φ 1 (y|x; β 2 ) be the conditional densities of two binary logistic regression experts, then the KL-divergence between them is given by

.

(5.24)

Analogously to the case of Gaussian regression experts, substituting the expression of KLdivergence between two logistic experts into (5.19), the objective function can be rewritten as a function of β by .25) where

ℓk (x), similarly to the case of Gaussian experts, is calculated by (5.18) with c( φ ℓ (•|x), φ (5.24). CHAPTER 5 DISTRIBUTED LEARNING FOR MIXTURES OF EXPERTS and the estimated parameters are well-ordered, i.e., satisfy the identification condition for ME model [START_REF] Jiang | On the identifiability of mixtures-of-experts[END_REF], θW is a potential candidate.

Evaluation metrics

We compare the estimators in terms of the following metrics

• Transportation distance, defined in (5.10) with c being KL divergence, between the true mixture and the estimated mixture. This metric is used for simulated data where we know the true mixture.

• Log likelihood value of the estimated parameter evaluated on the testing set.

• Mean squared error (MSE) between the true parameter and the estimated parameter, used for simulated data.

• Relative prediction error (RPE) on testing set, defined by RPE = n i=1 (y iy i ) 2 )/ n i=1 y 2 i , where y i and y i are, respectively, the true and the predicted responses of the ith observation in the testing set.

• Adjusted Rand Index (ARI) between the true clustering in the testing set and the clustering predicted by the estimated mixture. ARI is used for simulated data where we know prior labels of the clusters.

• Learning time, i.e., the time consumed to obtain the estimator.

Simulated data

Data generating process

In this simulation, we fixed K = 4 and d = 20. We considered datasets of different sizes, from moderate to very large, to illustrate the statistical performance of the estimators. The datasets were generated as follows. Firstly we specified a true mixture of K components, which to be estimated, and the centers

Then for each N ∈ {100000, 500000, 1000000}, we generated 100 random dataset D = {(x i , y i )} N i=1 by drawing, for each k, N/K points from a multivariate Gaussian distribution with mean x (k) and covariance matrix defined by

. Finally, the responses y i are generated by the standard generating process for ME model, i.e.,

where M(1, π) is the multinomial distribution with parameter vector π. Here, the true parameter vector of the model to be estimated is θ

Application to MMASH dataset

Data description and pre-processing

In this study, we investigate the prediction performance of all the learning approaches (i.e., global, distributed, median and weighted averaged) to the Multilevel Monitoring of Activity and Sleep in Healthy people (MMASH) dataset. This dataset consists of 24 hours of continuous inter-beat interval data (IBI), triaxial accelerometer data, sleep quality, physical activity and psychological characteristics of 22 healthy young males. During the 24-hours data recording, participants wore two devices that continuously recorded heartbeats and beat-to-beat interval (via a heart rate monitor), and actigraphy information such as accelerometer data, sleep quality and physical activity (via an Actigraph device). The data was recorded on a second-by-second basis for each subject.

More details about the experiment setup of MMASH are provided in [START_REF] Rossi | A public dataset of 24-h multi-levels psycho-physiological responses in young healthy adults[END_REF], [START_REF] Oyeleye | A predictive analysis of heart rates BIBLIOGRAPHY using machine learning techniques[END_REF].

While MMASH is a rich dataset to assess many relations between physical, psychological and physiological characteristics, in this study we consider the problem of predicting the upcoming heart rate (HR) given the most recent records from them. In particular, we applied the data preprocessing procedure as in [START_REF] Oyeleye | A predictive analysis of heart rates BIBLIOGRAPHY using machine learning techniques[END_REF] to extract a HR series dataset that contains HR series of all participants in the MMASH dataset. Then, for each participant, his/her HR series is reframed into a set of predictor-response data pairs, i.e., into supervised setting. More specifically, we will use the last HRs during 300 seconds, to predict the 5-minutes-ahead HR. Each predictor can be viewed as a sample of 300 discretized values U (t j ), j ∈ [300], of some HR function X(t), and the 5-minutes-ahead HR (to be predicted) is its scalar response. Therefore, the FME model proposed in Chapter 3 will be the choice in this context.

The recording times vary among the participants, so the number of pairs (X i (t), y i ) extracted from their HR series are different from each other. However, the differences are not large, and the average is about 66000 pairs for each participant. In total, our dataset contains 1453000 data pairs (U i (t), y i ). The dataset is split into training and testing data with a ratio of 80%-20% (to ensure the balance, this split is also made w.r.t. each participant). The size of our training set is N = 1162400, and the size of testing set is 290600. Figure 5.7 shows some randomly taken predictors U i (t), as well as their smooth versions using a 50-dimensional B-spline basis. This appendix presents the detailed calculations to obtain the penalized least squares estimator (2.27). The problem under consideration now is minimize The equation (A.2) can thus be continued as

in which the last equality can be verified by a simple vectorization step that uses the definition of the matrices.

The second term of the objective function in (A.1) can be written in terms of vec(B ⊤ ) as

where Λ = diag(λ 1 , . . . , λ p ). Gathering (A.4) and (A.5), the problem (A.1) is now reformulated to minimize

which takes the same form as the usual ridge regression. Therefore, by calculating the first order derivative with respect to vec(B ⊤ ) and setting it to zero, the solution to (A. 6) is established as given in (2.27).

A.1.2 Penalized LS estimation of function-on-function regression

This appendix presents the detailed calculations to obtain the penalized least squares estimator (2.37) for the function-on-function regression problem. First, by plugging (2.33) into (2.32) and using the vectorization representation of the penalty term, the problem under consideration is now rewritten as

where Y (t) and X ⋆ (t) are defined as in (2.31).

From the relations (2.31) and (2.34), the first term of the objective function in (A.7) can be expressed as

with the second equality was obtained using the same manner as in (A.2) and (A.3).

Appendix B B.1 EM for the FME model

The FME model can be fitted by iteratively maximizing the observed-data log-likelihood (3.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 loglikelihood (B.1), given the observed data D, and the current parameter vector Ψ (s) :

where τ (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 (3.15).

M-step

This step updates the value of the parameter vector Ψ by maximizing the Q-function (B.2) 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.

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

The EM-Lasso algorithm for the maximization of (3.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 (B.1). 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 (3.18)).

E-step. This step computes the expectation of the complete-data log-likelihood (B.4), 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 (3.15).

M-step. This step updates the value of the parameter vector Ψ by maximizing the Q-function (B.5) 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.2.1 Updating the gating network parameters

Updating the gating network parameters at the sth EM iteration consists of maximizing the following function -1) is the gating network parameter vector and Ψ (s) is the current estimation of model's parameters. One can see this is equivalent to solving 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 ξ (t+1) = (1ν)ξ (t) + ν ξ(t) 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.2 Updating the experts network parameters

The maximization step for updating the expert parameters θ k consists of maximizing the function

where

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 τ

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

•

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

where the latter is the optimal solution of (B.8). Finally, the update for σ 2 k is given by

Hence, the update for the value of parameters vector Ψ at M-step, i.e., the solution to problem (B.6), is

) where ξ (s+1) and θ

), respectively, using the algorithms described above. The EM algorithm monotonically increases (3.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].

The inequality in the opposite direction is also true. Indeed, let g ⋆ = arg inf g∈M K R c (g) be a minimizer of R c (g), i.e., we have R c (g ⋆ ) = inf g∈M K R c (g). Then, in the remainder we will show that inf

We denote the gating and expert functions of g ⋆ by {π ⋆ (x),

According to (5.15), the gating functions are calculated by

where P ℓk (g ⋆ , x) denotes the entry (ℓ, k) of P(g ⋆ , x) as defined in (5.13). The relation (C.1) means that the matrix P(g ⋆ , x) is belonging to Π x (•, π ⋆ ). On the other hand, by its definition P(g ⋆ , x) is obviously belonging to Π x ( π, •). Hence, P(g ⋆ , x) ∈ Π x ( π, π ⋆ ). Note that, this holds for all x ∈ X . Therefore, taking expectation we have inf

C APPENDIX C which completes the proof.

C.1.2 Proof of Proposition 5.3.2

Continuity. Firstly, we show that, conditional on x, I(g, x) is continuous in g. Indeed, because C(P , φ) is a continuous function, and Π x ( π, •) is a compact set for all x ∈ X , the function defined by φ → inf

C(P , φ) is continuous, moreover, the infimum is attained. By definition, I(g, x)

depends on g only through φ, it follows that I(g, x) is also continuous in g.

Then, since I(g, x) is bounded for all g and x by the assumption, from the Lebesgue's dominated convergence theorem we have that: for all sequence (g (t) ) t⩾1 such that g (t) → g,

Convexity. We observe that Π x ( π, •) is a convex set for all x ∈ X (it can be verified easily by definition). It follows that I(g, x) is convex in g because of the facts that C(P , φ) is convex and taking infimum over a convex set preserves convexity. Then, let g 1 ,

which is the definition for the convexity of R c (g).

C.1.3 Proof of Proposition 5.3.3

For each m ∈ [M ], the assumption that the local estimator θ m is consistent implies that its expert parameters β (m) 1 , . . . , β (m) K must be ordered and its gating parameter α (m) must be initialized as those of θ * . Therefore, their weighted average θW = M m=1 λ m θ m is also a consistent estimator of θ * , i.e., we have θW → θ * in probability. It follows that the expected transportation divergence between the two corresponding ME models, f W and f * , will tend to zero almost surely as N → ∞, i.e., we have T ( f W , f * ) → 0.

On the other hand, by definition, the ME model corresponding to the reduction estimator θR is the one, among all g ∈ M K , minimizes the expected transportation divergence from the true model f * . Therefore, we have

almost surely. As a consequence, βR k → β * k in probability, because any possibility that it does not

List of Publications and Communications

Papers and communications • Pham, T. N. (November 2019). Quantile regression in large-scale data sets. Talk at LMNO Lab.

Other scientific activities

• Participation to "Working Group on Model-Based Clustering Summer Session", Vienne, 14-20 July 2019.

• Research visit at LIPN, University Paris-Nord, October-November 2020.

• Participation to "Conference: Math & IA", 9 March 2021.

• Participation to "MiMo 2021: Workshop on Mixture Models", Rouen, 8-9 April 2021.

• Participation to conference MHC2021, 2-4 June 2021

• Participation to "Third Conference on Statistics and Data Science", 28-30 October 2021.

List of Tables