Statistical Inference in Mean-Field Variational Bayes

Key words:

Bootstrap; Mean-field approximation; Sampling algorithm; Uncertainty quantification; Variational inference.

Introduction

Variational inference is a popular computational approach for approximating complicated probability densities that often involve intractable integrals and many latent variables arising in complex Bayesian hierarchical models. In variational inference, the complicated target is approximated by a closest member relative to the Kullback-Leibler (KL) divergence in a pre-specified family of tractable densities. In many large-scale machine learning applications including clustering problems , image classification and topic models , variational inference can be orders of magnitude faster than the traditional sampling based approaches such as Markov Chain Monte Carlo (MCMC). In particular, by turning the integration, or sampling, problem into an optimization problem, variational inference can take advantage of modern optimization tools such as stochastic optimization techniques and distributed optimization architecture for further improving its efficiency.

Among various approximating schemes, mean-field approximation is the most common type of variational inference that is conceptually simple, implementation-wise easy and particularly suitable for problems involving large numbers of latent variables. The word “mean-field” is originated from the mean-field theory in physics where despite complex interactions among many particles in a many (infinite) body system, all interactions to any one particle can be approximated by a single averaged effect from a “mean-field”. In variational inference, by restricting the approximating family of the mean-field to be all density functions that are fully factorized over (blocks of) unknown variables, the associated optimization problem of finding a closest density can be efficiently solved via the (block) coordinate ascent algorithm . However, the ease of computation comes at a price of poor approximation as these fully factorized densities in the mean-field family fail to capture any dependence structure among the variables. As noticed in earlier studies , this disregard of dependence structure may lead to undesirable consequences such as under-estimating the uncertainties if the resulting mean-field densities are blindly used for constructing credible intervals for the parameters. Therefore, variational inference, including the mean-field approximation, are primarily used for rapidly obtaining point estimates in complex Bayesian hierarchical models where traditional methods such as EM algorithms and MCMC are either mathematically intractable ($\rm E$-step in the EM) or computationally inefficient (slow mixing in the MCMC).

Despite the great empirical success achieved by variational inference over the past decades, researchers have not developed much general theory explaining why variational approximation, in particular the mean-field approximation, works so well until recently. Some earlier threads of research characterize their statistical properties in specific problems such as Bayesian linear models , Poisson mixed effect models , stochastic block models and normal mixture models , among others. Many of these studies prove estimation consistency and derive convergence rate of a point estimator based on the variational proxy by explicitly analyzing the fixed point equation of the variational optimization problem, or directly analyzing the iterative algorithm for solving the optimization problem. In addition, these analyses require the strong conjugacy assumption on the priors of their models.

More recently, Wang and Blei prove that the KL minimizer in variational Bayes asymptotically approaches a normal limit in regular parametric models. Their proof uses the $\Gamma$ convergence technique, and is based on a crucial local asymptotic normality (LAN) assumption on the variational objective. This LAN assumption implicitly assumes the estimation consistency and may require a case-by-case verification. Three groups provide general conditions for deriving the contraction rate of variational approximation as a probability distribution towards the $\delta$-measure at the true parameter of the data generating model, which includes both regular parametric models and infinite-dimensional nonparametric models and implies estimation consistency. Specifically, focus on models that contain no latent variables, while the theory in can be applied to latent variable models such as normal mixture models. All these results justify the use of variational inference as a valid approach for rapidly obtaining rate-optimal point estimators in complex Bayesian models. However, it remains unclear how good the variational point estimator is when compared to some benchmark, such as the maximum likelihood estimator (MLE) and the posterior mean in regular parametric models—at least theoretically, since the MLE (posterior mean) may be computationally expensive to calculate when the $\rm E$-step (full conditional distributions) does not admit a closed form expression when applying the EM algorithm (Gibbs sampler). In addition, there is little work on how to conduct statistical inference, such as creating credible intervals and performing hypothesis testing in variational procedures.

In this work, we develop a new framework for studying theoretical properties of the mean-field variational approximation and for conducting statistical inference on the model parameters based on the mean-field estimator in parametric models involving latent variables. First, we prove a non-asymptotic result that provides an explicit upper bound on the KL divergence between the mean-field approximation to the marginal posterior of the model parameter and its normal approximation, where the center of the normal is precisely the MLE and the covariance matrix is the diagonal of the inverse of the observed data information matrix (which is the asymptotic covariance of the MLE) plus an extra latent variable information matrix (c.f. Section 2.3). The covaraince structure of the approximating normal limit implies that due to the neglect of the dependence between model parameters and latent variables in the mean-field approximation, the uncertainty under-estimation phenomenon is more severe in models with latent variables than models without (latent variable information is zero). As a direct consequence of the normal approximation, we show that the mean-field variational estimator, defined as the expectation of the model parameter under the mean-field approximation to the posterior distribution, matches the MLE (or posterior mean) up to a higher-order term relative to the root-$n$ convergence rate. In other words, there is no loss of efficiency (at least asymptotically) in terms of the mean squared error criterion in using the mean-field inference for point estimation. Second, we propose a new class of variational weighted likelihood Bootstrap (VWLB) methods for conducting statistical inference on the model parameter via perturbing with random weights the (joint) likelihood function in the mean-field inference in the same spirit as bootstrapping. Interestingly, the VWLB can also be viewed as a new sampling scheme that produces independent samples approximating the marginal posterior of the model parameter. In terms of methodology, VWLB extends the classical ideas of weighted likelihood bootstrap and Bayesian bootstrap to complex Bayesian latent variable models. In terms of computation, the VWLB does not suffer from the slowing mixing issue in MCMC due to the independence of the generated samples and is free of tuning. In addition, unlike the sequential nature of MCMC, sampling via VWLB can be conducted in an embarrassingly parallel manner that has the same time complexity as solving a single variational optimization problem via any distributed learning architecture.

A key ingredient in our proof is a relaxed “triangle inequality” around the projection of the limiting normal approximation to the posterior for the KL divergence when restricted to the mean-field family (c.f. Lemma 5). In particular, the mean-field family is not a convex family of distributions, and a strict triangle inequality around the projection of a distribution onto this family with leading factor one (e.g. Theorem 11.6.1 in ) is no longer true. In addition, previous results only show a slow polynomial decay on the tail probability of the variational approximation to the posterior (when away from the true parameter), while in order to control the KL-divergence between the variational approximation and its normal approximation, we prove a stronger sub-Gaussian type tail bound (square exponential decay) that uses essential structures of the mean-field family via a variational type analysis (c.f. Proof of Lemma 2 in Section 6.1).

Overall, our results reveal that despite uncertainty under-estimation, point estimators from the mean-field variational inference have essentially no loss of efficiency as the maximum likelihood estimator and also attains the Cramér-Rao lower bound in parametric models involving latent variables. In addition, by combining variational inference with bootstrap, the resulting VWLB has the potential of providing a principled and more efficient algorithm for sampling from the posterior in complicated Bayesian hierarchical models.

The rest of the paper is organized as follows. In Section 2, we briefly review the mean-field variational inference for approximating the posterior in a general class of Bayesian latent variable models, and present our theoretical results on the non-asymptotic properties of the mean-field approximation. Motivated by these theoretical developments, we propose in Section 3 a new class of variational weighted likelihood Bootstrap methods for statistical inference via the mean-field approximation, and show in Section 4 the estimation consistency in terms of approximating the target posterior of the model parameter. In Section 5, we provide two simulations studies, one with latent variable and one without for validating the theory and illustrating the method. Proofs of some selected results are provided in Section 6. Further details about the simulation and other proofs are deferred to appendices in the supplement material.

Non-asymptotic analysis of mean-field variational approximation

In this section, we begin with a brief review on the mean-field variational inference for a class of Bayesian latent variable models. Then we provide two perspectives for explaining the mechanism behind the mean-field approximation. After that, we state our main results in this section providing non-asymptotic analysis of the mean-field variational procedure. Our results imply the estimation consistency, and characterize the center and shape of the variational approximation to the exact posterior. In particular, our results reveal that although the mean-field approximation fails to capture the uncertainty, a point estimator obtained as the expectation with respect to the variational distribution matches that of the exact posterior up to high-order terms. This favorable property on the variational mean provides the basis of our inference procedure proposed in the next section.

where $p(X^{n}|\,\theta,\,S^{n})$ is the conditional density function of $X^{n}$ given $S^{n}$, and the joint density $p(X^{n}\,|\,\theta)$ of $S^{n}$ is also parametrized by (a subset of) $\theta$ under $P_{\theta}^{n}$. In other cases, a complex probability model, including the latent Dirichlet allocation and Bayesian hierarchical models, may itself be defined in a hierarchical fashion by first specifying the distribution of the data given latent variables and parameters, and then the latent variable distribution given parameters, as formulated in (1). Due to the negative result on the inconsistency of mean-field variatioanl approximation for general state-space models with non-independent observations with non-independent latent variables, we assume the observation latent variable pair $(X_{i},S_{i})$ to be mutually independent, that is,

where $p(S_{i}|\,\theta)$ denotes the marginal density function of $S_{i}$ parametrized by parameter $\theta$ under $P_{\theta}^{n}$.

In the Bayesian paradigm, we impose a prior distribution, denoted by $\Pi(\cdot)$, on the model parameter $\theta$ over $\Theta$, whose density function is denoted by $\pi(\cdot)$. Figure 1 provides a graphical representation of this Bayesian latent variable model considered in the paper. In this framework, all inference is based on the posterior probability $p(Z^{n}|\,X^{n})$ of the collection of latent variables $Z^{n}=(\theta,\,S^{n})$ given visible variables $X^{n}$. According to Bayes’ theorem, this posterior probability has the following form:

In particular, we are interested in the marginal posterior distribution $\Pi_{n}$ of the parameter $\theta$ by integrating out $S^{n}$ in the joint posterior,

where the marginal likelihood $L(\theta;\,X^{n})$, as a function of $\theta$, is

Unfortunately, in most cases $p(S^{n},\theta\,|\,X^{n})$ and $\Pi_{n}(\cdot)$ in equations (3) and (4) can be inconvenient to use for direct analysis due to the intractable normalization constant $p(X^{n})$ involving multi-dimensional integration. Sampling based procedures such as MCMC algorithms could be computationally inefficient due to the high computational cost and slow mixing. Alternatively, variational inference turns the integration problem into an optimization problem by approximating the target distribution $p(Z^{n}|\,X^{n})$ with a closest member $\widehat{q}_{Z^{n}}$ in a pre-specified family $\Gamma$. Formally, the variational approximation $\widehat{q}_{Z^{n}}$ to $p(Z^{n}|\,X^{n})$ is obtained by solving the following optimization problem,

In particular, we focus on the mean-field approximation where the variational family $\Gamma$ is composed of all fully factorized distributions as

Alternatively, one can apply a block mean-field approximation that preserves dependence structures within some multidimensional components, such as the $d$-dim $\theta=(\theta_{1},\ldots,\theta_{d})$ block, in $Z^{n}$, and our results can be readily applied to this less stringent scheme.

2 Two perspectives of the mean-field approximation

The following decomposition of the KL divergence in (6) reveals the interplay between the parameter $\theta$ and latent variables $S^{n}$ pertaining to the mean-field approximation,

where $p_{\theta,S^{n}}(\cdot\,|X^{n})$ stands for the joint posterior density of $(\theta,\,S^{n})$, $\pi_{n}$ the density function induced from the marginal posterior distribution $\Pi_{n}$ of $\theta$ as in (4), and $p_{S^{n}}(\cdot\,|\,\theta,\,X^{n})$ the conditional posterior density of $S^{n}$ given $\theta$. Therefore, jointly minimizing the KL-divergence over $(q_{\theta},q_{S^{n}})$ is equivalent to first profiling out the nuisance part $q_{S^{n}}$ by minimizing the second term for a fixed $q_{\theta}$, and then finding the primary quantity of interest $q_{\theta}$ that minimizes the resulting “profile divergence”. In particular, this semiparametric profiling perspective plays a crucial role in identifying the limiting center and shape of the variational approximation $\widehat{q}_{\theta}$ in its normal approximation presented in the following subsection. The lemma below provides an explicit expression for this profile divergence, whose proof is provided in Section B.1. Recall that $\pi_{n}$ is the marginal posterior density of $\theta$.

For each fixed density $q_{\theta}$ over $\Theta$, the profile divergence takes the following form:

where $r_{i}(s)=\int_{\Theta}\log p(s\,|\,\theta,\,X_{i})\,q_{\theta}(\theta)\,\differential\theta$ for all $s\in\mathcal{S}$.

Roughly speaking, as $q_{\theta}$ approaches the $\delta$-measure at $\theta^{\ast}$, $r_{i}(s)$ tends to $\log p(s\,|\,\theta^{\ast},X_{i})$ and the second term in the preceding display vanishes. Consequently, minimizing the KL-divergence between the joint distributions over $(\theta,S^{n})$ boils down to minimizing $D(q_{\theta}\,\|\,\pi_{n})$. However, the second term still contributes to the limiting shape of $\widehat{q}_{\theta}$ as we will see in the next subsection.

The original derivation of the variance inference provides an alternative interpretation via Jensen’s inequality. Precisely, using the concavity of $\log(x)$, we can obtain an manageable lower bound to the log normalization constant (called evidence) as,

where $L(q_{Z^{n}})$ is called the evidence lower bound (ELBO, ). In particular, the KL divergence

quantifies the discrepancy between the evidence $\log p(X^{n})$ and its lower bound approximation $L(q_{Z^{n}})$. Consequently, minimizing the KL divergence in optimization problem (6) is equivalent to finding a best $q_{Z^{n}}$ to maximize the ELBO. The KL minimization formulation (6) is convenient for our theoretical analysis, while the ELBO formulation leads to various computational algorithms for implementing the variational inference.

In this paper, our attention toward the model is inference on $\theta$, the model parameter, where our theory and methodology is centered on. Towards this goal, it is helpful to inspect a finer decomposition of the ELBO from ,

which consists of three terms: an integrated (relative to the variational distribution of $\theta$) log-marginal likelihood, the Jensen gap $\Delta_{J}$ due to the mean-field decomposition on latent variables $\{S_{i}\}_{i=1}^{n}$ in approximating the marginal likelihood $p(X^{n}|\,\theta)$ with $\widehat{p(X^{n}|\,\theta)}$, and the KL divergence between the variational distribution $q_{\theta}$ and the prior $\pi(\theta)$. When there is no likelihood approximation with latent variables, the Jensen gap $\Delta_{J}$ term vanishes, and maximizing of the ELBO value in decomposition (10) resembles a regularized $M$-estimation problem of minimizing an objective function composed of a goodness of fit term $\int_{\Theta}-\log p(X^{n}|\,\theta)\,Q_{\theta}(\differential\theta)$ plus a regularizing term $D(q_{\theta}\,||\,\pi_{\theta})$ over all distributions in the variational family $\Gamma$. This perspective is useful in proving the consistency and characterizing the contraction rate of the variational approximation $\widehat{q}_{\theta}$ towards the true parameter $\theta^{\ast}$, as describe in the next subsection.

3 Contraction of mean-field variational approximation and normal approximation

Assumption A1 (Prior continuity and growth): The prior density satisfies $\pi(\theta^{\ast})>0$. In addition, $\log\pi(\theta)$ is differentiable in a neighborhood of $\theta^{\ast}$, and satisfies

We define the following quantity of Hellinger bracketing entropy that provides a measure on the model space complexity.

Definition (Hellinger bracketing entropy): For a set $\mathcal{F}$ of functions over $\mathcal{X}$ and any $\varepsilon>0$, we call a set (of pairs of functions) $\{(f_{j}^{L},f_{j}^{U},j=1,\ldots,N)\}$ a (Hellinger) $\varepsilon$-bracketing of $\mathcal{F}$, if $\int_{\mathcal{X}}\big[(f_{j}^{L})^{1/2}(x)-(f_{j}^{U})^{1/2}(x)\big]^{2}\differential x\leq\varepsilon^{2}$ for $j=1,2,\ldots,N$ and for any $f\in\mathcal{F}$, there is a $j$ such that $f_{j}^{L}\leq f\leq f_{j}^{U}$. The (Hellinger) $\varepsilon$-bracketing metric entropy, denote by $H_{B}(\varepsilon,\mathcal{F})$, is defined as the logarithm of the smallest cardinality of such an $\varepsilon$-bracketing of $\mathcal{F}$.

Assumption A2 (Marginal likelihood regularity):

(Smoothness) The log-marginal likelihood function $l(\theta;x)=\log p(x\,|\,\theta)$ is thrice continuously differentiable with respect to $\theta$.

(Information matrix non-degeneracy) In addition, the order of taking expectation with respect to $P_{\theta^{\ast}}$ and differentiation at $\theta^{\ast}$ is valid so that

The $d\times d$ Fisher information matrix in this display, denoted by $I(\theta^{\ast})$, is positive definite.

(Euclidean metric equivalence) The squared Hellinger distance satisfies that for some constants $(c_{1},c_{2})$,

(Local metric entropy growth) There exists a constant $c_{3}$, such that the Hellinger entropy satisfies

The first three assumptions in A2 are standard regularity conditions for parametric models (c.f. Chapter 1.4 in ). The Euclidean metric equivalence assumption is also made in Theorem 5.1 in as one of their sufficient conditions for proving posterior contraction in parametric models. The last assumption on the local metric entropy assumption is adopted from , and often holds for parametric models (c.f. for examples).

In addition, the $d\times d$ latent (variable) information matrix $I_{s}(\theta)=E_{\theta^{\ast}}[\nabla^{2}l_{s}(\theta;S,X)]$ is locally Lipschitz in the neighborhood $\mathcal{B}(\theta^{\ast};\delta)$ of $\theta^{\ast}$, that is,

Assumption A3 includes regularity conditions on the conditional distribution of latent variables. In particular, by viewing the latent variable $S$ as missing data, and interpreting $I_{s}(\theta^{\ast})$ as the missing data information matrix and $I(\theta^{\ast})$ as the observed data information matrix , we can define the complete data information matrix as $I_{c}(\theta^{\ast})=I(\theta^{\ast})+I_{s}(\theta^{\ast})$. As we will see, the inverse of ${\rm diag}(I_{c}(\theta^{\ast}))$ characterizes the limiting shape of the variational approximation $\widehat{q}_{\theta}$, where the second term $I_{s}(\theta^{\ast})$ causes the extra variance reduction due to the neglect of the posterior dependence between $\theta$ and $S^{n}$.

The following lemma shows that with high probability, the marginal distribution $\widehat{Q}_{\theta}$ obtained from the mean-field variational approximation (6) has a sub-Gaussian tail probability outside an $\varepsilon$-ball centered at the truth $\theta^{\ast}$, for all $\varepsilon\geq\sqrt{\log n/n}$. This exponentially decaying tail behavior is essential for controlling the tail integrals in the proof of our next result that approximates $\widehat{Q}_{\theta}$ with a normal distribution. A proof is provided in Section 6.1.

Under Assumptions A1, A2 and A3, there exist constants $(C_{0},C_{1},C_{2},C_{3})$ such that for any $M\geq 1$ and $\varepsilon_{n}=C_{0}\sqrt{\frac{\log n}{n}}$ it holds with probability at least $1-C_{1}M^{-2}$ that the mean-field approximation $\widehat{Q}_{\theta}$ satisfies

Although we focus on the finite-dimensional parametric model $P_{\theta}$, the proof of this lemma is based on a general treatment under a similar setting as . Therefore, this result can also be extended to mean-field approximations for infinite-dimensional models, for which the same sub-Gaussian tail bound holds for all $\varepsilon$ greater than a benchmark contraction rate $\varepsilon_{n}$ slower than the parametric root-$n$ rate pertaining to the model by making certain assumptions (c.f. conditions (2.2)–(2.4) in ) on the prior thickness and the complexity of model space. In comparison, earlier results on the consistency and convergence rates of variational approximations, such as , only show a polynomially decay $C_{3}\varepsilon_{n}^{2}/\varepsilon^{2}$ on the tail probability $\widehat{Q}_{\theta}\big(\|\theta-\theta^{\ast}\|\geq C_{2}\varepsilon\big)$. The proof of our exponentially decaying bound utilizes the factorization structure (7) of the mean-field approximation, and it is still an interesting open problem whether similar sub-Gaussian type tail bounds hold for a broader class of variational approximations beyond the mean-field.

Let $\widehat{\theta}_{\mbox{\tiny MLE}}$ denote the maximum likelihood estimator of $\theta$,

The classical Bernstein von-Mises (BvM) theorem states that the marginal posterior distribution $\Pi_{n}$ of $\theta$ approaches in the total variation metric to $N\big(\widehat{\theta}_{\mbox{\tiny MLE}},\,[nI(\theta^{\ast})]^{-1}\big)$ as $n\to\infty$ (our Lemma 7 with $Q_{\theta}=\Pi_{n}$ gives a stronger KL divergence version of the BvM theorem). Our next theorem shows that the marginal variational distribution $\widehat{Q}_{\theta}$ can also be approximated by a normal distribution with the same center as $\widehat{\theta}_{\mbox{\tiny MLE}}$, but a different variance-covariance matrix, under the stronger KL-divergence. A proof is deferred to Section 6.2. Recall that $I_{c}(\theta^{\ast})=I(\theta^{\ast})+I_{s}(\theta^{\ast})$ is the complete data information matrix.

Under Assumptions A1, A2 and A3, there exist constants $(C_{4},C_{5})$ such that for any $M\geq 1$ it holds with probability at least $1-C_{4}M^{-2}$ that

where $Q^{\ast}_{VB}=N\big(\widehat{\theta}_{\mbox{\tiny MLE}},\,[nI_{VB}]^{-1}\big)$, and $I_{VB}=\mbox{Diag}\big(I_{c}(\theta^{\ast})\big)$.

Theorem 1 is a non-asymptotic result that applies to any sample size $n\geq 1$. The complementary probability $C_{4}M^{-2}$ decays polynomially in $M$ because we simply apply the Markov inequality with the second order moment assumption on the derivatives of the log-likelihood function. If we instead make a sub-Gaussian type assumption as in , then this remainder probability will be exponentially small in $M^{2}$ as $\exp\{-C_{4}M^{2}\}$. As a special when there is no latent variables ($I_{s}(\theta^{\ast})=0$), Theorem 1 shows that the mean-field approximation $\widehat{Q}_{\theta}$ tends to the normal distribution $N\big(\widehat{\theta}_{\mbox{\tiny MLE}},\,[nI_{VB}(\theta^{\ast})]^{-1}\big)$ whose covariance matrix simply removes all off-diagonal components in $I(\theta^{\ast})$, which is consistent with earlier results such as and explains the overly small variances exhibited by the mean-field approximation due to the neglect of the dependence among components of $\theta$. In the general case of Bayesian latent variable models, Theorem 1 shows that the overly small variances phenomenon is even more severe due to the neglect of the dependence between $\theta$ and $S^{n}$.

In practice, although this mismatch on the covariance structures between the variational approximation and the exact posterior is not a serious issue when doing point estimation, erroneous characterization of uncertainty can be produced. As a consequence, variational inference is widely used for rapidly obtaining a point estimator for the model parameter $\theta$. Let $\widehat{\theta}_{VB}$ to denote the variational posterior mean $\widehat{\theta}_{VB}=\int_{\Theta}\theta\,\widehat{q}_{\theta}(\theta)\,\differential\theta$. The following corollary, as a direct consequence of the normal approximation in Theorem 1, shows that although the shape of the exact posterior $\Pi_{n}$ is not properly captured by $\widehat{Q}_{\theta}$, their centers $\widehat{\theta}_{\mbox{\tiny MLE}}$ and $\widehat{\theta}_{VB}$ match up to $\mathcal{O}(n^{-3/4})$.

Under the conditions and the high-probability event of Theorem 1, there exists a constant $C_{6}$ such that

This corollary implies that there is essentially no loss of efficiency in using the mean-field approximation as a fast approach for obtaining a point estimator in low-dimensional parametric models. Moreover, this interesting finding suggests that we can also conduct statistical inference in mean-field approximation by Bootstrapping the point estimator $\widehat{\theta}_{VB}$.

Statistical inference in mean-field approximation

Motivated by results in the previous section, we propose an inferential framework for mean-field variational Bayes in Bayesian models with latent variables by borrowing the classical idea of weighted likelihood Bootstrap (WLB, ) for approximate Bayesian computation. We begin this section with a brief review on the original WLB as a way to simulate approximately from a posterior distribution when there is no latent variables. After that, we extend the WLB to incorporate latent variables, which further leads to our variational weighted likelihood Bootstrap (VWLB) for approximating the marginal posterior distribution of $\theta$ in the mean-field variational Bayes.

WLB is an extension of the Bayesian Bootstrap from nonparametric models to parametric and semiparametric models by approximating the exact posterior via a random sample of parameter values, each maximizing a weighted likelihood function with random weights. The original WLB proposed in directly operates on the marginal density function $p(\cdot\,|\,\theta)$ of the i.i.d. observations $\{X_{i}\}_{i=1}^{n}$ without introducing the latent variables $\{S_{i}\}_{i=1}^{n}$. More specifically, in WLB the $b$th random sample $\widetilde{\theta}^{(b)}$ of the parameter $\theta$, for $b=1,2,\ldots,B$, is produced by maximizing the following weighted likelihood function obtained from tilting the likelihood function $L(\theta;\,X^{n})$:

where the weights $W^{n,(b)}=(W_{1}^{(b)},W_{2}^{(b)}\ldots,W_{n}^{(b)})$ satisfy the following assumption.

Assumption W (Weight randomness): The weights $\{W_{i}^{(b)}:\,i=1,2,\ldots,n,\,b=1,2,\ldots,B\}$ are i.i.d. copies of a nonnegative random variable $W$ with $E[W]=\mbox{Var}(W)=1$. In addition, $W$ is sub-exponential, that is, there exist some constants $(c_{0},c_{1})$ such that $E[e^{\lambda(W-1)}]\leq e^{c_{0}\lambda^{2}/2}$ holds for all $|\lambda|\leq c_{1}$.

Unlike other commonly used sampling schemes such as MCMC, the random samples $\{\widetilde{\theta}^{(b)}\}_{b=1}^{B}$ from WLB are conditionally independent given the data $X^{n}$, where the extra randomness in $\widetilde{\theta}^{(b)}$ is induced by the distribution of the random weights. Under some mild conditions on the model, one can show that the conditional distribution of $\widetilde{\theta}^{(b)}$ given data $X^{n}$ approaches the exact posterior distribution $p(\theta\,|\,X^{n})$ of $\theta$ as $n\to\infty$ (readers may refer to for more details on the WLB and its accompanied theory). To accommodate this WLB idea to variational inference, it is helpful to also associate the weighted likelihood function $\widetilde{L}^{(b)}(\theta;\,X^{n})$ with a weighted posterior distribution

This weighted posterior can be viewed as a generalization of the fractional posterior by raising the probability density $p(X_{i}\,|\,\theta)$ of $X_{i}$ in the likelihood to a sample specific power $W_{i}^{(b)}$. Similar to the classical BvM theorem, the following proposition shows that the mean $\widetilde{\theta}^{(b)}_{B}$ of the weighted posterior distribution $\widetilde{\Pi}^{(b)}_{n}$ based on the marginal likelihood matches the maximum weighted likelihood estimator $\widetilde{\theta}^{(b)}$ up to a higher-order remainder term.

Under Assumptions A1, A2 and W, there exist constants $(C_{7},C_{8})$, such that for any $M\geq 1$, it holds with probability at least $1-C_{7}M^{-2}$ that

2 Variational weighted likelihood Bootstrap

In this subsection, we propose variational weighed likelihood Bootstrap (VWLB) as a variational approximation method for simulating random samples from the marginal posterior distribution $\Pi_{n}$ of $\theta$ in the Bayesian latent variable model (2), thereby facilitating statistical inference on parameter $\theta$. Motivated by the MLB method described in the previous subsection, we define the weighted joint likelihood function $\widetilde{p}(X^{n},S^{n}|\,\theta)$ and weighted joint posterior density $\widetilde{p}(S^{n},\theta\,|\,X^{n})$ that incorporate latent variables $S^{n}$ as,

where recall that $\widetilde{L}^{(b)}(\theta;\,X^{n})$ is the weighted (marginal) likelihood function defined in (11). Note that the “marginalization” of $S^{i}$ in the denominator of (14) is before raising to the power $W^{(b)}_{i}$. Therefore, the weighted joint posterior density is not properly normalized and strictly speaking, not a real density function. Similar to the optimization problem (6) of variational approximation, we define the weighted variational approximation $\widetilde{q}^{(b)}_{Z^{n}}=\widetilde{q}_{\theta}^{(b)}\otimes\bigotimes_{i=1}^{n}\widetilde{q}^{(b)}_{S_{i}}\in\Gamma$ to $\widetilde{p}^{(b)}(Z^{n}|\,X^{n})$ as

where $\widetilde{D}^{(b)}\big(q_{Z^{n}}(\cdot)\,||\,\widetilde{p}^{(b)}(\cdot\,|\,X^{n})\big)$ calculates the expectation with respect to $q_{Z^{n}}$ of the log-ratio between the tilted $q_{Z^{n}}$ and $\widetilde{p}^{(b)}(Z^{n}\,|\,X^{n})$. It is worthy noticing that in practical implementations, the denominator $\widetilde{p}^{(b)}(z^{n}|X^{n})$ in the preceding display simply contributes a constant independent of $q_{Z^{n}}$ in above objective function, and there is no need to explicitly compute this quantity. We include this term mainly for theoretical purposes—the $\widetilde{D}^{(b)}(\cdot\,||\,\cdot)$ with this term reduces to the usual KL divergence when $W_{i}^{(b)}\equiv 1$. In addition, the weighted variational approximation leads to the following decomposition that generalizes the unweighted KL decomposition formula (8) and plays an essential role in our theoretical analysis,

This identity implies the weighted variational objective function to be nonnegative. Accordingly, this decomposition formula enables us to divide the joint minimization problem (15) into two steps—first profiling out the nuisance part $q_{S^{n}}$ by minimizing the second term for a fixed $q_{\theta}$, and then minimizing the resulting “weighted profile divergence” as a function of $q_{\theta}$. In particular, the first step of minimizing over $q_{S^{n}}$ admits a convenient closed form expression for the theoretical analysis, as summarized in the following lemma.

For each fixed density $q_{\theta}$ over $\Theta$, the weighted profile divergence takes the form as:

where functions $\{r_{i}(\cdot)\}_{i=1}^{n}$ are defined in Lemma 1.

The same remark after Lemma 1 regarding the limiting behavior of $r_{i}$ and its implications as $q_{\theta}$ approaches to the $\delta$-measure at $\theta^{\ast}$ applies to the weighted case. Note that the decomposition formula in the lemma is convenient for deriving a normal approximation to $\widetilde{q}^{(b)}_{\theta}$ and cannot be directly used for practical computation since the weighted marginal posterior $\widetilde{\pi}^{(b)}_{n}$in the first term is computationally intractable. In addition, it is worthy mentioning that a similar weighted ELBO decomposition generalizing (10) holds by replacing the first goodness of fit term and the second Jensen gap term with their weighted counterparts. This weighed ELBO decomposition would be useful for deriving the contraction rate of the weighted variational approximation $\widetilde{q}^{(b)}_{\theta}$ beyond parametric models by adapting the proof techniques in .

Having obtained the weighted variational approximation $\widetilde{q}^{(b)}_{\theta}$ to the weighted posterior $\widetilde{\pi}_{n}$ of $\theta$ with the $b$th random weights $W^{n,(b)}$, we can then proceed as in the usual variational inference by using the variational mean $\widetilde{\theta}^{(b)}_{VB}=\int_{\Theta}\theta\,\widetilde{q}^{(b)}_{\theta}(\theta)\,\differential\theta$ as the $b$th random sample approximately drawn from the marginal posterior $\pi_{n}$, for $b=1,2,\ldots,B$. Unlike MCMC sampling algorithms, the random samples $\{\widetilde{\theta}^{(b)}_{VB}\}_{b=1}^{B}$ from the variational WLB are i.i.d. draws approximately from $\pi_{n}$ given data $X^{n}$. These random samples can be used for statistical inference, such as constructing credible sets and conducting hypothesis testing. Algorithm 1 below summarizes the pseudo-code for implementing the variational WLB.

3 Coordinate ascent algorithm for computation

In this subsection, we discuss computational aspects of the optimization problem (15) in the inner loop of Algorithm 1 via a weighted variant of the coordinate ascent. Coordinate ascent variational inference (CAVI, ) is a popular optimization algorithm tailored for solving (6) in the usual mean-field approximation that is scalable to large datasets. CAVI, as an optimization counterpart of Gibbs sampling, utilizes the special structure of the mean-field solution $\widehat{q}_{Z^{n}}$ to (6) that based upon the optimality, each factor in the decomposition (7) should be proportional to the exponential of the expected log of the joint posterior with respect to the rest factors, which under our notation, simplifies to

where the notation $E_{\widehat{q}_{-\theta_{j}}}$ stands for taking expectation with respect to all factors in (7) except for $\widehat{q}_{\theta_{j}}$ (similar notation convention applies to the weighted case below). CAVI iteratively updates each factor $q_{\theta_{j}}$ or $q_{Z_{i}}$ until convergence. Since the ELBO value is non-decreasing along the iterations, CAVI is guaranteed to converge to a local minimum. In practice, the convergence of the CAVI can be assessed by monitoring the ELBO value, and multiple random initializations can be deployed for finding the global minimum by picking one that yields the highest ELBO value.

Now we generalize the CAVI to its weighted counterpart for solving optimization (15). More specifically, the optimality condition of the optimization problem (15) with random weights $W^{n,(b)}=(W^{(b)}_{i})$ is

which only differs from the optimality condition (17) of the usual mean-field optimization in replacing the joint likelihood function $p(X^{n},S^{n}|\,\theta)$ with its weighted version $\widetilde{p}(X^{n},S^{n}|\,\theta)$. Similarly to the CAVI, we can repeatedly update each $q_{\theta_{j}}$ and $q_{S_{i}}$ until convergence, and a stopping criterion can be based on the change in the weighted ELBO. Again, the weighted CAVI reduces the the CAVI when are weights are identically one. Due to the similarity between (17) and the preceding display, only minor changes are needed in order to implement the weighted CAVI based on existing statistical softwares for the CAVI. Algorithm 2 below summarizes the pseudo-code for the weighted CAVI.

Non-asymptotic analysis of variational weighted likelihood Bootstrap

In this section, we develop theoretical justifications for the statistical inference procedure developed in Section 3. Our results show that unlike the mean-field variational approximation $\widehat{q}_{\theta}$ to the marginal posterior $\pi_{n}$ of $\theta$ that generally underestimates the variance (c.f. Theorem 1), the variational weighted likelihood Bootstrap generates independent random samples of $\theta$ given the data whose distribution approaches $\pi_{n}$ as $n\to\infty$. As a consequence, credible intervals constructed from these random samples of $\theta$ has frequentist coverage approaching to their nominal levels as $n\to\infty$. Our analysis is non-asymptotic and leads to explicit high probability error bounds on the discrepancies.

In this subsection, we investigate the theoretical properties of the weighted variational approximation $\widetilde{q}^{(b)}_{\theta}$ as the optimum of the weighted variational optimization (15). Similar to the study of the mean-field approximation $\widehat{q}_{\theta}$ in Section 2.3, we begin with he contraction property of the weighted posterior distribution $\widetilde{\pi}_{n}^{(b)}$ defined in (12) that leads to the contraction of $\widetilde{q}^{(b)}_{\theta}$ with a sub-Gaussian type tail bound. It proof is provided in Section B.4.

Under Assumptions A1, A2 and W, there exist constants $(C^{\prime}_{0},C^{\prime}_{1},C^{\prime}_{2},C^{\prime}_{3})$ such that for any $M\geq 1$ and $\widetilde{\varepsilon}_{n}=C^{\prime}_{0}\frac{\log n}{\sqrt{n}}$ it holds with probability at least $1-C^{\prime}_{1}M^{-2}$ that the weighted posterior $\widetilde{\Pi}^{(b)}_{n}$ satisfies

In addition, if Assumption A3 is also true, then there exist some constants $(C^{\prime\prime}_{2},C^{\prime\prime}_{3})$ such that under the same high probability event, the mean-field approximation $\widehat{Q}_{\theta}$ satisfies

The proof of (20) involves a uniform control on the weighted likelihood ratio via the bracket entropy, and uses the proof technique of Theorem 2 in . Similar to Lemma 2, the proof of Theorem 2 is not specific to parametric models where $\widetilde{\varepsilon}_{n}\asymp\sqrt{\log n/n}$ and can be extended to general cases such as infinite-dimensional models as long as the bracket entropy of the model space is properly controlled.

Recall that $\widetilde{\theta}^{(b)}$ is the maximizer of the weighted likelihood function,

Let $\widetilde{Q}^{\ast(b)}$ denote the normal distribution $N\big(\widetilde{\theta}^{(b)},\,[nI(\theta^{\ast})]^{-1}\big)$. The next theorem extends the classical BvM theorem to the normal approximation $\widetilde{Q}^{\ast(b)}$ of the weighted posterior distribution $\widetilde{\Pi}_{n}^{(b)}$ of $\theta$.

Under Assumptions A1, A2 and W, there exist constants $(C^{\prime}_{4},C^{\prime}_{5})$ such that for any $M\geq 1$, it holds with probability at least $1-C^{\prime}_{4}M^{-2}$ that

Theorem 3 is a non-asymptotic result providing an explicit upper bound to certain discrepancy measures between $\widetilde{\Pi}_{n}^{(b)}$ and its normal approximation $\widetilde{Q}^{\ast(b)}$. This theorem also extends the classical BvM type results from the total variational metric to the stronger KL-divergence, due to the strong sub-Gaussian tail bound (20) for controlling the expectation of the log-density ratio $\log(\wt\pi_n^{(b)}/\wt q^{\ast(b)})$ relative to $\widetilde{\Pi}_{n}^{(b)}$ outside a $\sqrt{\log n/n}$-neighborhood of $\widetilde{\theta}^{(b)}$.

2 Consistency of variational weighted likelihood Bootstrap

In this subsection, we discuss the consistency of variational weighted likelihood Bootstrap summarized in Algorithm 1 as a new sampling scheme for approximating the marginal posterior distribution $\Pi_{n}$ of $\theta$. First, we present a result on the normal approximation of the weighted mean-field approximation $\widetilde{Q}^{(b)}_{\theta}$.

Under Assumptions A1, A2, A3 and W, there exist constants $(C^{\prime}_{6},C^{\prime}_{7})$ such that for any $M\geq 1$, it holds with probability at least $1-C^{\prime}_{6}M^{-2}$ that

This result shows that the weighted version $\widetilde{Q}^{\ast{(b)}}_{VB}$ shares the same covariance structure as $\widetilde{Q}^{\ast}_{VB}$ in Theorem 1, but the center changes from the MLE $\widehat{\theta}$ to the weighted MLE $\widetilde{\theta}^{(b)}$. Since the diagonal matrix $I_{VB}$ not only ignores all off-diagonal entries in the information matrix $I(\theta^{\ast})$ but also inflates the diagonals by an extra additive term $I_{s}(\theta^{\ast})$ due to the mean-field approximation on the latent variables, statistical inference based on $\widetilde{Q}^{(b)}_{\theta}$ will be erroneous. Fortunately, as the center of $\widetilde{Q}^{(b)}_{\theta}$ approximates the weighted MLE $\widetilde{\theta}^{(b)}$, we may instead conduct inference based on this quantity. Formally, recall that $\widetilde{\theta}^{(b)}_{VB}=\int_{\Theta}\theta\,\widetilde{q}_{\theta}^{(b)}(\theta)\,\differential\theta$ is the weighted variational mean of $Q_{\theta}^{(b)}$ in the $b$th replicate of Algorithm 1.

Under the conditions and the high-probability event of Theorem 4, there exists a constant $C^{\prime}_{8}$ such that

This corollary is the weighted version of Corollary 1 that indicates the closeness between the weighted variational mean of $Q_{\theta}^{(b)}$ and the weighted MLE $\widetilde{\theta}^{(b)}$, whose conditional distribution given $X^{n}$ approximates the sampling distribution of the MLE $\widehat{\theta}_{\mbox{\tiny MLE}}$.

Under the assumptions of Theorem 4, there exist constants $(C^{\prime}_{9},C^{\prime}_{10},C^{\prime}_{11})$ such that for any $M\geq 1$, it holds with probability at least $1-C^{\prime}_{9}M^{-2}$ that

where $Z\sim N(0,I_{d})$ is the $d$-variate standard normal distribution, and the randomness of conditional probability $P(\cdot\,|\,X^{n})$ is on the random weights $(W^{(b)}_{1},W^{(b)}_{2},\ldots,W^{(b)}_{n})$.

The first display in Theorem 5 implies that the conditional cdf of $\sqrt{n}\,[I(\theta^{\ast})]^{1/2}\big(\widetilde{\theta}^{(b)}_{VB}-\widehat{\theta}_{\mbox{\tiny MLE}}\,\big)$ given $X^{n}$ uniformly converges to that of $N(0,I_{d})$ as $n\to\infty$, and the second display implies the uniform convergence of the conditional cdf of $\widetilde{\theta}^{(b)}_{VB}$ to that of the marginal posterior distribution $\Pi(\cdot\,|\,X^{n})$ of $\theta$. As a direct consequence of Theorem 5 and the classical BvM result on the posterior $\Pi(\cdot\,|\,X^{n})$, we may use sample quantiles of $\{\widetilde{\theta}^{(b)}_{VB}\}_{b=1}^{B}$ to construct a credible interval, whose frequentist coverage is at most $\mathcal{O}(\sqrt{\log n}/{n^{1/4}})$ away from its nominal level for sufficiently large $B$.

Numerical study

In this section, we provide two numerical examples, the Gaussian mixture model and the Bayesian linear regression, to evaluate the performance of the variational weighted likelihood Bootstrap for credible interval constructions. We will compare four types of credible intervals with level $\alpha=95\%$. The first interval is based on the sample quantiles of the draws from the Gibbs sampler for sampling from the posterior $\Pi_{n}$. The second interval is based on the quantile of the mean-field variational approximation $\widehat{Q}_{\theta}$. The third interval is based on the sample quantile of the draws $\{\widetilde{\theta}_{VB}^{(b)}\}_{b=1}^{B}$ from our VWLB procedure (Algorithm 1). The last interval is the usual bootstrap interval based on the sample quantile of $\{2\widetilde{\theta}_{VB}^{(b)}-\widehat{\theta}_{VB}\}_{b=1}^{B}$, since by Corollary 1 and Theorem 5, the conditional distribution of $\sqrt{n}(\widetilde{\theta}_{VB}^{(b)}-\widehat{\theta}_{VB})$ given $X^{n}$ provides a good approximation to the sampling distribution of $\sqrt{n}(\widehat{\theta}_{VB}-\theta^{\ast})$. Note that the last two types of intervals are asymptotically equivalent since the limiting distribution of $\sqrt{n}(\widehat{\theta}_{VB}-\theta^{\ast})$ is symmetric.

Gaussian Mixture Model(GMM) is a classical example of latent variable models. We consider the following one-dimensional GMM in this simulation. Assume $(X_{1},X_{2},\ldots,X_{n})$ as i.i.d. sample from the data generating model $P_{\theta}=\sum_{k=1}^{K}\frac{1}{K}N(\mu_{k},1)$, with $K=3$, where the parameter is the centers $\theta=\mu=(\mu_{1},\ldots,\mu_{K})$. This model has a latent variable representation by associating each $X_{i}$ with a latent assignment (variable) $c_{i}$ that follows the categorical distribution: $c_{i}\sim\text{Categorical}(1/K,\ldots,1/K)$. Conditioning on $c_{i}$ and $\mu$, the observation $X_{i}$ follows the normal distribution $N(c_{i}^{T}\mu,1)$. We specify the prior distribution of the parameter $\mu=(\mu_{1},\mu_{2},\mu_{3})$ as i.i.d. $\mu_{k}\sim N(0,\sigma^{2})$ with $\sigma=5$. The true parameter $\mu^{\ast}=(-\Delta,0,\Delta)$, with $\Delta$ ranging from $0$ to $8$.

We apply the mean-field approximation that approximates the joint posterior distribution $p(\{c_{i}\}_{i=1}^{n},\mu\,|\,X^{n})$ with a fully factorized distribution $q_{Z^{n}}=\bigotimes_{k=1}^{K}q_{\mu_{k}}\otimes\bigotimes_{i=1}^{n}q_{c_{i}}$. It turns out that the normal prior of $\mu$ is “conjugate” in the sense that the KL minimizer $\widehat{q}_{Z^{n}}=\bigotimes_{k=1}^{K}\widehat{q}_{\mu_{k}}\otimes\bigotimes_{i=1}^{n}\widehat{q}_{c_{i}}$ as in the optimization problem (6) must be from the same distribution family: each $\widehat{q}_{\mu_{k}}$ is a normal distribution, and each $\widehat{q}_{c_{i}}$ is a categorical distribution. Therefore, we can parametrize them by

for $k=1,2,3$ and $i=1,\ldots,n$. Solving the infinite-dimensional optimization problem (6) boils down to optimizing the objective function over these parameters $\{m_{k},s_{k}\}_{k=1}^{3}$ and $\{(\phi_{i1},\ldots,\phi_{iK})\}_{i=1}^{n}$. Implementation details of the algorithm is provided in Appendix A.1.

We compare the numeric performance of the four credible intervals. In the Gibbs sampler, we take out the first $10,000$ iterations as the burn-in, and fetch $500$ samples every $20$ iterations to reduce the auto-correlation. In our VWLB, we set the number of draws to be $B=500$. Figure 2 reports the averaged coverage probability over three parameters $(\mu_{1},\mu_{2},\mu_{3}\}$ in each type of intervals, under sample size $n=100$, $200$, $500$ and $1000$ respectively, and Table 1 reports the credible interval lengths under $\Delta\in\{0,1,3,5\}$ and $n=500$ (clusters $1$ and $3$ are symmetric, so we only report clusters $1$ and $2$). In Appendix A.1, we provide more details about the individual coverage probability for each of the three normal centers, and the respectively estimated posterior density curves. Interestingly, when $\Delta$ is near $0$, all methods except for the VWLB have degeneracy problem, as the three clusters are no longer distinguishable (our regularity Assumptions A2 and A3 are violated). As we can see, the degeneracy window width decreases as the sample size grows. In the good region where the cluster gap $\Delta$ is sufficiently large so that the normal components are statistically distinguishable and our theory applies, the two credible intervals based on VWLB tends to attain the nominal $95\%$ level as the Gibbs sampler. Moreover, the lengths of the credible intervals based on Gibbs and VWLB are similar. In contrast, the credible intervals directly constructed from the mean-field approximation $\widehat{Q}_{\theta}$ are shorter than those from the Gibbs sampling, and tend to under-estimate the uncertainty until the gap $\Delta$ exceeds 5 where the three normal components become nearly disjoint across all sample sizes.

2 Bayesian Linear Regression

In this example we consider the Bayesian linear regression of the following form:

where $\delta_{a}$ denotes the point mass measure at point $a$. The prior for $\gamma_{j}$’s is specified as i.i.d. Bernoulli$(\xi)$, where $\xi$ follows the Beta prior distribution as Beta$(a_{0},b_{0})$, and the prior for the noise variance $\sigma^{2}$ is the inverse Gamma distribution as IG$(\nu/2,\nu\lambda/2)$. In our simulation, we take $v_{1}=2$, $a_{0}=1$, $b_{0}=1$, $\nu=0.002$ and $\lambda=1$ for the hyperparameters. In this example, there is no latent variable and the parameter is $\theta=(\{(\gamma_{j},\beta_{j})\}_{j=1}^{p},\sigma^{2})$.

We apply the following block mean-field approximation for approximating the joint posterior distribution $p(\{(\gamma_{j},\beta_{j})\}_{j=1}^{p},\sigma^{2}\,|\,Y^{n})$ using the blockwise-factorized family $q_{\theta}=\bigotimes_{j=1}^{p}q_{\gamma_{j},\beta_{j}}\otimes q_{\sigma^{2}}$. It again turns out that the point mass mixture priors on $\beta$ and inverse gamma prior on $\sigma$ is the “conjugate” prior in the sense that in the KL minimizer $\widehat{q}_{\theta}=\bigotimes_{j=1}^{p}\widehat{q}_{\gamma_{j},\beta_{j}}\otimes\widehat{q}_{\sigma^{2}}$, each $\widehat{q}_{\gamma_{j},\beta_{j}}$ is a point mass mixture and $\widehat{q}_{\sigma^{2}}$ is an inverse Gamma. Therefore, we can parametrize them by

for $j=1,\ldots,p$. Therefore, finding $\widehat{Q}_{\theta}$ amounts to optimizing the objective function (6) over these parameters $\{\phi_{j},\mu_{j},\sigma_{j}^{2})\}_{j=1}^{p}$ and $\{c,d\}$ via the coordinate descent algorithm. Implementation details of the algorithm is provided in Appendix A.2, where we have adopted a slightly different but more efficient algorithm .

We generate the data in our simulation as follows. The number of observations is $n=1,000$ and number of covariates $p=10$, with ground truth parameter $\beta^{\ast}=(2,3,2,4,1,2,1,0,0,2)^{T}$. We generate $X$ conforming to an AR(1) process associated with a unit variance white noise process as below: for each $i=1,\ldots,n$,

where $\rho\in[0,1)$ is the auto-correlation, and each $X_{ij}$ has the same marginal distribution as $N(0,(1-\rho^{2})^{-1})$. In the simulation, we consider different settings of the auto-correlation as $\rho=0,0.05,\ldots,0.95$. We draw $B=1000$ samples from the Gibbs sampler and the VWLB sampler (Algorithm 1). When evaluating four credible intervals, we compute the coverage probability based on $1,000$ replicates for each setting. Figure 4 reports the trends of coverage probabilities for two coefficients $\beta_{1}$ and $\beta_{4}$ using the four aforementioned types of credible intervals.

As we can expect, the coverage probabilities of the credible intervals based on the mean-field approximation $\widehat{Q}_{\theta}$ rapidly fall below the nominal level $95\%$ when the auto-correlation is large, as $\widehat{Q}_{\theta}$ completely ignores the (high) correlations among $\{\beta_{j}\}_{j=1}^{p}$ in the joint posterior distribution $\Pi_{n}$. In comparison, the rest three methods exhibit similar patterns and nearly attain the nominal level (horizontal dotted line) across all $\rho$ values. We also report the lengths of the intervals reflecting the estimated uncertainty magnitudes in Table 2.

As we can infer from this table, the degree of uncertainty underestimation in the mean-field approximation $\widehat{Q}_{\theta}$ increases as the correlation among $\{\beta_{j}\}_{j=1}^{p}$ in their joint posterior increases, as the posterior covariance matrix of the regression coefficient vector $\beta$ is approximately proportional to the auto-covariance matrix of the AR(1) process with auto-correlation $\rho$. For example, when $\rho=0$, there is no visible uncertainty underestimation, while when $\rho=0.95$, the variational standard deviation of $\beta_{4}$ reduces to roughly one quarter of the true (marginal) posterior standard deviation. In comparison, the lengths from Gibbs sampler and our VWLB sampler are always close.

Proofs of the main results

In this section, we provide a selective proofs of the main results in the paper, and leave the rest and some technique results to the supplement. To simplify the presentation, we use letter $C$ to denote a generic constant whose value may change from one line to another throughout the proof.

Before showing that the variational approximation $\widehat{Q}_{\theta}$ has the desired sub-Gaussian type tail bound, we first show that the marginal posterior distribution $\Pi_{n}$ of $\theta$ has a similar type bound. In fact, Assumption A2 implies $E_{\theta^{\ast}}\log\frac{p(X\,|\,\theta^{\ast})}{p(X\,|\,\theta)}\leq C\|\theta-\theta^{\ast}\|^{2}$ and $E_{\theta^{\ast}}\big[\log\frac{p(X\,|\,\theta^{\ast})}{p(X\,|\,\theta)}\big]^{2}\leq C\|\theta-\theta^{\ast}\|^{2}$. Therefore, we can apply Theorem 5.1 in to obtain that for any $M\geq 1$, it holds with probability at least $1-CM^{-2}$ that the posterior $\Pi_{n}$ of $\theta$ satisfies

where we slightly strengthen their result by providing the explicit posterior tail bound and by extending it from a single $\varepsilon=M\varepsilon_{n}$ to all $\varepsilon\geq M\varepsilon_{n}$ simultaneously, due to the same argument as the proof of equation (6.10) in .

We will use the optimality of $\widehat{Q}=\widehat{Q}_{Z^{n}}=\widehat{Q}_{\theta}\otimes\widehat{Q}_{S^{n}}$ for optimization problem (6) to prove the desired result. Let $Q|_{A}$ denote the restriction of a probability measure $Q$ onto a set $A\subset\Theta\times\mathcal{S}^{n}$, that is, $Q|_{A}(B)=Q(A\cap B)/Q(A)$ for all $B\subset\Theta\times\mathcal{S}^{n}$. For a fixed $\varepsilon\geq M\varepsilon_{n}$ and each $j=1,2,\ldots,d$, we construct a sequence $\mathcal{Q}=\{Q^{\dagger}_{\lambda}:\,\lambda\in\}$ of distributions as

where $A_{n,j}=\{\theta:\,|\theta_{j}-\theta^{\ast}_{j}|\geq D\varepsilon\}\times\mathcal{S}^{n}$ for some sufficiently large constant $D\geq C$. It is easy to verify that $Q^{\dagger}_{\lambda}$ is a valid distribution belonging to the mean-field family $\Gamma$ for any $\lambda\in$. Let $\widehat{\lambda}=\widehat{Q}(A_{n,j})$, so that $\widehat{Q}=Q^{\dagger}_{\widehat{\lambda}}$. Due to the optimality of $\widehat{Q}$ for minimizing $D\big(Q\,||\,P(\cdot\,|\,X^{n})\big)$ when restricting to the smaller family $\mathcal{Q}$, we obtain $\widehat{\lambda}=\mathop{\rm argmin~}_{\lambda\in}D\big(Q^{\dagger}_{\lambda}\,||\,P(\cdot\,|\,X^{n})\big)$, where we can express

where the second equality is due to mutual orthogonality of $\widehat{Q}|_{A_{n,j}}$ and $\widehat{Q}|_{A_{n,j}^{c}}$. Due to a similar decomposition $P(\differential Z^{n}\,|\,X^{n})=P|_{A_{n,j}}(\differential Z^{n}\,|\,X^{n})\,P(A_{n,j}\,|\,X^{n})+P|_{A_{n,j}^{c}}(\differential Z^{n}\,|\,X^{n})\,P(A_{n,j}^{c}\,|\,X^{n})$ for $P(\differential Z\,|\,X^{n})$, the preceding display can be further decomposed into

where $\beta_{n}=P(A_{n,j}\,|\,X^{n})=\Pi_{n}(|\theta_{j}-\theta^{\ast}_{j}|\geq D\varepsilon)$, $\text{Ber}(\lambda)$ denotes a Bernoulli distribution with success probability $\lambda$, and two constants $(d_{n1},d_{n2})$ independent of $\lambda$ are

Since $\widehat{\lambda}$ minimizes $D\big(Q_{\lambda}^{\dagger}\,\|\,P(\cdot\,|\,X^{n})\big)$, by setting the derivative of (24) as a function of $\lambda$ to be zero, we obtain

where we have used the fact that $d_{n1}\geq 0$ in the last step. On the other hand, from the fact that (24) is equal to $D\big(Q_{\lambda}^{\dagger}\,\|\,P(\cdot\,|\,X^{n})\big)$ at $\lambda=\widehat{\lambda}$ and the non-negativeness of the three terms therein, we obtain $D\big(\widehat{Q}\,||\,P(\cdot\,|\,X^{n})\big)\geq(1-\widehat{\lambda})\,d_{n2}+D\big(\text{Ber}(\widehat{\lambda})\,\|\,\text{Ber}(\beta_{n})\big)$. Combining this with the preceding display, we can reach

where we have used inequality (23) so that $\beta_{n}\leq\Pi_{n}(\|\theta-\theta^{\ast}\|\geq D\varepsilon)\leq e^{-CD^{2}n\varepsilon^{2}}\leq n^{-CD^{2}}\leq 1/2$ for a sufficiently large $D$.

Now we invoke the following lemma for bounding $D\big(\widehat{Q}\,\|\,\Pi_{n})$ from above. Its proof is deferred to Appendix B.9 in the supplement.

Under Assumptions A1, A2 and A3, it holds with probability at least $1-CM^{-2}$ that

Using this lemma, the second display in inequality (25), the bound in inequality (23) on $\beta_{n}$, and the fact that $\lambda\log\lambda+(1-\lambda)\log(1-\lambda)\geq-\log 2$ for any $\lambda\in$, we can obtain that $\widehat{\lambda}\leq C/D^{2}\leq 1/2$ by choosing a sufficiently large $D$. Then a combination of the same Lemma 4, the inequality (23) on $\beta_{n}$, and the first display in inequality (25) on $\widehat{\lambda}$ implies

for a sufficiently large constant $D$, which yields the claimed result on the tail probability of $\widehat{Q}_{\theta}$ via a union bound over $j=1,\ldots,d$ (since $d^{-1/2}\|\theta-\theta^{\ast}\|\leq\max_{j}|\theta_{j}-\theta_{j}^{\ast}|$).

2 Proof of Theorem 1

We focus on illustrating the key proof idea of Theorem 1 in this subsection, and leave the proofs of technical lemmas to the appendix. In addition, the proof of Theorem 4 (weighted version of Theorem 1) will follow the same strategy, and details are also deferred to Appendix B.6.

One main difficulty of the proof lies in the fact that the KL divergence does not satisfy the triangle inequality (unless the approximate family is convex, which is not true in the mean-field case), when specialized to our problem, taking the form as

where recall that $Q^{\ast}_{VB}$ shares the same center $\widehat{\theta}_{\mbox{\tiny MLE}}$ as $N(\widehat{\theta}_{\mbox{\tiny MLE}},\,[nI_{c}(\theta^{\ast})]^{-1})$, but the precision matrix (inverse of covariance matrix) of the former is the diagonal part of the latter. This triangle inequality, if true, would imply the desired bound on $D(\widehat{Q}_{\theta}\,||\,Q^{\ast}_{VB})$. In fact, according to Lemma 5 below, $Q^{\ast}_{VB}$ minimizes $D\big(Q_{\theta}\,||\,N(\widehat{\theta}_{\mbox{\tiny MLE}},\,[nI_{c}(\theta^{\ast})]^{-1}\big)$ over all $Q_{\theta}$ within the mean-field family ($Q_{\theta}$ factorizes into $\bigotimes_{j=1}^{d}Q_{\theta_{j}}$). Moreover, the variational optimum $\widehat{Q}_{\theta}$ also approximately minimizes this divergence since according to Lemma 6, the KL divergence $D\big(Q_{\theta}\,||\,N(\widehat{\theta}_{\mbox{\tiny MLE}},\,[nI_{c}(\theta^{\ast})]^{-1})\big)$ is, for all $Q_{\theta}$ with a similar sub-Gaussian tail as in Lemma 2, close to the “profile divergence” defined in Lemma 1, for which $\widehat{Q}_{\theta}$ is the minimizer since by definition $(\widehat{Q}_{\theta},\widehat{Q}_{S^{n}})$ jointly minimizes the variational objective function (7). In other words, the preceding triangle inequality reveals the local strongly convexity structure of the profile variational objective function (after profiling out $Q_{S^{n}}$) around the populational level minimizer $Q^{\ast}_{VB}$. Unfortunately, such a triangle inequality for KL divergence is not true in general. Therefore, our first step of the proof is to establish a similar “triangle inequality” restricted on the mean-field family with respect to the KL-divergence around the populational level minimizer $Q^{\ast}_{VB}$. After that, we will formalize the above intuition that $\widehat{Q}_{\theta}$ effectively minimizes $D\big(Q_{\theta}\,||\,N(\widehat{\theta}_{\mbox{\tiny MLE}},\,[nI_{c}(\theta^{\ast})]^{-1})\big)$.