Consistency of Variational Bayes Inference for Estimation and Model Selection in Mixtures

Introduction

This paper studies the statistical properties of variational inference as a tool to tackle two problems of interest: estimation and model selection in mixture models. Mixtures are often used for modelling population heterogeneity, leading to practical clustering methods . Moreover they have enough flexibility to approximate accurately almost every density . Mixtures are used in many various areas such as computer vision , genetics , economics , transport data analysis and others. We refer the reader to for an account of the recent advances on mixtures. The most famous procedure for mixture density estimation in the frequentist literature is probably Expectation-Maximization , a maximum-likelihood algorithm that yields increasingly higher likelihood. At the same time, the Bayesian paradigm has raised great interest among researchers and practitioners, especially through the Variational Bayes (VB) framework which aims at maximizing a quantity referred to as Evidence Lower Bound on the marginal likelihood (ELBO). Variational Bayes inference is a useful tool for approximating intractable posteriors. It is known to work well in practice for mixture models: one of the most recent survey on VB chooses mixtures as an example of choice to illustrate the power of the method. Moreover states: "the [evidence lower] bound is a good approximation of the marginal likelihood, which provides a basis for selecting a model. Though this sometimes works in practice, selecting based on a bound is not justified in theory". The main contribution of this paper is to prove that VB is consistent for estimation in mixture models, and that the ELBO maximization strategy used in practice is consistent for model selection. Thus we solve the question raised by .

Variational Bayes is a method for computing intractable posteriors in Bayesian statistics and machine learning. Markov Chain Monte Carlo (MCMC) algorithms remain the most widely used methods in computational Bayesian statistics. Nevertheless, they are often too slow for practical uses when the dataset is very large. A more and more popular alternative consists in finding a deterministic approximation of the target distribution called Variational Bayes approximation. The idea is to minimize the Kullback-Leibler divergence of a tractable distribution $\rho$ with respect to the posterior, which is also equivalent to maximizing the ELBO. This optimization procedure is much faster and efficient than MCMC sampling with numerous applications in different fields: matrix completion for collaborative filtering , computer vision , computational biology and natural language processing , to name a few prominent examples.

However, variational inference is mainly used for its practical efficiency and only little attention has been put in the literature towards theoretical properties of the VB approximation until very recently. In the properties of variational approximations of Gibbs distributions used in machine learning are derived. The results are essentially valid for bounded loss functions, which makes them difficult to use beyond the problem of supervised classification. Based on some technical advances from , removed the boundedness assumption in , allowing to study more general statistical models. In , the authors extended the range of models covered by . This allowed them to study mixture of Gaussian distributions as an example. Many questions are still left unanswered: model selection, and the estimation of mixture of non-Gaussian distributions. For example mixture of multinomials are widely used in practice , as well as more intricated examples such as nonparametric mixtures . Note that all the results in are limited to so-called tempered posteriors, that is, where the likelihood is taken to some power $\alpha$. Still, the use of tempered posteriors is highly recommended by many authors as a way to overcome model misspecification, see and the references therein. Indeed some results in are valid in a misspecified setting. Alternative approaches were developed to study VB: established Bernstein-von-Mises type theorems on the variational approximation of the posterior. They provide very interesting results for parametric models but it is unclear whether these results can be extended to model selection or misspecified case. More recently, succeeded in adapting the now classical results of to Variational Bayes and showed that a slight modification in the three classical "prior mass and testing conditions" leads to the convergence of their variational approximations, again under the assumption that the model is true. With respect to these works, our contribution is a complete study of the consistency of VB for mixtures of general distributions. In particular, we explicit independent conditions on the prior on the weights, and on the prior on the parameters of the components. The study is done in the case $\alpha<1$ which allows to prove results in the misspecified case.

The other point addressed in this paper is model selection. This is a natural question which can be interpreted in this context as the determination of the number of components of the mixture. This point is crucial: indeed, too many components can lead to estimates with too large variances whereas with too few components, we may obtain mixtures which are not able to fit the data properly. This is a common issue and a lot of statisticians worked on this question. In the literature, criteria such as AIC and BIC are popular. It is well known that in some collections of models, AIC optimizes the prediction ability while BIC recovers with high probability the true model (when there is one). These two objectives are not compatible in general . Anyway, these results depend on asssumptions that are not satisfied by mixtures. It seems thus more natural to develop criteria suited to a given objective. For example, proposed a procedure to select a number of components that is the most relevant for clustering. A non-asymptotic theory of penalization has been developed during the last two decades using oracle inequalities . In the wake of those works, our paper studies mixture model selection based on the ELBO criterion. We prove a general oracle inequality. This result establishes the consistency of ELBO maximization when the primary objective is the estimation of the distribution of the data.

The rest of this paper is organized as follows. In Section 2 we introduce the background and the notations that will be adopted. Consistency of the Variational Bayes for estimation in a mixture model is studied in Section 3. First, we give the general results under a "prior mass" assumption, as well as a general form for the algorithm to compute the VB approximation (Subsection 3.1). We then apply these results to mixtures of multinomials (Subsection 3.2) and Gaussian mixtures (Subsection 3.3). In each case, we provide a rate of convergence of VB and discuss its numerical implementation. We extend the setting to the misspecified case in Subsection 3.4. Finally, we address the issue of selecting based on the ELBO in Section 4. We discuss possible extensions in Section 5, while Section 6 is dedicated to the proofs.

Background and notations

First, we consider the well-specified case, assuming that the true distribution belongs to the $K$-components mixture model. Thus, we define the true distribution $P^{0}$ from which data are sampled:

Hence, we want to estimate the true distribution $P_{\theta^{0}}$ using a Bayesian approach. Therefore, we define a prior $\pi=\pi_{p}\bigotimes_{j=1}^{K}\pi_{j}$ on $\theta$, $\pi_{p}\in\mathcal{M}_{1}^{+}(\mathcal{S}_{K})$ being a probability distribution on some measurable space $(\mathcal{S}_{K},\mathcal{A})$, and each $\pi_{j}\in\mathcal{M}_{1}^{+}(\Theta)$ a probability distribution on the measurable space $(\Theta,\mathcal{T})$. We will also consider in this paper the misspecified case where the true distribution does not belong to our statistical model i.e. is not necessarily a mixture, but the specific notations and framework will be described later.

The negative log-likelihood ratio $r_{n}$ between two distributions $P$ and $R$ is given by

(note that $r_{n}(\theta,\theta^{\prime})$ is used by many authors instead of $r_{n}(P_{\theta},P_{\theta^{\prime}})$ but our notation is more convenient for the extension to the misspecified case). The Kullback-Leibler (KL) divergence between two probability distributions $P$ and $R$ is given by

If some measure $\lambda$ dominates both $P$ and $R$ distributions represented here by their densities $f$ and $g$ with respect to this measure, we have

and we will use $K(P,R)$ or $\mathcal{K}(f,g)$ to denote this quantity, depending on the context.

We also remind that the $\alpha$-Renyi divergence between $P$ and $R$,

When for some $\lambda$ we have $f=\frac{dP}{d\lambda}$ and $g=\frac{dR}{d\lambda}$,

Some useful properties of Renyi divergences can be found in . In particular, the Renyi divergence between two probability distributions $P$ and $R$ can be related to the classical total variation $TV$ and Hellinger $H$ distances respectively defined as $TV(P,R)=\frac{1}{2}\int|dP-dR|$ and $H(P,R)^{2}={\frac{1}{2}\int(\sqrt{dP}-\sqrt{dR})^{2}}={1-e^{-\frac{1}{2}D_{1/2}(P,R)}}$ through:

The tempered Bayesian posterior $\pi_{n,\alpha}(.|X_{1}^{n})$, which is our target here, is defined for $0<\alpha\leq 1$ by

(it is also referred to as fractional posterior, for example in ). Note that when $\alpha=1$, then we recover the "true" Bayesian posterior, but the case $\alpha<1$ has many advantages: it is often more tractable from a computational perspective , it is consistent under less stringent assumptions than required for $\alpha=1$ and it is more robust to misspecification .

The mean-field approximation is very popular in the Variational Bayes literature. It is based on a decomposition of the space of parameters $\Theta_{K}$ as a product. Then $\mathcal{F}$ consists in compatible product distributions. Here, a natural choice is $\Theta_{K}=\mathcal{S}_{K}\times\Theta\times\dots\times\Theta$ and

We will work on this particular set in the following and we will often use $\rho$ instead of $\rho_{p}\bigotimes_{j=1}^{K}\rho_{j}$ to ease notation.

We end this section by recalling Donsker and Varadhan’s variational formula. Refer for example to for a proof (Lemma 1.1.3).

with the convention $\infty-\infty=-\infty$. Moreover, if $h$ is upper-bounded on the support of $\lambda$, then the supremum on the right-hand side is reached by the distribution of the form:

This technical lemma is one of the main ingredients for the proof of our results, but it is also very helpful to understand variational approximations. Indeed, for $\textbf{E}=\Theta_{K}$ and using the definition of ${\pi}_{n,\alpha}(.|X_{1}^{n})$, we get:

The quantity maximized in (1) is called the ELBO in the litterature (ELBO stands for Evidence Lower Bound), and many authors actually take this as the definition of VB .

In practice, the choice of $\alpha$ is not staightforward. Depending on the objective, some heuristic might be available: for example, proposed a nice method to calibrate $\alpha$ in order to get confidence intervals on a parameter of interest. More generally, cross-validation can give good results. We have to acknowledge that there is no universal method to calibrate $\alpha$. This could lead the reader to the idea that the proper Bayesian approach ($\alpha=1$) is simpler to use. We insist on the fact $\alpha=1$ can produce catastrophic results in case of misspecification , while we present below some results in the misspecified case with $\alpha<1$. We believe that the calibration of $\alpha$ is a very important research direction.

Variational Bayes estimation of a mixture

We start with a result for general mixtures. Later in this section we provide corollaries obtained by applying this theorem to special cases: mixture of multinomials and Gaussian mixtures.

As a special case, when there exists $r_{n,K}$ such that there is are distributions $\rho_{p,n}\in\mathcal{M}_{1}^{+}(\mathcal{S}_{K})$ and $\rho_{j,n}\in\mathcal{M}_{1}^{+}(\Theta)$ ($j=1,...,K$) such that for $j=1,...,K$

The proof is given in Section 6. This theorem provides the consistency of the Variational Bayes for mixture models as soon as (3) and (4) are satisfied. In , the authors use similar conditions ((3) and (4) in their Theorem 2.6), and show that they are strongly linked to the assumptions on the prior used by to derive concentration of the posterior. Thus they cannot be removed in general. Theorem 3.1 states that finding $r_{n,K}$ fulfilling (3) and (4) independently for the weights and for each component is sufficient to obtain the rate of convergence $Kr_{n,K}$ of the VB estimator towards the true distribution.

Clearly, the theorem cannot be directly extended to the case $\alpha=1$. As discussed above, the case $\alpha=1$ is studied in : it requires a testing condition in addition to the prior mass condition given by (3) and (4). In the case of mixtures, such a testing condition was studied in .

Note that there always exists a distribution $\rho_{p,n}\in\mathcal{M}_{1}^{+}(\mathcal{S}_{K})$ such that the two quantities corresponding to the weights $\int\mathcal{K}(p^{0},p)\rho_{p,n}(dp)$ and $\mathcal{K}(\rho_{p,n},\pi_{p})$ are bounded as required in Theorem 3.1 for $r_{n,K}=\frac{4\log(nK)}{n}$ when the chosen prior is a Dirichlet distribution $\pi_{p}=\mathcal{D}_{K}(\alpha_{1},...,\alpha_{K})$ under some minor restriction on $\alpha_{1},...,\alpha_{K}$. This result summarized below for any $K\geq 2$ helps find explicit rates of convergence for the VB approximation.

For $r_{n,K}=\frac{4\log(nK)}{n}$ and a prior $\pi_{p}=\mathcal{D}_{K}(\alpha_{1},...,\alpha_{K})\in\mathcal{M}_{1}^{+}(\mathcal{S}_{K})$ with $\frac{2}{K}\leq\alpha_{j}\leq 1$ for $j=1,...,K$, we can find a distribution $\rho_{p,n}\in\mathcal{M}_{1}^{+}(\mathcal{S}_{K})$ such that

Thus, conditions (3) and (4) concerning the mixture components are always satisfied for guaranteeing consistency and obtaining convergence rates of the Variational Bayes procedure.

When $K=1$, Lemma 3.2 does not apply as $\frac{2}{K}>1$. Nevertheless, as there is only one component, then any $p\in\mathcal{S}_{K}$ is equal to $1$ and the two conditions are immediately satisfied for any prior $\pi_{p}$ and any rate $r_{n,K}$ with $\rho_{p,n}=\pi_{p}$.

The central idea of the proof of Lemma 3.2 (given in details in Section 6) is to consider the ball $\mathcal{B}$ centered at $p^{0}$ of radius $Kr_{n,K}$ defined as:

Hence, when considering the restriction $\rho_{p,n}\in\mathcal{M}_{1}^{+}(\mathcal{S}_{K})$ of $\pi_{p}$ to $\mathcal{B}$, condition (3) is trivially satisfied and condition (4) is restricted to

This is a very classical assumption stated in many papers to study the concentration of the posterior . However, the computation of such a prior mass $\pi_{p}(\mathcal{B})$ is a major difficulty. Lemma 6.1 in treated the case of $L_{1}$-balls for Dirichlet priors. Since then, only a few papers in the literature addressed this issue. Our result extends the work in to KL-balls, which is of great interest in our study. Moreover, the range of Dirichlet priors for which Lemma 3.2 is applicable is the same as the one in .

We conclude Subsection 3.1 by a short discussion on the implementation of the VB approximation. Indeed, VB methods are meant to be practical objects, so there would be no point in proving the consistency of a VB approximation that would not be computable in practice. Many algorithms have been studied in the literature, with good performances – see and the references therein. In the case of mean-field approximation, the most popular method is to optimize iteratively with respect to all the independent components. Here this might seem difficult: it is indeed as difficult as maximizing the likelihood of a mixture. But a trick widely used in practice (see for example Section 7 in ) is to use the equality

This equality is once again a consequence of Lemma 2.1 (take $\textbf{E}=\{1,...,K\}$, $\lambda=(1/K,...,1/K)$ and $h(j)=\log(p_{j}q_{\theta_{j}}(X_{i}))$). This leads to the program:

This version can be solved by coordinate descent, see Algorithm 1. Update formulas once again follow from Lemma 2.1 (for instance, line 7 can be obtained with $\textbf{E}=\{1,...,K\}$, $\lambda=(1/K,...,1/K)$ and $h(j)=\int\log(p_{j})\rho_{p}(dp)+\int\log(q_{\theta_{j}}(X_{i}))\rho_{j}(d\theta_{j})$, more details are provided in Section 6). This algorithm is, in the case $\alpha=1$, exactly equivalent to the popular CAVI algorithm , where the $\omega^{i}_{j}$’s are interpreted as the posterior means of the latent variables $Z^{i}_{j}$’s. A very short numerical study is provided in the Supplementary Material but note that CAVI has already been extensively tested in practice .

2 Application to multinomial mixture models

The following corollary of Theorem 3.1 states that convergence of the VB approximation for the multinomial mixture model is achieved at rate $\frac{KV\log(nV)}{n}$ as soon as $V^{V}\geq K$, which is the case in many text mining models such as Latent Dirichlet Allocation for which the size of the vocabulary is very large:

The proof is in Section 6. We also specialize Algorithm 1 to the present setting (see Algorithm 2). Here $\psi$ denotes the Digamma function, $\psi(x)=\frac{{\rm d}}{{\rm d}x}\log[\Gamma(x)]$ where $\Gamma$ stands for the Gamma function $\Gamma(x)=\int_{0}^{\infty}\exp(-t)t^{x-1}{\rm d}t$.

3 Application to Gaussian mixture models

Let us now address the case of the Gaussian mixture model. This is one of the most popular mixture models for many applications including model based clustering and VB approximations have been studied in depth for this model . First, we will give rates of convergence of the VB approximation of the tempered posterior when the variance is known, and then when the variance is unknown.

Let us define $r_{n,K}=\frac{4\log(nK)}{n}\bigvee_{j=1}^{K}\frac{1}{n}\bigg[\frac{1}{2}\log\bigg(\frac{n}{2}\bigg)+\frac{V^{2}}{n\mathcal{V}^{2}}+\log\bigg(\frac{\mathcal{V}}{V}\bigg)+\frac{(\mu_{j}^{0})^{2}}{2\mathcal{V}^{2}}-\frac{1}{2}\bigg]$. Then, for any $\alpha\in(0,1)$,

One can see that for $n$ large enough, the convergence rate is $\frac{K\log(nK)}{n}$, which comes from the estimation of the weights of the mixture.

The Normal-Inverse-Gamma $\mathcal{NIG}(\mu,\theta^{2},a,b)$ is the distribution which density $w$ with respect to Lebesgue measure is defined by $w(x,y)=g(x|\mu,\frac{y}{\theta^{2}})h(y|a,b)$, where $g(.|\mu,\sigma^{2})$ is the density function of a Gaussian distribution of mean $\mu$ and variance $\sigma^{2}$, and $h(.|a,b)$ is the density distribution of an Inverse-Gamma of parameters $a$ and $b$.

For a Normal-Inverse-Gamma prior $\pi_{j}=\mathcal{NIG}(0,\mathcal{V}^{-2},1,\gamma^{2})$ for each $j=1,...,K$. With $r_{n,K}=\frac{4\log(nK)}{n}\bigvee_{j=1}^{K}\frac{1}{n}\bigg[2\log(n\sqrt{\mathcal{V}})+\frac{1}{2n\mathcal{V}^{2}}+\frac{(\mu_{j}^{0})^{2}}{2(\sigma_{j}^{0})^{2}\mathcal{V}^{2}}+\log\bigg(\frac{(\sigma_{j}^{0})^{2}}{\gamma^{2}}\bigg)+\frac{\gamma^{2}}{(\sigma_{j}^{0})^{2}}-\frac{1}{2}\log(2\pi)\bigg]$,

For the factorized prior $\pi_{j}=\mathcal{N}(0,\mathcal{V}^{2})\bigotimes\mathcal{IG}(1,\gamma^{2})$ for each $j=1,...,K$. With $r_{n,K}=\frac{4\log(nK)}{n}\bigvee_{j=1}^{K}\frac{1}{n}\bigg[2\log(n\sqrt{\mathcal{V}})+\frac{(\sigma_{j}^{0})^{2}}{2n\mathcal{V}^{2}}+\frac{(\mu_{j}^{0})^{2}}{2\mathcal{V}^{2}}+\frac{1}{2}\log\bigg(\frac{(\sigma_{j}^{0})^{2}}{\gamma^{4}}\bigg)+\frac{\gamma^{2}}{(\sigma_{j}^{0})^{2}}-\frac{1}{2}\log(2\pi)\bigg]$,

One can see that even when the variance has to be estimated, the convergence rate still achieves $\frac{K\log(nK)}{n}$ for $n$ large enough, whatever the form of the prior - factorized or not.

We give in Algorithm 3 a version of Algorithm 1 for unit-variance Gaussian mixtures with priors $\pi_{p}=\mathcal{D}_{K}(\alpha_{1},...,\alpha_{K})$ and $\pi_{j}=\mathcal{N}(0,\mathcal{V}^{2})$ where $\frac{2}{K}\leq\alpha_{j}\leq 1$ for $j=1,...,K$ and $\mathcal{V}^{2}>0$.

4 Extension to the misspecified case

From now we do not assume any longer that the true distribution $P^{0}$ belongs to the $K$-mixtures model. We still consider a prior $\pi=\pi_{p}\bigotimes_{j=1}^{K}\pi_{j}$ on $\theta\in\Theta_{K}$ for which $\pi_{p}\in\mathcal{M}_{1}^{+}(\mathcal{S}_{K})$ and $\pi_{j}\in\mathcal{M}_{1}^{+}(\Theta)$ for $j=1,...,K$.

For some value $r_{n,K}$, we introduce the set $\Theta_{K}(r_{n,K})$ of parameters $\theta^{*}\in\Theta_{K}$ such that:

there exists a set $\mathcal{A}_{n,K}\subset\mathcal{S}_{K}$ satisfying:

for each $p\in\mathcal{A}_{n,K}$, for each $j=1,...,K$, $\log\big(\frac{p_{j}^{*}}{p_{j}}\big)\leq Kr_{n,K}$,

$\pi_{p}(\mathcal{A}_{n,K})\geq e^{-nKr_{n,K}}$.

there are distributions $\rho_{j,n}\in\mathcal{M}_{1}^{+}(\Theta)$ ($j=1,...,K$) such that for $j=1,...,K$:

Let us discuss this definition. To begin with, the first item of the definition of $\Theta_{K}(r_{n,K})$ can seem quite restrictive. It is even a much more stronger assumption than (3) and (4). Nevertheless, the way to find the required measures $\rho_{p,n}$ in Lemma 3.2 in the well-specified case implies constructing in the proof such sets $\mathcal{A}_{n,K}$ for the true weight parameter $p^{0}$. As a consequence, it might seem reasonable to replace conditions (3) and (4) by the first part of the definition of $\Theta_{K}(r_{n,K})$. On the other hand, the condition given by (5) looks like those of Theorem 2.7 in . Once again, the difference is that inequalities must be satisfied here for each component. A condition on both the true distribution $P^{0}$ and the parameter $\theta^{*}$ considered is required through the expectation term. Besides, condition (5) is equivalent to (3) and (4) when the model is well-specified.

If there is no $r_{n,K}$ such that $\Theta_{K}(r_{n,K})$ is not empty, then the right-hand side is equal to infinity (by convention) for any value of $r_{n,K}$ and the inequality is useless. Nevertheless, this is not the case in models used in practice. We show an example below.

It is worth mentioning that even if this is not exactly an oracle inequality as the risk function in the left-hand side ($\alpha$-Renyi divergence) is lower than the right-hand side one (Kullback-Leibler divergence), but the theorem still remains of great interest. Indeed, when the minimizer of $\mathcal{K}(P^{0},P_{\theta})$ with respect to $\theta\in\Theta_{K}(r_{n,K})$ exists and is such that the corresponding Kullback-Leibler divergence is small, then our oracle inequality is informative as it gives a small bound on the expected risk of the Variational Bayes.

Variational Bayes model selection

with $\Theta_{K}=\mathcal{S}_{K}\times\Theta^{K}$, $\mathcal{S}_{K}=\{p_{K}=(p_{1,K},...,p_{K,K})\in^{K}/\sum_{j=1}^{K}p_{j,K}=1\}$ and the general notation $\theta_{K}=(p_{K},\theta_{1,K},...,\theta_{K,K})$. We would like to emphasize that the notations are slightly different as the size of each component parameter depends on the model complexity $K$. The entire parameter space $\Omega$ is the union of all parameter spaces $\Theta_{K}$ associated with each model index $K$: $\Omega=\cup_{K=1}^{\infty}\Theta_{K}$, and we can think of a whole statistical model $\mathcal{M}=\cup_{K=1}^{\infty}\mathcal{M}_{K}$ as the union of all collections $\mathcal{M}_{K}$. First, we can notice that different models $\mathcal{M}_{K}$ never overlap as parameters in each one do not have the same length. Nonetheless, parameters in complex models (models $\mathcal{M}_{K}$ with large $K$) can be sparse and therefore contain the "same information" as parameters in less complex ones, i.e. can lead to the same distribution $P_{\theta}$.

The prior specification is a crucial point. As mentioned above, each parameter depends on the number of components. Then, we specify a prior weight $\pi_{K}$ assigned to the model $\mathcal{M}_{K}$ and a conditional prior $\Pi_{K}(.)$ on $\theta_{K}\in\Theta_{K}$ given model $\mathcal{M}_{K}$. More precisely, we define our conditional prior on $\theta_{K}=(p_{K},\theta_{1,K},...,\theta_{K,K})$ as follows: given $K$, the weight parameter $p_{K}=(p_{1,K},...,p_{K,K})$ is supposed to follow a distribution $\pi_{p,K}$ on $\mathcal{M}_{1}^{+}(\mathcal{S}_{K})$; finally, given $K$, we set independent priors $\pi_{j,K}$ for the component parameters $\theta_{j,K}$ where each $\pi_{j,K}$ is a probability distribution on $\mathcal{M}_{1}^{+}(\Theta)$. In a nutshell:

We have to adapt the notations for the VB approximations. The tempered posteriors $\pi_{n,\alpha}^{K}(.|X_{1}^{n})$ on parameter $\theta_{K}\in\Theta_{K}$ given model $\mathcal{M}_{K}$, is defined again as

We recall that an alternative way to define the variational estimate is to use the Evidence Lower Bound via the optimization program (1):

where the function inside the argmax operator is the ELBO $\mathcal{L}(\rho_{K})$. For simplicity, we will just call ELBO $\mathcal{L}(K)$ the closest approximation to the log-evidence, i.e. the value of the lower bound evaluated in its maximum:

which is a penalized version of the ELBO. Note that taking $(\pi_{K})$ as uniform on a finite set $\{1,2,\dots,K_{\max}\}$ leads to the procedure described in . We discuss below the choice $\pi_{K}=2^{-K}$.

This oracle inequality shows that our variational distribution adaptively satisfies the best possible balance between bias (misspecification error) and variance (estimation error). If we assume that there is actually a $K_{0}$ and $\theta^{*}\in\Theta_{K_{0}}$ such that $P^{0}=P_{\theta^{*}}$ then the theorem will imply

The variance term is composed of two parts. The first one, $Kr_{n,K}$ up to a multiplicative constant, corresponds to the rate obtained when approximating the true distribution with mixtures of model $\mathcal{M}_{K}$. The second part of the overall rate can be interpreted as a complexity term over the different models reflecting our prior belief. For instance, if we want to penalize more complex models, we can take $\pi_{K}=2^{-K}$ and the corresponding term will be of order $K/n$. In practice, as soon as $\frac{1}{n}\lesssim r_{n,K}$, then this penalty term is negligible when compared to the approximating rate $Kr_{n,K}$: this means that this choice can be considered safe, as it does not interfere with the estimation rate.

Conclusion

Using variational inference, we studied consistency of variational approximations for estimation and model selection in mixtures. When considering tempered posteriors, we showed that Variational Bayes is consistent and we gave statistical guarantees to model selection based on the ELBO. For further investigation, it would be interesting to explore the case of Bayesian posteriors when $\alpha=1$. The recent work of Zhang and Gao gives the tools for tackling such an issue, and allows one to consider risk functions different from $\alpha$-Renyi divergence. But the conditions would be more stringent, and misspecification would be more problematic in this case.

Another point of interest is the study of the non-convex optimization program (2). Indeed, the proposed coordinate optimization can lead to a local extremum, and this implies that one needs to pay attention to initialization. The same problem also occurs in the Expectation-Maximization (EM) algorithm. In practice, users often run EM or CAVI several times with different initial distributions. Many practical ideas were proposed to target the global extrema more efficiently with EM and could be extended to CAVI. But the question of convergence remains open in theory.

Finally, note that our results are remarkable as there are almost no conditions on the mixtures considered. In this paper we have focused on estimating the true probability distribution $P^{0}$, even in the well-specified case. We have no results on the estimation of the parameters. In the case of mixtures, these results are extremely difficult to obtain even for Gaussian mixtures . They require restrictions on the parameters set and lead to different rates of convergence. The consistency of VB for the estimation of the parameters remains open.

Acknowledgements

We thank the Associate Editor and the anonymous Referee for their insightful comments on the paper.

References

Proofs

We provide in this section two useful lemmas required in many proofs below.

The lemma below was first stated by for mixtures of Gaussians, checked that the proof remains valid for general mixtures. It is a tool widely used in signal processing . We provide the proof for the sake of completeness.

Let $p,p^{0}\in\mathcal{S}_{K}$ and $\theta_{j},\theta_{j}^{0}\in\Theta$ for $j=1,...,K$. Then,

For any nonnegative numbers $\alpha_{1},...,\alpha_{K}$ and positive $\beta_{1},...,\beta_{K}$, we have:

$\begin{aligned} \left(\sum_{j=1}^{K}\alpha_{j}\right)\log\left(\frac{\sum_{j=1}^{K}\alpha_{j}}{\sum_{j=1}^{K}\beta_{j}}\right)&=\left(\sum_{j=1}^{K}\beta_{j}\right)\left(\frac{\sum_{j=1}^{K}\alpha_{j}}{\sum_{j=1}^{K}\beta_{j}}\right)\log\left(\frac{\sum_{j=1}^{K}\alpha_{j}}{\sum_{j=1}^{K}\beta_{j}}\right)\\ &=\left(\sum_{j=1}^{K}\beta_{j}\right)\left(\sum_{j=1}^{K}\frac{\beta_{j}}{\sum_{l=1}^{K}\beta_{l}}\frac{\alpha_{j}}{\beta_{j}}\right)\log\left(\sum_{j=1}^{K}\frac{\beta_{j}}{\sum_{l=1}^{K}\beta_{l}}\frac{\alpha_{j}}{\beta_{j}}\right)\\ &=\left(\sum_{j=1}^{K}\beta_{j}\right)f\left(\sum_{j=1}^{K}\frac{\beta_{j}}{\sum_{l=1}^{K}\beta_{l}}\frac{\alpha_{j}}{\beta_{j}}\right)\end{aligned}$

where $f$ is the convex function $x\longmapsto x\log(x)$. As $\sum_{j=1}^{K}\frac{\beta_{j}}{\sum_{l=1}^{K}\beta_{l}}=1$, then using Jensen’s inequality:

$\begin{aligned} \left(\sum_{j=1}^{K}\alpha_{j}\right)\log\left(\frac{\sum_{j=1}^{K}\alpha_{j}}{\sum_{j=1}^{K}\beta_{j}}\right)&=\left(\sum_{j=1}^{K}\beta_{j}\right)f\left(\sum_{j=1}^{K}\frac{\beta_{j}}{\sum_{l=1}^{K}\beta_{l}}\frac{\alpha_{j}}{\beta_{j}}\right)\\ &\leq\left(\sum_{j=1}^{K}\beta_{j}\right)\sum_{j=1}^{K}\frac{\beta_{j}}{\sum_{l=1}^{K}\beta_{l}}f\left(\frac{\alpha_{j}}{\beta_{j}}\right)\\ &=\left(\sum_{j=1}^{K}\beta_{j}\right)\sum_{j=1}^{K}\frac{\beta_{j}}{\sum_{l=1}^{K}\beta_{l}}\frac{\alpha_{j}}{\beta_{j}}\log\left(\frac{\alpha_{j}}{\beta_{j}}\right)\\ &=\sum_{j=1}^{K}\alpha_{j}\log\left(\frac{\alpha_{j}}{\beta_{j}}\right).\end{aligned}$

The inequality remains valid when some or all $\beta_{j}$’s are zero. Indeed, assume that $\beta_{j}=0$. If $\alpha_{j}\neq 0$, then the $j^{th}$ term of the sum in the right-hand side is $\alpha_{j}\log({\alpha_{j}}/{\beta_{j}})=+\infty$, and the result is obvious. Otherwise, $\alpha_{j}=0$, hence the $j^{th}$ term of each sum in the inequality is zero as $\alpha_{j}\log({\alpha_{j}}/{\beta_{j}})=0$, and the inequality can be obtained considering only the other numbers.

Thus, for $p,p^{0}\in\mathcal{S}_{K}$ and $\theta_{j},\theta_{j}^{0}\in\Theta$ for $j=1,...,K$:

$\begin{aligned} \mathcal{K}\left(\sum_{j=1}^{K}p_{j}^{0}q_{\theta_{j}^{0}},\sum_{j=1}^{K}p_{j}q_{\theta_{j}}\right)&=\int\left(\sum_{j=1}^{K}p_{j}^{0}q_{\theta_{j}^{0}}\right)\log\left(\frac{\sum_{j=1}^{K}p_{j}^{0}q_{\theta_{j}^{0}}}{\sum_{j=1}^{K}p_{j}q_{\theta_{j}}}\right)\\ &\leq\int\sum_{j=1}^{K}p_{j}^{0}q_{\theta_{j}^{0}}\log\left(\frac{p_{j}^{0}q_{\theta_{j}^{0}}}{p_{j}q_{\theta_{j}}}\right)\\ &=\int\sum_{j=1}^{K}p_{j}^{0}q_{\theta_{j}^{0}}\log\left(\frac{p_{j}^{0}}{p_{j}}\right)+\int\sum_{j=1}^{K}p_{j}^{0}q_{\theta_{j}^{0}}\log\left(\frac{q_{\theta_{j}^{0}}}{q_{\theta_{j}}}\right)\\ &=\sum_{j=1}^{K}p_{j}^{0}\log\left(\frac{p_{j}^{0}}{p_{j}}\right)\left(\int q_{\theta_{j}^{0}}\right)+\sum_{j=1}^{K}p_{j}^{0}\int q_{\theta_{j}^{0}}\log\left(\frac{q_{\theta_{j}^{0}}}{q_{\theta_{j}}}\right)\\ &=\mathcal{K}(p^{0},p)+\sum_{j=1}^{K}p_{j}^{0}\mathcal{K}(q_{\theta_{j}^{0}},q_{\theta_{j}}),\end{aligned}$

1.2 KL-divergence between Gaussian distributions and between Normal-Inverse-Gamma distributions

We give in this section the Kullback-Leibler divergence between 1-dimensional Gaussian distributions and between Normal-Inverse-Gamma distributions.

We denote $u$ and $v$ the density functions of the respective Gaussian distributions $\mathcal{N}(\mu_{u},\sigma^{2}_{u})$ and $\mathcal{N}(\mu_{v},\sigma^{2}_{v})$. Similarly, we denote $p$ and $q$ the two densities of $\mathcal{NIG}(\mu_{1},\theta_{1}^{2},a_{1},b_{1})$ and $\mathcal{NIG}(\mu_{2},\theta_{2}^{2},a_{2},b_{2})$. Then:

$\begin{aligned} \mathcal{K}(p,q)=\frac{1}{2}\log\bigg(\frac{\theta_{1}^{2}}{\theta_{2}^{2}}\bigg)&+\frac{\theta_{2}^{2}}{2\theta_{1}^{2}}+\frac{\theta_{2}^{2}(\mu_{2}-\mu_{1})^{2}}{2}\frac{a_{1}}{b_{1}}-\frac{1}{2}\\ &+(a_{1}-a_{2})\psi(a_{1})+\log\bigg(\frac{\Gamma(a_{2})}{\Gamma(a_{1})}\bigg)+a_{2}\log\bigg(\frac{b_{1}}{b_{2}}\bigg)+a_{1}\frac{b_{2}-b_{1}}{b_{1}}.\end{aligned}$

The first equality is extremely classical so we don’t provide the proof. For the second one,

Using the KL-divergence between Gaussians:

and using the KL-divergence between Inverse-Gamma distributions

where $\Gamma$ and $\psi$ are respectively the Gamma and Digamma functions, we have:

$\begin{aligned} \mathcal{K}(p,q)=\frac{1}{2}\log\bigg(\frac{\theta_{1}^{2}}{\theta_{2}^{2}}\bigg)&+\frac{\theta_{2}^{2}}{2\theta_{1}^{2}}+\frac{\theta_{2}^{2}(\mu_{2}-\mu_{1})^{2}}{2}\frac{a_{1}}{b_{1}}-\frac{1}{2}\\ &+(a_{1}-a_{2})\psi(a_{1})+\log\bigg(\frac{\Gamma(a_{2})}{\Gamma(a_{1})}\bigg)+a_{2}\log\bigg(\frac{b_{1}}{b_{2}}\bigg)+a_{1}\frac{b_{2}-b_{1}}{b_{1}}.\end{aligned}$

2 Proof of Theorem 3.1

This result relies on an application of Theorem 2.6 in to mixture models. The proof of Theorem 2.6 in itself relies mostly on a deviation inequality from and on PAC-Bayesian theory .

Fix $0<\alpha<1$. Theorem 2.6 from gives:

the last inequality being obtained thanks to Theorem 28 in . Gathering all the pieces together leads to

that is the result stated in Theorem 3.1. ∎

3 Proof of Lemma 3.2

Let us define $\rho_{p,n}\in\mathcal{M}_{1}^{+}(\mathcal{S}_{K})$ by the following formula $\rho_{p,n}(dp)\propto\mathbf{1}(p\in\mathcal{B})\pi_{p}(dp)$ with

where $A=\frac{2}{K}$ and $M_{p}^{0}=\max\{p_{j}^{0}/j=1,...,K\}$. We adopt the notation $S=\sum_{j=1}^{K}\alpha_{j}$ in the following. Recall that by assumption $K\geq 2$ and hence $A=\frac{2}{K}\leq 1$.

First, $\int\mathcal{K}(p^{0},p)\rho_{p,n}(dp)\leq Kr_{n,K}^{\prime}$.

Then, let us show that $\mathcal{K}(\rho_{p,n},\pi_{p})\leq Knr_{n,K}^{\prime}$. For that, let us define

where $K$ is such that $p_{K}^{0}=\max\{p_{j}^{0}/j=1,...,K\}$ (this assumption can always be fulfilled by reordering and relabelling the vector components). Then, $p_{K}^{0}\geq\frac{1}{K}$ (otherwise, the sum of the components of $p^{0}$ would be strictly lower than $1$ and the vector would not be included in $\mathcal{S}_{K}$). We will show that $\mathcal{A}\subset\mathcal{B}$ and that $\pi_{p}(\mathcal{A})\geq e^{-Knr_{n,K}^{\prime}}$. Then, we will conclude thanks to the following formula: $\mathcal{K}(\rho_{p,n},\pi_{p})=-\log(\pi_{p}(\mathcal{B}))$.

First, let us show that $\mathcal{A}\subset\mathcal{B}$.

Let $p\in\mathcal{A}$. As $p_{K}=1-\sum_{j=1}^{K-1}p_{j}$, we just need to check that $\mathcal{K}(p^{0},p)\leq Kr_{n,K}^{\prime}$ and that $p_{j}\geq 0$ for each $j=1,...,K$.

The first part can be proven using the definition of $\mathcal{A}$. According to the $K-1$ left-hand side inequalities in the definition of $\mathcal{A}$,

$\begin{aligned} \mathcal{K}(p^{0},p)=\sum_{j=1}^{K-1}p_{j}^{0}\log\bigg(\frac{p_{j}^{0}}{p_{j}}\bigg)+p_{K}^{0}\log\bigg(\frac{p_{K}^{0}}{p_{K}}\bigg)&\leq\sum_{j=1}^{K-1}p_{j}^{0}\log(e^{Kr_{n,K}^{\prime}})+p_{K}^{0}\log\bigg(\frac{p_{K}^{0}}{p_{K}}\bigg)\\ &=\sum_{j=1}^{K-1}p_{j}^{0}Kr_{n,K}^{\prime}+p_{K}^{0}\log\left(\frac{p_{K}^{0}}{p_{K}}\right)\\ &=(1-p_{K}^{0})Kr_{n,K}^{\prime}+p_{K}^{0}\log\left(\frac{p_{K}^{0}}{p_{K}}\right).\end{aligned}$

All we need to show now is that $\log\left(\frac{p_{K}^{0}}{p_{K}}\right)\leq Kr_{n,K}^{\prime}$. This comes from the following inequalities:

$\begin{aligned} \log\bigg(\frac{p_{K}^{0}}{p_{K}}\bigg)=\log\bigg(\frac{p_{K}^{0}}{1-\sum_{j=1}^{K-1}p_{j}}\bigg)&\leq\log\bigg(\frac{p_{K}^{0}}{1-\sum_{j=1}^{K-1}p_{j}^{0}e^{-Kr_{n,K}^{\prime}}-\frac{p_{K}^{0}}{n}}\bigg)\\ &=\log\bigg(\frac{p_{K}^{0}}{1-(1-p_{K}^{0})e^{-Kr_{n,K}^{\prime}}-\frac{p_{K}^{0}}{n}}\bigg)\\ &\leq\frac{p_{K}^{0}}{1-(1-p_{K}^{0})e^{-Kr_{n,K}^{\prime}}-\frac{p_{K}^{0}}{n}}-1\\ &=\frac{p_{K}^{0}-1+(1-p_{K}^{0})e^{-Kr_{n,K}^{\prime}}+\frac{p_{K}^{0}}{n}}{1-(1-p_{K}^{0})e^{-Kr_{n,K}^{\prime}}-\frac{p_{K}^{0}}{n}}\end{aligned}$

$\begin{aligned} \log\bigg(\frac{p_{K}^{0}}{p_{K}}\bigg)&\leq\frac{p_{K}^{0}-1+(1-p_{K}^{0})e^{-Kr_{n,K}^{\prime}}+\frac{p_{K}^{0}}{n}}{1-(1-p_{K}^{0})e^{-Kr_{n,K}^{\prime}}-\frac{p_{K}^{0}}{n}}&&\\ &=\frac{\frac{p_{K}^{0}}{n}-(1-p_{K}^{0})(1-e^{-Kr_{n,K}^{\prime}})}{p_{K}^{0}(1-\frac{1}{n})+(1-p_{K}^{0})(1-e^{-Kr_{n,K}^{\prime}})}&&\\ &=\frac{\frac{1}{n}-(\frac{1}{p_{K}^{0}}-1)(1-e^{-Kr_{n,K}^{\prime}})}{(1-\frac{1}{n})+(\frac{1}{p_{K}^{0}}-1)(1-e^{-Kr_{n,K}^{\prime}})}&&\\ &\leq\frac{\frac{1}{n}}{1-\frac{1}{n}}=\frac{1}{n-1}&&\\ &\leq Kr_{n,K}^{\prime}.&&\end{aligned}$

Hence $\hskip 5.69046pt\mathcal{K}(p^{0},p)\leq(1-p_{K}^{0})Kr_{n,K}^{\prime}+p_{K}^{0}\log\left(\frac{p_{K}^{0}}{p_{K}}\right)\leq(1-p_{K}^{0})Kr_{n,K}^{\prime}+p_{K}^{0}Kr_{n,K}^{\prime}=Kr_{n,K}^{\prime}$.

On the other hand, for $j=1,...,K-1$, $p_{j}\geq p_{j}^{0}e^{-Kr_{n,K}^{\prime}}\geq 0$ and:

$\begin{aligned} p_{K}=1-\sum_{j=1}^{K-1}p_{j}&\geq 1-\sum_{j=1}^{K-1}\big(p_{j}^{0}e^{-Kr_{n,K}^{\prime}}+\frac{p^{0}_{K}}{n(K-1)}\big)\\ &=1-\bigg((1-p_{K}^{0})e^{-Kr_{n,K}^{\prime}}+\frac{p_{K}^{0}}{n}\bigg)\\ &\geq 1-(1-p_{K}^{0})e^{-Kr_{n,K}^{\prime}}-p_{K}^{0}\\ &=(1-p_{K}^{0})(1-e^{-Kr_{n,K}^{\prime}})\\ &\geq 0.\end{aligned}$

Then, $p\in\mathcal{B}$, and finally $\mathcal{A}\subset\mathcal{B}$.

Now, let us show that $\pi_{p}(\mathcal{A})\geq e^{-Knr_{n,K}^{\prime}}$.

Let us denote $f$ the density of the $\pi_{p}=\mathcal{D}_{K}(\alpha_{1},...,\alpha_{K})$ Dirichlet distribution:

Thus, we can lower bound $\pi_{p}(\mathcal{A})$:

$\begin{aligned} \pi_{p}(\mathcal{A})=\int_{\mathcal{A}}f(p_{1},...,p_{K})dp&=\int_{\mathcal{A}}\frac{\Gamma\big(S\big)}{\prod\limits_{j=1}^{K}\Gamma(\alpha_{j})}\prod_{j=1}^{K}p_{j}^{\alpha_{j}-1}\hskip 2.84544pt\mathbf{1}(p\in\mathcal{S}_{K})dp\\ &\geq\frac{\Gamma\big(S\big)}{\prod\limits_{j=1}^{K}\Gamma(\alpha_{j})}\prod_{j=1}^{K-1}\int_{p_{j}^{0}e^{-Kr_{n,K}^{\prime}}}^{p_{j}^{0}e^{-Kr_{n,K}^{\prime}}+\frac{p_{K}^{0}}{n(K-1)}}p_{j}^{\alpha_{j}-1}\hskip 2.84544ptdp_{j}\end{aligned}$

as for $p\in\mathcal{A}$, $0\leq p_{j}^{0}e^{-Kr_{n,K}^{\prime}}\leq p_{j}\leq p_{j}^{0}e^{-Kr_{n,K}^{\prime}}+\frac{p_{K}^{0}}{n(K-1)}\leq 1$ for each $j=1,...,K-1$ (as $\mathcal{A}\subset\mathcal{B}$), and then $p_{j}^{\alpha_{j}-1}\geq 1$.

Then, by definition of $r_{n,K}^{\prime}$, $\frac{p_{K}^{0}}{n(K-1)}\geq\Gamma(A)^{\frac{K}{K-1}}e^{-nr_{n,K}^{\prime}}$, and using inequalities $\Gamma(A)\geq\Gamma(\alpha_{j})$ as $A\leq\alpha_{j}\leq 1$ and $\Gamma(S)\geq 1$ as $S\geq 2$,

$\begin{aligned} \pi_{p}(\mathcal{A})&\geq\frac{\Gamma\big(S\big)}{\prod\limits_{j=1}^{K}\Gamma(\alpha_{j})}\prod_{j=1}^{K-1}\int_{p_{j}^{0}e^{-Kr_{n,K}^{\prime}}}^{p_{j}^{0}e^{-Kr_{n,K}^{\prime}}+\Gamma(A)^{\frac{K}{K-1}}e^{-nr_{n,K}^{\prime}}}p_{j}^{\alpha_{j}-1}\hskip 2.84544ptdp_{j}\\ &\geq\frac{\Gamma\big(S\big)}{\prod\limits_{j=1}^{K}\Gamma(\alpha_{j})}\prod_{j=1}^{K-1}\int_{p_{j}^{0}e^{-Kr_{n,K}^{\prime}}}^{p_{j}^{0}e^{-Kr_{n,K}^{\prime}}+\Gamma(A)^{\frac{K}{K-1}}e^{-nr_{n,K}^{\prime}}}\hskip 2.84544ptdp_{j}\\ &=\frac{\Gamma\big(S\big)}{\prod\limits_{j=1}^{K}\Gamma(\alpha_{j})}\prod_{j=1}^{K-1}\Gamma(A)^{\frac{K}{K-1}}e^{-nr_{n,K}^{\prime}}\\ &=\frac{\Gamma\big(S\big)}{\prod\limits_{j=1}^{K}\Gamma(\alpha_{j})}\Gamma(A)^{K}e^{-n(K-1)r_{n,K}^{\prime}}\\ &\geq e^{-nKr_{n,K}^{\prime}}.\end{aligned}$

Hence, as $\mathcal{A}\subset\mathcal{B}$, $\pi_{p}(\mathcal{B})\geq\pi_{p}(\mathcal{A})\geq e^{-nKr_{n,K}^{\prime}}$, and finally, $\mathcal{K}(\rho_{p,n},\pi_{p})=-\log(\pi_{p}(\mathcal{B}))\leq Knr_{n,K}^{\prime}$.

We just proved the lemma but with the rate $r_{n,K}^{\prime}$ instead of the value $r_{n,K}$ used in the lemma. We can conlude by noticing that the result is valid for every $r$ such that $r_{n,K}^{\prime}\leq r$, and that in particular $r_{n,K}^{\prime}\leq r_{n,K}$. This last result comes from the inequality:

which is a direct application of the left-hand side of inequality (3.2) part 3 in with $x=\frac{A}{2}>0$ and $\lambda=\frac{A}{2}\in(0,1)$. As, $1+\frac{A}{2}\in$, then $\Gamma(1+\frac{A}{2})\leq 1$, and $\frac{1}{\left(\frac{A}{2}\right)^{1-\frac{A}{2}}}=K^{1-\frac{A}{2}}\leq K$. Thus:

and as $K\geq 2$ and $p^{0}_{K}\geq\frac{1}{K}$, it follows that

i.e. $r_{n,K}^{\prime}\leq\max(\frac{1}{K(n-1)},\frac{\log(nK^{4})}{n})\leq\max(\frac{1}{K(n-1)},\frac{4\log(nK)}{n})$. Besides, $\frac{n}{n-1}=1+\frac{1}{n-1}\leq 2$ implies $\frac{1}{K(n-1)}\leq\frac{1}{2(n-1)}\leq\frac{1}{n}\leq\frac{4\log(2)}{n}\leq\frac{4\log(nK)}{n}$, and finally $r_{n,K}^{\prime}\leq\frac{4\log(nK)}{n}=r_{n,K}$.

4 Proof of Corollary 3.3

According to Lemma 3.2, there exists a distribution $\rho_{p,n}\in\mathcal{M}_{1}^{+}(\mathcal{S}_{K})$ such that

Similarly, the same result states that there exists distributions $\rho_{j,n}\in\mathcal{M}_{1}^{+}(\mathcal{S}_{V})$ for $j=1,...,K$ such that