A non-asymptotic risk bound for model selection in a high-dimensional mixture of experts via joint rank and variable selection

TrungTin Nguyen; Dung Ngoc Nguyen; Hien Duy Nguyen; Faicel Chamroukhi

Abstract

We are motivated by the problem of identifying potentially nonlinear regression relationships between high-dimensional outputs and high-dimensional inputs of heterogeneous data. This requires regression, clustering, and model selection, simultaneously. In this framework, we apply the mixture of experts models which are among the most popular ensemble learning techniques developed in the field of neural networks. In particular, we consider a more general case of mixture of experts models characterized by multiple Gaussian experts whose means are polynomials of the input variables and whose covariance matrices have block-diagonal structures. More especially, each expert is weighted by a gating network that is a softmax function of a polynomial of the input variables. These models require several hyper-parameters, including the number of mixture components, the complexity of the softmax gating networks and Gaussian mean experts, and the hidden block-diagonal structures of the covariance matrices. We provide a non-asymptotic theory for model selection of such complex hyper-parameters using the slope heuristic approach in a penalized maximum likelihood estimation framework. Specifically, we establish a non-asymptotic risk bound on the penalized maximum likelihood estimation, which takes the form of an oracle inequality, given lower bound assumptions on the penalty function.

Keywords: Dimensionality reduction Low rank estimation Mixture of experts Finite mixture regression Non-asymptotic model selection Oracle inequality Variable selection.

1 Introduction

Mixture of experts (MoE) models, introduced by Jacobs et al. [16] are widely applied to decompose the prediction model through a combination of gating models and expert models, both of which depend on the input variables. These flexible models are specific instances of conditional computation [3] , where different model experts are responsible for different regions of the input space. Thus, by applying only a subset of parameters to each example, MoE can increase model capacity while keeping training and inference costs roughly constant. For reviews on this topic, we refer to [19] ; [25] . Furthermore, they have gained popularity due to universal approximation properties in various special cases, including mixture models [30] ; [28] , mixture of regression models [15] , and fully-parameterized mixture of experts models [26] ; [27] . In high-dimensional multivariate multiple regression for heterogeneous data, we refer to outputs $𝐘 \in 𝒴 \subset ℝ^{Q}$ as target or response variables, and inputs $𝐗 \in 𝒳 \subset ℝ^{P}$ as explanatory or predictor variables, where $Q$ and $P$ are both much larger than the sample size. Additionally, hidden interactions may exist in the graphical structure between response variables. In such cases, regression, clustering, and model selection need to be performed simultaneously. Consequently, we employ MoE models to identify potential non-linear relationships between output and input variables in the high-dimensional heterogeneous data. We assume that $𝐘$ , conditional on $𝐗$ , follows a distribution with the true but unknown probability density function $s_{0} (\cdot ∣ 𝐗 = 𝐱)$ . Motivated by universal approximation theorems for MoE models, $s_{0}$ can be estimated by

	$s_{𝝍_{K}} (𝐲 ∣ 𝐱)$	$= \sum_{k = 1}^{K} g_{k} (𝐰 (𝐱)) ϕ (𝐲, 𝐯_{k} (𝐱), 𝚺_{k} (B_{k})), with$
	$g_{k} (𝐰 (𝐱))$	$= \frac{\exp (w_{k} (𝐱))}{\sum_{l = 1}^{K} \exp (w_{l} (𝐱))}, for k = 1, \dots, K .$		(1)

Here, on each cluster $k \in {1, \dots, K}$ , $g_{k}$ is called a softmax gating network corresponding to the weight functions, $𝐰 (𝐱) = (w_{1} (𝐱), \dots, w_{K} (𝐱))$ , of $𝐱$ , and $ϕ (\cdot, 𝐯_{k} (𝐱), 𝚺_{k} (B_{k}))$ is a Gaussian expert with the mean function $𝐯_{k} (𝐱)$ and covariance matrix $𝚺_{k} (B_{k})$ depending on the block-diagonal structure $B_{k}$ . We call $s_{𝝍_{K}} (𝐲 ∣ 𝐱)$ , defined as in (1), the softmax-gated block-diagonal MoE (SGaBloME) models with unknown functional parameters $𝝍_{K} = (w_{k}, 𝐯_{k}, 𝚺_{k} (B_{k}))_{k \in {1, \dots, K}}$ . Furthermore, when the weights and the means of the SGaBloME model $s_{𝝍_{K}}$ are the functions depending on polynomials of the input variables $𝐱$ which are specified, for $k \in {1, \dots, K}$ , respectively, as

	$w_{k} (𝐱)$	$= ω_{k 0} + \sum_{d = 1}^{D_{W}} 𝝎_{k d}^{T} 𝐱^{d}, with ω_{k 0} \in ℝ, 𝝎_{k d} \in ℝ^{P},$		(2)
	$𝐯_{k} (𝐱)$	$= 𝝊_{k 0} + \sum_{d = 1}^{D_{V}} 𝚼_{k d} 𝐱^{d}, with 𝝊_{k 0} \in ℝ^{Q}, 𝚼_{k d} \in ℝ^{Q \times P} .$		(3)

Here, $𝝎_{0} = (ω_{k 0})_{k \in {1, \dots, K}}$ , $𝝎 = (𝝎_{k 1}, \dots, 𝝎_{k D_{W}})_{k \in {1, \dots, K}}$ and $𝝊_{0} = (𝝊_{k 0})_{k \in {1, \dots, K}}$ , $𝚼 = (𝚼_{k 1}, \dots, 𝚼_{k D_{V}})_{k \in {1, \dots, K}}$ are $K$ -tuples of unknown coefficients with the maximum degrees $D_{W}$ and $D_{V}$ of polynomials for the weight and mean functions, respectively, and $𝐱^{d} = (x_{1}^{d}, \dots, x_{P}^{d})$ is a vector of all components of $𝐱$ with power $d$ . Then we call an SGaBloME model $s_{𝝍_{K}}$ defined by (1) with the weights and mean experts specified as in (2)-(3) the polynomial SGaBloME model.

Motivation for block-diagonal covariance matrices.

It is worth mentioning that the block-diagonal covariance matrices $𝚺 (𝐁) = (𝚺_{k} (B_{k}))_{k \in {1, \dots, K}}$ depend on the block structures $𝐁 = (B_{k})_{k \in {1, \dots, K}}$ that are the partitions of the outputs’ index set ${1, \dots, Q}$ for each cluster. This structure is not only a trade-off between the model complexity and sparsity but is also motivated by some real-world applications, where one wishes to perform prediction on data sets with heterogeneous observations and graph-structured hidden interactions between the outputs. A relevant example is the gene expression data where, subject to phenotypic response, genes interact with only a few other genes, there are small modules of correlated genes, see e.g. [14] for more details.

Motivation for polynomial regression.

To solve the high-dimensional regression problem, some authors applied SGaBloME models with certain simplifying assumptions. More specifically, Devijver [13] focused on a mixture of Gaussian linear regression models where the gating networks do not depend on the input variables. On the other hand, Chamroukhi et al. [7] considered MoE for multiple regression models with the univariate output variable, however, the weights and means are linear functions of the inputs and thus the capacity of MoE models is limited. In fact, in the context of convolutional neural networks, Chen et al. [8] have empirically found that the mixture of linear experts performs better than a single expert, but is still significantly worse than the mixture of non-linear experts. Within this framework, we are motivated to integrate nonlinearities into SGaBloME models by defining the weights and mean experts as linear combinations of bounded functions (LinBo) whose coefficients belong to a compact set. Such a general setting may include the polynomial basis with a bounded input domain, the suitable re-normalized wavelet dictionaries, or the Fourier basis on an interval. If the dimensions of the inputs and outputs are not too large, it is not necessary to select relevant variables and/or use rank sparse models. Then we can work on the softmax-gated MoE models with the linear combinations of bounded functions for weight and mean functions as in [24] . However, to deal with high-dimensional data and simplify the interpretation of sparsity, we consider a special case of LinBo-SGaBloME models to explore the presence of nonlinearities that is the class of polynomial SGaBloME models defined by (1)-(3). On the convergence rates of polynomial SGaBloME models, we refer to [23] for a discussion of the optimal convergence rate of an MoE model where each expert is associated with a polynomial regression model.

Model selection for polynomial SGaBloME models.

The estimation of SGaBloME models can be performed by using a well-known expectation-maximization (EM) algorithm [11] , which obtains the global convergence in regression mixture models [18] . However, it crucially requires data-driven hyper-parameter choices, including the number of mixture components, the degree of complexity of each softmax gating network and each Gaussian expert mean function, and the hidden block-diagonal structures of the covariance matrices. Hyper-parameter choices from the data-driven learning algorithms belong to the class of model selection problems that select the model with the lowest risk from the data. Typically, penalization is one of the main strategies proposed for model selection that minimizes the sum of the empirical risk with a term of penalty so that the model can be fitted to data while avoiding the overfitting problem.

Related works.

Typically, model selection for MoE models is performed using the asymptotic criteria [31] ; [4] , whose uses in small samples are limited. Birgé et al. [5] proposed a novel approach, called slope heuristic, supported by a non-asymptotic oracle inequality via a general model selection theorem, see [2] and the references therein for recent reviews. This method leads to an optimal data-driven choice of multiplicative constants for penalties. In fact, oracle inequalities for the least absolute shrinkage and selection operator (Lasso) [32] and general penalized maximum likelihood estimators were established in the spirit of the methods based on concentration inequalities developed by [20] . These results include work on the simplified assumptions of MoE models such as high-dimensional Gaussian graphical models [14] , Gaussian mixture model selection [21] , and finite mixture regression models [12] ; [13] , or softmax-gated MoE models with linear combinations of bounded functions for weight and mean functions without consideration of variable selection in the high-dimensional setting [24] .

1.1 Main contributions

In this work, we established an oracle inequality for model selection, as shown in Theorem 1, under lower bound assumptions on penalty terms. This allows us to obtain non-asymptotic risk bounds in the form of weak oracle inequalities allowing the numbers of predictor and response variables that grow or are even much larger than the sample size. More concretely, the constructed oracle inequality shows that the performance of our penalized maximum likelihood estimations is comparable to that of oracle models with sufficiently large constant multiples of the penalties. The forms of these constants are only known up to multiplicative constants and are proportional to the dimensions of the models. Moreover, the flexibility of polynomial SGaBloME models requires the hyper-parameters comprising 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 output. Therefore, the aforementioned theoretical justifications for the penalty shapes motivate the use of the heuristic slope criterion to select these hyper-parameters of the models under consideration.

1.1.1 Notations.

For any matrix $𝐀$ with the elements $A_{i j}$ , we denote $| | | 𝐀 | | |_{\infty} = \max_{i, j} | A_{i j} |$ the max-norm, and $𝐀_{\cdot j}$ the $j^{t h}$ column, of $𝐀$ . Furthermore, the smallest and largest eigenvalues of $𝐀$ are denoted by $eig (𝐀)$ and $Eig (𝐀)$ , respectively. The notation $𝐀 ≻ 0$ indicates that $𝐀$ is positive definite. For any vector $𝐚$ , we denote $∥ 𝐚 ∥_{p}$ the $l_{p}$ -norm of $𝐚$ for $0 < p \leq \infty$ . We call $dim (𝒮)$ the total number of parameters to be estimated or the dimension of a parametric model $𝒮$ . If $S$ is a finite set, we denote $c a r d (S)$ the cardinality, $𝒫 (S)$ the set of all subsets, and $ℬ (S)$ the set of all partitions, of $S$ . The set of all natural numbers without zero is denoted by $ℕ^{*}$ . For any $K \in ℕ^{*}$ , the notation $[K]$ corresponds to the set ${1, \dots, K}$ . Finally, we refer to $a \land b$ as $\min {a, b}$ for $a, b \in ℝ$ .

1.1.2 Paper organization.

The rest of the paper is organized as follows. Section 2 is devoted to the construction of a collection of polynomial SGaBloME models for high-dimensional heterogeneous data. In Section 3, we state the main theoretical results of oracle inequality for the penalized maximum likelihood estimations under some conditions on the parameter space and input domain of the models. Finally, Section 4 contains concluding remarks and future directions.

2 Collection of polynomial SGaBloME models

For high-dimensional data, it is necessary to work with parsimonious models by combining two well-known approaches: selection of relevant variables and rank sparse models. Within this framework, the collection of polynomial SGaBloME models is then constructed.

2.1 Variable selection via selecting relevant variables

In this section, we introduce the index sets for the input and output variables so that they are related to each other. This facilitates the variable selection of the models in a highly dimensional framework. In particular, for every $p \in [P], q \in [Q]$ , we call a couple $(X_{p}, Y_{q})$ irrelevant if the elements $(𝚼_{k d})_{q p} = 0$ and $(𝝎_{k l})_{p} = 0$ for all cluster $k \in [K]$ and degrees $d \in [D_{V}]$ , $l \in [D_{W}]$ . Therefore, the variables $(X_{p}, Y_{q})$ are relevant if they are not irrelevant. Formally, we denote $I = {(p, q) \in [P] \times [Q] : (X_{p}, Y_{q}) is irrelevant}$ the set of indices of irrelevant couples, and the complement of $I$ , called $J = ([P] \times [Q]) ∖ I$ , is thus the set of indices of relevant couples with $J \in 𝒫 ([P] \times [Q])$ . In addition, we also denote $J_{i n} = {p \in [P] : \exists q \in [Q], (p, q) \in J}$ the set of indices of input variables that are relevant to the outputs so that $J_{i n} \subseteq [P]$ .

We notice that, for every cluster $k \in [K]$ and degree $d \in [D_{V}]$ , all entries of $𝚼_{k d}$ belonging to columns indexed by $[P] ∖ J_{i n}$ equal to $0$ , in other words, $𝚼_{k d}$ has the relevant columns indexed by $J_{i n}$ . Hence, the matrix $𝚼_{k d}$ will have $Q \times c a r d (J_{i n})$ coefficients to be estimated, which are smaller than $Q \times P$ when all variables are considered. The number of parameters in the regression matrices is therefore considerably reduced when the cardinality of $J_{i n}$ is much smaller than the number of input variables $P$ . The subsets $J$ or $J_{i n}$ can be constructed by the Lasso [32] and has been extended to deal with multiple multivariate regression models for column sparsity using the Group-Lasso [33] .

2.2 Variable selection via rank sparse models

Anderson et al. [1] introduced rank sparse models in the regression framework which is if regression matrices have low rank or at least can be well approximated by low-rank matrices, then the corresponding regression models are said to be rank sparse. In the polynomial SGaBloME models, we assume that for every cluster $k \in [K]$ and degree $d \in [D_{V}]$ , the matrix $𝚼_{k d}$ has the associated rank $R_{k d}$ and therefore it is completely determined by $R_{k d} \times (P - (Q - R_{k d}))$ coefficients, which can be less than $Q \times P$ . Combined with the selection of relevant variables method, we denote a rank matrix by $𝐑 = (R_{k d})_{k \in [K], d \in [D_{V}]}$ with the element $R_{k d} \in [c a r d (J_{i n}) \land Q]$ for each $k \in [K], d \in [D_{V}]$ .

2.3 Collection of polynomial SGaBloME models

So far, each polynomial SGaBloME model defined by (1)-(3) can be characterized by the set $𝐦 = (K, D_{W}, D_{V}, 𝐁, J, 𝐑)$ where $K$ is the number of clusters, $D_{W}$ and $D_{V}$ are the maximum degrees of polynomials of the weight and mean functions, respectively, $𝐁$ is the set of the block-diagonal structures of the covariance matrices, $J$ is the set of relevant variables, and $𝐑$ is the rank matrix of coefficient matrices. Let $𝒮_{𝐦}$ be a class of (conditional) densities of polynomial SGaBloME models with respect to $𝐦$ , which is specified as

$𝒮_{𝐦}$	$= {s_{𝝍_{𝐦}} \equiv s_{𝝍_{K}} defined by (1)-(3) with$
	$𝝍_{𝐦} = (𝝎_{0}, 𝝎, 𝝊_{0}, 𝚼, 𝚺 (𝐁)) \in 𝚿_{𝐦},$
	$𝚿_{𝐦} = ℝ^{K} \times 𝐖_{J}^{K \times D_{W}} \times ℝ^{K \times Q} \times 𝐕_{J, 𝐑}^{K \times D_{V}} \times 𝛀_{𝐁}^{K},$	(4)
	$𝐖_{J} = {𝜶 = (α_{1}, \dots, α_{P}) \in ℝ^{P} : α_{j} = 0, \forall j \in ([P] ∖ J_{i n})},$
	$𝐕_{J, 𝐑} = {𝐀 \in ℝ^{Q \times P} : 𝐀_{\cdot j} = 𝟎, \forall j \in [P] ∖ J_{i n} and rank(𝐀) = R \in 𝐑}$
	$𝛀_{𝐁} = {𝚺 (B) \in ℝ^{Q \times Q} : 𝚺 (B) ≻ 0 and 𝚺 (B) depends on block B \in 𝐁}},$

for every $𝐦 \in ℕ^{*} \times ℕ^{*} \times ℕ^{*} \times ℬ ([Q])^{K} \times 𝒫 ([P] \times [Q]) \times [c a r d (J_{i n}) \land Q]^{K \times D_{V}}$ . The collection of polynomial SGaBloME models defined in (4) is generally large and therefore not feasible in practice. Therefore, we restrict the set of $(K, D_{W}, D_{V})$ to a finite set of $(𝒦, 𝒟_{W}, 𝒟_{V})$ where, for $K^{*}, D_{W}^{*}, D_{V}^{*} \in ℕ^{*}$ , $𝒦 = [K^{*}], 𝒟_{W} = [D_{W}^{*}], 𝒟_{V} = [D_{V}^{*}]$ . Accordingly, the collection of polynomial SGaBloME models on the deterministic set of hyper-parameters can be defined as

	$𝒮$	$= {𝒮_{𝐦} defined by (4) such that 𝐦 \in ℳ,$		(5)
		$ℳ = 𝒦 \times 𝒟_{W} \times 𝒟_{V} \times ℬ ([Q])^{K} \times 𝒫 ([P] \times [Q]) \times [c a r d (J_{i n}) \land Q]^{K \times D_{V}}} .$

Furthermore, because the block structures are specified by the partitions of the index set ${1, \dots, Q}$ , the number of such structures follows the so-called Bell number, which grows exponentially even for a moderate number of variables $Q$ and clusters $K$ . Therefore, it is infeasible to consider an exhaustive exploration of the combination of all the partitions to detect the block structures for covariance matrices. Motivated by the recent work of [14] , for a set of thresholds $ℰ$ and on each cluster $k \in [K]$ , we restrict our attention to the sub-collection $ℬ_{k, ℰ} = (ℬ_{k, ϵ})_{ϵ \in ℰ}$ of $ℬ ([Q])$ , where $ℬ_{k, ϵ}$ is the partition of the output variables corresponding to the block-diagonal structure of the adjacency matrix $𝐄_{k, ϵ}$ based on the thresholded absolute values of the sample covariance matrix $𝐒_{k}$ . More formally, for each $ϵ \in ℰ$ , $(𝐄_{k, ϵ})_{q q^{'}} = 1$ if $| (𝐒_{k})_{q q^{'}} | > ϵ$ , otherwise it is equal to $0$ for $q, q^{'} \in [Q]$ . In fact, Mazumder et al. [22] have shown that the class of block-diagonal structures detected by the graphical Lasso algorithm is identical to the block-diagonal structures detected by the thresholding of the sample covariance matrices, which supports our motivation for this restriction.

For the set of relevant variables, we focus on a random subset $𝒥$ of $𝒫 ([P] \times [Q])$ with the controlled size of $𝒥$ required in the high-dimension case. Accordingly, the number of possible vectors of ranks is reduced by working on a random subset of $[c a r d (J_{i n}) \land Q]^{K \times D_{V}}$ , which is denoted by $ℛ_{(K, J, D_{V})}$ depending on $J \in 𝒥$ with the dimension of $K \in 𝒦$ and $D_{V} \in 𝒟_{V}$ . As a result, the collection of polynomial SGaBloME models based on a random sub-collection of hyper-parameters can be specified as

	$\tilde{𝒮}$	$= {𝒮_{𝐦} defined by (4) such that 𝐦 \in \tilde{ℳ},$		(6)
		$\tilde{ℳ} = 𝒦 \times 𝒟_{W} \times 𝒟_{V} \times (ℬ_{k, ℰ})_{k \in [K]} \times 𝒥 \times ℛ_{(K, J, D_{V})}} .$

3 Main theoretical results

In this section, we begin by introducing conditions on the parameter space of the models and give an overview of loss functions that are useful for comparing two (conditional) probability density functions. A general principle of penalized maximum likelihood estimation is also derived. Next, we show a finite sample oracle inequality used to ensure that if we penalize the log-likelihood in an approximate approach, we are able to select a model that is as good as the oracle.

3.1 Boundedness conditions on the parameter space

By motivation of integrating the nonlinearities into SGaBloME models discussed in Section 1, we consider the class of linear combinations of bounded functions for the weights and mean experts whose coefficients belong to compact sets with a bounded input domain. More specifically, we let $(𝐗_{[N]}, 𝐘_{[N]}) = ((𝐗_{1}, 𝐘_{1}), \dots, (𝐗_{N}, 𝐘_{N}))$ be $N$ pairs of real-valued random variables $(𝐗, 𝐘)$ where the covariates $𝐗$ are assumed to belong to a hypercube, that is $𝒳 = [0, 1]^{P}$ . Then, there exist the constants $C_{𝝎}, C_{𝚼}, c_{𝚺}, C_{𝚺} > 0$ such that, for every $k \in [K]$ ,

∥ 𝝎_{k d} ∥_{\infty} \leq C_{𝝎}, | | | 𝚼_{k l} | | |_{\infty} \leq C_{𝚼}, for every d \in [D_{W}], l \in [D_{V}],

(7)

moreover, the eigenvalues of the block-diagonal covariances of the Gaussian experts lie on a positive interval, that is

0 < c_{𝚺} \leq eig (𝚺_{k} (B_{k})) \leq Eig (𝚺_{k} (B_{k})) \leq C_{𝚺} .

(8)

This setting can be applied to the case of polynomial functions for the weights of the softmax gates and the means of the Gaussian experts as we described in (2)-(3). More generally, the oracle inequality provided by Theorem 1 still holds for monomials of weights, allowing for the interaction between different inputs.

3.2 Loss function

To evaluate the maximum likelihood estimate, the Kullback-Leibler (KL) divergence is the most natural loss function, which is generally defined by

K L (s, t) = {\begin{matrix} \int_{ℝ^{D}} \ln (\frac{s (x)}{t (x)}) s (x) d y x & if s d x is absolutely continuous w.r.t. t d x, \\ + \infty & otherwise, \end{matrix}

where $s (\cdot)$ and $t (\cdot)$ are two density functions. In our work, we will apply the tensorized KL divergence to capture the structure of the density functions conditional on the random variables $𝐗$ , that is

K L^{\otimes N} (s, t) = 𝔼_{𝐗_{[N]}} [\frac{1}{N} \sum_{n = 1}^{N} K L (s (\cdot ∣ 𝐗_{n}), t (\cdot ∣ 𝐗_{n}))] .

Another case of the tensorized KL divergence is the tensorized Jensen-KL divergence [9] , which is given, for any $ρ \in (0, 1)$ , by

J K L_{ρ}^{\otimes N} (s, t) = 𝔼_{𝐗_{[N]}} [\frac{1}{N} \sum_{n = 1}^{N} \frac{1}{ρ} K L (s (\cdot | 𝐗_{n}), (1 - ρ) s (\cdot ∣ 𝐗_{n}) + ρ t (\cdot ∣ 𝐗_{n}))] .

A relationship between the tensorized KL and the tensorized Jensen-KL divergence can be found in (10, Proposition 1).

3.3 Penalized maximum likelihood estimation (PMLE)

In the context of maximum likelihood estimation, given a collection $S_{𝐦}$ , we aim to estimate $s_{0}$ by the conditional density $\hat{s}_{𝐦}$ that minimizes the negative log-likelihood (NLL) as

\hat{s}_{𝐦} = a r g m i n_{s_{𝐦} \in S_{𝐦}} \sum_{n = 1}^{N} - \ln [s_{𝐦} (𝐘_{n} ∣ 𝐗_{n})] .

It is important to us to look for almost minimizer of this quantity and thereby define an $η$ -log-likelihood minimizer (LLM) that satisfies

\sum_{n = 1}^{N} - \ln [\hat{s}_{𝐦} (𝐘_{n} ∣ 𝐗_{n})] \leq i n f_{s_{𝐦} \in 𝒮_{𝐦}} \sum_{n = 1}^{N} - \ln [s_{𝐦} (𝐘_{n} ∣ 𝐗_{n})] + η,

(9)

where the error term $η > 0$ is added to avoid any existence issue such as the infimum may not be reached. See (6, Chapter 2), [10] ; [9] ; [24] for more details of this literature. However, this approach underestimates the risk of the estimation and leads to the selection of overly complex models. Therefore, a trade-off between good data fit and model complexity can be found by adding an appropriate penalty term $pen (𝐦)$ . More concretely, for a given choice of $pen (𝐦)$ , the selected model $𝒮_{\hat{𝐦}}$ is chosen as the one whose index $\hat{𝐦}$ is a $η^{'}$ -minimizer of the sum of the NLL and penalty function, that is

\sum_{n = 1}^{N} - \ln [\hat{s}_{\hat{𝐦}} (𝐘_{n} | 𝐗_{n})] + pen (\hat{𝐦}) \leq i n f_{𝐦 \in ℳ} {\sum_{n = 1}^{N} - \ln [\hat{s}_{𝐦} (𝐘_{n} | 𝐗_{n})] + pen (𝐦)} + η^{'},

(10)

for $η^{'} > 0$ . We then call $\hat{s}_{\hat{𝐦}}$ the $η^{'}$ -PMLE that depends on both error terms $η$ and $η^{'}$ . From now on, the term selected model or best data-driven model is used to indicate the model that satisfies (10).

3.4 Oracle inequality

In this section, we provide the construction of an oracle inequality that guarantees a non-asymptotic theory for model selection in high-dimensional polynomial SGaBloME models.

Theorem 1.

Let $(𝐗_{[N]}, 𝐘_{[N]})$ be a random sample of $(𝐗, 𝐘)$ where $𝐘 | 𝐗$ arises from the unknown conditional density $s_{0}$ . For every $𝐦 = (K, D_{W}, D_{V}, 𝐁, J, 𝐑) \in ℳ$ , the model $𝒮_{𝐦}$ can be specified by (4). Assume that there exists $τ > 0$ and $ϵ_{K L} > 0$ such that, for all $𝐦 \in ℳ$ , one can find $\bar{s}_{𝐦} \in 𝒮_{𝐦}$ such that

K L^{\otimes N} (s_{0}, \bar{s}_{𝐦}) \leq i n f_{s \in 𝒮_{𝐦}} K L^{\otimes N} (s_{0}, s) + \frac{ϵ_{K L}}{N}, and \bar{s}_{𝐦} \geq e^{- τ} s_{0} .

(11)

Furthermore, we construct a random sub-collection $\tilde{𝒮}$ of $𝒮$ such that every model of $\tilde{𝒮}$ depends on the sets of $\tilde{ℳ} \subset ℳ$ as in (5)-(6). Then, there is a constant $C$ such that, for any $ρ \in (0, 1)$ and $C_{1} > 1$ , there are two constants $κ$ and $C_{2}$ depending only on $ρ$ and $C_{1}$ such that, for every $𝐦 \in ℳ$ , $ξ_{𝐦} \in ℝ^{+}$ , $Ξ = \sum_{𝐦 \in ℳ} e^{- ξ_{𝐦}} < \infty$ and

pen (𝐦) \geq κ [(C + \ln N) dim (𝒮_{𝐦}) + (1 \lor τ) ξ_{𝐦}],

and the $η^{'}$ -PMLE $\hat{s}_{\hat{𝐦}}$ defined in (10) on the subset $\tilde{ℳ} \subset ℳ$ satisfies

	$𝔼_{𝐗_{[N]}, 𝐘_{[N]}} [J K L_{ρ}^{\otimes N} (s_{0}, \hat{s}_{\hat{𝐦}})]$	$\leq C_{1} 𝔼_{𝐗_{[N]}, 𝐘_{[N]}} [i n f_{𝐦 \in \tilde{ℳ}} (i n f_{s \in 𝒮_{𝐦}} K L^{\otimes N} (s_{0}, s) + 2 \frac{pen (𝐦)}{N})]$
		$+ C_{2} (1 \lor τ) \frac{Ξ^{2}}{N} + \frac{η^{'} + η}{N} .$		(12)

Remarks.

Theorem 1 guarantees that a penalized criterion leads to a good model selection and that the penalty is only known up to multiplicative constants $κ$ , and is proportional to the dimensions of the models $dim (𝒮_{𝐦})$ . In particular, in the small and finite sample setting, these multiplicative constants can be calibrated using the slope heuristic approach. We notice that (11) is not a strong assumption and is satisfied in the case $s_{0}$ is bounded with compact support. This oracle inequality compares the performance of our PMLE with the best model in the collection. However, Theorem 1 allows us to approximate well a rich class of conditional PDFs if we use polynomials of weights and Gaussian expert means of sufficient degrees, or enough clusters due to the universal approximation of MoE models. This results in the term on the right of (1) being small, for $𝒟_{W}, 𝒟_{V}$ and $𝒦$ well chosen. It should be emphasized that Theorem 1 extends the main result of [24] , which is only valid for a full collection of LinBoSGaME models in the low-dimensional setting. Furthermore, in the context of MoE models, our non-asymptotic oracle inequality for SGaME models in Theorem 1 can be seen as a complementary result to a classical asymptotic theory (17, Theorems 1, 2, and 3), and an $l_{1}$ oracle inequality that focuses on the properties of the Lasso estimator rather than the model selection procedure [29] .

Main challenges on the proof of Theorem 1.

To prove Theorem 1 it is inspired by [24] for handling LinBo-SGaME models, however, our method and most of the technical details differ. This is because their approach is not directly applicable to our high-dimensional SGaBloME models, due to restrictions on relevant predictor variables and rank reduction, and Gaussian experts with block-diagonal covariance matrices. In particular, the main difficulty in proving our oracle inequality lies in bounding the bracketing entropy for the collections of SGaBloME models. This requires several regularity assumptions, which are not easy to verify due to the complexity of SGaBloME models and technical reasons. Therefore, our proofs require the development of several new ideas. Furthermore, unlike [24] , which uses a model selection theorem for a deterministic collection of models from [10] ; [9] , we need to find a way to use the model selection theorem for MLE among a random sub-collection (cf. (12, Theorem 5.1) and (14, Theorem 7.3)). We refer readers to Sections S-1 and S-2 in the supplementary materials for a sketch of proof and detailed proof of Theorem 1.

4 Conclusion and perspectives

We have studied PMLEs for polynomial SGaBloME models in high-dimensional heterogeneous data. Our main contribution is to establish a non-asymptotic risk bound in the form of an oracle inequality, provided that lower bounds of the penalty hold. The future direction is to empirically evaluate our oracle inequality and to extend the current oracle inequality to more general settings where Gaussian experts are replaced by elliptic distributions.

References

(1) Anderson, C.W., Stolz, E.A., Shamsunder, S.: Multivariate autoregressive models for classification of spontaneous electroencephalographic signals during mental tasks. IEEE Transactions on Biomedical Engineering 45(3), 277–286 (1998) ↩ 2.2
(2) Arlot, S.: Minimal penalties and the slope heuristics: a survey. Journal de la Société Française de Statistique 160(3), 1–106 (2019) ↩ 1
(3) Bengio, Y.: Deep Learning of Representations: Looking Forward. In: Statistical Language and Speech Processing. pp. 1–37. Springer Berlin Heidelberg, Berlin, Heidelberg (2013) ↩ 1
(4) Biernacki, C., Celeux, G., Govaert, G.: Assessing a mixture model for clustering with the integrated completed likelihood. IEEE Transactions on Pattern Analysis and Machine Intelligence 22(7), 719–725 (2000) ↩ 1
(5) Birgé, L., Massart, P.: Minimal penalties for Gaussian model selection. Probability Theory and Related Fields 138(1), 33–73 (2007) ↩ 1
(6) Borwein, J.M., Zhu, Q.J.: Techniques of Variational Analysis. Springer (2004) ↩ 3.3
(7) Chamroukhi, F., Huynh, B.T.: Regularized maximum-likelihood estimation of mixture-of-experts for regression and clustering. In: 2018 International Joint Conference on Neural Networks (IJCNN). pp. 1–8 (2018) ↩ 1
(8) Chen, Z., Deng, Y., Wu, Y., Gu, Q., Li, Y.: Towards Understanding the Mixture-of-Experts Layer in Deep Learning. In: NeurIPS (2022) ↩ 1
(9) Cohen, S.X., Le Pennec, E.: Partition-based conditional density estimation. ESAIM: Probability and Statistics 17, 672–697 (2013) ↩ 3.2 3.3 3.4
(10) Cohen, S., Le Pennec, E.: Conditional density estimation by penalized likelihood model selection and applications. Technical report, INRIA (2011) ↩ 3.2 3.3 3.4
(11) Dempster, A.P., Laird, N.M., Rubin, D.B.: Maximum Likelihood from Incomplete Data via the EM Algorithm. Journal of the Royal Statistical Society. Series B 39(1), 1–38 (1977) ↩ 1
(12) Devijver, E.: Finite mixture regression: a sparse variable selection by model selection for clustering. Electronic Journal of Statistics 9(2), 2642–2674 (2015) ↩ 1 3.4
(13) Devijver, E.: Joint rank and variable selection for parsimonious estimation in a high-dimensional finite mixture regression model. Journal of Multivariate Analysis 157, 1–13 (2017) ↩ 1
(14) Devijver, E., Gallopin, M.: Block-diagonal covariance selection for high-dimensional Gaussian graphical models. Journal of the American Statistical Association 113(521), 306–314 (2018) ↩ 1 2.3 3.4
(15) Ho, N., Yang, C.Y., Jordan, M.I.: Convergence Rates for Gaussian Mixtures of Experts. Journal of Machine Learning Research 23(323), 1–81 (2022) ↩ 1
(16) Jacobs, R.A., Jordan, M.I., Nowlan, S.J., Hinton, G.E.: Adaptive mixtures of local experts. Neural Computation 3(1), 79–87 (1991) ↩ 1
(17) Khalili, A.: New estimation and feature selection methods in mixture-of-experts models. Canadian Journal of Statistics 38(4), 519–539 (2010) ↩ 3.4
(18) Kwon, J., Qian, W., Caramanis, C., Chen, Y., Davis, D.: Global convergence of the EM algorithm for mixtures of two component linear regression. In: COLT. vol. 99, pp. 2055–2110. PMLR (2019) ↩ 1
(19) Masoudnia, S., Ebrahimpour, R.: Mixture of experts: a literature survey. Artificial Intelligence Review 42(2), 275–293 (2014) ↩ 1
(20) Massart, P.: Concentration Inequalities and Model Selection: Ecole d’Eté de Probabilités de Saint-Flour XXXIII-2003. Springer (2007) ↩ 1
(21) Maugis, C., Michel, B.: A non asymptotic penalized criterion for Gaussian mixture model selection. ESAIM: Probability and Statistics 15, 41–68 (2011) ↩ 1
(22) Mazumder, R., Hastie, T.: Exact covariance thresholding into connected components for large-scale graphical lasso. The Journal of Machine Learning Research 13(1), 781–794 (2012) ↩ 2.3
(23) Mendes, E.F., Jiang, W.: On convergence rates of mixtures of polynomial experts. Neural Computation 24(11), 3025–3051 (2012) ↩ 1
(24) Montuelle, L., Le Pennec, E., et al.: Mixture of Gaussian regressions model with logistic weights, a penalized maximum likelihood approach. Electronic Journal of Statistics 8(1), 1661–1695 (2014) ↩ 1 3.3 3.4
(25) Nguyen, H.D., Chamroukhi, F.: Practical and theoretical aspects of mixture-of-experts modeling: An overview. Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery 8(4), e1246 (2018) ↩ 1
(26) Nguyen, H.D., Chamroukhi, F., Forbes, F.: Approximation results regarding the multiple-output Gaussian gated mixture of linear experts model. Neurocomputing 366, 208–214 (2019) ↩ 1
(27) Nguyen, H.D., Nguyen, T., Chamroukhi, F., McLachlan, G.J.: Approximations of conditional probability density functions in Lebesgue spaces via mixture of experts models. Journal of Statistical Distributions and Applications 8(1), 13 (2021) ↩ 1
(28) Nguyen, T., Chamroukhi, F., Nguyen, H.D., McLachlan, G.J.: Approximation of probability density functions via location-scale finite mixtures in Lebesgue spaces. Communications in Statistics - Theory and Methods pp. 1–12 (2022) ↩ 1
(29) Nguyen, T., Nguyen, H.D., Chamroukhi, F., McLachlan, G.J.: An $l_{1}$ -oracle inequality for the Lasso in mixture-of-experts regression models. arXiv preprint arXiv:2009.10622 (2020) ↩ 3.4
(30) Nguyen, T., Nguyen, H.D., Chamroukhi, F., McLachlan, G.J.: Approximation by finite mixtures of continuous density functions that vanish at infinity. Cogent Mathematics & Statistics 7(1), 1750861 (2020) ↩ 1
(31) Schwarz, G., et al.: Estimating the dimension of a model. The Annals of Statistics 6(2), 461–464 (1978) ↩ 1
(32) Tibshirani, R.: Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B 58(1), 267–288 (1996) ↩ 1 2.1
(33) Yuan, M., Lin, Y.: Model selection and estimation in regression with grouped variables. Journal of the Royal Statistical Society: Series B 68(1), 49–67 (2006) ↩ 2.1

[bib.bib1] (1) Anderson, C.W., Stolz, E.A., Shamsunder, S.: Multivariate autoregressive models for classification of spontaneous electroencephalographic signals during mental tasks. IEEE Transactions on Biomedical Engineering 45(3), 277–286 (1998) ↩ 2.2

[bib.bib2] (2) Arlot, S.: Minimal penalties and the slope heuristics: a survey. Journal de la Société Française de Statistique 160(3), 1–106 (2019) ↩ 1

[bib.bib3] (3) Bengio, Y.: Deep Learning of Representations: Looking Forward. In: Statistical Language and Speech Processing. pp. 1–37. Springer Berlin Heidelberg, Berlin, Heidelberg (2013) ↩ 1

[bib.bib4] (4) Biernacki, C., Celeux, G., Govaert, G.: Assessing a mixture model for clustering with the integrated completed likelihood. IEEE Transactions on Pattern Analysis and Machine Intelligence 22(7), 719–725 (2000) ↩ 1

[bib.bib5] (5) Birgé, L., Massart, P.: Minimal penalties for Gaussian model selection. Probability Theory and Related Fields 138(1), 33–73 (2007) ↩ 1

[bib.bib6] (6) Borwein, J.M., Zhu, Q.J.: Techniques of Variational Analysis. Springer (2004) ↩ 3.3

[bib.bib7] (7) Chamroukhi, F., Huynh, B.T.: Regularized maximum-likelihood estimation of mixture-of-experts for regression and clustering. In: 2018 International Joint Conference on Neural Networks (IJCNN). pp. 1–8 (2018) ↩ 1

[bib.bib8] (8) Chen, Z., Deng, Y., Wu, Y., Gu, Q., Li, Y.: Towards Understanding the Mixture-of-Experts Layer in Deep Learning. In: NeurIPS (2022) ↩ 1

[bib.bib9] (9) Cohen, S.X., Le Pennec, E.: Partition-based conditional density estimation. ESAIM: Probability and Statistics 17, 672–697 (2013) ↩ 3.2 3.3 3.4

[bib.bib10] (10) Cohen, S., Le Pennec, E.: Conditional density estimation by penalized likelihood model selection and applications. Technical report, INRIA (2011) ↩ 3.2 3.3 3.4

[bib.bib11] (11) Dempster, A.P., Laird, N.M., Rubin, D.B.: Maximum Likelihood from Incomplete Data via the EM Algorithm. Journal of the Royal Statistical Society. Series B 39(1), 1–38 (1977) ↩ 1

[bib.bib12] (12) Devijver, E.: Finite mixture regression: a sparse variable selection by model selection for clustering. Electronic Journal of Statistics 9(2), 2642–2674 (2015) ↩ 1 3.4

[bib.bib13] (13) Devijver, E.: Joint rank and variable selection for parsimonious estimation in a high-dimensional finite mixture regression model. Journal of Multivariate Analysis 157, 1–13 (2017) ↩ 1

[bib.bib14] (14) Devijver, E., Gallopin, M.: Block-diagonal covariance selection for high-dimensional Gaussian graphical models. Journal of the American Statistical Association 113(521), 306–314 (2018) ↩ 1 2.3 3.4

[bib.bib15] (15) Ho, N., Yang, C.Y., Jordan, M.I.: Convergence Rates for Gaussian Mixtures of Experts. Journal of Machine Learning Research 23(323), 1–81 (2022) ↩ 1

[bib.bib16] (16) Jacobs, R.A., Jordan, M.I., Nowlan, S.J., Hinton, G.E.: Adaptive mixtures of local experts. Neural Computation 3(1), 79–87 (1991) ↩ 1

[bib.bib17] (17) Khalili, A.: New estimation and feature selection methods in mixture-of-experts models. Canadian Journal of Statistics 38(4), 519–539 (2010) ↩ 3.4

[bib.bib18] (18) Kwon, J., Qian, W., Caramanis, C., Chen, Y., Davis, D.: Global convergence of the EM algorithm for mixtures of two component linear regression. In: COLT. vol. 99, pp. 2055–2110. PMLR (2019) ↩ 1

[bib.bib19] (19) Masoudnia, S., Ebrahimpour, R.: Mixture of experts: a literature survey. Artificial Intelligence Review 42(2), 275–293 (2014) ↩ 1

[bib.bib20] (20) Massart, P.: Concentration Inequalities and Model Selection: Ecole d’Eté de Probabilités de Saint-Flour XXXIII-2003. Springer (2007) ↩ 1

[bib.bib21] (21) Maugis, C., Michel, B.: A non asymptotic penalized criterion for Gaussian mixture model selection. ESAIM: Probability and Statistics 15, 41–68 (2011) ↩ 1

[bib.bib22] (22) Mazumder, R., Hastie, T.: Exact covariance thresholding into connected components for large-scale graphical lasso. The Journal of Machine Learning Research 13(1), 781–794 (2012) ↩ 2.3

[bib.bib23] (23) Mendes, E.F., Jiang, W.: On convergence rates of mixtures of polynomial experts. Neural Computation 24(11), 3025–3051 (2012) ↩ 1

[bib.bib24] (24) Montuelle, L., Le Pennec, E., et al.: Mixture of Gaussian regressions model with logistic weights, a penalized maximum likelihood approach. Electronic Journal of Statistics 8(1), 1661–1695 (2014) ↩ 1 3.3 3.4

[bib.bib25] (25) Nguyen, H.D., Chamroukhi, F.: Practical and theoretical aspects of mixture-of-experts modeling: An overview. Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery 8(4), e1246 (2018) ↩ 1

[bib.bib26] (26) Nguyen, H.D., Chamroukhi, F., Forbes, F.: Approximation results regarding the multiple-output Gaussian gated mixture of linear experts model. Neurocomputing 366, 208–214 (2019) ↩ 1

[bib.bib27] (27) Nguyen, H.D., Nguyen, T., Chamroukhi, F., McLachlan, G.J.: Approximations of conditional probability density functions in Lebesgue spaces via mixture of experts models. Journal of Statistical Distributions and Applications 8(1), 13 (2021) ↩ 1

[bib.bib28] (28) Nguyen, T., Chamroukhi, F., Nguyen, H.D., McLachlan, G.J.: Approximation of probability density functions via location-scale finite mixtures in Lebesgue spaces. Communications in Statistics - Theory and Methods pp. 1–12 (2022) ↩ 1

[bib.bib29] (29) Nguyen, T., Nguyen, H.D., Chamroukhi, F., McLachlan, G.J.: An $l_{1}$ -oracle inequality for the Lasso in mixture-of-experts regression models. arXiv preprint arXiv:2009.10622 (2020) ↩ 3.4

[bib.bib30] (30) Nguyen, T., Nguyen, H.D., Chamroukhi, F., McLachlan, G.J.: Approximation by finite mixtures of continuous density functions that vanish at infinity. Cogent Mathematics & Statistics 7(1), 1750861 (2020) ↩ 1

[bib.bib31] (31) Schwarz, G., et al.: Estimating the dimension of a model. The Annals of Statistics 6(2), 461–464 (1978) ↩ 1

[bib.bib32] (32) Tibshirani, R.: Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B 58(1), 267–288 (1996) ↩ 1 2.1

[bib.bib33] (33) Yuan, M., Lin, Y.: Model selection and estimation in regression with grouped variables. Journal of the Royal Statistical Society: Series B 68(1), 49–67 (2006) ↩ 2.1