Frequentist Consistency of Variational Bayes

TL;DR

This paper establishes the theoretical foundations of variational Bayes (VB), proving its consistency and asymptotic normality, and connects VB to classical frequentist estimators, thus enhancing its credibility for scalable Bayesian inference.

Contribution

We provide the first rigorous proof of frequentist consistency and asymptotic normality for VB methods, linking them to variational Bernstein-von Mises theorems and applying results to various Bayesian models.

Findings

VB posterior converges to the KL minimizer centered at the truth

Variational expectation of parameters is consistent and asymptotically normal

Theoretical results are supported by simulation studies

Abstract

A key challenge for modern Bayesian statistics is how to perform scalable inference of posterior distributions. To address this challenge, variational Bayes (VB) methods have emerged as a popular alternative to the classical Markov chain Monte Carlo (MCMC) methods. VB methods tend to be faster while achieving comparable predictive performance. However, there are few theoretical results around VB. In this paper, we establish frequentist consistency and asymptotic normality of VB methods. Specifically, we connect VB methods to point estimates based on variational approximations, called frequentist variational approximations, and we use the connection to prove a variational Bernstein-von Mises theorem. The theorem leverages the theoretical characterizations of frequentist variational approximations to understand asymptotic properties of VB. In summary, we prove that (1) the VB posterior converges to the Kullback-Leibler (KL) minimizer of a normal distribution, centered at the truth and (2) the corresponding variational expectation of the parameter is consistent and asymptotically normal. As applications of the theorem, we derive asymptotic properties of VB posteriors in Bayesian mixture models, Bayesian generalized linear mixed models, and Bayesian stochastic block models. We conduct a simulation study to illustrate these theoretical results.

Figures

Tables

Table 1. Table 1 : Glossary of terms

Name	Definition
Variational log likelihood	$M_n(\theta ;x):={sup}_{q\left(z\right)\in {𝒬}^n}\int q\left(z\right)\log \frac{p(x,z|\theta )}{q\left(z\right)}\text{d}z$
Variational frequentist estimate (vfe)	${\hat{\theta }}_n:={\arg \max }_{\theta }{M}_n(\theta ;x)$
vb ideal	${\pi }^{\ast }\left(\theta |x\right):=\frac{p\left(\theta \right)\exp \left\{M_n(\theta ;x)\right\}}{\int p\left(\theta \right)\exp \left\{M_n(\theta ;x)\right\}\text{d}\theta }$
Evidence Lower Bound (elbo)	elbo ($q(\theta ,z)):=\int \int q(\theta \left)q\right(z)\log \frac{p(x,z,\theta )}{q\left(\theta \right)q\left(z\right)}\text{d}\theta \text{d}z$
vb posterior	$q^{\ast }\left(\theta \right):={\arg \max }_{q\left(\theta \right)\in {𝒬}^d}{sup}_{q\left(z\right)\in {𝒬}^n}\text{<span>elbo</span>}\left(q(\theta ,z)\right)$
vb estimate (vbe)	${\hat{\theta }}_n^{\ast }:=\int \theta \cdot q^{\ast }\left(\theta \right)\text{d}\theta$

Equations

$$p(\theta ,x,z)=p(\theta )\prod _{i=1}^{n}p(z_i|\theta )p(x_i|z_i,\theta ).$$

$${\mathcal{Q}}^{n+d}=\left\{q:q(\theta ,z)={\prod }_{i=1}^d{q}_{{\theta }_i}({\theta }_i){\prod }_{j=1}^n{q}_{z_j}(z_j)\right\}.$$

$$q^{\ast }(\theta ,z)=\underset{q(\theta ,z)\in {\mathcal{Q}}^{n+d}}{\mathrm{argmin}}\text{lx@glossaries@gls@link}acronymKL\text{leavevmode}\text{textsc}kl(q(\theta ,z)||p(\theta ,z\text{leavevmode}|\text{leavevmode}x)).$$

$$\text{lx@glossaries@gls@link}acronymELBO\text{leavevmode}\text{textsc}elbo(q(\theta ,z))=-\int q(\theta ,z)\log \frac{q(\theta ,z)}{p(\theta ,x,z)}\text{d}\theta \text{d}z.$$

$${\hat{\theta }}_n^{\ast }:=\int \theta \cdot q^{\ast }(\theta )\text{d}\theta .$$

$$M_n(\theta ;x)=\underset{q(z)}{\max }{\mathbb{E}}_{q(z)}\left[\log p(x,z|\theta )-\log q(z)\right].$$

$$q^{†}(z)=\arg \underset{q}{\min }\text{KL}(q(z)||p(z|x,\theta )).$$

$${\hat{\theta }}_n=\arg \underset{\theta }{\max }{M}_n(\theta ;x).$$

$${\pi }^{\ast }(\theta |x)=\frac{p(\theta )\exp \{M_n(\theta ;x)\}}{\int p(\theta )\exp \{M_n(\theta ;x)\}\text{d}\theta }.$$

$${\pi }^{\ast }(\theta \text{leavevmode}|\text{leavevmode}x)=\frac{p(\theta )\exp (M_n(\theta ;x))}{\int p(\theta )\exp (M_n(\theta ;x))\mathrm{d}\theta },$$

$$\ell (\theta ;x)\propto \exp (M_n(\theta ;x)).$$

$$\int {\varphi }_n(x)p(x,z\text{leavevmode}|\text{leavevmode}{\theta }_0)\text{d}z\text{d}x\to 0$$

$$\underset{\theta :||\theta -{\theta }_0||\ge ϵ}{\sup }\int (1-{\varphi }_n(x))\frac{\ell (\theta ;x)}{\ell ({\theta }_0;x)}p(x,z\text{leavevmode}|\text{leavevmode}{\theta }_0)\text{d}z\text{d}x\to 0,$$

$$\underset{h\in K}{\sup }|M_n(\theta +{\delta }_{n}h;x)-M_n(\theta ;x)-h^{\top }{V}_{{\theta }_0}{\Delta }_{n,{\theta }_0}+\frac{1}{2}{h}^{\top }{V}_{{\theta }_0}h|\stackrel{P_{{\theta }_0}}{\to }0,$$

$$\tilde{\theta }={\delta }_n^{-1}(\theta -{\theta }_0).$$

$${\pi }_{\tilde{\theta }}^{\ast }(\tilde{\theta }\text{leavevmode}|\text{leavevmode}x)={\pi }^{\ast }({\theta }_0+{\delta }_n\tilde{\theta }\text{leavevmode}|\text{leavevmode}x)\cdot |\det ({\delta }_n)|.$$

$$||{\pi }_{\tilde{\theta }}^{\ast }(\cdot \text{leavevmode}|\text{leavevmode}x)-\mathcal{N}(\cdot ;{\Delta }_{n,{\theta }_0},V_{{\theta }_0}^{-1})|{|}_{\text{lx@glossaries@gls@link}acronymTV\text{leavevmode}\text{textsc}tv}\stackrel{P_{{\theta }_0}}{\to }0.$$

$${\int }_{|\tilde{\theta }|>M_n}{\pi }_{\tilde{\theta }}^{\ast }(\tilde{\theta }\text{leavevmode}|\text{leavevmode}x)\mathrm{d}\tilde{\theta }\stackrel{P_{{\theta }_0}}{\to }0$$

$$\underset{q(\theta )\in {\mathcal{Q}}^d}{\mathrm{argmin}}\text{lx@glossaries@gls@link}acronymKL\text{leavevmode}\text{textsc}kl(q(\theta )||{\pi }^{\ast }(\theta \text{leavevmode}|\text{leavevmode}x))\stackrel{d}{\to }{\delta }_{{\theta }_0}\text{in Pθ0-probability}.$$

$${\Vert \underset{q\in {\mathcal{Q}}^d}{\mathrm{argmin}}\text{lx@glossaries@gls@link}acronymKL\text{leavevmode}\text{textsc}kl(q(\cdot )||{\pi }_{\tilde{\theta }}^{\ast }(\cdot \text{leavevmode}|\text{leavevmode}x))-\underset{q\in {\mathcal{Q}}^d}{\mathrm{argmin}}\text{lx@glossaries@gls@link}acronymKL\text{leavevmode}\text{textsc}kl(q(\cdot )||\mathcal{N}(\cdot ;{\Delta }_{n,{\theta }_0},V_{{\theta }_0}^{-1}))\Vert }_{\text{lx@glossaries@gls@link}acronymTV\text{leavevmode}\text{textsc}tv}\stackrel{P_{{\theta }_0}}{\to }0.$$

$$q(\theta ,z)=q(\theta )q(z),$$

$$\text{lx@glossaries@gls@link}acronymELBO\text{leavevmode}\text{textsc}elbo(q(\theta )q(z))$$

$$=\int \int q(\theta )q(z)\log \frac{p(\theta )p(x,z\text{leavevmode}|\text{leavevmode}\theta )}{q(\theta )q(z)}\mathrm{d}\theta \mathrm{d}z$$

$$=\int q(\theta )\log \frac{p(\theta )}{q(\theta )}\mathrm{d}\theta +\int q(\theta )\int q(z)\log \frac{p(x,z|\theta )}{q(z)}\mathrm{d}\theta \mathrm{d}z.$$

$$q^{\ast }(\theta )=\underset{q(\theta )}{\mathrm{argmax}}\underset{q(z)}{\sup }\int q(\theta )\left(\log \left[p(\theta )\exp \left\{\int q(z)\log \frac{p(x,z\text{leavevmode}|\text{leavevmode}\theta )}{q(z)}\mathrm{d}z\right\}\right]-\log q(\theta )\right)\mathrm{d}\theta .$$

$$\text{lx@glossaries@gls@link}acronymELBO{\text{leavevmode}\text{textsc}elbo}_p(q(\theta )):=\underset{q(z)}{\sup }\int q(\theta )\left(\log \left[p(\theta )\exp \left\{\int q(z)\log \frac{p(x,z\text{leavevmode}|\text{leavevmode}\theta )}{q(z)}\mathrm{d}z\right\}\right]-\log q(\theta )\right)\mathrm{d}\theta .$$

$$\text{lx@glossaries@gls@link}acronymELBO{\text{leavevmode}\text{textsc}elbo}_p(q(\theta ))=-\text{lx@glossaries@gls@link}acronymKL\text{leavevmode}\text{textsc}kl(q(\theta )||{\pi }^{\ast }(\theta \text{leavevmode}|\text{leavevmode}x))+o_P(1).$$

$$-\text{lx@glossaries@gls@link}acronymKL\text{leavevmode}\text{textsc}kl(q(\theta )||{\pi }^{\ast }(\theta \text{leavevmode}|\text{leavevmode}x))$$

$$=$$

$$q^{\ast }(\theta )\stackrel{d}{\to }{\delta }_{{\theta }_0}\text{in Pθ0-probability}.$$

$${\Vert q_{\tilde{\theta }}^{\ast }(\cdot )-\underset{q\in {\mathcal{Q}}^d}{\mathrm{argmin}}\text{lx@glossaries@gls@link}acronymKL\text{leavevmode}\text{textsc}kl(q(\cdot )||\mathcal{N}(\cdot ;{\Delta }_{n,{\theta }_0},V_{{\theta }_0}^{-1}))\Vert }_{\text{lx@glossaries@gls@link}acronymTV\text{leavevmode}\text{textsc}tv}\stackrel{P_{{\theta }_0}}{\to }0.$$

Full text

Keywords: Bernstein–von Mises, Bayesian inference, variational methods, consistency, asymptotic normality, statistical computing

1 Introduction

Bayesian modeling is a powerful approach for discovering hidden patterns in data. We begin by setting up a probability model of latent variables and observations. We incorporate prior knowledge by setting priors on latent variables and a functional form of the likelihood. Finally we infer the posterior, the conditional distribution of the latent variables given the observations.

For many modern Bayesian models, exact computation of the posterior is intractable and statisticians must resort to approximate posterior inference. For decades, Markov chain Monte Carlo (mcmc) sampling (Hastings:1970; gelfand1990sampling; Robert:2004) has maintained its status as the dominant approach to this problem. mcmc algorithms are easy to use and theoretically sound. In recent years, however, data sizes have soared. This challenges mcmc methods, for which convergence can be slow, and calls upon scalable alternatives. One popular class of alternatives is variational Bayes (vb) methods.

To describe vb, we introduce notation for the posterior inference problem. Consider observations $x=x_{1:n}$. We posit local latent variables $z=z_{1:n}$, one per observation, and global latent variables $\theta=\theta_{1:d}$. This gives a joint,

[TABLE]

The posterior inference problem is to calculate the posterior $p(\theta,z\,|\,x)$.

This division of latent variables is common in modern Bayesian statistics.111In particular, our results are applicable to general models with local and global latent variables (Hoffman:2013). The number of local variables $z$ increases with the sample size $n$; the number of global variables $\theta$ does not. We also note that the conditional independence of Equation 1 is not necessary for our results. But we use this common setup to simplify the presentation. In the Bayesian Gaussian mixture model (gmm) (roberts1998bayesian), the component means, covariances, and mixture proportions are global latent variables; the mixture assignments of each observation are local latent variables. In the Bayesian generalized linear mixed model (glmm) (breslow1993approximate), the intercept and slope are global latent variables; the group-specific random effects are local latent variables. In the Bayesian stochastic block model (sbm) (Wiggins:2008), the cluster assignment probabilities and edge probabilities matrix are two sets of global latent variables; the node-specific cluster assignments are local latent variables. In the latent Dirichlet allocation (lda) model  (blei2003latent), the topic-specific word distributions are global latent variables; the document-specific topic distributions are local latent variables. We will study all these examples below.

vb methods formulate posterior inference as an optimization (jordan1999introduction; wainwright2008graphical; Blei:2016). We consider a family of distributions of the latent variables and then find the member of that family that is closest to the posterior.

Here we focus on mean-field variational inference (though our results apply more widely). First, we posit a family of factorizable probability distributions on latent variables

[TABLE]

This family is called the mean-field family. It represents a joint of the latent variables with $n+d$ (parametric) marginal distributions, $\{q_{\theta_{1}},\ldots,q_{\theta_{d}},q_{z_{1}},\ldots,q_{z_{n}}\}$.

vb finds the member of the family closest to the exact posterior $p(\theta,z\,|\,x)$, where closeness is measured by kl divergence. Thus vb seeks to solve the optimization,

[TABLE]

In practice, vb finds $q^{*}(\theta,z)$ by optimizing an alternative objective, the evidence lower bound (elbo),

[TABLE]

This objective is called the elbo because it is a lower bound on the evidence $\log p(x)$. More importantly, the elbo is equal to the negative KL plus $\log p(x)$, which does not depend on $q(\cdot)$. Maximizing the elbo minimizes the kl (jordan1999introduction).

The optimum $q^{*}(\theta,z)=q^{*}(\theta)q^{*}(z)$ approximates the posterior, and we call it the vb posterior.222For simplicity we will write $q(\theta,z)=\prod^{d}_{i=1}q(\theta_{i})\prod_{j=1}^{n}q(z_{j})$, omitting the subscript on the factors $q(\cdot)$. The understanding is that the factor is indicated by its argument. Though it cannot capture posterior dependence across latent variables, it has hope to capture each of their marginals. In particular, this paper is about the theoretical properties of the vb posterior $q^{*}(\theta)$, the vb posterior of $\theta$. We will also focus on the corresponding expectation of the global variable, i.e., an estimate of the parameter. It is

[TABLE]

We call $\theta^{*}$ the variational Bayes estimate (vbe).

vb methods are fast and yield good predictive performance in empirical experiments (Blei:2016). However, there are few rigorous theoretical results. In this paper, we prove that (1) the vb posterior converges in total variation (tv) distance to the kl minimizer of a normal distribution centered at the truth and (2) the vbe is consistent and asymptotically normal.

These theorems are frequentist in the sense that we assume the data come from $p(x,z\,;\,\theta_{0})$ with a true (nonrandom) $\theta_{0}$. We then study properties of the corresponding posterior distribution $p(\theta\,|\,x)$, when approximating it with variational inference. What this work shows is that the vb posterior is consistent even though the mean field approximating family can be a brutal approximation. In this sense, vb is a theoretically sound approximate inference procedure.

1.1 Main ideas

We describe the results of the paper. Along the way, we will need to define some terms: the variational frequentist estimate (vfe), the variational log likelihood, the vb posterior, the vbe, and the vb ideal. Our results center around the vb posterior and the vbe. (Table 1 contains a glossary of terms.)

The variational frequentist estimate (vfe) and the variational log likelihood.  The first idea that we define is the variational frequentist estimate (vfe). It is a point estimate of $\theta$ that maximizes a local variational objective with respect to an optimal variational distribution of the local variables. (The vfe treats the variable $\theta$ as a parameter rather than a random variable.) We call the objective the variational log likelihood,

[TABLE]

In this objective, the optimal variational distribution $q^{\dagger}(z)$ solves the local variational inference problem,

[TABLE]

Note that $q^{\dagger}(z)$ implicitly depends on both the data $x$ and the parameter $\theta$.

With the objective defined, the vfe is

[TABLE]

It is usually calculated with variational expectation maximization (em) (wainwright2008graphical; ormerod2010explaining), which iterates between the E step of Equation 5 and the M step of Equation 6. Recent research has explored the theoretical properties of the vfe for stochastic block models (bickel2013asymptotic), generalized linear mixed models (hall2011asymptotic), and Gaussian mixture models (westling2015establishing).

We make two remarks. First, the maximizing variational distribution $q^{\dagger}(z)$ of Equation 5 is different from $q^{*}(z)$ in the vb posterior: $q^{\dagger}(z)$ is implicitly a function of individual values of $\theta$, while $q^{*}(z)$ is implicitly a function of the variational distributions $q(\theta)$. Second, the variational log likelihood in Equation 4 is similar to the original objective function for the em algorithm (Dempster:1977). The difference is that the em objective is an expectation with respect to the exact conditional $p(z\,|\,x)$, whereas the variational log likelihood uses a variational distribution $q(z)$.

Variational Bayes and ideal variational Bayes.  While earlier applications of variational inference appealed to variational em and the vfe, most modern applications do not. Rather they use vb, as we described above, where there is a prior on $\theta$ and we approximate its posterior with a global variational distribution $q(\theta)$. One advantage of vb is that it provides regularization through the prior. Another is that it requires only one type of optimization: the same considerations around updating the local variational factors $q(z)$ are also at play when updating the global factor $q(\theta)$.

To develop theoretical properties of vb, we connect the vb posterior to the variational log likelihood; this is a stepping stone to the final analysis. In particular, we define the vb* ideal posterior* $\pi^{*}(\theta\,|\,x)$,

[TABLE]

Here the local latent variables $z$ are constrained under the variational family but the global latent variables $\theta$ are not. Note that because it depends on the variational log likelihood $M_{n}(\theta\,;\,x)$, this distribution implicitly contains an optimal variational distribution $q^{\dagger}(z)$ for each value of $\theta$; see Equations 4 and 5.

Loosely, the vb ideal lies between the exact posterior $p(\theta\,|\,x)$ and a variational approximation $q(\theta)$. It recovers the exact posterior when $p(z\,|\,\theta,x)$ degenerates to a point mass and $q^{\dagger}(z)$ is always equal to $p(z\,|\,\theta,x)$; in that case the variational likelihood is equal to the log likelihood and Equation 7 is the posterior. But $q^{\dagger}(z)$ is usually an approximation to the conditional. Thus the vb ideal usually falls short of the exact posterior.

That said, the vb ideal is more complex that a simple parametric variational factor $q(\theta)$. The reason is that its value for each $\theta$ is defined by the optimization within $M_{n}(\theta\,;\,x)$. Such a distribution will usually lie outside the distributions attainable with a simple family.

In this work, we first establish the theoretical properties of the vb ideal. We then connect it to the vb posterior.

Variational Bernstein–von Mises.  We have set up the main concepts. We now describe the main results.

Suppose the data come from a true (finite-dimensional) parameter $\theta_{0}$. The classical Bernstein–von Mises theorem says that, under certain conditions, the exact posterior $p(\theta\,|\,x)$ approaches a normal distribution, independent of the prior, as the number of observations tends to infinity. In this paper, we extend the theory around Bernstein–von Mises to the variational posterior. Here we summarize our results.

• •

Lemma 1 shows that the vb ideal $\pi^{*}(\theta\,|\,x)$ is consistent and converges to a normal distribution around the vfe. If the vfe is consistent, the vb ideal $\pi^{*}(\theta\leavevmode\nobreak\ |\leavevmode\nobreak\ x)$ converges to a normal distribution whose mean parameter is a random vector centered at the true parameter. (Note the randomness in the mean parameter is due to the randomness in the observations $x$.)

• •

We next consider the point in the variational family that is closest to the vb ideal $\pi^{*}(\theta\leavevmode\nobreak\ |\leavevmode\nobreak\ x)$ in kl divergence. Lemma 2 and Lemma 3 show that this kl minimizer is consistent and converges to the kl minimizer of a normal distribution around the vfe. If the vfe is consistent (bickel2013asymptotic; hall2011asymptotic) then the kl minimizer converges to the kl minimizer of a normal distribution with a random mean centered at the true parameter.

• •

Lemma 4 shows that the vb posterior $q^{*}(\theta)$ enjoys the same asymptotic properties as the kl minimizers of the vb ideal $\pi^{*}(\theta\,|\,x)$.

• •

Theorem 5 is the variational Bernstein–von Mises theorem. It shows that the vb posterior $q^{*}(\theta)$ is asymptotically normal around the vfe. Again, if the vfe is consistent then the vb posterior converges to a normal with a random mean centered at the true parameter. Further, Theorem 6 shows that the vbe $\hat{\theta}_{n}^{*}$ is consistent with the true parameter and asymptotically normal.

• •

Finally, we prove two corollaries. First, if we use a full rank Gaussian variational family then the corresponding vb posterior recovers the true mean and covariance. Second, if we use a mean-field Gaussian variational family then the vb posterior recovers the true mean and the marginal variance, but not the off-diagonal terms. The mean-field vb posterior is underdispersed.

Related work. This work draws on two themes. The first is the body of work on theoretical properties of variational inference. you2014variational and ormerod2014variational studied variational Bayes for a classical Bayesian linear model. They used normal priors and spike-and-slab priors on the coefficients, respectively. wang2004convergence studied variational Bayesian approximations for exponential family models with missing values. wang2005inadequacy and wang2006convergence analyzed variational Bayes in Bayesian mixture models with conjugate priors. More recently, zhang2017theoretical studied mean field variational inference in stochastic block models (sbms) with a batch coordinate ascent algorithm: it has a linear convergence rate and converges to the minimax rate within $\log n$ iterations. sheth2017excess proved a bound for the excess Bayes risk using variational inference in latent Gaussian models. ghorbani2018instability studied a version of latent Dirichlet allocation (lda) and identified an instability in variational inference in certain signal-to-noise ratio (snr) regimes. zhang2017convergence characterized the convergence rate of variational posteriors for nonparametric and high-dimensional inference. pati2017statistical provided general conditions for obtaining optimal risk bounds for point estimates acquired from mean field variational Bayes. alquier2016properties and alquier2017concentration studied the concentration of variational approximations of Gibbs posteriors and fractional posteriors based on PAC-Bayesian inequalities. yang2017alpha proposed $\alpha$-variational inference and developed variational inequalities for the Bayes risk under the variational solution.

On the frequentist side, hall2011theory; hall2011asymptotic established the consistency of Gaussian variational em estimates in a Poisson mixed-effects model with a single predictor and a grouped random intercept. westling2015establishing studied the consistency of variational em estimates in mixture models through a connection to M-estimation. celisse2012consistency and bickel2013asymptotic proved the asymptotic normality of parameter estimates in the sbm under a mean field variational approximation.

However, many of these treatments of variational methods—Bayesian or frequentist—are constrained to specific models and priors. Our work broadens these works by considering more general models. Moreover, the frequentist works focus on estimation procedures under a variational approximation. We expand on these works by proving a variational Bernstein–von Mises theorem, leveraging the frequentist results to analyze vb posteriors.

The second theme is the Bernstein–von Mises theorem. The classical (parametric) Bernstein–von Mises theorem roughly says that the posterior distribution of $\sqrt{n}(\theta-\theta_{0})$ “converges”, under the true parameter value $\theta_{0}$, to $\mathcal{N}(X,1/I(\theta_{0}))$, where $X\sim\mathcal{N}(0,1/I(\theta_{0}))$ and $I(\theta_{0})$ is the Fisher information (ghosh2003bayesian; van2000asymptotic; le1953some; le2012asymptotics). Early forms of this theorem date back to Laplace, Bernstein, and von Mises (laplace1809memoire; bernstein; vonmises). A version also appeared in lehmann2006theory. kleijn2012bernstein established the Bernstein–von Mises theorem under model misspecification. Recent advances include extensions to extensions to semiparametric cases (murphy2000profile; kim2006bernstein; de2009bernstein; rivoirard2012bernstein; bickel2012semiparametric; castillo2012semiparametric; castillo2012semiparametric2; castillo2014bernstein; panov2014critical; castillo2015bernstein; ghosal2017fundamentals) and nonparametric cases (cox1993analysis; diaconis1986consistency; diaconis1997consistency; diaconis1998consistency; freedman1999wald; kim2004bernstein; boucheron2009bernstein; james2008large; johnstone2010high; bontemps2011bernstein; kim2009bernstein; knapik2011bayesian; leahu2011bernstein; rivoirard2012bernstein; castillo2012nonparametric; castillo2013nonparametric; spokoiny2013bernstein; castillo2014bayesian; castillo2014bernstein; ray2014adaptive; panov2015finite; lu2017bernstein). In particular, lu2016gaussapprox proved a Bernstein–von Mises type result for Bayesian inverse problems, characterizing Gaussian approximations of probability measures with respect to the kl divergence. Below, we borrow proof techniques from lu2016gaussapprox. But we move beyond the Gaussian approximation to establish the consistency of variational Bayes.

This paper.  The rest of the paper is organized as follows. Section 2 characterizes theoretical properties of the vb ideal. Section 3 contains the central results of the paper. It first connects the vb ideal and the vb posterior. It then proves the variational Bernstein–von Mises theorem, which characterizes the asymptotic properties of the vb posterior and vb estimate. Section 4 studies three models under this theoretical lens, illustrating how to establish consistency and asymptotic normality of specific vb estimates. Section 5 reports simulation studies to illustrate these theoretical results. Finally, Section 6 concludes with paper with a discussion.

2 The vb ideal

To study the vb posterior $q^{*}(\theta)$, we first study the vb ideal of Equation 7. In the next section we connect it to the vb posterior.

Recall the vb ideal is

[TABLE]

where $M_{n}(\theta\,;\,x)$ is the variational log likelihood of Equation 4. If we embed the variational log likelihood $M_{n}(\theta\,;\,x)$ in a statistical model of $x$, this model has likelihood

[TABLE]

We call it the frequentist variational model. The vb ideal $\pi^{*}(\theta\leavevmode\nobreak\ |\leavevmode\nobreak\ x)$ is thus the classical posterior under the frequentist variational model $\ell(\theta\,;\,x)$; the vfe is the classical maximum likelihood estimate (mle).

Consider the results around frequentist estimation of $\theta$ under variational approximations of the local variables $z$ (bickel2013asymptotic; hall2011asymptotic; westling2015establishing). These works consider asymptotic properties of estimators that maximize $M_{n}(\theta\,;\,x)$ with respect to $\theta$. We will first leverage these results to prove properties of the vb ideal and their kl minimizers in the mean field variational family $\mathcal{Q}^{d}$. Then we will use these properties to study the vb posterior, which is what is estimated in practice.

This section relies on the consistent testability and the local asymptotic normality (lan) of $M_{n}(\theta\,;\,x)$ (defined later) to show the vb ideal is consistent and asymptotically normal. We will then show that its kl minimizer in the mean field family is also consistent and converges to the kl minimizer of a normal distribution in tv distance.

These results are not surprising. Suppose the variational log likelihood behaves similarly to the true log likelihood, i.e., they produce consistent parameter estimates. Then, in the spirit of the classical Bernstein–von Mises theorem under model misspecification (kleijn2012bernstein), we expect the vb ideal to be consistent as well. Moreover, the approximation through a factorizable variational family should not ruin this consistency— point masses are factorizable and thus the limiting distribution lies in the approximating family.

2.1 The vb ideal

The lemma statements and proofs adapt ideas from ghosh2003bayesian; van2000asymptotic; bickel1967asymptotically; kleijn2012bernstein; lu2016gaussapprox to the variational log likelihood. Let $\Theta$ be an open subset of $\mathbb{R}^{d}$. Suppose the observations $x=x_{1:n}$ are a random sample from the measure $P_{\theta_{0}}$ with density $\int p(x,z\leavevmode\nobreak\ |\leavevmode\nobreak\ \theta=\theta_{0})\mathop{}\!\mathrm{d}z$ for some fixed, nonrandom value $\theta_{0}\in\Theta$. $z=z_{1:n}$ are local latent variables, and $\theta=\theta_{1:d}\in\Theta$ are global latent variables. We assume that the density maps $(\theta,x)\mapsto\int p(x,z\leavevmode\nobreak\ |\leavevmode\nobreak\ \theta)\mathop{}\!\mathrm{d}z$ of the true model and $(\theta,x)\mapsto\ell(\theta\,;\,x)$ of the variational frequentist models are measurable. For simplicity, we also assume that for each $n$ there exists a single measure that dominates all measures with densities $\ell(\theta\,;\,x),\theta\in\Theta$ as well as the true measure $P_{\theta_{0}}$.

Assumption 1**.**

We assume the following conditions for the rest of the paper:

(Prior mass) The prior measure with Lebesgue-density $p(\theta)$ on $\Theta$ is continuous and positive on a neighborhood of $\theta_{0}$. There exists a constant $M_{p}>0$ such that $|(\log p(\theta))^{\prime\prime}|\leq M_{p}e^{|\theta|^{2}}$. 2. 2.

(Consistent testability) For every $\epsilon>0$ there exists a sequence of tests $\phi_{n}$ such that

[TABLE]

and

[TABLE] 3. 3.

(Local asymptotic normality (lan)) For every compact set $K\subset\mathbb{R}^{d}$, there exist random vectors $\Delta_{n,\theta_{0}}$ bounded in probability and nonsingular matrices $V_{\theta_{0}}$ such that

[TABLE]

where $\delta_{n}$ is a $d\times d$ diagonal matrix. We have $\delta_{n}\rightarrow 0$ as $n\rightarrow\infty$. For $d=1$, we commonly have $\delta_{n}=1/\sqrt{n}$.

These three assumptions are standard for Bernstein–von Mises theorem. The first assumption is a prior mass assumption. It says the prior on $\theta$ puts enough mass to sufficiently small balls around $\theta_{0}$. This allows for optimal rates of convergence of the posterior. The first assumption further bounds the second derivative of the log prior density. This is a mild technical assumption satisfied by most non-heavy-tailed distributions.

The second assumption is a consistent testability assumption. It says there exists a sequence of uniformly consistent (under $P_{\theta_{0}}$) tests for testing $H_{0}:\theta=\theta_{0}$ against $H_{1}:||\theta-\theta_{0}||\geq\epsilon$ for every $\epsilon>0$ based on the frequentist variational model. This is a weak assumption. For example, it suffices to have a compact $\Theta$ and continuous and identifiable $M_{n}(\theta\,;\,x)$. It is also true when there exists a consistent estimator $T_{n}$ of $\theta$. In this case, we can set $\phi_{n}:=1\{T_{n}-\theta\geq\epsilon/2\}.$

The last assumption is a local asymptotic normality assumption on $M_{n}(\theta\,;\,x)$ around the true value $\theta_{0}$. It says the frequentist variational model can be asymptotically approximated by a normal location model centered at $\theta_{0}$ after a rescaling of $\delta_{n}^{-1}$. This normalizing sequence $\delta_{n}$ determines the optimal rates of convergence of the posterior. For example, if $\delta_{n}=1/\sqrt{n}$, then we commonly have $\theta-\theta_{0}=O_{p}(1/\sqrt{n})$. We often need model-specific analysis to verify this condition, as we do in Section 4. We discuss sufficient conditions and general proof strategies in Section 3.4.

In the spirit of the last assumption, we perform a change-of-variable step:

[TABLE]

We center $\theta$ at the true value $\theta_{0}$ and rescale it by the reciprocal of the rate of convergence $\delta_{n}^{-1}.$ This ensures that the asymptotic distribution of $\tilde{\theta}$ is not degenerate, i.e., it does not converge to a point mass. We define $\pi^{*}_{\tilde{\theta}}(\cdot\leavevmode\nobreak\ |\leavevmode\nobreak\ x)$ as the density of $\tilde{\theta}$ when $\theta$ has density $\pi^{*}(\cdot\leavevmode\nobreak\ |\leavevmode\nobreak\ x)$:

[TABLE]

Now we characterize the asymptotic properties of the vb ideal.

Lemma 1**.**

The vb ideal converges in total variation to a sequence of normal distributions,

[TABLE]

Proof sketch of lemma 1. This is a consequence of the classical finite-dimensional Bernstein–von Mises theorem under model misspecification (kleijn2012bernstein). Theorem 2.1 of kleijn2012bernstein roughly says that the posterior is consistent if the model is locally asymptotically normal around the true parameter value $\theta_{0}$. Here the true data generating measure is $P_{\theta_{0}}$ with density $\int p(x,z\leavevmode\nobreak\ |\leavevmode\nobreak\ \theta=\theta_{0})\mathop{}\!\mathrm{d}z$, while the frequentist variational model has densities $\ell(\theta\,;\,x),\theta\in\Theta$.

What we need to show is that the consistent testability assumption in 1 implies assumption (2.3) in kleijn2012bernstein:

[TABLE]

for every sequence of constants $M_{n}\rightarrow\infty.$. To show this, we mimic the argument of Theorem 3.1 of kleijn2012bernstein, where they show this implication for the iid case with a common convergence rate for all dimensions of $\theta$. See Appendix A for details. ∎

This lemma says the vb ideal of the rescaled $\theta$, $\tilde{\theta}=\delta_{n}^{-1}(\theta-\theta_{0})$, is asymptotically normal with mean $\Delta_{n,\theta_{0}}$. The mean, $\Delta_{n,\theta_{0}}$, as assumed in 1, is a random vector bounded in probability. The asymptotic distribution $\mathcal{N}(\cdot;\Delta_{n,\theta_{0}},V^{-1}_{\theta_{0}})$ is thus also random, where randomness is due to the data $x$ being random draws from the true data generating measure $P_{\theta_{0}}$. We notice that if the vfe, $\hat{\theta}_{n}$, is consistent and asymptotically normal, we commonly have $\Delta_{n,\theta_{0}}=\delta_{n}^{-1}(\hat{\theta}_{n}-\theta_{0})$ with $\mathbb{E}(\Delta_{n,\theta_{0}})=0$. Hence, the vb ideal will converge to a normal distribution with a random mean centered at the true value $\theta_{0}$.

2.2 The KL minimizer of the vb ideal

Next we study the kl minimizer of the vb ideal in the mean field variational family. We show its consistency and asymptotic normality. To be clear, the asymptotic normality is in the sense that the kl minimizer of the vb ideal converges to the kl minimizer of a normal distribution in tv distance.

Lemma 2**.**

The kl minimizer of the vb ideal over the mean field family is consistent: it converges weakly to a point mass in $P_{\theta_{0}}$-probability,

[TABLE]

Proof sketch of lemma 2. The key insight here is that point masses are factorizable. Lemma 1 above suggests that the vb ideal converges in distribution to a point mass. We thus have its kl minimizer also converging to a point mass, because point masses reside within the mean field family. In other words, there is no loss, in the limit, incurred by positing a factorizable variational family for approximation.

To prove this lemma, we bound the mass of $B^{c}(\theta_{0},\eta_{n})$ under $q(\theta)$, where $B^{c}(\theta_{0},\eta_{n})$ is the complement of an $\eta_{n}$-sized ball centered at $\theta_{0}$ with $\eta_{n}\rightarrow 0$ as $n\rightarrow\infty$. In this step, we borrow ideas from the proof of Lemma 3.6 and Lemma 3.7 in lu2016gaussapprox. See Appendix B for details. ∎

Lemma 3**.**

The kl minimizer of the vb ideal of $\tilde{\theta}$ converges to that of $\mathcal{N}(\cdot\,;\,\Delta_{n,\theta_{0}},V^{-1}_{\theta_{0}})$ in total variation: under mild technical conditions on the tail behavior of $\mathcal{Q}^{d}$ (see 2 in Appendix C),

[TABLE]

Proof sketch of lemma 3. The intuition here is that, if the two distribution are close in the limit, their kl minimizers should also be close in the limit. Lemma 1 says that the vb ideal of $\tilde{\theta}$ converges to $\mathcal{N}(\cdot;\Delta_{n,\theta_{0}},V^{-1}_{\theta_{0}})$ in total variation. We would expect their kl minimizer also converges in some metric. This result is also true for the (full-rank) Gaussian variational family if rescaled appropriately.

Here we show their convergence in total variation. This is achieved by showing the $\Gamma$-convergence of the functionals of $q$: $\lx@glossaries@gls@link{acronym}{KL}{\leavevmode\textsc{kl}}(q(\cdot)||\pi^{*}_{\tilde{\theta}}(\cdot\leavevmode\nobreak\ |\leavevmode\nobreak\ x))$ to $\lx@glossaries@gls@link{acronym}{KL}{\leavevmode\textsc{kl}}(q(\cdot)||\mathcal{N}(\cdot\,;\,\Delta_{n,\theta_{0}},V^{-1}_{\theta_{0}}))$, for parametric $q$’s. $\Gamma$-convergence is a classical tool for characterizing variational problems; $\Gamma$-convergence of functionals ensures convergence of their minimizers (dal2012introduction; braides2006handbook). See Appendix C for proof details and a review of $\Gamma$-convergence. ∎

We characterized the limiting properties of the vb ideal and their kl minimizers. We will next show that the vb posterior is close to the kl divergence minimizer of the vb ideal. Section 3 culminates in the main theorem of this paper – the variational Bernstein–von Mises theorem – showing the vb posterior share consistency and asymptotic normality with the kl divergence minimizer of vb ideal.

3 Frequentist consistency of variational Bayes

We now study the vb posterior. In the previous section, we proved theoretical properties for the vb ideal and its kl minimizer in the variational family. Here we first connect the vb ideal to the vb posterior, the quantity that is used in practice. We then use this connection to understand the theoretical properties of the vb posterior.

We begin by characterizing the optimal variational distribution in a useful way. Decompose the variational family as

[TABLE]

where $q(\theta)=\prod_{i=1}^{d}q(\theta_{i})$ and $q(z)=\prod_{i=1}^{n}q(z_{i})$. Denote the prior $p(\theta)$. Note $d$ does not grow with the size of the data. We will develop a theory around vb that considers asymptotic properties of the vb posterior $q^{*}(\theta)$.

We decompose the elbo of Equation 3 into the portion associated with the global variable and the portion associated with the local variables,

[TABLE]

The optimal variational factor for the global variables, i.e., the vb posterior, maximizes the elbo. From the decomposition, we can write it as a function of the optimized local variational factor,

[TABLE]

One way to see the objective for the vb posterior is as the elbo profiled over $q(z)$, i.e., where the optimal $q(z)$ is a function of $q(\theta)$ (Hoffman:2013). With this perspective, the elbo becomes a function of $q(\theta)$ only. We denote it as a functional $\lx@glossaries@gls@link{acronym}{ELBO}{\leavevmode\textsc{elbo}}_{p}(\cdot)$:

[TABLE]

We then rewrite Equation 9 as $q^{*}(\theta)=\operatorname*{arg\,max}_{q(\theta)}\lx@glossaries@gls@link{acronym}{ELBO}{\leavevmode\textsc{elbo}}_{p}(q(\theta))$. This expression for the vb posterior is key to our results.

3.1 kl minimizers of the vb ideal

Recall that the kl minimization objective to the ideal vb posterior is the functional $\lx@glossaries@gls@link{acronym}{KL}{\leavevmode\textsc{kl}}(\cdot||\pi^{*}(\theta\leavevmode\nobreak\ |\leavevmode\nobreak\ x))$. We first show that the two optimization objectives $\lx@glossaries@gls@link{acronym}{KL}{\leavevmode\textsc{kl}}(\cdot||\pi^{*}(\theta\leavevmode\nobreak\ |\leavevmode\nobreak\ x))$ and $\lx@glossaries@gls@link{acronym}{ELBO}{\leavevmode\textsc{elbo}}_{p}(\cdot)$ are close in the limit. Given the continuity of both $\lx@glossaries@gls@link{acronym}{KL}{\leavevmode\textsc{kl}}(\cdot||\pi^{*}(\theta\leavevmode\nobreak\ |\leavevmode\nobreak\ x))$ and $\lx@glossaries@gls@link{acronym}{ELBO}{\leavevmode\textsc{elbo}}_{p}(\cdot)$, this implies the asymptotic properties of optimizers of $\lx@glossaries@gls@link{acronym}{KL}{\leavevmode\textsc{kl}}(\cdot||\pi^{*}(\theta\leavevmode\nobreak\ |\leavevmode\nobreak\ x))$ will be shared by the optimizers of $\lx@glossaries@gls@link{acronym}{ELBO}{\leavevmode\textsc{elbo}}_{p}(\cdot)$.

Lemma 4**.**

The negative kl divergence to the vb ideal is equivalent to the profiled elbo in the limit: under mild technical conditions on the tail behavior of $\mathcal{Q}^{d}$ (see for example 3 in Appendix D), for $q(\theta)\in\mathcal{Q}^{d},$

[TABLE]

Proof sketch of Lemma 4. We first notice that

[TABLE]

Comparing Equation 13 with Equation 10, we can see that the only difference between $-\lx@glossaries@gls@link{acronym}{KL}{\leavevmode\textsc{kl}}(\cdot||\pi^{*}(\theta\leavevmode\nobreak\ |\leavevmode\nobreak\ x))$ and $\lx@glossaries@gls@link{acronym}{ELBO}{\leavevmode\textsc{elbo}}_{p}(\cdot)$ is in the position of $\sup_{q(z)}$. $\lx@glossaries@gls@link{acronym}{ELBO}{\leavevmode\textsc{elbo}}_{p}(\cdot)$ allows for a single choice of optimal $q(z)$ given $q(\theta)$, while $-\lx@glossaries@gls@link{acronym}{KL}{\leavevmode\textsc{kl}}(\cdot||\pi^{*}(\theta\leavevmode\nobreak\ |\leavevmode\nobreak\ x))$ allows for a different optimal $q(z)$ for each value of $\theta$. In this sense, if we restrict the variational family of $q(\theta)$ to be point masses, then $\lx@glossaries@gls@link{acronym}{ELBO}{\leavevmode\textsc{elbo}}_{p}(\cdot)$ and $-\lx@glossaries@gls@link{acronym}{KL}{\leavevmode\textsc{kl}}(\cdot||\pi^{*}(\theta\leavevmode\nobreak\ |\leavevmode\nobreak\ x))$ will be the same.

The only members of the variational family of $q(\theta)$ that admit finite $-\lx@glossaries@gls@link{acronym}{KL}{\leavevmode\textsc{kl}}(q(\theta)||\pi^{*}(\theta\leavevmode\nobreak\ |\leavevmode\nobreak\ x))$ are ones that converge to point masses at rate $\delta_{n}$, so we expect $\lx@glossaries@gls@link{acronym}{ELBO}{\leavevmode\textsc{elbo}}_{p}(\cdot)$ and $-\lx@glossaries@gls@link{acronym}{KL}{\leavevmode\textsc{kl}}(\cdot||\pi^{*}(\theta\leavevmode\nobreak\ |\leavevmode\nobreak\ x))$ to be close as $n\rightarrow\infty.$ We prove this by bounding the remainder in the Taylor expansion of $M_{n}(\theta\,;\,x)$ by a sequence converging to zero in probability. See Appendix D for details. ∎

3.2 The vb posterior

Section 2 characterizes the asymptotic behavior of the vb ideal $\pi^{*}(\theta\leavevmode\nobreak\ |\leavevmode\nobreak\ x)$ and their kl minimizers. Lemma 4 establishes the connection between the vb posterior $q^{*}(\theta)$ and the kl minimizers of the vb ideal $\pi^{*}(\theta\leavevmode\nobreak\ |\leavevmode\nobreak\ x)$. Recall $\operatorname*{arg\,min}_{q(\theta)\in\mathcal{Q}^{d}}\lx@glossaries@gls@link{acronym}{KL}{\leavevmode\textsc{kl}}(q(\theta)||\pi^{*}(\theta\leavevmode\nobreak\ |\leavevmode\nobreak\ x))$ is consistent and converges to the kl minimizer of a normal distribution. We now build on these results to study the vb posterior $q^{*}(\theta)$.

Now we are ready to state the main theorem. It establishes the asymptotic behavior of the vb posterior $q^{*}(\theta)$.

Theorem 5**.**

(Variational Bernstein-von-Mises Theorem)**

The vb posterior is consistent:

[TABLE] 2. 2.

The vb posterior is asymptotically normal in the sense that it converges to the kl minimizer of a normal distribution:

[TABLE]

Here we transform $q^{*}(\theta)$ to $q_{\tilde{\theta}}(\tilde{\theta})$, which is centered around the true $\theta_{0}$ and scaled by the convergence rate; see Equation 8. When $\mathcal{Q}^{d}$ is the mean field variational family, then the limiting vb posterior is normal:

[TABLE]

where $V^{\prime}_{\theta_{0}}$ is diagonal and has the same diagonal terms as $V_{\theta_{0}}$.

Proof sketch of Theorem 5. This theorem is a direct consequence of Lemma 2, Lemma 3, Lemma 4. We need the same mild technical conditions on $\mathcal{Q}^{d}$ as in Lemma 3 and Lemma 4. Equation 15 can be proved by first establishing the normality of the optimal variational factor (see Section 10.1.2 of bishop2006pattern for details) and proceeding with Lemma 8. See Appendix E for details. ∎

Given the convergence of the vb posterior, we can now establish the asymptotic properties of the vbe.

Theorem 6**.**

(Asymptotics of the vbe)

Assume $\int|\theta|^{2}\pi(\theta)\mathop{}\!\mathrm{d}\theta<\infty$. Let $\hat{\theta}_{n}^{*}=\int\theta\cdot q_{1}^{*}(\theta)\mathop{}\!\mathrm{d}\theta$ denote the vbe.

The vbe is consistent:

[TABLE] 2. 2.

The vbe is asymptotically normal in the sense that it converges in distribution to the mean of the kl minimizer:333The randomness in the mean of the kl minimizer comes from $\Delta_{n,\theta_{0}}$. if $\Delta_{n,\theta_{0}}\stackrel{{\scriptstyle d}}{{\rightarrow}}X$ for some $X$,

[TABLE]

Proof sketch of Theorem 6. As the posterior mean is a continuous function of the posterior distribution, we would expect the vbe is consistent given the vb posterior is. We also know that the posterior mean is the Bayes estimator under squared loss. Thus we would expect the vbe to converge in distribution to squared loss minimizer of the kl minimizer of the vb ideal. The result follows from a very similar argument from Theorem 2.3 of kleijn2012bernstein, which shows that the posterior mean estimate is consistent and asymptotically normal under model misspecification as a consequence of the Bernsterin–von Mises theorem and the argmax theorem. See Appendix E for details. ∎

We remark that $\Delta_{n,\theta_{0}}$, as in 1, is a random vector bounded in $P_{\theta_{0}}$ probability. The randomness is due to $x$ being a random sample generated from $P_{\theta_{0}}$.

In cases where vfe is consistent, like in all the examples we will see in Section 4, $\Delta_{n,\theta_{0}}$ is a zero mean random vector with finite variance. For particular realizations of $x$ the value of $\Delta_{n,\theta_{0}}$ might not be zero; however, because we scale by $\delta_{n}^{-1}$, this does not destroy the consistency of vb posterior or the vbe.

3.3 Gaussian vb posteriors

We illustrate the implications of Theorem 5 and Theorem 6 on two choices of variational families: a full rank Gaussian variational family and a factorizable Gaussian variational family. In both cases, the vb posterior and the vbe are consistent and asymptotically normal with different covariance matrices. The vb posterior under the factorizable family is underdispersed.

Corollary 7**.**

Posit a full rank Gaussian variational family, that is

[TABLE]

with $\Sigma$ positive definite. Then

$q^{*}(\theta)\stackrel{{\scriptstyle d}}{{\rightarrow}}\delta_{\theta_{0}}$* in $P_{\theta_{0}}$-probability.* 2. 2.

$||q^{*}_{\tilde{\theta}}(\cdot)-\mathcal{N}(\cdot\,;\,\Delta_{n,\theta_{0}},V^{-1}_{\theta_{0}})||_{\lx@glossaries@gls@link{acronym}{TV}{\leavevmode\textsc{tv}}}\stackrel{{\scriptstyle P_{\theta_{0}}}}{{\rightarrow}}0.$** 3. 3.

$\hat{\theta}_{n}^{*}\stackrel{{\scriptstyle P_{\theta_{0}}}}{{\rightarrow}}\theta_{0}.$** 4. 4.

$\delta_{n}^{-1}(\hat{\theta}_{n}^{*}-\theta_{0})-\Delta_{n,\theta_{0}}=o_{P_{\theta_{0}}}(1)$.

Proof sketch of corollary 7. This is a direct consequence of Theorem 5 and Theorem 6. We only need to show that Lemma 3 is also true for the full rank Gaussian variational family. The last conclusion implies $\delta_{n}^{-1}(\hat{\theta}_{n}^{*}-\theta_{0})\stackrel{{\scriptstyle d}}{{\rightarrow}}X$ if $\Delta_{n,\theta_{0}}\stackrel{{\scriptstyle d}}{{\rightarrow}}X$ for some random variable $X$. We defer the proof to Appendix F. ∎

This corollary says that under a full rank Gaussian variational family, vb is consistent and asymptotically normal in the classical sense. It accurately recovers the asymptotic normal distribution implied by the local asymptotic normality of $M_{n}(\theta\,;\,x)$.

Before stating the corollary for the factorizable Gaussian variational family, we first present a lemma on the kl minimizer of a Gaussian distribution over the factorizable Gaussian family. We show that the minimizer keeps the mean but has a diagonal covariance matrix that matches the precision. We also show the minimizer has a smaller entropy than the original distribution. This echoes the well-known phenomenon of vb algorithms underestimating the variance.

Lemma 8**.**

The factorizable kl minimizer of a Gaussian distribution keeps the mean and matches the precision:

[TABLE]

where $\Sigma_{1}^{*}$ is diagonal with $\Sigma^{*}_{1,ii}=((\Sigma_{1}^{-1})_{ii})^{-1}$ for $i=1,2,...,d$. Hence, the entropy of the factorizable kl minimizer is smaller than or equal to that of the original distribution:

[TABLE]

Proof sketch of Lemma 8. The first statement is consequence of a technical calculation of the kl divergence between two normal distributions. We differentiate the kl divergence over $\mu_{0}$ and the diagonal terms of $\Sigma_{0}$ and obtain the result. The second statement is due to the inequality of the determinant of a positive matrix being always smaller than or equal to the product of its diagonal terms (amir1969product; beckenbach2012inequalities). In this sense, mean field variational inference underestimates posterior variance. See Appendix G for details. ∎

The next corollary studies the vb posterior and the vbe under a factorizable Gaussian variational family.

Corollary 9**.**

Posit a factorizable Gaussian variational family,

[TABLE]

where $\Sigma$ positive definite and diagonal. Then

$q^{*}(\theta)\stackrel{{\scriptstyle d}}{{\rightarrow}}\delta_{\theta_{0}}$* in $P_{\theta_{0}}$-probability.* 2. 2.

$||q^{*}_{\tilde{\theta}}(\cdot)-\mathcal{N}(\cdot\,;\,\Delta_{n,\theta_{0}},V^{\prime-1}_{\theta_{0}})||_{\lx@glossaries@gls@link{acronym}{TV}{\leavevmode\textsc{tv}}}\stackrel{{\scriptstyle P_{\theta_{0}}}}{{\rightarrow}}0,$**

where $V^{\prime}$ is diagonal and has the same diagonal entries as $V_{\theta_{0}}$. 3. 3.

$\hat{\theta}_{n}^{*}\stackrel{{\scriptstyle P_{\theta_{0}}}}{{\rightarrow}}\theta_{0}.$** 4. 4.

$\delta_{n}^{-1}(\hat{\theta}_{n}^{*}-\theta_{0})-\Delta_{n,\theta_{0}}=o_{P_{\theta_{0}}}(1)$.

Proof of corollary 9. This is a direct consequence of Lemma 8, Theorem 5, and Theorem 6. ∎

This corollary says that under the factorizable Gaussian variational family, vb is consistent and asymptotically normal in the classical sense. The rescaled asymptotic distribution for $\tilde{\theta}$ recovers the mean but underestimates the covariance. This underdispersion is a common phenomenon we see in mean field variational Bayes.

As we mentioned, the vb posterior is underdispersed. One consequence of this property is that its credible sets can suffer from under-coverage. In the literature on vb, there are two main ways to correct this inadequacy. One way is to increase the expressiveness of the variational family $\mathcal{Q}$ to one that accounts for dependencies among latent variables. This approach is taken by much of the recent vb literature, e.g. tran2015copula; tran2015variational; ranganath2016hierarchical; ranganath2016operator; liu2016stein. As long as the expanded variational family $\mathcal{Q}$ contains the mean field family, Theorem 5 and Theorem 6 remain true.

Alternative methods to handling underdispersion center around sensitivity analysis and bootstrap. giordano2017covariances identified the close relationship between Bayesian sensitivity and posterior covariance. They estimated the covariance with the sensitivity of the vb posterior means with respect to perturbations of the data. chen2017use explored the use of bootstrap in assessing the uncertainty of a variational point estimate. They also studied the underlying bootstrap theory. giordano2017measuring assessed the clutering stability in Bayesian nonparametric models based on an approximation to the infinitesimal jackknife.

3.4 The lan condition of the the variational log likelihood

Our results rest on 1.3, the lan expansion of the variational log likelihood $M_{n}(\theta\,;\,x)$. For models without local latent variables $z$, their variational log likelihood $M_{n}(\theta\,;\,x)$ is the same as their log likelihood $\log p(x\leavevmode\nobreak\ |\leavevmode\nobreak\ \theta)$. The lan expansion for these models have been widely studied. In particular, iid sampling from a regular parametric model is locally asymptotically normal; it satisfies 1.3 (van2000asymptotic). When models do contain local latent variables, however, as we will see in Section 4, finding the lan expansion requires model-specific characterization.

For a certain class of models with local latent variables, the lan expansion for the (complete) log likelihood $\log p(x,z\leavevmode\nobreak\ |\leavevmode\nobreak\ \theta)$ concurs with the expansion of the variational log likelihood $M_{n}(\theta\,;\,x)$. Below we provide a sufficient condition for such a shared lan expansion. It is satisfied, for example, by the stochastic block model (bickel2013asymptotic) under mild identifiability conditions.

Proposition 10**.**

The log likelihood $\log p(x,z\leavevmode\nobreak\ |\leavevmode\nobreak\ \theta)$ and the variational log likelihood $M_{n}(\theta\,;\,x)$ will have the same lan expansion if:

The conditioned nuisance posterior is consistent under $\delta_{n}$-perturbation at some rate $\rho_{n}$ with $\rho_{n}\downarrow 0$ and $\delta_{n}^{-2}\rho_{n}\rightarrow 0$:

For all bounded, stochastic $h_{n}=O_{P_{\theta_{0}}}(1)$, the conditional nuisance posterior converges as

[TABLE]

where $D^{c}(\theta,\rho_{n})=\{z:d_{H}(z,z_{\text{profile}})\geq\rho_{n}\}$ is the Hellinger ball of radius $\rho_{n}$ around $z_{\text{profile}}=\operatorname*{arg\,max}_{z}p(x,z\leavevmode\nobreak\ |\leavevmode\nobreak\ \theta)$, the maximum profile likelihood estimate of $z$. 2. 2.

$\rho_{n}$* should also satisfy that the likelihood ratio is dominated:*

[TABLE]

where the expectation is taken over $x$.

Proof sketch of Proposition 10. The first condition roughly says the posterior of the local latent variables $z$ contracts faster than the global latent variables $\theta$. The second conidtion is a regularity condition. The two conditions together ensure the log marginal likelihood $\log\int p(x,z\leavevmode\nobreak\ |\leavevmode\nobreak\ \theta)\mathop{}\!\mathrm{d}z$ and the complete log likelihood $\log p(x,z\leavevmode\nobreak\ |\leavevmode\nobreak\ \theta)$ share the same lan expansion. This condition shares a similar flavor with the condition (3.1) of the semiparametric Bernstein–von Mises theorem in bickel2012semiparametric. This implication can be proved by a slight adaptation of the proof of Theorem 4.2 in bickel2012semiparametric: We view the collection of local latent variables $z$ as an infinite-dimensional nuisance parameter.

This proposition is due to the following key inequality. For simplicity, we state the version with only discrete local latent variables $z$:

[TABLE]

The continuous version of this inequality can be easily adapted. The lower bound is due to

[TABLE]

and

[TABLE]

The upper bound is due to the Jensen’s inequality. This inequality ensures that the same lan expansion for the leftmost and the rightmost terms would imply the same lan expansion for the middle term, the variational log likelihood $M_{n}(\theta\,;\,x)$. See Appendix H for details. ∎

In general we can appeal to Theorem 4 of le2012asymptotics to argue for the preservation of the lan condition, showing that if it holds for the complete log likelihood then it holds for the variational log likelihood too. In their terminology, we need to establish the vfe as a “distinguished” statistic.

4 Applications

We proved consistency and asymptotic normality of the variational Bayes (vb) posterior (in total variation (tv) distance) and the variational Bayes estimate (vbe). We mainly relied on the prior mass condition, the local asymptotic normality of the variational log likelihood $M_{n}(x\,;\,\theta)$ and the consistent testability assumption of the data generating parameter.

We now apply this argument to three types of Bayesian models: Bayesian mixture models (bishop2006pattern; murphy2012machine), Bayesian generalized linear mixed models (mcculloch2001generalized; jiang2007linear), and Bayesian stochastic block models (wang1987stochastic; snijders1997estimation; mossel2012stochastic; abbe2015community; Wiggins:2008). For each model class, we illustrate how to leverage the known asymptotic results for frequentist variational approximations to prove asymptotic results for vb. We assume the prior mass condition for the rest of this section: the prior measure of a parameter $\theta$ with Lebesgue density $p(\theta)$ on $\Theta$ is continuous and positive on a neighborhood of the true data generating value $\theta_{0}$. For simplicity, we posit a mean field family for the local latent variables and a factorizable Gaussian variational family for the global latent variables.

4.1 Bayesian Mixture models

The Bayesian mixture model is a versatile class of models for density estimation and clustering (bishop2006pattern; murphy2012machine).

Consider a Bayesian mixture of $K$ unit-variance univariate Gaussians with means $\mu=\{\mu_{1},...,\mu_{K}\}$. For each observation $x_{i},i=1,...,n$, we first randomly draw a cluster assignment $c_{i}$ from a categorical distribution over $\{1,...,K\}$; we then draw $x_{i}$ randomly from a unit-variance Gaussian with mean $\mu_{c_{i}}$. The model is

[TABLE]

For a sample of size $n$, the joint distribution is

[TABLE]

Here $\mu$ is a $K$-dimensional global latent vector and $c_{1:n}$ are local latent variables. We are interested inferring the posterior of the $\mu$ vector.

We now establish asymptotic properties of vb for Bayesian Gaussian mixture model (gmm).

Corollary 11**.**

Assume the data generating measure $P_{\mu_{0}}$ has density $\int p(\mu_{0},c,x)\mathop{}\!\mathrm{d}c$. Let $q^{*}(\mu)$ and $\mu^{*}$ denote the vb posterior and the vbe. Under regularity conditions (A1-A5) and (B1,2,4) of westling2015establishing, we have

[TABLE]

and

[TABLE]

where $\mu_{0}$ is the true value of $\mu$ that generates the data. We have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

*The diagonal matrix $V_{0}(\mu_{0})$ satisfies $(V_{0}(\mu_{0})^{-1})_{ii}=(A(\mu_{0}))_{ii}$. The specification of Gaussian mixture model is invariant to permutation among $K$ components; this corollary is true up to permutations among the $K$ components.

Proof sketch for Corollary 11. The consistent testability condition is satisfied by the existence of a consistent estimate due to Theorem 1 of westling2015establishing. The local asymptotic normality is proved by a Taylor expansion of $m(\mu\,;\,x)$ at $\mu_{0}$. This result then follows directly from our Theorem 5 and Theorem 6 in Section 3. The technical conditions inherited from westling2015establishing allow us to use their Theorems 1 and 2 for properties around variational frequentist estimate (vfe). See Appendix I for proof details. ∎

4.2 Bayesian Generalized linear mixed models

Bayesian generalized linear mixed models (glmms) are a powerful class of models for analyzing grouped data or longitudinal data (mcculloch2001generalized; jiang2007linear).

Consider a Poisson mixed model with a simple linear relationship and group-specific random intercepts. Each observation reads $(X_{ij},Y_{ij}),1\leq i\leq m,1\leq j\leq n$, where the $Y_{ij}$’s are non-negative integers and the $X_{ij}$’s are unrestricted real numbers. For each group of observations $(X_{ij},Y_{ij}),1\leq j\leq n$, we first draw the random effect $U_{i}$ independently from $N(0,\sigma^{2})$. We follow by drawing $Y_{ij}$ from a Poisson distribution with mean $\exp(\beta_{0}+\beta_{1}X_{ij}+U_{i})$. The probability model is

[TABLE]

The joint distribution is

[TABLE]

We establish asymptotic properties of vb in Bayesian Poisson linear mixed models.

Corollary 12**.**

Consider the true data generating distribution $P_{\beta_{0}^{0},\beta_{1}^{0},(\sigma^{2})^{0}}$ with the global latent variables taking the true values $\{\beta_{0}^{0},\beta_{1}^{0},(\sigma^{2})^{0}\}$. Let $q^{*}_{\beta_{0}}(\beta_{0})$, $q^{*}_{\beta_{1}}(\beta_{1})$, $q^{*}_{\sigma^{2}}(\sigma^{2})$ denote the vb posterior of $\beta_{0},\beta_{1},\sigma^{2}$. Similarly, let $\beta_{0}^{*},\beta_{1}^{*},(\sigma^{2})^{*}$ be the vbe s accordingly. Consider $m=O(n^{2})$. Under regularity conditions (A1-A5) of hall2011asymptotic, we have

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here $\phi(\cdot)$ is the moment generating function of $X$.

Also,

[TABLE]

Proof sketch for Corollary 12. The consistent testability assumption is satisfied by the existence of consistent estimates of the global latent variables shown in Theorem 3.1 of hall2011asymptotic. The local asymptotic normality is proved by a Taylor expansion of the variational log likelihood based on estimates of the variational parameters based on equations (5.18) and (5.22) of hall2011asymptotic. The technical conditions inherited from hall2011asymptotic allow us to leverage their Theorem 3.1 for properties of the vfe. The result then follows directly from Theorem 5 and Theorem 6 in Section 3. See Appendix J for proof details. ∎

4.3 Bayesian stochastic block models

Stochastic block models are an important methodology for community detection in network data (wang1987stochastic; snijders1997estimation; mossel2012stochastic; abbe2015community).

Consider $n$ vertices in a graph. We observe pairwise linkage between nodes $A_{ij}\in\{0,1\},1\leq i,j\leq n$. In a stochastic block model, this adjacency matrix is driven by the following process: first assign each node $i$ to one of the $K$ latent classes by a categorical distribution with parameter $\pi$. Denote the class membership as $Z_{i}\in\{1,...,K\}$. Then draw $A_{ij}\sim\textrm{Bernoulli}(H_{Z_{i},Z_{j}})$. The parameter $H$ is a symmetric matrix in $[0,1]^{K\times K}$ that specifies the edge probabilities between two latent classes; the parameter $\pi$ are the proportions of the latent classes. The Bayesian stochastic block model is

[TABLE]

The dependence in stochastic block model is more complicated than the Bayesian gmm or the Bayesian glmm.

Before establishing the result, we reparameterize $(\pi,H)$ by $\theta=(\omega,\nu)$, where $\omega\in\mathbb{R}^{K-1}$ is the log odds ratio of belonging to classes $1,...,K-1$, and $\nu\in\mathbb{R}^{K\times K}$ is the log odds ratio of an edge existing between all pairs of the $K$ classes. The reparameterization is

[TABLE]

The joint distribution is

[TABLE]

where

[TABLE]

We now establish the asymptotic properties of vb for stochastic block models.

Corollary 13**.**

Consider $\nu_{0},\omega_{0}$ as true data generating parameters. Let $q^{*}_{\nu}(\nu),q^{*}_{\omega}(\omega)$ denote the vb posterior of $\nu$ and $\omega$. Similarly, let $\nu^{*},\omega^{*}$ be the vbe. Then

[TABLE]

where $\lambda_{0}=\mathbb{E}_{P_{\nu_{0},\omega_{0}}}$(degree of each node), $(\log n)^{-1}\lambda_{0}\rightarrow\infty.$ $Y_{1}$ and $Y_{2}$ are two zero mean random vectors with covariance matrices $\Sigma_{1}$ and $\Sigma_{2}$, where $\Sigma_{1},\Sigma_{2}$ are known functions of $\nu_{0},\omega_{0}$. The diagonal matrix $V(\nu_{0},\omega_{0})$ satisfies $V^{-1}(\nu_{0},\omega_{0})_{ii}=\text{diag}(\Sigma_{1},\Sigma_{2})_{ii}$. Also,

[TABLE]

The specification of classes in stochastic block model (sbm) is permutation invariant. So the convergence above is true up to permutation with the $K$ classes. We follow bickel2013asymptotic to consider the quotient space of $(\nu,\omega)$ over permutations.

Proof sketch of Corollary 13. The consistent testability assumption is satisfied by the existence of consistent estimates by Lemma 1 of bickel2013asymptotic. The local asymptotic normality,

[TABLE]

for $(\nu_{0},\omega_{0})\in\mathcal{T}$ for compact $\mathcal{T}$ with $\rho_{n}=\frac{1}{n}\mathbb{E}(\text{degree of each node})$, is established by Lemma 2, 3 and Theorem 3 of bickel2013asymptotic. The result then follows directly from our Theorem 5 and Theorem 6 in Section 3. See Appendix K for proof details. ∎

5 Simulation studies

We illustrate the implication of Theorem 5 and Theorem 6 by simulation studies on Bayesian glmm (mccullagh1984generalized). We also study the vb posteriors of latent Dirichlet allocation (lda) (blei2003latent). This is a model that shares similar structural properties with sbm but has no consistency results established for its vfe.

We use two automated inference algorithms offered in Stan, a probabilistic programming system (carpenter2015stan): vb through automatic differentiation variational inference (advi) (kucukelbir2016automatic) and Hamiltonian Monte Carlo (hmc) simulation through No-U-Turn sampler (nuts) (hoffman2014nuts). We note that optimization algorithms used for vb in practice only find local optima.

In both cases, we observe the vb posteriors get closer to the truth as the sample size increases; when the sample size is large enough, they coincide with the truth. They are underdispersed, however, compared with hmc methods.

5.1 Bayesian Generalized Linear Mixed Models

We consider the Poisson linear mixed model studied in Section 4. Fix the group size as $n=10$. We simulate data sets of size $N=$ (50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000). As the size of the data set grows, the number of groups also grows; so does the number of local latent variables $U_{i},1\leq i\leq m$. We generate a four-dimensional covariate vector for each $X_{ij},1\leq i\leq m,1\leq j\leq n$, where the first dimension follows i.i.d $\mathcal{N}(0,1)$, the second dimension follows i.i.d $\mathcal{N}(0,25)$, the third dimension follows i.i.d Bernoulli$(0.4)$, and the fourth dimension follows i.i.d. Bernoulli$(0.8)$. We wish to study the behaviors of coefficient efficients for underdispersed/overdispersed continuous covariates and balanced/imbalanced binary covariates. We set the true parameters as $\beta_{0}=5$, $\beta_{1}=(0.2,-0.2,2,-2)$, and $\sigma^{2}=2$.

Figure 1 shows the boxplots of vb posteriors for $\beta_{0},\beta_{1},$ and $\sigma^{2}$. All vb posteriors converge to their corresponding true values as the size of the data set increases. The box plots present rather few outliers; the lower fence, the box, and the upper fence are about the same size. This suggests normal vb posteriors. This echoes the consistency and asymptotic normality concluded from Theorem 5. The vb posteriors are underdispersed, compared to the posteriors via hmc. This also echoes our conclusion of underdispersion in Theorem 5 and Lemma 8.

Regarding the convergence rate, vb posteriors of all dimensions of $\beta_{1}$ quickly converge to their true value; the vb posteriors center around their true values as long as $N\geq 1000$. The convergence of vb posteriors of slopes for continuous variables ($\beta_{11},\beta_{12}$) are generally faster than those for binary ones ($\beta_{13},\beta_{14}$). The vb posterior of $\sigma^{2}$ shares a similarly fast convergence rate. The vb posterior of the intercept $\beta_{0}$, however, struggles; it is away from the true value until the data set size hits $N=20000$. This aligns with the convergence rate inferred in Corollary 12, $\sqrt{mn}$ for $\beta_{1}$ and $\sqrt{m}$ for $\beta_{0}$ and $\sigma^{2}.$

Computation wise, vb takes orders of magnitude less time than hmc. The performance of vb posteriors is comparable with that from hmc when the sample size is sufficiently large; in this case, we need $N=20000.$

5.2 Latent Dirichlet Allocation

Latent Dirichlet Allocation (lda) is a generative statistical model commonly adopted to describe word distributions in documents by latent topics.

Given $M$ documents, each with $N_{m},m=1,...,M$ words, composing a vocabulary of $V$ words, we assume $K$ latent topics. Consider two sets of latent variables: topic distributions for document $m$, $(\theta_{m})_{K\times 1}$, $m=1,...,M$ and word distributions for topic $k$, $(\phi_{k})_{V\times 1}$, $k=1,...,K$. The generative process is

[TABLE]

The first two rows are assigning priors assigned to the latent variables. $w_{m,j}$ denotes word $j$ of document $m$ and $z_{m,j}$ denotes its assigned topic.

We simulate a data set with $V=100$ sized vocabulary and $K=10$ latent topics in $M=$ (10, 20, 50, 100, 200, 500, 100) documents. Each document has $N_{m}$ words where $N_{m}\stackrel{{\scriptstyle iid}}{{\sim}}$Poi(100). As the number of documents $M$ grows, the number of document-specific topic vectors $\theta_{m}$ grows while the number of topic-specific word vectors $\phi_{k}$ stays the same. In this sense, we consider $\theta_{m},m=1,...,M$ as local latent variables and $\phi_{k},k=1,...,K$ as global latent variables. We are interested in the vb posteriors of global latent variables $\phi_{k},k=1,...,K$ here. We generate the data sets with true values of $\theta$ and $\phi$, where they are random draws from $\theta_{m}\stackrel{{\scriptstyle iid}}{{\sim}}$Dir$((1/K)_{K\times 1})$ and $\phi_{k}\stackrel{{\scriptstyle iid}}{{\sim}}$Dir$((1/V)_{V\times 1})$.

Figure 2 presents the Kullback-Leibler (kl) divergence between the $K=10$ topic-specific word distributions induced by the true $\phi_{k}$’s and the fitted $\phi_{k}$’s by vb and hmc. This kl divergence equals to kl (Mult($\phi_{k}^{0}$)||Mult($\hat{\phi}_{k}))=\sum^{V}_{i=1}\phi_{ki}^{0}(\log\phi_{ki}^{0}-\log\hat{\phi}_{ki})$, where $\phi_{ki}^{0}$ is the $i$th entry of the true $k$th topic and $\hat{\phi}_{ki}$ is the $i$th entry of the fitted $k$th topic.

Figure 2a shows that vb posterior (dark blue) mean kl divergences of all $K=10$ topics get closer to 0 as the number of documents $M$ increase, faster than hmc (light blue). We become very close to the truth as the number of documents $M$ hits 1000. Figure 2b444We only show boxplots for Topic 2 here. The boxplots of other topics look very similar. shows that the boxplots of vb posterior mean kl divergences get closer to 0 as $M$ increases. They are underdispersed compared to hmc posteriors. These align with our understanding of how vb posterior behaves in Theorem 5.

Computation wise, again vb is orders of magnitude faster than hmc. In particular, optimization in vb in our simulation studies converges within 10,000 steps.

6 Discussion

Variational Bayes (vb) methods are a fast alternative to Markov chain Monte Carlo (mcmc) for posterior inference in Bayesian modeling. However, few theoretical guarantees have been established. This work proves consistency and asymptotic normality for variational Bayes (vb) posteriors. The convergence is in the sense of total variation (tv) distance converging to zero in probability. In addition, we establish consistency and asymptotic normality of variational Bayes estimate (vbe). The result is frequentist in the sense that we assume a data generating distribution driven by some fixed nonrandom true value for global latent variables.

These results rest on ideal variational Bayes and its connection to frequentist variational approximations. Thus this work bridges the gap in asymptotic theory between the frequentist variational approximation, in particular the variational frequentist estimate (vfe), and variational Bayes. It also assures us that variational Bayes as a popular approximate inference algorithm bears some theoretical soundness.

We present our results in the classical vb framework but the results and proof techniques are more generally applicable. Our results can be easily generalized to more recent developments of vb beyond Kullback-Leibler (kl) divergence, $\alpha$-divergence or $\chi$-divergence for example (li2016renyi; dieng2017variational). They are also applicable to more expressive variational families, as long as they contain the mean field family. We could also allow for model misspecification, as long as the variational loglikelihood $M_{n}(\theta\,;\,x)$ under the misspecified model still enjoys local asymptotic normality.

There are several interesting avenues for future work. The variational Bernstein–von Mises theorem developed in this work applies to parametric and semiparametric models. One direction is to study the vb posteriors in nonparametric settings. A second direction is to characterize the finite-sample properties of vb posteriors. Finally, we characterized the asymptotics of an optimization problem, assuming that we obtain the global optimum. Though our simulations corroborated the theory, vb optimization typically finds a local optimum. Theoretically characterizing these local optima requires further study of the optimization loss surface.

**Acknowledgements. **  We thank the associate editor and two anonymous reviewers for their constructive comments. We thank Adji Dieng, Prateek Jaiswal, Christian Naesseth, and Dustin Tran for their valuable feedback on our manuscript. We also thank Richard Nickl for pointing us to a key reference. This work is supported by ONR N00014-11-1-0651, DARPA PPAML FA8750-14-2-0009, the Alfred P. Sloan Foundation, and the John Simon Guggenheim Foundation.