Concentration of tempered posteriors and of their variational approximations

Introduction

In many applications of Bayesian statistics, the posterior is not tractable. Markov Chain Monte Carlo algorithms (MCMC) were developed to allow the statistician to sample from the posterior distribution even in situations where a closed-form expression is not available. MCMC methods were successfully used in many applications, and are still one of the most valuable tools in the statistician’s toolbox. However, many modern applications of statistics and machine learning involve such massive datasets that sampling schemes such as MCMC have become impractical. In order to allow the use of Bayesian approaches with these datasets, it is actually much faster to compute variational approximations of the posterior by using optimization algorithms. Variational Bayes (VB) has indeed become a corner stone algorithm for fast Bayesian inference.

VB has been applied to many challenging problems: matrix completion for collaborative filtering , NLP on massive datasets , video processing , classification with Gaussian processes , among others. Chapter 10 in is a good introduction to VB and provides an exhaustive survey.

Despite its practical success very little attention has been put towards theoretical guaranties for VB. Asymptotic results in exponential models were provided in . More recently, proposed a very nice asymptotic study of approximations in parametric models. The main problem with these results is that by nature they cannot be applied to high-dimensional or nonparametric models, or to model selection. In the machine learning community, also studied VB approximations. In a distribution-free setting, there is actually no likelihood, but a pseudo-likelihood can be defined through a suitable loss function and thus it is possible to define a pseudo-posterior. Thanks to PAC-Bayesian inequalities from , derived rates of convergence for VB approximation of this pseudo-posterior. However, the tools used in are valid for bounded loss functions, so there is no direct way to adapt this method to study VB approximations when the log-likelihood is unbounded.

In this paper, we propose a general way to derive concentration rates for approximations of fractional posteriors. Concentration rates are the most natural way to assess “frequentist guarantees for Bayesian estimators”: the objective is to prove that the posterior is asymptotically highly concentrated around the true value of the parameter. This approach is now very well understood, we refer the reader to the milestone paper , an account of recent advances can be found in in . Recently, studied the situation where the likelihood $L(\theta)$ is replaced by $L^{\alpha}(\theta)$ for $0<\alpha<1$, leading to what is usually called a fractional or tempered posterior. They proved that concentration of the fractional posterior requires actually fewer hypothesis than concentration of the (true) posterior. Extending the technique of , we analyze the concentration of VB approximations of (fractional) posteriors. Especially, we derive a condition for the VB approximation to concentrate at the same rate as the fractional posterior.

2 Definitions and notations

Let $\alpha\in(0,1)$. Let $P$ and $R$ be two probability measures. Let $\mu$ be any measure such that $P\ll\mu$ and $R\ll\mu$, for example $\mu=P+R$. The $\alpha$-Rényi divergence and the Kullback-Leibler (KL) divergence between two probability distributions $P$ and $R$ are respectively defined by

We remind the reader of a few properties proven in . First, it is obvious that $D_{\alpha}(P,R)$ does actually not depend on the choice of the reference measure $\mu$. This is sometimes made explicit by the (informal) statement $D_{\alpha}(P,R)=(1/(\alpha-1))\log\int({\rm d}P)^{\alpha}({\rm d}R)^{1-\alpha}$. The measures $P$ and $R$ are mutually singular if and only if $D_{\alpha}(P,R)=(\frac{1}{\alpha-1})\log(0)=+\infty$.

We have $\lim_{\alpha\rightarrow 1}D_{\alpha}(P,R)=\mathcal{K}(P,R)$ which gives ground to the notation $D_{1}(P,R)=\mathcal{K}(P,R)$. For $\alpha\in(0,1]$, $(\alpha/2)d^{2}_{TV}(P,R)\leq D_{\alpha}(P,R)$, $d_{TV}$ being the total variation distance – for $\alpha=1$ this is Pinsker’s inequality. The map $\alpha\mapsto D_{\alpha}(P,R)$ is nondecreasing. Also, the authors of note that the $\alpha$-Rényi divergences are all equivalent for $0<\alpha<1$, through the formula $\frac{\alpha}{\beta}\frac{1-\beta}{1-\alpha}D_{\beta}\leq D_{\alpha}\leq D_{\beta}$ for $\alpha\leq\beta$. Additivity holds: $D_{\alpha}(P_{1}\otimes P_{2},R_{1}\otimes R_{2})=D_{\alpha}(P_{1},R_{1})+D_{\alpha}(P_{2},R_{2})$, thus $D_{\alpha}(P^{\otimes n},R^{\otimes n})=nD_{\alpha}(P,R)$; $D_{1/2}(P,R)\geq 2[1-\exp(-(1/2)D_{1/2}(P,R))]=H^{2}(P,R)$ the squared Hellinger distance.

The fractional posterior, that will be our ideal estimator, is given by

Let $\mathcal{F}\subset\mathcal{M}_{1}^{+}(\Theta)$,

In Sections 3, 4 and 5 we apply our general results in various settings. In Section 3 we study the parametric family of Gaussian approximations

Main results

For any $\alpha\in(0,1)$, for any $\varepsilon\in(0,1)$,

It is tempting to minimize the right-hand side (r.h.s) of the inequality in order to ensure a good estimation. The minimizer of the r.h.s can actually be explicitly given. In order to do this, let us recall Donsker and Varadhan’s variational inequality (Lemma 1.1.3 in ).

Using Lemma 2.2 with $h(\theta)=-\alpha r_{n}(\theta,\theta_{0})$ and the definition of $\pi_{n,\alpha}$ we obtain

so the minimizer of the r.h.s of Theorem 2.1 is actually $\pi_{n,\alpha}({\rm d}\theta|X_{1}^{n})$.

Theorem 2.1 can be used to study other approximations of the posterior. For example, as suggested by one of the Referees, we can use it to study distributions centered around the maximum a posteriori (MAP) or the maximum likelihood estimate (MLE). For example, Laplace approximations are Gaussian distributions centered at the MLE. However, there are models $(P_{\theta},\theta\in\Theta)$ where the MLE and the MAP are not defined, while the posterior and some variational approximations are consistent. Such an example is provided in the Supplementary Material.

2 Concentration of VB approximations

We specialize the above results to the variational approximation. Elementary calculations show that

As a consequence, we obtain the following corollary of Theorem 2.1.

For any $\alpha\in(0,1)$ and $\varepsilon\in(0,1)$, with probability at least $1-\varepsilon$,

Fix $\mathcal{F}\subset\mathcal{M}_{1}^{+}(\Theta)$. Assume that a sequence $\varepsilon_{n}>0$ is such that there is a distribution $\rho_{n}\in\mathcal{F}$ such that

Then, for any $\alpha\in(0,1)$, for any $(\varepsilon,\eta)\in(0,1)^{2}$,

This theorem is a consequence of Corollary 2.3, its proof is provided in Section 7. Let us now discuss the main consequences of this theorem.

Note that the assumption involving a distribution $\rho_{n}$ is not standard. This requires some explanations. Consider first the case $\mathcal{F}=\mathcal{M}_{1}^{+}(\Theta)$. Define $B(r)$, for $r>0$, as

Then the choice $\rho_{n}=\pi_{|B(\varepsilon_{n})}$, i.e. $\pi$ restricted to $B(\varepsilon_{n})$, ensures immediately (2.1), and (2.2) can be rewritten

This assumption is standard to study concentration of the posterior, see Theorem 2.1 page 503 in or Subsection 3.2 in . Our message is that in the studies of concentration of the posterior, the choice $\rho_{n}=\pi_{|B(\varepsilon_{n})}$ was hidden. Other choices might lead to easier calculations in some situations. More importantly, in the relevant case $\mathcal{F}\subsetneq\mathcal{M}_{1}^{+}(\Theta)$, $\pi_{|B(\varepsilon_{n})}\notin\mathcal{F}$ in general. Thus $-\log\pi(B(\varepsilon_{n}))\leq n\varepsilon_{n}$ is no longer sufficient, and (2.1) and (2.2) are natural extensions of this assumption to study VB. They provide an explicit condition on the family $\mathcal{F}$ in order to ensure concentration of the approximation.

Choosing $\eta=\frac{1}{n\varepsilon_{n}}$ and $\varepsilon=\exp(-n\varepsilon_{n})$ we obtain a more readable concentration result. It shows that, as soon as $(1/n)\ll\varepsilon_{n}\ll 1$, the sequence $\varepsilon_{n}$ gives a concentration rate for VB.

Under the same assumptions as in Theorem 2.4,

As a special case, when $\alpha=1/2$, the theorem leads to a concentration result in terms of the more classical Hellinger distance

Also, with a general $\alpha\in(0,1)$, from the properties recalled in Remark 1.1, we have, for $0<\beta\leq\alpha$,

3 A simpler result in expectation

It is possible to simplify the assumptions at the price of stating a result in expectation instead of concentration.

Fix $\mathcal{F}\subset\mathcal{M}_{1}^{+}(\Theta)$. Then

Assume that $\varepsilon_{n}>0$ is such that there is distribution $\rho_{n}\in\mathcal{F}$ such that

4 Extension of the result in expectation to the misspecified case

In this section we do not assume any longer that the true distribution is in $\{P_{\theta},\theta\in\Theta\}$. In order not to change all the notations we define an extended parameter set $\Theta\cup\{\theta_{0}\}$ where $\theta_{0}\notin\Theta$ and define $P_{\theta_{0}}$ as the true distribution. Theorem 2.6 can be applied to this setting, and we obtain:

Now, rewriting, for $\theta^{*}\in\Theta$,

Assume that, for $\theta^{*}=\arg\min_{\theta\in\Theta}\mathcal{K}(P_{\theta_{0}},P_{\theta})$, there is $\varepsilon_{n}>0$ and $\rho_{n}\in\mathcal{F}$ with

In the well-specified case, $\theta^{*}=\theta_{0}$ and we recover Theorem 2.6. Otherwise, this result takes the form of an oracle inequality. It is not a sharp oracle inequality as that the risk measure used in the l.h.s and the r.h.s are not the same, but remains informative when $\mathcal{K}(P_{\theta_{0}},P_{\theta^{*}})$ is small. For example, in Section 5 below, we provide a nonparametric example where $\mathcal{K}(P_{\theta_{0}},P_{\theta^{*}})$ and Theorem 2.7 leads to the minimax rate of convergence.

Gaussian variational Bayes

thus the algorithm will consist in projecting onto the set of Gaussian distributions. Depending on the hypotheses made on the covariance matrix we can build different approximations. For instance define:

We have by definition $\mathcal{F}^{\Phi}_{id}\subseteq\mathcal{F}^{\Phi}_{diag}\subseteq\mathcal{F}^{\Phi}$.

The remarkable fact of Gaussian VB is that it allows to recast integration as a finite dimension optimization problem. The choice of a specific Gaussian is a trade off between accuracy and computational complexity. We will show in the following that, under some assumption on the likelihood, the integrated $\alpha$-Rényi divergence is convergent for most of the approximations.

To simplify the exposition of the results we will restrict our study to the case of Gaussian priors: $\pi=\mathcal{N}(0,\vartheta^{2}I_{p})$. One can readily see that in Theorem 2.4 the prior appears only in the condition $\frac{1}{n}\mathcal{K}(\rho,\pi)\leq\varepsilon_{n}$, many other distribution could be used, providing different rates.

In the rest of the section we assume that the density is log Lipschitz.

There is a measurable real function $M(\cdot)$ such that

An example is logistic regression, see Subsection 3.2 below.

Let the approximation family be $\mathcal{F}$ with $\mathcal{F}^{\Phi}_{id}\subset\mathcal{F}$ as defined above and that the model satisfies Assumption 3.1. We put

Then for any $\alpha\in(0,1)$, for any $\eta,\epsilon$

In many cases the model is not conjugate, i.e. the VB objective does not have a closed-form solution. We can however use a full Gaussian approximation and a stochastic gradient descent on the objective function defined by the KL divergence. This approach has been studied in .

We may write our variational bound as the following minimisation problem

In the authors suggest using a parametrization of the problem where we replace the optimization over $\Sigma$ by a minimization over the matrix $C$ where $CC^{t}=\Sigma$. To simplify the notations in this section define

to be the objective of the minimization problem (3.1), where $\xi\sim\mathcal{N}(0,I_{d})$ and

In order to be able to state non-asymptotic results on the stochastic gradient algorithms, we restrict the parameter space to an Euclidean ball, that is (3.1) is transformed into

The objective can now be replaced by a Monte Carlo estimate and we can use stochastic gradient descent as described in Algorithm 1.

Assume that $f$, as defined in (3.2), is convex in its first component $x$ and that it has $L$-Lipschitz gradients.

On most examples the gradient is a sum of at least $n$ components. If each term is Lipschitz with constant $L_{i}$, an estimate of the constant will be $L\leq n\max_{i}L_{i}$. The additional term of the bound is therefore of the order $(2B\max L_{i}/(nT))^{1/2}/(1-\alpha)$, hence a good choice is $T=O(\sqrt{n})$ to mitigate the impact of the numerical approximation on the rate.

2 Example: logistic regression

We consider the case of a binary regression model. Although estimation of parameters is relatively simple for small datasets , it remains challenging when the size of the dictionary is large. Furthermore usual deterministic methods do not come with theoretical guarantees as would a gradient descent algorithm for maximum likelihood. Note that the logistic regression is not conjugate in the sense that we cannot find an iterative scheme based on a mean field approximation, as will be done for the matrix completion example in Section 4.

We will prove results in the case of random design where we suppose that the distribution of $Z_{1}^{n}$ does not depend on the parameter.

then for any $\alpha\in(0,1)$, for any $\eta,\epsilon$

Note that the only assumption on the distribution of $X_{1}$ is that $K_{2}<\infty$. Still, it is interesting to compute $K_{1}$ and $K_{2}$ on some examples. For example, when $X_{1}$ is uniform on the unit sphere, $K_{1}\leq 2$ and $K_{2}\leq 4$. When $X_{1}\sim\mathcal{N}(0,s^{2}I_{d})$ then $K_{2}=4s^{2}d$ and $K_{1}\leq 2\sqrt{s^{2}d}$. In both cases, the terms in $K_{1}$ and $K_{2}$ do not deterioriate the parametric rate of convergence $d/n$. Furthermore the Lipschitz constant can be bounded explicitly under additional assumptions on the design matrix (e.g. bounded singular value) and leads to $L=\mathcal{O}\left(nd+\frac{d}{\psi}\right)$. Hence taking $\psi=1/(n\sqrt{d})$ one would get a bound in $\mathcal{O}\left(\sqrt{\frac{d^{3/2}}{nT}}+(d/n)\log{nd}\right)$. We can take $T$ of the order $\frac{n}{d^{1/4}}$ in order not to deteriorate the rate.

Application to matrix completion

where the $(i_{k},j_{k})$ are i.i.d $\mathcal{U}(\{1,\dots,m\}\times\{1,\dots,p\})$. For the sake of simplicity we will assume that the $\xi_{k}$ are i.i.d $\mathcal{N}(0,\sigma^{2})$, and that $\sigma^{2}$ is known, so we only have to estimate $M$. Note that for $\alpha\leq 1$, $\mathcal{D}_{\alpha}(\mathcal{N}(\mu_{1},\sigma^{2}),\mathcal{N}(\mu_{2},\sigma^{2}))=\alpha(\mu_{1}-\mu_{2})^{2}/(2\sigma^{2})$, see (10) page 3800 in . Thus, for $0<\alpha<1$,

which depends only on $\alpha$, $\sigma_{2}$ and the matrices $M$ and $N$ so we will use the notation $d_{\alpha,\sigma}(M,N)=\mathcal{D}_{\alpha}(P_{M},P_{N})$. In the case $\alpha=1$,

where $\|\cdot\|_{F}$ denotes the Frobenius norm. In the noiseless case $\sigma^{2}=0$, proved that it is possible to recover exactly $M$ under the assumption that its rank is small enough. Various extensions to noisy settings, approximately low-rank matrices, or other loss functions can be found in . The main message of these papers is that the minimax rate of convergence is $(m+p){\rm rank}(M)/n$, possibly up to log terms. Bayesian estimators were proposed in using factorized Gaussian priors. Convergence of the posterior mean was proven in for a bounded prior, excluding the Gaussian prior used in practice. Similarly, proves concentration of a truncated version of the posterior. For very large datasets the MCMC algorithm proposed in is too slow, a VB approximation was proposed in with very good results on the Netflix dataset. This approximation was re-used and extended by many authors including . But the consistency of the Bayesian estimator with Gaussian priors and of its variational approximations are opened questions.

First, we will recall the Gaussian prior and the VB approximation . We will then prove the concentration of the VB approximation, and as a consequence the concentration of the tempered posterior.

2 Definition of the prior and of the VB approximation

Fix $K\in\{1,\dots,m\wedge p\}$. The main idea of factorized priors is that, when ${\rm rank}(M)\leq K$ then we have

for some matrices $U$ of dimension $p\times K$ and $V$ of dimension $m\times K$. Thus, we can define a prior on $M$ by specifying priors on $U$ and $V$. A usual choice is that the entries $U_{i,k}$ and $V_{j,k}$ are independent $\mathcal{N}(0,\gamma_{k})$ and finally $\gamma_{k}$ is inverse gamma, that is $1/\gamma_{k}\sim\Gamma(a,b)$. These choices ensure conjugacy: put $\gamma=(\gamma_{1},\dots,\gamma_{K})$, it is then possible to compute the conditional posteriors of $U|V,\gamma$, of $V|U,\gamma$ and $\gamma|U,V$. This allows to use the Gibbs sampler . For large datasets, proposed mean-field VB with $\mathcal{F}$ given by

The minimization of the VB program is shown in many cited papers, see and all the references therein. Shortly: $\rho_{U_{i}}$ is $\mathcal{N}(\mathbf{m}_{i,\cdot}^{t},\mathcal{V}_{i})$, $\rho_{V_{j}}$ is $\mathcal{N}(\mathbf{n}_{j,\cdot}^{t},\mathcal{W}_{j})$ and $\rho_{\gamma_{k}}$ is $\Gamma(a+(m_{1}+m_{2})/2,\beta_{k})$ for some $m\times K$ matrix $\mathbf{m}$ whose rows are denoted by $\mathbf{m}_{i,\cdot}$, some $p\times K$ matrix $\mathbf{n}$ whose rows are denoted by $\mathbf{n}_{j,\cdot}$ and some vector $\beta=(\beta_{1},\dots,\beta_{K})$. The parameters are updated iteratively through the formulae

(where $(\mathcal{V}_{i})_{k,k}$ denotes the $(k,k)$-th entry of the matrix $\mathcal{V}_{i}$ and $(\mathcal{W}_{j})_{k,k}$ denotes the $(k,k)$-th entry of the matrix $\mathcal{W}_{j}$).

3 Concentration of the posteriors

Fix $a$ as any constant. There is a small enough $b>0$ such that

where the constant $\mathcal{C}(a)=\log(8\sqrt{\pi}\Gamma(a)2^{10a+1})+3$. In particular, the result holds for the choice $b=B^{2}/\{512(nmp)^{4}[(m\vee p)K]^{2}\}$.

In practice, it is important that $b$ is small to ensure a good approximation of low-rank matrices . We don’t claim that $b=B^{2}/\{512(nmp)^{4}[(m\vee p)K]^{2}\}$ is the optimal value, recommends cross-validation to tune $b$.

Note as a special case that when $M=\bar{U}\bar{V}^{t}$ for $(\bar{U},\bar{V})\in\mathcal{M}(r,B)$ then we have exactly

This result is the first consistency result for the VB approximation with Gaussian priors, that is used in practice. Still, it is stated for a “weak” distance criterion $d_{\alpha,\sigma}(M,M_{0})$. Under additional assumptions, it is actually possible to relate this criterion to the standard Frobenius norm. Assume that there is a known $C$ such that $\max_{i,j}|(M_{0})_{i,j}|\leq C$. This assumption is satisfied in many applications like collaborative filtering: in the Netflix data the entries are between $1$ and $5$. Then it is natural to project any estimator to the set of matrices with bounded entries. Precisely, define for any $M$ the matrix ${\rm clip}_{C}(M)$ its $(i,j)$-th entry: $\min(\max(M_{i,j},-C),C)$. A simple study of $d_{\alpha,\sigma}$, detailed in the proofs section, leads to the following result.

Under the assumptions of Theorem 4.1, and when in addition $\max_{i,j}|(M_{0})_{i,j}|\leq C$, then

Note that once the Gaussian approximation of the posterior is known, it is easy to sample from it and to clip the samples to approximate the posterior mean of ${\rm clip}(M)$. So under the boundedness assumption we have a bound based on the Frobenius norm for an effective procedure based on VB. It is known that for the squared Frobenius norm, the rate $r(m+p)/n$ is minimax optimal – maybe up to log terms .

Still assuming that $M=\bar{U}\bar{V}^{t}$ for $(\bar{U},\bar{V})\in\mathcal{M}(r,B)$ it is also possible to state a proper concentration result as an application of Corollary 2.5. We omit the proof as it is exactly similar to the one of Theorem 4.1.

Assume $M=\bar{U}\bar{V}^{t}$ for $(\bar{U},\bar{V})\in\mathcal{M}(r,B)$ and take $b$ as in Theorem 4.1. Then

where for some explicit constant $\mathcal{D}(a,\sigma^{2},B)$,

Nonparametric regression estimation

In this section, we provide a nonparametric example. Thus, the parameter will actually be a function $f$. We assume that $X_{1}=(W_{1},Y_{1}),\dots,X_{n}=(W_{n},Y_{n})$ are i.i.d from a distribution $P_{f_{0}}$, and the model $(P_{f})$ is given by: $W_{i}\sim\mathcal{U}()$ and

where $\xi_{i}\sim\mathcal{N}(0,1)$. We will provide a prior and a mean-field approximation of the posterior. We will show that we estimate the functions $f$ belonging to a Sobolev ellipsoid $\mathcal{W}\left(r,C^{2}\right)$ at the minimax rate of convergence, up to $\log$ terms (the definitions of the ellipsoids will be reminded below). The reader might think that this example is not the most striking application of VB. On the other hand, it is an illustration of the generality of our method. We will estimate $f$ using projections on the Fourier basis and the choice of the number of coefficients will be done by model selection. It appears that in this case, model selection can be seen as a variational approximation where the constraint on the posterior is to give all its mass to only one model. This leads to adaptation of the estimator, in the sense that it is not required to know $r$ nor $C$ to compute the estimator.

First, we recall the definition of the trigonometric basis $(\varphi_{k})_{k=1}^{\infty}$:

We now define a prior distribution $\pi$ by describing how to draw from $\pi$: we first draw $K$ from a geometric distribution, $\pi(K=k)=2^{-k}$. We then draw $\beta_{1},\dots,\beta_{K}$ i.i.d from a $\mathcal{N}(0,1)$ distribution. We finally put

Note that when $f_{0}(\cdot)=\sum_{k=1}^{\infty}\beta_{k}^{0}\varphi_{k}(\cdot)$ and all the $\beta_{k}^{0}$’s are non-zero, such a function is never “produced” by the prior. Still, we will see that the prior gives enough mass to functions in the neighborhood of $f_{0}$, ensuring consistency.

2 Construction of the variational approximation

Note that the support of $\pi(\cdot|K)$ has dimension $K$, but the support of $\pi$ is infinite-dimensional. Thus, we can expect the support of the tempered posterior $\pi_{n,\alpha}$ to be also infinite-dimensional, and $\pi_{n,\alpha}$ to be intractable. We define a variational approximation that will fix these problems.

First, for $K\geq 1$ define $\mathcal{F}_{K}$ as the set of probability measures $\rho_{\mathbf{m},s^{2}}$ where $\mathbf{m}=(m_{1},\dots,m_{K})$ on functions $f(\cdot)=\sum_{k=1}^{K}\beta_{k}\varphi_{k}(\cdot)$ such that under $\rho_{\mathbf{m},s^{2}}$, the $\beta_{k}$’s are independent and $\beta_{k}\sim\mathcal{N}(m_{k},s^{2})$. We put $\mathcal{F}=\bigcup_{k=1}^{\infty}\mathcal{F}_{k}$. Note that the choice of a constant variance $s^{2}$ was motivated by the fact that the estimator of $\beta_{k}$ studied for example in , $\hat{\beta}_{k}=(1/n)\sum_{i=1}^{n}Y_{i}\varphi_{k}(X_{i})$, satisfies $\hat{\beta}_{k}\sim\mathcal{N}(\beta_{k},\sigma^{2}/n)$. Then

(that is, the approximated posterior mean is simply a ridge regression estimator).

3 Nonparametric rates of convergence

We remind the definition of the Sobolev ellipsoid given (see e.g. Chapter 1 in ) for $C>0$ and $r\geq 2$:

Fix $\alpha\in(0,1)$. Assume that there is an $r\in[2,\infty[$ and a $C>0$ such that $f_{0}\in\mathcal{W}(r,C^{2})$. Then

The proof is in Section 7. Note that on the contrary to previous sections, we only provide an asymptotic statement here. However, from the proof of Theorem 5.1, it is clear that it is possible to provide a non-asymptotic statement as well (with cumbersome constants).

Here again, note that the distance criterion used in the left-hand side is not standard. We actually have:

However, when $f_{0}\in\mathcal{W}_{r,C^{2}}$, $f$ is bounded by a constant that depends on $r$ and $C^{2}$. If we moreover assume that $f_{0}$ is bounded by a known constant $c_{0}$, we can as in Section 4 define a clip operator: ${\rm clip}_{c_{0}}(f)(x)=\min(\max(-c_{0},f(x)),c_{0})$ and obtain:

The rate $1/n^{2r/(2r+1)}$ is known to be minimax optimal on $\mathcal{W}(r,C^{2})$ for the squared $\|\cdot\|_{2}$-norm . The additional $\log$ term is sometimes referred to as “the price to pay for adaptation”. In the case of the $\|\cdot\|_{2}$-norm this is misleading as it is actually possible to build an adaptive estimator that reaches the minimax rate without the additional $\log$, but up to our knowledge this is not possible with a fully Bayesian estimator.

Conclusion

Based on PAC-Bayesian inequalities, we introduced a generic method to study the concentration of variational Bayesian approximations. This is a very general approach that can be applied to many models. We studied applications to logistic regression, matrix completion and density estimation. Still, some questions remain open. From a theoretical perspective, the oracle inequality in Theorem 2.7 compares a Rényi divergence to a Kullback-Leibler divergence. It would be very interesting to obtain a result with the Kullback divergence in the left-hand side. This is probably more difficult, if possible at all. We believe that tools from could be of some help, but some work is needed to make explicit the assumptions of this paper in our context.

Also, since the first version of this work was submitted, extensions were proven by other authors: extended our results to models with hidden variables, such as mixture models, and proved results in the case $\alpha=1$ and study many nonparametric examples. Note that while $\alpha=1$ remains the most popular choice in practice, these results require much stronger assumptions and cannot in general be extended to the misspecified case .

An important open issue is the choice of the parameter $\alpha$. It is clear that our results are not helpful to solve this issue. Some previous work proposes to use cross-validation , but this is computationaly expensive. Moreover, no theoretical guarantees are known in this case. In the misspecified case, proposed an online adaptive tuning of this parameter. However, it is not clear if this method could work in our context. This should be the object of a future work.

Finally, it would be nice to get rid of the extra log in the rates. Catoni’s localization technique is a nice tool to remove extra log factors in PAC-Bayesian bounds, but its adaptation to our setting is not direct. It could be the object of future works.

Acknowledgements

We would like to thank Badr-Eddine Chérief-Abdellatif, as well as the Associate Editor and the anonymous Referees, for their helpful comments and suggestions on the paper.

Proofs

We adapt the proof given in . Fix $\alpha\in(0,1)$, and $\theta\in\Theta$. It’s immediate to check that

The key argument here, introduced by , is to use Lemma 2.2. Note that almost surely with respect to the sample, we know that

Multiply both sides by $\varepsilon$ to get

2 Proof of Theorem 2.4

Now apply take the union bound of this inequality and of the inequality in Corollary 2.3. We obtain, for any $\alpha\in(0,1)$, for any $\varepsilon\in(0,1)$, with probability at least $1-\varepsilon-\eta$,

where in the last step we use the assumptions on $\rho_{n}$.

3 Proof of Theorem 2.6

The beginning is as for Theorem 2.1. Fix $\alpha\in(0,1)$, then

This is where things change: we now use Jensen’s inequality to obtain

4 Proof of Theorem 3.1

We start by defining a sequence $\rho_{n}(d\theta):=\Phi(d\theta;\theta_{0},\sigma_{n}^{2}I)\in\mathcal{F}_{\Phi}^{id}$ indexed by a positive scalar $\sigma^{2}_{n}$ to be later defined. As before by proving the result on the smallest family of distribution, it will remain true on larger ones using the fact that $\min_{\mathcal{F}^{id}}\leq\min_{\mathcal{F}^{diag}}\leq\min_{\mathcal{F}^{full}}$. Under Assumption 3.1 we can check the hypotheses on the KL between the likelihood terms as required in Theorem 2.4. We have

When integrating with respect to $\rho_{n}$ we have

To apply Theorem 2.4 it remains to compute the KL between the approximation of the pseudo-posterior and the prior,

To obtain an estimate of the rate $\varepsilon_{n}$ of Theorem 2.4 we put together those bounds. Choosing $\sigma^{2}_{n}=\frac{1}{n\sqrt{d}}$ we can apply it with

5 Proof of Theorem 3.2

Following the rest of the proof of 2.6 we get

To bound the first term of the right hand-side we use Assumption 3.1 and the proof of Theorem 3.1. In particular notice that $\Phi(d\theta;\theta_{0},\frac{1}{n\sqrt{d}}I_{d})\in\mathcal{F}^{\Phi}_{B}$, we get straight away

We now study the term inside the brackets on the right hand-side.

Divide by $T$, take expectation with respect to $(\xi_{t})_{t}$

Notice that $x_{t}$ belongs to the $\sigma$-algebra generated by $(x_{1},\cdots,x_{t-1})$. By a multiple use of the tower property we get,

Putting everything together concludes the proof.

6 Proof of Corollary 3.3

Direct calculation shows that the log-likelihood is $2\|X\|$-Lipschitz hence satisfying Assumption 3.1. We conclude using Theorem 3.1 and the assumption on the design matrix.

7 Proof of Corollary 3.4

Start by noticing that we can take $f$ as

where $\rho(.)=\Phi(.;m,CC^{t})$ the likelihood part is convex with Lipschitz gradient as a composition of a convex and gradient Lipschitz function with a affine map. The Lipschitz constant for this term is bounded by $\sum_{i=1}^{n}\|x_{i}x_{i}^{t}\|$. The KL part can be written as $\mathcal{K}(\rho,\pi)=\frac{\|m\|^{2}}{2\vartheta}+\left(\frac{1}{2\vartheta}\text{trace}(CC^{t})-\log|C|\right)$ which is convex for positive semi-definite $C$. We need to check that the gradients of the objectives are also Lipschitz, the only problematic term is $\log\text{det}(C)$. Denote $(\lambda_{i})$ the eigen values of $\Sigma=CC^{t}$

To apply Theorem 3.2 we also need to check that the new constraint contains the Gaussian distribution used in the proof. This is the case as long as $\psi\leq\sigma^{2}=\frac{1}{n\sqrt{d}}$.

The supplementary material contains the toy example mentioned in Remark 2.1 above.

The remaining proofs, that is, the proofs of Theorems 4.1 and 5.1 and of Corollary 4.2, are also provided in the supplementary material.

References

Supplementary material

In this subsection, we provide the toy example announced in the paper, where the MLE (and thus the MAP) are not defined. Then, we show that there is a variational approximation that leads to a consistent estimator. This also implies that the tempered posterior is also consistent in this case.

The prior $\pi$ is given by: $m\sim\mathcal{N}(0,1)$ and $\sigma^{2}\sim\mathcal{U}(0,1)$ (uniform distribution).

8.2 Non-existence of the MLE

It is easy to check that when $X_{1},\dots,X_{n}$ are i.i.d from $P_{(m_{0},\sigma_{0}^{2})}$ then the likelihood function

Thus, the MLE is not defined. For the same reason, the MAP does not exist either.

8.3 A variational approximation family

Note that the family $\mathcal{F}$ is inspired by Catoni’s point of view to use a “perturbed MLE” in PAC-Bayesian bounds. An application of Theorem 2.6 leads to the following result.

As a corollary, we also have that the tempered posterior $\pi_{n,\alpha}(\cdot|X_{1}^{n})$ satisfies the same inequality.

8.4 Proof of Proposition 7.1

Assume that $m\in[m_{0}-\sqrt{\delta\sigma_{0}^{2}},m_{0}+\sqrt{\delta\sigma_{0}^{2}}]$ and $\sigma^{2}\in[\sigma^{2}_{0}-\delta\sigma^{2}_{0},\sigma^{2}_{0}]$ for some $0<\delta<1$. Then:

This implies that for any $\delta\in(0,1)$,

The value $\delta=1/(2n)$ gives, using $(1-\delta)>1/2$ to simplify things,

9 Proof of Theorem 4.1

Fix $B>0$, $r\geq 1$ and any pair $(\bar{U},\bar{V})\in\mathcal{M}_{r,B}$ and define for $\delta\in(0,B)$ that will be chosen later,

Note that it can be factorized so it belongs to the family $\mathcal{F}$.

We adapt the calculations from but simplify a lot. First, note that

and that for any $(U,V)$ in the support of $\rho_{n}$ we have

with the choice $\delta=B/[8(nmp)^{2}]$ which satisfies $0<\delta<B$. Then, we derive

for any event $E$. We actually take $E=\{\gamma_{1},\dots,\gamma_{r}\in[B^{2},2B^{2}],\gamma_{r+1},\dots,\gamma_{K}\in[s,2s]\}$ and $s\in(0,B^{2})$ is to be chosen later. Then note that

with $s=\frac{1}{2}\left(\frac{\delta}{2(m\vee p)K}\right)^{2}$ which satisfies $0<s<B^{2}$. Then, for $k\leq r$,

where we replaced $\delta$ and $s$ by their respective value. In order to keep the expressions as simple as possible we can use $K\leq m\vee p\leq m+p\leq r(m+p)$ and $2\leq m+p\leq r(m+p)$ to get

We are now in position to apply Theorem 2.6. Then

10 Proof of Corollary 4.2

We start from (4.1). Under the boundedness assumption on $M_{0}$ it is obvious that $\forall M$, $d_{\alpha,\sigma}(M,M_{0})\geq d_{\alpha,\sigma}({\rm clip}_{C}(M),M_{0})$, so

Fix $M$ and for short, put $N={\rm clip}_{C}(M)$. We have:

By assumption, $(N_{i,j}-(M_{0})_{i,j})^{2}/(2\sigma^{2})\leq(2C)^{2}/(2\sigma^{2})=2C^{2}/\sigma^{2}$. Straightforward derivations show that for any $x\in[0,2C^{2}/\sigma^{2}]$ we have

Pluging this into (7.4) gives the result claimed.

11 Proof of Theorem 5.1

Let $(\beta_{k}^{0})$ denote the coefficients of $f^{0}$: $f_{0}=\sum_{k=1}^{\infty}\beta_{k}^{0}\varphi_{k}$. Theorem 2.6 gives:

The choice $(m_{1},\dots,m_{K})=(\beta_{1},\dots,\beta_{K})$ gives:

From Chapter 1 in we know that $f_{0}\in\mathcal{W}(r,C^{2})$ implies $\sum_{k=K+1}^{\infty}(\beta_{k}^{0})^{2}\leq\Lambda(r,C)K^{-2r}$ for some $\Lambda(k,C)$. Moreover, $\sum_{k=1}^{K}(\beta_{k}^{0})^{2}\leq\sum_{k=1}^{\infty}k^{2r}(\beta_{k}^{0})^{2}\leq C^{2}$. So finally:

The choice $K=\lceil(n/\log(n))^{1/(2r+1)})\rceil$ leads to the result.