Conceptio › Archive › arXiv CS

arXiv CSopen access

On Bayesian Softmax-Gated Mixture-of-Experts Models

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

neural-networks

machine learning, deep learning, neural networks

On Bayesian Softmax-Gated Mixture-of-Experts Models Nicola Bariletto1

Huy Nguyen1

Nhat Ho

Alessandro Rinaldo

The University of Texas at Austin April 23, 2026

arXiv:2604.20551v1 [stat.ML] 22 Apr 2026

Abstract Mixture-of-experts models provide a flexible framework for learning complex probabilistic input-output relationships by combining multiple expert models through an input-dependent gating mechanism. These models have become increasingly prominent in modern machine learning, yet their theoretical properties in the Bayesian framework remain largely unexplored. In this paper, we study Bayesian mixture-of-experts models, focusing on the ubiquitous softmaxbased gating mechanism. Specifically, we investigate the asymptotic behavior of the posterior distribution for three fundamental statistical tasks: density estimation, parameter estimation, and model selection. First, we establish posterior contraction rates for density estimation, both in the regimes with a fixed, known number of experts and with a random learnable number of experts. We then analyze parameter estimation and derive convergence guarantees based on tailored Voronoi-type losses, which account for the complex identifiability structure of mixture-of-experts models. Finally, we propose and analyze two complementary strategies for selecting the number of experts. Taken together, these results provide one of the first systematic theoretical analyses of Bayesian mixture-of-experts models with softmax gating, and yield several theory-grounded insights for practical model design.

Introduction

Mixture-of-experts (MoE) models, originally introduced by Jacobs et al. (1991), extend classical mixture models to supervised learning settings such as regression and classification. In these models, predictions are produced by combining several specialized submodels, called experts, whose contributions are weighted through a learnable gating mechanism that depends on the input variables. As a result, the model adaptively allocates different regions of the input space to different experts, allowing complex conditional relationships to be represented through a collection of simpler components. This architecture is particularly appealing when the conditional distribution of the response varies substantially across the input domain. Jacobs et al. (1991) introduced the original formulation and showed that such models can automatically decompose prediction tasks into simpler subtasks handled by specialized experts. The MoE framework has since been extended in several directions. Jordan and Jacobs (1994) introduced hierarchical mixtures-of-experts, in which experts are organized in a tree structure that enables further specialization across the input space. More recently, large-scale machine learning applications have renewed interest in these architectures. In particular, sparse routing mechanisms have been proposed to activate only a subset of experts for each input, enabling models with very large capacity while maintaining manageable computational costs (Shazeer et al., 2017). Such architectures have been employed in a wide range of modern deep learning applications, including language modeling (Jiang et al., 2024; DeepSeek-AI et al., 2024), computer vision (Riquelme et al., 1

Equal contribution.

2021), multi-task learning (Hazimeh et al., 2021), multi-modal learning (Han et al., 2024), domain adaptation (Li et al., 2023; Nguyen et al., 2025), continual learning (Le et al., 2024), and parameterefficient fine-tuning (Le et al., 2025; Diep et al., 2025). Despite their empirical success, these models are typically estimated using maximum likelihood methods, often based on the Expectation– Maximization algorithm, which may lead to overfitting and requires the model complexity to be fixed in advance (Masoudnia and Ebrahimpour, 2014). Uncertainty quantification also remains an important open challenge, particularly in scientific domains where reliable assessment of predictive and parameter uncertainty is essential. Addressing this issue within the maximum likelihood paradigm is often difficult, since likelihood-based procedures typically produce point estimates and do not naturally provide a coherent characterization of uncertainty in complex models such as MoE. Bayesian formulations of MoE models provide a natural framework for addressing these limitations by enabling uncertainty quantification and principled model selection. Early contributions include the fully Bayesian treatment of Peng et al. (1996), who used Markov chain Monte Carlo methods to infer both expert parameters and gating parameters. In parallel, Waterhouse et al. (1995) proposed a variational approach based on the minimization of variational free energy, providing a scalable approximation to the posterior distribution. These works highlighted the advantages of combining the modular structure of MoE models with Bayesian probabilistic inference, allowing uncertainty in both predictions and model structure to be quantified. Subsequent work has explored Bayesian approaches to model selection in MoE architectures. Ueda and Ghahramani (2002) proposed a variational Bayesian framework that allows the number of experts to be inferred from data through optimization of an evidence lower bound. Nonparametric Bayesian extensions have also been considered. For example, Rasmussen and Ghahramani (2002) proposed an infinite MoE model combining Dirichlet processes with Gaussian process experts, allowing the number of mixture components to grow with the size of the data. More broadly, several Bayesian nonparametric constructions have been developed to introduce covariate-dependent mixture weights, including dependent Dirichlet processes (MacEachern, 1999), hierarchical Dirichlet processes (Teh et al., 2006), nested Dirichlet processes (Rodriguez et al., 2008), and logistic stick-breaking processes (Ren et al., 2011), among many others; indeed, although many of these models were originally developed to capture dependence across a finite set of disjoint groups, they can also be interpreted as mechanisms for introducing input-dependent mixture weights, establishing a close conceptual connection with MoE models. Posterior inference in Bayesian MoE models remains challenging due to the complexity of the resulting posterior distributions, particularly in nonconjugate or hierarchical settings. Markov chain Monte Carlo methods provide asymptotically exact inference but may become computationally demanding in high-dimensional or large-scale problems. Variational inference methods provide a computationally attractive alternative by approximating the posterior distribution through optimization. For example, Bishop and Svensén (2003) developed a variational treatment of MoE models that yields efficient approximate inference while enabling model comparison through the evidence lower bound. In this paper, we focus on a widely used parametric formulation of MoE models, in which the gating mechanism is given by softmax weights and the experts are Gaussian distributions with input-dependent means. This formulation directly addresses regression tasks (see Section 7 for an extension to classification), which are among the most fundamental applications of MoE models and statistical learning in general. At the same time, concentrating on regression keeps the analysis tractable while still yielding insights that may extend to other settings, including unsupervised tasks 2

where MoE models are also commonly employed. Against this backdrop, we study the theoretical properties of Bayesian MoE models with softmax gating. In particular, we investigate the asymptotic behavior of the posterior distribution with respect to density estimation, parameter estimation, and model selection. In fact, although Bayesian MoE models have been widely used in practice, their theoretical properties remain comparatively less understood. Therefore, our goal is to provide one of the first systematic analyses of these models from a Bayesian perspective, establishing conditions under which posterior consistency and related guarantees can be obtained and, as a consequence, deriving practical insights for model design and specification. Plan of the aricle. The remainder of the article is organized as follows. In Section 2, we provide an overview of the literature related to our work. Section 3, we introduce the notation and mathematical framework used throughout the paper, together with some preliminary results on the identifiability of MoE models. Section 4 studies posterior convergence for density estimation with MoE models, while Section 5 focuses on parameter estimation. Section 6 presents theoretical results and empirical evidence on model selection criteria for MoE models. Although most of the paper concentrates on MoE models for regression tasks, Section 7 briefly discusses extensions to classification. Section 8 concludes the article. The proofs of all theoretical results and additional experimental details are collected in Appendix A and Appendix B, respectively.

Related literature

Beyond its connections to the literature on MoE models (see Gormley and Frühwirth-Schnatter, 2019, for a recent review in the statistical context), our work is also closely related to a well-established body of research in Bayesian asymptotics. In particular, our analysis builds on foundational results on posterior consistency (Schwartz, 1965; Barron et al., 1999; Ghosal et al., 1999; Walker, 2004) and posterior contraction rates (Ghosal et al., 2000; Shen and Wasserman, 2001; Walker et al., 2007). Following these contributions, a substantial literature has studied consistency and contraction properties for mixture models (Ghosal and van der Vaart, 2001, 2007b; Lijoi et al., 2005, among others). Our results parallel these analyses, although in a setting that is both more complex, due to the presence of covariates, and somewhat simpler, since we focus on finite mixtures rather than infinite-dimensional ones. Our analysis of parameter estimation is also related to recent developments on the estimation of mixing measures in mixture models, where Wasserstein distances have played a central role (Nguyen, 2013; Ho and Nguyen, 2016; Guha et al., 2021). In the presence of covariates, however, standard Wasserstein-based approaches become less effective. In particular, Ho et al. (2022) used a generalized Wasserstein distance to characterize the convergence behavior of parameter estimation in Gaussian MoE with covariate-free gating and the number of experts being over-specified. They showed that the rates for estimating exactly-specified parameters, which were fitted by one component, shared the same order as those for estimating over-specified parameters, which were fitted by at least two components. However, in practice, the former should be faster than the latter. This motivated the use of more refined tools based on Voronoi cells, following the recent work of Manole and Ho (2022), where they distinguished the estimation rates of exactly-specified parameters from those of over-specified parameters in the context of finite mixture models. Subsequently, Nguyen et al. (2023) generalized the use of Voronoi-based losses to a more challenging setting of Gaussian MoE with covariate-dependent gating function, showing that exactly-specified parameters enjoy significantly 3

faster estimation rates than over-specified parameters. More recently, Nguyen et al. (2026) examined the convergence of parameter estimation in a regression framework where the regression function took either the form of a softmax gating MoE or its hierarchical version. In addition to Wasserstein distance and Voronoi loss function, the Kullback-Leibler (KL) divergence was also utilized to study the maximum likelihood estimator for parameters of the MoE with each expert being a polynomial regression model by Mendes and Jiang (2012) studied, leading to useful insights on how many experts should be chosen. Finally, our study of model selection is closely connected to recent work on the behavior of posterior distributions when inferring the number of mixture components. This includes settings where consistency can be established (Miller, 2023; Ascolani et al., 2023), as well as well-known examples demonstrating that Dirichlet process mixtures can be inconsistent for this task (Miller and Harrison, 2014). Similarly, our findings on over-specified MoE models relate to the asymptotic behavior of overfitted mixture models, where the posterior tends to empty redundant components (Rousseau and Mengersen, 2011). As we will clarify, our preliminary empirical evidence supporting model selection via variational inference is also motivated by the recent theoretical advances of Zhang and Yang (2024); Wang and Yang (2024).

Formal setup and identifiability

To begin our formal analysis, assume that we observe a X × R-valued sample (Yi , Xi )ni=1 of inputresponse pairs, where X := [−1, 1]d and Y ∈ R. The class of models we consider is that of Softmax-gated Mixture of Gaussian Experts (SMoGE) densities, from which we assume the data are generated: iid (Yi , Xi ) ∼ gG⋆ , gG⋆ (Y, X) := fG⋆ (Y | X) × p(X), where p(·) is the uniform distribution2 on X and ⋆

⋆

G :=

K X

⋆ exp(α0j )δ(α⋆ ,β ⋆ ,σ⋆2 ) , 1j

j=1 ⋆

fG⋆ (Y | X) :=

K X

⋆ + X ⊤ α⋆ ) exp(α0j 1j N (Y | E(X, βj⋆ ), σj⋆2 ). PK ⋆ ⋆ ⋆ ) ⊤ exp(α + X α ℓ=1 0ℓ 1ℓ j=1

Here N (· | µ, σ 2 ) denotes the normal distribution with mean µ ∈ R and variance σ 2 > 0, E : Rd × Rp → R is the expert function describing the covariate-dependent mean of each component density, parametrized by β ∈ Rp , and G⋆ is a positive mixing measure that completely characterizes the distribution of the data. For any K ∈ N, denote by GK the set of positive measures with K support points ωj = (α1j , βj , σj2 ) ∈ Ω ⊂ Rd+p+1 and weights exp(α0j ), where α0j ∈ A ⊂ R. Let Θ := Ω × A. The specification of a Bayesian model in this setting starts with the choice of (i) a set of integers S K, and (ii) a prior Π supported on k∈K Gk . The posterior distribution is then obtained via Bayes’ rule n Y n Π(dG | (Yi , Xi )i=1 ) ∝ Π(dG) × fG (Yi | Xi ). (1) i=1 2

Since our focus is on the conditional distribution of Y given X, we fix the marginal distribution of X to a simple choice, although more general distributions could also be considered.

It is important to note S that, while more compact notationally, the convention of placing S a prior andk deriving a posterior on k∈K Gk is equivalent to placing a prior on the Euclidean space k∈K (Ω × A) of SMoGE parameters. At times we will switch to this alternative notation.

3.1

Losses for density estimation

Since we study the asymptotic behavior of the posterior distribution around the data-generating density, we first introduce several notions of divergence between SMoGE densities. In particular, for any two probability densities f and f ′ defined on a Borel set S ⊆ Rm , define the squared Hellinger distance Z p p 2 d2H (f, f ′ ) := f − f ′ dλ, S

where λ denotes a σ-finite measure on S. For any two SMoGE densities gG and gG′ with G ∈ GK , G′ ∈ GK ′ this definition specializes to Z 2 p p 2 dH (gG , gG′ ) = gG (Y, X) − gG′ (Y, X) d(Y, X) X×R ≡ EX∼p d2H (fG (· | X), fG′ (· | X)) . Also define the Kullback-Leibler (KL) divergence Z ′ K(f, f ) := log(f /f ′ )f dλ, S

which yields the following important notion of KL neighborhood (restricted to GK ⋆ ): B(G⋆ , ε) := G ∈ GK ⋆ : K(gG⋆ , gG ) ≤ ε2 for any ε > 0. Finally, denote by N (ε, M, d) the ε-covering number of a set M endowed with the metric d, while by ≳ and ≲ we will mean inequalities up to constants.

3.2

Losses for parameter estimation

In the previous subsection we considered losses that quantify discrepancies between SMoGE, or more general, density functions. In our theoretical analysis, however, we will also study the asymptotic behavior of the posterior distribution in terms of the SMoGE parameters themselves. This requires losses that are tailored to parameter estimation. SMoGE models present several challenges that require care in the construction of such losses. First, since K ⋆ is unknown and the statistical model may contain a different number of experts, we need a way to compare SMoGE parameters belonging to models with different numbers of experts. Second, even when the number of experts coincides, several identifiability issues arise, for instance permutation invariance. These observations, which will be discussed in detail later, prevent the use of simple Euclidean-type metrics between parameter vectors as meaningful measures of discrepancy. To solve this issue, we adopt a refined approach based on Voronoi cells (Manole and Ho, 2022; Nguyen et al., 2023). In particular, given any mixing measure G ∈ GK with K ≥ K ⋆ experts, we partition its components into the following Voronoi cells generated by the components of the ground-truth mixing measure G⋆ : Cj ≡ Cj (G) := {i ∈ {1, 2, . . . , K} : ∥ωi − ωj⋆ ∥ ≤ ∥ωi − ωℓ⋆ ∥, ∀ℓ ̸= j}, 5

(2)

⋆ , β ⋆ , σ ⋆2 ) for any j ∈ {1, 2, . . . , K ⋆ }. Based on this construcwhere ωi := (α1i , βi , σi2 ) and ωj⋆ := (α1j j j tion, we define the following Voronoi losses, which quantify the discrepancy between the parameters of G and G⋆ and will play a central role in our posterior asymptotic analysis: ⋆

⋆

L1 (G, G ) :=

K X X

⋆ exp(α0i ) − exp(α0j )

j=1 i∈Cj

K⋆ X X

⋆ exp(α0i ) ∥α1i − α1j ∥ + ∥βi − βj⋆ ∥ + |σi2 − σj⋆2 |

(3)

j=1 i∈Cj

and ⋆

⋆

L2 (G, G ) :=

K X X

⋆ exp(α0i ) − exp(α0j )

j=1 i∈Cj

⋆ exp(α0i )(∥α1i − α1j ∥ + ∥βi − βj⋆ ∥ + |σi2 − σj⋆2 |)

j∈[K ⋆ ]:|Cj |=1 i∈Cj

⋆ 2 exp(α0i )(∥α1i − α1j ∥ + ∥βi − βj⋆ ∥2 + |σi2 − σj⋆2 |2 ),

(4)

j∈[K ⋆ ]:|Cj |>1 i∈Cj

where [m] := {1, . . . , m} for any m ∈ N. Although the geometry induced by these losses will be discussed in detail in Section 5, it is clear that both L1 and L2 measure discrepancies between the component parameters of G and the closest component parameters of G⋆ . In the context of statistical estimation, this can be interpreted as comparing the parameters of the statistical model implied by G with the parameters of the component of the data-generating model implied by G⋆ that they estimate. As a preview of our analysis, L1 will be useful both for characterizing the support of the prior and for analyzing posterior convergence at the level of the SMoGE parameters when the model contains exactly the same number of components as G⋆ . In contrast, L2 will be crucial for characterizing parameter convergence when the statistical model contains potentially more experts than the ground-truth model.

3.3

Basic assumptions and model identifiability

As a final step before proceeding with the analysis, we state several assumptions that will be used throughout the paper and discuss some of their immediate consequences. Assumption 1. The parameter space Θ := Ω × A is compact, and the variance space is bounded away from 0. Assumption 2. There exists L > 0 such that, for all X ∈ X and β, β ′ ∈ Rp in the parameter space, |E(X, β) − E(X, β ′ )| ≤ L∥β − β ′ ∥. Assumption 3. The expert function E : [0, 1]d × Rp → R is identifiable. That is, for each K ∈ N the ′ ) such following holds. If there exist distinct parameter collections (β1 , β2 , . . . , βK ) and (β1′ , β2′ , . . . , βK that, for almost every x ∈ X , there exists a permutation function σx : {1, 2, . . . , K} → {1, 2, . . . , K} ′ }. satisfying E(X, βσx (i) ) = E(X, βi′ ) for all 1 ≤ i ≤ K, then {β1 , β2 , . . . , βK } ≡ {β1′ , β2′ , . . . , βK 6

⋆ − α⋆ , for 1 ≤ i < j ≤ K ⋆ , are Assumption 4. The discrepancies between gating parameters α1i 1j pairwise distinct. In addition, expert parameters (βj⋆ , σj⋆2 ), for 1 ≤ j ≤ K ⋆ , are also pairwise different. The same holds for the parameter configuration associated to Π-almost every mixing measure G.3

The first two assumptions are crucial to ensure expert parameter estimation consistency. As for the third assumption, it can be verified that, for example, if φ : R → R is any injective function, then an expert function of the form E(X, β) = φ(β ⊤ X), with β ∈ Rd , is identifiable. This assumption and the fourth one have important consequences for the identifiability of mixing measures, as shown in the following result. Proposition 1. Assume G is a mixing measure G=

K X

exp(α0j )δ(α1j ,βj ,σ2 ) j

j=1

such that gG (y, x) = gG⋆ (y, x) for almost every pair (y, x) ∈ R × X. Assume also that Assumption 3 holds, and Assumption 4 also holds for both G⋆ and G. Then K = K ⋆ and ⋆

G ≡ G⋆t0 ,t1 :=

K X

⋆ exp(α0j + t0 )δ(α⋆ +t1 ,β ⋆ ,σ⋆2 ) , 1j

j=1

for some t0 ∈ R and t1 ∈ Rd . Proposition 1, whose proof is in Appendix A.8, shows that SMoGE densities are identifiable up to translation of the gating parameters. This dependence arises from the softmax gating structure, ⋆ and α⋆ are translated to α⋆ + t as softmax weights remain unchanged if the gating parameters α0j 0 1j 0j ⋆ d and α1j + t1 , respectively, for some t0 ∈ R and t1 ∈ R . Formally, ⋆ + X ⊤ α⋆ ) exp(α0j 1j = , PK ⋆ PK ⋆ ⋆ ⋆ ⋆ ⋆ ⊤ ⊤ ℓ=1 exp(α0ℓ + t0 + X (α1ℓ + t1 )) ℓ=1 exp(α0ℓ + X α1ℓ ) ⋆ + t + X ⊤ (α⋆ + t )) exp(α0j 0 1 1j

(5)

for all j = 1, 2, . . . , K ⋆ . ⋆ ⋆ In light of this, if one fixes a pair of gating parameters to zero, for example α0K ⋆ = 0 and α1K ⋆ = 0, ⋆ ⊤ ⋆ then the numerator of the last softmax weight becomes equal to one, that is, exp(α0K ⋆ +X α1K ⋆ ) = 1. In this case, the softmax weights are no longer invariant to translations of the gating parameters and, as a consequence, the SMoGE density becomes identifiable, which is crucial for parameter estimation. In light of this, we formulate the following additional assumption. Assumption 5. Assume that the last pair of ground-truth gating parameters are zero, that is, ⋆ ⋆ α0K ⋆ = 0 and α1K ⋆ = 0.

Density estimation

We are now ready to tackle one of the core objectives of the paper, which is to study the asymptotic properties of the posterior distribution in the SMoGE model we have described. Recall that we 3

This holds, for instance, as long as the prior Π is of product form across experts and diffuse.

assume that the observed data are generated iid according to a fixed, yet unknown SMoGE model gG⋆ . A natural question is therefore whether, and at what rate, the posterior concentrates around G⋆ , and the losses introduced in Section 3 allow us to formalize the notion of being “close” to G⋆ . In this section, we begin by asking whether the posterior distribution in Equation (1) asymptotically concentrates within arbitrary Hellinger neighborhoods of the data-generating density gG⋆ . Two classical approaches address this question. The first studies posterior consistency in the Hellinger metric, that is, whether the posterior concentrates its mass in every neighborhood of the form {G : dH (gG , gG⋆ ) < ε} for all fixed ε > 0 (Barron et al., 1999; Ghosal et al., 1999; Walker and Hjort, 2001; Walker, 2004; Ghosal and Van der Vaart, 2017; Bariletto et al., 2025). The second studies contraction within shrinking neighborhoods of the form {G : dH (gG , gG⋆ ) < εn } for a decreasing contraction rate (εn )n∈N (Ghosal et al., 2000; Shen and Wasserman, 2001; Walker et al., 2007; Ghosal and van der Vaart, 2007a; Ghosal and Van der Vaart, 2017). Our aim is to investigate these properties under two specifications of interest for the integer set K that determines the support of the prior on the number of experts: the exactly-specified case where K = {K ⋆ }, and the over-specified case where K = {1, . . . , K} for some large enough K ≥ K ⋆ (featuring a random-a-priori number of experts).

4.1

Exactly-specified number of experts

We first consider the exactly-specified case, that is, when the set of SMoGE models supported by the prior comprises SMoGE densities with exactly K ⋆ components. While precise knowledge of K ⋆ is not necessarily a tenable assumption in practice, first analyzing this setting provides a useful benchmark to compare our subsequent results of more realistic scenarios, especially in terms of parameter estimation. To this end, the following theorem deals with posterior consistency in arbitrary Hellinger neighborhoods of the true density. Theorem 1. Assume that Π(GK ⋆ ) = 1. If Π ({G ∈ GK ⋆ : K(gG⋆ , gG ) < δ}) > 0 for all δ > 0, then for all ε > 0 lim Π ({G ∈ GK ⋆ : dH (gG , gG⋆ ) > ε} | (Xi , Yi )ni=1 ) = 0

n→∞

∞ a.s.-gG ⋆.

(6)

If in addition Assumptions 1 and 2 hold, then Equation (6) holds for all ε > 0 provided that ∀t > 0,

Π ({G ∈ GK ⋆ : L1 (G⋆ , G) < t}) > 0.

(7)

Theorem 1 shows that a minimal support condition on the prior, expressed either in terms of the KL divergence (Schwartz, 1965) or in terms of the Voronoi loss L1 , is sufficient to ensure that the posterior concentrates within arbitrary Hellinger neighborhoods of the true density. Inspection of the proof in Appendix A.1 reveals that the argument combines Schwartz’s KL support theorem with the observation that, despite their flexibility, SMoGE models satisfy the simple sequential identifiability condition introduced by Bariletto et al. (2025), which ensures Hellinger consistency. We now turn to contraction rates. Before stating the main result, we introduce an assumption describing a form of prior concentration around the true mixing measure. Assumption 6. There exists c > 0 such that the following holds for all sufficiently small t > 0: Π ({G ∈ GK ⋆ : L1 (G⋆ , G) ≤ ct}) ≥ t. 8

(8)

Assumption 6 can be interpreted as a quantitative version of the support condition stated in Equation (7). While that condition only requires every L1 neighborhood of G⋆ to be charged positive prior mass, Assumption 6 specifies how this mass must scale with the neighborhood radius. In particular, it requires that L1 (G⋆ , G) be approximately uniformly distributed near zero under the prior, ensuring that the latter does not place too little mass near G⋆ . With the aid of Assumption 6, we formulate the next theorem dealing with posterior contraction rates in Hellinger metric. Theorem 2. Consider a SMoGE model with Π(GK ⋆ ) = 1 and satisfying Assumptions 1, 2, and 6. Then for any sequence (Mn )n∈N such that limn→∞ Mn = ∞, o n p lim Π G ∈ GK ⋆ : dH (gG , gG⋆ ) ≥ Mn log(n)/n | (Xi , Yi )ni=1 = 0 n→∞

n -probability. in gG ⋆

Theorem 2 shows that, for density estimation, the posterior achieves a parametric contraction rate up to a logarithmic factor. The proof, reported in Appendix A.2, relies on the general theory of posterior contraction developed in Ghosal et al. (1999) (see also Ghosal and Van der Vaart, 2017). The argument proceeds by constructing an explicit cover of the model class and linking the KL neighborhood B(G⋆ , ε) with a corresponding neighborhood defined through the Voronoi loss L1 .

4.2

Over-specified random-a-priori number of experts

An appealing feature of Bayesian methods is that, in principle, they allow inference on any quantity for which a prior can be specified. In SMoGE models, a key quantity influencing both density and parameter estimation is the number of experts κ. In the previous subsection we assumed that κ = K ⋆ was fixed a priori. In practice, however, one rarely has exact knowledge of the data-generating number of experts, and may therefore wish to consider a range of possible values for κ ∈ K, so to achieve greater model flexibility or to hedge against model misspecification. The Bayesian framework is uniquely suited to accommodate this situation, as we show next. Let Πk denote the prior on the mixing measure G considered in Theorem 2, when the number of ⋆ experts κ equals k ∈ N (in Theorem SK 2, k = K ). To allow for multiple values of κ, we enlarge the support of the prior to OK := j=1 Gj and place a prior πκ directly on κ: G | κ = k ∼ Πk , κ ∼ πκ ,

(9)

with πκ ({1, . . . , K}) = 1 for some K ∈ N, chosen large enough so that K ≥ K ⋆ . That is, instead of requiring exact knowledge of the unknown parameter K ⋆ as in the exactly-specified case, we now only assume to be able to choose an large (but finite) upper-bound K. The resulting Parbitrarily K prior on the mixing measure is Π = k=1 πκ ({k})Πk , which induces the posterior Π(· | (Xi , Yi )ni=1 ). This formulation also allows inference on κ itself, although a thorough asymptotic analysis of this aspect is deferred to Section 6. We now deal with density estimation, and specifically posterior contraction rates, in this random-κ setting.4 4

Although we omit a formal statement, a result analogous to Theorem 1 can also be obtained in this setting, ensuring posterior consistency under a simple KL or L1 support condition for ΠK ⋆ .

Theorem 3. Consider a SMoGE model with a prior structure as in Equation (9), such that Assumptions 1, 2 and 6 hold. Then for any sequence (Mn )n∈N such that limn→∞ Mn = ∞, o n p lim Π G ∈ OK : dH (gG , gG⋆ ) ≥ Mn log(n)/n | (Xi , Yi )ni=1 = 0 n→∞

n -probability. in gG ⋆

According to Theorem 3, allowing the number of experts to be random a priori does not affect the Hellinger contraction rate of the posterior around the true density, provided that the prior on the number of experts has finite support, and a parametric rate (up to a logarithmic factor) is recovered in this setting as well.

Parameter and expert estimation

Having tackled density estimation, this section is devoted to the analysis of posterior convergence in terms of parameter estimation. To further motivate the need for such a refined analysis, we recall that, during the training of MoE architectures in deep learning, individual experts are meant to specialize in different subproblems or latent regions of the input distribution. For instance, in large language models built upon the MoE framework (Jiang et al., 2024; Liu et al., 2024; Google Gemini Team, 2025), some experts may learn syntactic patterns, others may focus on semantic representations, while the remaining ones capture broader contextual information. In our Bayesian setting, expert specialization can be formalized as posterior concentration on models whose distinct experts are close to unknown heterogeneous regimes in the underlying data-generating process. From a statistical perspective, this makes the rate of specialization a critical subject of study (Dai et al., 2024; Oldfield et al., 2024; Nguyen et al., 2024c). To operationalize this problem, an effective strategy is to focus on parameter estimation convergence, which in turn indicates how fast the parametrized experts specialize. In terms of sample complexity, faster convergence rates indicate that experts need less data to achieve a given degree of estimation accuracy or, in MoE jargon, specialization. Consequently, analyzing the asymptotics of parameter estimation not only provides theoretical insights into the dynamics of expert learning, but also informs strategies to improve the design of sample-efficient MoE models. There are two notable existing works on parameter estimation in Gaussian MoE models. First, Ho et al. (2022) study the effects of mean and variance experts on parameter estimation rates in Gaussian MoE with covariate-free gating. In particular, they show that if mean and variance experts meet some specific algebraic independence conditions, then the ground-truth parameters can be estimated at rates of order (log(n)/n)1/4 . Otherwise, the rates turn out to be substantially slower due to the interaction between mean and variance parameters expressed in terms of some underlying partial differential equations (PDEs). Second, Nguyen et al. (2023) consider a more challenging SMoGE model (in a maximum likelihood estimation framework) but limited to linear experts. They capture a PDE-type interaction between gating and expert parameters and demonstrate that parameter convergence rates depend on the solvability of a complex system of polynomial equations derived from that interaction. In this work, we further investigate the SMoGE model by characterizing expert structures that lead to improved parameter estimation rates compared to those in Nguyen et al. (2023), also connecting them to our overarching Bayesian approach. Before delving into the main findings of this section, it is worth noting that the posterior convergence on the density space established in Theorems 2 and 3 is a key component in characterizing 10

the convergence behavior of parameter estimation. More specifically, given the results of Theorems 2 and 3, the convergence analysis for parameter estimation boils down to constructing a discrepancy measure between the model parameters and their ground-truth counterparts, say L(G, G⋆ ), which is bounded above by the Hellinger distance between model and ground-ground truth densities. That is, more concretely: dH (gG , gG⋆ ) ≳ L(G, G⋆ ).

(10)

A bound of the above type, combined with Theorems 2 and 3, immediately translates into convergence rates for parameter estimation. In the sequel, we will conduct our convergence analysis of parameter estimation in two different settings, namely the exactly-specified and over-specified settings considered in the context of density estimation in Section 4.

5.1

Exactly-specified setting

To begin with, we investigate the exactly-specified setting in which the number of ground-truth experts K ⋆ is known. In doing so, we face some technical challenges, which we briefly discuss. In our proofs, a key step in characterizing parameter estimation rates is to decompose the density discrepancy gG − gG⋆ into a combination of linearly independent terms. This can be achieved by applying a Taylor expansion to the function F (Y |X, α, β, σ 2 ) := exp(X ⊤ α)N (Y |E(X, β), σ 2 ). In the process, it is important to ensure that the function F and its partial derivatives with respect to its parameters are linearly independent. If that is the case, when the Hellinger distance goes to zero, the coefficients of these terms in the decomposition, which are represented as parameter discrepancies, will also converge to zero. For that purpose, we need to introduce a non-trivial algebraic independence condition on the expert function E(X, β), as in Definition 1 below, which we refer to as first-order strong identifiability under the exactly-specified setting. Definition 1. An expert function X 7→ E(X, n β) is said to be o first-order strongly identifiable if it is ∂E(X,β) differentiable with respect to β, and the set : u ∈ [d] is linearly independent with respect to ∂β (u) X, where we denote β = (β (u) )du=1 . That is, for any β, if there exist tu ∈ R, for 1 ≤ u ≤ d, such that d X u=1

t(u)

∂E(X, β) = 0, ∂β (u)

for almost every X, then we have t(u) = 0 for all 1 ≤ u ≤ d. Example 1. It can be verified that a linear expert E(X, β1 , β0 ) = X ⊤ β1 + β0 satisfies the first-order strong identifiability condition. On the other hand, a constant expert E(X, β) = c is not first-order strongly identifiable. PK ⋆ ⋆ Next, it should be noted that the mixing measure G⋆ = j=1 exp(α0j )δ(α⋆1j ,βj⋆ ,σj⋆2 ) is only identifiable up to a permutation of parameters. In particular, if two mixing measures are equivalent, ⋆ ⋆ ⋆ , σ ⋆2 ) for some permutation G = G′ , we can only deduce that (α0j , α1j , βj , σj2 ) = (α0ϕ(j) , α1ϕ(j) , βϕ(j) ϕ(j) ϕ of the set {1, 2, . . . , K ⋆ }. To resolve this permutation invariance, we employ a Voronoi cell-based approach, which was briefly introduced in Section 3 and which we recall and expand on here. Given 11

an arbitrary mixing measure G with K ⋆ components, we partition its components into the following Voronoi cells, which are generated by the components of ground-truth mixing measure G⋆ : Cj ≡ Cj (G) := {i ∈ {1, 2, . . . , K ⋆ } : ∥ωi − ωj⋆ ∥ ≤ ∥ωi − ωℓ⋆ ∥, ∀ℓ ̸= j},

(11)

⋆ , β ⋆ , σ ⋆2 ) for any j ∈ {1, 2, . . . , K ⋆ }. Notably, when the where ωi := (α1i , βi , σi2 ) and ωj⋆ := (α1j j j components of G are sufficiently close to those of G⋆ , then each of these Voronoi cells will contain one element. As a result, we can deal with the aforementioned permutation-invariance of mixing measures with respect to parameters. The corresponding Voronoi-based loss capturing parameter discrepancies is given by ⋆

⋆

L1 (G, G ) :=

K X X

⋆ exp(α0i ) − exp(α0j )

j=1 i∈Cj ⋆

K X X

⋆ exp(α0i ) ∥α1i − α1j ∥ + ∥βi − βj⋆ ∥ + |σi2 − σj⋆2 | .

j=1 i∈Cj

Above, if a Voronoi cell Cj is empty, then we let the corresponding summation term be zero. Additionally, it can be checked that L1 (G, G⋆ ) = 0 if and only if G ≡ G⋆ . Thus, when L1 (G, G⋆ ) is ⋆ , β −β ⋆ , and σ 2 −σ ⋆2 are also small. This property indicates sufficiently small, the differences α1i −α1j i j i j ⋆ that L1 (G, G ) is an appropriate loss function for measuring parameter discrepancies. However, since the loss L1 (G, G⋆ ) is not symmetric, it is not a proper metric. The Voronoi loss function L1 (G, G⋆ ) also enjoys computational benefits, as it features a low computational complexity, namely O(K ⋆2 ). We finally note that the key benefits of Voronoi-based loss functions are more evident in the overspecified setting where the rates for estimating over-specified parameters (parameters of G⋆ fitted by at least two parameters) are slower than those for estimating exactly-specified parameters (parameters of G⋆ fitted by one parameter). In particular, Nguyen et al. (2023) have shown that Voronoi-based loss functions are capable of distinguishing these rates, while the generalized Wasserstein loss used in (Ho et al., 2022) cannot capture the rate difference. Given the Voronoi-based loss in Equation (3), we are ready to present our results on the posterior convergence rates for parameter estimation in exactly-fitted SmoGE models. Theorem 4. Suppose that the expert function X 7→ E(X, β) is first-order strongly identifiable. Under Assumptions 1, 2, 3, 4, 5, and 6, the lower bound dH (gG , gG⋆ ) ≳ L1 (G, G⋆ ) holds for any mixing measure G ∈ GK ⋆ . Then, the following holds for any sequence (Mn )n∈N such that limn→∞ Mn = ∞: n o p lim Π G ∈ GK ⋆ : L1 (G, G⋆ ) ≥ Mn log n/n | (Xi , Yi )ni=1 = 0. (12) n→∞

n -probability. in gG ⋆

The proof of Theorem 4 is in Appendix A.4. A few remarks regarding the results of Theorem 4 are in order. (i) Parameter estimation rates: The result in Equation (12), together with the structure of the Voronoi-based loss L1 (G, G⋆ ) in Equation (3), indicates p that the posterior rates for estimating ⋆ ⋆ ⋆ ⋆2 parameters exp(α0j ), α1j , βj , and σj are of the same log(n)/n order, which is parametric on the sample size n. 12

(ii) Expert estimation rates: Recall that the expert function E(X, β) is Lipschitz continuous with respect to its parameter β. Thus, we have |E(X, βi ) − E(X, βj⋆ )| ≤ L∥βi − βj⋆ ∥,

(13)

for all i ∈ Cj and j ∈ [K ⋆ ]. Then, the estimation rates for parameters βj⋆ in part (i) indicate that the estimators of experts E(X, βj∗ ) also enjoy the posterior contraction rates of the p parametric order log(n)/n, for j ∈ {1, 2, . . . , K ⋆ }. Therefore, our results reveal that one needs a polynomial number O(ϵ−2 ) of data points to estimate an expert E(X, βj⋆ ) with a given error tolerance ϵ > 0.

5.2

Over-specified setting

We now turn to the over-specified setting, in which the number of ground-truth experts is unknown and allowed to range up to an over-specified number K > K ⋆ . Similar to the exactly-specified setting, we also need to introduce a strong identifiability condition for the expert function in this setting to guarantee linear independence among terms in the Taylor expansion of F (Y |X, α, β, σ 2 ) = exp(X ⊤ α)N (Y |E(X, β), σ 2 ) when decomposing the density discrepancy gG − gG⋆ . However, since the number of ground-truth experts is unknown, the first-order strong identifiability condition in Definition 1 is not sufficient. Instead, it is necessary to control all derivatives of the expert function up to the second order. For this purpose, let us formally present the following second-order strong identifiability condition. Definition 2. An expert function X 7→ E(X, β) is said to be second-order strongly identifiable if it is twice differentiable with respect to β, and the following set is linearly independent with respect to X: ( ) ∂E(X, β) ∂ 2 E(X, β) (u) ∂E(X, β) , , X : u, v ∈ [d] . ∂β (u) ∂β (u) ∂β (v) ∂β (v) (u)

(uv)

That is, for any β, if there exist t1 , t2

(uv)

, t3

∈ R such that

2 ∂E(X, β) (uv) ∂ E(X, β) (uv) (u) ∂E(X, β) + t + t3 X (u) = 0, 2 (u) (u) (v) ∂β ∂β ∂β ∂β (v)

(u)

(uv)

for almost every X, then we have t1 = t2

(uv)

= t3

= 0 for all 1 ≤ u, v ≤ d.

Example 2. We can validate that the function E(X, β) = exp(X ⊤ β)/(1+exp(X ⊤ β)) is second-order strongly identifiable. By contrast, linear experts E(X, β1 , β0 ) = X ⊤ β1 + β0 are not second-order 1 ,β0 ) 1 ,β0 ) strongly identifiable because of the PDE relation ∂E(X;β = X ⊤ ∂E(X;β , which leads to the ∂β1 ∂β0 undesired linear dependence ∂2F ∂F = . ∂β1 ∂α∂β0

(14)

It can be seen that second-order strong identifiability implies the first-order version in Definition 1. In addition, this condition has substantive technical and intuitive interpretations. Technically, it helps eliminate any linear dependence among terms in the aforementioned Taylor expansion of the 13

density discrepancy, which is crucial for our proofs. Intuitively, second-order strong identifiability helps prevent potential adverse interactions among parameters, namely, where gating parameters α interact with expert parameters β1 and β0 as in Equation (14), thereby negatively parameter estimation efficiency. In fact, when the expert function is of linear form, which violates the first-order strong identifiability condition, Nguyen et al. (2023) showed that parameter estimation rates became significantly slow since, they hinge on the solvability of a system of polynomial equations originating precisely in such linear dependence interactions. Next, by relying on the following Voronoi loss (recall it from Section 3) ⋆

⋆

L2 (G, G ) :=

K X X

⋆ exp(α0i ) − exp(α0j )

j=1 i∈Cj

⋆ exp(α0i )(∥α1i − α1j ∥ + ∥βi − βj⋆ ∥ + |σi2 − σj⋆2 |)

j∈[K ⋆ ]:|Cj |=1 i∈Cj

⋆ 2 exp(α0i )(∥α1i − α1j ∥ + ∥βi − βj⋆ ∥2 + |σi2 − σj⋆2 |2 ),

(15)

j∈[K ⋆ ]:|Cj |>1 i∈Cj

we are ready to capture the convergence behavior of parameter estimation in Theorem 5. Theorem 5. Suppose that the expert function X 7→ E(X, β) is second-order strongly identifiable. Under Assumptions 1, 2, 3, 4, 5, and 6, the lower bound dH (gG , gG⋆ ) ≳ L2 (G, G⋆ ) holds for any mixing measure G ∈ OK . As a result, we obtain n o p lim Π G ∈ OK : L2 (G, G⋆ ) ≥ Mn log n/n | (Xi , Yi )ni=1 = 0. n→∞

n -probability. in gG ⋆

The proof of Theorem 5 is in Appendix A.5. There are three main implications of the above result which are worth highlighting. Before going into their details, let us introduce some preliminary ⋆ ), α⋆ , β ⋆ , σ ⋆2 and experts E(X, β ⋆ ) exactlynotions. We call the ground-truth parameters exp(α0j 1j j j j specified if they are fitted by one component, that is, |Cj | = 1, and over-specified if they are fitted by more than one component, that is, |Cj | > 1. (i) Estimation rates for exactly-specified parameters and experts: Combining the result of Theorem 5 and the formulation of the Voronoi loss L2 (G, G⋆ ), we deduce p that the rates for estimating ⋆ ), α⋆ , β ⋆ , σ ⋆2 are of order exactly-specified parameters exp(α0j log(n)/n, which is identical 1j j j to the rate uncovered under the exactly-specified setting. As a result, from the inequality in Equation (13), it follows that exactly-specified experts E(X, βj⋆ ) enjoy posterior estimation rates of the same order. (ii) Estimation rates for over-specified parameters and experts: On the other hand, when fitted ⋆ ), by more than one component, the rates for estimating over-specified parameters exp(α0j ⋆ , β ⋆ , σ ⋆2 become slower, in particular of order (log(n)/n)−1/4 , and the same rate holds α1j j j for estimating the over-specified experts E(X, βj⋆ ). Consequently, our results suggest that one needs a faster-growing number of data points to estimate over-specified experts than exactly-specified experts with the same error ϵ > 0: O(ϵ−4 ) compared to O(ϵ−2 ). 14

(iii) Sample efficiency of second-order strongly identifiable experts vs. linear experts: From the above remarks, we p see that the estimation rates for second-order strongly identifiable experts are either of order log(n)/n or (log(n)/n)1/4 . These rates are significantly faster than those for linear experts studied by Nguyen et al. (2023), which violate the second-order strong identifiability condition. In particular, Nguyen et al. (2023) showed that when fitted by two and three components, the rates for estimating linear experts are substantially slower, of orders (log(n)/n)1/8 and (log(n)/n)1/12 , respectively. Therefore, our analysis uncovers that using second-order strongly identifiable experts is more sample-efficient than linear experts. We conclude this section by remarking that, while expressed in the framework of Bayesian posterior convergence, our parameter convergence results, because they arise from a fundamental inequality between Hellinger and Voronoi losses, may be applied to any other estimation procedure (e.g., maximum likelihood) that may be analyzed through the lenses of density convergence.

Model selection

So far we have focused on density and parameter estimation. Another fundamental statistical task is selecting the model class that best describes the data-generating process. In the SMoGE setting, perhaps the most important quantity to determine is the number of experts κ. Indeed, if κ is too small, the model is misspecified, while on the other hand, as shown in Section 5, an over-specified number of experts negatively affects posterior convergence for parameter estimation, as reflected by the need to rely on the loss L2 rather than the more favorable L1 . Aside from theoretical considerations, model selection may also be relevant in applications where interpretability is important. For example, the number of experts may correspond to the number of latent populations across which input-response relationships differ, and identifying this number can therefore be of direct scientific interest. It should be noted from the outset that, in our context, the model selection problem may be illdefined or meaningful only within specific applications. Unless one assumes that the data-generating process itself belongs to a SMoGE family with K ⋆ components, a “true” number of experts need not exist. Nevertheless, one may still wish to select the number of experts that yields the best predictive performance, or particular applications may motivate ad hoc but practically effective criteria. For instance, Ludziejewski et al. (2024) study how the configuration of MoE models should scale with model size and training data in large language models. Rather than directly selecting the number of experts as a fixed architectural hyperparameter, the paper introduces a granularity parameter that controls the size of each expert relative to the feed-forward layer in the Transformer architecture (Vaswani et al., 2017). Increasing granularity effectively splits a standard expert into multiple smaller experts, thereby increasing the total number of experts while keeping the overall parameter budget comparable. Using extensive experiments, the authors derive scaling laws that relate model loss to the number of training tokens, total model parameters, and granularity, enabling the computation of compute-optimal configurations for a given training budget. Their results show that the commonly used design—setting expert size equal to the feed-forward layer width—is generally suboptimal; instead, the optimal number of experts should be determined jointly with model size and data scale through these scaling laws, balancing the benefits of finer expert specialization against the additional routing overhead. In what follows, instead, we adopt the same theoretical perspective used throughout the paper and define model selection relative to the data-generating process, investigating conditions under 15

which model selection consistency can be achieved. As in the rest of our analysis, this requires assuming that a true value K ⋆ exists, and therefore that the data-generating distribution is of SMoGE type. We view this as a natural first step toward understanding model selection in this flexible class of models. Before presenting the two strategies we consider, we note that both approaches are inherently Bayesian. In particular, both (i) learning κ directly from the data and (ii) using the evidence lower bound obtained from variational inference as a measure of model fit do not readily translate to the frequentist procedures commonly used to estimate mixture-of-expert models. This highlights a distinctive advantage of Bayesian methods, where coherent inference can be carried out for any parameter of interest, including parameters that determine the model size, and can therefore be naturally adapted to the task of model selection.

6.1

Model selection via a random-a-priori number of experts

Our first approach, similarly to Subsection 4.2, treats the number of experts κ as an unknown parameter and places a prior πκ on it. In contrast to the earlier setting, we now allow πκ to have potentially full support on N, so to enable model selection over an unbounded number of components. The parameter space5 on which we place the prior Π is therefore [ Θ∞ := Θk . k∈N

An application of Doob’s posterior consistency theorem (Doob, 1949; Miller, 2018), together with arguments inspired by the analysis of traditional mixtures in Nobile (1994); Miller (2023), yields the following result on posterior consistency for the number of experts. Theorem 6. Let Assumptions 1 and 2 hold. Moreover, assume the following: 1. The model is identifiable, meaning that gG = gG′ 6 implies G = G′ for all G, G′ ∈

k∈N Gk ;

2. All priors are continuous and factorize across experts. Then there exists a set Θ⋆ ⊆ Θ∞ such that Π(Θ⋆ ) = 1 and, for all θ⋆ ∈ Θ⋆ , lim Π(K = K(θ⋆ ) | (Xi , Yi )ni=1 ) = 1

n→∞

∞ a.s.-gG(θ , where, for any θ ∈ Θ∞ , K(θ) and G(θ) denote respectively the number of experts and the ⋆) mixing measure associated with θ.

Theorem 6 implies that, whenever the data-generating parameter belongs to a set Θ⋆ of full prior measure, the posterior asymptotically recovers the true number of experts. In other words, the Bayesian procedure is consistent for model selection for such data-generating processes. The caveat, which is typical of results based on Doob’s consistency theorem (Doob, 1949), is that the set Θ⋆ is 5

Without loss of generality, and as anticipated in Section 3, we now adopt notation in which priors and posteriors are distributions on Θ∞ rather than on the space of mixing measures. S 6 With a slight abuse of notation, for any mixing measure G ∈ j∈N Gj , we denote by gG both the resulting SMoGE joint density and the corresponding probability distribution. In the statement of the theorem we mean equality of the two probability measures, which implies equality of the densities almost everywhere. See Proposition 1 above for a characterization of the scope of this assumption.

not characterized explicitly. Rather, it arises from the proof, which relies on martingale almost sure convergence arguments, and as a consequence, verifying whether a specific parameter value belongs to this set is generally impossible. This limitation is the price paid for obtaining a result of such generality and strength in terms of model selection consistency (Nobile, 1994; Miller, 2023). We conclude with a natural corollary of Theorem 6, which shows that, under the same setting with a prior on κ having potentially full support on N, posterior Hellinger contraction also holds for those parameter values in Θ⋆ . Corollary 1. Let θ⋆ ∈ Θ⋆ , G⋆ = G(θ⋆ ), and K ⋆ = K(θ⋆ ), with Θ⋆ as defined in Theorem 6. Consider the SMoGE model with a prior on the number of experts having full support on N. Then, under the assumptions of Theorem 2 and Theorem 6, the following holds for any sequence Mn → ∞:      [ p lim Π  G ∈ Gj : dH (gG , gG⋆ ) ≥ Mn log n/n | (Xi , Yi )ni=1  = 0. n→∞   j∈N

n -probability. in gG ⋆

The parameter estimation results of Section 5 may also be extended to this setting, but we do not do so explicitly here as such extensions are minimal in light of Corollary 1.

6.2

Model selection by variational inference

While computation is not the primary focus of this paper, fitting Bayesian SMoGE models would typically rely on two main approaches: Markov chain Monte Carlo (MCMC), which constructs a Markov chain whose stationary measure is the target posterior (Robert and Casella, 2004; Gelman et al., 2013), and Variational Inference (VI), which approximates the posterior by finding the closest distribution within a tractable family (Jordan et al., 1999; Blei et al., 2017). As we detail below, VI naturally provides a mechanism for performing model selection as a byproduct of the fitting process. Specifically, VI proceeds by selecting a variational family Q = {qµ : µ ∈ M} of distributions over the parameter space, and optimizing the evidence lower bound (ELBO), E(µ) := Eθ∼qµ [log p ((Xi , Yi )ni=1 , θ) − log qµ (θ)] , with respect to the variational parameters µ ∈ M ⊆ Rdµ ; here, p ((Xi , Yi )ni=1 , θ) denotes the joint distribution of the data and parameters under the Bayesian prior-likelihood model. For simpler models and specific choices of Q, E(µ) can often be maximized via closed-form coordinate ascent updates (Bishop, 2006; Wainwright and Jordan, 2008). However, for more complex likelihoods like the SMoGE models considered here, optimization typically requires stochastic gradient-based methods. In this setting, the gradient of the ELBO is approximated using Monte Carlo samples from qµ , often in conjunction with variance reduction techniques (Ranganath et al., 2014; Kucukelbir et al., 2017). Returning to the problem of model selection, the ELBO provides a highly practical criterion. Because the ELBO bounds the log marginal likelihood, it serves as a computationally tractable proxy for model evidence and can guide the selection of model complexity, such as the number of experts in a SMoGE architecture. In our framework, we advocate complementing the prior over the number of experts with an ELBO-based heuristic: by fitting several SMoGE models across a 17

1.0

0.8

Proportion of simulations with highest ELBO

1.0

0.6 0.4 0.2 0.0

50 Sample size K=1

Model size K=2 K=3

0.6 0.4 0.2 0.0 100

100

K=4

500

K=1

(a) Data-generating process with K ⋆ = 2, d = 2

1000 Sample size K=2

Model size K=3 K=4

2000

K=5

K=6

(b) Data-generating process with K ⋆ = 4, d = 6

Figure 1: Model selection by VI-derived ELBO maximization. As the sample size increases, the procedure selects the correct number of experts with high probability. Win proportions are approximated over 50 independent simulations per sample size. plausible range of expert counts, κ = K, we select the model that maximizes the ELBO as the one exhibiting the highest approximate evidence. While employing the ELBO for model selection in mixtures of experts is not a new strategy (Waterhouse et al., 1995; Ueda and Ghahramani, 2002), rigorous analyses establishing the frequentist consistency of this procedure have only recently emerged for non-singular models (Zhang and Yang, 2024) and specific mixture specifications (Wang and Yang, 2024). Extending these highly technical consistency results from simple mixtures to our more complex SMoGE formulation is beyond the scope of this paper. Nevertheless, we view this as a highly promising strategy for future theoretical validation, primarily because VI scales significantly better than the transdimensional MCMC algorithms (Green, 1995) that would otherwise be required to fit SMoGE models with a random number of experts. To support this direction, we conclude the section by reporting preliminary empirical evidence validating the ELBO-based model selection procedure. We conducted a series of stylized synthetic experiments, whose results we describe next. Figure 1 illustrates the procedure’s performance across two scenarios: a simpler setup with K ⋆ = 2 true linear experts in d = 2 dimensions (Figure 1(a)), and a more complex K ⋆ = 4 linear expert architecture in d = 6 dimensions (Figure 1(b)). For both setups, we evaluate a range of candidate model sizes by maximizing the ELBO using black-box variational inference (Ranganath et al., 2014; Kucukelbir et al., 2017) over 50 independent simulations. For each sample size, we record the proportion of trials in which each candidate K achieved the highest final evidence lower bound. The results provide empirical evidence for consistency: at small sample sizes, the ELBO’s inherent complexity penalty conservatively favors fewer experts, but as the sample size increases, the selection procedure overwhelmingly identifies the true number of experts K ⋆ . It is worth noting that the true datagenerating mechanisms in these simulations, as described in Appendix B, are not strictly SMoGE models, as they employ deterministic expert assignments that introduce sharp discontinuities in the conditional mean structure. Nonetheless, model selection consistency, in this looser sense of identifying the correct number of underlying linear regimes, appears to hold robustly. To further validate the robustness of the ELBO-based selection, we also evaluate its performance on datasets generated with a logit-based gating mechanism that closely resembles the softmax structure of true SMoGE models. In this setup, we fix the sample size to n = 500 and vary the true number of regimes (K ⋆ ∈ {1, 2, 3}), the dimensionality (d ∈ {2, 4}), and the “sharpness” of 18

Table 1: Proportion of simulations in which each candidate K attained the highest ELBO over 100 independent trials. The true data-generating process relies on a max-logit gating structure evaluated at two different separation scales. Bold values indicate the most frequently selected model size. (a) Low gating sharpness (separation = 5.0)

K⋆

K=1

K=2

K=3

K=4

K=5

K=6

K=7

2 2 4

1 2 3

1.00 0.00 0.00

0.00 0.76 0.00

0.00 0.20 0.73

0.00 0.04 0.20

0.00 0.00 0.06

0.00 0.00 0.01

0.00 0.00 0.00

(b) High gating sharpness (separation = 10.0)

K⋆

K=1

K=2

K=3

K=4

K=5

K=6

K=7

2 2 4

1 2 3

1.00 0.00 0.00

0.00 0.84 0.00

0.00 0.14 0.78

0.00 0.02 0.21

0.00 0.00 0.01

0.00 0.00 0.00

the boundary separating the expert regions. Table 1 reports the proportion of 100 independent simulations in which each candidate K ∈ {1, . . . , 7} attained the highest ELBO. The procedure reliably identifies the correct K ⋆ as the most frequent winner across all configurations. Notably, selection accuracy improves when the gating sharpness is higher (Table 1b), as the clearer separation between expert regimes naturally reduces the posterior overlap and makes the true model complexity easier to distinguish. Furthermore, we observe that increasing the true number of regimes K ⋆ predictably slows down the convergence of the algorithm, as the variational optimization landscape becomes inherently more complex to navigate.

Extension to mixture-of-experts for classification

Throughout the article, we have focused on SMoGE models for regression tasks and have analyzed their theoretical properties in depth. Another fundamental supervised learning problem, for which Bayesian MoE models may be effectively employed to increase model flexibility, is classification. Although we have not treated this setting explicitly, it can be incorporated into our framework with only minor modifications, as follows. Suppose that the response variable takes values in a finite set of classes {1, . . . , S}. In this case, the data may be assumed to be drawn according to the following Softmax Mixture of Multinomial Experts (SMoME) model: ⋆

⋆

G :=

K X

⋆ exp(α0j )δ(α⋆1j ,βj⋆ ) ,

j=1 ⋆

fG⋆ (Y = s | X) :=

K X j=1

⋆ + X ⊤ α⋆ ) exp(α0j 1j

PK ⋆

⋆ ⊤ ⋆ ℓ=1 exp(α0ℓ + X α1ℓ )

exp Es (X, βj⋆ )

× PS

⋆ r=1 exp Er (X, βj )

for s = 1, . . . , S, where Es (X, β) denotes the expert function associated with class s, parametrized 19

by β. In this case, fG⋆ is a density with respect to the counting measure on {1, . . . , S}, and the usual Bayesian pipeline involving prior specification and posterior computation can be implemented without difficulty. While we omit the details for brevity, extending much of our theoretical analysis to this setting presents no particular difficulty. In particular, posterior contraction for density estimation may still be established using the general theory developed by Ghosal et al. (2000), while parameter estimation results based on Voronoi losses can be adapted from the recent contribution of Nguyen et al. (2024a).

Discussion

In this paper, we have presented a theoretical study of Bayesian MoE models for Gaussian regression with softmax gating. We focused on three key aspects: density estimation, parameter estimation, and model selection. For density estimation, we derived posterior contraction rates for both fixed and random numbers of experts. For parameter estimation, we established convergence guarantees using tailored Voronoi-type losses, which account for the complex identifiability structure of MoE models in both correctly specified and over-specified settings. Finally, we studied model selection using a random-a-priori number of experts and variational inference.

8.1

Implications of our work

Our results have several practical implications for working with MoE models. In particular, we have highlighted the role of Voronoi-type losses in characterizing prior support conditions that ensure Hellinger contraction and in guiding parameter convergence; this provides a useful tool to design priors with good asymptotic behavior and to assess convergence speed as a function of the sample size. Moreover, regarding model selection, our analysis shows that placing a prior with sufficiently large support on the number of experts can lead to model selection consistency, illustrating the robustness of the Bayesian approach for SMoGE models. At the same time, our parameter estimation results warrant some caution: overly diffuse priors on the number of experts can harm performance, as parameter contraction rates deteriorate when too many redundant experts are included. In light of these considerations, our preliminary empirical results indicate that using the ELBO from VI procedures may offer a practical alternative for model selection. This approach avoids the computational overhead of treating the number of experts as random, while remaining effective provided that one focuses on a reasonable range of candidate values when running VI.

8.2

Future directions

Several directions for future work remain. First, while we studied densely activated MoE models—where all experts are active for each data point—modern architectures often use sparsely activated experts (Shazeer et al., 2017; Nguyen et al., 2024b) to increase capacity while keeping computation manageable. Extending our theoretical results to sparse MoE models would help provide a solid understanding of these widely used, scalable architectures. Second, all of our results are based on an inherent well-specification assumption, which may be hard to justify in practice. Given the remarkable practical success of MoE models in capturing complex data-generating processes, it would be interesting to investigate their behavior under 20

misspecification (e.g., in an under-specified setting with fewer components than K ⋆ ), and to study whether this class of models offers benefits, in terms of improved approximation capacity, compared to more traditional models like mixtures with input-free gating. Third, the connections between MoE models and Bayesian nonparametric models with covariate dependence are still not fully explored. It would be interesting to see if the tools we developed here for finite mixtures can be applied to more complex, infinite-dimensional settings. Finally, while we focused on traditional Bayesian inference, recent approaches to uncertainty quantification—sometimes called “post-Bayesian” or “generalized Bayesian” methods (Bissiri et al., 2016; Knoblauch et al., 2022; Fong et al., 2023; Fortini and Petrone, 2024)—offer alternative ways to perform inference using loss functions and predictive rules rather than the standard prior-likelihood framework. These methods have been applied successfully to mixture models (Rodríguez et al., 2025) and deep learning pipelines (Lee et al., 2023; Wu and Williamson, 2024; Ng et al., 2025; Bariletto and Walker, 2026), which are closely connected to modern MoE architectures. Extending such approaches to MoE models could provide scalable and robust alternatives to conventional posterior inference.

Appendices A

Proofs

In this appendix, we report the proofs of the theoretical results presented in the main body of the article.

A.1

Proof of Theorem 1

As for the first statement, the assumption it requires is equivalent to the classical KL support condition by Schwartz (1965). Moreover, one easily checks that, if a sequence of SMoGE models gG1 , gG2 , . . . converges weakly to gG⋆ , it must converge in Hellinger distance as well. In the terminology of Bariletto et al. (2025), the model is sequentially identifiable and consistency immediately follows (see also Bariletto and Walker, 2025; Walker et al., 2005). As for the second statement, if Assumptions 1 and 2 also hold, then the upcoming proof of Theorem 3 shows that any set of the form {G ∈ GK ⋆ : K(gG⋆ , gG ) < δ} always includes a set of the form {G ∈ GK ⋆ : L1 (G⋆ , G) < t} for some t > 0, as long as δ > 0 is small enough. Therefore, by the assumption stated in the theorem, we have Π({G ∈ GK ⋆ : K(gG⋆ , gG ) < δ}) ≥ Π({G ∈ GK ⋆ : L1 (G⋆ , G) < t}) > 0. The arguments presented in the last paragraph conclude the proof of the second statement in the theorem.

A.2

Proof of Theorem 2

Clearly, the statement of Theorem 2 is a particular case of the statement of Theorem 3 (setting πκ ({K ⋆ }) = 1), so we refer the reader to the next proof for a more general treatment.

A.3

Proof of Theorem 3

We rely on the following result, which is an adaptation of Theorem 8.11 in Ghosal and Van der Vaart (2017). Theorem 7 (Ghosal and Van der Vaart (2017)). Let (εn )n∈N and (ε̄n )n∈N be such that εn , ε̄n ≥ n−1/2 for all large n ∈ N. Moreover assume that, for all large j ∈ N and some sequence of sets O n ⊆ OK , the following conditions are satisfied: 1. supε≥εn log N (ε/2, {G ∈ O n : dH (gG , gG⋆ ) ≤ 2ε}, dH ) ≤ nε2n ; Π(OK \O n ) −Dn nε̄2n for some sequence D → ∞; 2. Π(B(G ⋆ ,ε̄ )) = o e n n n

: jεn <dH (gG ,gG⋆ )≤2jεn }) 3. Π({G∈O ≤ enεn j /16 . Π({G∈OK : K(g ⋆ ,gG )≤ε2 }) 2 2

Then, for any positive sequence Mn → ∞, we have lim Π {G ∈ OK : dH (gG , gG⋆ ) ≥ Mn εn } | (Xi , Yi )ni=1 = 0.

n→∞ n -probability. in gG ⋆

Preliminarily, we note that the statement of Theorem 8.11 in Ghosal and Van der Vaart (2017) is based on distances more general than dH , as long as they satisfy certain testing conditions for some parameters ξ and K (in the notation of the book). Nevertheless, on page 197, the authors show that for the Hellinger metric these conditions are satisfied for ξ = 1/2 and K = 1/8, which we directly incorporate in the above version of the result. In what follows, we show that Conditions 1, 2, and 3 in Theorem 7 are satisfied by our SMoGE model under Assumptions 1, 2 and 6 and the following choices: p • εn = ε̄n = M log n/n for some large enough M > 0 • O n ≡ OK for all n ∈ N Condition 1. We begin by noting that sup log N (ε/2, {G ∈ O n : dH (gG , gG⋆ ) ≤ 2ε}, dH ) ≤ sup log N (ε/2, OK , dH )

ε≥εn

= log N (εn /2, OK , dH ), so it is enough to show

log N (εn /2, GK , dH ) ≤ nε2n .

Moreover, it is easily checked that, by choosing M large enough, the above condition is implied by log N (ε, OK , dH ) ≲ log(1/ε) for all ε > 0. Finally, a further simplification comes from noting that d2H (f, f ′ ) ≤ ∥f − f ′ ∥L 1 for any two density functions f, f ′ , where ∥ · ∥L 1 denotes the L 1 norm on the space of integrable functions. In particular, this implies that we are left to show log N (ε, OK , ∥ · ∥L 1 ) ≲ log(1/ε). We now proceed as follows. We will explicitly construct an ε-cover of GK in ∥ · ∥L 1 , whose cardinality NK we show to be O(ε−1 ). Because this O(ε−1 ) result does not depend on theSspecific choice of K ∈ N, it holds for all k ∈ {1, . . . , K}. Therefore, we will have shown that OK = K k=1 Gk PK −1 itself has an ε-cover in ∥ · ∥L 1 of cardinality k=1 Nk = O(ε ). Therefore, Condition 1 will be satisfied. So let ∆ := {(α0 , α1 ) ∈ Rd+1 : ∃(β, σ 2 ) ∈ Rp+1 s.t. (α0 , α1 , β, σ 2 ) ∈ Θ}, Γ := {(β, σ 2 ) ∈ Rp+1 : ∃(α0 , α1 ) ∈ Rd+1 s.t. (α0 , α1 , β, σ 2 ) ∈ Θ} be the projections of Θ onto the subspaces of gate parameters and expert parameters, respectively. Because Θ is compact in the Euclidean norm ∥ · ∥, so are ∆ and Ω. This implies that one can cover them with at most C/εd+1 and C ′ /εp+1 Euclidean balls of radius ε > 0, where C, C ′ > 0 are scaling constants (Wainwright, 2019). Denote by O∆ and OΓ the sets of centers of these ε-covers and, for each G ∈ GK such that G=

K X

exp(α0j )δ(α1j ,βj ,σ2 ) j

K for some (α0j , α1j , βj , σj2 )K j=1 ∈ Θ ,

j=1

define Ḡ ∈ GK as G1 :=

K X

exp(ᾱ0j )δ(ᾱ1j ,βj ,σ2 ) , j

j=1

where

(ᾱ0j , ᾱ1j ) ∈ arg

min (a0 ,a1 )∈O∆

∥(a1 , a2 ) − (α0j , α1j )∥.

for all j = 1, . . . , K. In turn, for the same G ∈ GK , define G2 :=

K X

exp(ᾱ0j )δ(ᾱ1j ,β̄j ,σ̄2 ) , j

j=1

where

(β̄j , σ̄j2 ) ∈ arg

∥(b, s2 ) − (βj , σj2 )∥.

min

(b,s2 )∈O

for all j = 1, . . . , K. Intuitively, G1 replaces the gate parameter vectors with the closest ones among the centers O∆ , and G2 further replaces the expert parameter vectors with the closest ones among the centers OΓ . Now denote SoftMax : RK → [0, 1]K , z 7→

exp(z1 ) PK

exp(zK )

ℓ=1 exp(zℓ )

, · · · , PK

ℓ=1 exp(zℓ )

which is a 1-Lipschitz continuous function in the Euclidean norm (Gao and Pavel, 2017), meaning that ∥SoftMax(z) − SoftMax(z ′ )∥ ≤ ∥z − z ′ ∥ for all z, z ′ ∈ RK . From this and the boundedness of X, we deduce Z Z ∥gG − gG1 ∥L 1 = |fG (Y | X) − fG1 (Y | X)|dY p(X)dX X

Z X K

exp(ᾱ0j + X ⊤ ᾱ1j ) exp(α0j + Xi⊤ α1j ) − PK p(X)dX PK ⊤ ⊤ X j=1 ℓ=1 exp(α0ℓ + X α1ℓ ) ℓ=1 exp(ᾱ0ℓ + X ᾱ1ℓ ) Z ⊤ ⊤ K = SoftMax (α0j + X ⊤ α1j )K − SoftMax (ᾱ + X ᾱ ) 1K p(X)dX 0j 1j j=1 j=1 X √ Z ⊤ K SoftMax (α0j + X ⊤ α1j )K p(X)dX ≤ K j=1 − SoftMax (ᾱ0j + X ᾱ1j )j=1 X Z √ ⊤ K ≤ K (α0j + X ⊤ α1j )K j=1 − (ᾱ0j + X ᾱ1j )j=1 p(X)dX X √ Z ⊤ K ≤ K (α0j − ᾱ0j )K j=1 + (X (α1j − ᾱ1j ))j=1 p(X)dX v X v  Z u Z u K K X uX u √ 2 t = K t |α0j − ᾱ0j |2 p(X)dX + |X ⊤ (α1j − ᾱ1j )| p(X)dX  ≤

j=1

≲ ε. 24

j=1

Moreover (16) ∥gG1 − gG2 ∥L 1 Z X K exp(ᾱ0j + X ⊤ ᾱ1j ) ≤ ∥N (· | E(X, βj ), σj2 ) − N (· | E(X, β̄j ), σ̄j2 )∥L 1 p(X)dX. PK ⊤ X j=1 ℓ=1 exp(ᾱ0ℓ + X ᾱ1ℓ ) For all j = 1, . . . , K and X ∈ X, ∥N (· | E(X, βj ), σj2 ) − N (· | E(X, β̄j ), σ̄j2 )∥L 1 ≤ ∥N (· | E(X, βj ), σj2 ) − N (· | E(X, β̄j ), σj2 )∥L 1 + ∥N (· | E(X, β̄j ), σj2 ) − N (· | E(X, β̄j ), σ̄j2 )∥L 1 . Using a first-order Taylor expansion on the mean of the normal density around E(X, β̄j ) and invoking the boundedness of X and the Lipschitz property of E(X, ·), we get ∥N (· |E(X, βj ), σj2 ) − N (· | E(X, β̄j ), σj2 )∥L 1 Z ∂ ≤ E(X, βj ) − E(X, β̄j ) N (· | µ, σj2 ) dY µ=µ⋆ R ∂µ ≲ ε. In a similar fashion, we obtain ∥N (· | E(X, β̄j ), σj2 ) − N (· | E(X, β̄j ), σ̄j2 )∥L 1 ≲ ε and therefore, by Equation (16), ∥gG1 − gG2 ∥L 1 ≲ ε. Thus, for all G, we have shown the existence of a measure G2 ∈ GK such that ∥gG − gG2 ∥L 1 ≤ ∥gG − gG1 ∥L 1 + ∥gG1 − gG2 ∥L 1 ≲ ε, with G2 belonging to the set

Ξ :=

 K X

exp(a0j )δ(a1j ,bj ,s2 ) ∈ GK : (a0j , a1j ) ∈ O∆ , (bj , s2j ) ∈ OΓ

 





j=1

Therefore, for some large enough C1 > 0, Ξ yields a C1 ε-cover of GK with respect to ∥·∥L 1 . Moreover, recalling that the cardinalities of O∆ and OΓ are upper-bounded respectively by C/εd+1 and C ′ /εp+1 , the cardinality of Ξ is at most C2 K/εd+p+2 for some constant C2 > 0. This yields the desired result log N (ε, GK , ∥ · ∥L 1 ) ≲ log(1/ε). Condition 2. Because OK = O n for all n ∈ N, Condition 2 is trivially satisfied. 25

Condition 3. We begin by noting that Condition 3 is implied by 2 Π( G ∈ OK : K(gG⋆ , gG ) ≤ ε2n ) ≥ e−nεn /16 , and moreover that Π( G ∈ OK : K(gG⋆ , gG ) ≤ ε2n ) ≥ Π( G ∈ GK ⋆ : K(gG⋆ , gG ) ≤ ε2n ). Therefore, we focus on proving that 2 Π( G ∈ GK ⋆ : K(gG⋆ , gG ) ≤ ε2n ) ≥ e−nεn /16 Notice that, for any G ∈ GK ⋆ , we have Z K(gG⋆ , gG ) ≡

K(fG⋆ (· | X), fG (· | X)) p(X)dX,

where, for each fixed X, fG⋆ (· | X) and fG (· | X) are Gaussian mixture densities with K ⋆ components. The means and variances lie in compact sets (and the variances are bounded away from 0, so say they belong to [v, V ] with 0 < v < V < ∞). We now aim to show that there exists δ ∈ (0, 1] (depending on the parameter space and G⋆ , but not on G nor X) such that Z +∞

fG⋆ (Y | X)

−∞

fG⋆ (Y | X) fG (Y | X)

δ dY < ∞

Because of the compactness of the space of means and the Gaussianity of each component, there Y2

exists Y1 > 0 such that, for all Y > Y1 , we have supX∈X fG⋆ (Y | X) ≲ e− 2V . Similarly, there exists Y2

for some constant C > 0, where the bound on the first addendum comes from compactness of the parameter space and truncation. The second addendum is then made finite by choosing δ ∈ (0, min{1, [V − v]−1 v}). Because these inequalities are uniform over values of X, we also obtain Z

Z +∞

gG⋆ (Y, X)

−∞

gG⋆ (Y, X) gG (Y, X)

δ dY dX ≤ Mδ < ∞

where Mδ is constant across G. The above result allows us to invoke Theorem 5 in Wong and Shen (1995), which yields that, for all G and ε such that d2H (gG , gG⋆ ) ≤ ε2 < 21 (1 − e−1 ), we have K(gG⋆ , gG ) ≤ Cε2 log(1/ε) for some 26

p C > 0 depending on the parameter space. Now define the function H(ε) := ε/ 2 log(1/ε), which is such that H 2 (ε) log(1/H(ε)) ≤ ε2 for all ε > 0 small enough. Therefore, for n large enough, B(G⋆ , εn ) = G ∈ GK ⋆ : K(gG⋆ , gG ) ≤ ε2n ⊇ G ∈ GK ⋆ : K(gG⋆ , gG ) ≤ H 2 (εn ) log(1/H(εn )) ⊇ G ∈ GK ⋆ : d2H (gG⋆ , gG ) ≲ H 2 (εn ) ⊇ G ∈ GK ⋆ : ∥gG⋆ − gG ∥L 1 ≲ H 2 (εn ) ⊇ G ∈ GK ⋆ : L(G⋆ , G) ≲ H 2 (εn ) , where the last inclusion holds for any loss function (17)

L(G⋆ , G) ≳ ∥gG⋆ − gG ∥L 1 . In that case, we get Π(B(G⋆ , εn )) ≥ Π

G ∈ GK ⋆ : L(G, G⋆ ) ≲ H 2 (εn )

So it suffices to require Π

G ∈ GK ⋆ : L(G, G⋆ ) ≲ H 2 (εn )

or equivalently, given our choice εn = M

≥ e−nεn /16 ,

p log(n)/n,

M 2 /16 log n 1 DM 2 ⋆ ≥ , Π G ∈ GK ⋆ : L(G, G ) ≤ n log n − log log n − 2 log M n where D is an explicit constant. Because the fraction is of constant order, we can simply update D and require M 2 /16 DM 2 1 ⋆ ⋆ Π G ∈ GK : L(G, G ) ≤ ≥ . n n Moreover, choosing M large enough and requiring n c o 1 ≥ , G ∈ GK ⋆ : L(G, G⋆ ) ≤ n n

Π for some c > 0, we would obtain Π

DM 2 G ∈ GK ⋆ : L(G, G ) ≤ n ⋆

M 2 /16 n c o 1 1 ⋆ ≥ Π G ∈ GK ⋆ : L(G, G ) ≤ ≥ ≥ . n n n

Hence, replacing t with 1/n, it is sufficient to require Π ({G ∈ GK ⋆ : L(G, G⋆ ) ≤ ct}) ≥ t for some c > 0 and all t > 0 small enough. 27

We are left to show that Equation (17) holds by taking L = L1 as defined in Equation (3). To that end, it suffices to show that the relation holds for ∥fG⋆ (· | X) − fG (· | X)∥L 1 for all X. For any j = 1, . . . , K ⋆ and i = 1, . . . , K ⋆ , define ! ⋆ ) ⋆ ))2 (Y − E(X, β exp(α 1 0j j ⋆ , f˜j⋆ (Y | X) := q exp − + X ⊤ α1j p̃⋆j (X) := PK ⋆ ⋆2 ⋆ + X ⊤ α⋆ ) 2σ ⋆2 exp(α j 2πσj ℓ=1 1ℓ 0ℓ exp(α0i ) (Y − E(X, βi ))2 1 ⊤ ˜ p̃i (X) := PK ⋆ exp − + X α1i , , fi (Y | X) := q ⊤ 2σi2 2πσi2 ℓ=1 exp(α0ℓ + X α1ℓ ) so that ⋆

K X

Z ∥fG⋆ (· | X) − fG (· | X)∥L 1 =





p̃⋆j (X)f˜j⋆ (Y | X) −

R j=1 K⋆

≤

X Z

f˜j⋆ (Y | X)dY

≲

i∈Cj

p̃⋆j (X) −

p̃i (X) +

Z p̃i (X)

f˜j⋆ (Y | X) − f˜i (Y | X) dY

j=1 i∈Cj

i∈Cj

⋆

K X

p̃i (X)f˜i (Y | X) dY

K⋆

j=1

⋆

p̃⋆j (X) −

j=1

p̃i (X) +

K X X

⋆ exp(α0i ) ∥α1i − α1j ∥ + ∥βi − βj⋆ ∥ + |σi2 − σj⋆2 | ,

j=1 i∈Cj

i∈Cj

P ⋆ ⊤ where the last inequality follows from the boundedness of 1/ K ℓ=1 exp(α0ℓ + X α1ℓ ) and a first R order Taylor expansion of f˜⋆ (Y | X) − f˜i (Y | X) dY around (α1i , βi , σ 2 ) = (α⋆ , β ⋆ , σ ⋆2 ). As for R

the first term, the triangle inequality yields ⋆

K X

⋆

p̃⋆j (X) −

j=1

p̃i (X) ≲

K X

⋆ exp(α0j )−

j=1

i∈Cj

exp(α0i )

i∈Cj

⋆

K X X j=1 i∈Cj

exp(α0i ) PK ⋆

⋆ ⊤ ⋆ ℓ=1 exp(α0ℓ + X α1ℓ )

− PK ⋆

ℓ=1 exp(α0ℓ + X

⊤α

1ℓ )

As for the last term above, using the Lipschitz property of t 7→ 1/t over a domain bounded away ⋆ ), one shows that is it is from 0 and an appropriate triangle inequality replacing X ⊤ α1ℓ with X ⊤ α1ℓ less than (up to a multiplicative constant) ⋆

K X j=1

⋆

⋆ exp(α0j )−

exp(α0i ) +

K X X

⋆ exp(α0i )∥α1i − α1j ∥.

j=1 i∈Cj

i∈Cj

Combining all of the above, we obtain ⋆

∥fG⋆ (· | X) − fG (· | X)∥L 1 ≲

K X X

⋆ exp(α0i ) − exp(α0j )

j=1 i∈Cj ⋆

K X X

⋆ exp(α0i ) ∥α1i − α1j ∥ + ∥βi − βj⋆ ∥ + |σi2 − σj⋆2 | ,

j=1 i∈Cj

as desired. 28

A.4

Proof of Theorem 4

Overview. Following from the result of Theorem 2, it suffices to prove that dH (gG , gG⋆ ) ≳ L1 (G, G⋆ ),

(18)

for any mixing measure G ∈ GK ⋆ . For that purpose, we first demonstrate that lim

inf

ε↘0 G∈GK ⋆ :L1 (G,G⋆ )≤ε

dH (gG , gG⋆ ) > 0. L1 (G, G⋆ )

(19)

Assume that the above equation holds for now. Then, there exists a positive constant ε′ such that inf G∈GK ⋆ :L1 (G,G⋆ )≤ε′ dH (gG , gG⋆ )/L1 (G, G⋆ ) > 0. As a consequence, we can complete the proof by showing that inf

G∈GK ⋆ :L1 (G,G⋆ )>ε′

dH (gG , gG⋆ ) > 0. L1 (G, G⋆ )

(20)

Given these arguments, we will provide the proofs for Equation (19) and Equation (20), respectively, in the sequel. Proof of Equation (19). Note that the Hellinger distance is bounded below by the L 1 -norm, that is, dH (g, g ′ ) ≳ ∥g − g ′ ∥L 1 . Therefore, it is sufficient to show that lim

inf

ε↘0 G∈GK ⋆ :L1 (G,G⋆ )≤ε

∥gG − gG⋆ ∥L 1 > 0. L1 (G, G⋆ )

(21)

We will use a proof-by-contradiction method to prove this result. In particular, assume that Equation (21) not hold true. Then, there exists a sequence of mixing measures (Gn ) of the P does ⋆ ⋆ form Gn = K exp(α n,0i )δ(αn,1i ,βn,i ,σ 2 ) that satisfies L1 (Gn , G ) → 0 and i=1 n,i

∥gGn − gG⋆ ∥L 1 /L1 (Gn , G⋆ ) → 0,

(22)

n over the range as n → ∞. Note that there are only a finite number of distinct sets C1n × . . . × CK ⋆ of n ∈ N. Thus, up to replacing Gn by its subsequence, we may assume without loss of generality that Cj = Cjn does not change with n for all j ∈ [K ⋆ ]. Additionally, since K ∗ is known under the exactly-specified setting and L1 (Gn , G⋆ ) → 0 as n → ∞, each Voronoi cell Cj has only one element for any j ∈ [K ∗ ]. Without loss of generality, we assume that Cj = {j} for simplicity. Thus, the Voronoi loss L1n := L1 (Gn , G⋆ ) can be written as ⋆

L1n =

K X

⋆ exp(αn,0j ) − exp(α0j )

j=1 ⋆

K X

⋆ ⋆ 2 exp(αn,0j )(∥αn,1j − α1j ∥ + ∥βn,j − β1j ∥ + |σn,j − σj⋆2 |),

(23)

j=1 ⋆ ) and (α 2 ⋆ ⋆ ⋆2 Since L1n → 0 as n → ∞, we have exp(αn,0j ) → exp(α0j n,1j , βn,j , σn,j ) → (α1j , βj , σj ) as ⋆ n → ∞, for all j ∈ [K ]. Next, we divide the rest of this proof into three main steps.

Step 1 - Density Decomposition. In this step, we decompose the density discrepancy gGn (Y, X) − gG⋆ (Y, X) = [fGn (Y | X) − fG⋆ (Y | X)]p(X) into a combination of linearly independent terms. In particular, let us denote   K⋆ X ⋆ ⋆ Dn :=  exp(α0j + X ⊤ (α1j )) · [gGn (Y, X) − gG⋆ (Y, X)], j=1

F (Y | X; α1 , β, σ ) := exp(X ⊤ α1 )N (Y | E(X, β), σ 2 ), 2

Hn (Y | X; α1 ) := exp(X ⊤ α1 )fGn (Y | X). Then, we can decompose Dn as Dn = An − Bn + Cn , where ⋆

An :=

K X

2 ⋆ exp(αn,0j )[F (Y | X; αn,1j , βn,j , σn,j ) − F (Y | X; α1j , βj⋆ , σj⋆2 )]p(X),

j=1 ⋆

Bn :=

K X

⋆ exp(αn,0j )[Hn (Y | X; αn,1j ) − Hn (Y | X; α1j )]p(X),

j=1 ⋆

Cn :=

K X

⋆ ⋆ ⋆ exp(αn,1j ) − exp(α1j ) [F (Y | X; α1j , βj⋆ , σj⋆2 ) − Hn (Y | X; α1j )]p(X).

j=1 ⋆ , ∆β ⋆ Next, we will further decompose An and Bn . Let ∆αn,1j := αn,1j − α1j n,j := βn,j − βj , and 2 := σ 2 − σ ⋆2 . By means of a first-order Taylor expansion, we have ∆σn,j n,j j ⋆

An =

K X

exp(αn,0j )

j=1

X 1 2 ρ3 (∆αn,1j )ρ1 (∆βn,j )ρ2 (∆σn,j ) ρ!

ρ:|ρ|=1

∂ |ρ| F (Y | E(X, βj⋆ ), σj⋆2 )p(X) + Rn,1 (Y | X)p(X). ∂α1ρ1 ∂β ρ2 ∂(σ 2 )ρ3 ¯

Above, we denote ρ = (ρ1 , ρ2 , ρ3 ) ∈ Nd × Np × N. Additionally, for any vectors a = (a1 , . . . , ad¯) ∈ Rd ¯ b and b = (b1 , . . . , bd¯) ∈ Nd , we let ab = ab11 . . . ad¯d¯, |a| = a1 + . . . + ad¯, and b! = b1 ! . . . bd¯!. Lastly, Rn,1 (Y | X)is a Taylor remainder such that Rn,1 (Y | X)/L1n → 0 as n → ∞. By taking the first-order partial derivatives of F (Y | X; α1 , β, σ 2 ) with respect to its parameters, we have ∂F (u) ∂α1

(Y | X; α1 , β, σ 2 ) = X (u) exp(X ⊤ α1 )N (Y | E(X, β), σ 2 ),

∂F ∂E ∂N 2 ⊤ (Y | E(X, β), σ 2 ), ′ ) (Y | X; α1 , β, σ ) = ′ ) (X, β) exp(X α1 ) (u (u ∂E ∂β ∂β ∂F 1 ∂2N 2 ⊤ (Y | X; α , β, σ ) = exp(X α ) (Y | E(X, β), σ 2 ), 1 1 ∂σ 2 2 ∂E 2 for all u ∈ [d] and u′ ∈ [d′ ]. Based on these derivatives, we can rewrite An,1 as ⋆

An =

K X 2 X j=1 γ=0

A(j) n,γ

∂γ N (Y | E(X, β), σ 2 )p(X) + Rn,1 (Y | X)p(X), ∂E γ 30

(24)

where (j) An,0 := exp(αn,0j )

d X

⋆ (∆αn,0ij )(u) X (u) exp(X ⊤ α1j ),

u=1 ′

d X

(j) An,1 := exp(αn,0j )

′

(∆βn,j )(u )

u′ =1

∂E ⋆ (X, βj⋆ ) exp(X ⊤ α1j ), ∂β (u′ )

1 (j) 2 ⋆ An,2 := exp(αn,0j ) (∆σn,j ) exp(X ⊤ α1j ), 2 for all j ∈ [K ⋆ ]. Next, by applying a first-order Taylor expansion, we rewrite Bn as ⋆

Bn =

K X j=1

d X ⋆ exp(αn,0j ) (∆αn,1j )(u) X (u) Hn (Y | X; α1j ) + Rn,2 (Y | X), u=1

where Rn,2 (Y | X) is a Taylor remainder such that Rn,2 (Y | X)/L1n → 0 as n → ∞. Step 2 - Non-vanishing coefficients. In this step, we prove by contradiction that at least one among the coefficients in the representations of [An − Rn,1 (Y | X)]p(X)/L1n , [Bn − Rn,2 (Y | X)]p(X)/L1n , Cn /L1n does not go to zero as n → ∞. Assume that all these coefficients converge to zero. From the coefficients of the terms: ⋆ )N (Y |X; α⋆ , β ⋆ , σ ⋆2 ) for j ∈ [K ⋆ ], we get • exp(X ⊤ α1j 1j j j K⋆

1 X ⋆ · exp(αn,0j ) − exp(α1j ) → 0; L1n j=1

⋆ )N (Y |X; α⋆ , β ⋆ , σ ⋆2 ) for j ∈ [K ⋆ ] and u ∈ [d], we get • X (u) exp(X ⊤ α1j 1j j j K⋆

1 X · exp(αn,0j )∥∆αn,1j ∥ → 0; L1n j=1

⋆ ′ ′ ⊤ ⋆ ∂N ⋆ ⋆ ⋆2 • ∂β∂E (u′ ) exp(X α1j ) ∂E (Y |X; α1j , βj , σj ) for j ∈ [K ] and u ∈ [d ], we get K⋆

1 X · exp(αn,0j )∥∆βn,j ∥ → 0; L1n j=1

⋆ ) ∂ N (Y |X; α⋆ , β ⋆ , σ ⋆2 ) for j ∈ [K ⋆ ], we get • exp(X ⊤ α1j j 1j j ∂E 2 K⋆

1 X 2 · exp(αn,0j )|∆σn,j | → 0; L1n j=1

By taking the sum of these limits, we have 1 = L11n · L1n → 0 as n → ∞, which is a contradiction. Therefore, not all the coefficients in the representations of [An −Rn,1 (Y | X)]p(X)/L1n , [Bn −Rn,2 (Y | X)]p(X)/L1n , Cn /L1n converge to zero as n → ∞. Step 3 - Fatou’s argument. In this step, we will point out a contradiction to the results of Step 2. To begin, we denote by mn the maximum of the absolute values of the coefficients in the representations of [An − Rn,1 (Y | X)]p(X)/L1n , [Bn − Rn,2 (Y | X)]p(X)/L1n , Cn /L1n . Then, the results of Step 2 indicate that 1/mn ̸→ ∞ as n → ∞. Next, let us denote 1 ⋆ · exp(αn,0j ) − exp(α0j ) → t0,j , mn L1n 1 ′ (u′ ) · exp(αn,0j )(∆βn,j )(u ) → t2,j , mn L1n

1 (u) · exp(αn,0j )(∆αn,1j )(u) → t1,j , mn L1n 1 2 · exp(αn,0j )(∆σn,j ) → t3,j , mn L1n

for all j ∈ [K ⋆ ] as n → ∞. It follows from Step 2 that not all these limits are zero. Next, recall that we have ∥gGn − gG⋆ ∥L 1 /L1 (Gn , G⋆ ) → 0 as n → ∞. Then, by means of the Fatou’s lemma, we have Z ∥gGn − gG⋆ ∥L 1 |gGn (Y, X) − gG⋆ (Y, X)| lim ≥ lim inf d(Y, X), n→∞ n→∞ mn L1n mn L1n which implies that [gGn (Y, X) − gG⋆ (Y, X)]/[mn L1nP ] → 0 as n → ∞ along a subsequence (focus on ⋆ ⋆ ⊤ ⋆ such a subsequence from now on). Since the term K j=1 exp(α0j + X α1j ) is bounded, we deduce Dn /[mn L1n ] → 0 as n → ∞, for almost every (Y, X). From the decomposition of Dn , it follows that 1 · (An − Bn + Cn ) → 0. mn L1n

(25)

⋆ ) ∂ N (Y | X; α⋆ , β ⋆ , σ ⋆2 ) and H (Y | For ease of presentation, let us denote Fγ,j := exp(X ⊤ α1j j 1j j j ∂E γ ⋆ X) := limn→∞ Hn (Y | X; α1j ), for all 0 ≤ γ ≤ 2 and j ∈ [K ⋆ ]. Then, we have ′

⋆

K hX d d X X An (u) (u′ ) ∂E lim = t1,j X (u) F0,j (Y |X) + t2,j (X, βj⋆ )F1,j (Y |X) (u′ ) n→∞ mn L1n ∂β j=1 u=1 u′ =1 i 1 + t3,j F2,j (Y |X) p(X), 2 K⋆ X d X Bn (u) lim = t1,j X (u) Hj (Y |X)p(X), n→∞ mn L1n j=1 u=1 K⋆

X Cn lim = t0,j [F0,j (Y |X) − Hj (Y |X)]p(X). n→∞ mn L1n j=1

n o Note that for almost every X, the set Fγ,j (Y |X), Hj (Y |X) : 0 ≤ γ ≤ 2, j ∈ [K ⋆ ] is linearly independent with respect to Y . Therefore, it follows that the coefficients of these terms in the limit in Equation (25) become zero. For j ∈ [K ⋆ ], by considering the coefficients of P (u) • F0,j (Y |X), we have t0,j + du=1 t1,j X (u) = 0, for almost every X. Then, we deduce t0,j = (u)

t1,j = 0 for all u ∈ [d]; 32

(u′ ) ∂E ⋆ u′ =1 t2,j ∂β (u′ ) (X, βj ), for almost every X. As the expert function E is (u′ ) first-order strongly identifiable, we get t2,j = 0 for all u′ ∈ [d′ ];

• F1,j (Y |X), we have

Pd′

• F2,j (Y |X), we have t3,j = 0. (u)

(u′ )

Putting the above results together, we have t0,j = t1,j = t2,j = t3,j = 0 for all j ∈ [K ⋆ ], u, v ∈ [d] and u′ , v ′ ∈ [d′ ]. This contradicts to the fact that at least one among them is non-zero. Therefore, we achieve the result in Equation (19). Proof of Equation (20). Assume by contrary that Equation (20) does not hold. Then, we can find a sequence of mixing measure (Gn ) such that L1 (Gn , G⋆ ) > ε′ and ∥gGn −gG⋆ ∥L 1 /L1 (Gn , G⋆ ) → 0 as n → ∞. These two properties imply that ∥gGn − gG⋆ ∥L 1 → 0. Since the parameter space Θ is compact, we can substitute the sequence (Gn ) with its subsequence (G′n ) that converges to some mixing measure G′ . Recall that L1 (Gn , G⋆ ) > ε′ , then we also have L1 (G′ , G⋆ ) ≥ ε′ . On the other hand, by the Fatou’s lemma, we get Z 0 = lim ∥gG′n − gG⋆ ∥L 1 ≥ lim inf |gG′n (Y, X) − gG⋆ (Y, X)|d(Y, X) n→∞ n→∞ Z = |gG′ (Y, X) − gG⋆ (Y, X)|d(Y, X). This inequality indicates that gG′ (Y, X) = gG⋆ (Y, X) for almost every (Y, X). Since the SMoGE model is identifiable, we deduce G′ ≡ G⋆ . As a result, we get L1 (G′ , G⋆ ) = 0, which contradicts the previous result that L1 (G′ , G⋆ ) > ε′ > 0. Hence, the proof is completed.

A.5

Proof of Theorem 5

Overview. Following the result of Theorem 3 and using the same arguments for the proof of Theorem 4 in Appendix A.4, it suffices to establish the following inequality ∥gG − gG⋆ ∥L 1 > 0. ε↘0 G∈OK :L2 (G,G⋆ )≤ε L2 (G, G⋆ ) lim

inf

(26)

Assume by contrary that Equation (26) is not true, that is, we can find a sequence (Gn ) such that L2 (Gn , G⋆ ) → 0 and ∥gGn − gG⋆ ∥L 1 /L2 (Gn , G⋆ ) → 0, as n → ∞. Since Kn ≤ K for all n, there exists a subsequence of Gn such that Kn does not change with n. Therefore, up to replacing Gn by this subsequence, we may assume that Kn = K ′ ≤ K for all n. Similarly, since there are only a n over the range of n ∈ N, we may assume without loss of finite number of distinct sets C1n × . . . × CK ⋆ n generality that Cj = Cj does not change with n for all j ∈ [K ⋆ ]. Then, we can rewrite the Voronoi loss L2n := L2 (Gn , G⋆ ) as ⋆

L2n =

K X X

⋆ exp(αn,0i ) − exp(α0j )

j=1 i∈Cj

⋆ ⋆ 2 exp(αn,0i )(∥αn,1i − α1j ∥ + ∥βn,i − β1j ∥ + |σn,i − σj⋆2 |)

j∈[K ⋆ ]:|Cj |=1 i∈Cj

⋆ 2 2 ⋆ 2 exp(αn,0i )(∥αn,1i − α1j ∥ + ∥βn,i − β1j ∥ + |σn,i − σj⋆2 |2 ),

j∈[K ⋆ ]:|Cj |>1 i∈Cj

(27)

P ⋆ Recall that we have L2n → 0 as n → ∞, which implies i∈Cj exp(αn,0i ) → exp(α0j ) and 2 ) → (α⋆ , β ⋆ , σ ⋆2 ) as n → ∞, for all i ∈ C and j ∈ [K ⋆ ]. Subsequently, we sep(αn,1i , βn,i , σn,i j 1j j j arate the rest into three main steps. Step 1 - Density Decomposition. First, we decompose the density discrepancy gGn (Y, X) − gG⋆ (Y, X) = [fGn (Y | X) − fG⋆ (Y | X)]p(X) into a combination of linearly independent terms through the quantity   K⋆ X ⋆ ⋆ Dn :=  exp(α0j + X ⊤ (α1j )) · [gGn (Y, X) − gG⋆ (Y, X)]. j=1

In particular, let F (Y | X; α1 , β, σ 2 ) := exp(X ⊤ α1 )N (Y | E(X, β), σ 2 ) and Hn (Y | X; α1 ) := exp(X ⊤ α1 )fGn (Y | X). Then, the quantity Dn can be represented as Dn = An − Bn + Cn , where we define ⋆

An :=

K X X

2 ⋆ exp(αn,0i )[F (Y | X; αn,1i , βn,i , σn,i ) − F (Y | X; α1j , βj⋆ , σj⋆2 )]p(X),

j=1 i∈Cj ⋆

Bn :=

K X X

⋆ exp(αn,0i )[Hn (Y | X; αn,1i ) − Hn (Y | X; α1j )]p(X),

j=1 i∈Cj

  K⋆ X X ⋆  ⋆ ⋆  ) [F (Y | X; α1j , βj⋆ , σj⋆2 ) − Hn (Y | X; α1j )]p(X). Cn := exp(αn,1i ) − exp(α1j j=1

i∈Cj

Next, we rewrite An as An = An,1 + An,2 , where X X 2 ⋆ ⋆ An,1 = exp(αn,0i )[F (Y | X; αn,1i , βn,i , σn,i ) − F (Y | X; α1j , β0j , βj⋆ , σj⋆2 )]p(X), j∈[K ⋆ ]:|Cj |=1 i∈Cj

An,2 =

2 ⋆ ⋆ exp(αn,0i )[F (Y | X; αn,1i , βn,i , σn,i ) − F (Y | X; α1j , β0j , βj⋆ , σj⋆2 )]p(X).

j∈[K ⋆ ]:|Cj |>1 i∈Cj ⋆ , ∆β ⋆ 2 2 ⋆2 Let us denote ∆αn,1ij := αn,1i − α1j n,ij := βn,i − βj , and ∆σn,ij := σn,i − σj . By applying first-order and second-order Taylor expansions, we get

An,1 =

exp(αn,0i )

j∈[K ⋆ ]:|Cj |=1 i∈Cj

ρ:|ρ|=1

∂ |ρ| F

(Y | E(X, βj⋆ ), σj⋆2 )p(X) + Rn,1 (Y | X)p(X), ∂α1ρ1 ∂β ρ2 ∂(σ 2 )ρ3 X X 1 2 exp(αn,0i ) (∆αn,1ij )ρ1 (∆βn,ij )ρ2 (∆σn,ij )ρ3 ρ! ×

An,2 =

X 1 2 (∆αn,1ij )ρ1 (∆βn,ij )ρ2 (∆σn,ij )ρ3 ρ!

j∈[K ⋆ ]:|Cj |>1 i∈Cj

ρ:|ρ|=2

∂ |ρ| F ∂α1ρ1 ∂β ρ2 ∂(σ 2 )ρ3

(Y | E(X, βj⋆ ), σj⋆2 )p(X) + Rn,2 (Y | X)p(X), ¯

Rn,ℓ (Y | X), for ℓ ∈ {1, 2}, are Taylor remainders such that Rn,ℓ (Y | X)/L2n → 0 as n → ∞. Recall that the first-order partial derivatives of F (Y | X; α1 , β, σ 2 ) with respect to its parameters have been provided in Equation (24). Meanwhile, the second-order partial derivatives of F (Y | X; α1 , β, σ 2 ) with respect to its parameters are given by ∂2F (u) (v) ∂α1 ∂α1 ∂2F

(Y | X; α1 , β, σ 2 ) = X (u) X (v) exp(X ⊤ α1 )N (Y | E(X, β), σ 2 ),

∂2E ∂N ⊤ (Y | E(X, β), σ 2 ) ′) ′) ′ ) (X, β) exp(X α1 ) (u (u (v ∂E ∂β ∂β ∂β ∂β ∂E ∂2N ∂E ⊤ (X, β) (X, β) exp(X α ) (Y | E(X, β), σ 2 ), + 1 ′ ′ ∂E 2 ∂β (u ) ∂β (v ) ∂2F 1 ∂4N 2 ⊤ (Y | X; α , β, σ ) = exp(X α ) (Y | E(X, β), σ 2 ), 1 1 ∂ν 2 4 ∂E 4 (Y | X; α1 , β, σ 2 ) = (v ′ )

and ∂2F (u) ∂α1 ∂β (v′ ) ∂2F

(Y | X; α1 , β, σ 2 ) = X (u)

∂N ∂E ⊤ (Y | E(X, β), σ 2 ), ′ ) (X, β) exp(X α1 ) (v ∂E ∂β

∂2N 1 (u) ⊤ 2 X exp(X α ) (Y | E(X, β), σ 2 ), (Y | X; α , β, σ ) = 1 1 2 (u) 2 ∂E ∂α ∂ν 1

∂2F 1 ∂E ∂3N 2 ⊤ (Y | X; α , β, σ ) = (X, β) exp(X α ) (Y | E(X, β), σ 2 ), 1 1 2 ∂β (u′ ) ∂E 3 ∂β (u′ ) ∂ν for all i ∈ [d] and u′ ∈ [d′ ]. From these derivatives, An,1 can be written as 2 X

An,1 =

(j)

An,1,γ

j∈[K ⋆ ]:|Cj |=1 γ=0

∂γ N (Y | E(X, β), σ 2 )p(X) + Rn,1 (Y | X)p(X), ∂E γ

where (j) An,1,0 :=

exp(αn,0i )

d X

⋆ (∆αn,0ij )(u) X (u) exp(X ⊤ α1j ),

i∈Cj

u=1

d X

′

(j) An,1,1 :=

exp(αn,0i )

u′ =1

i∈Cj (j)

An,1,2 :=

X i∈Cj

′

(∆βn,ij )(u )

∂E ⋆ (X, βj⋆ ) exp(X ⊤ α1j ), ∂β (u′ )

1 2 ⋆ exp(αn,0i ) (∆σn,ij ) exp(X ⊤ α1j ), 2

for all j ∈ [K ⋆ ] : |Cj | = 1. Analogously, we can write An,2 as An,2 =

4 X

j∈[K ⋆ ]:|Cj |>1 γ=0

(j)

An,1,γ

∂γ N (Y | E(X, β), σ 2 )p(X) + Rn,2 (Y | X)p(X), ∂E γ 35

where " d # d X X (∆αn,1ij )(u) (∆αn,1ij )(v) (u) (v) (u) (u) ⋆ exp(αn,0i ) (∆αn,1ij ) X + X X exp(X ⊤ α1j ), 1 + 1{u=v} u=1 u,v=1 i∈Cj " d′ X X ∂E ′ (j) An,2,1 := exp(αn,0i ) (∆βn,ij )(u ) (u′ ) (X, βj⋆ ) ∂β ′ i∈C

(j) An,2,0 :=

u =1

′

′ ′ d X (∆βn,ij )(u ) (∆βn,ij )(v ) ∂2E (X, βj⋆ ) + (u′ ) ∂β (v ′ ) ′ =v ′ } 1 + 1 ∂β {u u′ ,v ′ =1 # d X d′ X (u) (v ′ ) (u) ∂E ⋆ ⋆ + (∆αn,1ij ) (∆βn,ij ) X (X, βj ) exp(X ⊤ α1j ), (v ′ ) ∂β u=1 ′

v =1

′

′ ′ d X (∆βn,ij )(u ) (∆βn,ij )(v ) ∂E 1 ∂E 2 ⋆ exp(αn,0i ) (∆σn,ij ) + (X, βj⋆ ) ′ ) (X, βj ) (u (v ′ ) ′ =v ′ } 2 1 + 1 ∂β ∂β {u i∈Cj u′ ,v ′ =1 # d X 1 (u) 2 (u) ⋆ (∆αn,1ij ) (∆σn,ij )X exp(X ⊤ α1j ), + 2

(j) An,2,2 :=

u=1

′

(j)

An,2,3 :=

exp(αn,0i )

u′ =1

i∈Cj (j)

An,2,4 :=

X i∈Cj

d X 1

′

2 (∆βn,ij )(u ) (∆σn,ij )

∂E ⋆ (X, βj⋆ ) exp(X ⊤ α1j ), ∂β (u′ )

1 2 ⋆ )2 exp(X ⊤ α1j ), exp(αn,0i ) (∆σn,ij 8

for all j ∈ [K ⋆ ] : |Cj | > 1. Next, we rewrite Bn as Bn = Bn,1 + Bn,2 , where X X ⋆ Bn,1 := exp(αn,0i )[Hn (Y | X; αn,1i ) − Hn (Y | X; α1j )]p(X), j∈[K ⋆ ]:|Cj |=1 i∈Cj

Bn,2 :=

⋆ exp(αn,0i )[Hn (Y | X; αn,1i ) − Hn (Y | X; α1j )]p(X).

j∈[K ⋆ ]:|Cj |>1 i∈Cj

By applying first-order and second-order Taylor expansions, we have Bn,1 =

exp(αn,0i )

j∈[K ⋆ ]:|Cj |=1 i∈Cj

Bn,2 =

⋆ (∆αn,1ij )(u) X (u) Hn (Y | X; α1j ) + Rn,3 (Y | X),

u=1 d hX ⋆ exp(αn,0i ) (∆αn,1ij )(u) X (u) Hn (Y | X; α1j ),

j∈[K ⋆ ]:|Cj |>1 i∈Cj

d X

u=1

d X (∆αn,1ij )(u) (∆αn,1ij )(v) u,v=1

1 + 1{u=v}

i ⋆ X (u) X (v) Hn (Y | X; α1j ) + Rn,4 (Y | X),

where Rn,ℓ (Y | X), for ℓ ∈ {3, 4}, are Taylor remainders such that Rn,ℓ (Y | X)/L2n → 0 as n → ∞. 36

Step 2 - Non-vanishing coefficients. In this step, we prove by contradiction that at least one among the coefficients in the representations of [An,1 − Rn,1 (Y | X)]p(X)/L2n , [An,2 − Rn,2 (Y | X)]p(X)/L2n , [Bn,1 − Rn,3 (Y | X)]p(X)/L2n , [Bn,2 − Rn,4 (Y | X)]p(X)/L2n , Cn /L2n does not go to zero as n → ∞. Assume that all these coefficients converge to zero. From the coefficients of the terms: ⋆ )N (Y |X; α⋆ , β ⋆ , σ ⋆2 ) for j ∈ [K ⋆ ], we get • exp(X ⊤ α1j 1j j j ⋆

K 1 X X ⋆ · exp(αn,0i ) − exp(α1j ) → 0; L2n j=1 i∈Cj

⋆ )N (Y |X; α⋆ , β ⋆ , σ ⋆2 ) for j ∈ [K ⋆ ] : |C | = 1 and u ∈ [d], we get • X (u) exp(X ⊤ α1j j 1j j j

1 · L2n

exp(αn,0i )∥∆αn,1ij ∥ → 0;

j∈[K ⋆ ]:|Cj |=1 i∈Cj

⋆ ⋆ ⋆2 ⋆ ′ ′ ⊤ ⋆ ∂N • ∂β∂E (u′ ) exp(X α1j ) ∂E (Y |X; α1j , βj , σj ) for j ∈ [K ] : |Cj | = 1 and u ∈ [d ], we get

1 · L2n

exp(αn,0i )∥∆βn,ij ∥ → 0;

j∈[K ⋆ ]:|Cj |=1 i∈Cj

⋆ ) ∂ N (Y |X; α⋆ , β ⋆ , σ ⋆2 ) for j ∈ [K ⋆ ] : |C | = 1, we get • exp(X ⊤ α1j j 1j j j ∂E 2

1 · L2n

2 exp(αn,0i )|∆σn,ij | → 0;

j∈[K ⋆ ]:|Cj |=1 i∈Cj

⋆ )N (Y |X; α⋆ , β ⋆ , σ ⋆2 ) for j ∈ [K ⋆ ] : |C | > 1 and u, v ∈ [d], we get • X (u) X (v) exp(X ⊤ α1j j 1j j j

1 · L2n •

∂2E

′ ∂β (u )

(X, βj⋆ )

exp(αn,0i )∥∆αn,1ij ∥2 → 0;

j∈[K ⋆ ]:|Cj |>1 i∈Cj

⋆ ) ∂ N (Y |X; α⋆ , β ⋆ , σ ⋆2 ) for j ∈ [K ⋆ ] : |C | > 1 and u′ ∈ [d′ ], we get exp(X ⊤ α1j j 1j j j ∂E 2

1 · L2n

exp(αn,0i )∥∆βn,ij ∥ → 0;

j∈[K ⋆ ]:|Cj |>1 i∈Cj

⋆ ) ∂ N (Y |X; α⋆ , β ⋆ , σ ⋆2 ) for j ∈ [K ⋆ ] : |C | > 1, we get • exp(X ⊤ α1j j j j 1j ∂E 4

1 · L2n

j∈[K ⋆ ]:|Cj |>1 i∈Cj

2 exp(αn,0i )|∆σn,ij |2 → 0;

By taking the sum of these limits, we have 1 = L12n · L2n → 0 as n → ∞, which is a contradiction. Therefore, not all the coefficients in the representations of [An,1 − Rn,1 (Y | X)]p(X)/L2n , [An,2 − Rn,2 (Y | X)]p(X)/L2n , [Bn,1 − Rn,3 (Y | X)]p(X)/L2n , [Bn,2 − Rn,4 (Y | X)]p(X)/L2n , Cn /L2n converge to zero as n → ∞. Step 3 - Fatou’s argument. In this step, we will show a contradiction to the results of Step 2. To begin with, we denote by mn the maximum of the absolute values of the coefficients in the representations of [An,1 − Rn,1 (Y | X)]p(X)/L2n , [An,2 − Rn,2 (Y | X)]p(X)/L2n , [Bn,1 − Rn,3 (Y | X)]p(X)/L2n , [Bn,2 − Rn,4 (Y | X)]p(X)/L2n , Cn /L2n . Then, the results of Step 2 indicate that 1/mn ̸→ ∞ as n → ∞. Next, let us denote X 1 ⋆ · exp(αn,0i ) − exp(α0j ) → t0,j , mn L2n i∈Cj

1 · mn L2n 1 · mn L2n

′

(u )

exp(αn,0i )(∆βn,ij )(u ) → t2,j ,

i∈Cj

X 1 2 exp(αn,0i )(∆σn,ij · ) → t3,j , mn L2n i∈Cj

(uv)

exp(αn,0i )(∆αn,1ij )(u) (∆αn,1ij )(v) → t4,j ,

i∈Cj

X ′ ′ 1 (u′ v ′ ) · exp(αn,0i )(∆βn,ij )(u ) (∆βn,ij )(v ) → t5,j , mn L2n i∈Cj

1 · mn L2n 1 · mn L2n

X 1 (u) · exp(αn,0i )(∆αn,1ij )(u) → t1,j , mn L2n

2 exp(αn,0i )(∆σn,ij )2 → t6,j ,

i∈Cj

X ′ 1 (uv ′ ) · exp(αn,0i )(∆αn,1ij )(u) (∆βn,ij )(v ) → t7,j , mn L2n i∈Cj

(u)

2 exp(αn,0i )(∆αn,1ij )(u) (∆σn,ij ) → t8,j ,

i∈Cj

X ′ 1 (u′ ) 2 · exp(αn,0i )(∆βn,ij )(u ) (∆σn,ij ) → t9,j , mn L2n i∈Cj

for all j ∈ [K ⋆ ] as n → ∞. The results of Step 2 indicate that not all these limits are zero. According to the Fatou’s lemma, we have

∥gGn − gG⋆ ∥L 1 ≥ n→∞ mn L2n lim

|gGn (Y, X) − gG⋆ (Y, X)| d(Y, X), n→∞ mn L2n

lim inf

Since ∥gGn − gG⋆ ∥L 1 /L2 (Gn , G⋆ ) → 0 as n → ∞, we get [gGn (Y, X) − gG⋆ (Y, X)]/[mn L2n ] → 0 as n → ∞ along a subsequence every (Y, X) (work on this subsequence from now on). P ⋆ for almost ⋆ + X ⊤ α⋆ ) is bounded, it follows that D /[m L ] → 0 as Furthermore, as the term K exp(α n n 2n 0j 1j j=1 n → ∞, for almost every (Y, X). From the decomposition of Dn in Step 1, we have

1 · [An,1 + An,2 − Bn,1 − Bn,2 + Cn ] → 0. mn L2n

(28)

⋆ ) ∂ N (Y | X; α⋆ , β ⋆ , σ ⋆2 ) and H (Y | X) := For ease of presentation, let Fγ,j := exp(X ⊤ α1j j 1j j j ∂E γ

⋆ ), for all 0 ≤ γ ≤ 4 and j ∈ [K ⋆ ]. Then, we have limn→∞ Hn (Y | X; α1j

An,2 n→∞ mn L2n lim

d′ X

d hX

′

An,1 = lim n→∞ mn L2n

j∈[K ⋆ ]:|Cj |>1

(u) t1,j X (u) F0,j (Y |X) +

u=1

u,v=1

′

(u′ ) t2,j

u′ =1

d X ∂E ∂2E (u′ v ′ ) ⋆ (X, β ) + (X, βj⋆ ) t ′ j 5,j (u′ ) ∂β (v ′ ) ∂β (u ) ∂β ′ ′ u ,v =1

d X d′ X ∂E (uv ′ ) + t7,j X (u) (v′ ) (X, βj⋆ ) F1,j (Y |X) ∂β u=1 ′ v =1

1 2

t3,j +

(u′ v ′ ) t5,j

u′ ,v ′ =1 d′

d′ X

X 1 u′ =1

(u′ )

t9,j

d X ∂E ∂E 1 (u) (u) ⋆ ⋆ F2,j (Y |X) t X (X, β ) (X, β ) + ′ ′ j j 2 8,j ∂β (u ) ∂β (v ) u=1

i ∂E 1 ⋆ (X, β ) F (Y |X) + t F (Y |X) p(X), 3,j 6,j 4,j j 8 ∂β (u′ )

and Bn,1 = n→∞ mn L2n

lim

Bn,2 lim = n→∞ mn L2n

d X (u) t1,j X (u) Hj (Y |X)p(X),

j∈[K ⋆ ]:|Cj |=1 u=1

d X

j∈[K ⋆ ]:|Cj |>1

u=1

(u) t1,j X (u) +

d X

(uv) t4,j X (u) X (v) Hj (Y |X)p(X),

u,v=1

K⋆

X Cn = t0,j [F0,j (Y |X) − Hj (Y |X)]p(X). n→∞ mn L2n lim

j=1

n o Note that for almost every X, the set Fγ,j (Y |X), Hj (Y |X) : 0 ≤ γ ≤ 4, j ∈ [K ⋆ ] is linearly independent with respect to Y . Therefore, it follows that the coefficients of these terms in the limit in Equation (28) become zero. For j ∈ [K ⋆ ] such that |Cj | = 1, by considering the coefficients of P (u) • F0,j (Y |X), we have t0,j + du=1 t1,j X (u) = 0, for almost every X. Then, we deduce t0,j = (u)

t1,j = 0 for all u ∈ [d]; (u′ ) ∂E ⋆ u′ =1 t2,j ∂β (u′ ) (X, βj ), for almost every X. As the expert function E is (u′ ) second-order strongly identifiable, we get t2,j = 0 for all u′ ∈ [d′ ];

• F1,j (Y |X), we have

Pd′

• F2,j (Y |X), we have t3,j = 0. 39

For j ∈ [K ⋆ ] such that |Cj | > 1, by considering the coefficients of • F0,j (Y |X), we have t0,j + (u)

(u) (u) Pd (uv) + u,v=1 t4,j X (u) X (v) = 0, for almost every X. Then, u=1 t1,j X

(uv)

we get t0,j = t1,j = t4,j for all u, v ∈ [d]. • F1,j (Y |X), we have ′

′

d d X X ∂2E (u′ ) ∂E (u′ v ′ ) ⋆ t2,j (X, β ) + (X, βj⋆ ) t ′ j 5,j (u ) (u′ ) ∂β (v ′ ) ∂β ∂β ′ ′ ′

u =1

u ,v =1

d X d′ X

(uv ′ )

t7,j X (u)

u=1 v ′ =1

∂E (X, βj⋆ ) = 0, ∂β (v′ )

for almost every X. As the expert function E meets the second-order strong identifiability (u′ ) (u′ v ′ ) (uv ′ ) condition, we get t2,j = t5,j = t7,j = 0 for all u′ , v ′ ∈ [d′ ] and u ∈ [d]; • F2,j (Y |X), we have ′

d d X X ∂E 1 (u) (u) 1 (u′ v ′ ) ∂E ⋆ ⋆ t3,j + t5,j t8,j X = 0, (X, β ) (X, β ) + j j (u′ ) (v ′ ) 2 2 ∂β ∂β ′ ′ u=1

u ,v =1

(u′ v ′ )

for almost every X. Since t5,j

= 0 for all u′ , v ′ ∈ [d′ ], it follows that 12 t3,j +

1 (u) (u) = u=1 2 t8,j X

(u)

0, for almost every X. Then, we get t3,j = t8,j = 0 for all u′ , v ′ ∈ [d′ ] and u ∈ [d]; ′ 1 (u ) ∂E ⋆ u′ =1 2 t9,j ∂β (u′ ) (X, βj ) = 0, for almost every X. As the expert function (u′ ) E is strongly identifiable, we get t9,j for all u′ ∈ [d′ ];

• F3,j (Y |X), we have

Pd′

• F4,j (Y |X), we have t6,j = 0. (u)

(u′ )

(uv)

(u′ v ′ )

Combining the above results, we deduce t0,j = t1,j = t2,j = t3,j = t4,j = t5,j (u) (u′ ) t8,j = t9,j = 0 for all j ∈ [K ⋆ ], u, v ∈ [d] and u′ , v ′ ∈ [d′ ].

(uv ′ )

= t6,j = t7,j

This contradicts the fact that at least one among them is non-zero. Therefore, we achieve the result in Equation (26). Hence, the proof is completed.

A.6

Proof of Theorem 6

Following the proof strategy of Miller (2023) for simple mixture models, we prove the result by leveraging the famous posterior consistency theorem by Doob (1949).7 However, the latter requires a stronger notion of identifiability with respect to SMoGE parameters. This, despite the identifiability of SMoGE densities with respect to mixing measures, does not hold on Θ∞ because mixing measures are invariant under (i) permutations of expert labels, and (ii) mergers and separations of equal mixing measure atoms. Therefore, we first construct a restricted space Θ̃∞ where the assumptions of Doob’s theorem (measurability of the model and parameter identifiability) hold; then we show that 7

Throughout the proof, we refer to the formulation of Doob’s theorem presented in Miller (2018) (Theorem 2.4), requiring measurability of the model with respect to parameters and parameter identifiability.

consistency on that space implies consistency on the original space Θ∞ , thanks to the assumptions on the prior. Before delving into the proof, here is some useful notation. Given any θ ∈ Θk ⊆ Θ∞ and any permutation ρ : {1, . . . , k} → {1, . . . , k}, define θ[ρ] := (θρ(1) , . . . , θρ(k) ). Moreover, recall the decomposition Θ = Ω × A, where Ω denotes the space of parameters ω = (α1 , β, σ 2 ), and A is the space of gating bias parameters α0 . Therefore, for all k ∈ N, Θk can be identified with Θk ≡ Ωk × Ak . Also, for all θ, θ′ ∈ Θ∞ , let ( min{∥θ − θ′ ∥2 , 1} if K(θ) = K(θ′ ), dΘ∞ (θ, θ ) := 1 otherwise. ′

By Propositions A.1 and A.2 of Miller (2023), this definition makes (Θ∞ , dΘ∞ ) a Borel measurable subset of a Polish space, which we endow with the corresponding Borel sigma-algebra. Construction of the restricted model. For all k ∈ N, define the restricted expert space Ωk := {(ω1 , . . . , ωk ) ∈ Ωk : ω1 ≺ · · · ≺ ωk }, where, for all ω, ω ′ ∈ Ω, ω ≺ ω ′ means that ω precedes ω ′ lexicographically, and ω ̸= ω ′ . Now S k define Θk := Ωk × A and Θ̃∞ := ℓ∈N Θℓ . For all k ∈ N and θ ∈ Θk , let T (θ) = θ[ρ], where the permutation ρ is chosen such that θ[ρ] ∈ Θk if possible (that is, if no two or more experts in θ share the same ω parameters), otherwise θ[ρ] = θ. Then, denoting θ ∼ Π, we obtain Π(T (θ) ∈ Θ̃∞ ) = 1 because the subset of Ωk where two or more experts share the same parameters has prior probability zero by the assumptions of the Theorem. Denoting B[ρ] = {θ ∈ Θk : θ[ρ] ∈ B} for all B ⊆ Θk , the definition of T implies [ T −1 (B) = {θ ∈ Θk : T (θ) ∈ B} = B[ρ], (29) ρ∈Rk

where Rk denotes the space of all permutations of {1, . . . , k}. Letting Q̃ denote the pushforward of Π through T (i.e., the distribution of T (θ) when θ ∼ Π) restricted to Θ̃∞ , we obtain the restricted model 8 iid (Xi , Yi ) | T (θ) ∼ gG(T (θ)) , i = 1, . . . , n (30) T (θ) ∼ Q̃ by Theorem 10.2.1 of Dudley (2002). This is because gG(θ) = gG(T (θ)) and, for all measurable C ⊆ (X × R)n and B ⊆ Θ̃∞ , the following holds: P((Xi , Yi )ni=1 ∈ C, T (θ) ∈ B) = P((Xi , Yi )ni=1 ∈ C, θ ∈ T −1 (B)) =

Z B

n gG(θ) (C) Q̃(dθ),

(n)

where measurability of θ 7→ gG(θ) (C) for all measurable C ⊆ (X × R)n follows from the measurability of θ 7→ gG(θ) (C ′ ) for all measurable C ′ ⊆ X × R (shown in the next paragraph) and from an application of Lemma 5.2 of Miller (2018). 8

We denote by P the joint probability measure over infinite data sequences and parameters implied by this model.

Proof of measurability. Let C ⊆ X × R be measurable and let CX denote its section at X ∈ X. Then for any k ∈ N and θ ∈ Θk , by Fubini’s theorem we can write Z Z fG(θ) (Y | X)dY p(X)dX. θ 7→ gG(θ) (C) = (31) X

By the compactness of X and the Gaussianity of the expert densities, fG(θ) (Y | X) can be bounded above by a function of Y (constant across X) that is integrable on R (e.g., a large enough multiple of a Gaussian density with high enough variance). So the dominated convergence theorem and the continuity of θ 7→ fG(θ) (Y | X) for all (X, Y ) ∈ X × R imply the continuity (hence measurability) of θ 7→ gG(θ) (C) as a function on Θk . Therefore, this mapping is also measurable as a function on Θk and, as a consequence, it is measurable as a function on Θ̃∞ (since the pre-image of a measurable subset of R is a union of measurable subsets of Θ1 , Θ2 , . . . , and is thus measurable by Proposition A.2 of Miller (2023)). Proof of parameter identifiability. Choose θ ∈ Θk ⊂ Θ̃∞ , θ′ ∈ Θk′ ⊂ Θ̃∞ such that gG(θ) = gG(θ′ ) , where θ = (α1j , βj , σj2 , α0j )kj=1 ≡ (ωj , α0j )kj=1 ′

′

′ ′ k ′ k )j=1 ≡ (ωj′ , α0j )j=1 θ′ = (α1j , βj′ , σj′2 , α0j

By the identifiabilty assumption in the statement of the Theorem, this implies G(θ) = G(θ′ ). By the definition of Θ̃∞ , we have ωj ̸= ωℓ and ωj ′ = ̸ ωℓ′ for all j ̸= ℓ and j ′ ̸= ℓ′ (with j, ℓ ∈ {1, . . . , k} and j ′ , ℓ′ ∈ {1, . . . , k ′ }). Moreover, the boundedness of A implies that all weights of G(θ) and G(θ′ ) ′ ′ are strictly positive and finite. This implies that k = k ′ , (α01 , . . . , α0k ) = (α0ρ(1) , . . . , α0ρ(k) ), and ′ ′ ′ (ω1 , . . . , ωk ) = (ωρ(1) , . . . , ωρ(k′ ) ) for some ρ ∈ Rk . Further, because ω1 ≺ · · · ≺ ωk and ω1 ≺ · · · ≺ ωk′ by the definition of Θ̃∞ , it must be the case that ρ is the identity function, implying θ = θ′ . Therefore, parameter identifiability holds on the restricted space Θ̃∞ . Application of Doob’s theorem on the restricted model. We have shown that the restricted model in Equation (30) satisfies the conditions of Doob’s theorem. Therefore, there exists Θ̃⋆ ⊆ Θ̃∞ such that P(T (θ) ∈ Θ̃⋆ ) = 1 and the restricted model is consistent at all T (θ⋆ ) ∈ Θ̃⋆ . Equivalently, ∞ for any neighborhood B ⊆ Θ̃∞ of T (θ⋆ ), we have P(T (θ) ∈ B | (Xi , Yi )ni=1 ) → 1 a.s.-gG(T (θ0 )) . Now S S define Θ⋆ := k∈N ρ∈Rk (Θ̃⋆ ∩ Θk )[ρ]. Then, by Equation (29), Π(Θ⋆ ) ≡ P(θ ∈ Θ⋆ ) = P(T (θ) ∈ Θ̃⋆ ) = 1. iid

Implications for the unrestricted model. Let θ⋆ ∈ Θ⋆ and define K ⋆ = K(θ⋆ ). Let (Xi , Yi ) ∼ g for i ∈ N, and define B := {θ ∈ Θ̃∞ : dΘ∞ (θ, T (θ⋆ )) < ε} ⊆ ΘK ⋆ for ε ∈ (0, 1). Observe that SG(θ⋆ ) ρ∈RK ⋆ B[ρ] ⊆ B̃(θ⋆ , ε), where B̃(θ⋆ , ε) :=

[

{θ ∈ Θ∞ : dΘ∞ (θ, θ⋆ [ρ]) < ε}.

ρ∈RK ⋆

So by Equation (29), Π(θ ∈ B̃(θ⋆ , ε) | (Xi , Yi )ni=1 ) ≡ P(θ ∈ B̃(θ⋆ , ε) | (Xi , Yi )ni=1 ) [ n ≥P θ∈ B[ρ] | (Xi , Yi )i=1 ρ∈RK ⋆

(32)

= P(T (θ) ∈ B | (Xi , Yi )ni=1 ) → 1 as n → ∞, ∞ a.s.-gG(θ , since gG(θ⋆ ) = gG(T (θ⋆ )) and the restricted model is consistent at all T (θ⋆ ) ∈ Θ̃⋆ . Finally, ⋆) consistency for the number of experts follows immediately from Equation (32), since ε < 1 implies ⋆ B̃(θ⋆ , ε) ⊆ ΘK and therefore ⋆

Π(K = K ⋆ | (Xi , Yi )ni=1 ) ≡ P(θ ∈ ΘK | (Xi , Yi )ni=1 ) ≥ P(θ ∈ B̃(θ⋆ , ε) | (Xi , Yi )ni=1 ) → 1 as n → ∞, ∞ a.s.-gG(θ ⋆)

A.7

Proof of Corollary 1

Note that      [ p Π G ∈ Gj : dH (gG , gG⋆ ) ≥ Mn log n/n | (Xi , Yi )ni=1    j∈N o X n p = Π G ∈ Gj : dH (gG , gG⋆ ) ≥ Mn log n/n | (Xi , Yi )ni=1 Π(K = j | (Xi , Yi )ni=1 ) j∈N

≤Π

n o p G ∈ GK ⋆ : dH (gG , gG⋆ ) ≥ Mn log n/n | (Xi , Yi )ni=1

+ Π(K ̸= K ⋆ | (Xi , Yi )ni=1 ), where the second addendum asymptotically vanishes almost surely by Theorem 6, while the first term can be dealt with exactly as in the proof of Theorem 2.

A.8

Proof of Proposition 1

Firstly, we expand equation gG (Y, X) = gG⋆ (Y, X) for almost every (Y, X) as follows: K X

exp(α0j + X ⊤ α1j ) N (Y | E(X, βj ), σj2 ) PK ⊤α ) exp(α + X 0ℓ 1ℓ ℓ=1 j=1 ⋆

K X

⋆ + X ⊤ α⋆ ) exp(α0j 1j N (Y | E(X, βj⋆ ), σj⋆2 ). PK ⋆ ⋆ + X ⊤ α⋆ ) exp(α ℓ=1 0ℓ 1ℓ j=1

(33)

As the mixture of location-scale Gaussian distributions is identifiable (Teicher, 1963), it follows that K = K ⋆ and ) ) ( ( ⋆ + X ⊤ α⋆ ) exp(α0j exp(α0j + X ⊤ α1j ) 1j : j ∈ [K] , : j ∈ [K] = PK ⋆ PK ⋆ ⊤ ⋆ ⊤ ℓ=1 exp(α0ℓ + X α1ℓ ) ℓ=1 exp(α0ℓ + X α1ℓ ) for almost every X. Without loss of generality, we may assume that ⋆ + X ⊤ α⋆ ) exp(α0j 1j , = PK PK ⋆ ⋆ ⋆ ⊤ ⊤ ℓ=1 exp(α0ℓ + X α1ℓ ) ℓ=1 exp(α0ℓ + X α1ℓ )

exp(α0j + X ⊤ α1j )

for almost every X and for all j ∈ [K]. Since the softmax function is invariant to translation, we ⋆ + t and α = α⋆ + t , for some t ∈ R and t ∈ Rd , for all j ∈ [K]. Thus, we deduce α0j = α0j 0 1j 1 0 1 1j can rewrite Equation (33) as K X

⋆

exp(α0j )F (Y |X; α1j , βj , σj2 ) =

j=1

K X

exp(α0j )F (Y |X; α1j , βj⋆ , σj⋆2 ),

(34)

j=1

for almost every (Y, X), where we define F (Y |X; α1 , β, σ 2 ) := exp(X ⊤ α1 )N (Y |E(X, β), σ 2 ). Next, we partition the index set [K] into m subsets J1 , J2 , . . . , Jm such that exp(α0j ) = exp(α0j ′ ) for any j, j ⋆ ∈ Ji and for i ∈ [m]. Meanwhile, if j and j ⋆ are not in the same subset, then we let exp(α0j ) ̸= exp(α0j ′ ). Then, we rewrite Equation (34) as m X X

exp(α0j )F (Y |X; α1j , βj , σj2 ) =

i=1 j∈Ji

m X X

exp(α0j )F (Y |X; α1j , βj⋆ , σj⋆2 ),

i=1 j∈Ji

for almost every (Y, X). The above equation implies that {(E(X, βj ), σj2 ) : j ∈ Ji } = {(E(X, βj⋆ ), σj⋆2 ) : j ∈ Ji }, for almost every X for all i ∈ [m]. Since the expert function E is identifiable, it follows that {(βj , σj2 : j ∈ Ji } = {(βj⋆ , σj⋆2 ) : j ∈ Ji } for all i ∈ [m]. As a result, we have G=

m X X

exp(α0j )δ(α1j ,βj ,σ2 ) =

m X X

i=1 j∈Ji

⋆ exp(α0j + t0 )δ(α⋆ +t1 ,β ⋆ ,σ⋆2 ) = G⋆t0 ,t1 . 1j

i=1 j∈Ji

Hence, the proof is completed.

Experiment Details

This appendix details the experimental design and variational inference methodology used to generate the model selection results presented in Section 6 of the main text. 44

B.1

Black-Box Variational Inference Implementation

To fit the candidate models, we utilize a black-box variational inference (BBVI) framework across all experiments. Variational family: We specify a fully factorized (mean-field) Gaussian variational family Q over the transformed parameter space. Specifically, we place Gaussian variational distributions over the gating intercepts α0 , gating slopes α1 , regression coefficients β, and the log-variances log σ 2 . The variational parameters µ therefore consist of the means and log-standard-deviations for each of these independent Gaussian distributions. Optimization: The ELBO is maximized using a pathwise (reparameterization) gradient estimator. During training, we draw standard normal samples ϵ ∼ N (0, 1) and apply the location-scale transformation (e.g., θ = µθ + ϵ exp(log σθ )) to obtain differentiable samples from the variational posterior. We perform full-batch optimization using the Adam optimizer. Hyperparameters and evaluation: Zero-mean Gaussian priors are placed on the gating and regression parameters (with prior variances τ1 = τ2 = τ3 = 10.0), alongside an Inverse-Gamma(2, 2) prior for the expert variances. Upon completion of the optimization routine, the final ELBO for each candidate model is estimated using 200 Monte Carlo samples from the optimized variational posterior qµ⋆ . The model attaining the highest estimated ELBO is selected as the winner for that simulation round.

B.2

Figure 1a experiment

For the first experiment, we simulate covariates uniformly from a two-dimensional square, Xi ∼ U(−1, 1)2 . The true gating mechanism assigns observations to the first expert if Xi1 > Xi2 , and to the second expert otherwise. Conditional on the expert assignment, the responses Yi are drawn from a Gaussian linear regression model. The true parameters for the first expert are set to intercept β0,1 = 2.0, slopes β1,1 = [1.0, 1.0]⊤ , and variance σ12 = 1.0. For the second expert, the parameters are β0,2 = −2.0, slopes β1,2 = [−1.0, −1.0]⊤ , and variance σ22 = 2.0. This creates a sharp discontinuity in the conditional mean surface. We evaluate candidate models K ∈ {1, 2, 3, 4} across sample sizes n ∈ {10, 25, 50, 100}, running 50 independent replications per size. For this setup, the Adam optimizer is run for 50,000 iterations. To ensure stable convergence, learning rates are tailored to each (n, K) pair, ranging from 0.0036 to 0.015.

B.3

Figure 1b experiment

For the more complex experiment, we increase the dimensionality to d = 6 and simulate covariates uniformly, Xi ∼ U(−1, 1)6 . The true model relies on K ⋆ = 4 experts. The gating mechanism is deterministic and assigns each observation based on the maximum value among the first four covariates: Zi = argmaxk∈{1,2,3,4} Xi,k . Conditional on the assignment Zi = k, the responses Yi are drawn from a Gaussian linear regression model. The true intercepts alternate in sign, set to β0,k = 2(−1)k−1 . The slope vectors β1,k are predominantly diagonal, with the k-th component set to 2(−1)k−1 , while all components are additionally perturbed by independent Gaussian noise N (0, 0.22 ). The expert variances σk2 are linearly spaced between 1.0 and 2.0. We evaluate this configuration across larger sample sizes n ∈ {100, 500, 1000, 2000}, running 50 independent replications per size. Candidate models are extended to K ∈ {1, . . . , 6}. Due to 45

the increased complexity and larger sample sizes, the Adam optimizer is run for 10,000 iterations per model, and the base learning rates are scaled linearly with the sample size to ensure stable convergence.

B.4

Table 1 experiment

We simulate n = 500 covariate vectors uniformly, Xi ∼ U(−1, 1)d , and evaluate three combinations of dimensionality and true expert counts: (d = 2, K ⋆ = 1), (d = 2, K ⋆ = 2), and (d = 4, K ⋆ = 3). The gating assignments are determined by evaluating a set of linear logits, Xi W ⊤ + b. The weight matrix W is strictly diagonal for the first K ⋆ dimensions, with Wk,k set to a predefined ‘separation‘ constant. The biases are set to bk = −0.2 × separation × k. The final expert assignment relies on a hard-max over these logits: Zi = argmaxk (Xi W ⊤ + b)k . We test two separation scales: 5.0 (low sharpness, Table 1a) and 10.0 (high sharpness, Table 1b). Conditional on the assignment Zi = k, the responses Yi are drawn from a Gaussian linear regression model. The true intercepts are linearly spaced in [−2.0, 2.0]. The slopes β1,k are given by a base diagonal value of 2(−1)k−1 , corrupted by independent Gaussian noise N (0, 0.32 ) across all dimensions. The true variances are fixed at σk2 = 0.8 for all experts. Candidate models evaluated range from K = 1 to K = 7. BBVI is run for 4,000 iterations per model (10,000 for models evaluated against the K ⋆ = 1 baseline to ensure adequate convergence). The Adam learning rate is set to an adaptive schedule, 0.1 + 0.000015n + 0.001K, with the exception of the K ⋆ = 1 baseline models, which are fixed at a rate of 0.06. Each configuration is run over 100 independent replications.

References Ascolani, F., Lijoi, A., Rebaudo, G., and Zanella, G. (2023). Clustering consistency with Dirichlet process mixtures. Biometrika, 110(2):551–558. (Cited on page 4.) Bariletto, N., Flores, B., and Walker, S. G. (2025). Posterior consistency in parametric models via a tighter notion of identifiability. arXiv preprint arXiv:2504.11360. (Cited on pages 8 and 22.) Bariletto, N. and Walker, S. G. (2025). On A Necessary Condition For Posterior Inconsistency: New Insights From A Classic Counterexample. arXiv preprint arXiv:2510.18126. (Cited on page 22.) Bariletto, N. and Walker, S. G. (2026). Scalable uncertainty quantification for black-box density-based clustering. arXiv preprint arXiv:2603.03188. (Cited on page 21.) Barron, A., Schervish, M., and Wasserman, L. (1999). The consistency of posterior distributions in nonparametric problems. The Annals of Statistics, 27:536–561. (Cited on pages 3 and 8.) Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer. (Cited on page 17.) Bishop, C. M. and Svensén, M. (2003). Bayesian hierarchical mixtures of experts. In Proceedings of the 19th Conference on Uncertainty in Artificial Intelligence (UAI-2003), pages 57–64. Morgan Kaufmann. (Cited on page 2.) Bissiri, P. G., Holmes, C. C., and Walker, S. G. (2016). A general framework for updating belief distributions. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 78(5):1103–1130. (Cited on page 21.) Blei, D. M., Kucukelbir, A., and McAuliffe, J. D. (2017). Variational inference: A review for statisticians. Journal of the American Statistical Association, 112(518):859–877. (Cited on page 17.) Dai, D., Deng, C., Zhao, C., Xu, R., Gao, H., Chen, D., Li, J., Zeng, W., Yu, X., Wu, Y., Xie, Z., Li, Y., Huang, P., Luo, F., Ruan, C., Sui, Z., and Liang, W. (2024). DeepSeekMoE: Towards ultimate expert specialization in mixture-of-experts language models. In Proceedings of the 62nd Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), pages 1280–1297. Association for Computational Linguistics. (Cited on page 10.) DeepSeek-AI et al. (2024). Deepseek-v3 technical report. arXiv preprint arXiv:2412.19437. (Cited on page 1.)

Diep, N. T., Nguyen, H., Nguyen, C., Le, M., Nguyen, D. M. H., Sonntag, D., Niepert, M., and Ho, N. (2025). On zero-initialized attention: Optimal prompt and gating factor estimation. In Proceedings of the ICML. (Cited on page 2.) Doob, J. L. (1949). Application of the theory of martingales. Le calcul des probabilités et ses applications, pages 23–27. (Cited on pages 16 and 40.) Dudley, R. M. (2002). Real Analysis and Probability. Cambridge Studies in Advanced Mathematics. Cambridge University Press, 2 edition. (Cited on page 41.) Fong, E., Holmes, C., and Walker, S. G. (2023). Martingale posterior distributions. Journal of the Royal Statistical Society Series B: Statistical Methodology, 85(5):1357–1391. (Cited on page 21.) 47

Fortini, S. and Petrone, S. (2024). Exchangeability, prediction and predictive modeling in Bayesian statistics. Statistical Science. In press. arXiv:2402.10126. (Cited on page 21.) Gao, B. and Pavel, L. (2017). On the properties of the softmax function with application in game theory and reinforcement learning. arXiv preprint arXiv:1704.00805. (Cited on page 24.) Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2013). Bayesian Data Analysis. CRC press, 3rd edition. (Cited on page 17.) Ghosal, S., Ghosh, J. K., and Ramamoorthi, R. V. (1999). Posterior Consistency of Dirichlet Mixtures in Density Estimation. The Annals of Statistics, 27(1):143–158. (Cited on pages 3, 8, and 9.) Ghosal, S., Ghosh, J. K., and Van Der Vaart, A. W. (2000). Convergence rates of posterior distributions. Annals of Statistics, pages 500–531. (Cited on pages 3, 8, and 20.) Ghosal, S. and van der Vaart, A. (2007a). Convergence rates of posterior distributions for noniid observations. The Annals of Statistics, 35(1):192 – 223. (Cited on page 8.) Ghosal, S. and Van der Vaart, A. (2017). Fundamentals of nonparametric Bayesian inference, volume 44. Cambridge University Press. (Cited on pages 8, 9, 22, and 23.) Ghosal, S. and van der Vaart, A. W. (2001). Entropies and rates of convergence for maximum likelihood and bayes estimation for mixtures of normal densities. The Annals of Statistics, 29(5):1233–1263. (Cited on page 3.) Ghosal, S. and van der Vaart, A. W. (2007b). Posterior convergence rates of dirichlet mixtures at smooth densities. The Annals of Statistics, 35(2):697–723. (Cited on page 3.) Google Gemini Team (2025). Gemini 2.5: Pushing the frontier with advanced reasoning, multimodality, long context, and next generation agentic capabilities. arxiv preprint arxiv 2507.06261. (Cited on page 10.)

Gormley, I. C. and Frühwirth-Schnatter, S. (2019). Mixture of experts models. In Handbook of mixture analysis, pages 271–307. Chapman and Hall/CRC. (Cited on page 3.) Green, P. J. (1995). Reversible jump markov chain monte carlo computation and bayesian model determination. Biometrika, 82(4):711–732. (Cited on page 18.) Guha, A., Ho, N., and Nguyen, X. (2021). On posterior contraction of parameters and interpretability in bayesian mixture modeling. Bernoulli, 27(4):2159–2188. (Cited on page 3.) Han, X., Nguyen, H., Harris, C., Ho, N., and Saria, S. (2024). Fusemoe: Mixture-of-experts transformers for fleximodal fusion. In Advances in Neural Information Processing Systems. (Cited on page 2.)

Hazimeh, H., Zhao, Z., Chowdhery, A., Sathiamoorthy, M., Chen, Y., Mazumder, R., Hong, L., and Chi, E. (2021). DSelect-k: Differentiable Selection in the Mixture of Experts with Applications to Multi-Task Learning. In Ranzato, M., Beygelzimer, A., Dauphin, Y., Liang, P. S., and Vaughan, J. W., editors, Advances in Neural Information Processing Systems, volume 34, pages 29335–29347. Curran Associates, Inc. (Cited on page 2.) 48

Ho, N. and Nguyen, X. (2016). On strong identifiability and convergence rates of parameter estimation in finite mixtures. Electronic Journal of Statistics, 10(1):271–307. (Cited on page 3.) Ho, N., Yang, C.-Y., and Jordan, M. I. (2022). Convergence rates for Gaussian mixtures of experts. Journal of Machine Learning Research, 23(323):1–81. (Cited on pages 3, 10, and 12.) Jacobs, R. A., Jordan, M. I., Nowlan, S. J., and Hinton, G. E. (1991). Adaptive mixtures of local experts. Neural Computation, 3(1):79–87. (Cited on page 1.) Jiang, A. Q., Sablayrolles, A., Roux, A., Mensch, A., Savary, B., Bamford, C., Chaplot, D. S., Casas, D. d. l., Hanna, E. B., Bressand, F., et al. (2024). Mixtral of experts. arXiv preprint arXiv:2401.04088. (Cited on pages 1 and 10.) Jordan, M. I., Ghahramani, Z., Jaakkola, T. S., and Saul, L. K. (1999). An introduction to variational methods for graphical models. Machine Learning, 37(2):183–233. (Cited on page 17.) Jordan, M. I. and Jacobs, R. A. (1994). Hierarchical mixtures of experts and the EM algorithm. Neural Computation, 6(2):181–214. (Cited on page 1.) Knoblauch, J., Jewson, J., and Damoulas, T. (2022). An optimization-centric view on bayes’ rule: Reviewing and generalizing variational inference. Journal of Machine Learning Research, 23(132):1–109. (Cited on page 21.) Kucukelbir, A., Tran, D., Ranganath, R., Gelman, A., and Blei, D. M. (2017). Automatic differentiation variational inference. Journal of Machine Learning Research, 18(14):1–45. (Cited on pages 17 and 18.)

Le, M., Nguyen, C., Nguyen, H., Tran, Q., Le, T., and Ho, N. (2025). Revisiting prefix-tuning: Statistical benefits of reparameterization among prompts. In The Thirteenth International Conference on Learning Representations. (Cited on page 2.) Le, M., The, A. N., Nguyen, H., Vu, T. T. N., Pham, H. T., Van, L. N., and Ho, N. (2024). Mixture of experts meets prompt-based continual learning. In The Thirty-eighth Annual Conference on Neural Information Processing Systems. (Cited on page 2.) Lee, H., Yun, E., Nam, G., Fong, E., and Lee, J. (2023). Martingale posterior neural processes. In The Eleventh International Conference on Learning Representations. (Cited on page 21.) Li, B., Shen, Y., Yang, J., Wang, Y., Ren, J., Che, T., Zhang, J., and Liu, Z. (2023). Sparse Mixture-of-Experts are Domain Generalizable Learners. In The Eleventh International Conference on Learning Representations. (Cited on page 2.) Lijoi, A., Prünster, I., and Walker, S. G. (2005). On consistency of nonparametric normal mixtures for Bayesian density estimation. Journal of the American Statistical Association, 100(472):1292–1296. (Cited on page 3.)

Liu, A., Feng, B., Xue, B., Wang, B., Wu, B., Lu, C., Zhao, C., Deng, C., Zhang, C., Ruan, C., et al. (2024). Deepseek-v3 technical report. arXiv preprint arXiv:2412.19437. (Cited on page 10.) 49

Ludziejewski, J., Krajewski, J., Adamczewski, K., Pióro, M., Krutul, M., Antoniak, S., Ciebiera, K., Król, K., Odrzygóźdź, T., Sankowski, P., Cygan, M., and Jaszczur, S. (2024). Scaling laws for fine-grained mixture of experts. In Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pages 33270–33288. PMLR. (Cited on page 15.)

MacEachern, S. N. (1999). Dependent nonparametric processes. In ASA proceedings of the section on Bayesian statistical science, volume 1, pages 50–55. Alexandria, Virginia. Virginia: American Statistical Association; 1999. (Cited on page 2.) Manole, T. and Ho, N. (2022). Refined convergence rates for maximum likelihood estimation under finite mixture models. In Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 14979–15006. PMLR. (Cited on pages 3 and 5.)

Masoudnia, S. and Ebrahimpour, R. (2014). Mixture of experts: a literature survey. Artificial Intelligence Review, 42(2):275–293. (Cited on page 2.) Mendes, E. F. and Jiang, W. (2012). On convergence rates of mixtures of polynomial experts. Neural Computation, 24(11):3025–3051. (Cited on page 4.) Miller, J. W. (2018). A detailed treatment of Doob’s theorem. arXiv preprint arXiv:1801.03122. (Cited on pages 16, 40, and 41.)

Miller, J. W. (2023). Consistency of mixture models with a prior on the number of components. Dependence Modeling, 11(1):20220150. (Cited on pages 4, 16, 17, 40, 41, and 42.) Miller, J. W. and Harrison, M. T. (2014). Inconsistency of pitman-yor process mixtures for the number of components. Journal of Machine Learning Research, 15(1):3333–3370. (Cited on page 4.) Ng, K., Fong, E., Frazier, D. T., Knoblauch, J., and Wei, S. (2025). Tabmgp: Martingale posterior with tabpfn. arXiv preprint arXiv:2510.25154. (Cited on page 21.) Nguyen, H., Akbarian, P., Nguyen, T., and Ho, N. (2024a). A general theory for softmax gating multinomial logistic mixture of experts. In Proceedings of the 41st International Conference on Machine Learning. (Cited on page 20.) Nguyen, H., Akbarian, P., Pham, T., Nguyen, T., Zhang, S., and Ho, N. (2025). Statistical advantages of perturbing cosine router in mixture of experts. In International Conference on Learning Representations. (Cited on page 2.) Nguyen, H., Akbarian, P., Yan, F., and Ho, N. (2024b). Statistical Perspective of Top-K Sparse Softmax Gating Mixture of Experts. In The Twelfth International Conference on Learning Representations. (Cited on page 20.) Nguyen, H., Han, X., Harris, C. W., Saria, S., and Ho, N. (2024c). On expert estimation in hierarchical mixture of experts: Beyond softmax gating functions. arxiv preprint arxiv 2410.02935. (Cited on page 10.)

Nguyen, H., Ho, N., and Rinaldo, A. (2026). Convergence rates for softmax gating mixture of experts. IEEE Transactions on Information Theory, 72(2):1276–1304. (Cited on page 4.) 50

Nguyen, H., Nguyen, T., and Ho, N. (2023). Demystifying softmax gating function in Gaussian mixture of experts. In Advances in Neural Information Processing Systems. (Cited on pages 3, 5, 10, 12, 14, and 15.)

Nguyen, X. (2013). Convergence of latent mixing measures in finite and infinite mixture models. The Annals of Statistics, 41(1):370 – 400. (Cited on page 3.) Nobile, A. (1994). Bayesian Analysis of Finite Mixture Distributions. PhD thesis, Carnegie Mellon University, Pittsburgh, PA. (Cited on pages 16 and 17.) Oldfield, J., Georgopoulos, M., Chrysos, G. G., Tzelepis, C., Panagakis, Y., Nicolaou, M. A., Deng, J., and Patras, I. (2024). Multilinear mixture of experts: Scalable expert specialization through factorization. In Advances in Neural Information Processing Systems. (Cited on page 10.) Peng, F., Jacobs, R. A., and Tanner, M. A. (1996). Bayesian inference in mixtures-of-experts and hierarchical mixtures-of-experts models with an application to speech recognition. Journal of the American Statistical Association, 91(434):953–960. (Cited on page 2.) Ranganath, R., Gerrish, S., and Blei, D. (2014). Black box variational inference. In Artificial intelligence and statistics, pages 814–822. PMLR. (Cited on pages 17 and 18.) Rasmussen, C. E. and Ghahramani, Z. (2002). Infinite mixtures of Gaussian process experts. In Dietterich, T. G., Becker, S., and Ghahramani, Z., editors, Advances in Neural Information Processing Systems 14, pages 881–888. MIT Press. (Cited on page 2.) Ren, L., Du, L., Dunson, D. B., et al. (2011). Logistic stick-breaking process. Journal of Machine Learning Research, 12(1). (Cited on page 2.) Riquelme, C., Puigcerver, J., Mustafa, B., Neumann, M., Jenatton, R., Pint, A. S., Keysers, D., and Houlsby, N. (2021). Scaling vision with sparse mixture of experts. In Advances in Neural Information Processing Systems, volume 34, pages 8583–8595. Curran Associates, Inc. (Cited on page 1.)

Robert, C. P. and Casella, G. (2004). Monte Carlo Statistical Methods. Springer. (Cited on page 17.) Rodriguez, A., Dunson, D. B., and Gelfand, A. E. (2008). The nested dirichlet process. Journal of the American Statistical Association, 103(483):1131–1154. (Cited on page 2.) Rodríguez, C. E., Mena, R. H., and Walker, S. G. (2025). Martingale posterior inference for finite mixture models and clustering. Journal of Computational and Graphical Statistics, pages 1–10. (Cited on page 21.)

Rousseau, J. and Mengersen, K. (2011). Asymptotic behaviour of the posterior distribution in overfitted mixture models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(5):689–710. (Cited on page 4.) Schwartz, L. (1965). On Bayes procedures. Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete, 4(1):10–26. (Cited on pages 3, 8, and 22.) 51

Shazeer, N., Mirhoseini, A., Maziarz, K., Davis, A., Le, Q., Hinton, G., and Dean, J. (2017). Outrageously large neural networks: The sparsely-gated mixture-of-experts layer. In In International Conference on Learning Representations. (Cited on pages 1 and 20.) Shen, X. and Wasserman, L. (2001). Rates of convergence of posterior distributions. The Annals of Statistics, 29(3):687–714. (Cited on pages 3 and 8.) Teh, Y. W., Jordan, M. I., Beal, M. J., and Blei, D. M. (2006). Hierarchical dirichlet processes. Journal of the American Statistical Association, 101(476):1566–1581. (Cited on page 2.) Teicher, H. (1963). Identifiability of finite mixtures. Annals of Statistics, 34:1265–1269. (Cited on page 44.)

Ueda, N. and Ghahramani, Z. (2002). Bayesian model search for mixture models based on optimizing variational bounds. Neural Networks, 15(10):1223–1241. (Cited on pages 2 and 18.) Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, L. u., and Polosukhin, I. (2017). Attention is all you need. In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc. (Cited on page 15.) Wainwright, M. J. (2019). High-dimensional statistics: A non-asymptotic viewpoint, volume 48. Cambridge university press. (Cited on page 23.) Wainwright, M. J. and Jordan, M. I. (2008). Graphical models, exponential families, and variational inference. Foundations and Trends® in Machine Learning, 1(1–2):1–305. (Cited on page 17.) Walker, S. G. (2004). New approaches to Bayesian consistency. The Annals of Statistics, 32(5):2028 – 2043. (Cited on pages 3 and 8.) Walker, S. G. and Hjort, N. L. (2001). On Bayesian consistency. Journal of the Royal Statistical Society, Series B, 63:811–821. (Cited on page 8.) Walker, S. G., Lijoi, A., and Prünster, I. (2005). Data tracking and the understanding of bayesian consistency. Biometrika, 92(4):765–778. (Cited on page 22.) Walker, S. G., Lijoi, A., and Prünster, I. (2007). On rates of convergence for posterior distributions in infinite-dimensional models. The Annals of Statistics, 35(2):738–746. (Cited on pages 3 and 8.) Wang, C. and Yang, Y. (2024). Estimating the number of components in finite mixture models via variational approximation. arXiv preprint arXiv:2404.16746. (Cited on pages 4 and 18.) Waterhouse, S., MacKay, D., and Robinson, T. (1995). Bayesian methods for mixtures of experts. In Advances in Neural Information Processing Systems, volume 8. (Cited on pages 2 and 18.) Wong, W. H. and Shen, X. (1995). Probability inequalities for likelihood ratios and convergence rates of sieve MLEs. The Annals of Statistics, pages 339–362. (Cited on page 26.) Wu, L. and Williamson, S. A. (2024). Posterior uncertainty quantification in neural networks using data augmentation. In International conference on artificial intelligence and statistics, pages 3376–3384. PMLR. (Cited on page 21.) 52

Zhang, Y. and Yang, Y. (2024). Bayesian model selection via mean-field variational approximation. Journal of the Royal Statistical Society Series B: Statistical Methodology, 86(3):742–770. (Cited on pages 4 and 18.)

Record · ID 124070 · SHA-256 8277561ee70e5362
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.