A non-asymptotic approach for model selection via penalization in high-dimensional mixture of experts models

Mixture of experts (MoE) models, originally introduced as neural network architectures in [48, 50], are flexible models that generalize the classical finite mixture models as well as finite mixtures of regression models [66, Sec. 5.13]. Their flexibility comes from allowing the mixture weights (or the gating functions) to depend on the explanatory variables, along with the component densities (or experts). In the context of regression, MoE models with Gaussian experts and softmax or normalized Gaussian gating functions are the most popular choices and are powerful tools for modeling more complex nonlinear relationships between outputs (responses) and inputs (predictors) that arise from different subpopulations. The popularity of these conditional mixture density models arise largely due to their universal approximation properties, which have been studied in [49, 81, 74, 45, 73, 75], and which improve upon approximation capabilities of unconditional finite mixture models, as studied in [41, 85, 80, 43, 44, 78, 77]. Detailed reviews on the practical and theoretical aspects of MoE models can be found in [101, 59, 72, 76].

In this paper, we wish to investigate MoE models with Gaussian experts and normalized Gaussian gating functions for clustering and regression, first introduced by [99], which extended the original MoE models of [48]. From hereon in, we refer to these models as Gaussian-gated localized MoE (GLoME) models and block-diagonal covariance for localized MoE (BLoME) models, and we provide their precise definitions in Section 2.1. Furthermore, the original MoE models with softmax gating functions will be referred to as linear-combination-of-bounded-functions softmax-gated MoE (LinBoSGaME) regression models. In particular, we simply refer to affine instances of LinBoSGaME models as softmax-gated MoE (SGaME) regression models. It is worth pointing out that the BLoME models generalize GLoME models by utilizing a parsimonious covariance structure, via block-diagonal structures for covariance matrices in the Gaussian experts. It is also interesting to point out that supervised Gaussian locally-linear mapping (GLLiM) and block-diagonal covariance for Gaussian locally-linear mapping (BLLiM) models in [26] and [34] are affine instances of GLoME and BLoME models, respectively, where the latter pair consider general linear combination of bounded functions instead of affine for mean functions of Gaussian experts, in the prior pair.

One of the main disadvantages of LinBoSGaME models is the difficulty of applying an expectation–maximization (EM) algorithm [27, 65], which requires an internal iterative numerical optimization procedure (e.g., iteratively-reweighted least squares, Newton-Raphson algorithm) to update the softmax parameters. To overcome this problem, we instead use the Gaussian gating network that enables us to link BLoME with hierarchical Gaussian mixture model. Then, the maximization with respect to the parameters of the gating network can be solved analytically with the EM algorithm framework, which decreases the computational complexity of the estimation routine. Furthermore, we can then also make use of well-established theoretical results for finite mixture models.

We note here that both GLoME and BLoME models have been thoroughly studied in the statistics and machine learning literatures in many different guises, including localized MoE [86, 87, 69, 15], normalized Gaussian networks [90], MoE modeling of priors in Bayesian nonparametric regression [83, 82], cluster-weighted modeling [47], deep mixture of linear inverse regressions [55], hierarchical Gaussian locally linear mapping structured mixture (HGLLiM) model [95], multiple-output Gaussian gated mixture of linear experts [73], and approximate Bayesian computation with surrogate posteriors using GLLiM [39].

Our main goal is to learn potentially nonlinear relationships between a multivariate output and a high-dimensional input issued from a heterogeneous population. This involves performing regression, clustering and model selection, simultaneously. While estimation can be performed using standard EM algorithms, it crucially depends and requires data-driven hyperparameter choices, including the number of mixture components (or clusters), the degree of complexity of each Gaussian expert’s mean function, and the hidden block-diagonal structures of the covariance matrices. Recall that hyperparameter choices of data-driven learning algorithms belong to the class of model selection problems, which have attracted a lot of attention in statistics and machine learning over the past 50 years [1, 58, 2, 60, 3]. Specifically, given a set of models, how do we select the one with the lowest possible risk from the data? Note that penalization is one of the main strategies proposed for model selection. It suggests choosing the estimator that minimizes the sum of its empirical risk and some penalty terms corresponding to the fit of the model to the data while avoiding overfitting.

In general, model selection for MoE models is often performed using the Akaike information criterion (AIC) [1], the Bayesian information criterion (BIC) [91] or the BIC-like approximation of integrated classification likelihood (ICL-BIC) [10]. An important limitation of these criteria, however, is that they are only valid asymptotically. This implies that there is no finite-sample guarantee when using AIC, BIC or ICL-BIC, to choose between different levels of complexity. Their use in small samples is, therefore, ad hoc. To overcome such difficulties, [11] proposed a novel approach, called the slope heuristics, supported by a non-asymptotic oracle inequality. This method leads to an optimal data-driven choice of multiplicative constants for the penalties. Recent reviews and practical issues concerning the slope heuristic can be found in [8, 3] and the references given therein.

It should be emphasized that a general model selection result, originally established by [60, Theorem 7.11], ensures a penalized criterion leads to a good model selection and that the penalty is known only up to multiplicative constants and is proportional to the dimensions of models. In particular, these multiplicative constants can be calibrated by the slope heuristic approach in a finite-sample setting. Then, in the spirit of the methods based on concentration inequalities developed in [60, 61, 22, 23], a number of finite-sample oracle type inequalities have been established for the least absolute shrinkage and selection operator (LASSO) [94] and general penalized maximum likelihood estimators (PMLE). These results include the works for high-dimensional Gaussian graphical models [33], Gaussian mixture model selection [62, 63, 64], finite mixture regression models [68, 28, 29, 31, 32], and LinBoSGaME models, outside the high-dimensional setting [70].

To the best of our knowledge, no attempt has been made in the literature to develop a finite-sample oracle inequality for the framework of MoE regression models for high-dimensional data. In this paper, our first original contribution is to provide finite-sample oracle inequalities, Theorems 3.1 and 3.2 for high-dimensional GLoME and BLoME models, respectively. More specifically, we establish non-asymptotic risk bounds that take the form of weak oracle inequalities, provided that the lower bounds on the penalties are true, in high-dimensional GLoME and BLoME models, based on an inverse regression strategy and block-diagonal structures of the Gaussian covariance matrices experts. Unlike traditional criteria such as AIC, BIC, or ICL-BIC, which are based on asymptotic theory or a Bayesian approach, our contributed non-asymptotic risk bounds allow the number $n$ of observations to be fixed while the dimensionality and cardinality of the models, characterized by the number of covariates and the size of the response, are allowed to grow with respect to $n$ and can be much larger than $n$.

In particular, our oracle inequalities show that the Jensen–Kullback–Leibler loss performance of our PMLEs is comparable to that of oracle models if we take sufficiently large constants in front of the penalties, whose shapes are only known up to multiplicative constants and proportional to the dimensions of the models. These theoretical justifications for the shapes of the penalties motivate us to make use of the slope heuristic criterion to select several hyperparameters, including the number of mixture components, the degree of polynomial mean functions, and the potential hidden block-diagonal structures of the covariance matrices of the multivariate predictors.

Moreover, in Theorem 3.1, we extends a corollary of [70, Theorem 1], which can be verified via Lemma 1 from [75], which makes explicit the relationship between softmax and Gaussian gating classes.

Another significant contribution is the numerical experiments for simulated and real data sets in Section 5 that support our theoretical results, and the statistical study of non-asymptotic model selection in GLoME models. To the best of our knowledge, instead of using classical asymptotic approaches for model selection such as AIC, BIC and ICL-BIC, we are the first to illustrate that the slope heuristic works well for MoEs with Gaussian experts and normalized Gaussian gating networks such as GLLiM, BLLiM, GLoME, and BLoME models. Note that our main objective here is to investigate how well the empirical tensorized Kullback–Leibler divergence between the true model ($s_0^{\ast }$) and the selected model ${\widehat{s}}_{\widehat{𝐦}}^{\ast }$ follows the finite-sample oracle inequality of Theorem 3.1, as well as the convergence rate $𝒪(n^{-1})$ of the error upper bound. Therefore, we focus on $1$-dimensional data sets. Beyond the statistical estimation and model selection purposes considered here, the dimensionality reduction capability of GLLiM and BLLiM models in high-dimensional regression data, can be found in [26, Section 6] and [34, Sections 3 and 4], respectively.

Specifically, besides the important theoretical issues regarding the tightness of upper bounds, how to integrate a priori information, and a minimax analysis of our proposed PMLEs, our finite-sample oracle inequalities and the corresponding illustrative numerical experiments help to partially answer the following two important questions raised in the field of MoE regression models: (1) What is the number of mixture components $K$ to choose, given the sample size $n$; and (2) Whether it is better to use a few complex experts or a combination of many simple experts, given the total number of parameters. Note that such problems are also considered in the work of [67, Proposition 1], where the authors provided some qualitative guidance and only suggested a practical method for choosing $K$ and $d$, involving a complexity penalty or cross-validation. Furthermore, their model is only non-regularized maximum-likelihood estimation, and thus is not suitable for the high-dimensional setting.

Lastly, we emphasize that although the finite-sample oracle inequalities compare performances of our estimators with the best model in the collection, they also allow us to approximate a rich class of conditional densities if we take enough clusters and/or a high enough degree $d$ of polynomials in the mean of the functions of Gaussian experts, in the context of mixture of Gaussian experts [49, 67, 74, 45, 75]. This leads to the upper bounds of the risks being small.

For the Gaussian gating parameters, the technique for handling the logistic weights in the SGaME models of [70] is not directly applicable to the GLoME or BLoME framework, due to the quadratic form of the canonical link. Therefore, we have to propose a reparameterization trick¹¹1Note that we only use this nomenclature to perform a change of variables of the Gaussian gating parameters space of GLoME models via the logistic weights of SGaME models. This reparameterization trick does not stand for the well-known one of variational autoencoders (VAEs) in deep learning literature (see [52], for more details). to bound the metric entropy of the Gaussian gating parameters space; see (4.1.3) and Lemma 4.5 for more details.

To work with conditional density estimation in the BLLiM and BLoME models, it is natural to make use of a general conditional density model selection theorem from [22, 23], see also Theorem 4.1 for its adaption to our context. However, it is worth mentioning that since the collection of models constructed by the BLLiM model [34] or by some appropriate procedures for BLoME models in practice is usually random, we have to utilize a model selection theorem for MLE among a random subcollection (cf. [29, Theorem 5.1] and [33, Theorem 7.3]).

In particular, our collections of BLoME models must satisfy certain regularity assumptions, see Section 4.2.3 for details, which are not easy to verify due to the complexity of BLoME models and technical reasons. For BLoME models, the main difficulty in proving our oracle inequality lies in bounding the bracketing entropy of the Gaussian gating functions and Gaussian experts with block-diagonal covariance matrices. To solve the first problem, we need to use our previous trick of reparametrizing the Gaussian gating parameter space. For the second one, based on some ideas of Gaussian mixture models from [41, 62], we provide a novel extension for standard MoE models with Gaussian gating networks. Note that our important contributions also extend the recent result on block-diagonal covariance matrices in [33], which is only developed for Gaussian graphical models.

In addition, there are two possible ways to decompose the bracketing entropy of collection models based on different distances. As mentioning in Appendix B.2.1 from [70], their decomposition boils down to assuming that the predictor domain is bounded. And this limiting assumption can then be relaxed by using a smaller distance, namely a tensorized extension of the Hellinger distance ${\mathrm{d}}^{\otimes \mathrm{n}}$ in (2.24), but bounding the corresponding bracketing entropy becomes much more difficult. Nevertheless, it is important to point out that we are able to weaken this limiting assumption by using the smaller distance: ${\mathrm{d}}^{\otimes \mathrm{n}}$, for the bracketing entropy of our collection of BLoME models, see Lemma 4.10 for more details. This reinforces our original contributions concerning the control of bracketing entropy of BLoME models.

The rest of this paper is organized as follows. In Section 2, we present the notation and framework for GLoME, and BLoME models and their special cases, GLLiM, and BLLiM models. In Section 3, we state the main results of this paper: finite-sample oracle inequalities satisfied by the PMLEs. Section 4 is devoted to the proofs and main mathematical challenges in establishing our new finite-sample oracle inequalities. We experimentally evaluate our new results in simulated and real datasets in Section 5. Some conclusions and perspectives are provided in Section 6. The proofs of technical lemmas can be found in Appendix A.

Motivated by an inverse regression framework, where the role of predictor and response variables will be exchanged such that $𝐘$ becomes the input and $𝐗$ plays the role of a multivariate output, we consider GLoME and BLoME models with inverse conditional probability density functions (PDFs) of the form (2.1). This construction goes back to the work of [57, 99, 26, 84], and is very useful in a high-dimensional regression context, where typically $D\gg L$. In this way, we define PDF of $𝐗$ conditional on $𝐘=𝐲$ as follows:

Here, $g_k(\cdot ;𝝎)$ and ${\varphi }_D(\cdot ;{𝝊}_{k,d}(\cdot ),{𝚺}_k)$, $k\in [K]$, $K\in {\mathbb{N}}^{\star }$, $d\in {\mathbb{N}}^{\star }$, are called normalized Gaussian gating functions (or Gaussian gating networks) and Gaussian experts, respectively. Furthermore, we decompose the parameters of the model as follows: ${𝝍}_{K,d}=(𝝎,{𝝊}_d,𝚺)\in {𝛀}_K\times {𝚼}_{K,d}\times {𝐕}_K=:{𝚿}_{K,d}$, $𝝎=(𝝅,\mathrm{c},𝚪)\in ({𝚷}_{K-1}\times {𝐂}_K\times {\mathrm{V}}_K^{\prime })=:{𝛀}_K$, $𝝅={\left({\pi }_k\right)}_{k\in [K]}$, $\mathrm{c}={\left({𝐜}_k\right)}_{k\in [K]}$, $𝚪={\left({𝚪}_k\right)}_{k\in [K]}$, ${𝝊}_d={\left({𝝊}_{k,d}\right)}_{k\in [K]}\in {𝚼}_{K,d}$, and $𝚺={\left({𝚺}_k\right)}_{k\in [K]}\in {𝐕}_K$. Note that ${𝚷}_{K-1}=\left\{{\left({\pi }_k\right)}_{k\in [K]}\in {\left({ℝ}^{+}\right)}^K,{\sum }_{k=1}^K{\pi }_k=1\right\}$ is a $K-1$ dimensional probability simplex, ${𝐂}_K$ is a set of $K$-tuples of mean vectors of size $L\times 1$, ${\mathrm{V}}_K^{\prime }$ is a set of $K$-tuples of elements in ${𝒮}_L^{++}$, where ${𝒮}_L^{++}$ denotes the collection of symmetric positive definite matrices on ${ℝ}^L$, ${𝚼}_{K,d}$ is a set of $K$-tuples of mean functions from ${ℝ}^L$ to ${ℝ}^D$ depending on a degree $d$ (e.g.,  a degree of polynomials), and ${𝐕}_K$ is a set containing $K$-tuples from ${𝒮}_D^{++}$.

Recall that GLLiM and BLLiM models are affine instances of GLoME and BLoME models, and are especially useful for high-dimensional regression data, since there exist link functions between the inverse and forward conditional density; see Figure 1 for comprehensive classification and nomenclature of standard MoE regression models with Gaussian gating networks. Note that the principle of inverse regression is only useful when the functions ${𝝊}_{k,d}(𝐲)$ are linear, because there is no explicit way to express the law of $𝐘\mid 𝐗$ from that of $𝐗\mid 𝐘$ for higher degree of polynomials. However, to have more consistent notations with the previous affine results of GLLiM, BLLiM models from [26, 34], we decide to use the inverse regression frameworks instead of the forward one.

In the BLoME model, we wish to make use of block-diagonal covariance structures by replacing ${𝚺}_k$ and ${𝐕}_K$ with ${𝚺}_k\left({𝐁}_k\right)$ and ${𝐕}_K\left(𝐁\right)$, defined in (2.5), respectively (see, e.g.,  [34, 33]). This block-diagonal structures for covariance matrices are not only used to trade-off between complexity and sparsity but is also motivated by some real applications, where we want to perform prediction on data sets with heterogeneous observations and hidden graph-structured interactions between covariates; for instance, for gene expression data sets in which conditionally on the phenotypic response, genes interact with few other genes only, i.e., there are small modules of correlated genes (see [34, 33] for more details). To be more precise, for $k\in [K]$, we decompose ${𝚺}_k\left({𝐁}_k\right)$ into $G_k$ blocks, $G_k\in {\mathbb{N}}^{\star }$, and we denote by $d_k^{[g]}$ the set of variables into the $g$th group, for $g\in \left[G_k\right]$, and by $card\left(d_k^{[g]}\right)$ the number of variables in the corresponding set. Then, we define ${𝐁}_k={\left(d_k^{[g]}\right)}_{g\in \left[G_k\right]}$ to be a block structure for the cluster $k$, and $𝐁={\left({𝐁}_k\right)}_{k\in [K]}$ to be the covariate indexes into each group for each cluster. In this way, to construct the block-diagonal covariance matrices, up to a permutation, we make the following definition: ${𝐕}_K\left(𝐁\right)={\left({𝐕}_k\left({𝐁}_k\right)\right)}_{k\in [K]}$, for every $k\in [K]$,

Here, ${𝐏}_k$ corresponds to the permutation leading to a block-diagonal matrix in cluster $k$. It is worth pointing out that outside the blocks, all coefficients of the matrix are zeros and we also authorize reordering of the blocks: e.g., $\{(1,3);(2,4)\}$ is identical to $\{(2,4);(1,3)\}$, and the permutation inside blocks: e.g., the partition of $4$ variables into blocks $\{(1,3);(2,4)\}$ is the same as the partition $\{(3,1);(4,2)\}$.

It is interesting to point out that the BLLiM framework aims to model a sample of high-dimensional regression data issued from a heterogeneous population with hidden graph-structured interaction between covariates. In particular, the BLLiM model is considered as a good candidate for performing model-based clustering and for predicting the response in situations affected by the “curse of dimensionality” phenomenon, where the number of parameters could be larger than the sample size. Indeed, to deal with high-dimensional regression problems, the BLLiM model is based on an inverse regression strategy, which inverts the role of the high-dimensional predictor and the multivariate response. Therefore, the number of parameters to estimate is drastically reduced. More precisely, BLLiM utilizes GLLiM, described in [25, 26], in conjunction with a block-diagonal structure hypothesis on the residual covariance matrices, to make a trade-off between complexity and sparsity. This prediction model is fully parametric and highly interpretable. For instance, it can be used for the analysis of transcriptomic data in molecular biology to classify observations or predict phenotypic states; see e.g.,  [42, 71, 56]. Indeed, if the predictor variables are gene expression data measured by microarrays or by the RNA-seq technologies, and the response is a phenotypic variable, the BLLiM model not only provides clusters of individuals based on the relation between gene expression data and the phenotype, but also implies a gene regulatory network specific to each cluster of individuals (see [34] for more details).

In order to establish our finite-sample oracle inequalities, Theorems 3.1 and 3.2, we need to explicitly impose some classical boundedness conditions on the parameter space.

The GLLiM and BLLiM models, as originally introduced in [26, 34], are used to capture the nonlinear relationship between the response and the set of covariates from a high-dimensional regression data, imposed by a potential hidden graph structured interaction, typically in the case when $D\gg L$, by the $K$ locally affine mappings:

Here, $𝕀$ is an indicator function and $Z$ is a latent variable capturing a cluster relationship, such that $Z=k$ if $𝐘$ originates from cluster $k\in [K]$. Cluster specific affine transformations are defined by matrices ${𝐀}_k^{\ast }\in {ℝ}^{L\times D}$ and vectors ${𝐛}_k^{\ast }\in {ℝ}^L$. Furthermore, ${𝐄}_k^{\ast }$ are error terms capturing both the reconstruction error due to the local affine approximations and the observation noise in ${ℝ}^L$.

Following the common assumption that ${𝐄}_k^{\ast }$ is a zero-mean Gaussian variable with covariance matrix ${𝚺}_k^{\ast }\in {ℝ}^{L\times L}$, it holds that

where we denote by ${𝝍}_K^{\ast }$ the vector of model parameters and ${\mathrm{\Phi }}_L$ is the PDF of a Gaussian distribution of dimension $L$. In order to allow the affine transformations to be local, $𝐗$ is defined as a mixture of $K$ Gaussian components as follows:

where ${𝐜}_k^{\ast }\in {ℝ}^D,{𝚪}_k^{\ast }\in {ℝ}^{D\times D}$, ${𝝅}^{\ast }={\left({𝝅}_k^{\ast }\right)}_{k\in [K]}\in {𝚷}_{K-1}^{\ast }$, and ${𝚷}_{K-1}^{\ast }$ is the $K-1$ dimensional probability simplex. Then, according to formulas for conditional multivariate Gaussian variables and the following hierarchical decomposition

we obtain the following forward conditional density [26]:

where ${𝝍}_K^{\ast }=({𝝅}^{\ast },{𝜽}_K^{\ast })\in {\mathrm{\Pi }}_{K-1}\times {𝚯}_K^{\ast }=:{𝚿}_K^{\ast }$. Here, ${𝜽}_K^{\ast }={({𝐜}_k^{\ast },{𝚪}_k^{\ast },{𝐀}_k^{\ast },{𝐛}_k^{\ast },{𝚺}_k^{\ast })}_{k\in [K]}$ and

Without assuming anything on the structure of parameters, the dimension of the model, denoted by $dim(\cdot )$, is defined as the total number of parameters that has to be estimated, as follows:

It is worth mentioning that $dim\left({𝚿}_K\right)$ is typically very large compared to the sample size (see, e.g., [26, 34, 84] for more details in their real data sets) whenever $D$ is large and $D\gg L$. Furthermore, it is more realistic to make assumption on the residual covariance matrices ${𝚺}_k^{\ast }$ of error vectors ${𝐄}_k^{\ast }$ rather than on ${𝚪}_k^{\ast }$ (cf. [26, Section 3]). This justifies the use of the inverse regression trick from [26], which leads a drastic reduction in the number of parameters to be estimated.

More specifically, in (2.14), the roles of input and response variables were exchanged such that $𝐘$ becomes the covariates and $𝐗$ plays the role of the multivariate response. Therefore, its corresponding inverse conditional density is defined as a Gaussian locally-linear mapping (GLLiM) model, based on the previous hierarchical Gaussian mixture model, as follows:

where ${𝚺}_k$ is a $D\times D$ covariance structure (usually diagonal in GLLiM models, chosen to reduce the number of parameters) automatically learned from data and ${𝝍}_K$ is the set of parameters, denoted by ${𝝍}_K=(𝝅,{𝜽}_K)\in {𝚷}_{K-1}\times {𝚯}_K=:{𝚿}_K$. It is important to note that the BLLiM model imposes block-diagonal structures on ${\left({𝚺}_k\right)}_{k\in [K]}$ to make a trade-off between complexity and sparsity.

The following interesting feature of both GLLiM and BLLiM models is proved in Section A.1.

For GLoME, we only need to choose the degree of polynomials $d$ and the number of components $K$ among finite sets ${𝒟}_{𝚼}=\left[d_{\max }\right]$ and $𝒦=\left[K_{\max }\right]$, respectively, where $d_{\max }\in {\mathbb{N}}^{\star }$ and $K_{\max }\in {\mathbb{N}}^{\star }$ may depend on the sample size $n$. We wish to estimate the unknown inverse conditional density $s_0$ by inverse conditional densities belonging to the following collection of models ${\left(S_{𝐦}\right)}_{𝐦\in ℳ}$. Here, we define $ℳ=\{(K,d):K\in 𝒦,d\in {𝒟}_{𝚼}\}$, and

where $s_{{𝝍}_{K,d}}$, ${\widetilde{𝛀}}_K$, ${𝚼}_{K,d}$ and ${𝐕}_K$ are define previously in (2.1), (2.1.1), (2.8) (or more general (2.7)) and (2.1.2), respectively.

While for the BLoME model, we also need to select $𝐁$ from a list of candidate structures ${\left({ℬ}_k\right)}_{k\in [K]}\equiv {\left(ℬ\right)}_{k\in [K]}$, where $ℬ$ denotes the set of all possible partitions of the covariables indexed by $[D]$, for each cluster of individuals.The collection of BLoME models is defined as follows: ${\left(S_{𝐦}\right)}_{𝐦\in ℳ}$, $ℳ=\{(K,d,𝐁):K\in 𝒦,d\in {𝒟}_{𝚼},𝐁\in {\left(ℬ\right)}_{k\in [K]}\}$,

where $s_{{𝝍}_{K,d,𝐁}}$, ${\widetilde{𝛀}}_K$, ${𝚼}_{K,d}$, and ${\widetilde{𝐕}}_K\left(𝐁\right)$ are defined in (2.1), (2.1.1), (2.8) (or more generally (2.7)), and (2.1.2), respectively. In theory, we would like to consider the whole collection of models ${\left(S_{𝐦}\right)}_{𝐦\in ℳ}$. However, the cardinality of $ℬ$ is large; its size is a Bell number. Even for a moderate number of variables $D$, it is not possible to explore the set $ℬ$, exhaustively. We restrict our attention to a random subcollection ${ℬ}^R$ of moderate size. For example, we can consider the BLLiM procedure from [34, Section 2.2].

In the maximum likelihood approach, the Kullback–Leibler divergence is the most natural loss function, which is defined for two densities $s$ and $t$ by

However, to take into account the structure of conditional densities and the random covariates $\left({𝐘}_{[n]}\right)$, we consider the tensorized Kullback–Leibler divergence ${\mathrm{KL}}^{\otimes \mathrm{n}}$, defined as:

if $sdy$ is absolutely continuous w.r.t. $tdy$, and $+\mathrm{\infty }$ otherwise. Note that if the predictors are fixed, this divergence is the classical fixed design type divergence in which there is no expectation. We refer to our result as a weak oracle inequality, because its statement is based on a smaller divergence, when compared to ${\mathrm{KL}}^{\otimes \mathrm{n}}$, namely the tensorized Jensen-Kullback–Leibler divergence:

with $\rho \in (0,1)$. We note that ${\mathrm{JKL}}_{\rho }^{\otimes \mathrm{n}}$ was first used in [22, 23]. However, a version of this divergence appears explicitly with $\rho =\frac{1}{2}$ in [60], and it is also found implicitly in [12]. This loss is always bounded by $\frac{1}{\rho }\ln \frac{1}{1-\rho }$ but behaves like ${\mathrm{KL}}^{\otimes \mathrm{n}}$, when $t$ is close to $s$. The main tools in the proof of such a weak oracle inequality are deviation inequalities for sums of random variables and their suprema. These tools require a boundedness assumption on the controlled functions which is not satisfied by $-\ln \frac{s_{𝐦}}{s_0}$, and thus also not satisfied by ${\mathrm{KL}}^{\otimes \mathrm{n}}$. Therefore, we consider instead the use of ${\mathrm{JKL}}_{\rho }^{\otimes \mathrm{n}}$. In particular, in general, it holds that $C_{\rho }{\mathrm{d}}^{2\otimes \mathrm{n}}\le {\mathrm{JKL}}_{\rho }^{\otimes \mathrm{n}}\le {\mathrm{KL}}^{\otimes \mathrm{n}}$, where $C_{\rho }=\frac{1}{\rho }\min (\frac{1-\rho }{\rho },1)\left(\ln \left(1+\frac{\rho }{1-\rho }\right)-\rho \right)$ (see [22, Prop. 1]) and ${\mathrm{d}}^{2\otimes \mathrm{n}}$ is a tensorized extension of the squared Hellinger distance ${\mathrm{d}}^{2\otimes \mathrm{n}}$, defined by

Moreover, if we assume that, for any $𝐦\in ℳ$ and any $s_{𝐦}\in S_{𝐦},s_{0}d\lambda \ll s_{𝐦}d\lambda$, then (see [70, 22])

In the context of MLE, given the collection of conditional densities $S_{𝐦}$, we aim to estimate $s_0$ by the conditional density ${\widehat{s}}_{𝐦}$ that maximizes the likelihood (conditionally to ${\left({𝐲}_i\right)}_{i\in [n]}$) or equivalently that minimizes the negative log-likelihood (NLL), which we can write as:

We should work with almost minimizer of this quantity and define a $\eta$-log-likelihood minimizer as any ${\widehat{s}}_{𝐦}$ that satisfies:

where the error term $\eta$ is necessary to avoid any existence issue, e.g.,  the infimum may not be reached. Roughly speaking, the Ekeland variational principle asserts that, for any extended-valued lower semicontinuous function which is bounded below, one can add a small perturbation to ensure the existence of the minimum, see e.g.,  [14, Chapter 2]. This framework is also used in [22, 23, 70].

As always, using the NLL of the estimate in each model as a criterion is not sufficient. It is an underestimation of the risk of the estimate and this leads to choosing models that are too complex. In the context of PMLE, by adding a suitable penalty $\text{pen}(𝐦)$, one hopes to create a trade-off between good data fit and model complexity. For a given choice of $\text{pen}(𝐦)$, the selected model $S_{\widehat{𝐦}}$ is chosen as the one whose index is an ${\eta }^{\prime }$-almost minimizer of the sum of the NLL and this penalty:

Note that ${\widehat{s}}_{\widehat{𝐦}}$ is then called the ${\eta }^{\prime }$-penalized likelihood estimate and depends on both the error terms $\eta$ and ${\eta }^{\prime }$. From hereon in, the term selected model (estimate) or best data-driven model (estimate) is used to indicate that it satisfies the definition in (2.27).

Theorem 3.1 provides a lower bound on the penalty function, $\text{pen}(𝐦)$, which guarantees that the PMLE selects a best data-driven model that performs almost as well as the best model. We highlight our main contribution and briefly establish the proof of Theorem 3.1 in Section 4.1.3. A more detail technical proof can be found in Sections A.2 and A.3.

It is worth noting that Theorem 3.1 extends a corollary of [70, Theorem 1], which can be verified via Lemma 1 from [75], which makes explicit the relationship between softmax and Gaussian gating classes.

Note that we only aim to construct upper bounds and do not focus on the important question of the existence of the corresponding lower bounds or minimax adaptive estimation. However, as a special case of GLoME models, [62, 64] consider a collection of univariate Gaussian mixture models having the same and known component variance. They show that the PMLE ${\widehat{s}}_{\widehat{𝐦}}$ is minimax adaptive to the regularity $\beta$, $\beta >0$, of univariate density classes ${\mathscr{H}}_{\beta }$ whose logarithm of the elements is locally $\beta$-Hölder, with convergence rate $n^{-\frac{2\beta }{2\beta +1}}$ up to a power of $\ln (n)$. Although this result is stated for the Hellinger risk, it remains valid for the Kullback–Leibler divergence if we further assume that $\ln \left({\Vert s/t\Vert }_{\mathrm{\infty }}\right)$ is uniformly bounded on ${\cup }_{𝐦\in ℳ}{S}_{𝐦}$, see e.g.,  Lemma 7.23 in [60]. Note that this assumption guarantees that the Kullback–Leibler divergence and the Hellinger distance are equivalent.

A special case of GLoME models, namely model-selection (nearly-D-sparse) aggregation in mixture models, is considered by [17, 9, 89, 18, 24] and is related to our model selection results. More precisely, the authors of [89] did not consider PMLEs with Kullback-Leibler loss but only with $L_2$-loss or aggregation of a finite number of densities. Similarly, the results from [17, 9] dealt with the $L_2$-loss and investigate the Lasso and the Dantzig estimators, respectively, suitably adapted to the problem of density estimation. With respect to the Kullback–Leibler loss, [18] established model selection type oracle inequalities with high probability rather than in expectation. In particular, a bound in expectation with Kullback–Leibler loss in [24] is perhaps the most relevant reference to our result. They established exact oracle inequalities with a rate-optimal remainder term ${\left(\left(\ln K\right)/n\right)}^{1/2}$, up to a possible logarithm correction, in the problem of convex aggregation when the number $K$ of components is larger than $n^{1/2}$. It is worth noting that they did not consider PMLEs as we do or as [62, 64], see e.g.,  (2.27). Instead, the constraint that the weight vector in the mixture model belongs to the probabilistic simplex acts as a sparsity-inducing penalty. However, in their collection of mixture models, the mixture components are deterministic and chosen from a dictionary obtained on the basis of previous experiments or expert knowledge rather estimated from data. The adjective exact refers to the fact that the “bias term” ${\mathrm{KL}}^{\otimes \mathrm{n}}(s_0,s_{𝐦})$ is not multiplied by factor strictly larger than one as in our Theorem 3.1.

To the best of our knowledge, providing a minimax analysis for the PMLEs or the problem of model-selection (nearly-D-sparse) aggregation in the context of standard MoE models is still an open question. In particular, it is not trivial to extend such optimal risk bounds regarding Gaussian mixtures, see e.g.,  [64, Theorem 2.8] or [24, Theorems 2.3 and 3.1], to standard MoE models. However, we wish to provide such a minimax analysis in future research.

In contract to Theorem 3.1, in Theorem 3.2, we consider a random subcollection of models $\widetilde{ℳ}$, included in the whole collection $ℳ$. This is particularly useful in a high-dimensional context where we cannot test all the models. Therefore, we have to restrict ourselves to a smaller subcollection of models, which is then potentially random.

Note that the constructed collection of models with block-diagonal structures for each cluster of individuals is designed, for example, by the BLLiM procedure from [34], where each collection of partitions is sorted by sparsity level. Nevertheless, our finite-sample oracle inequality, Theorem 3.2, still holds for any random subcollection of $ℳ$, which is constructed by some suitable tools in the framework of BLoME regression models. We highlight our main contribution and briefly establish the proof of Theorem 3.2 in Section 4.2. A more detailed technical proof can be found in Sections A.4, 4.2.3 and A.5.

It is worth mentioning that the comparison of the two inequalities in Theorems 3.1 and 3.2 involves the following major differences. First, to control the random subcollection, in Theorem 3.2, we further assume condition (3.3). However, this is not a strong assumption and is satisfied, for example, if $s_0$ is bounded, with a compact support. Assumption (3.3) is needed because we consider a random subcollection from the whole collection. Thanks to this assumption, we can utilize Bernstein’s inequality to control the additional randomness. It is important to stress that the parameter $\tau$ depends on the true unknown density $s_0$ and cannot be explicitly determined for this reason. It appears not only in the oracle type inequality (3.5), but also in the penalty term (3.4). However, in some special cases of BLoME models, for instance finite mixture regression models, we can construct a larger penalty independent of $\tau$ to obtain an oracle type inequality but the price to pay for getting rid of $\tau$ in the risk bound is the increased value of error upper bound, see Appendix C in [32] for more details. Second, the constant $\mathrm{\Xi }$ associated to the Kraft-type inequality for the collection appears squared in the upper bound of the oracle inequality in Theorem 3.2. This is because of the random subcollection $\widetilde{ℳ}$ of $ℳ$, if the model collection is fixed, we get a linear bound as in Theorem 3.1.

It is important to mention that Theorems 3.1 and 3.2 provide some theoretical justification regarding the penalty shapes when using the slope heuristic for our collection of MoE models, including GLLiM, GLoME, BLLiM, BLoME models.

More precisely, our oracle inequalities show that the performance in ${\mathrm{JKL}}_{\rho }^{\otimes \mathrm{n}}$ loss of our PMLEs are roughly comparable to that of oracle models if we take large enough constants in front of the penalties, whose forms are only known up to multiplicative constants and proportional to the dimensions of models. This motivate us to make use of the slope heuristic criterion to select our data-driven hyperparameters, including the number of mixture components, the degree of polynomial Gaussian expert’s mean functions, and the potential hidden block-diagonal structures of the covariance matrices of the multivariate predictor, see more in Section 5.

To be more precise, we consider the condition (3.6) for GLoME models. As shown in the proof of Theorem 3.1, in fact we can replace the assumption on $\text{pen}(𝐦)$ by a milder one. More precisely, given a constant $ℭ$, which we specify later, there is a constant $\kappa$ depending only on $\rho$ and $C_1$, such that for every index $𝐦\in ℳ$, $\text{pen}(𝐦)$ is bounded from below by

In our numerical experiment, namely, Figures 3 and 4, we confirm the validity of the linear penalty shape assumption, which supports the use of penalties proportional to the dimension. This implies that the logarithm terms are not detected in practice as shown in Figures 3 and 4 and thus only the preponderant term in $dim\left(S_{𝐦}\right)$ is retained in the penalty form.

In particular, based on Appendix B.4 from [22], we can make explicit the dependence of the theoretical constant $\kappa$, with respect to $\rho$ and $C_1$ as follows. For instance, in Theorem 3.1, given any $\rho \in (0,1)$ and $C_1>1$, define ${ϵ}_{\text{pen}}=1-\frac{1}{C_1}$. Then $\kappa$ is determined by

where ${ϵ}_d$ is a given positive constant and

For example, if $\rho =\frac{1}{2}$, $C_1=2$, ${ϵ}_d=1$, ${\kappa }_3^{\prime }=27$, then ${ϵ}_{\text{pen}}=1-\frac{1}{C_1}=\frac{1}{2}$, ${\kappa }_0^{\prime }=3$, ${\kappa }_1^{\prime }=56+162\sqrt{2}$, ${\kappa }_2^{\prime }=84+\frac{\sqrt{3}}{2}$, and