On Statistical Optimality of Variational Bayes
Debdeep Pati Anirban Bhattacharya Yun Yang
math.ST stat.ML stat.TH






1]Debdeep Patidebdeep@stat.tamu.edu
1]Anirban Bhattacharyaanirbanb@stat.tamu.edu
2]Yun Yang yyang@stat.fsu.edu
[2]Department of Statistics, Florida State University
[1]Department of Statistics, Texas A&M University

On Statistical Optimality of Variational Bayes[To appear in AISTATS 2018]
    [
    

=========================================================================




The article  addresses a long-standing open problem on the justification of using variational Bayes methods for parameter estimation.  We provide general conditions for obtaining optimal risk bounds for point estimates acquired from mean-field variational Bayesian inference. The conditions pertain to the existence of certain test functions for the distance metric on the parameter space and minimal assumptions on the prior.  A general recipe for verification of the conditions is outlined which is  broadly applicable to existing Bayesian models with or without latent variables. As illustrations, 
specific applications to Latent Dirichlet Allocation and Gaussian mixture models are discussed. 

Keywords: Bayes risk; Evidence lower bound; Latent variable models; optimality; Variational inference.




§ INTRODUCTION


Variational inference <cit.> is now a well-established tool to approximate intractable posterior distributions in hierarchical multi-layered Bayesian models. The traditional Markov chain Monte Carlo (MCMC; <cit.>) approach of approximating distributions with intractable normalizing constants draws (correlated) samples according to a discrete-time Markov chain whose stationary distribution is the target distribution. Despite their success and popularity, MCMC methods can be slow to converge and lack scalability in big data problems and/or problems involving very many latent variables, which has fueled search for alternatives. 

In contrast to the sampling approach of MCMC, variational inference approaches the problem from an optimization viewpoint. First, a class of analytically tractable distributions, referred to as the variational family, is identified for the problem at hand. For example, in mean-field approximation, the set of parameters and latent variables is divided into blocks and the variational distribution is assumed to be independent across blocks. The distribution in the variational family closest to the target distribution relative to the Kullback–Leibler (KL) divergence is then used as a proxy to the target. Implementation-wise, the above optimization is commonly solved using coordinate descent or alternating minimization. A comprehensive review of various aspects of variational inference can be found in the recent article by <cit.>.

Variational inference has arguably found its most potent applications in latent variable models such as mixture models, hidden Markov models, graphical models, topic models, and neural networks; see <cit.> for a flavor of this enormous literature. Due to the fast convergence properties of the variational objective, variational inference algorithms are typically orders of magnitude faster in big data problems compared to MCMC approaches <cit.>. However, in spite of their tremendous empirical success, a general statistical theory qualifying the statistical properties of a variational solution is lacking. Existing results operate in a case-by-case manner, often directly analyzing the specific iterative algorithm to characterize properties of its limit; see Section 5.2 of <cit.> for a comprehensive review. These analyses typically require sufficient tractability of the successive steps of the iterative algorithm, and can be difficult to adapt to minor changes in the prior. More recently <cit.> modified the objective function to introduce an inverse-temperature parameter, and obtained general guarantees for the variational solution under this modified objective function.

In this article, we aim to address the general question as to whether point estimates obtained from usual variational approximations share the same statistical accuracy as those from the actual posterior. We clarify at the onset that we operate in a frequentist setting assuming the existence of a true data generating parameter. Our novel contribution is to relate the Bayes risk relative to the variational solution for a general distance metric (defined on the parameter space) to (i) the existence of certain test functions for testing the true parameter against complements of its neighborhood under the error metric, and (ii) the size of the variational objective function. As an important consequence of such Bayes risk bounds, risk bounds for variational point estimates can be readily derived when the distance metric is convex. If the risk of the variational point estimate coincides with the contraction rate of the exact posterior at the true parameter, it can be argued that there is no loss of statistical efficiency, at least asymptotically, in using the variational approximation. 

We identify a number of popularly used models where the conditions can be satisfied and the variational point estimate attains the minimax rate. Since variational Bayes is primarily used for point estimation, our theory suggests that variational Bayes successfully achieves its desiderata. As vignettes, we present two non-trivial examples in the form of density estimation in Latent Dirichlet Allocation (LDA) for topic modeling, and estimating component specific parameters in Gaussian mixture models. 






§ BACKGROUND


In this section, we introduce notations and offer some background to set up our theoretical results. Let h(p || q) =[∫(√(p) - √(q))^2dμ]^1/2 and D(p || q) = ∫ plog (p/q)dμ denote the Hellinger distance and Kullback–Leibler divergence, respectively, between two probability density functions p and q relative to a a common dominating measure μ. Also denote V(p || q) = ∫ p log^2(p/q) dμ. For a set A, we use I_A to denote its indicator function. For any vector μ and positive semidefinite matrix Σ, we use N(μ, Σ) to denote the normal distribution with mean μ and covariance matrix Σ, and use N(θ; μ, Σ) to denote its pdf at θ. We use w.h.p. to abbreviate “with high probability", when the probability is evident from the context. Throughout, C denotes a constant independent of everything else whose value may change from one line to the other. We write a ≲ b to denote a ≤ C b for some constant C> 0. Similarly, a ≳ b. 

Suppose we have n observations Y_1, …, Y_n and a probabilistic model P_θ^(n) for the joint distribution of the data Y^n=(Y_1,…,Y_n), with a density p_θ^(n) relative to the Lebesgue measure. Here, θ is the unknown parameter to be estimated from the data which lives in some parameter space Θ. Our formulation does not require the Y_is to be identically distributed or even independent. In a Bayesian setup, uncertainty regarding the parameter is quantified through a prior distribution p_θ on Θ, which upon observing the data is updated to the posterior distribution using Bayes' theorem: 

    p(θ∈ B | Y^n) = ∫_B[p(Y^n | θ)]  p_θ(θ)  dθ/∫_Θ[p(Y^n | θ)]  p_θ(θ)  dθ,

for any measurable subset B ⊂Θ and p(Y^n | θ) := p_θ(Y^n).  

In a wide variety of practical problems, the likelihood function p(Y^n  | θ) may be intractable or difficult to analyze directly. For example, in a 2-component Gaussian mixture model, p(Y^n  | θ) is a combinatorial sum over 2^n terms. A widely used trick in such situations is to introduce latent variables to simplify the conditional likelihood. Specifically, assume one can decompose

    p(Y^n  | θ) = ∑_s^n p(Y^n  |  S^n = s^n, μ)  π_s^n,

where S^n = (S_1, …, S_n) denotes a collection of discrete latent variables, with S_i∈{1,2,…,K} the latent variable for the ith observation. We have assumed the parameter θ can be decomposed as θ = (μ, π), with p(Y^n  |  S^n = s^n, θ) = p(Y^n | S^n = s^n, μ) and 
and π_s^n : = p(S^n = s^n  | θ) denotes the probability of the latent vector taking on the value s^n. In the 2-component mixture model example, S_i ∈{1, 2} denotes the latent membership indicator for the ith observation.  We assume discrete latent variables for notational convenience and note that our results generalize to continuous latent variables in a straightforward fashion. 

Let Z^n = (θ, S^n). The augmented posterior p(Z^n  |  Y^n) assumes the form 

    p(θ, s^n  |  Y^n) ∝ p(Y^n  | μ, s^n)  π_s^n  p_θ(θ),

whose constituent terms are typically more tractable. 

Variational inference in this setup proceeds by first identifying a variational family comprising of distributions Γ on Z^n and finding the closest member in this family to p(Z^n | Y^n) relative to the KL divergence 

    _Z^n   :  = _q_Z^n∈Γ D[ q_Z^n(·) ||  p(· |  Y^n) ] = _q_Z^n∈Γ L(q_Z^n),

where 

    L(q_Z^n) = ∫ q_Z^n(z^n) logp(Y^n | z^n)  p_Z^n(z^n)/q_Z^n(z^n)  dz^n

is the evidence lower bound (ELBO) which gives a lower bound to the log marginal likelihood log p(Y^n). If Γ is completely unrestricted, then _Z^n coincides with the posterior distribution p(Z^n | Y^n). In practice, the choice of the variational family Γ is dictated by a trade-off between flexibility and computational tractability. For example, in mean-field variational approximation, it is common to assume independence among the parameters and the latent variables in the variational family, whence q_Z^n decomposes as 

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

We shall assume the mean-field decomposition (<ref>) throughout, so that the minimizer _Z^n in (<ref>) necessarily is of the form 
_Z^n = _θ⊗_S^n.
In many situations, the constituent terms may be further decomposed as q_S^n(s^n) = ⊗_i=1^n q_S_i(s_i) and q_θ(θ) = q_μ(μ) ⊗ q_π(π). Such a decomposition may be either due to computational reasons or implied by the conditional independence structure of the model.

Under the mean-field decomposition (<ref>) and using (<ref>), we have, after some simplification, 

    D[ q_Z^n(·) ||  p(· |  Y^n) ] = D(q_θ || p_θ) 
        - ∫_Θ[ ∑_s^n q_S^n(s^n)  logp(Y^n | μ, s^n)  π_s^n/q_S^n(s^n)] _ℓ_n(θ)  q_θ(dθ).

The quantity ℓ_n(θ) is an approximation to the log likelihood ℓ_n(θ) : = log p(Y^n |θ) in terms of the latent variables. To see this, multiply and divide the right hand side of (<ref>) by q_S^n(s^n) and apply Jensen's inequality to the concave function x ↦log x to conclude that ℓ_n(θ) ≥ℓ_n(θ). Replacing ℓ_n(θ) by ℓ_n(θ) in (<ref>) and adjusting the error term, we obtain

    D[ q_Z^n(·) ||  p(· |  Y^n) ] = D(q_θ || p_θ) - ∫_Θℓ_n(θ) q_θ(dθ) 
        + ∫_Θ[ ℓ_n(θ) - ℓ_n(θ) ]    q_θ(dθ)_Δ(q_θ, q_S^n).

The quantity Δ is an average error due to the likelihood approximation and is clearly nonnegative. In the specific situation where no latent variables are present in the model, Δ≡ 0. However, in general, Δ is a strictly positive quantity. 



§ VARIATIONAL RISK BOUNDS


We are now prepared to state a general Bayes risk bound for the variational distribution. We shall operate in a frequentist framework and assume the existence of a true data generating parameter θ^∗∈Θ. In other words, we assume that the data Y^n is distributed according to P_θ^∗^(n). Let d(·, ·) be a distance metric on the parameter space Θ which quantifies distance between two putative parameter values. For example, if Θ designates a space of densities so that each θ can be identified with a density function, d can be chosen as the Hellinger or total variation distance. We are interested in obtaining bounds on 

    ∫ d^κ(θ, θ^∗)  _θ(θ)

for some κ > 0, that hold with high probability with respect to P_θ^∗^(n). In particular, if d^κ is convex in its first argument, then by Jensen's inequality, 

    d^κ(θ_VB, θ^∗) ≤∫ d^κ(θ, θ^∗)  _θ(θ),
 where θ_VB = ∫θ _θ(dθ) is the mean of the variational distribution and a surrogate for the posterior mean. We are specifically interested in obtaining sufficient conditions for the variational point estimate θ_VB to contract at the same rate as posterior mean. Since variational approaches are overwhelmingly used for rapidly obtaining point estimates, such a result will indicate that at least in terms of rates of convergence, there is no loss of statistical accuracy in using a variational approximation. Moreover, the negative result from <cit.> shows that the spread of the variational distribution _θ is typically “too small" compared with that of the sampling distribution of the maximum likelihood estimator. This fact combined with the Bernstein von-Mises theorem (Chapter 10 of <cit.>) implies the inadequacy of using _θ for approximating the true posterior distribution, and a rate optimal variational point estimator is the best one can hope for in general.

Define ℓ_n(θ, θ^∗) :  = ℓ_n(θ) - ℓ_n(θ^∗) to be the log-likelihood ratio between θ and θ^∗. We can replace ℓ_n(θ) with ℓ_n(θ, θ^∗) inside the integrand in (<ref>) without affecting the minimization problem - this is done for purely theoretical reasons to harness the structure of the log-likelihood ratio. Let us call the equivalent objective function Ω, so that 

    Ω(q_θ, q_S^n) : = 
        - ∫_Θℓ_n(θ, θ^∗)   q_θ(dθ) + Δ(q_θ, q_S^n) + D(q_θ || p_θ).


We are now ready to state the main assumption. 

Assumption T:  (existence of tests) Let ε_n be a sequence satisfying ε_n → 0 and n ε_n^κ→∞ for some κ≥ 2. Let ϕ_n ≡ϕ_n(Y^n) ∈ (0, 1) be a sequence of test functions for testing 

    H_0: θ = θ^∗ versus.  H_1: d(θ, θ^∗)  > ε_n

with type-I and II error rates satisfying

    E_θ^∗ [ϕ_n] ≤ e^-2 n ε_n^κ, 
        E_θ [1 - ϕ_n] ≤ e^-C n d^κ(θ, θ^∗),

for any θ∈Θ with d(θ, θ^∗) > ε_n, where E_θ denotes an expectation with respect to P_θ^(n). While κ = 2 appears naturally in most problems, we provide a non-standard example (estimating component specific means in a mixture model) in Section 4 with κ = 4. 

In models satisfying the monotone likelihood ratio property <cit.> such
as exponential families, one can construct such tests (with κ = 2) from the generalized likelihood ratio test (GLRT) statistic when d corresponds to the Euclidean metric on the natural parameter. 
A general recipe <cit.> to construct such tests when Θ is compact relative to d is to (i) construct an ε_n/2-net N = {θ_1, …,  θ_N} such that for any θ with d(θ, θ^∗) > ε_n, there exists θ_j ∈ N with d(θ, θ_j) < ε_n/2, (ii) construct a test ϕ_n,j for H_0: θ = θ^∗ versus H_1: θ = θ_j with type-I and II error rates as in Assumption T, and (iii) set ϕ_n = max_1≤ j ≤ Nϕ_n, j. 
The type-II error of ϕ_n retains the same upper bound, while the type-I error can be bounded by N   e^-2 n ε_n^κ. Since N can be further bounded by N(Θ, ε_n/2, d), the covering number of Θ by d-balls of radius ε_n/2, it suffices to show that N(Θ, ε_n/2, d) ≲ e^n ε_n^κ. For example, when Θ is a compact subset of ℝ^d, then N(Θ, ε_n/2, d) ≲ε_n^-d≲ e^n ε_n^2 as long as ε_n ≲√(log n/n). More generally, if Θ is a space of densities and d the Hellinger/L_1 metric, then construction of the point-by-point tests in (i) (with κ = 2) from the LRT statistics follows from the classical Birgé-Lecam testing theory <cit.>; see also <cit.>. 

We are now ready to state our first theorem on the variational risk bound, which relates the Bayes risk under the variational solution to the size of the objective function Ω. 


Suppose there exists a sequence of test functions {ϕ_n} for the metric d satisfying Assumption T. Then, it holds with P_θ^∗^(n) probability at least (1 - e^- C n ε_n^κ) that 

    ∫_Θ d^κ(θ, θ^∗)  _θ(dθ) ≤Ω(q_θ, q_S^n)/n + ε_n^κ,

for any q_θ≪ p_θ and any probability distribution q_S^n on S^n which is nowhere zero. 

Theorem <ref> implies that minimizing the Bayes risk is equivalent to minimizing the objective function Ω in (<ref>). <cit.> obtained a similar result for a modified variational objective function and d being limited to Rényi divergence measures. Theorem <ref> instead allows any metric d as long as the testing condition in Assumption T can be satisfied. 

A detailed proof of Theorem <ref> (as well as proofs of other results) is provided in the Section <ref>. We sketch some of the key steps to highlight the main features of our argument. Our first key step is to show using the testing assumption T that w.h.p. (w.r.t. P_θ^∗^(n)), 

    ∫_Θ[ exp{ℓ_n(θ, θ^∗) + n d^κ(θ, θ^∗) }_ξ(θ, θ^∗)]  p_θ(dθ) ≤ e^C n ε_n^κ.

To show this, write the integral in (<ref>) as T_1 + T_2, where T_1 and T_2 respectively split the integral  over {d(θ, θ^∗) ≤ε_n} and {d(θ, θ^∗) > ε_n}. Using Markov's inequality along with the fact that E_θ^∗ [e^ℓ_n(θ, θ^∗)] = 1, it can be shown that T_1 ≤ e^C n ε_n^κ w.h.p. 
To tackle T_2, write T_2 = T_21 + T_22 by decomposing ξ(θ, θ^∗) = ξ(θ, θ^∗) (1 - ϕ_n) + ξ(θ, θ^∗) ϕ_n, where ϕ_n is the test from Assumption T. Using Markov's inequality and the fact that ϕ_n has small type-II error, it can be shown that T_21≤ e^-C n ε_n^κ w.h.p. The bound on the type-I error of ϕ_n along with Markov's inequality yields ϕ_n ≤ e^-2 n ε_n^κ w.h.p., which yields T_22≤ e^-2 n ε_n^κ T_2 w.h.p. Combining, one gets T_2 ≤ e^-C n ε_n^κ w.h.p. 



Once (<ref>) is established, the next step is to link the integrand in (<ref>) with the latent variables. To that end, observe that  

    ξ(θ, θ^∗) = ∑_s^nexp{h(θ, s^n)} _S^n(s^n),

where 

    h(θ, s^n) = logp(Y^n | μ, s^n)  π_s^n/p(Y^n | θ^∗)  _S^n(s^n) + n d^κ(θ, θ^∗).

Combining the above with (<ref>), we have, w.h.p.,

    ∫_Θ∑_s^nexp{ h(θ, s^n) } _S^n(s^n)   p_θ(d θ) ≤ e^C n ε_n^κ.

Next, use a well-known variational/dual representation of the KL divergence (see, e.g., Corollary 4.15 of <cit.>) which states that for any probability measure μ and any measurable function h with e^h ∈ L_1(μ), 

    log∫ e^h(η) μ(d η) = sup_ρ[ ∫ h(η)  ρ(d η) -  D(ρ || μ) ],

where the supremum is over all probability measures ρ≪μ. In the present context, setting η = (θ, s^n), μ : = _S^n⊗ p_θ, and ρ = _θ⊗_S^n, it follows from the variational lemma (<ref>) and some rearrangement of terms that w.h.p.

    n ∫_Θ d^κ(θ, θ^∗)  _θ(d θ) ≤ n ε_n^κ + D(_θ || p_θ)  - ∫_Θ∑_s^n h(θ, s^n)  _θ(d θ).

It then follows from (<ref>)–(<ref>) that the right hand side of the above display equals n ε_n^κ + Ω(_θ, _S^n). The proof of the theorem then follows, since by definition, Ω(_θ, _S^n) ≤Ω(q_θ, q_S^n) for any (q_θ, q_S^n) in the variational family Γ. 



Now we discuss choices of good variational distributions q_θ and q_S^n for minimizing Ω(q_θ, q_S^n), the stochastic component of the variational upper bound  in Theorem <ref>. We make some additional assumptions first on the augmented likelihood and prior in (<ref>) for the subsequent development. 
First, assume independent priors on μ and π so that p_θ = p_μ⊗ p_π. 
Next, assume p(Y^n | μ, s^n) = ∏_i=1^n p(Y_i | μ, s_i) and π_s^n = ∏_i=1^n π_s_i splits into independent components. This implies p(Y^n | θ) = ∏_i=1^n p(Y_i | θ), where p(Y_i | θ) = ∑_s_i p(Y_i | μ, s_i) π_s_i. For the variational distribution q_S^n, we additionally assume a mean-field decomposition q_S^n(s^n) =∏_i=1^n q_S_i(S_i). 


Recall that the objective function decomposes as Ω(q_θ, q_S^n) = - ∫_Θℓ_n(θ, θ^∗)   q_θ(dθ) + Δ(q_θ, q_S^n) + D(q_θ || p_θ).
The first model-fit term is an averaged (with respect to the variational distribution) log-likelihood ratio which tends to get small as the variational distribution q_θ places more mass near the true parameter θ^∗. The second term arising from the approximation of the likelihood function ℓ_n(θ) by ℓ_n(θ) will become small under some proper choice of q_S^n, as we will illustrate in the proof of Theorem <ref> below. The last regularization or penalty term prevents over-fitting to the data by constricting the KL divergence between the variational solution and the prior. Consequently, a good variational distribution q_θ should put all its mass into an appropriately small neighborhood around the truth θ^∗ so that
the first two terms in V become small; on the other hand, the neighborhood has to be large enough so that the last regularization term is not too large.

Motivated by the above discussion, we follow the development of <cit.> by defining two KL neighborhoods around (π^∗, μ^∗) with radius (ε_π, ε_μ) as

    B_n(π^∗, ε_π)    = {D(π^∗  || π) ≤ε_π^2,  V(π^∗  || π) ≤ε_π^2}, 
     B_n(μ^∗, ε_μ)    =  {sup_sD[ p(· | μ^∗,s) ||   p(· | μ,s)]≤ε_μ^2, 
           sup_s V[ p(· | μ^∗,s) ||   p(· | μ,s)] ≤ε_μ^2},

where we used the shorthand D(π^∗  || π)=∑_s π^∗_s log(π^∗_s/π_s) to denote the KL divergence between categorical distributions with parameters π^∗∈ S_K and π∈ S_K in the K-dim simplex S_K. Consistent with notation introduced at the beginning of Section 2, V(π^∗  || π)  = ∑_s π^∗_s log^2(π^∗_s/π_s). 


For any fixed (ε_π, ε_μ) ∈ (0, 1)^2, with _θ^∗ probability at least 1 - 5/{(D-1)^2  n   (ε_π^2+ε_μ^2)}, it holds that

    ∫ d^κ(θ, θ^∗) _θ(θ) dθ≤ε_n^κ + (ε_π^2+ε_μ^2) + 
       { - log P_π[ B_n(π^∗, ε_π)]/n}  + 
     { - log P_μ[B_n(μ^∗, ε_μ) ] /n}.


The proof follows Theorem 4.5 of <cit.>, and is omitted; we provide a sketch here. As discussed above, we first make a good choice of q_S^n as follows. Let _S^n be a probability distribution over S^n defined as 

    _S^n(s^n) = ∏_i=1^n _S_i(s_i)=∏_i=1^n π^∗_s_i  p(Y_i | μ^∗,  s_i)/p(Y_i | θ^∗).

Intuitively,  takes the full conditional distribution of S^n |θ, Y^n, and replaces θ by the true parameter θ^∗. With this choice, Δ(q_θ, _S^n) simplifies to 

    - ∫_Θ  q_θ(θ) ∑_i=1^n ∑_s_i_S_i(s_i) logp(Y_i | μ,s_i)  π_s_i/p(Y_i | μ^∗,s_i)  π^∗_s_i  dθ
        +  ∫_Θ q_θ(θ)  logp(Y^n | θ)/p(Y^n | θ^∗) dθ.

It now remains to choose q_θ. The first term in the above display naturally suggests choosing 
q_θ as the restriction of p_θ into B_n(π^∗, ε_π)× B_n(μ^∗, ε_μ). For this choice, the second term can also be controlled w.h.p., leading to the conclusion in Theorem <ref>. 







§ APPLICATIONS

As described in Section 3, variational risk bounds for the parameter of interest depend on the existence of appropriate test functions which characterizes the ability of the likelihood to identify the parameter.  
 
Developing test functions for studying convergence rates of estimators in classical and Bayesian statistics dates back to <cit.>, with renewed attention in the Bayesian context due to <cit.>. Specific tests have been constructed for nonparametric density estimation <cit.>, semi/non-parametric regression <cit.>, 
convergence of latent mixing measures <cit.>, high-dimensional problems <cit.>, and empirical Bayes methods <cit.>, among others. As long as the prior distributions are supported on compact subsets of the parameter space, these existing tests can be used to prove minimax optimality of the variational estimate θ̂_VB in each of these problems. We skip the details for space constraints. 




In this section, we focus on two novel examples involving latent variables where variational methods are commonly used and no theoretical guarantee is available for the variational solutions. The first one is the Latent Dirichlet Allocation (LDA; <cit.>), a generative probabilistic model for topic modeling. The second example is concerned with estimating the component specific parameters in Gaussian mixture models.






 §.§ Latent Dirichlet allocation

We first consider LDA <cit.>, a conditionally conjugate probabilistic topic model <cit.> for learning the latent “topics” contained in a collection of documents. Starting from the original paper <cit.>, the mean-field variational Bayes approximation has become a routine approach for implementing LDA. However, theoretical guarantees for the variational solution is still an open problem despite its empirical success. In this subsection, we show the rate optimality of the estimate from the mean-field approximation to LDA.

In LDA, each document is assumed to contain multiple topics, where a topic is defined as a
distribution over words in a vocabulary.
Our presentation of the model follows the notation of <cit.>. Let K be the total number of topics, V the vocabulary size, D the total number of documents, and N the number of words in each document (for simplicity, we assume the same number of words across documents). Recall that we use the notation 𝒮_d to denote the d-dim simplex. LDA contains two parameters: word distribution matrix B=[β_1, β_2, …, β_K]∈ℝ^V × K and topic distribution matrix T_γ = [γ_1, γ_2,…, γ_D]∈ℝ^K× D, where β_k∈ S_V is the word proportion vector of the kth topic, k=1,2,…,K, and γ_d∈ S_K is the topic proportion vector of the dth document, d=1,2,…,D.
Given parameters B and T_γ, the data generative model of LDA is:
    
      
  * for each document in d=1,…,D, draw a topic assignment z_dn∼(γ_d), then
      
  * for each word in n=1,…, N, draw a word w_dn∼(β_z_dn).
  
 where (π) stands for categorical distribution with probability π. 
Here, z_dn∈{1,…,K} is the latent class variable over topics so that
z_dn = k indicates the nth word in document d is assigned to the
kth topic. Similarly, w_dn∈{1,…,V} is the latent class variable over the words in the vocabulary so that w_dn = v indicates
that the nth word in document d is the vth word in the vocabulary. A common prior distribution over the parameters B and T_γ is:
  
    
  * For each topic in k=1,…,K, word proportion vector has prior β_k ∼_V(η_β),
    
  * For each document in d=1,…,D, topic proportion vector has prior γ_d ∼_K(η_γ).
  
Here, η_β∈ R_+ is a hyper-parameter of the symmetric
Dirichlet prior on the topics β, and η_γ∈ R_+^K
are hyper-parameters of the Dirichlet prior on the topic proportions for each document.
To facilitate adaptation to sparsity using Dirichlet distributions when V,  K≫1, we choose η_β=1/V^c and η_γ=1/K^c for some fixed number c>1 <cit.>.

To apply our theory, we view N as the sample size, and D as the “dimension" of the parameters in the model.
Under our vanilla notation, we are interested in learning parameters θ=(π, μ), with π=Γ and μ=D, from the posterior distribution P(π, μ, z |  Y^n), where S^N={S_n: n=1,…,N} with S_n={z_dn:  d=1,…, D} are latent variables, Y^N={Y_n: n=1,…,N} with Y_n={w_dn:  d=1,…, D} are the data, and the priors for (π, μ) are independent Dirichlet distributions _V(η_β) and _K(η_γ) whose densities are denoted by p_π and p_μ. The conditional distribution p(Y^N | μ, S^N) of the observation given the latent variable is 

    (w_dn | μ, z_dn) ∼(β_z_dn), d∈[D]  n∈[N].

We consider the following mean-field approximation <cit.> by decomposing the variational distribution into

    q(μ, π, S^N)
         =q_π(π)   q_μ(μ) ∏_n=1^Nq_S_n(S_n)
       =
      ∏_k=1^K q_β_k(β_k)
      ∏_d=1^D
      (
      q_θ_d(θ_d)
      ∏_n=1^N
      q_z_dn(z_dn)
      )

for approximating the joint posterior distribution of (μ, π,  S^N). Since for LDA, each observation Y_n is composed of D independent observations, it is natural to present the variational oracle inequality with respect to the average squared Hellinger distance D^-1∑_d=1^D h^2[p_d(· | θ) || p_d(· | θ^∗)], where p_d(· | θ) denotes the likelihood function of the dth observation w_d· in Y_·. We make the following assumption.




  
Assumption S: (sparsity and regularity condition) Suppose for each k, β^∗_k is d_k≪ V sparse, and for each d, γ^∗_d is e_d≪ K sparse. Moreover, there exists some constant δ_0>0, such that each nonzero component of β^∗_k or γ^∗_d is at least δ_0. 



Suppose Assumption S holds. If (∑_d=1^De_d+∑_k=1^K d_k)/(DN)→ 0 as N→∞, then it holds with probability tending to one that as N→∞

    ∫ D^-1∑_d=1^D h^2[p_d(· | θ) || p_d(· | θ^∗)]  _θ(θ) dθ
         ≲∑_d=1^De_d/DN log(DKN) + ∑_k=1^K d_k/DN log(KVN).



Theorem <ref> implies the estimation consistency as long as the “effective" dimensionality ∑_d=1^De_d+∑_k=1^K d_k of the model is o(DN) as the “effective sample size" DN→∞.
In addition, the upper bound depends only logarithmically on the vocabulary size V due to the sparsity assumption.


 


 §.§ Gaussian mixture models

Variational inference methods are routinely used in conjugate exponential-family mixtures <cit.> to speed up computation and perform inference on component-specific parameters, with Gaussian mixtures constituting an important special case. Traditional MCMC methods face difficulties in inferring component-specific parameters due to label-switching <cit.>. 
It has been empirically verified <cit.> that variational inference for Gaussian mixtures provides accurate estimates of the true density as well as the labels (up to permutation of the indices). However, theoretical guarantees of such a phenomenon is an open problem till date.  In this section, we close this gap and provide an affirmative answer under reasonable assumptions on the true mixture density.  
In particular, we show that variational techniques using mean field approximation provide optimal estimates for the component specific parameters up to some permutation of the labels. 
 
Suppose the true data generating model is the d-dimensional Gaussian mixture model with K components,

    Y_i |μ, π∼∑_k=1^K π_k   N(μ_k,   I_d),   i=1, …, n,

where μ_k  is the mean vector associated with the kth component and π=(π_1,…,π_K)∈ S^K-1 is a vector of probabilities lying in the K-dimensional simplex S^K-1.  Assume that μ_k∈ C_μ where C_μ is a compact  subset of ℝ^d.  Set μ = (μ_1, …, μ_K) and θ = (μ, π) as before.  Although we assume the covariance matrix of each Gaussian component to be I_d, it is straightforward to extend our results to diagonal covariance matrices. 

Introducing independent latent variables S_i ∈{1, …, K} for i=1, …, n such that π_k = P(S_i = k) for k=1, …, K, (<ref>) can be re-written as 

    Y_i |μ, π, S^n ∼ N(μ_S_i,   I_d),  
       S_i |θ∼(1; π),   i=1, …, n,


where S^n = (S_1, …, S_n) and (1; π) denotes a discrete distribution with support {1, …, K} and probabilities π_1, …, π_K. For simplicity, we assume independent priors p_μ⊗ p_π for (μ, π).

We apply the mean field approximation by using the family of density functions of the form

    q(π, μ, S^n) = q_π(π)  q_μ(μ)  q_S^n(s^n) = q_π(π)  q_μ(μ) ∏_i=1^n q_S_i(s_i)

to approximate the joint posterior distribution of (π, μ,  S^n).  

Since we are interested in studying accuracy in estimating the component specific parameters (π, μ),  we turn our attention to relating the discrepancy in estimating f with that of (π, μ).  We work with the Wasserstein metric <cit.>  between the mixing measures associated with the density.  Note that (<ref>)-(<ref>) can also be written in terms of the mixing measure P = ∑_k=1^K π_k δ_μ_k as  

    Y_i | P ∼ f (·| P) = ∫ N( · ; μ,   I_d)  dP(μ).

Henceforth {μ_k:k = 1, …, K} will be referred to as the atoms of P.  Let 𝒫 denote the class of all such mixing measures 

    P = {P = ∑_k=1^K π_k δ_μ_k : π∈ S^K-1,  μ_k ∈ C_μ ∀ k},

such that the atoms lie in a compact subset C_μ of ℝ^d.  Define the L_r-Wasserstein distance, denoted W_r, between two mixing measures P_1 = ∑_k=1^K π_1,kδ_μ_1,k and P_2 = ∑_k=1^K π_2,kδ_μ_2,k  in 𝒫 as 

    W_r^r(P_1, P_2) = inf_C ∈ C_X_1X_2X_1 - X_2^r,

where X_i ∼ P_i for i = 1, 2, and C_X_1 X_2 is the set of all possible couplings, i.e. joint distributions of X_1 and X_2 with marginals P_1 and P_2 respectively.  One can write (<ref>) in terms of (π, μ) as 

    W_r^r(P_1, P_2) = inf_q ∈ C_P_1P_2∑_k, k' q_kk'μ_1,k -  μ_2,k'^r,

where q varies over C_P_1P_2, the set of joint probability mass functions over {1, …, K}^2 satisfying ∑_k q_kk' = π_2,k' and ∑_k' q_kk' = π_1,k.   

In the following, we will consider r=1 to work with the W_1 metric.  It is known <cit.> that 𝒫 is compact with respect to W_1.  Let h{ f(·| P_1)  ||  f(·| P_2)} denote the Hellinger distance between the densities with corresponding mixing measures P_1 and P_2.  Denote by N(ε, 𝒫, W_1) and N(ε, ℱ,h) the ε-covering numbers of 𝒫 and the corresponding space of densities ℱ with respect to W_1 and h respectively. 


Following <cit.>, we investigate the minimum separation between the densities in (<ref>) in Hellinger distance when the corresponding mixing measures are separated by at least ε in W_1.   Fix   P^∗ = ∑_k=1^K π_k^∗δ_μ_k^∗∈𝒫 and P ∈𝒫. 




From Birgé-Lecam theory <cit.>, there exists a sequence test functions Φ_n based on observations Y_1, …, Y_n such that 

    _P^∗Φ_n ≤ N(ε, ℱ, h)  e^-C_1n ε^2
    
       _P [ 1- Φ_n]  ≤ e^-C_2n h^2[f(· | P) || f(· |  P^∗)],
 
 for any P with h[f(· | P) || f(· |  P^∗)] > ε. <cit.> discusses construction of test functions in the W_1 metric. Using their Theorem 3.1 along with (<ref>)-(<ref>), one obtains a test function Ψ_n such that 

    _P^∗ [Ψ_n] ≤  N(ε, P, W_1)  e^-C_1 n ε^2
       _P [ 1- Ψ_n]  ≤ e^-C_2 n W_1^2(P, P^∗),

for any P ∈ P with W_1(P, P^∗) > ε. We show in the supplement that N(ε, P, W_1) ≲ e^C n ε^2. Hence Assumption T is satisfied with κ = 2. We remark that for the W_2 metric, it is possible to construct a test with κ = 4. 

To state the risk bounds in the W_1 metric, we require the following assumptions on the component specific weights  associated with the true mixing measure P^∗ and the prior densities p_μ and p_π respectively. 


  
Assumption R: (lower bound on component weights)  There exists some constant δ∈ (0, 1), such that the true π_k^∗ > δ for all k = 1, …, K. 




  
Assumption P: (Prior thickness)  Assume that the prior densities for μ and π satisfy 
p_π(π^∗)>0 and p_μ(μ^∗)>0. Additionally, assume that p_μ is supported on the compact set C_μ.



Suppose Assumption R holds, and the prior densities p_μ and p_π satisfy Assumption P.  Then if  d K/n→ 0 as n→∞,  it holds with probability tending to one as n→∞ that

    ∫ h^2[p(· | θ) ||  p(· | θ^∗)]  _θ(θ) dθ≲d K/n log (dn), 
    ∫ W_1^2(P, P^∗)  _θ(θ) dθ≲d K/n log (K n).

 
The convergence rates obtained in (<ref>)-(<ref>) are minimax upto logarithmic terms  <cit.>.
An important consequence of Theorem <ref> is that it allows us to study the accuracy of the estimates of the component specific parameters (π_k, μ_k) for k=1, …, K. The following lemma relates the accuracy in estimating the mixing measures with respect to the W_1 distance and that for the component specific weights and atoms.  The convergence of the weights requires an additional assumption in terms of separability of atoms in P^∗ which prohibits two distinct atoms of the true mixing distribution P^∗ to be vanishingly close.   



  
Assumption S: (separability of atoms)  There exists some ζ>0, such that  inf_k≠ k'μ_k^∗ -   μ_k'^∗≥ζ.


 Fix P^∗∈𝒫 such that P^∗ satisfies Assumption R.  Given ε > 0, if P = ∑_k=1^K π_k δ_μ_k satisfies W_1(P^∗, P) < ε, then  for ε < ζδ, 

    μ_k - μ_k^∗ < ε/δ,  ∀ k=1, …, K,

upto a permutation of the indices.  
If in addition, P^∗ satisfies  both Assumption R and S, then,  

    ∑_k=1^K |π_k - π^∗_k| <  ε/ (ζ - ε/δ).














Theorem <ref> together with Lemma <ref> implies the following corollary about parameter convergence rates.

 
 Let  M_n ↑∞ be any sequence of positive numbers and  ε_n =  √((dK /n)log (K n)). If 
 P^∗ satisfies Assumption R,   
 then _θ( μ_k^∗ -  μ_k < M_n ε_n,  ∀  k) → 1 
 If P^∗ satisfies both assumptions R and S, then  
_θ( μ_k^∗ -  μ_k < M_n ε_n,  ∀  k,  ∑_k=1^K |π_k^∗ - π_k| < M_nε_n ) → 1
 
 
According to <cit.>, n^-1/2 is the minimax rate of estimating μ_k's when the number of components are correctly specified to be K (π_k's are not estimable if some μ_k^∗'s are the same). The second part of this corollary implies n^-1/2 convergence rates (up to log n factors) for both μ and π under the strong identifiability assumption S.



§ CONCLUSION


We have provided general purpose tools to verify the statistical accuracy of point estimates obtained from mean-field variational inference. Our analysis incorporates latent variables commonly augmented to simplify variational inference, and hence our theory is applicable to a broad range of existing algorithms in exactly the way they are practically implemented. The theory does not require prior conjugacy or analytical tractability of the iterative algorithms. The only two ingredients are the existence of appropriate tests which quantify the ability of the likelihood to separate points in the parameter space relative to the distance metric d, and the prior mass assigned to appropriate KL neighborhoods of the true parameter. Since similar quantities also appear in quantifying the convergence rate of the actual posterior <cit.>, one can leverage on this large body of literature along with our theory to offer insights into variational point estimates. Future work will involve relaxing the prior compactness assumption and considering more general variational families.  



§ PROOFS
 


 §.§ Convention

Equations in the main document are cited as (1), (20 etc., retaining their numbers, while new equations defined in this document are numbered (S1), (S2) etc. 



 §.§ Proof of Theorem <ref>

As in the proof sketch in the main document, our first step is to show that under the testing assumption T, 

    ∫_Θξ(θ, θ^∗) p_θ(dθ) ≤ e^C n ε_n^κ,

w.h.p. (w.r.t. P_θ^∗^(n)), where recall logξ(θ, θ^∗) = ℓ_n(θ, θ^∗) + n d^κ(θ, θ^∗). We first establish (<ref>). Define 

    T_1 = ∫_d(θ, θ^∗) ≤ε_nξ(θ, θ^∗) p_θ(dθ), 
        T_2 = ∫_d(θ, θ^∗) > ε_nξ(θ, θ^∗) p_θ(dθ).

Let us first tackle T_1. Since E_θ^∗ [e^ℓ_n(θ, θ^∗] = 1, we have, 

    E_θ^∗ T_1 = ∫_d(θ, θ^∗) ≤ε_n e^n d^κ(θ, θ^∗) p_θ(dθ) ≤ e^n ε_n^κ.

Hence, by Markov's inequality, T_1 ≤ e^C n ε_n^κ with probability at least 1 - e^- C n ε_n^κ. 

Let us now focus on T_2. Write T_2 = T_21 + T_22, where 

    T_21 = ∫_d(θ, θ^∗) > ε_n (1 - ϕ_n)  ξ(θ, θ^∗) p_θ(dθ), 
        T_22 = ∫_d(θ, θ^∗) > ε_nϕ_n  ξ(θ, θ^∗) p_θ(dθ),

where ϕ_n is the test function from Assumption T. Focus on T_21 first. Observe 

    E_θ^∗ T_21    = ∫_d(θ, θ^∗) > ε_n E_θ[1 - ϕ_n]   e^n d^κ(θ, θ^∗) p_θ(d θ) 
       ≤ e^- C n ε_n^κ.

This implies, by Markov's inequality, than T_21≤ e^- C n ε_n^κ with probability at least 1 - e^- C n ε_n^κ. 

Finally, focus on T_22. Since E_θ^∗ [ϕ_n] ≤ e^-  n ε_n^κ, it follows from Markov's inequality that ϕ_n ≤ e^- C n ε_n^κ with probability at least 1 - e^- C n ε_n^κ. Hence, T_22≤ e^-C n ε_n^κ T_2 w.h.p. Adding the w.h.p. bound for T_21, we obtain, w.h.p.,

    T_2≤ e^-C n ε_n^κ T_2 + e^-C n ε_n^κ.

Rearranging, T_2 ≤ e^- C n ε^κ with probability at least 1 - e^- C n ε_n^κ. Combining with the bound for T_1, (<ref>) is established. 

Once (<ref>) is established, the next step is to link the integrand in (<ref>) with the latent variables. To that end, observe that  

    ξ(θ, θ^∗) = ∑_s^nexp{h(θ, s^n)} _S^n(s^n),

where 

    h(θ, s^n) = logp(Y^n | μ, s^n)  π_s^n/p(Y^n | θ^∗)  _S^n(s^n) + n d^κ(θ, θ^∗).

Combining the above with (<ref>), we have, w.h.p.,

    ∫_Θ∑_s^nexp{ h(θ, s^n) } _S^n(s^n)   p_θ(d θ) ≤ e^C n ε_n^κ.

Next, use a well-known variational/dual representation of the KL divergence (see, e.g., Corollary 4.15 of <cit.>) which states that for any probability measure μ and any measurable function h with e^h ∈ L_1(μ), 

    log∫ e^h(η) μ(d η) = sup_ρ[ ∫ h(η)  ρ(d η) -  D(ρ || μ) ],

where the supremum is over all probability measures ρ≪μ. In the present context, setting η = (θ, s^n), μ : = _S^n⊗ p_θ, and ρ = _θ⊗_S^n, it follows from the variational lemma (<ref>) and some rearrangement of terms that w.h.p.

    n ∫_Θ d^κ(θ, θ^∗)  _θ(d θ) ≤ n ε_n^κ + D(_θ || p_θ)  - ∫_Θ∑_s^n h(θ, s^n)  _θ(d θ).

From (<ref>)–(<ref>) (in the main document), it follows that the right hand side of the above display equals n ε_n^κ + Ω(_θ, _S^n). The proof of the theorem then follows, since by definition, Ω(_θ, _S^n) ≤Ω(q_θ, q_S^n) for any (q_θ, q_S^n) in the variational family Γ. 





 §.§ Proof of Lemma <ref>


Since W_1(P^∗, P) < ε,  there exists a coupling q such that  ∑_k, k' q_kk'μ_k^∗ -  μ_k' < ε.   Then ∑_kπ_k^∗inf_k'μ_k^∗ -  μ_k' < ε. Since π_k^∗≥δ, we have  inf_k'μ_k^∗ -  μ_k'≤ε/δ for all k= 1, …, K.  
 This means for any k, there exists a k' such that μ_k^∗ -  μ_k' < ε/δ.  Without loss of generality, let k' = k.   This proves the first part of the assertion.   To prove the second part, observe that for k ≠ k', μ_k^∗ -  μ_k'≥ζ -  μ_k'^∗ -  μ_k'≥κ - ε/δ.  Then 
 
    ε > W_1(P^∗, P) ≥inf_q∑_k≠ k' q_kk'μ_k^∗ -  μ_k'
          ≥ (ζ - ε/δ)
     inf_C ∈ C_XY(X ≠ Y) 
           = (ζ - ε/δ) ∑_k=1^K |π_k^∗ - π_k|,

 implying ∑_k=1^K |π_k^∗ - π_k| ≤ε / (ζ - ε/δ).   




 §.§ Proof of Theorem <ref>

We first ensure the existence of the test functions Φ_n, and Ψ_n as described in 
(<ref>)-(<ref>). 

First, we find find the covering numbers N(ε, 𝒫, W_1) and N(ε, ℱ, h) to upper bound the Type I and II errors of the test functions Φ_n and Ψ_n.  Note that 

    h^2[f(· | P_1) || f(· |  P_2)] ≤∑_k=1^K |π_1,k - π_2,k| + 
    ∑_k=1^k π_1,kμ_1,k - μ_2,k.

Hencd N(ε, ℱ, h) ≤ N(ε^2/2, 𝒮^K-1, || · ||_1) ×{N(ε^2/2, C_μ, || · ||)}^K where || · ||_1 denotes the L_1 norm between
 two probability vectors and  || · || denotes the Euclidean norm.  From Lemma A.4 of <cit.>, we obtain N(ε^2/2, 𝒮^K-1, || · ||_1) ≤ (10/ε^2)^K-1.  Also,  {N(ε^2/2, C_μ, || · ||)}≤  (2C_U/ ε^2)^d for a global constant C_U is the diameter of the set C_μ. Then 
N(ε, ℱ, h) ≤ (C/ε^2)^dK for some constant C > 0.  
To obtain an upper bound for N(ε, 𝒫, W_1), we note that 

    W_1(P_1, P_2) ≤∑_k=1^K max{π_1,k, π_2,k}μ_1,k - μ_2,k
     + C_U∑_k=1^K |π_1,k - π_2,k|.

Hence N(ε, 𝒫, W_1) ≤ N(ε/(2 C_U), 𝒮^K-1, || · ||_1) ×{N(ε/(2K), C_μ, || · ||)}^K ≤( CK/ε)^dK (10/ε)^K-1. 
Hence   log N(ε, ℱ, h) ≲ dK log (1/ε) and log N(ε, 𝒫, W_1) ≲ dK log (K/ε). 
Then, we have from (<ref>)-(<ref>)

    _P^∗Φ_n ≤ e^-C_1n ε^2 + dK log (1/ε)
    _P [ 1- Φ_n]  ≤ e^-C_2n h^2[f(· | P) || f(· |  P^∗)],
 
for any P with h[f(· | P) || f(· |  P^∗)] > ε. 
In this case, we choose ε≡ε_n to be as constant multiple of {(dK/n) log n}^1/2.  
Also, we  have from  (<ref>)–(<ref>)

    _P^∗Ψ_n ≤ e^-C_1 n ε^2 +  dK log (K/ε)
    _P [ 1- Ψ_n]  ≤ e^-C_2 n W_1^2(P, P^∗),

for any P with W_1(P, P^∗) > ε. 
In this case, we choose ε≡ε_n to be as constant multiple of {(dK/n) log (Kn)}^1/2.  



Recall the two KL neighborhoods around (π^∗, μ^∗) with radius (ε_π, ε_μ) as

    B_n(π^∗, ε_π)    = {D(π^∗  || π) ≤ε_π^2,  V(π^∗  || π) ≤ε_π^2}, 
     B_n(μ^∗, ε_μ)    =  {sup_sD[ p(· | μ^∗,s) ||   p(· | μ,s)]≤ε_μ^2, 
           sup_s V[ p(· | μ^∗,s) ||   p(· | μ,s)] ≤ε_μ^2},

where we used the shorthand D(π^∗  || π)=∑_s π^∗_s log(π^∗_s/π_s) to denote the KL divergence between multinomial distributions with parameters π^∗, π∈ S_K. We choose q_θ as the restriction of p_θ into B_n(π^∗, ε_π)× B_n(μ^∗, ε_μ). 

















It is easy to verify that under Assumption R, there exists some constant C_1 depending only on δ_0 such that 
B_n(π^∗, √(K) ε)⊃{π: max_k|π_k-π^∗_k|≤ C_1 ε} (by using the inequality D(p || q) ≥ 2 h^2(p || q)). In addition, for Gaussian mixture model, it is easy to verify that the KL neighborhood B_n(μ^∗, ε)  contains the set {μ: max_kμ_k-μ_k^∗≤ 2 ε}. As a consequence,  with ε_π=√(K) ε and ε_μ=ε yields (using the prior thickness assumption and the fact that the volumes of {π: max_k|π_k-π^∗_k|≤ C_1 ε} and {μ: maxμ_k-μ^∗_k≤ C_2 ε} are at least O(ε^-K) and O((√(d)/ε)^dK) respectively). Then we have from   Theorem <ref>,with probability tending to one as n→∞,

    ∫{h^2[f(· | θ) || f(· | θ^∗)] } _θ(θ) dθ   ≲d K/nlog n+K ε^2  
        + d K/n logd/ε.

Choosing ε = √(d/n) in the above display yields the claimed bound.

Also, we have with high probability 

    ∫{W_1^2[f(· | θ) || f(· | θ^∗)] } _θ(θ) dθ   ≲d K/nlog (Kn) 
        + K ε^2+ d K/n logd/ε.

Choosing ε = √(d/n) in the above display yields the claimed bound noting that the first term in the right hand side of the preceding display is dominant. 





 §.§ Proof of Theorem <ref>

Under the notation in the paper, for each n=1,…,N, the latent variable S_n={z_dn:  d=1,…,D}. We use Theorem <ref> with d=h (Hellinger metric) and view each latent variable S_n per observation in the theorem as a block of D independent latent variable per observation. The existence of the test is automatic <cit.> with the Hellinger metric (parameter space is compact). This leads to that with probability tending to one as N→∞,

    ∫∑_d=1^D h^2[p_d(· | θ) || p_d(· | θ^∗)] dθ≤(∑_d=1^Dε_γ_d^2+∑_k=1^Kε_β_k^2)
            + { - 1/N∑_d=1^Dlog P_γ_d[ B_N(γ^∗_d, ε_γ_d)] }
           + { - 1/N∑_k=1^Klog P_β_k[B_N(β^∗_k, ε_β_k) ] },

where KL neighborhoods B_N(γ_d^∗; ε_γ_d): ={D(γ_d^∗ || γ_d) ≤ε_γ_d^2, V(γ_d^∗ || γ_d)≤ε_γ_d^2}, for d=1,…,D, and  B_N(β^∗_k, ε_β_k) = {max_k D[p(· | β_k, k) ||  p(· | β_k, k)]≤ε_β_k^2, max_S_n V[p(· | β_k, k) ||  p(· | β_k, k)]≤ε_β_k^2}.

Let S^β_k denote the index set corresponding to the non-zero components of β_k for k=1,…,K, and S^γ_d the index set corresponding to the non-zero components of γ_d for d=1,…,D.
Under Assumption S, it is easy to verify that for some sufficiently small constants c_1,c_2>0, it holds for all d=1,…,D that B_N(γ^∗_d, ε_γ_d)⊃{(γ_d)_(S^γ_d)^c_1 ≤ c_1 ε_γ_d,  (γ_d)_S^γ_d- (γ^∗_d)_S^γ_d_∞≤ c_1 ε_γ_d}, and for all k=1,…,K that B_N(β^∗_k, ε_β_k)⊃{(β_k)_(S^β_k)^c_1 ≤ c_2 ε_β_k,  (β_k)_S^β_k- (β^∗_k)_S^β_d_∞≤ c_2 ε_β_d}.
Applying Theorem 2.1 in <cit.>, we obtain the following prior concentration bounds for high-dimensional Dirichlet priors 

    P_γ_d{(γ_d)_(S^γ_d)^c_1 ≤ c_1 ε_γ_d,  (γ_d)_S^γ_d- (γ^∗_d)_S^γ_d_∞≤ c_1 ε_γ_d}
                   ≳exp{-C e_d logK/ε_γ_d}, d=1,…,D;
       P_β_k{(β_k)_(S^β_k)^c_1 ≤ c_2 ε_β_k,  (β_k)_S^β_k- (β^∗_k)_S^β_k_∞≤ c_2 ε_β_d}
                   ≳exp{-C d_k logV/ε_β_k},  k=1,…,K,

for some constant C>0.

Putting pieces together, we obtain

    ∫∑_d=1^D h^2[p_d(· | θ) || p_d(· | θ^∗)] dθ≲(∑_d=1^Dε_γ_d^2+∑_k=1^Kε_β_k^2) 
       + 1/N∑_d=1^D e_d logK/ε_γ_d +
      1/N∑_k=1^K d_k logV/ε_β_k,

which leads to the desired bound by optimally choosing ε_γ_d's and ε_β_k's.



plain



