Convergence Rates of Variational Posterior Distributions

Introduction

Variational Bayes inference is a popular technique to approximate difficult-to-compute probability posterior distributions. Given a posterior distribution $\Pi(\cdot|X^{(n)})$, and a variational family $\mathcal{S}$, variational Bayes inference seeks a $\widehat{Q}\in\mathcal{S}$ that best approximates $\Pi(\cdot|X^{(n)})$ under the Kullback-Leibler divergence. Though it is not exact Bayes inference, the variational class $\mathcal{S}$ often gives computational advantage and leads to algorithms such as coordinate ascent that can be efficiently implemented on large-scale data sets. Researchers in many fields have used variational Bayes inference to solve real problems. Successful examples include statistical genetics , natural language processing , computer vision , and network analysis , to name a few. We refer the readers to an excellent recent review on this topic.

The goal of this paper is to study the variational posterior distribution $\widehat{Q}$ from a theoretic perspective. We propose general conditions on the prior, the likelihood and the variational class to characterize the convergence rate of the variational posterior to the true data generating process.

Before discussing our results, we give a brief review on the theory of convergence rates of the posterior distributions in the literature. In order that the posterior distribution concentrates around the true parameter with some rate, the “prior mass and testing” framework requires three conditions on the prior and the likelihood: a) The prior is required to put a minimal amount of mass in a neighborhood of the true parameter; b) Restricted to a subset of the parameter space, there exists a testing function that can distinguish the truth from the complement of its neighborhood; c) The prior is essentially supported on the subset described above. Rigorous statements of these three conditions can be found in seminal papers . Earlier versions of these conditions go back to . We also mention another line of work that established posterior rates of convergence using other approaches.

In this paper, we show that under almost the same three conditions, the variational posterior $\widehat{Q}$ also converges to the true parameter, and the rate of convergence is given by

The first term $\epsilon_{n}^{2}$ is the rate of convergence of the posterior distribution $\Pi(\cdot|X^{(n)})$. The second term is the variational approximation error with respect to the class $\mathcal{S}$ under the data generating process $P_{0}^{(n)}$. Since we are able to generalize the “prior mass and testing” theory with the same old conditions, many well-studied problems in the literature can now be revisited under our framework of variational Bayes inference with very similar proof techniques. This will be illustrated with several examples considered in the paper.

Remarkably, for a special class of prior distributions and a corresponding variational class, the second term of (1) will be automatically dominated by $\epsilon_{n}^{2}$ under a modified “prior mass” condition. We illustrate this result by a prior distribution of product measure

As long as there exists a subset $\otimes_{j}\widetilde{\Theta}_{j}\subset\left\{\theta:D_{\rho}\left(P_{0}^{(n)}\|P_{\theta}^{(n)}\right)\leq C_{1}n\epsilon_{n}^{2}\right\}$, such that the prior mass condition

holds together with the testing conditions, then the variational posterior distribution $\widehat{Q}$ converges to the true parameter with the rate $\epsilon_{n}^{2}$. In other words, the variational approximation error term in (1) is dominated under this stronger prior mass condition (2). This is the result of Theorem 2.4. Here, $D_{\rho}(\cdot\|\cdot)$ stands for a Rényi divergence with some $\rho>1$. The implication of the condition (2) is important. It says that as long as the prior satisfies a “prior mass” condition that is coherent with the structure of the variational class, the resulted variational approximation error will always be small compared with the statistical error from the true posterior. Therefore, the condition (2) offers a practical guidance on how to choose a good prior for variational Bayes inference. In addition, as a condition only on the prior mass, (2) is usually very easy to check. This mathematical simplicity is not just for independent priors and the mean-field class. In Section 4, a more general condition is proposed that includes the setting of (2) as a special case.

Besides the general formulation of conditions to ensure convergence of the variational posteriors, several interesting aspects of variational Bayes inference are also discussed in the paper. We show that for a general likelihood with a sieve prior, its mean-field variational approximation of the posterior distribution has an interesting relation to an empirical Bayes procedure. We also show that the empirical Bayes procedure is exactly a variational Bayes procedure using a specially designed variational class. This connection between empirical Bayes and variational Bayes is interesting, and may suggest similar theoretical properties of the two.

Finally, we would like to remark that the general rate (1) for variational posteriors is only an upper bound. It is not always true that the variational posterior has a slower convergence rate than the true posterior. Sometimes the variational posterior may not be a good approximation to the true posterior, but it can still contract faster to the true parameter if additional regularity is imposed by the variational class $\mathcal{S}$. We construct examples in Section 5.2 to illustrate this point.

Statistical properties of variational posterior distributions have also been studied in the literature. A recent work by established Bernstein-von Mises type of results for parametric models. We refer the readers to for other related references on theories for parametric variational Bayes inference. For nonparametric and high-dimensional models, recent work by studied variational approximation to tempered posteriors, where the likelihood $dP_{\theta}^{(n)}/dP_{0}^{(n)}$ is replaced by $\left(dP_{\theta}^{(n)}/dP_{0}^{(n)}\right)^{\alpha}$ for some $\alpha\in(0,1)$. Just as the convergence of tempered posteriors , the convergence of the variational approximation can also be established under generalizations of the prior mass condition. In addition, the paper also studied convergence rates under model misspecification, and the paper considered a more general setting that can handle latent variables, which is quite useful to analyze mixture models. We would like to point out that these results do not apply to the usual posterior distributions with $\alpha=1$. After the first version of our paper was posted, similar results on $\alpha=1$ have also been obtained independently by Some extensions of the results of were later added in the revised version of by the same authors.. An early related work on this topic is by , where the results cover both posterior distributions and their variational approximations. However, the conditions in are rather abstract and are not easy to check in applications.

Organization

The rest of the paper is organized as follows. In Section 2, we formulate the problem and introduce the general conditions that characterize convergence rates of variational posteriors. This section also includes results for the mean-field variational class, where the variational approximation error can be explicitly analyzed. In Section 3, we apply our general theory to three examples that use three different variational classes. Then, in Section 4, for a general class of prior distributions and a mean-field class under a model selection setting, we propose a new prior mass condition that leads to an automatic control of the variational approximation error. In Section 5, we discuss the relation between variational Bayes and empirical Bayes. We also discuss possible situations where the variational posterior outperforms the true posterior in this section. An extension of the main results under model misspecification is also discussed in Section 5. All the proofs will be given in the Appendix.

Notation

Main Results

We start this section by introducing a class of divergence functions.

Let $\rho>0$ and $\rho\neq 1$. The $\rho$-Rényi divergence between two probability measures $P_{1}$ and $P_{2}$ is defined as

The relations between the Rényi divergence and other divergence functions are summarized below.

When $\rho\rightarrow 1$, the Rényi divergence converges to the Kullback-Leibler divergence, defined as

When $\rho=1/2$, the Rényi divergence is related to the Hellinger distance by

When $\rho=2$, the Rényi divergence is related to the $\chi^{2}$-divergence by

and the $\chi^{2}$-divergence is defined as

The total variation distance between two probability measures $P_{1}$ and $P_{2}$ is defined as

The relation among the divergence functions defined above is given by the following proposition (see ).

With the above definitions, the following inequalities hold,

Moreover, the Rényi divergence $D_{\rho}(P_{1}\|P_{2})$ is a non-decreasing function of $\rho$.

Now we are ready to introduce the variational posterior distribution. Given a statistical model $P_{\theta}^{(n)}$ parametrized by $\theta$, and a prior distribution $\theta\sim\Pi$, the posterior distribution is defined by

To address possible computational difficulty of the posterior distribution, variational approximation is a way to find the closest object in a class $\mathcal{S}$ of probability measures to $\Pi(\cdot|X^{(n)})$. The class $\mathcal{S}$ is usually required to be computationally or analytically tractable. The most popular mathematical definition of variational approximation is given through the KL-divergence.

Let $\mathcal{S}$ be a family of distributions. The variational approximation of the posterior is defined as

Just like the posterior distribution $\Pi(\cdot|X^{(n)})$, the variational posterior $\widehat{Q}$ is a data-dependent measure that summarizes information from both the prior and the data. For a variational set $\mathcal{S}$, the corresponding variational posterior can be regarded as the projection of the true posterior onto $\mathcal{S}$ under KL-divergence. When $\mathcal{S}$ is the set of all distributions, $\widehat{Q}$ turns out to be the true posterior $\Pi(\cdot|X^{(n)})$. The choice of the class $\mathcal{S}$ usually determines the difficulty of the optimization (3). In this paper, our main goal is to study the statistical property of the data-dependent measure $\widehat{Q}$ for a general $\mathcal{S}$.

2 Results for General Variational Posteriors

Assume the observation $X^{(n)}$ is generated from a probability measure $P_{0}^{(n)}$, and $\widehat{Q}$ is the variational posterior distribution driven by $X^{(n)}$. The goal of this paper is to analyze $\widehat{Q}$ from a frequentist perspective. In other words, we study statistical properties of $\widehat{Q}$ under $P_{0}^{(n)}$. The first theorem gives conditions that guarantee convergence of the variational posterior $\widehat{Q}$.

Suppose $\epsilon_{n}$ is a sequence that satisfies $n\epsilon_{n}^{2}\geq 1$. Consider a loss function $L(\cdot,\cdot)$, such that for any two probability measures $P_{1}$ and $P_{2}$, $L(P_{1},P_{2})\geq 0$. Let $C,C_{1},C_{2},C_{3}>0$ be constants such that $C>C_{2}+C_{3}+2$. We assume

For any $\epsilon>\epsilon_{n}$, there exists a set $\Theta_{n}(\epsilon)$ and a testing function $\phi_{n}$, such that

For any $\epsilon>\epsilon_{n}$, the set $\Theta_{n}(\epsilon)$ above satisfies

Then for the variational posterior $\widehat{Q}$ defined in (3), we have

for some constant $M$ only depending on $C_{1},C$ and $\rho$, where the quantity $\gamma_{n}^{2}$ is defined as

Conditions (C1)-(C3) resemble the three conditions of “prior mass and testing” in . Interestingly, Theorem 2.1 shows that with a slight modification, these three conditions also lead to the convergence of the variational posterior. The testing conditions (C1) and (C2) are required to hold for all $\epsilon>\epsilon_{n}$. In the prior mass condition (C3), the neighborhood of $P_{0}^{(n)}$ is defined through a Rényi divergence with a $\rho>1$, compared with the KL-divergence used in . According to Proposition 2.1, $D_{\rho}(P_{1}\|P_{2})\geq D(P_{1}\|P_{2})$ for $\rho>1$, so the condition (C3) in our paper is slightly stronger than that in . This stronger “prior mass” condition ensures that the loss $L(P_{\theta}^{(n)},P_{0}^{(n)})$ is exponentially integrable under the true posterior $\Pi(\cdot|X^{(n)})$, which is a key step in the proof of Theorem 2.1. In all the examples considered in this paper, we will check (C3) with $D_{2}(P_{0}^{(n)}\|P_{\theta}^{(n)})$, which turns out to be a very convenient choice.

The convergence rate is the sum of two terms, $\epsilon_{n}^{2}$ and $\gamma_{n}^{2}$. The first term $\epsilon_{n}^{2}$ is the convergence rate of the true posterior $\Pi(\cdot|X^{(n)})$. The second term $\gamma_{n}^{2}$ characterizes the approximation error given by the variational set $\mathcal{S}$. A larger $\mathcal{S}$ means more expressive power given by the variational approximation, and thus the rate of $\gamma_{n}^{2}$ is smaller.

It is worth mentioning that we characterize the convergence of the variational posterior $\widehat{Q}$ through the expected loss $P_{0}^{(n)}\widehat{Q}L(P_{\theta}^{(n)},P_{0}^{(n)})$. Bounds for this quantity are also obtained by independently with a stronger testing condition on the entire space. We remark that convergence in $P_{0}^{(n)}\widehat{Q}L(P_{\theta}^{(n)},P_{0}^{(n)})$ automatically implies that the entire variational posterior distribution concentrates in a neighborhood of the true distribution $P_{0}^{(n)}$ with a radius of the same rate. When the loss function is convex, it also implies the existence of a point estimator that enjoys the same convergence rate. We summarize these results in the next corollary.

Under the same setting of Theorem 2.1, for any diverging sequence $M_{n}\rightarrow\infty$, we have

Furthermore, if the loss $L(P_{\theta}^{(n)},P_{0}^{(n)})$ is convex respect to $\theta$, then the variational posterior mean $\widehat{\theta}=\widehat{Q}\theta$ satisfies

The first result is an application of Markov’s inequality

The second result is directly implied by Jensen’s inequality that

To apply Theorem 2.1 to specific problems, we need to analyze the variational approximation error $\gamma_{n}^{2}=\frac{1}{n}\inf_{Q\in\mathcal{S}}P_{0}^{(n)}D(Q||\Pi(\cdot|X^{(n)}))$ in each individual setting. However, this task may not be trivial for many problems. Now we borrow a technique in to get a useful upper bound for $\gamma_{n}^{2}$. For any $Q\in\mathcal{S}$, we have

where $P_{\Pi}^{(n)}=\int P_{\theta}^{(n)}d\Pi(\theta)$. Then, we obtain the upper bound

Now, it is easy to see that a sufficient condition for the variational posterior to converge at the same rate as the true posterior is

We incorporate this condition into the next theorem.

Suppose $\epsilon_{n}$ is a sequence that satisfies $n\epsilon_{n}^{2}\geq 1$, for which the conditions (C1), (C2), (C3), (C4) hold. Then, for the variational posterior $\widehat{Q}$ that is defined in (3), we have

We would like to remark that the quantity $\inf_{Q\in\mathcal{S}}R(Q)$ is easier to analyze compared with the original definition of $\gamma_{n}^{2}$. According to its definition given by (5), it is sufficient to find a distribution $Q\in\mathcal{S}$, such that

These are exactly the two conditions formulated by as a natural extension of the prior mass condition. The relation between the prior mass condition and (7) has also been discussed in .

One way to construct such a distribution $Q$ that satisfies the above two inequalities is to focus on those whose supports are within the set $\mathcal{C}=\{\theta:D(P_{0}^{(n)}\|P_{\theta}^{(n)})\leq Cn\epsilon_{n}^{2}\}$ for some constant $C>0$. We summarize this method into the following theorem.

Suppose there exist constants $C_{1},C_{2}>0$, such that

where $\mathcal{E}=\{Q:{\rm supp}(Q)\subset\mathcal{C}\}$ with $\mathcal{C}=\{\theta:D(P_{0}^{(n)}\|P_{\theta}^{(n)})\leq C_{2}n\epsilon_{n}^{2}\}$. Then, we have

3 Results for Mean-Field Variational Posteriors

A special choice of $\mathcal{S}$ is the mean-field class of distributions. Not only does this class leads to computationally efficient algorithms such as coordinate ascent, but in this section, we will also show that the structure of this class leads to a convenient convergence analysis. We begin with its definition.

For parameters in a product space that can be written as $\theta=(\theta_{1},\theta_{2},...,\theta_{m})$ with some $1\leq m\leq\infty$, the mean-field variational family is defined as

The following theorem can be viewed as an application of Theorem 2.3 to the mean-field class.

Suppose there exists a $\widetilde{Q}\in\mathcal{S}_{\rm MF}$ and a subset $\otimes_{j=1}^{m}\widetilde{\Theta}_{j}$, such that

for some constants $C_{1},C_{2},C_{3}>0$. Then, we have

Note that the condition (9) can also be written as

In other words, Theorem 2.4 gives an interesting “distribution mass” type of characterization for $\inf_{Q\in\mathcal{S}}R(Q)$. Checking (9) is very similar to checking the “prior mass” condition (C3), and is usually not hard in many examples. We only need to make sure that $\widetilde{Q}$ is not too far away from the prior $\Pi$ in the sense of (8). In fact, if the prior $\Pi$ belongs to the class $\mathcal{S}_{\rm MF}$, then one can take $\widetilde{Q}=\Pi$, and the conditions of Theorem 2.4 simply become a “prior mass” condition $\Pi\left(\otimes_{j=1}^{m}\widetilde{\Theta}_{j}\right)\geq\exp\left(-C_{3}n\epsilon_{n}^{2}\right)$, with the choice of $\otimes_{j=1}^{m}\widetilde{\Theta}_{j}$ being a subset of the KL-neighborhood $\left\{\theta:D(P_{0}^{(n)}\|P_{\theta}^{(n)})\leq C_{1}n\epsilon_{n}^{2}\right\}$. A more general characterization of the variational approximation error under model selection setting through a prior mass condition will be studied in Section 4.

Applications

In this section, we consider several examples to illustrate the theory developed in Section 2.

Consider observations generated by a Gaussian sequence model,

We use the notation $P_{\theta}^{(n)}=\otimes_{j}N(\theta_{j},n^{-1})$ for the distribution above. Our goal is to use variational Bayes methods to estimate the true parameter $\theta^{*}$ that belongs to the following Sobolev ball,

Here, the smoothness $\alpha>0$ and the radius $B>0$ are considered as constants throughout the paper. The loss function for this problem is $L(P_{\theta}^{(n)},P_{\theta^{*}}^{(n)})=n\|\theta-\theta^{*}\|^{2}$, which is a natural choice for the Gaussian sequence model.

The prior distribution $\theta\sim\Pi$ is described through the following sampling process.

Conditioning on $k$, sample $\theta_{j}\sim f_{j}$ for all $j\in[k]$, and set $\theta_{j}=0$ for all $j>k$.

In other words, the prior on $\theta$ is a mixture of product measures,

Priors of similar forms are also considered in . Direct calculation implies that the posterior is also in the form of a mixture of product measures.

Consider the variational posterior $\widehat{Q}$ defined by (3) with $\mathcal{S}=\mathcal{S}_{\rm MF}$. That is, we seek a data-dependent measure in a more tractable form of a product measure. In most cases, the variational posterior does not have a closed form and needs to be solved by coordinate ascent algorithms . However, for the Gaussian sequence model (10) with the prior distribution (12), one can write down the exact form of the mean-field variational posterior distribution.

Consider the variational posterior $\widehat{Q}$ induced by the likelihood (10), the prior (12) and the mean-field variational set $\mathcal{S}_{\rm MF}$. The distribution $\widehat{Q}$ is a product measure with the density of each coordinate specified by

The number $\pi(k|Y)$ is the posterior probability of the model dimension, and according to Bayes formula, it is

In other words, the mean-field variational posterior $\widehat{Q}$ is nearly equivalent to a thresholding rule. It estimates all $\theta_{j}^{*}$ by $0$ after $\widetilde{k}$ and applies the usual posterior distribution for each coordinate before $\widetilde{k}$. A mixed strategy is applied to the $\widetilde{k}$th coordinate. The effective model dimension $\widetilde{k}$ is found in a data-driven way through (14).

Next, we will show that even though the posterior itself is not a product measure, using $\widehat{Q}$ from the mean-field class still gives us a rate-optimal contraction result. The conditions on the prior distributions are summarized below.

There exist some constants $C_{1},C_{2}>0$ such that

There exist some constants $C_{3},C_{4}>0$ such that for $k_{0}=\left\lceil\left(\frac{n}{\log n}\right)^{\frac{1}{2\alpha+1}}\right\rceil$,

These three conditions on $\Pi$ include a large class of prior distributions. We remark that even though (17) involves $\alpha$, it does not mean that one needs to know $\alpha$ when defining the prior $\Pi$. For example, the choice that $\pi(k)\propto e^{-\tau k}$ and $f_{j}$ being $N(0,\sigma^{2})$ for some constants $\tau,\sigma^{2}>0$ easily satisfies all the three conditions (15)-(17).

Conditions (15)-(17) will be used to derive the four conditions in Theorem 2.2. To be specific, (C1) and (C2) are consequences of (15) (see Lemma B.7 in the appendix), and (C3) and (C4) can be derived from (16) and (17) (see Lemma B.8 in the appendix). Then, by Theorem 2.2, we obtain the following result.

Consider the prior $\Pi$ that satisfies (15)-(17). Then, for any $\theta^{*}\in\Theta_{\alpha}(B)$, we have

where $\widehat{Q}$ is the variational posterior defined by (3) with $\mathcal{S}=\mathcal{S}_{\rm MF}$.

It is well known that the minimax rate of estimating $\theta^{*}$ in $\Theta_{\alpha}(B)$ is $n^{-\frac{2\alpha}{2\alpha+1}}$ . Using a mean-field variational posterior, we achieve the minimax rate up to a logarithmic factor. In fact, the following proposition demonstrates that this rate cannot be improved for a very general class of priors.

Consider the prior $\Pi$ specified in (12). Assume that $\max_{j}\|f_{j}\|_{\infty}\leq a$ and $\pi(k)$ is nonincreasing over $k$. Then, we have

where $\widehat{Q}$ is the variational posterior defined by (3) with $\mathcal{S}=\mathcal{S}_{\rm MF}$.

On the other hand, the extra logarithmic factor can actually be removed by a rescaling of the prior. Details of this improvement are given in Appendix A.1.

2 Infinite Dimensional Exponential Families

In this section, we study another interesting variational family. The Gaussian mean-field family is defined as

This class offers better interpretability of the results because every distribution in $\mathcal{S}_{G}$ is fully determined by a sequence of mean and variance parameters. Note that we allow $\sigma_{j}^{2}$ to be zero and $N(\mu_{j},0)$ is understood as the delta measure $\delta_{\mu_{j}}$ on $\mu_{j}$.

The application of $\mathcal{S}_{G}$ is illustrated by an infinite dimensional exponential family model. We define the probability measure $P_{\theta}$ by

Since $\phi_{0}(x)=1$ and $\theta_{0}$ can take arbitrary values without changing $P_{\theta}$, we simply set $\theta_{0}=0$. In other words, $P_{\theta}$ is fully parameterized by $\theta=(\theta_{1},\theta_{2},...)$. Given i.i.d. observations from $P_{\theta^{*}}^{n}$, our goal is to estimate $P_{\theta^{*}}$, where $\theta^{*}$ is assumed to belong to the Sobolev ball $\Theta_{\alpha}(B)$ defined in (11). The loss function is chosen as $n$ times the squared Hellinger distance $L(P_{\theta}^{n},P_{\theta^{*}}^{n})=nH^{2}(P_{\theta},P_{\theta^{*}})$.

We consider a prior distribution $\Pi$ that is similar to the one used in Section 3.1. Its sampling process is described as follows.

Conditioning on $k$, sample $\theta_{j}\sim f_{j}$ for all $j\in[k]$, and set $\theta_{j}=0$ for all $j>k$.

We impose the following conditions on the prior $\Pi$.

There exist some constants $C_{1},C_{2}>0$ such that

There exist some constants $C_{3},C_{4}>0$ such that for $k_{0}=\left\lceil\left(\frac{n}{\log n}\right)^{\frac{1}{2\alpha+1}}\right\rceil$

The conditions (20)-(23) are satisfied by a large class of prior distributions. For example, one can choose $k\sim\text{Poisson}(\tau)$ and $f_{j}$ being the density of $N(0,\sigma^{2})$ for some constants $\tau,\sigma^{2}>0$, and then the four conditions are easily satisfied.

Consider the prior $\Pi$ that satisfies (20)-(23). Then, for any $\theta^{*}\in\Theta_{\alpha}(B)$ with some $\alpha>1/2$, we have

where $\widehat{Q}$ is the variational posterior defined by (3) with $\mathcal{S}=\mathcal{S}_{\rm G}$.

The theorem shows that the Gaussian mean-field variational posterior is able to achieve the minimax rate $n^{-\frac{2\alpha}{2\alpha+1}}$ up to a logarithmic factor. We remark that the same result also holds for the mean-field variational posterior defined with $\mathcal{S}_{\rm MF}$. This is because $\mathcal{S}_{\rm G}\subset\mathcal{S}_{\rm MF}$, and thus $\inf_{Q\in\mathcal{S}_{\rm MF}}R(Q)\leq\inf_{Q\in\mathcal{S}_{\rm G}}R(Q)$. Compared with the class $\mathcal{S}_{\rm MF}$, the objective function using the parametric family $\mathcal{S}_{\rm G}$ can be optimized by algorithms such as stochastic gradient descent over the parameters $(\mu_{j},\sigma_{j}^{2})$. The objective function can be greatly simplified according to the general mean-field solution given in Theorem 5.1.

3 Piecewise Constant Model

The previous two sections consider examples of the mean-field variational set and its variant. In this section, we use another example to illustrate a situation where the mean-field variational set only gives a trivial rate. On the other hand, we show that alternative variational classes with appropriate dependence structures are able to achieve the optimal rate.

We consider the following piecewise constant model,

where $Z_{i}\sim N(0,1)$ independently for all $i\in[n]$. We assume $n\geq 2$ throughout the section. The true parameter $\theta^{*}$ is assumed to belong to the class $\Theta_{k^{*}}(B)=\left\{\theta\in\Theta_{k^{*}}:\|\theta\|_{\infty}\leq B\right\}$, where for a general $k\in[n]$,

Here for any two integers $a<b$, we use $(a:b]$ to denote all integers from $a+1$ to $b$. We assume both $B>0$ and $\sigma^{2}>0$ are constants throughout this section. A vector $\theta^{*}\in\Theta_{k^{*}}(B)$ is a piecewise constant signal with at most $k^{*}$ pieces. We use $P_{\theta}^{(n)}$ to denote the probability distribution of $N(\theta,\sigma^{2}I_{n})$ in this section.

We put a prior distribution $\Pi$ on the parameter $\theta$. Consider $\Pi$ that has the following sampling process.

Sample $w\sim\text{Beta}(\alpha_{0},\beta_{0})$;

Conditioning on $w$, sample $z_{i}\sim\text{Bernoulli}(w)$ for $i=2,3,...,n$;

Conditioning on $(z_{2},...,z_{n})$, sample $\theta_{1}\sim g$, and then for $i=2,3,...,n$, sample $\theta_{i}$ according to $\theta_{i}\sim g$ if $z_{i}=1$ and $\theta_{i}=\theta_{i-1}$ if $z_{i}=0$.

We first consider variational inference via the mean-field class, defined as

We also define $\mathcal{S}=\mathcal{S}_{\rm MF}^{\rm joint}$ on the joint distribution of $(w,z,\theta)$ by

The variational posteriors $\widehat{Q}_{\rm MF}$ and $\widehat{Q}_{\rm MF}^{\rm joint}$ are given by (3) with variational classes defined above respectively To be rigorous, the posterior distribution $\Pi(\cdot|X^{(n)})$ used in $D(Q\|\Pi(\cdot|X^{(n)}))$ are the marginal posterior of $\theta$ and the joint posterior of $(w,z,\theta)$, respectively.. Interestingly, for the piecewise constant model, both $\widehat{Q}_{\rm MF}$ and $\widehat{Q}_{\rm MF}^{\rm joint}$ give a trivial rate.

For the prior $\Pi$ specified above with any $g$ absolutely continuous with respect to the Lebesgue measure, we have

for any $k^{*}\in[n]$, where $\widehat{Q}_{\rm MF}$ and $\widehat{Q}_{\rm MF}^{\rm joint}$ are the variational posteriors defined by (3) with $\mathcal{S}=\mathcal{S}_{\rm MF}$ and $\mathcal{S}=\mathcal{S}_{\rm MF}^{\rm joint}$, respectively.

The result of Theorem 3.4 shows that the mean-field variational posteriors $\widehat{Q}_{\rm MF}$ and $\widehat{Q}_{\rm MF}^{\rm joint}$ are unable to achieve a better rate than simply estimating $\theta^{*}$ by the naive estimator $\widehat{\theta}=X$. The proof, given in Appendix B.5, reveals the reason of this phenomenon. Since the independence structure of the two classes fails to capture the underlying dependence structure of the parameter space $\Theta_{k^{*}}(B)$, the variational posterior distributions are equivalent to the posterior distribution induced by the prior $\Pi=\otimes_{i=1}^{n}g$, and therefore the condition (C4) is violated. Note that this is the first negative result in the literature on the statistical convergence of the mean-field approximation.

In order to achieve the minimax rate of the space $\Theta_{k^{*}}(B)$, it is necessary to introduce some dependence structure in the variational class. One of the simplest classes of dependent distributions is the class of first-order Markov chains, defined by

The class $\mathcal{S}_{\rm MC}$ introduces a natural dependence structure for the piecewise constant model, and it is compatible with the prior distribution $\Pi$, because conditioning on the change point pattern $z$, the prior distribution of $\theta|z$ belongs to the class $\mathcal{S}_{\rm MC}$. We also introduce a similar variational class on the joint distribution of $(w,z,\theta)$, defined by

Besides the distribution of $\theta$ restricted to $\mathcal{S}_{\rm MC}$, the distributions of $w$ and $z$ are both in the mean-field classes.

In order to derive the rates for the variational posterior distributions induced by $\mathcal{S}_{\rm MC}$ and $\mathcal{S}_{\rm MC}^{\rm joint}$, we impose the following conditions on the prior distribution $\Pi$.

There exist some constants $C_{2}>C_{1}>1$ such that

According to Theorem 2.2, we get the following result.

Consider a prior distribution $\Pi$ that satisfies (26) and (27). Then, for any $\theta^{*}\in\Theta_{k^{*}}(B)$, we have

where $\widehat{Q}_{\rm MC}$ and $\widehat{Q}_{\rm MC}^{\rm joint}$ are the variational posterior distributions defined by (3) with $\mathcal{S}=\mathcal{S}_{\rm MC}$ and $\mathcal{S}=\mathcal{S}_{\rm MC}^{\rm joint}$, respectively.

Theorem 3.5 shows that both $\widehat{Q}_{\rm MC}$ and $\widehat{Q}_{\rm MC}^{\rm joint}$ are able to achieve the minimax rate of the problem. This example illustrates the importance of the choice of the variational class. According to Theorem 2.1, the rate of a variational posterior is upper bounded by $\epsilon_{n}^{2}$, the rate of the true posterior, plus $\gamma_{n}^{2}$, the variational approximation error. The choice of $\mathcal{S}_{\rm MF}$ for the piecewise constant model leads to a very large $\gamma_{n}^{2}$, and thus a trivial rate in Theorem 3.4. On the other hand, the variational approximation errors given by the classes $\mathcal{S}_{\rm MC}$ and $\mathcal{S}_{\rm MC}^{\rm joint}$ are small, which are dominated by the minimax rate.

Though the statistical properties of the two classes $\mathcal{S}_{\rm MC}$ and $\mathcal{S}_{\rm MC}^{\rm joint}$ are both satisfactory, the class $\mathcal{S}_{\rm MC}^{\rm joint}$ enjoys a computational advantage, and the solution $\widehat{Q}_{\rm MC}^{\rm joint}$ can be computed exactly via dynamic programming. In order to characterize the solution $\widehat{Q}_{\rm MC}^{\rm joint}$, we consider the following discrete optimization problem:

The solution of (28) is denoted as the sequence $0=\widehat{a}_{0}<\widehat{a}_{1}<\cdots<\widehat{a}_{\widehat{k}}=n$. We remark that under the condition (26), the penalty term of (28) comes from the fact that

Let the maximizer of (28) be $(\widehat{a}_{0},\widehat{a}_{1},...,\widehat{a}_{\widehat{k}})$. For

the distributions $\widehat{Q}^{(w)}$, $\widehat{Q}^{(z)}$ and $\widehat{Q}^{(\theta)}$ are specified as follows.

Under $\widehat{Q}^{(z)}$, $z_{\widehat{a}_{j}+1}=1$ for $j=1,...,\widehat{k}-1$, and $z_{i}=0$ elsewhere with probability $1$.

We have $\widehat{Q}^{(w)}=\text{Beta}(\widehat{k}+\alpha_{0}-1,n-\widehat{k}+\beta_{0})$.

We have $d\widehat{Q}^{(\theta)}(\theta)=d\widehat{Q}_{1}^{(\theta)}(\theta_{1})\prod_{i=2}^{n}d\widehat{Q}_{i}^{(\theta)}(\theta_{i}|\theta_{i-1})$, where

By Theorem 3.6, in order to get $\widehat{Q}_{\rm MC}^{\rm joint}$, it is sufficient to solve (28). This can be done through a dynamic programming given in Algorithm 1. To simplify the notation, we define

We note that the computational cost of the dynamic programming above is $O(n^{3})$ (see ), and for any integers $0\leq a<b\leq n$, (29) has a closed form as long as we use a conjugate $g(\cdot)$.

Variational Bayes with Model Selection

In this section, we consider a general form of probability models

Here, the probability $P_{k,\theta^{(k)}}^{(n)}$ is determined by an index $k$ and a parameter $\theta^{(k)}$. We assume that the set $\mathcal{K}$ is either countable or finite. For a given $k$, the probability $P_{k,\theta^{(k)}}^{(n)}$ is parametrized by a $\theta^{(k)}$ in a parameter space $\Theta^{(k)}$ that is indexed by this $k$. Without loss of generality, we assume that the parameter $\theta^{(k)}$ can be written in a blockwise structure

Note that the dimension of $\theta^{(k)}$ may vary with $k$.

The model $\mathcal{M}$ is very natural for many applications. One can think of $k$ as a model dimension index, which determines the complexity of the parameter space $\Theta^{(k)}$. A leading example is the mixture density model, where $k$ stands for the number of components.

To model the hierarchical structure of $(k,\theta^{(k)})$, one naturally uses a hierarchical prior distribution, which is specified through the following sampling process:

Firstly, sample $k\sim\pi$ from $\mathcal{K}$;

Conditioning on $k$, sample $\theta^{(k)}$ from the probability measure $\Pi^{(k)}$, and $\Pi^{(k)}$ has a product structure

For variational inference, we consider a mean-field class that naturally takes advantage of the structure of the prior distribution. For a given $k\in\mathcal{K}$, the corresponding mean-field class is defined as

In order to select the best model from the data, we consider optimizing the evidence lower bound (ELBO). With the notation $p(X^{(n)}|\theta^{(k)})$ standing for the joint likelihood function, the marginal likelihood given a model $k\in\mathcal{K}$ is defined by

Then, a straightforward model selection procedure is to maximize $\log\left(p(X^{(n)}|k)\pi(k)\right)$ over $k\in\mathcal{K}$. In order to overcome the intractability of the integral (32), we instead optimize a lower bound, which is given by

which can be derived by a direct application of Jensen’s inequality. Denote the right hand side of (33) by $F(Q^{(k)},k)$, and we will solve the following optimization problem,

Finally, the solution to (34) leads to the variational posterior distribution $\widehat{Q}=\widehat{Q}^{(\widehat{k})}$ that we use in a model selection context. A similar variational approximation to the tempered posterior in the model selection setting was studied by .

2 Convergence Rates

Assume the observation $X^{(n)}$ is generated from a probability measure $P_{0}^{(n)}$, and $\widehat{Q}=\widehat{Q}^{(\widehat{k})}$ is the variational posterior that is a solution to (34). For the general settings described above, we show that the variational approximation error can be automatically controlled by a prior mass condition. Let $\Pi$ be the prior distribution on $P_{k,\theta^{(k)}}$ induced by the sampling process of $(k,\theta^{(k)})$.

Suppose $\epsilon_{n}$ is a sequence that satisfies $n\epsilon_{n}^{2}\geq 1$. Let $\rho>1$ be a constant and $C_{2},C_{3}>0$ be constants. We assume that there exists a $k_{0}\in\mathcal{K}$ and a subset $\Theta^{(k_{0})}=\otimes_{j=1}^{m_{k_{0}}}\Theta_{j}^{(k_{0})}\subset\left\{\theta^{(k_{0})}:D_{\rho}\left(P_{0}^{(n)}\|P_{k_{0},\theta^{(k_{0})}}^{(n)}\right)\leq C_{3}n\epsilon_{n}^{2}\right\}$, such that

where $\pi(k_{0})$ and $\Pi_{j}^{(k_{0})}$ are defined in the prior sampling procedure. Moreover, assume that the conditions (C1) and (C2) hold for all $\epsilon>\epsilon_{n}$ with respect to prior procedure $\Pi$ and some constant $C>C_{2}+C_{3}+2$. Then for the variational posterior $\widehat{Q}^{(\widehat{k})}$ defined as the solution of (34), we have

Theorem 4.1 characterizes the convergence rate of mean-field variational posterior with model selection using the conditions (C1), (C2) and (C3*). Given the structure of the prior distribution, an equivalent way of writing (C3*) is

for the factorized structure of $\Theta^{(k_{0})}$. Therefore, our three conditions (C1), (C2) and (C3*) still fall into the “prior mass and testing” framework, and directly correspond to the three conditions in for convergence rates of the true posterior.

An interesting special case is when the set $\mathcal{K}$ is a singleton. Then, for a product prior measure and the mean-field variational class, the condition (C3*) is reduced to (2) discussed in Section 1.

3 Density Estimation via Location-Scale Mixtures

In this section, we consider the location-scale mixture model as an application of the theory. The location-scale mixture density is defined as

for some positive even integer $p$. The kernel $\psi_{\sigma}(\cdot)$ has a pre-specified form, for example, Gaussian density when $p=2$, while the parameters $k$ and $\theta^{(k)}=(w,\mu,\sigma)$ are to be learned from the data.

Given i.i.d. observations $X_{1},...,X_{n}$ sampled from some density function $f_{0}$, our goal is to estimate the density $f_{0}$ through the location-scale mixture model (36). We denote the probability distribution of the mixture density $p(x|k,\theta^{(k)})$ as $P_{k,\theta^{(k)}}$ and a probability distribution with a general density $f$ as $P_{f}$. In the paper , a Bayesian procedure is proposed and a nearly minimax optimal convergence rate is derived for the true posterior distribution. We will follow the same setting in , but analyze the variational posterior.

We first specify the prior distribution $\Pi$ through the following sampling process:

Sample the number of mixtures $k\sim\pi$;

Conditioning on $k$, sample the location parameters $\mu_{1},\cdots,\mu_{k}$ independently from $p_{\mu}$, sample the weights $w=(w_{1},\cdots,w_{k})$ from $p_{w}^{(k)}$, and then sample the precision parameter $\tau=\sigma^{-2}$ from $p_{\tau}$.

In order to optimize (34) in the variational Bayes framework, we specify the blockwise structure (31) in this case as

Note that we do not factorize $dQ_{w}^{(k)}(w)$ because of the constraint $\sum_{j=1}^{k}w_{j}=1$. The variational posterior distribution is defined as $\widehat{Q}=\widehat{Q}^{(\widehat{k})}$ that solves (34). The loss function here is chosen as $n$ times squared Hellinger distance, i.e., $L(P_{f}^{n},P_{f_{0}}^{n})=nH^{2}(P_{f},P_{f_{0}})$.

In order that $\widehat{Q}$ enjoys a good convergence rate, we need conditions on the prior distribution and the true density function $f_{0}$. We first list the conditions on the prior.

There exist constants $C_{1},C_{2}>0$, such that

for all $m>0$. There exist constants $t,C_{3},C_{4}>0$, such that

for all $n^{\frac{1}{2\alpha+1}}\leq k_{0}\leq n^{\frac{1}{2\alpha+1}+t}$.

There exist constants $c_{1},c_{2},c_{3}>0$, such that

for all $x_{0}>0$ and constants $c_{4},c_{5},c_{6}$, such that

There exist constants $t,d_{1},d_{2},d_{3}>0$, such that

for all $w_{0}\in\Delta_{k_{0}}$ and $n^{\frac{1}{2\alpha+1}}\leq k_{0}\leq n^{\frac{1}{2\alpha+1}+t}$, where $\Delta_{k_{0}}(w_{0},\epsilon)=\{w\in\Delta_{k_{0}}:\|w-w_{0}\|_{1}\leq\epsilon\}$.

There exist constants $b_{0},b_{1},b_{2},b_{3}>0$, such that

for all $\tau_{0}>0$. There exist constants $b_{4},b_{5}>0$ and a constant $b_{6}\in(0,1]$ that satisfy

The conditions on the prior distribution are quite general. For example, one can choose $k\sim\text{Poisson}(\xi_{0})$, $\mu_{j}\sim N(0,\sigma_{0}^{2})$, $w\sim\rm{Dir}(\alpha_{0},\alpha_{0},\cdots,\alpha_{0})$ and $\tau\sim\Gamma(a_{0},b_{0})$ for some positive constants $\xi_{0},\sigma_{0},\alpha_{0},a_{0},b_{0}$. Then, the conditions above are all satisfied.

Next, we list the conditions on the true density function $f_{0}$:

(Smoothness) The logarithmic density function $\log f_{0}$ is assumed to be locally $\alpha$-Hölder smooth. In other words, for the derivative $l_{j}(x)=\frac{d^{j}}{dx^{j}}\log f_{0}(x)$, there exists a polynomial $L(\cdot)$ and a constant $\gamma>0$ such that,

for all $x,y$ that satisfies $|x-y|\leq\gamma$. Here, the degree and the coefficients of the polynomial $L(\cdot)$ are all assumed to be constants. Moreover, the derivative $l_{j}(x)$ satisfies the bound $\int|l_{j}(x)|^{\frac{2\alpha+\epsilon}{j}}f_{0}(x)dx<s_{\max}$ for all $j=1,...,{\left\lfloor{\alpha}\right\rfloor}$ with some constants $\epsilon,s_{\max}>0$.

(Tail) There exist positive constants $T$, $\xi_{1}$, $\xi_{2}$, $\xi_{3}$ such that

(Monotonicity) There exist constants $x_{m}<x_{M}$ such that $f_{0}$ is nondecreasing on $(-\infty,x_{m})$ and is nonincreasing on $(x_{M},\infty)$. Without loss of generality, we assume $f_{0}(x_{m})=f_{0}(x_{M})=c$ and $f_{0}(x)\geq c$ for all $x_{m}<x<x_{M}$ with some constant $c>0$.

These conditions are exactly the same as in and similar conditions are also considered in . The conditions allow a well-behaved approximation to the true density by a location-scale mixture. There are many density functions that satisfy the conditions (B1)-(B3), for which we refer to .

The convergence rate of the variational posterior is given by the following theorem.

Consider i.i.d. observations generated by $P_{f_{0}}^{n}$, and the density function $f_{0}$ satisfies conditions (B1)-(B3). For the prior that satisfies (40)-(46), we have

where $\widehat{Q}=\widehat{Q}^{(\widehat{k})}$ is the solution of (34), and $r=\frac{p}{\min\{p,\xi_{3}\}}+\max\{d_{3}+1,\frac{c_{6}}{\min\{p,\xi_{3}\}}\}$, with $p,\xi_{3},c_{6},d_{3}$ defined in (37), (48), (43) and (44), respectively.

The proof of Theorem 4.2 largely follows the arguments in that are used to establish the corresponding result for the true posterior distribution, thanks to the fact that Theorem 4.1 requires three very similar “prior mass and testing” conditions to that of . The only difference is that function approximations via location-scale mixtures need to be analyzed under a stronger divergence $D_{\rho}(\cdot\|\cdot)$ for some $\rho>1$. For this reason, the proof of Theorem 4.2 relies on the construction of a surrogate density function $\widetilde{f}_{0}$. We first apply Theorem 4.1 and establish a convergence rate under $\widetilde{f}_{0}$. Then, the conclusion is transferred to $f_{0}$ with a change-of-measure argument. Details of the proof are given in Appendix B.6.

4 Dealing with Latent Variables

For the mixture model considered in Section 4.3, we discuss a variation of the variational Bayes approach (34) by including latent variables. This facilitates computation and leads to a simple coordinate ascent algorithm that has closed-form updates. In the setting of mixture model, our approach is adaptive to the unknown number of components, and can be regarded as an extension of for variational inference with latent variables.

Since $p(X^{(n)}|k,\theta^{(k)})=\prod_{i=1}^{n}\sum_{j=1}^{k}w_{j}\psi_{\sigma}(X_{i}-\mu_{j})$ with $\theta^{(k)}=(\mu,w,\sigma)$, we can write

where $p(X^{(n)}|z^{(k)},\theta^{(k)})=\prod_{i=1}^{n}\prod_{j=1}^{k}\psi_{\sigma}(X_{i}-\mu_{j})^{\mathbf{1}{\left\{z_{i}^{(k)}=j\right\}}}$, and the probability of $z_{i}^{(k)}=j$ is $w_{j}$ under $w^{(k)}(\cdot)$. We use the notation $\bar{\Pi}^{(k)}$ for the joint distribution of $(z^{(k)},\theta^{(k)})$, and then the marginal likelihood (32) can be written as

Similar to (33), the evidence lower bound with the latent variables is given by

The right hand side of (49) is shorthanded by $\bar{F}(\bar{Q}^{(k)},k)$. Define

Then, we solve the following optimization problem,

The solution to (50) leads to the variational posterior distribution $\widehat{Q}=\widehat{Q}^{(\widehat{k})}_{\rm latent}$. It is worth noting that even though $\widehat{Q}$ is a joint distribution of $(z,\mu,w,\sigma)$, the posterior inference only relies on the marginal of $(\mu,w,\sigma)$, since the parametrization of the density $f(\cdot)$ in (36) does not depend on the latent variables. The existence of the latent variables only facilitates computation.

Consider i.i.d. observations generated by $P_{f_{0}}^{n}$, and the density function $f_{0}$ satisfies conditions (B1)-(B3). For the prior that satisfies (40)-(46), we have

where $\widehat{Q}=\widehat{Q}^{(\widehat{k})}_{\rm latent}$ is the solution to (50), and $r=\frac{p}{\min\{p,\xi_{3}\}}+\max\{d_{3}+1,\frac{c_{6}}{\min\{p,\xi_{3}\}}\}$, with $p,\xi_{3},c_{6},d_{3}$ defined in (37), (48), (43) and (44), respectively.

Theorem 4.3 shows that the variational posterior with latent variables achieves the same contraction rate as in Theorem 4.2. In fact, the two variational lower bounds (33) and (49) satisfy the following relation,

which implies that the introduction of latent variables makes the variational approximation looser. On the other hand, Theorem 4.3 shows that the worse variational approximation does not compromise the statistical convergence rate. Moreover, with the help of latent variables, $\widehat{Q}^{(\widehat{k})}_{\rm latent}$ can be computed via standard variational inference algorithms. Details of the computational issues are given in Appendix A.2.

Discussion

Conditioning on $k$, sample $\theta_{j}\sim f_{j1}$ for all $j\in[k]$, and sample $\theta_{j}\sim f_{j2}$ for all $j>k$.

An empirical Bayes procedure maximizes $e^{m_{k}(X^{(n)})}\pi(k)$ The canonical form of empirical Bayes has a flat prior on $k$., where

is the logarithm of marginal likelihood. With the maximizer $\widehat{k}$, the empirical Bayes posterior is defined as

Compared with a hierarchical Bayes approach, the empirical Bayes procedure does not need to evaluate the posterior distribution of $k$, and thus in many cases is easier to implement.

We also study mean-field approximation of the posterior distribution. In order to characterize its form, we need a few definitions. For any $g=(g_{j})_{j=1}^{\infty}$, define

for any $g$. We also define the density classes $\mathcal{G}_{j1}=\left\{g\geq 0:\int g=\int_{\Theta_{j1}}g=1\right\}$ and $\mathcal{G}_{j2}=\left\{g\geq 0:\int g=\int_{\Theta_{j2}}g=1\right\}$. The next theorem gives the exact form of the mean-field variational posterior.

Consider the variational posterior $\widehat{Q}_{\rm VB}$ induced by the sieve prior and the mean-field variational set $\mathcal{S}_{\rm MF}$. The distribution $\widehat{Q}_{\rm VB}$ is a product measure with the density of each coordinate specified by

where for each given $k$, $(\widetilde{g}_{j1}^{(k)})_{j=1}^{k}$ and $(\widetilde{g}_{j2}^{(k)})_{j=k}^{\infty}$ maximize the following objective function,

under the constraints that $g_{j1}\in\mathcal{G}_{j1}$ and $g_{j2}\in\mathcal{G}_{j2}$ for all $j$, $\widetilde{k}$ maximizes

The result of Theorem 5.1 also applies to the class $\mathcal{S}_{\rm G}$ discussed in Section 3.2 with $\mathcal{G}_{j1}$ replaced by the Gaussian class. We note that Theorem 5.1 can be viewed as an extension of Theorem 3.1. In fact, if the likelihood function can be factorized over each coordinate of $\theta$, the form of $\widehat{Q}_{\rm VB}$ can be greatly simplified.

Under the same setting of Theorem 5.1, if we further assume that $p(X^{(n)}|\theta)=\prod_{j=1}^{\infty}p(X_{j}^{(n)}|\theta_{j})$, then we will have

In light of Theorem 5.1, we can compare the variational Bayes approach and the empirical Bayes approach, especially the definitions of $\widetilde{k}$ and $\widehat{k}$. The empirical Bayes chooses the best model by maximizing $e^{m_{k}(X^{(n)})}\pi(k)$, or equivalently $\pi(k|X^{(n)})$, while the variational Bayes maximizes (54). There are two major differences. The first difference is that empirical Bayes uses the exact marginal likelihood function $m_{k}(X^{(n)})$ and variational Bayes uses a mean-field approximation of $m_{k}(X^{(n)})$. We remark that in the case of likelihood that can be factorized, the mean-field approximation is exact, which leads to (55). The second difference is that empirical Bayes maximizes the posterior probability of the $k$th model, but the variational Bayes maximizes the sum of the posterior probabilities (or their mean-field approximations) of the $(k-1)$th and the $k$th models.

Despite the two differences, the empirical Bayes approach and the variational Bayes approach have a lot in common. Both are random probability distributions that summarize the information in data and prior. Both select a sub-model according to very similar criteria. To close this section, we show that with a special variational class, the induced variational posterior is exactly the empirical Bayes posterior.

Then, the empirical Bayes posterior $\widehat{Q}_{\rm EB}$ defined by (51) is the variational posterior induced by the sieve prior and the variational class $\mathcal{S}_{\rm EB}$.

The result of Theorem 5.2 shows that for sieve priors, one can view the empirical Bayes approach as a variational Bayes approach, which suggests that it may be possible to unify the theoretical analysis in this paper and the analysis of empirical Bayes procedures in .

2 Variational Approximation as Regularization

According to Theorem 2.1, the convergence rate of the posterior is determined by the sum of $\epsilon_{n}^{2}$, the rate of the true posterior, and $\gamma_{n}^{2}$, the variational approximation error. Since $\epsilon_{n}^{2}+\gamma_{n}^{2}\geq\epsilon_{n}^{2}$, it seems that the convergence rate of variational posterior is always no faster than that of the true posterior. However, Theorem 2.1 just gives an upper bound. In this section, we give two examples, and we show that it is possible for a variational posterior to have a faster convergence rate than that of the true posterior.

We consider the setting of Gaussian sequence model (10). The true signal $\theta^{*}$ that generates the data is assumed to belong to the Sobolev ball $\Theta_{\alpha}(B)$. The prior distribution is specified as

Note that a similar Gaussian process prior is well studied in the literature . We force all the coordinates after $n$ to be zero, so that the variational approximation through Kullback-Leibler divergence will not explode. For the specified prior, the posterior contraction rate is $n^{-\frac{2(\alpha\wedge\beta)}{2\beta+1}}$, and when $\beta=\alpha$, the optimal minimax rate $n^{-\frac{2\alpha}{2\alpha+1}}$ is achieved.

for a given integer $k$. It is easy to see that the variational posterior $\widehat{Q}_{[k]}$ defined by (3) with $\mathcal{S}=\mathcal{S}_{[k]}$ can be written as

In other words, the class $\mathcal{S}_{[k]}$ does not put any constraint on the first $k$ coordinates and shrink all the coordinates after $k$ to zero. Ideally, one would like to use $\delta_{0}$ for the coordinates after $k$. However, that would lead to $D\left(Q\|\Pi(\cdot|Y)\right)=\infty$ for all $Q\in\mathcal{S}_{[k]}$ given that the support of $\delta_{0}$ is a singleton. That is why we use $N(0,e^{-jn})$ instead. The rate of $\widehat{Q}_{[k]}$ for each $k$ is given by the following theorem.

For the variational posterior $\widehat{Q}_{[k]}$, we have

where $\widehat{Q}_{[k]}$ is the variational posterior defined by (3) with $\mathcal{S}=\mathcal{S}_{[k]}$.

Note that Theorem 5.3 gives both upper and lower bounds for $\widehat{Q}_{[k]}$. This makes the comparison between variational posterior and true posterior possible. Observe that when $k=\infty$, we have $\widehat{Q}_{[\infty]}=\Pi(\cdot|Y)$, and the result is reduced to the posterior contraction rate $n^{-\frac{2(\alpha\wedge\beta)}{2\beta+1}}$ in .

Depending on the values of $\alpha,\beta$ and $k$, the rate for $\widehat{Q}_{[k]}$ can be better than that of the true posterior. For example, when $\beta<\alpha$, the choice $k=n^{\frac{1}{2\alpha+1}}$ leads to the minimax rate $n^{-\frac{2\alpha}{2\alpha+1}}$, which is always faster than $n^{-\frac{2(\alpha\wedge\beta)}{2\beta+1}}$. This is because for a $\beta<\alpha$, the true posterior distribution undersmooths the data, but the variational class $\mathcal{S}_{[k]}$ with $k=n^{\frac{1}{2\alpha+1}}$ helps to reduce the extra variance resulted from undersmoothing by thresholding all the coordinates after $k$. On the other hand, when $\beta\geq\alpha$, an improvement through the variational class $\mathcal{S}_{[k]}$ is not possible. In this case, the true posterior has already overly smoothed the data, and the information loss cannot be recovered by the variational class.