Estimating the Number of Components in Finite Mixture Models via Variational Approximation
Chenyang Wang Yun Yang





#1 1








Estimating the Number of Components in Finite Mixture Models via Variational Approximation
    Chenyang Wangcw80@illinois.edu  1, Yun Yangyy84@illinois.edu  1

  1Department of Statistics, University of Illinois Urbana-Champaign
    
=========================================================================================================================================




1



This work introduces a new method for selecting the number of components in finite mixture models (FMMs) using variational Bayes, inspired by the large-sample properties of the Evidence Lower Bound (ELBO) derived from mean-field (MF) variational approximation. Specifically, we establish matching upper and lower bounds for the ELBO without assuming conjugate priors, suggesting the consistency of model selection for FMMs based on maximizing the ELBO. As a by-product of our proof, we demonstrate that the MF approximation inherits the stable behavior (benefited from model singularity) of the posterior distribution, which tends to eliminate the extra components under model misspecification where the number of mixture components is over-specified. This stable behavior also leads to the n^-1/2 convergence rate for parameter estimation, up to a logarithmic factor, under this model overspecification. Empirical experiments are conducted to validate our theoretical findings and compare with other state-of-the-art methods for selecting the number of components in FMMs.




Keywords:    Evidence lower bound; Mean-field approximation; Mixture models; Model selection; Singular models.


1.15 


§ INTRODUCTION

One of the fundamental problems in statistics and machine learning is model selection, which attempts to select the most parsimonious model to fit the observed data from a list of candidate models. To address this problem, various model selection criteria have been proposed in the literature. Among them, two well-known criteria are the Akaike Information Criterion (AIC), proposed by <cit.>, and the Bayesian Information Criterion (BIC), introduced by <cit.>. Under suitable conditions, AIC is known to lead to optimal prediction performance but tends to overestimate the model size; in comparison, BIC, which is derived as an asymptotic approximation to the log-marginal likelihood (a.k.a. model evidence) under a Bayesian framework, is able to consistently identify the true model (the model with the smallest size that contains the true data-generating process) as the sample size grows to infinity under the frequentist perspective, a property known as model selection consistency. The validity of BIC in approximating the model evidence requires the considered candidate model to be non-singular, meaning that its Fisher information matrix needs to be invertible at the true model parameter so that the maximum likelihood estimator (MLE) enjoys the asymptotic normality property. However, in many practical applications, especially when the considered model is over-specified, the Fisher information matrix is non-invertible, a situation referred to as singular models. Here, an over-specified model means that this model contains the true model, but the number of parameters in the model exceeds that of the true model; an accompanying notion is an under-specified model, where the model does not contain the true model. 

In this paper, we consider the problem of model selection for a common and representative type of singular models, the finite mixture model (FMM), whose distribution family consists of all finite mixtures of probability distributions from a shared parametric family, with a fixed number of mixture components. 
FMMs have been extensively studied in the statistical literature <cit.> and are widely applied in various fields, including agriculture, astronomy, biology, economics, engineering, and psychology, among others.
In an FMM, the model size is determined by the number of components, which is typically unknown in practice. Thus, the model selection problem essentially becomes the task of identifying the minimal number of components required in the FMM to adequately fit the data. Over-specified FMMs are singular due to the non-identifiability of model parameters; in other words, different parameters can correspond to the same probability distribution.
Over-specifying the number of components in an FMM also has the drawback of deteriorating the statistical efficiency of parameter estimation. For example, in a two-component Gaussian mixture model that is over-specified with the true model being a single Gaussian, the convergence rate of the MLE for estimating the location parameter decreases from the parametric rate n^-1/2 to n^-1/4, as shown by <cit.>; more generally, this rate of convergence further slows down with the inclusion of additional components <cit.>.

Extensions of BIC-type model selection criteria to singular models date back to a series of seminal works <cit.> by Watanabe. In these works, tools from algebraic geometry, particularly those based on the resolution of singularities, are employed to determine the leading terms in the asymptotic expansion of the model evidence. The leading terms now depend on the real log canonical threshold (RLCT) and its multiplicity of the model; see Section <ref> for a brief review or the recent book <cit.> for further information. However, unlike the usual BIC 
that only depends on the considered model, the RLCT also requires knowledge of the true parameter value and therefore cannot be used in practice for model selection. Recently, <cit.> proposed a singular Bayesian information criterion (sBIC) for model selection in the context of singular models. This method involves solving a system of equations, the unique solution of which is shown to provide valid asymptotic approximations to the model evidences of candidate models. However, this approach requires knowledge of the real log canonical threshold (RLCT) and its multiplicity for all candidate models, the determination of which poses a highly challenging problem even in field of algebraic geometry. For example, for multi-variate Gaussian mixture models, their complete characterizations remain unresolved.

In this paper, we propose and analyze a new model selection approach for singular models based on variational Bayes (VB). VB has become a popular alternative to traditional sampling-based algorithms like MCMC for approximate Bayesian computation, owing to its computational scalability and good generalization performance in many applications.
Our proposed approach draws inspiration from the recent work of <cit.>, which shows that for non-singular models, the evidence lower bound (ELBO) derived from the mean-field (MF) approximation in variational Bayes can recover the leading terms in model evidence. Notably, the model selected by the ELBO tends to converge asymptotically with the one chosen by the BIC as the sample size increases. Furthermore, the calculation of the ELBO is a straightforward by-product of the Coordinate Ascent Variational Inference (CAVI) algorithm <cit.> for implementing the MF approximation. In particular, this approach eliminates the need for a case-by-case analysis to determine the effective number of parameters required in computing the BIC for non-singular models. In this work, we also consider the use of the ELBO from MF approximation for model selection in FMMs. Theoretically, we find that while the ELBO can no longer recover the leading term involving the RLCT in the model evidence, it still contains a term that penalizes over-specified models, leading to consistent model selection. 




 §.§ Related literature

In this subsection, we briefly review more recent works closely related to ours and defer a more extensive literature review and discussions to Appendix <ref>.
For Bayesian finite mixture models, <cit.> studied the asymptotic behavior of posterior distributions in over-specified Bayesian FMMs, and proved that under certain conditions on the prior, the posterior tends to empty the extra components. <cit.> and <cit.> derived the contraction rate of posterior distributions on the parameters of interest in Bayesian FMMs under model misspecification; more reviews of literature on mixture model estimation can also be found therein. For variational Bayes (VB), recent works <cit.> have provided general conditions under which VB leads to estimation consistency of the model parameter or data distribution, and derived the corresponding convergence rate. In particular, the theoretical framework developed in <cit.> and <cit.> covers Bayesian latent variable models, including FMMs. Moreover, the theory in <cit.> implies an n^-1/2 convergence rate, up to a logarithmic factor, for over-specified FMMs under the Hellinger metric, a result that will also play a key role in our proof.
In the context of applying VB to Bayesian FMMs, <cit.> derived upper and lower bounds for the ELBO (also called stochastic complexity therein) from the MF approximation for Bayesian Gaussian mixture models; also see Section <ref> for further details. However, the gap between these two bounds does not ensure model selection consistency based on the ELBO. <cit.> extended this analysis to general FMMs, but the issue of the gap remains unresolved.
More recently, <cit.> has shown that the ELBO from MF approximation recovers the leading term involving RLCT for a singular model in its normal-crossing form. However, while Hironaka’s theorem in algebraic geometry guarantees the existence of a reparametrization that can locally transform the posterior distribution of any singular model into its normal-crossing form, explicitly finding such a transformation is typically challenging. 





 §.§ Our contributions

In this study, we investigate the large sample properties of the ELBO derived from the MF variational approximation for FMMs under model misspecification. Specifically, we establish matching upper and lower bounds for the ELBO for a general class of FMMs whose constituting mixture components are exponential families. Our refined result will imply the consistency of model selection based on the ELBO. Moreover, unlike earlier analysis relying on conjugate priors to simplify the analysis, we only require the smoothness of the priors and employ Laplace approximation to analyze the conditional posterior distribution of parameters given latent class indicator variables. 
In our proof, we first establish a lower bound for the ELBO by constructing specific feasible solutions to the VB optimization problem, similar to those used by <cit.>. Our main technical contribution comes from our upper bound analysis. We utilize the variational risk bound under the Hellinger metric given by <cit.> to identify certain important behaviors of the mixing weights, yielding a substantially improved upper bound that matches the lower bound up to high-order terms. Interestingly, our findings indicate that, unlike singular models in normal-crossing forms <cit.>, directly applying MF to Bayesian FMMs does not recover the leading term in the model evidence, although it still results in consistent model selection. 
In addition, as a by-product of our analysis, we prove a corollary demonstrating that the MF variational approximation to the posterior distribution in Bayesian FMMs inherits the same interesting and stable behavior of tending to empty the extra components under model over-specification <cit.>.  This stable behavior also leads to the n^-1/2 convergence rate for parameter estimation, up to a logarithmic factor, for those statistically identifiable components that do not empty out (see Section <ref>). It is also worth mentioning that our results do not require strong identifiability, as assumed by <cit.> and <cit.>, and are applicable to weakly identifiable FMMs <cit.>.
Through experiments conducted on a location Gaussian mixture model and a real dataset, we empirically validate our theoretical findings and demonstrate that our model selection method is computationally efficient and tends to have higher power than existing state-of-the-art methods — it can correctly select the true number of components in FMMs using fewer samples.



The rest of the paper is organized as follows. Section <ref> outlines the notation, provides background knowledge, formulates the problems we aim to address in this work and introduce our model selection method in finite mixture models. Section <ref> lists our main assumptions and presents our main theoretical results about model selection via ELBO maximization and parameter estimation via mean-field variational approximation. Section <ref> showcases simulated experiments of the Gaussian mixture model under various settings to complement our theoretical findings and reports results from a real data analysis. Finally, Section <ref> concludes the paper with a discussion. Proofs and additional literature reviews are deferred to the appendices in the supplementary material.





§ BACKGROUND AND NEW METHOD
 

In this section, we begin with an overview of the notation, followed by a review of necessary background, and then introduce the problem formulation.




 §.§ Notation

Throughout the paper, random variables are denoted by capital letters, while their realizations are represented by corresponding lowercase letters. For two probability measures P and Q, let p and q be their Radon-Nikodym derivatives with respect to a σ-finite measure μ, usually chosen as the Lebesgue measure or the counting measure. We denote the total variation and the (squared) Hellinger distance between p and q respectively as d_ TV(p,q) = 1/2∫ |p-q|  dμ and h^2(p,q) = ∫(√(p) - √(q))^2  dμ. Moreover, their Kullback-Leibler (KL) divergence is denoted as D_ KL(pq) = ∫ plog(p/q)  dμ.

We use 𝒩(μ,Σ) to denote a (multivariate) normal distribution with mean μ∈ℝ^d and variance-covariance matrix Σ∈ℝ^d× d.

Additionally, for a vector a∈ℝ^d, we use |a| to denote its Euclidean ℓ_2-norm. Throughout this paper, we use C to denote a generic positive constant whose value may vary from line to line. We also write a≳ b to mean a≥ Cb, a≲ b to mean a≤ Cb and a≍ b to mean a≳ b and a≲ b.





 §.§ Singular models

In statistical learning, many classical techniques and conclusions are applicable primarily to regular models, such as the Laplace approximation, and the convergence rate and asymptotic normality of the maximum likelihood estimations (MLE). A statistical model ℳ={p(· | θ) | θ∈Θ} with parameter space Θ is called regular (or non-singular) if its parameter θ is identifiable and p(· | θ) has a well-defined and positive definite Fisher information matrix I(θ) = E_θ[-∇^2 log p(X | θ)] at each parameter θ∈Θ.  In contrast, we say a model is singular when at least one of the conditions for regularity is not met. Many common statistical models are singular, such as finite mixture models, hidden Markov models, factor models, and neural networks. Mixture models, characterized by their ability to represent complex data distributions through convex combinations of simpler distributions from a shared parametric family, are a representative type of singular model. In this article, we focus exclusively on finite mixture models.

Concretely, let us start from a simple parametric family G = {g(· ; η) | η∈Ω} parameterized by parameter η∈Ω, whose parameter space Ω is a subset of ℝ^d, such as a location-scale family. For an integer K, a mixture model with K components in G is model ℳ_K ={p(· | θ) | θ=(w_k,η_k)_k=1^K,  η_k∈Ω, w_k≥ 0, k∈[K] ∑_k=1^Kw_k=1} with

    p(x | θ) = ∑_k=1^K w_k g(x; η_k),  θ =(w_k,η_k)_k=1^K ∈ R^dK+K,

where w =(w_1,w_2,…,w_K)∈ [0,1]^K is called the mixing weight (vector) parameter, and, by slightly abusing the notation, we use η=(η_1,η_2,…,η_K)∈Ω^K to denote the parameters associated with the K mixture components.
We also denote the corresponding (discrete) mixing distribution as G(θ) = ∑_k=1^K w_k δ_η_k,
where δ_η_k denotes the Dirac mass at η_k.
Generally, mixture models may lose identifiability, since different parameters may lead to the same data distribution p(· | θ). However, G(θ) is always identifiable, that is, there is a one-to-one correspondence between G(θ) and p(· | θ). This serves as one of the key assumptions in our theory (see Section <ref>). For example, the following two (Gaussian) mixture models correspond to the same standard normal distribution: p(x | θ_1) = 1·𝒩(0,1) + 0·𝒩(1,1) and p(x | θ_2) = 1/2·𝒩(0,1) + 1/2·𝒩(0,1). Fortunately, they still have the same mixing distribution G(θ_1) = G(θ_2)=δ_0.

The loss of identifiability may lead to a non-invertible Fisher information matrix, resulting in decreased statistical efficiency — even when an identifiable submodel is considered, the best point estimator may suffer from a slow convergence rate as the sample size n increases. For example, <cit.> demonstrated such a phenomenon when the mixing distribution is restricted to the form G(θ) = 2/3δ_-h + 1/3δ_2h for h∈ℝ, and g satisfies certain regular conditions to make the reduced parameter h identifiable. However, even for this identifiable submodel, the MLE for h still converges at the slow rate of n^-1/4 when the true parameter h^∗ = 0. Therefore, it is crucial to correctly specify the number K of mixture components before conducting statistical inference for finite mixture models, which motivates the problem of model selection for K in this paper.





 §.§ Model selection criteria for singular models

The Bayesian Information Criterion (BIC) is widely used for model selection for regular parametric models. Let X^n = (X_1, ..., X_n) be n i.i.d. observations from the true model p^*(x), which may or may not belong to the considered model M.
Under this setup, the BIC for a model M is derived as an asymptotic approximation to the so-called model evidence, or the log-marginal likelihood in a Bayesian framework, defined as

    F(ℳ) = log∫_Θ_ℳπ(θ | ℳ) p(X^n | θ,ℳ)   dθ =  log p(X^n | θ_ M,ℳ)- d/2log n_ M + 𝒪_P(1),

where Θ_ℳ∈ R^d denotes its parameter space, π(θ |  M) the prior density, p(X^n | θ,ℳ) the likelihood function, and 
θ_ M=_θ∈Θ_ Mp(X^n | θ,ℳ) is the MLE of θ under model M. The BIC in the preceding display consists of two terms: a model goodness-of-fit term and an additional term of -d/2log n resulting from the integration over ℝ^d, which penalizes larger models that may lead to overfitting. It is a classical result <cit.> that selecting a model based on maximizing BIC enjoys the model selection consistency property.
Unfortunately, for singular models with non-invertible Fisher information matrices, the asymptotic expansion (<ref>) for the model evidence no longer holds, as the effective dimension of the parameter space is no longer equal to the number d of parameters.

In a seminal work,  <cit.> employed algebraic geometry tools to show that for a singular model M, the asymptotic behavior of the log-marginal likelihood can be characterized through

    F(ℳ) = log p^*(X^n) - λ_ Mlog n + (m_ M-1)loglog n + 𝒪_P(1),

where λ_M is called the real log canonical threshold (RLCT) of model M that determines the effective number of parameters, and m_ M is its multiplicity. Leveraging this general result, <cit.> proposed a model selection method by maximizing a singular Bayesian information criterion (sBIC), which is computed by solving a system of equations whose unique solution approximates the leading terms in F(ℳ). However, this approach requires knowledge of the real log canonical threshold (RLCT) and its multiplicity for all candidate models, the determination of which is a highly nontrivial problem even in the field of algebraic geometry within pure mathematics.




 §.§ Selecting number of components in mixture models via variational approximation

To handle the analytically intractable integral in F(ℳ), variational inference attempts to approximate F(ℳ) defined in (<ref>) from below by utilizing Jensen's inequality (to -log(t)),

    log∫_Θ_ℳπ(θ | ℳ) p(X^n | θ,ℳ)/q_θ(θ) q_θ(θ)   dθ≥∫_Θ_ℳ q_θ(θ) logπ(θ | ℳ) p(X^n | θ,ℳ)/q_θ(θ) dθ: = L(q_θ),

where q_θ(θ) is a distribution from a tractable distribution family Γ for approximating the posterior distribution p(θ | X^n, M). We will focus on the commonly used mean-field family Γ_ MF in this work, where q_θ factorizes over the (blocks of) components of θ. The evidence lower bound (ELBO) is then defined as the best approximation to F(ℳ) by maximizing the lower bound L(q_θ) over q_θ∈Γ_ MF. In the rest of the paper, when we consider a fixed model ℳ, we may suppress ℳ in our notation for simplicity.

Now let us specialize the above general framework to finite mixture models (<ref>).
As a common practice, we will augment the mixture model into a hierarchical model via introducing latent (class indicator) variables to facilitate efficient computation.  In fact, a key property utilized in our theoretical analysis is that the conditional distribution of observation given latent variable is a regular parametric model. Let S^n = (S_1, ..., S_n) denoting the collection of all latent variables, with each S_i being discrete and taking values from {1,..,K}, indicating which mixture component that X_i comes from. Then we have the hierarchical form of the model as [X_i | S_i=k]iid∼ g(x ; η_k) and S_i are i.i.d. with P(S_i=k) = w_k, for i∈[n] and k∈[K]. Moreover, we can express
the joint distribution of X^n given θ as

    p(X^n | θ) = ∑_s^n∈[K]^n p(X^n | s^n,θ) p(s^n | θ) = ∏_i=1^n∑_s_i∈[K] p(X_i | s_i,η) p(s_i | w).

Let Z^n = (θ,S^n) denote the collection of all hidden variables, then the augmented posterior distribution p(Z^n | X^n) satisfies p(θ,s^n | X^n) p(X^n,s^n | θ) π(θ). Now the variational inference will use q_Z^n∈Γ to approximate the joint posterior distribution p(Z^n | X^n) (by replacing θ with Z^n therein). 
To reduce the approximation error, we consider the block MF approximation where Γ_ MF={q_Z^n = q_θ⊗ q_S^n}. Let q_Z^n = q_θ⊗q_S^n := _q_Z^n∈Γ_ MF L(q_Z^n) denote the best variational density that approximates p(Z^n | X^n), and the corresponding L(q_Z^n) is the ELBO for the mixture model with K components. Then the best variational approximation to its posterior p(Z^n | X^n) is defined as the minimizer of the following KL divergence,

    q_Z^n := inf_q_Z^n∈Γ D_ KL(q_Z^n(·)  p(· |X^n)).

Note that for any q_Z^n, we have the following relationship 

    D_ KL(q_Z^n(·)  p(· |X^n)) = log∫ p(X^n|θ) π(θ)  dθ - ℒ(q_Z^n)≥ 0,

where the first term is the evidence and the second term 

    ℒ(q_Z^n) = ∫ q_Z^n(z^n)logp(X^n,s^n|θ) π(θ)/q_Z^n(z^n)d z^n.

is the corresponding ELBO function, since the KL divergence is non-negative. Therefore, maximizing the ELBO over q_Z^n is equivalent to minimizing the KL divergence. In addition, the ELBO valued at its maximizer q_Z^n(z^n) is a potential surrogate to the model evidence in singular model selection, as long as this approximation retains some characteristics of the evidence. As aforementioned, a commonly used variational family Γ is the (block) mean-field (MF) approximation, where q_Z^n(z^n) is factorized as

    q_Z^n(z^n) = q_θ(θ)⊗ q_S^n(s^n).

We will focus on the MF, and use q_Z^n to denote the maximizer of the ELBO function ℒ(·), which,  according to the MF definition, can be further decomposed into q_Z^n = q_θ⊗q_S^n. We will refer to q_Z^n as the variational posterior (distribution). If we use ℒ_K = ℒ_K(q_Z^n) to denote the (maximal) ELBO value under the FMM with K components, then our selected number of components K via ELBO maximization can be defined via

    K = _Kℒ_K.



Regarding the block MF approximation, the following theorem gives a necessary condition for q_Z^n to be the maximizer of the ELBO, as stated in the Theorem 2.1 by <cit.>. This theorem provides a system of distributional equations determining q_θ(θ) and q_S^n(s^n), which plays an important role and will be repeatedly used in our proofs.

    Under the block MF approximation, the variational posterior q_Z^n = q_θ⊗q_S^n satisfies
    
    q_θ(θ) = 1/C_rπ(θ)exp{∑_s^n∈ [K]^nq_S^n(s^n)log p(X^n, s^n | θ)}
         q_S^n(s^n) = 1/C_Qexp{∫q_θ(θ)log p(X^n, s^n | θ)  d θ},

    where C_r and C_Q are the normalization constants.



Under the above settings, <cit.> derived upper and lower bounds for ℒ(q_Z^n) under a symmetric Dirichlet prior Diri(ϕ_0,…,ϕ_0), ϕ_0 ≥ 0, on the mixing weights w and a conjugate prior for the parameters η associated with mixture components from an exponential family. Specifically, if we denote the variational posterior mean (estimator) of θ as θ:= ∫θ q_θ(θ)  dθ, then <cit.> proved that (after some notation adaptions)

    log p^*(X^n) - λlog n + 𝒪_P(1) ≤ℒ(q_Z^n)≤log p^*(X^n) -λlog n + logp(X^n | θ)/p^*(X^n) + 𝒪_P(1)

as sample size n tends to infinity, and the leading coefficients λ and λ are given by

    λ= 
            (K-K^*)ϕ_0 + dK^*+K^*-1/2,    ϕ_0≤d+1/2,
    dK+K-1/2,    ϕ_0>d+1/2;
                λ= 
            (K-1)ϕ_0 + d/2,    ϕ_0≤d+1/2,
    dK+K-1/2,    ϕ_0>d+1/2,
 
where K^* denotes the true number of components underlying the data-generating process.
When ϕ_0<(d+1)/2, there exists a gap between λ and λ. It should be noted that (d+1)/2 also serves as the threshold for ϕ_0, below which the desirable stable behavior of the posterior distribution and its variational approximation (see Corollary <ref>) to empty the extra components can be maintained. This stable behavior is also shown to lead to improved parameter estimation efficiency (see Theorem <ref>). In contrast, when ϕ_0>(d+1)/2 and K is over-specified, then two or more components can be very close with non-negligible weights each, leading to unstable and less accurate parameter estimation (for example, see our numerical results in Section <ref>). In practice, a choice of ϕ_0 less than one is also typically preferred in the literature. For example, a default choice of ϕ_0 in topic models <cit.> is K^-1, which encourages a sparse mixing weight w for easier interpretation. Conversely, a larger ϕ_0 tends to promote a denser mixing weight; in particular, Diri(ϕ_0,…,ϕ_0) tends to concentrate around the uniform weights (K^-1,…,K^-1) as ϕ_0 approaches infinity. Therefore, it is an important theoretical question whether this gap can be eliminated for small ϕ_0, so that model selection based on maximizing the ELBO is provably consistent. In this work, we close this gap by proving matching lower and upper bounds for the ELBO under any positive values of ϕ_0. The discussion following our Theorem <ref>, Section <ref>, and our numerical results (e.g., Section <ref>) both suggest that a small ϕ_0 value is indeed preferred for accurate model selection and parameter estimation.




§ THEORETICAL RESULTS

In this section, we present our main theoretical results. Here, we assume that the dataset X^n = (X_1, ..., X_n) consists of n observations that are i.i.d. samples from a fixed true data-generating process p^*(x) = ∑_k=1^K^* w_k^*g(x;η_k^*),
where K^* denotes the true component number, and θ^* =  (w_k^*, η_k^*)_k=1^K^*∈ℝ^dK^* + K^* denotes the (unique) true parameters under the mixture model (<ref>) with K^* components, where all η_k's are distinct and min_k∈[K^*] w_k^∗>0. Let G^* denote the corresponding mixing measure for the truth. Throughout this section, we consider a finite mixture model with K≥ K^* components. In addition, we follow the literature by assuming the parametric family G = {g(· ; η) | η∈Ω} for the mixture components to be an exponential family in the canonical form, that is, g(x;η) = exp{η^T T(x)-T_0(x)-A(η)}, where η∈ R^d is the natural parameter, T(x) is the sufficient statistic and A(η) is the log-partition function. The Fisher information matrix for this exponential family is I(η) = ∇^2 A(η), which implies the convexity of A(η) over its domain ℝ^d. Consequently, the log-likelihood function log g(x;η) is concave with respect to η for all x. One important assumption of our theory (see Section <ref> below) is the relationship between total variation and Wasserstein-type distances in FMMs, which transforms the closeness of mixture models in density functions into closeness in their mixing measure or model parameters. Such a relationship can be verified for most commonly used FMMs; see Section <ref> for concrete examples. An r-th order Wasserstein distance (or simply r-Wasserstein distance) between two mixing measures G(θ) = ∑_k=1^K w_k δ_η_k and G(θ') = ∑_k=1^K' w'_k δ_η'_k is defined as:

    W_r(G(θ),G(θ')) = (inf_q_ij∑_i,jq_ij|η_i - η_j'|^r)^1/r,

where the infimum is taken over all nonnegative probability masses {q_ij:  i∈[K],  j∈[K']} satisfying two marginal constraints ∑_i=1^K q_ij = w'_j for j∈[K'] and ∑_j=1^K' q_ij = w_i for i∈[K].




 §.§ Assumptions
 
Recall that we have the symmetric Dirichlet distribution Diri(ϕ_0,…,ϕ_0), ϕ_0 > 0 as our prior for the mixing weights w, whose density takes the form of
π( w) = Γ(Kϕ_0)/Γ(ϕ_0)^K∏^K_k=1 w_k^ϕ_0 - 1, for all w=(w_1,…,w_K)∈[0,1]^K such that ∑_k=1^Kw_k=1, where Γ(·) is the Gamma function.
We also consider independent priors on η, that is, π(η) = ∏_k=1^Kπ(η_k), η_k∈Ω, where for each k, the density π(η_k) is twice continuously differentiable on Ω. In addition, we make the following assumptions.

Assumption A: Recall that p^* is the true data-generating probability with mixing measure G^* and has K^* components.

    
  * (Parameter space compactness) The parameter space Ω of mixture components is a compact subset of ℝ^d, i.e., there exists a constant L<∞ such that sup_η,η'∈Ω|η-η'|≤ L.
    
  * (Regularity of mixture component family) There exist two positive constants a and b such that y^T I(η) y∈[a,b] over all |y| = 1 and all η∈Ω, where I(η) = ∇^2A(η) is the Fisher information matrix for the mixture components.
    
  * (Identifiability of mixing measure) The total variation distance between any FMM p(x|θ) with K(≥ K^∗) components and the true model p^∗ satisfies d_ TV(p(· |θ),p^*(·)) = 0 if and only if its mixing measure G(θ) satisfies G(θ) = G^*. Moreover, there exist an positive integer r and positive constants ε and c_0 such that 
    d_ TV(p(· |θ),p^*(·))≥ c_0 W_r^r(G(θ),G^*) for any p(· |θ) with mixing measure G(θ) (with K components) satisfying W_r(G(θ),G^*)<ε.
 


Assumption A1 is commonly made in the literature for technical simplicity when analyzing singular models. Assumption A2 is a mild condition with Ω being compact and is satisfied by all non-singular (i.e., not curved) exponential families. Assumption A3 requires certain partial identifiability for the true model p^* with K^* components. It essentially states that for any mixing measure G(θ) sufficiently close to G^*, the total variation of the density functions is lower bounded by some Wasserstein-type distance between the mixing measures. For example, if the family g(x;η), η∈Ω is identifiable in the second-order (see Definition 3.2 in <cit.>), in other words, g(x;η_k) and their derivatives with respect to η_k up to the second order are linearly independent, then Assumption A3 is satisfied with r = 2. For families that are not identifiable in the second-order, for instance, the location-scale Gaussian distributions, verification of Assumption A3 can be done in a case-by-case manner, where r may also depend on K^* and K. See Section <ref> for more examples that satisfy these assumptions.

It is also worth noting that, with Assumptions A1 and A3, the (localized) inequality in the second part of Assumption A3 can be made global. Specifically, there exists a constant c_1 such that d_ TV(p(· |θ), p^*(·)) ≥ c_1 W_r^r(G(θ), G^*) for all possible θ. The proof is straightforward. Let 𝒜:={θ: W_r(G(θ), G^*) ≥ε, η∈Ω^K}. Since W_r(G(θ), G^*) is continuous in θ, 𝒜 is a compact set and d_ TV(p(· |θ), p^*(·)) > 0 on 𝒜 by the first part of Assumption A3. For the considered exponential family, the total variation is also a continuous function with respect to θ. Therefore, it attains a strictly positive infimum on 𝒜, which we denote as τ := inf_θ∈𝒜 d_ TV(p(· |θ), p^*(·)) > 0. Meanwhile, from Assumption A1, we have W_r^r(G(θ), G^*) ≤ L^r. Therefore, on 𝒜, we also have d_ TV(p(· |θ), p^*(·)) ≥ (τ/L^r) W_r^r(G(θ), G^*). Letting c_1 = min{c_0, τ/L^r}, it then follows from above and the second part of Assumption A3 that for any θ, we have:

    d_ TV(p(· |θ), p^*(·))≥ c_1 W_r^r(G(θ),G^*),

which provides a global inequality relating the two distance metrics over any compact set. The above argument allows us to avoid assuming a reverse inequality
d_ TV(p(· |θ),p^*(·)) ≤ c_0 W_r^α(G(θ),G^*) for some α > 0 and for any p(· |θ) with mixing measure G(θ) in order to extend the local inequality into a global one, as was done in Corollary 3.1 of <cit.>.




 §.§ Large sample properties of ELBO
 
Now we are ready to present our first main result as follows, which provides two-sided bounds to the ELBO ℒ(q_Z^n) with matching leading terms. Recall that we use θ:= ∫θ q_θ(θ)  dθ to denote the variational point estimator for θ.


   Suppose Assumptions A1–A3 hold and K≥ K^*. Then there exist some positive constants C_1, C_2 and c independent of n such that with probability at least 1 - n^-c, we have
    
    log p^*(X^n)-λlog n + C_1≤ℒ(q_Z^n) ≤log p^*(X^n)-λlog n +  logp(X^n | θ)/p^*(X^n) +  C_2,

    where C_1=KlogΓ(ϕ_0) - Γ(Kϕ_0) - C and C_2 = KlogΓ(ϕ_0) - Γ(Kϕ_0) + C for some constant C independent of (n, ϕ_0), and
    λ is given by
    
    λ = 
            (K-K^*)ϕ_0 + dK^*+K^*-1/2,    ϕ_0≤d+1/2,
    dK+K-1/2,    ϕ_0>d+1/2.


We note that KlogΓ(ϕ_0) - Γ(Kϕ_0) is an increasing function with respect to both K and ϕ_0. In particular, a large ϕ_0 may cause the constants C_1 and C_2 in the ELBO to rapidly grow with the model size K; as a result, a larger sample size is needed to correctly select the true model by ELBO maximization. See our numerical results in Section <ref> for an empirical comparison, where at a sample size of n=50, using ELBO with ϕ_0=1 selects the correct model, while using ELBO with ϕ_0=6 (d = 6) leads to the selection of an over-specified model. Therefore, our theoretical result suggests that in practice a small ϕ_0 value is preferred for accurate model selection with limited data. In Section <ref> below, we will also demonstrate the theoretical benefits of using a small ϕ_0 < (d+1)/2 for parameter estimation.

There is still a gap Δ_n=log p(X^n | θ) -log p^*(X^n) in Theorem <ref>, which is related to the likelihood ratio test statistic LRT = log p(X^n | θ ) -log p^*(X^n) through Δ_n≤ LRT, with θ=_θ∈Θlog p(X^n | θ) denoting the MLE of θ. 
By utilizing the singular learning theory <cit.>, we show in Appendix <ref> that this LRT weakly converges to a limiting random variable defined through the supreme of a Gaussian process as n→∞, by using the invariance property of MLE and a resolution of singularity technique developed in <cit.> which transform the log-likelihood function of the singular model into its canonical normal crossing form. As a result, we obtain the following corollary, ensuring the gap in our ELBO upper and lower bounds to be at most 𝒪_P(1).


    Under Assumptions A1–A3, the asymptotic expansion of ELBO for K≥ K^* has the following form:
    
    ℒ(q_Z^n) = log p^*(X^n)-λlog n + KlogΓ(ϕ_0) - Γ(Kϕ_0) + 𝒪_P(1),  n→∞.



As a direct consequence, we obtain the following corollary about model selection consistency based on maximizing ELBO. Here, recall that K^* denotes the (true) number of components in the data-generating process, and K=_Kℒ_K denotes the estimated number of components via ELBO maximization, with L_K being the ELBO value under the FMM with K components. 

    Suppose Assumptions A1–A3 hold, then the estimated number of components K  in (<ref>) via ELBO maximization satisfies 
    
    ℙ(K= K^*| X^n) p→ 1,  n→∞.

    This means that model selection via ELBO maximization in FMMs is (asymptotically) consistent.




 §.§ Parameter estimation in finite mixture models via variational approximation

From the Theorem <ref>, we can also see that the large sample behavior of ELBO exhibits two different patterns, depending on the value of the prior hyperparameter ϕ_0 for the mixing weight w. For large ϕ_0 that exceeds (d+1)/2, the effective number of parameters in the variational approximation q_Z^n is (dK+K-1), which matches the degrees of freedom of parameter θ = ( w,η) (we lose one degree of freedom due to the constraint ∑_k=1^K w_k=1). Therefore, this case can be interpreted as the non-singular regime since ELBO and the regular BIC asymptotically match each other. The more interesting case occurs when ϕ_0 < (d+1)/2, which corresponds to the singular regime, as the effective number of parameters λ is strictly smaller than the degrees of freedom (dK+K-1) in the parameterization θ. In fact, the following corollary as a by-product of our proof suggests that under this singular regime, the variational approximation q_ w to the posterior distribution on the mixing weight will have an interesting stable behavior of emptying out the (K-K^*) redundant components when K is overspecified. 


    Suppose Assumptions A1–A3 holds and there exists some constant C_3 such that Δ_n≤ C_3loglog n with probability 1-o(1). Then it holds with probability at least 1 - n^-c - o(1) that: 1. when ϕ_0<(d+1)/2 and for ρ_1 = 2(C_3 + 1)/(d + 1-2ϕ_0),
    
    inf_σ∈𝒮_K∑_k = K^*+1^K w_σ(k)< (K-K^*)[(log n)^ρ_1 + ϕ_0]/n,

    where 𝒮_K is the set of all permutations over [K]; 2. when ϕ_0>(d+1)/2 and for ρ_2 = 2(C_3 + 1)/(2ϕ_0 - d -1),
    
    inf_k∈[K]w_k ≥1/(log n)^ρ_2+ϕ_0/n.


Our condition Δ_n≤ C_3loglog n requires a finite sample analysis of the likelihood ratio test statistic under singularity, and can be verified in many concrete examples <cit.>; also see Appendix <ref> for a brief review.
This corollary indicates that for small ϕ_0, the redundant mixture components are emptied out at the rate of n^-1 (up to log n factors), a rate faster than the usual parametric rate of n^-1/2. This phenomenon is again a distinctive feature of the singular regime, where some components of the posterior distribution may exhibit a faster rate of convergence, depending on the behavior of the prior distribution near the boundary of the parameter space <cit.>. In comparison, under the regular regime where ϕ_0 > (d+1)/2, the mixing weight point estimates tend to spread out over all K components and are not stable, and the “super-efficiency” property disappears; this finding is also supported by our numerical results in Section <ref>. A similar result regarding this robust behavior for the posterior distribution p(w | X^n) of mixing weights w is also proved in <cit.>; however, their obtained convergence rate under small ϕ_0 is n^-1/2 up to logarithm terms and requires a stronger identifiability condition on the model (i.e., identifiability in the second-order). Therefore, Corollary <ref> indicates that variational approximation also inherits the stability and (possibly) super-efficiency property of the posterior distribution in finite mixture models. 

When ϕ_0 < (d+1)/2, due to the rapid emptying out speed of redundant components, only K^* statistically identifiable components survive. It is then natural to ask the question of whether the convergence rate of the parameters associated with the true components (i.e., those K^* components whose weights do not empty out) can reach the parametric rate n^-1/2 (possibly modulo logarithmic terms). The answer is affirmative. Before presenting the result, we need to first introduce a precise notion of identifiability on the parameter estimation.


    The family {g(x;η),η∈Ω} is said to be weakly identifiable, or identifiable in the first order, if g(x;η) is differentiable in η and for any finite K different {η_k}_k=1^K, for some α_k∈ℝ and β_k∈ℝ^d, the following equation holding for almost every x,
    
    ∑_k=1^K[α_k g(x;η_k) + β_k^T∇_η g(x;η_k) ]= 0,

    implies α_k = 0 and β_k = 0 for all 1≤ k≤ K.


This condition of weak identifiability, or identifiability in the first order, is the least stringent form of identifiability condition for parameter estimation in FMMs in the literature (e.g., see <cit.>), and is satisfied by most of the commonly used distributions. The term “weak” is used in contrast to the so-called strong identifiability condition (e.g., see <cit.>), which is also called identifiability in higher orders (with orders greater than one). For example, a one-dimensional family {g(x;η), η∈Ω} is said to be identifiable in m-th order (or m-strongly identifiable) if for any finite K different {η_k}_k=1^K, for some α_k∈ℝ and β_k^(j)∈ℝ^d, j∈[m], the following equation holding for almost all x,

    ∑_k=1^K[α_k g(x;η_k) + ∑_j=1^mβ_k^(j) g^(j)(x;η_k) ]= 0,

implies α_k = 0 and β_k^(j) = 0 for all 1≤ k≤ K and j=1,…,m. Here, g^(j)(x;η_k) denotes the j-th order derivative of g(x;η) with respect to η. Note that the strong identifiability always implies the weak identifiability.
A strong identifiability condition (with order two) is commonly adopted in the literature to characterize the convergence rate of location or scale mixture models in the literature <cit.>. However, for weakly identifiable models, the convergence rate of a general FMM must be examined on a case-by-case basis (e.g., see <cit.>) by proving an inequality similar to our Assumption A3, which relates the total variation distance to a Wasserstein-type distance. Moreover, the derived convergence rate also depends on the order r in the inequality (see Assumption A3) and is typically slower than the square root-n rate for parameters. In contrast, our result below indicates that, thanks to the stability behavior of mixing weights in the singular regime (Corollary <ref> under ϕ_0 < (d+1)/2), weak identifiability, coupled with Assumption A3 under any order r, leads to a parametric rate for the parameters associated with those identifiable mixture components. From Section <ref> below, we will see that some simple FMMs, such as the location-scale Gaussian family, are only weakly identifiable and not strongly identifiable. Therefore, our results imply that the variational posterior estimator of θ can achieve the n^-1/2 convergence rate for parameter estimation, up to logarithmic factors, for weakly identifiable FMMs such as the location-scale Gaussian mixture model. This is a result that cannot be inferred from the existing literature (even for the MLE).



    Suppose the same assumptions of Corollary <ref> hold and the family {g(x;η),η∈Ω} is weakly identifiable. If ϕ_0<(d+1)/2, then there exists a large constant M such that it holds with probability at least 1-n^-c-o(1) that
    
    inf_σ∈𝒮_Ksup_1≤ k≤ K^*|η_σ(k) - η_k^*| ≤Mlog n/√(n),

     where 𝒮_K is the set of all permutations over [K].

Theorem <ref> can be interpreted as a consequence of Corollary <ref> where the mixing weights of redundant components decrease to zero at a sufficiently fast rate. Consequently, despite the overspecification, we can effectively treat the working models as well-specified models with the correct number of components. In such well-specified cases, <cit.> demonstrated in their Theorem 3.1 that d_ TV(p(· |θ), p^*(·)) ≥ c_0 W_1(G(θ), G^*) when W_1 is small, provided that the family is identifiable in the first order. This leads to a nearly parametric convergence rate for estimating the mixture component parameters, due to the nearly parametric convergence rate under the total variation metric. This heuristic argument illustrates how the stability behavior of mixing weights in the singular regime enhances parameter estimation efficiency. The numerical results in Section <ref> also demonstrate the advantages of using a small ϕ_0 to induce this stability behavior.



 §.§ Applications to concrete examples

In this subsection, we verify Assumption A and the weak identifiability condition for several representative finite mixture models, and determine their effective dimensions, λ, arising in the ELBO from the MF approximation.

[Location Gaussian mixture model]
    In this example, we consider the location Gaussian mixture model where the density function of each component is from the univariate location family of Gaussian distributions {𝒩(μ,σ^2);μ∈Ω}, where σ^2 is assumed to be known and the canonical parameter in this exponential family is simply η : = μ. Moreover, the Fisher information matrix is 1/σ^2 times the identity matrix. For simplicity, if we assume σ^2 = 1 without loss of generality, then
    
    g(x;μ) = 1/√(2π)exp{-x^2/2 + μ x -μ^2/2}.

    According to Theorem 3 in <cit.>, g(x;μ) is identifiable in the second order, which also implies the weak identifiability; in fact, it is also identifiable in any finite order according to Theorem 2.4 in <cit.>. Therefore, Assumption A3 is satisfied with r = 2 by Theorem 3.2 (a) from <cit.>, and our Theorem <ref> applies with λ being 
    
    λ = 
            (K-K^*)ϕ_0 + 2K^*-1/2,    ϕ_0≤1/2,
    2K-1/2,    ϕ_0>1/2.




[Scale exponential mixture model]
    In this example, we take g(x;η) = η e^-η x = exp{η x + logη} as the (scale) exponential distribution family with rate parameter η (or scale parameter η^-1), whose fisher information is 1/η^2. According to Theorem 2.4 in <cit.>, g(x;μ) is also identifiable in any finite order. Therefore, Assumption A3 is satisfied with r = 2 and the weak identifiability also holds. As a result, Theorem <ref> also applies to the scale exponential mixture model with the same λ as Example <ref>.


[Location-scale Gaussian mixture model]
    We consider the classical form of the univariate location-scale Gaussian distribution family where
    
    g(x;μ,σ^2) = exp{-x^2/2σ^2 + xμ/σ^2 - μ^2/2σ^2-1/2log (2πσ^2)}.

    This family is no longer strongly identifiable due to the algebraic relationship between the two partial derivatives (up to the second order) 
    
    ∂^2 g/∂μ^2(x;μ,σ^2) = 2∂ g/∂σ^2(x;μ,σ^2).

    For this parametric family, Proposition 2.2 from <cit.> shows that if G(θ) has at most K components, then for any G(θ) such that W_r(G(θ),G^*) is sufficiently small, we have
    
    d_ TV(p(· | θ),p^*(·))≳ W_r^r(G(θ),G^*),

    for some constant r≥ 2 depending only on K-K^*. This verifies Assumption A3 with r=r. Hence, we can apply Theorem <ref> to conclude that the effective dimension λ in the ELBO is (d=2)
    
    λ = 
            (K-K^*)ϕ_0 + 3K^*-1/2,    ϕ_0≤3/2,
    3K-1/2,    ϕ_0>3/2.

  Moreover, although the location-scale Gaussian distribution family is not strongly identifiable, it is weakly identifiable, as demonstrated in <cit.>. Consequently, all results from Section <ref> regarding parameter estimation via mean-field approximation remain applicable.



[Multinomial mixture model] In the final example, we consider a mixture of discrete distributions, specifically, the multinomial mixture model that is widely used in topic modeling. In particular, we take g(x;η) = Mx_1, …, x_d∏_j = 1^d η_j^x_j, x_j∈[M], ∑_j=1^d x_i=M. This is the probability mass function of the multinomial distribution Mult(M,(η_1,…,η_d)) with a known M and parameters (η_1,…,η_d), where ∑_j=1^dη_j =1 and has an effective dimension of d-1. A sufficient condition for this multinomial mixture model with K components to be strongly identifiable is 3K -1≤ M according to Corollary 1 in <cit.>. Therefore, Assumption A3 holds with r=2 and our Theorem <ref> applies with 
    
    λ = 
            (K-K^*)ϕ_0 + dK^*-1/2,    ϕ_0≤d/2,
    dK-1/2,    ϕ_0>d/2,

   where one degree of freedom is lost due to the linear constraint ∑_i  = 1^dη_i = 1. The weak identifiability is again implied by the strong identifiability.





 §.§ Technical highlights

In this subsection, we highlight some key technical steps and results in our proofs that could be interesting in their own right.
The main goal of introducing Assumption A3 is to obtain the following Lemma <ref>, which implies that when the total variation distance between a K component mixture model and the true model with K^* components is very close to zero, it is always possible to merge its K components into K^* groups I_1,…,I_K^*⊂ [K], so that the aggregated mixing weights {∑_j∈ I_kw_j:  k∈[K^*]} are close to the true mixing weights w^*={w^*_1,…,w^*_K^*}.



   Under Assumptions A2 and A3, there exists a positive constant c so that for any θ, we have
    
    d_ TV(p(· |θ),p^*(·))≥ cinf_I_1,...I_K^*∑_k=1^K^*|∑_j∈ I_kw_j - w_k^*|,

    where the infimum is taken over all index sets I_1, I_2, …, I_K^* that are disjoint subsets of [K]. 

Note that {I_1, I_2, …, I_K^*} is not necessarily a partition of [K]; in other words, the union of these subsets is not necessarily [K]. Lemma <ref> establishes a relationship between the total variation distance between two FMMs and their component weights. Since model selection is concerned only with the number of components and not with the parameter estimation of η_k's, this inequality plays a key role and is also sufficient for proving our main conclusions regarding the model selection consistency. Compared to a typical theoretical analysis on the estimation consistency, an advantage of adopting this lemma is that we no longer need to require the mixture model to be strongly identifiable <cit.>; instead, we only need to assume some weaker identifiable conditions as in Assumption A3. A proof of this lemma is provided in Appendix <ref>.

Finally, we remark that in order to improve the ELBO upper bound analysis from <cit.>, one of the key steps in our proof is to utilize the following result about the consistency of variational Bayes, which is adapted from Theorem 3.2 in <cit.> and is also applicable to singular models. In particular, this result, combined with Lemma <ref> above, implies the consistency of the estimation of the mixing weights w (after properly merging similar redundant components) in the mean-field variational approximation.


    Under Assumption A1, there exist constants (c', C_4) such that it holds with probability at least 1-n^-c',
    
    ∫_Θ h^2(p(x | θ), p^*(x))  q_θ(θ)  dθ≤ C_4 ( log p^*(X^n) - ℒ(q_Z^n)/n + log n/n).


This lemma is a direct consequence of the result from <cit.> since their testing condition is automatically satisfied for the Hellinger distance with a compact parameter space; so we omit the proof. By using the inequality d_ TV(p(x | θ), p^*(x)) ≤ h(p(x | θ), p^*(x)), we can then obtain that, under Assumption A3, at least K^* of the mixing weights are bounded away from 0 with high probability under the variational posterior distribution q_w. This property leads to an improved upper bound of ℒ(q_Z^n) (which matches the lower bound) compared to the one obtained in <cit.> (reviewed at the end of Section <ref>).




§ NUMERICAL STUDIES

In this section, we further demonstrate our theoretical findings about ELBO through simulation experiments, and illustrate and compare the two predicted regimes, namely, the singular regime where ϕ_0<(d+1)/2 and the regular regime where ϕ_0>(d+1)/2, in terms of their model selection and parameter estimation performances. Additionally, we compare our proposed method of selecting the number of components in finite mixture models via ELBO maximization with other state-of-the-art selection methods and apply them to a real dataset.





 §.§ Impact of ϕ_0 on ELBO in finite mixture models

To study the effectiveness of using the ELBO ℒ(q_Z^n) analyzed in Theorem <ref> for selecting the true number of components, and to examine how different choices of ϕ_0 may impact the selection results, we conduct a numerical experiment using a location Gaussian mixture model under different settings. The ELBO is computed using the Coordinate Ascent Variational Inference (CAVI) algorithm <cit.>. The simulated data X^n=(X_1,…,X_n) is generated from a d-dimensional Gaussian mixture model with K^* = 2 components and d = 6. We use a common variance-covariance matrix for the mixture components, which is the identity matrix I and treated as known. The true mixing weights are w_1^* = w_2^* = 1/2 and the two true location parameters are η_1^* = -2/√(d)·1 and η_2^* = 2/√(d)·1, where 1 denotes the d-dimensional all one vector. We employ a symmetric Dirichlet prior Diri(ϕ_0,…,ϕ_0) for the mixing weight w and an independent N(0, I) prior for the location parameters {η_k}_k=1^K. We consider K_max=5 candidate Gaussian mixture models with component number K∈{1,2,…,K_max}. For our implementations, we found that the CAVI algorithm can be sensitive to the initialization of mixing weights w, especially in the regular regime when ϕ_0>(d+1)/2.
Therefore, in every simulation, the following initialization scheme is applied to enhance the computational stability: we randomly and evenly assign all data points X^n into K_ int groups at the beginning to initialize the latent class indicators S^n for each K_ int∈{1,2,…,K}, and then iteratively update q_S^n(s^n) and q_θ(θ) until the ELBO value converges. For each K, we select the largest ELBO value obtained among different random initializations (i.e., under different K_ int) as our final output of ℒ_K, the ELBO value in the K components mixture model. 


In the first setting, we fix the sample size at n∈{10^4, 10^5} and let K and ϕ_0 change. We choose ϕ_0 ∈{1,2,3.5,4.5,6}, where (d+1)/2=3.5 corresponds to the critical threshold distinguishing the singular and regular regimes, as predicted in Theorem <ref>. In particular, our theory predicted that for K>K^*, the ELBO value L_K satisfies

    L_K -  L_K^* =𝒪_P(1) - min{ϕ_0, 3.5}  (K-K^*) log n.

Consequently, when plotting ( L_K -  L_K^*)/log n against K, the resulting curve will approximate a straight line with a slope of -min{ϕ_0, 3.5}. The numerical results shown in Figure <ref> confirm this theoretical prediction. The two graphs on the left show the ELBO values overall considered K, while the two graphs on the right zoom in on the over-specified range of K≥ K^*. As we can see, for all cases, ELBO is maximized at the true model with K=K^*=2. Additionally, under a relatively small sample size, as in Figure <ref>(b), the slope under ϕ_0=6 has a clear deviation from the theoretical value of 3.5 based on equation (<ref>); this situation improves under the relatively larger sample size, as in Figure <ref>(d). This observation can be explained by our discussion after Theorem <ref>, which suggests that using a large ϕ_0 may not be a good choice for small sample sizes due to higher approximation errors (C_1,C_2).




To formally study the impact of ϕ_0 on model selection, we consider a second setting where we fix ϕ_0 and vary n∈{10,10· 5,10· 5^2,10· 5^3,10· 5^4} and K∈{2,3,4,5}. Figure <ref> reports the results, where we plot [ L_K - log p^*(X^n)] against log n with a fixed ϕ_0∈{1,3.5,6}. As can be seen from Figure <ref>, for ϕ_0=1, the ELBO values L_K is already peaked at K^* at a sample size as small as n=10. However, for ϕ_0=6, a much larger sample size (around e^4.5 = 90) is required for the ELBO to peak at the true model with K^*=2. This empirical result again suggests that a smaller value of ϕ_0 is preferred for accurate model selection with limited data.





We also use the setting of a fixed K and varied n to numerically verify the theoretical dependence of the coefficients of log n in the ELBO on ϕ_0. According to Theorem <ref>, if we plot ℒ(q_Z^n) - log p(X^n|θ) against log n, we expect to see a linear relationship with a slope of -λ, which equals to (K-K^*)ϕ_0 + (dK^* + K^* -1)/2 for ϕ_0<3.5 and (dK + K -1)/2 for ϕ_0>3.5. The corresponding results are summarized in Figure <ref>. As we can see, the slopes of the lines for larger ϕ_0 ={3.5, 4.5, 6} exceeding the critical threshold (d+1)/2=3.5 remains essentially unchanged across different K. For ϕ_0<3.5, the lines becomes less steep compared to ϕ = 3.5 as K increases. Those empirical observations again align well with the prediction from our theory. 

Finally, to examine the stability behavior of the variational approximation q_w which tends to empty out the redundant mixture components under over-specified K>K^*, we conduct a last simulation study to report the variational point estimators of the mixing weights (expectations under q_w), denoted as w, under K=5 and K^∗=2.
In Table <ref>, we report the (sorted) estimated weights w_k's for each setting. We repeated every setting ten times and recorded the means and standard deviations. The behavior of w_k is consistent with our predictions from Corollary <ref>: for small ϕ_0<3.5, the mixing weights tend to only concentrate on the first two components; while for large ϕ> 3.5, the mixing weights tend to spread out over all K=5 mixture components. Moreover, for large ϕ≥ 3.5, the estimated weights tend to be fairly sensitive to the initialization so the standard deviations are higher, which again suggests that a small ϕ_0 should be used in order to enhance the parameter estimation robustness and accuracy. In addition, in Table <ref> we report the estimated Wasserstein distances between the estimated mixing distribution G(θ) and the true mixing distribution G^* to compare the convergence rate of component parameters. 
From the table, we can observe that the estimation errors in terms of W_1(G(θ), G^*) (and standard deviations) under both n=1000 and n=10000 are substantially smaller under a small ϕ_0<3.5 (singular regime) compared to those under a large ϕ_0> 3.5 (regular regime).










 §.§ Comparison with other model selection methods

Recently, a method for model selection of finite mixture models (FMMs) named Group-Sort-Fuse (GSF) has been proposed by <cit.>, which is shown to enjoy higher consistency and efficiency than existing approaches such as the Merge-Truncate-Merge (MTM) algorithm provided by <cit.>. Specifically, GSF is a penalized likelihood method to simultaneously estimate the number of components and the mixing measure G(θ). Three commonly used penalties, the Smoothly Clipped Absolute Deviation (SCAD), the Minimax Concave Penalty (MCP), and the Adaptive Lasso (ALasso), which meet the required theoretical conditions, are advocated in their paper. In <cit.>, they already compared their method to the MTM algorithm, which has lower accuracy and incurs high computational costs since it is based on sampling from the real posterior distribution. Therefore, in this simulation study, we will only study and compare the three methods, ELBO (under ϕ=1 and ϕ=5), BIC, and GSF with the three different penalties (denoted as GSF_S, GSF_M and GSF_A) for selecting the number of components in a location Gaussian mixture model. For the model setting, we adopt a similar data-generating model as <cit.>. Specifically, we consider the six mixture models as summarized in Table <ref> (two within each group) in the Appendix <ref> with true number of components K^* = 2, 3, 5 and mixture component parameter dimensions d = 2, 4, 6. The component parameters η_k^* (location parameters) in Models 1, 3, and 5 are the same as those from <cit.>, but we change the mixing weights w^*=(w_1^*,…,w_5^*) to be unevenly distributed. Models 2, 4, and 6 are less separated (i.e., with closer location parameters {η_k^∗}_k=1^5), with the same component numbers and dimensions as in Models 1, 3, and 5, respectively. For each model, we assume the identity covariance matrix as known. We randomly generate the data X^n 100 times with varying sample sizes n = 200, 400, 600, 800 for each model, then apply each method and record the model selection accuracy (percentages of selecting the correct model).


The findings are illustrated in Figure <ref>. A notable observation is that ELBO demonstrates superior performance compared to all other methods in scenarios where the separation between mixture components is relatively small (Models 2, 4, and 6). Particularly in Model 6, ELBO achieves an accuracy exceeding 0.875 with a sample size of 800, while alternative methods yield accuracies of 0. The tendency of other methods to select smaller models is attributed to the small (cluster) separation, often resulting in the selection of only one component. Moreover, in Models 3, and 5, where component weights are unevenly distributed, GSF exhibits suboptimal performance compared to the findings presented in <cit.>, as GSF imposes a large penalty when some component weights approach zero, leading to model underestimation. Additionally, for ELBO, similar to our earlier results, larger values of ϕ_0 (regular regime) yield lower accuracy compared to smaller ϕ_0 (singular regime), especially when the sample size is small, as demonstrated in Model 2, 4, and 6. However, overall, even ELBO with a large ϕ_0 still outperforms other methods. In summary, our simulation results for Gaussian mixture models indicate that ELBO generally outperforms GSF under the following conditions: (1) large number of components and parameter dimensions; (2) uneven distribution of mixing weights; and (3) similar mixture components or small cluster separations.






 §.§ ELBO and evidence

In this subsection, we empirically compare and calculate the difference between the ELBO and the true model evidence. <cit.> derived the real log canonical threshold (RLCT) for the univariate location Gaussian distribution. Note that the RLCT determines the (asymptotic) leading coefficient in front of log n for the true model evidence. The RLCT values in (<ref>) are given as

    λ_ M = K^* - 1 + j + j^2 + 2(K-K^*+1)/4(j+1)    j = max{i: i+i^2≤2(K-K^*+1)}.

We visualize these metrics by plotting the ELBO, which includes the empirical evidence (denoted as Evidence_ E) derived from MCMC sampling and the theoretical evidence (denoted as Evidence_ T) by plugging-in the theoretical RLCT value; the data-generating model is p^*(x) = 1/2𝒩(-2,1) + 1/2𝒩(2,1) with a sample size of n = 10000. The empirical evidence is obtained via the Monte Carlo method using samples generated from the posterior distribution using a Metropolis-Hastings algorithm, with a Monte Carlo sample size of N = 20000.   From the results in Figure <ref>, we observe that the ELBO under a relatively small ϕ=1/4 (singular regime) tends to be closer to both evidence values, while that under a relatively large ϕ=1 (regular regime) resembles the BIC, which is consistent with our theoretical predictions.
Moreover, despite the slight discrepancies between ELBO values and evidence values, they remain close and have the same overall trend as model size grows, ensuring that the ELBO inherits the consistency of model selection from using the true evidence.





 §.§ Real Data Analysis



In this subsection, we apply our method to the old Faithful geyser eruption data, which contains 272 two-dimensional data points. By examining the scatter plot (Section 4.2.2 in <cit.>), it appears that two components are sufficient to represent the data.
We preprocess the data to equalize the variances of the two dimensions by dividing the waiting time by 15 and setting σ^2 = 0.25. We apply each method using 5%, 7.5%, 10%, 15%, and 25% of the entire dataset to determine the minimal sample size required for each method to correctly select the true number of components. For each fraction, we repeat the subsampling 100 times and report the averaged number of components selected by different methods. As shown in Figure <ref>, when ϕ_0 = 1, the (averaged) number of components selected by the ELBO is already very close to 2 with just 14 data points (5% fraction) and remains unchanged when the fraction exceeds 7.5%. However, with ϕ_0 = 4, a preference for a larger model (6 is the upper limit we set) is observed when the sample size is small. All three penalties of the GSF algorithm and the BIC show a trend toward correctly selecting the true number of components as the sample size increases. However, when the sample size is relatively small, BIC tends to incorrectly select an overly large model compared to the ELBO. This phenomenon could again be explained by the stability behavior of the ELBO (with a small ϕ), which makes the variational posterior distribution under an overspecified model θ more accurate than the maximum likelihood estimator θ (see Theorem <ref> and Table <ref>), resulting in better power to distinguish true signals (i.e., separation between mixture components) from estimation errors due to a small size. A deterioration of the ELBO model selection performance observed in Figure <ref> using a larger ϕ_0=4 also empirically supports this explanation.
A similar overestimation of the model size occurs with the three GSF methods, which could also be due to the same reason.




§ DISCUSSION

In this paper, we advocate the use of ELBO from mean-field variational Bayes for model selection in finite mixture models. We analyze the large-sample properties of ELBO by proving matching lower and upper bounds that improve upon existing results by utilizing the consistency property of variational Bayes. As a direct consequence, we have established an asymptotic expansion for the ELBO and proved the consistency of model selection through ELBO maximization. As a by-product of our proof, we show that the stable behavior of the posterior distribution, which benefits from model singularity, can be inherited by the mean-field variational approximation to the posterior. We also derived a parametric convergence rate of log n/√(n) for component parameters.

As one future direction, we would like to extend the current development to other singular models such as hidden Markov models, neural networks, and factor models. Another future direction is to develop a finite sample analysis of the ELBO by providing a more refined characterization of the likelihood ratio test statistic under singularity, as we briefly commented after Corollary <ref>.



apalike









    Supplementary Materials for “Model Selection for Finite Mixture Models via Variational Approximation”





§ MORE LITERATURE REVIEW

In this appendix, we offer additional details on some of the literature reviewed in the main paper, along with information on several more related works.

Variational Bayes.
Variational inference was introduced by <cit.> for probability density approximation with intractable integrals and point estimation for parameter determination. In variational inference, the posterior distribution is approximated by the closest member relative to the Kullback-Leibler (KL) divergence in a specified family. Among the various approximating schemes, mean-field approximation, where the variational family adopts a factorized form, emerges as the most prevalent type of variational inference. It is characterized by its conceptual simplicity, ease of implementation, and particular suitability for addressing problems that involve a high number of latent variables. <cit.> relates the Bayes risk to the variational solution for a general distance metric. In addition, for non-singular models, which include finite mixture models with a known component number, <cit.> proved that the convergence rate of point estimators based on the variational posterior is n^-1/2, and they also provide the asymptotic normality of the optimal mean-field approximation centered at the maximum likelihood estimation.

Asymptotic properties of likelihood ratio tests in mixture models.
The use of maximum likelihood estimator (MLE) for fitting mixture models has received a lot of attention from statisticians since the 1960s. <cit.> introduced a general iterative approach, the expectation–maximization (EM) algorithm, for computing the MLE in latent variable models. The convergence properties of the MLE for the mixture problem were then theoretically established. Meanwhile, the testing of homogeneity, in other words, determining whether the mixture model has only one component or more, became an important research question. For non-singular models, when under the null hypothesis the parameter is constrained in a subspace of the whole parameter space, the limit distribution of the likelihood ratio testing statistic (LRTS) Λ_n is a chi-square distribution. Unlike in non-singular model testing, <cit.> observed that under homogeneity, the statistic diverges to infinity if the whole parameter space is unrestricted. <cit.> developed the asymptotic theory for the distribution of the LRTS under this setting. They showed that in the limit, 2logΛ_n is distributed as the square of the supremum of a centered Gaussian process W_S:S∈𝒮 admitting continuous sample paths. Also, <cit.> investigated the null behavior of the LRTS for this model. At the beginning of the 21st century, more mixture models were considered, and the divergence rate of the LRTS with an unconstrained parameter space was verified. For example, normal mixture models in <cit.> and gamma mixture models in <cit.> were shown to diverge at the rate 𝒪(loglog n). On the other hand, mixture models with restricted parameter spaces were also explored. While the limit distribution is still related to a Gaussian process, this time the index set 𝒮 is compact, so the LRTS is bounded in probability. The general results for models with restricted parameter spaces under loss of identifiability are provided by <cit.>, and the application to finite mixture models can also be found in Section 4 of the same paper. Although their method may test the true model against the alternative, it requires prior knowledge about the true component number to set the null hypothesis, and the exact form of the Gaussian process still requires a case-by-case analysis.



§ TABLE <REF>







§ PROOFS OF MAIN RESULTS

In this appendix, we collect all proofs of the theoretical results in the paper.



 §.§ Proof of Lemma <ref>

Consider the two mixing distributions G = ∑_k=1^K w_k δ_η_k and G^* = ∑_k=1^K^∗ w^*_k δ_η^*_k.
Let δ:= min_j≠ k|η^*_j - η^*_k| and I_k' := {j:|η_j - η_k^*|<δ/2} for k∈[K^*], then I_k''s are disjoint by definition. Recall that W_r^r(G,G^*) = inf_q_ij∑_i=1^K∑_j=1^K^* q_ij|η_i-η_j^*|^r, where the infimum ranges over all q_ij satisfying ∑_i∈[K]q_ij =  w^*_j and ∑_j∈[K^*]q_ij =  w_i.

For each k∈[K^*] and feasible q_ij, we will show 

    W_r^r(G,G^*) ≥|∑_i∈ I_k'w_i - w_k^*| δ^r.

To see this, we consider two cases. In the first case when ∑_i∈ I_k'w_i > w_k^*, we have
    
    W_r^r(G,G^*)    ≥inf_q_ij∑_i∈ I_k'∑_j=1^K q_ij|η_i-η_j^*|^r≥inf_q_ij(∑_i∈ I_k'q_ik|η_i-η_k^*|^r + ∑_i∈ I_k'∑_j≠ kq_ij|η_i - η_j^*|^r)
        ≥ 0 + inf_q_ij∑_i∈ I_k'∑_j≠ kq_ij(δ/2)^r.

    Since ∑_i∈ I_k'∑_j∈[K^*]q_ij = ∑_i∈ I_k'w_i and ∑_i∈ I_k'q_ik≤∑_i∈[K] q_ik = w_k^* for all feasible q_ij, we obtain ∑_i∈ I_k'∑_j≠ k q_ij≥∑_i∈ I_k'w_i - w_k^*. This inequality combined with (<ref>) leads to bound (<ref>).
   In the second case when ∑_i∈ I_k'w_i ≤ w_k^*, we have
    
    W_r^r(G,G^*) ≥inf_q_ij∑_i=1^K q_ik|η_i - η_k^*|^r = inf_q_ij( ∑_i∈ I_k'q_ik|η_i - η^*_k|^r + ∑_i∉ I_k'q_ik|η_i - η^*_k|^r).

    Since ∑_i∈[K]q_ik = w^*_k and ∑_i∈ I_k'q_ik≤∑_i∈ I_k'∑_j∈[K^*]q_ij = ∑_i∈ I_k'w_i, we can obtain that ∑_i∉ I_k'q_ik≥ w_k^* - ∑_i∈ I_k'w_i, which further implies bound (<ref>) due to the preceding display and the definition of I_k.

   Finally, since K^* is finite, we obtain from inequality (<ref>) that
    
    W_r^r(G,G^*)    ≥ (K^*)^-1∑_k=1^K^*|∑_i∈ I_k'w_j - w_k^*| δ^r 
       ≳∑_k=1^K^*|∑_i∈ I_k'w_j - w_k^*| δ^r ≥inf_I_1,...,I_K^*∑_k=1^K^*|∑_i∈ I_kw_j - w_k^*| δ^r.

   The claimed bound then follows by combing above with inequality (<ref>).





 §.§ Proof of Theorem <ref>


The proof is divided into three steps. Firstly, we decompose the ELBO into two parts and obtain some analytically manageable approximations to both parts. Secondly, we constructively prove the claimed lower bound by using two instances of q_Z^n, since for any q_Z^n, L(q_Z^n) provides a lower bound to the ELBO. Finally, we prove the claimed upper bound based on the key Lemma <ref> and our Assumption A2.

Step 1: For the complete-data {X^n,S^n}, we denote the posterior probability of X_i coming from the kth component as p_ik = q_S^n(S_i = k), and let n_k = ∑_i=1^n p_ik. Then according to formula (<ref>), the variational posterior distribution q_θ(θ) = q_ w( w) ⊗q_η(η) can be written as
    
    q_ w( w) = Γ(n + Kϕ_0)/∏^K_k=1Γ(n_k+ϕ_0)∏^K_k=1 w_k^(n_k+ϕ_0) - 1,

    and
    
    q_η(η)∝∏^K_k=1π(η_k)exp[∑_i=1^n p_ik(η_k^TT(x_i) - T_0(x_i) - A(η_k))].

   From equations (<ref>) and (<ref>), we observe that q_w(w) and q_η(η) are parameterized by (in other words, only depend on) {p_ik: i∈[n],k∈[K]}. Additionally, q_S^n(s^n) is also paramaterized by {p_ik} due to formula (<ref>). Therefore, we only need to analyze those q_θ(θ) = q_w(w)⊗ q_η(η) taking the form of
    
    q_ w( w) = Γ(n + Kϕ_0)/∏^K_k=1Γ( n_k+ϕ_0)∏^K_k=1 w_k^( n_k+ϕ_0) - 1,

    and
    
    q_η(η)∝∏^K_k=1π(η_k)exp[∑_i=1^n p_ik(η_k^TT(X_i) - T_0(X_i) - A(η_k))].

    where n_k = ∑_i=1^n p_ik; and given a q_θ(θ), the optimal q_S^n(s^n) takes the following form with a normalization constant C_Q (also parameterized by {p_ik}),
    
    q_S^n(s^n) = 1/C_Qexp{∫q_θ(θ)log p(X^n, s^n | θ)d θ}.

    In the rest of the proof, we only consider those q_ w( w), q_η(η) and q_S^n(s^n) having the forms as (<ref>), (<ref>) and (<ref>). Based on these characterizations, we know that maximizing ℒ(q_Z^n) over q_Z^n(z^n) = q_θ(θ)⊗ q_S^n(S^n) is equivalent to maximizing ℒ(q_Z^n) over all possible combinations of {p_ik: i∈[n],k∈[K]} subject to 0≤ p_ik≤1 and ∑_k=1^K p_ik = 1 for all i. 
    
    We now characterize the ELBO function. According to the definition of ℒ(q_Z^n) in (<ref>), we can decompose the ELBO into two parts as
    
    ℒ(q_Z^n)     = ∫ q_Z^n(z^n)logp(X^n,s^n | θ)π(θ)/q_Z^n(z^n)d z^n
        = ∫ q_θ(θ)logπ(θ)/q_θ(θ)dθ + ∑_s^n q_S^n(s^n)∫ q_θ(θ) logp(X^n,s^n | θ)/q_S^n(s^n) dθ,

    where the first part is just -D_ KL(q_θ(θ)π(θ)) and the second term is log C_Q. Then we can reformulate ℒ(q_Z^n) as
    
    ℒ(q_Z^n) = max_{p_ik}{-D_ KL( q_θ(θ)π(θ)) + 
        log C_Q}.

    We consider the KL divergence part first. Recall that priors on w and η are independent, which leads to
    
    D_ KL(q_θ(θ)π(θ)) = D_ KL(q_ w( w)π( w)) + D_ KL(q_η(η)π(η)).

    For the D_ KL(q_ w( w)π( w)) term, we can use the closed forms of the prior and the variational posterior of w to explicitly calculate
    
    ∫q_ w( w) log w_k dw = Ψ( n_k + ϕ_0) - Ψ(n + Kϕ_0),

    where Ψ(x) = Γ'(x)/Γ(x) is the so-called di-gamma function. This further leads to
    
        D_KL(q_w(w)π(w)) = ∑^K_k=1 (n_kΨ(n_k +     ϕ_0) -  logΓ(n_k + ϕ_0)) - nΨ(n+Kϕ_0)
      + logΓ(n + Kϕ_0) + logΓ(ϕ_0)^K/Γ(Kϕ_0),
    
    where the last term causes the undesirable behavior under large ϕ_0 as discussed after Theorem <ref>. To further simplify this expression, we resort to the following two inequalities: for any x>0 <cit.>, we have
    
    1/2x < log x - Ψ(x) < 1/x,

    and
    
    0≤logΓ(x) - ((x-1/2)log x - x +1/2log 2π) ≤1/12x.

    Applying these two inequalities to (<ref>), we can obtain
    
    |D_ KL(q_ w( w)π( w)) - (Kϕ_0 - 1/2)log n + (ϕ_0 -1/2)∑_k=1^K log (n_k + ϕ_0)| < C,

   for some constant C independent of n and ϕ_0.
   
    As to the D_ KL(q_η(η)π(η)) term, since both the prior π and the variational posterior of η_k are factorized under our setup, we have
    
    D_ KL(q_η(η)π(η)) = ∑_k=1^K D_ KL(q_η_k(η_k)π(η_k)).

    For each fixed k∈[K], we denote the variational posterior mode (i.e., maximizer of its density function) as
    
    η_k = max_η_k ∑_i=1^n p_ik(η_k^TT(x_i) + T_0(x_i) - A(η_k))_ℓ_nk(η_k),

    Note that the critical point η_k satisfies the first order condition ∇_η_kℓ_nk(η_k) = 0. Under this notation, we can decompose D_ KL(q_η_k(η_k)π(η_k)) as
    
    D_ KL(q_η_k(η_k)π(η_k))     = ∫ q_η_k(η_k)logq_η_k(η_k)/π(η_k)
        = ∫ q_η_k(η_k) logπ(η_k)exp{ℓ_nk(η_k) - ℓ_nk(η_k)}/π(η_k)∫π(η_k)exp{ℓ_nk(η_k) - ℓ_nk(η_k)}dη_k
        =∫ q_η_k(η_k)(ℓ_nk(η_k) - ℓ_nk(η_k))dη_k 
             - log∫π(η_k)exp{ℓ_nk(η_k) - ℓ_nk(η_k)}dη_k.

    For the second term, applying Taylor expansion of ℓ_nk(η_k) at η_k up to the second order, we obtain
    
    ℓ_nk(η_k) - ℓ_nk(η_k) = -1/2∑_i=1^n p_ik(η_k -η_k)^TI(η_ik)   (η_k - η_k),

    where η_ik lies on the segment between η_k and η_k such that η_ik = η_k + t_i(η_k - η_k),t_i∈[0,1]. Then a standard Laplace's approximation for integrals, e.g., equation (6.4.35) in <cit.>, implies that with n_k large enough, the second integral in (<ref>) satisfies
    
    |∫π(η_k)exp{ℓ_nk(η_k) - ℓ_nk(η_k)}dη_k/(2π/n_k)^d/2-π(η_k)/√(I(η_k))|≤C/n_k.

    Therefore, for large n_k, we have
    
    |- log∫π(η_k)exp{ℓ_nk(η_k) - ℓ_nk(η_k)}dη_k + d/2log n_k| ≤ C.

    For the first term, by plugging in the explicit form of q(η) from (<ref>) into the integral, we obtain
    
    ∫ q_η_k(η_k)(ℓ_nk(η_k) - ℓ_nk(η_k))dη_k 
            =     ∫π(η_k)(ℓ_nk(η_k) - ℓ_nk(η_k))exp{ℓ_nk(η_k) - ℓ_nk(η_k)}dη_k/∫π(η_k)exp{ℓ_nk(η_k) - ℓ_nk(η_k)}dη_k,

    whose denominator is already analyzed. Using Laplace's approximation again to the numerator, as well as the fact that the value of the integrand at η_k vanishes, we obtain
    
    |∫π(η_k)(ℓ_nk(η_k) - ℓ_nk(η_k))exp{ℓ_nk(η_k) - ℓ_nk(η_k)}dη_k/(2π/n_k)^d/2-0|≤C/n_k.

    Therefore, by noticing ℓ_nk(η_k) ≤ℓ_nk (η_k), we have
    
    -C/n_k≤∫ q_η_k(η_k)(ℓ_nk(η_k) - ℓ_nk(η_k))dη_k≤ 0.

    Combining (<ref>), (<ref>) and (<ref>) together, we finally obtain that for any sufficiently large n_k,
    
    |D_ KL(q_η_k(η_k)π(η_k)) - d/2log n_k|≤ C.


    Step 2: Next we prove a lower bound for ℒ(q_Z^n). According to the characterization of ELBO in (<ref>), every instance of {p_ik} gives a lower bound to ℒ(q_Z^n). Let us consider the following two special constructions:


    Construction (I): This construction applies to the case when ϕ_0>(d+1)/2.  For 1≤ k≤ K^*-1, we can always choose {p_ik: i∈[n],k∈[K^*-1]} such that
    
    n_k = ∑_i=1^n p_ik = w_k^*(n + Kϕ_0)     -ϕ_0, i.e. n_k+ϕ_0/n+Kϕ_0 = w_k^*,
      ∑_i = 1^n p_ikT(x_i) = n_k  ∇_η_k    A(η_k^*), i.e. η_k = η^*_k,

    and for K^*≤ k≤ K, we choose {p_ik: i∈[n],K^*≤ k≤ K} such that
    
    n_k = w^*_K^*n+Kϕ_0/K-K^*+1-ϕ_0, i.e.     (n_k+ϕ_0)(K-K^*+1)/n+Kϕ_0 = w_K^*^*
      ∑_i = 1^n p_ikT(x_i)  = n_k∇_η_k    A(η_K^*^*),  i.e. η_k   = η^*_K^*.

    With this construction, all n_k's are proportional to n, so we have from (<ref>) that
    
    D_ KL(q_η(η)π(η)) ≤d/2∑_k=1^K log n_k + C.

    Combining the above with the approximation (<ref>) to D_ KL(q_ w( w)π( w)), we have the following upper bound to the KL divergence between q_θ(θ) and π(θ),
    
    D_ KL(q_θ(θ)π(θ))     ≤ (Kϕ_0 - 1/2)log n + d/2∑_k=1^K log n_k + (1/2-ϕ_0) ∑_k=1^K log (n_k + ϕ_0) + C
        ≤ (Kϕ_0 - 1/2)log n + (d + 1/2 - ϕ_0) ∑_k=1^K log n_k + C
        = (Kϕ_0 - 1/2)log n + (d + 1/2 - ϕ_0) K log n + C
        = dK+K-1/2log n + C.

    As for the lower bound of log C_Q, we can first rewrite C_Q as
    
        C_Q     = ∏_i=1^n∑_s_iexp∫q_θ(θ)logp(x_i, s_i|θ) dθ

            = ∏_i=1^n∑_k=1^Kexp(∫q_w_k(w_k)logw_k dw_k + ∫q_η_k(η_k)logg(x_i;η_k) dη_k).
    
    Using equations (<ref>) and (<ref>) again, we obtain
    
    ∫q_w_k(w_k) log w_k dw_k = Ψ(n_k + ϕ_0) - Ψ(n + Kϕ_0)≥logn_k+ϕ_0/n+Kϕ_0-1/n_k+ϕ_0.

    Additionally, since η_k^* = η_k in this construction, we can apply a Taylor expansion at η_k^*. Concretely, we can write q_η_k(η_k) explicitly according to (<ref>) as,
    
    ∫q_η_k(η_k)log g(x_i;η_k) dη_k     = ∫π(η_k)log g(x_i;η_k)exp{ℓ_nk(η_k) - ℓ_nk(η_k^*)}dη_k/∫π(η_k)exp{ℓ_nk(η_k) - ℓ_nk(η_k^*)}dη_k.

    Under Assumption A1 and using the fact that n_k≍ n, we can apply Laplace's approximation to approximate the numerator and denominator respectively as,
    
    |∫π(η_k)log g(x_i;η_k)exp{ℓ_nk(η_k) - ℓ_nk(η_k^*)}dη_k/(2π/n_k)^d/2·√(I(η_k^*))/π(η_k^*)log g(x_i;η_k^*) - 1 |≤C/n_k,

    and
    
    |∫π(η_k)exp{ℓ_nk(η_k) - ℓ_nk(η_k^*)}dη_k/(2π/n_k)^d/2·√(I(η_k^*))/π(η_k^*) - 1|≤C/n_k.

    Therefore, with the choice of w_k^* = (n_k +ϕ_0)/(n + Kϕ_0) for all k<K^*, we have
    
    ∑_k=1^K^*-1exp(∫q_w_k(w_k) log w_k dw_k + ∫q_η_k(η_k)log g(x_i;η_k) dη_k)
    ≥    ∑_k=1^K^*-1exp(log w_k^* -1/n_k+ϕ_0 +log g(x_i;η_k^*)-log(1+C/n_k))
    ≥    ∑_k=1^K^*-1 w_k^* g(x_i;η_k^*)(1-C/n).

    The same argument can applies to k≥ K^*, where by using the fact that n_k≍ n, we can obtain that for all K≥ K^*, 
    
    ∑_k=K^*^K exp(∫q_w_k(w_k) log w_k dw_k + ∫q_η_k(η_k)log g(x_i;η_k) dη_k)
    ≥    ∑_k=K^*^K exp(logw_K^*^*/K-K^*+1 -1/n_k+ϕ_0 +log g(x_i;η_K^*^*)-log(1+C/n_k))
    ≥    w_K^*^* g(x_i;η_K^*^*) (1-C/n).

    According to (<ref>), the constant log C_Q can be bounded from below by
    
    log C_Q ≥∑_i=1^nlog(∑_k=1^K^* w_k^* g(x_i;η_k^*)(1-C/n)) = log p^*(X^n) - C.

    Finally, by noticing the formulation of ℒ(q_Z^n) as in (<ref>), we can add (<ref>) and (<ref>) together to get a lower bound when ϕ_0>(d+1)/2 as
    
    ℒ(q_Z^n) - log p^*(X^n)≥ -dK+K-1/2log n -C.


     Construction (II): This construction applies to the case when  ϕ_0≤(d+1)/2. For each k such that 1≤ k≤ K^*, we can always choose  {p_ik: i∈[n],k∈[K^*]} such that
    
    n_k = ∑_i=1^n p_ik = w_k^*(n + ϕ_0)    -ϕ_0,    i.e. n_k +ϕ_0/n+Kϕ_0 = w_k^*,
      ∑_i = 1^n p_ikT(x_i)    = n_k∇_η_kA(η_k^*);

    and for K^*< k≤ K, we set n_k = 0. Then for each k>K^*, q_η_k(η_k) is the same as prior density π(η_k), so we have D_ KL(q_η_k(η_k)π(η_k)) = 0 for all k>K^*. For each k≤ K^*, we have n_k≍ n, so (<ref>) leads to
    
    D_ KL(q_η(η)π(η))  ≤∑_k=1^K^*d/2log n_k + C.

   Combining the preceding display with (<ref>), we obtain
    
    D_ KL(q_θ(θ)π(θ)) ≤(dK^* + K^*-1/2 + (K-K^*)ϕ_0)log n + C.

    Now since the summation over k∈[K] for the C_Q in (<ref>) is dominated by its first K^* terms, we can lower bound this summation to the sum of the first K^* terms. Now, since n_k≍ n for k≤ K^*, we can apply (<ref>) and (<ref>) to obtain
    
    ∑_k= 1^K^*exp(∫ q_θ_k(θ_k) log(w_kg(x_i;η_k))dθ_k )
            ≥    ∑_k=1^K^*w_k^*g(x_i;η^*_k)(1-C/n).

    Putting all the pieces together, we obtain
    
    log C_Q ≥∑_i=1^nlog(∑_k=1^K^* w^*_kg(x_i;η^*_k)) + nlog(1-C/n) = log p^*(X^n) - C + 𝒪(1/n).

    Adding (<ref>) and (<ref>) together, we can obtain a lower bound to ℒ(q_Z^n(z^n)) for all ϕ_0≤(d+1)/2 as
    
    ℒ(q_Z^n) - log p^*(X^n)≥ -(dK^*+K^*-1/2 + (K-K^*)ϕ_0)log n -C.

    Finally, combining (<ref>) and (<ref>) leads to the claimed lower bound of ℒ(q_Z^n).

    Step 3:
    Now we derive an upper bound of the ELBO in (<ref>). We first tackle the logC_Q term, recall that C_Q is the normalization constant for the optimal distribution q_S^n(s^n) in (<ref>). Note that log w_k is concave with respect to w_k. For log g(x;η_k), the second order derivative with respect to η_K is the negative fisher information matrix for the mixture component
    
    ∇^2_η_klog g(x;η_k) = -∇^2_η_kA(η_k) = -I(η_k),

    which is negative definite according to Assumption A1. Therefore, log g(x;η) is concave with respect to η_k. By applying Jensen's inequality, we can then obtain
    
    ∫q_w_k(w_k) log w_k  dw_k + ∫q_η_k (η_k) log g(x;η_k)  dη_k≤log (w_kg(x;η_k)),

    which combined with (<ref>) leads to
    
    logC_Q ≤∑_i=1^nlog(∑_k=1^K w_kg(x_i;η_k)) = log p(X^n | θ).

    The KL divergence D_ KL(q_w(w)π(w)) is the same as (<ref>). As to the D_ KL(q_η(η)π(η)) term, note that for a large n_k, we can apply the approximation (<ref>). Consequently, we define a set ℱ:= {k  |  1 ≤ k ≤ K, n_k ≥ M }, where M is a sufficiently large (truncating) constant, to select those indices k for which the approximation is accurate. This approximation, when combined with the non-negativity of the KL divergence, gives
    
    D_ KL(q_η(η)π(η)) = ∑_k=1^K D_ KL(q_η_k(η_k)π(η_k))≥d/2∑_k∈ℱlogn_k -C.

    Using approximation (<ref>) again, we can bound the entire KL divergence term from below as
    
    D_ KL(q_θ(θ)π(θ))     ≥ (Kϕ_0-1/2)log n + (d+1/2-ϕ_0)∑_k ∈ℱlogn_k 
                                 + (1/2-ϕ_0)∑_k ∉ℱlogn_k -C 
        = (Kϕ_0-1/2)log n + (d+1/2-ϕ_0)∑_k ∈ℱlogn_k - C.


    In the last part of this proof, we will show that at least K^* of the n_k's are proportional to n, thus are in the set ℱ. Using this fact, we can obtain that for each ϕ_0<(d+1)/2,
    
    D_ KL(q_θ(θ)π(θ))     ≥(Kϕ_0-1/2)log n + (d+1/2-ϕ_0)K^*log n - C
        =(dK^*+K^*-1/2 + ϕ_0(K-K^*))log n -C.
 
    When ϕ_0 > (d+1)/2, by using the fact that 0<n_k≤ n and the lower bound in (<ref>), we can obtain,
    
    D_ KL(q_θ(θ)π(θ))     ≥ (Kϕ_0-1/2)log n + (d+1/2-ϕ_0)∑_k=1^Klogn_k - C
        ≥ (Kϕ_0-1/2)log n + (d+1/2-ϕ_0)Klog n - C
        = dK+K-1/2log n - C.

    Combining the two lower bounds of D_ KL(q_θ(θ)π(θ)) obtained above and the bound in (<ref>), we can get an upper bound to ℒ(q_Z^n) - log p^*(X^n) as 
    
    ℒ(q_Z^n) - log p(X^n | θ) ≤ -λlog n + C,

    where λ is given by equation (<ref>) in the theorem statement.


    It remains to show that at least K^* of the n_k's are of the same magnitude of n. Noticing that the Hellinger distance
    
    h^2(p(x | θ), p^*(x))     = ∫((∑_k=1^K w_kg(x;η_k))^1/2 - √(p^*(x)))^2  d x
        = ∫{∑_k=1^K w_kg(x;η_k) + p^*(x) - 2(∑_k=1^K w_kg(x;η_k))^1/2√(p^*(x))} dx

    is a convex function with respect to every w_k, then using (<ref>), Lemma <ref>, Jensen's inequality and the inequality d_ TV(p(x | θ),  p^*(x)) ≤ h(p(x | θ),  p^*(x)), we can get
    
    (inf_{I_1,...I_K^*}∑^K^*_k=1|∑_j∈ I_kw_j - w_k^*|)^2    ≲∫_Ω^K h^2(p(x | w,η), p^*(x)) q(η) dη
        ≤∫_Θ h^2(p(x |  w,η), p^*(x)) q(θ) dθ.
        ≤ C_4 ( log p^*(X^n) - ℒ(q_Z^n)/n + log n/n),

    where the convexity is used in the second inequality. Recall that, in Step 2, we proved a lower bound to ℒ(q_Z^n) - log p^*(X^n), which is -λlog n. Thus we can obtain
    
    (inf_{I_1,...I_K^*}∑_k = 1^K^*|∑_j∈ I_kw_j - w^*_k|)^2 ≲λlog n + log n/n.

    Since the cardinality of I_k is at most K, we have that for every k≤ K^*, there is at least one w_j,j∈ I_k such that w_j≥ (w^*_k - C√(log n/n))/K. Therefore, we can obtain
    
    n_j = w_j(n+Kϕ_0) - ϕ_0≥(w^*_k - C√(log n/n))(n + Kϕ_0)/K -ϕ_0≍ n,

    which completes the proof.





 §.§ Proof of Corollary <ref>

According to the classical Hironaka's Theorem for the resolution of singularity in algebraic geometry, which is Theorem 2.3 in <cit.>, the KL divergence K(θ) = D_ KL(p^*(·)  p(· | θ)) for a parametric family {p(· | θ): θ∈Θ} of a singular model, as a function of θ, can be transformed (via an invertible map) into its normal crossing form (uniquely determined by analytical properties of K(θ)) around any of its zero point θ^*, that is, any point θ^* such that K(θ^*)=0; see Theorem 2.3 in <cit.> for the precise statement. Moreover, Theorem 6.2 in <cit.> characterized the asymptotic behavior of the likelihood ratio statistic under a standard normal crossing form. Note that their Fundamental condition (I) (Definition 6.1) is satisfied by our considered mixtures of the exponential family under Assumptions A1–A3. 

Moreover, we know that the maximum likelihood estimator is invariant under any (one-to-one) reparameterization of the model (i.e., the invariance property of MLE). 
Therefore, the likelihood ratio statistic LRT defined after Theorem <ref> coincides with the likelihood ratio statistic under the normal crossing form analyzed in their Theorem 6.2, which is shown to converge in distribution as n→∞ to the supremum of a Gaussian process that is bounded in probability. These facts imply a bound of LRT=𝒪_P(1), which combined with Theorem <ref> leads to Corollary <ref>.




 §.§ Proof of Corollary <ref>

    According to Theorem <ref> and the bound Δ_n ≤ LRT = 𝒪_P(1) as n→∞, we can conclude that for all sufficiently large n and all K≥ K^*, the probability that the ELBO value ℒ_K strictly decreases with K approaches one as n→∞. Particularly, for K = K^*, we have
    
    log p^*(X^n) - dK^* + K^* - 1/2log n≤ℒ(q_Z^n)-C_1≤log p(X^n | θ) - dK^* + K^* - 1/2log n+C_2.

    On the other hand, when K<K^* and G(θ) = ∑_j=1^Kw_jδ_η_j, we will show below that the Wasserstein distance will be bounded from below by a constant order; then we can apply Lemma <ref> to show that the ELBO value ℒ_K is also strictly less than ℒ_K^*. Concretely, recall that in the proof of Lemma <ref>, we defined δ = min_j≠ k|η^*_j - η^*_k|>0 and I_k = {j:|η_j - η_k^*|<δ/2} for k∈[K]. Since K<K^* and I_k's are disjoint, there is at least one I_k that is empty. Without loss of generality, we may assume I_1 = ∅. This then leads to (according to the proof of Lemma <ref>)
    
    W_r^r(G,G^*)≥inf_q_i1∑_i=1^K q_i1|η_i-η_1^*|^r ≥inf_q_i1∑_i=1^K q_i1(δ/2)^r = w_1^*(δ/2)^r.

    Now an application of Lemma <ref>, Assumption A3 and the inequality h(p(x | θ), p^*(x))≥ d_ TV(p(x | θ), p^*(x)) leads to
    
    log p^*(X^n) - ℒ(q_Z^n)/n + log n/n≳∫_Θ d_ TV^2(p(x | θ), p^*(x))  q_θ(θ)≳w_1^*^2(δ/2)^2r.

    Therefore, we have
    
    ℒ(q_Z^n)≤log p^*(X^n) + log n - c nw_1^*^2(δ/2)^2r

    for some positive constant c. On the other hand, the ELBO value ℒ_K^* is at least log p^*(X^n) - 𝒪_P(log n), which  implies that the ELBO value ℒ_K at K<K^* is strictly smaller than that at K=K^* for all sufficiently large n.




 §.§ Proof of Corollary <ref>

    Under the assumption that T_n(θ)≤ C_3loglog n holds with high probability, we have
    
    -λlog n - C_3loglog n - C_1≤ℒ(q_Z^n) - log p(X^n | θ)≤-λlog n + C_2.


    Let us start with the case of ϕ_0<(d+1)/2.
    In the proof of Theorem <ref>, we showed that at least K^* of the n_k's are proportional to n when ϕ_0<(d+1)/2. We can then without loss of generality assume that it is the first K^* of n_k's. 
    We will prove the claimed bound by contradiction below. Now suppose there exists some k'>K^* such that w_k'≥ [(log n)^ρ_1 + ϕ_0]/n, then n_k' = w_k'(n + Kϕ_0) - ϕ_0 ≥ (log n)^ρ_1, implying k'∈ℱ (defined in the proof of Theorem <ref>). Now using inequality (<ref>), we obtain
    
    D_ KL(q_θ(θ)π(θ))     ≥ (Kϕ_0-1/2)log n + (d+1/2-ϕ_0)∑_k ∉ Flog n_k -C
        ≥[(K-K^*)ϕ_0 + dK^*+K-1/2]log n + (C_3+1)loglog n- C.

    Therefore, by further using (<ref>) and (<ref>), we obtain 
    
    ℒ(q_Z^n) - log p(X^n | θ) ≤ -  D_ KL(q_θ(θ)π(θ)) + C ≤ - λ_1 log n - (C_3+1)loglog n - C,

    where λ_1 = (K-K^*)ϕ_0 + (dK^*+K^*-1)/2. Note that this upper bound is smaller than the corresponding lower bound as in inequality (<ref>) for all sufficiently large n, which is a contradiction. Therefore, we must have w_k ≤ [(log n)^ρ_1 + ϕ_0]/n for all k>K^*, which leads to
    
    inf_σ∈𝒮_k∑_k = K^*+1^K w_σ(k)≤∑_k = K^*+1^K w_k <(K-K^*)[(log n)^ρ_1 + ϕ_0]/n.


    For the other case of ϕ_0>(d+1)/2, suppose there exists some k” such that w_k≤ 1/(log n)^ρ_2 + ϕ_0/n, then we have n_k”≤ n/(log n)^ρ_2. By further using inequality (<ref>), we get
    
    D_ KL(q_θ(θ)π(θ))     ≥ (Kϕ_0-1/2)log n + (d+1/2-ϕ_0)(Klog n - ρ_2loglog n) -C,

    which implies that, with λ = (dK+K-1)/2, we have
    
    ℒ(q_Z^n) - log p(X^n | θ)≤ - λlog n - (C_3 + 1)loglog n + C.

    This upper bound is again smaller than its corresponding lower bound from (<ref>), which is a contradiction. Therefore, for all k∈{1,...,K}, we must have w_k≥ 1/(log n)^ρ_2 + ϕ_0/n. This completes the proof.




 §.§ Proof of Theorem <ref>
 

Step 1: In the first step, we establish a relationship between the total variation distance between two finite mixture distributions and a discrepancy metric between their respective parameters, similar to the proof of Theorem 3.1 in <cit.>. Concretely, we consider a general finite mixture measure p(x | θ) = ∑_k=1^K^* w_kg(x;η_k), with θ=(w_1,…,w_K^*,η_1,…,η_K^*), and w = (w_1,…,w_K^*), η = (η_1,…,η_K^*) where the mixing weights satisfy w_k≥ 0 and ∑_k=1^K^*w_k ≤ 1 and every η_k is from the compact set Ω with diameter L<∞. Note that here p(x|θ) is not necessarily a density function since its integral can be smaller than 1. Also recall that p^*(x) = ∑_k=1^K^* w_k^* g(x;η_k^*) is the true finite mixture distribution. We consider the following discrepancy metric between the parameters of p(· | θ) and p^*(·),
    
    D(θ) := inf_σ∈𝒮_K^*∑_k=1^K^*(|w_σ(k) - w_k^*| + w_σ(k)|η_σ(k) - η_k^*|).

    We will show below that there exists positive constants ε and c such that for any θ such that D(θ) < ε, it holds 
    
    1/2∫|p(x|θ) - p^*(x)|dx > cD(θ).

    
    
    We will prove this inequality by contradiction. If this is not the case, then there exists some sequence {θ^m: m≥ 1} such that D(θ^m)→ 0 as m→∞ and 
    
    lim_m→∞1/ D(θ^m) ∫ |p(x|θ^m) - p^*(x)|dx = 0.

    Since D(θ^m)→ 0 as m→∞ and the number of possible permutations of [K^*] is finite, there is a sub-sequence {θ^m_ℓ} and a permutation σ'∈𝒮_K^* such that w^m_ℓ_σ'(k)→ w^*_k and η^m_ℓ_σ'(k)→η^*_k for all k = 1,..., K^* as ℓ→∞. Without loss of generality, we may assume that σ' is the identity permutation and substitute the sequence with its sub-sequence, so that we can write D(θ^m) = ∑_k=1^K^*(|w_k - w_k^*| + w_k|η_k - η_k^*|) for all sufficiently large m. Now we can apply a Taylor expansion as follows to obtain:
    
    p(x|θ^m) - p^*(x)     = ∑_k=1^K^* w_k^m (g(x;η^m_k) - g(x;η_k^*)) + ∑_k=1^K^*(w_k^m - w_k^*)g(x;η_k^*)
        = ∑_k=1^K^* w_k^m [(η^m_k - η_k^*)∇ g(x;η_k^*) + O(η^m_k - η_k^*)^2 ] + ∑_k=1^K^*(w_k^m - w_k^*)g(x;η_k^*)
        = ∑_k=1^K^* w_k^m (η^m_k - η_k^*)∇ g(x;η_k^*) + ∑_k=1^K^*(w_k^m - w_k^*)g(x;η_k^*) + R_m,

   where the remainder term satisfies R_m/D(θ^m)→ 0 as m→∞. Note that equation (<ref>) will then imply
    
    lim_m→∞p(x|θ^m) - p^*(x)/D(θ^m) = lim_m→∞{∑_k=1^K^*w_k^m (η^m_k - η_k^*)/D(θ^m) ∇ g(x;η_k^*) + ∑_k=1^K^*w_k^m - w_k^*/D(θ^m) g(x;η_k^*)} = 0

    for almost every x. 
    Note that the term inside the second limit is a linear combination of finite many g(x;η_k^*) and ∇ g(x;η_k^*) terms. Since the parameter space is compact under Assumption A1, each linear combination coefficient has a convergent sub-sequence. Without loss of generality, we may substitute the sequence with its sub-sequence, and denote the respective limits as
    
    α_k = lim_m→∞w_k^m (η^m_k - η_k^*)/D(θ^m)    β_k = lim_m→∞w_k^m - w_k^*/D(θ^m).

    Note that at least one of these limiting coefficients {α_k}_k=1^K^* and {β_k}_k=1^K^* does not vanish, since otherwise, we will have
    
    1= lim_m→∞D(θ^m)/D(θ^m) = lim_m→∞∑_k=1^K^*w_k^m|η^m_k - η_k^*| + ∑_k=1^K^*|w_k^m - w_k^*|/D(θ^m) = ∑_k=1^K^*(|α_k| + |β_k|) =  0,

    which is a contradiction. Therefore, we constructed some {α_k}_k=1^K^* and β=(β_1,…,β_K^*) such that they are not all zero, and 
    
    ∑_k = 1^K^*[α_k g(x;η_k^*) + β^T ∇ g(x;η_k^*)] = lim_m→∞p(x|θ^m) - p^*(x)/D(θ^m) = 0

    for almost every x. Note that the preceding display and the weak identifiability condition together would imply that α_k = 0, β_k = 0⃗ for all k = 1,..., K^*, which contradicts the fact that at least one of the coefficients is not zero. Therefore, (<ref>) holds. Since if 1/2∫ |p(x|θ) - p^*(x)|dx = 0, then ∑_k=1^K^*w_k = 1 and by Assumption A3, D(θ) = 0; if D(θ) = 0, it also implies ∑_k=1^K^*w_k = 1 and from Assumption A3 again we have 1/2∫ |p(x|θ) - p^*(x)|dx = 0. Therefore, 1/2∫ |p(x|θ) - p^*(x)|dx = 0 if and only if D(θ) = 0. Because both D(θ) and 1/2∫ |p(x|θ) - p^*(x)|dx are continuous in θ, we can further obtain (following the same argument from which we derived (<ref>)) from equation (<ref>) that for each possible θ,
    
    1/2∫ |p(x|θ) - p^*(x)|dx ≳ D(θ).

  

    Step 2: Recall that we have already shown in Theorem <ref> and Corollary <ref> that when ϕ_0<(d+1)/2, there are K^* of the mixing weights w_j's being bounded away from 0, and the rest K-K^* of the mixing weights w_j's will empty out at the rate (log n)^ρ_1/n. Without loss of generality, we may assume that the first K^* of the weights w_j's do not vanish and denote w_K^* = (w_j)_j = 1^K^* and η_K^* = (η_j)_j = 1^K^*. Then for any fixed ϕ_0<(d+1)/2, we have
    
    d_ TV(p(· |  w,η),p^*(·))     = 1/2∫|∑_j=1^K w_jg(x;η_j) - ∑_k=1^K^* w_k^* g(x;η_k^*) | dx
        ≥1/2∫|∑_j=1^K^*w_jg(x;η_j) - ∑_k=1^K^* w_k^* g(x;η_k^*) | dx - 1/2∑_j = K^* + 1^K w_j
        ≥1/2∫|∑_j=1^K^*w_jg(x;η_j) - ∑_k=1^K^* w_k^* g(x;η_k^*) | dx - C(log n)^ρ_1/n
        ≳ D( w_K^*,η_K^*) - C(log n)^ρ_1/n.

    where the last step is due to inequality (<ref>). Recall that we have sup_η,η'∈Ω|η-η'|≤ L from Assumption A1, which leads to the bound D( w_K^*,η_K^*)≤∑_k=1^K^*(w_k + w_K^* + w_kL)≤ 2 + L. Now since (C(log n)^ρ_1/n)^2>0, we can then get 
    d^2_TV(p(· |  w,η),p^*(·))≳  D^2( w_K^*,η_K^*) - 2C(L + 2)(log n)^ρ_1/n.


    Step 3: In this step, we show that with high probability, the variational posterior distribution of each non-emptying out mixture component parameter η_k is concentrated within a neighborhood of η_k with radius Mn^-1/2. 
     Recall that from Corollary <ref>, for each non-emptying out component k, we have n_k≥ c_kn for some positive constant c_k. In this step, we may again assume without loss of generality that the first K^* of the weights w_j's do not vanish. 
     Recall that from Theorem <ref>, we have that
    
    q_η_k(η_k)∝π(η_k)exp[∑_i=1^n p_ik(η_k^TT(x_i) - T_0(x_i) - A(η_k))],

    where ∑_i=1^n p_ik = n_k. Recall that η_k, defined in equation (<ref>), denotes the η_k that maximizes the exponent in the preceding display. We will first show below that 
    
    min_1≤ j≤ K^*|η_j - η^*_k|≲log n/√(n_k)

    holds for each 1≤ k≤ K^*. 
    
    From equations (<ref>) and (<ref>), and Assumption A2, we can get
    
    π(η_k)exp[∑_i=1^n p_ik(η_k^TT(x_i) - T_0(x_i) - A(η_k))- ℓ_nk(η_k)]≥π(η_k)exp[-bn_k/2(η_k-η_k)^2],

    which implies that for any positive number C,
    
    ∫_{|η_k-η_k|_∞≤ C(bn_k)^-1/2}q_η_k(η_k)dη_k ≥inf_{|η_k-η_k|_∞≤ C(bn_k)^-1/2}π(η_k)/π(η_k)(a(2Φ(C) - 1)/b)^d/2,

    where Φ denotes the CDF of the standard normal distribution. In particular, by choosing C = Φ^-1(3/4), the integral above can be bounded from below by c_π(a/2b)^d/2, for some positive constant c_π when n is sufficiently large. Therefore, for constant C_k = C(d/bc_k)^1/2, we have the following inequality,
    
    ∫_{|η_k-η_k|_∞≤ C_kn^-1/2}q_η_k(η_k)dη_k ≥ c_π(a/2b)^d/2.

 We can now prove (<ref>) by contradiction.
    Suppose min_1≤ j≤ K^*|η_j - η^*_k'|≥max_1≤ j≤ K^*C_jlog n /n^1/2 holds for some k'≤ K^*, then each η_j in the set ℬ := ⋂_j=1^K^*{η:|η_j - η_j|≤ C_jn^-1/2} cannot be close to η_k^* by the triangle inequality, i.e.,
    
    min_1≤ j≤ K^*|η_j - η^*_k|≥min_1≤ j≤ K^*(|η_j - η^*_k'| - |η_j - η^*_k|)≥min_1≤ j≤ K^*C_j(log n - 1)/n^1/2

    holds for any η=(η_j: j∈[K^*]) ∈ℬ. Since from Step 1, we obtain d^2_TV(p(x| w,η),p^*(x))≳ D^2( w_K^*,η_K^*) - 2C(L+2)(log n)^ρ_1/n, which furhter implies
    
    ∫_Ω^K h^2(p(x |  w,η), p^*(x)) q(η) dη    ≥∫_ℬ d_ TV^2(p(x |  w,η), p^*(x)) q(η) dη
        ≳∫_ℬ D^2( w_K^*,η_K^*) q(η) dη - 2CL(log n)^ρ_1/n
        ≳ (log n -1 )n^-1/2Q(ℬ) - 2CL(log n)^ρ_1/n
        ≥ (log n -1 )  n^-1/2(c_π(a/2b)^d/2)^K^* - 2CL(log n)^ρ_1/n.

    This is a contradiction with Lemma <ref>, which implies the left-hand side of the preceding display to be at most 𝒪_P(n^-1). This proves inequality (<ref>). 
    
    Lastly, we show that the variational posterior mode η_k is close to the variational posterior mean η_k for each k∈[K^*]. In fact, since η_k is given by
    
    η_k = ∫η_kπ(η_k) exp{ℓ_nk(η_k) - ℓ_nk(η_k)}dη_k/∫π(η_k) exp{ℓ_nk(η_k) - ℓ_nk(η_k)}dη_k,

    we may use the Laplace formula again to obtain
    
    |∫η_kπ(η_k) exp{ℓ_nk(η_k) - ℓ_nk(η_k)}dη_k/(2π/n_k)^d/2-η_kπ(η_k)/√(I(η_k))|≤C/n_k.

    This inequality combined with (<ref>) leads to 
    |η_k - η_k| ≤ C/n_k; and the conclusion follows by combining this with (<ref>) and using the fact that n_k≍ n for k∈[K^*].

