Adaptive variational Bayes:
Optimality, computation and applications
Adaptive variational Bayes

Ilsang<EMAIL_ADDRESS>Lizhen<EMAIL_ADDRESS>Department of Applied and Computational Mathematics and Statistics, The University of Notre Dame

I. Ohn and L. Lin

In this paper, we explore adaptive inference based on variational Bayes. Although a number of studies have been conducted to analyze contraction properties of variational posteriors, there is still a lack of a general and computationally tractable variational Bayes method that can achieve adaptive optimal contraction of the variational posterior. We propose a novel variational Bayes framework, called adaptive variational Bayes, which can operate on a collection of models with varying dimensions and structures. The proposed framework combines variational posteriors over individual models with certain weights to obtain a variational posterior over the entire model. It turns out that this combined variational posterior minimizes the Kullback-Leibler divergence to the original posterior distribution. We show that the proposed variational posterior achieves optimal contraction rates adaptively under very general conditions and attains model selection consistency when the true model structure exists. We apply the general results obtained for the adaptive variational Bayes to several examples including deep learning models and derive some new and adaptive inference results. Moreover, we consider the use of quasi-likelihood in our framework. We formulate conditions on the quasi-likelihood to ensure the adaptive optimality and discuss specific applications to stochastic block models and nonparametric regression with sub-Gaussian errors.

Variational Bayes
Adaptive inference
Posterior contraction rates
Model selection consistency
Deep neural networks
Particle-based variational inference

§ INTRODUCTION

Bias-variance trade-off, a fundamental principle of statistical inference, suggests that the structure of a probability model needs to be appropriately chosen in order to efficiently gain useful information from data. For example, in high-dimensional linear regression, it is almost necessary to assume sparsity of the regression coefficients in order to effectively reduce the dimension of the data. Another example is nonparametric regression/density estimation, where complexity of a model needs to be selected or specified according to smoothness of the true regression function/density in order to optimize the prediction or estimation risk. To determine which model structure is appropriate requires knowledge of the underlying data-generating process. However, this knowledge is rarely available in practice, so accurate inference requires the statistical method to be adaptive, i.e., to be able to learn from data the relevant features of the data-generating process so that appropriate model structure determinations can be made.

Bayesian inference is appealing from this adaptive inference viewpoint, because it has the ability to automatically learn the model structure through a hierarchical design of the prior distribution, which consists of a distribution on the model structure and a conditional distribution of the parameter given the model. Indeed, a number of studies [Belitser and Ghosal, 2003, Lember and van der Vaart, 2007, Ghosal et al., 2008, Yang and Pati, 2017, Gao et al., 2020, Han, 2021] have shown that with carefully designed hierarchical prior distributions, one can conduct optimal inference adaptively via the  (original) posterior distributions. However, posterior inference based on sampling the posterior distribution over models with varying structure is usually computationally demanding and often very inefficient, especially when moves between models are poorly proposed in a Markov chain Monte Carlo (MCMC) algorithm.

High-dimensional and/or huge data sets are common in practical applications, so instead of using algorithms that exactly approximate original posterior distributions such as MCMC, variational Bayes that offers fast approximation of posterior distributions, has been  widely used as an attractive alternative. Along with the increasing popularity of variational Bayes, theory for contraction properties of variational posteriors has emerged recently. [Alquier and Ridgway, 2020] investigated theoretical properties of the variational approximation to the fractional posterior where the likelihood function is replaced by its fractional version with  a positive exponent less than 1.  Pati et al., 2018 and Zhang and Gao, 2020 extended the result of [Alquier and Ridgway, 2020] to the variational approximation to the original posterior. Yang et al., 2020 studied the variational posterior for models with latent variables, which differentiates it from the related work. These  aforementioned studies   provided general theoretical frameworks for characterizing  contraction rates of the variational posterior or its fractional version and proved their optimality in many applications. Other works studied asymptotic properties of a variational posterior distribution for a specific statistical model, for example, mixture models [Chérief-Abdellatif and Alquier, 2018], sparse (Gaussian) linear regression  [Ray and Szabó, 2021, Yang and Martin, 2020], sparse logistic linear regression [Ray et al., 2020] and sparse principal component analysis [Ning, 2021].

However, adaptivity of the variational posteriors remains an important and largely open problem. Although the existing general theoretical frameworks  can lead to adaptive inferences in some specific examples, they require specially designed prior distributions to achieve adaptivity, which limits applicability of the theory. Another problem concerns about  computation. Construction of a computationally tractable variational Bayes method over a collection of models with different structure is not an easy task in general. Exceptionally, Zhang and Gao, 2020 provided a potentially adaptive variational Bayes procedure with the model selection step, which we call model selection variational Bayes. They showed that after selecting the best model among the collection of the models with different structure and conducting variational inference over the selected model, the resulting variational posterior has optimality under mild assumptions. It seems that this procedure can attain adaptive optimality in many problems, which however was not elaborated. Furthermore, their variational posterior does not fully preserve information of the original posterior distribution, in particular, does not provide posterior probabilities over the model structure.

In this paper, we propose a new variational Bayes method or framework, called adaptive variational Bayes. Instead of selecting the best model and using the individual variational posterior over that model for consequent inference tasks, the proposed framework aggregates all the individual variational posteriors with certain weights to obtain a variational posterior over the entire model. By doing so, the adaptive variational Bayes provides a posterior probability of each model unlike the model selection variational Bayes. Moreover, the adaptive variational Bayes can attain optimal contraction rates adaptively for a wide range of statistical problems under mild conditions on the prior and the variational family. To the best of our knowledge, our framework is the first general recipe for showing adaptive contraction of variational posteriors.

We summarized our contributions as follows.

* Computational tractability.  We show that the variational posterior which is obtained by minimizing the Kullback-Leibler divergence to the original posterior distribution over the entire model can be represented by a certain “average” of the variational posteriors over individual models. This enables us to avoid the optimization over the entire model, for which a tractable algorithm may not exist. This implies that, our proposed variational Bayes is computationally efficient as  variational inference methods for individual models are, while the computation of the original posterior distributions with hierarchical priors may not be easy and often intractable as we mentioned previously.

* Adaptive contraction rate and model selection consistency. We formulate very mild conditions under which the variational posterior in our adaptive variational Bayes framework can attain optimal contraction rates adaptively. We also provide some easily verifiable sufficient conditions for our theoretical assumptions. Moreover, when the “true” model structure exists, we show that  model selection consistency can be achieved by the proposed variational posterior. We quantify the trade-off between exact model selection consistency and optimal contraction, which is represented as an amount of the hurt in the contraction rate.

* Adaptive variational deep learning. We apply our general theory to (deep) neural network models and derive some new  and adaptive results. In the first example, we consider nonparametric regression with neural networks. We prove that the adaptive variational posterior over “non-sparse” neural networks with various depth and width attains the optimal contraction rates adaptively in some regression problems. To our best knowledge, this is the first result on adaptive inference with non-sparse neural networks. The second example considers binary classification problem using neural networks and we derive the convergence rate for estimation of the conditional probability function. The last example is about intensity estimation for Poisson point processes, which is completely new  in the deep learning literature.

* Theoretical and empirical superiority over the  model selection variational Bayes.  We compare the proposed adaptive variational Bayes with the  variational posterior with  model selection  by [Zhang and Gao, 2020]. Theoretically, we prove that the adaptive variational Bayes always gives a tighter risk bound for the variational posterior than that of the model selection variational Bayes. We also conduct numerical experiments that show superiority of the adaptive variational Bayes.

* Particle-based construction of prior and variational family. We propose a general particle-based approach for a construction of a prior distribution and variational family. The particle-based approach simplifies theoretical conditions as well as provides an optimization algorithm applicable in general situations.

* Theoretical justification of the use of quasi-likelihood. We consider the use of quasi-likelihood in our framework. We formulate conditions on quasi-likelihood to ensure the contraction rates remain the same. We provide two examples of using a variational posterior based on the quasi-likelihood. The first example considers stochastic block models, where the sample is given as a set of connections between nodes of a graph, which are encoded as 0 or 1. We propose to use Gaussian quasi-likelihood which solves a problematic behavior of the Bernoulli likelihood. The second example is regression with sub-Gaussian errors. We show that when we use the Gaussian likelihood with sufficiently large variance parameter, we can still attain optimal contraction of the variational posterior without specifying the error distribution exactly.

The rest of the paper is organized as follows. We first introduce some notation in the rest of this section.  In <ref>, we develop a general framework for adaptive variational Bayes inference. Concretely, we provide a design of the prior distribution and variational family, and describe a general and simple scheme for computation of the proposed variational posterior. In <ref>, we study contraction properties of the proposed variational posterior including  oracle rates, adaptivity and model selection consistency.  In <ref>, we apply our general results to variational deep learning for a number of statistical problems. In <ref>, we compare our adaptive variational posterior with the  model selection variational posterior proposed  by [Zhang and Gao, 2020] both theoretically and empirically. In <ref>, we propose a particle-based approach for a construction of a prior distribution and variational family, which simplifies theoretical conditions to derive contraction properties. In <ref>, we study theoretical properties of variational quasi-posteriors where the usual likelihood function is replaced with an alternative quasi-likelihood.

§.§ Notation

Let $\R$, $\R_+$, $\R_{\ge0}$, $\bZ$, $\bN_0$ and $\bN$  be the sets of real numbers, positive numbers, nonnegative numbers, integers, nonnegative integers and natural numbers, respectively. We denote by $\ind(\cdot)$ the indicator function.  For two integers $z_1,z_2\in\bZ$ with $z_1\le z_2$, we let $[z_1:z_2]:=\cbr{z\in\bZ:z_1\le z\le z_2}$ and if $z_2\ge2$, we use the shorthand $[z_2]:=[1:z_2]$ and $\bN_{\ge z_2}:=\bN\setminus[z_2-1]=\cbr{z\in\bN:z\ge z_2 }$. For $d\in\bN$, let $\zero_d$ and $\one_d$ denote the $d$-dimensional vectors of 0's and of 1's, respectively and let $\I_d$ denote the $d\times d$-dimensional identity matrix. For a $d$-dimensional vector $\x:=(x_1,\dots, x_d)\in\R^d$, we denote $|\x|_p:=\del[0]{\sum_{j=1}^d|x_j|^p}^{1/p}$ for  $p\ge1$, which are the usual Euclidean $p$-norm, and denote $|\x|_0:=\sum_{j=1}^d\ind(x_j\neq0)$ and $|\x|_\infty:=\max_{j\in[d]}|x_j|$.  Let $\Delta_d:=\cbr[0]{\balpha\in[0,1]^d:|\balpha|_1=1}$. For a real number $x\in\R$, we define $\floor{x}:=\max\{z\in\bZ:z\le x\}$,  and $\ceil{x}:=\min\{z\in\bZ:z\ge x\}$. For two real numbers $a,b\in\R$, we write $a\vee b:=\max\{a,b\}$ and $a\wedge b:=\min\{a,b\}$. Moreover, for two real vectors  $\a:=(a_j)_{j\in[d]}\in\R^d$ and $\b:=(b_j)_{j\in[d]}\in\R^d$ with the same dimension, we write $\a\vee\b:=(a_j\vee b_j)_{j\in[d]}$ and $\a\wedge\b:=(a_j\wedge b_j)_{j\in[d]}$. For an arbitrary set $B$, we denote by $B^\complement$ its complement and by $|B|$ its cardinality.  For two positive sequences $(a_n)_{n\in \mathbb{N}}$ and $(b_n)_{n\in \mathbb{N}}$, we write $a_n\lesssim b_n$ or  $b_n\gtrsim a_n$, if there exists a positive constant $C>0$ such that $a_n\le Cb_n$ for any $n\in \mathbb{N}$. Moreover, we write $a_n\asymp b_n$ if both $a_n\lesssim b_n$  and  $a_n\gtrsim b_n$ hold. We use the little $\sco$ notation, that is, we write $a_n=\sco(b_n)$ if $\lim_{n\to\infty}a_n/b_n=0.$

For a measurable space $(\fX,\cX)$, we denote by $\bP(\fX)$ the set of all probability measures supported on $\fX$. For a probability measure $\P\in\bP(\fX)$ and a function $g$ on $\fX$, we write $\P g:=\int_{\fX} g\d\P$, i,e., $\P g$ denotes the expectation of $g$ with respect to the measure $\P$. For two probability distributions $\P_1\in\bP(\fX)$ and $\P_2\in\bP(\fX)$, we denote by $\kl(\P_1, \P_2)$ the Kullback-Leibler (KL) divergence from $\P_2$ to $\P_1$, which is defined by

\begin{equation*}
\kl(\P_1,\P_2)
\int \log(\frac{\d\P_1}{\d\P_2})\d\P_1 & \mbox{ if }\P_1\ll\P_2,\\
+\infty & \mbox{ otherwise}.
\end{cases}
\end{equation*}

For simplicity, we slightly abuse a notation to denote $\kl(\balpha_1,\balpha_2):=\kl(\Cat(\balpha_1), \Cat(\balpha_2))$ for any $\balpha_1,\balpha_2\in\Delta_m$ and any $m\in\bN$, where $\Cat(\balpha)$ stands for the categorical distribution with probability vector $\balpha.$ For $\rho\in(0,1)\cup(1,\infty)$, the $\rho$-Rényi divergence from $\P_2$ to $\P_1$ is defined as

\begin{equation*}
\sD_\rho(\P_1,\P_2)
\frac{1}{\rho-1}\log(\int \del{\frac{\d\P_1}{\d\P_2}}^{\rho-1}\d\P_1)& \mbox{ if }\P_1\ll\P_2,\\
+\infty & \mbox{ otherwise}.
\end{cases}
\end{equation*}

Note that $\sD_\rho(\P_1,\P_2)$ is a nondecreasing function of $\rho$ and  $\sD_{\rho_1}(\P_1,\P_2)\le \kl(\P_1,\P_2)\le \sD_{\rho_2}(\P_1,\P_2)$ for any $\rho_1\in(0,1)$ and $\rho_2\in(1,\infty)$. Let $\|\P_1-\P_2\|_{\textsc{TV}}:=\sup_{B\in\cX}|\P_1(B)-\P_2(B)|$ be the total variation distance between $\P_1$ and $\P_2$. For $x\in\fX$, let $\delta(\cdot;x)$ be a Dirac-delta measure at $x$ such that $\delta(B;x)=1$ if $x\in B$ and $\delta(B;x)=0$ otherwise.

For a metric space $(\Theta,\scd)$, let $\cB_{\Theta,\scd}(\btheta_0,\zeta):=\cbr{\btheta\in\Theta:\scd(\btheta,\btheta_0)\le \zeta}$ be a $\scd$-(closed) ball of radius $\zeta>0$ centered at $\btheta_0\in\Theta$. For $\zeta>0$, we say that $\Theta^\bullet$ is a $\zeta$-net of $\Theta$ (with respect to (w.r.t.) the metric $\scd$), if $\Theta\subset\cup_{\btheta^\bullet\in\Theta^\bullet}\cB_{\Theta,\scd}(\btheta^\bullet,\zeta)$. We define the $\zeta$-covering number of $\Theta$ with respect to the metric $\scd$ by $\scN(\zeta, \Theta,\scd):=\min\cbr{|\Theta^\bullet|:\mbox{$\Theta^\bullet$ is a $\zeta$-net of $\Theta$ w.r.t. $\scd$}}$.
Furthermore, we denote a projection of an element $\btheta\in\Theta$ to a subset $\Theta'\subset\Theta$ by $ \proj_{\Theta',\scd}(\btheta)$, i.e., $\proj_{\Theta',\scd}(\btheta)\in\argmin_{\btheta'\in\Theta'}\scd(\btheta',\btheta).$

§ ADAPTIVE VARIATIONAL BAYES

In this section, we develop  a novel variational Bayes framework called adaptive variational Bayes. We first describe a statistical setup, then introduce an estimation procedure in the proposed framework as well as a general and simple scheme for its computation.

§.§ Statistical experiment and models

We describe our setup for a statistical experiment, which, in this paper, is defined as a pair of a sample space and a set of some distributions on the sample space. For each sample size $n\in\bN$,  suppose that we observe the $\fY_n$-valued sample $\Y^{(n)}$, where $\fY_n$ is a measurable sample space equipped with a reference $\sigma$-finite measure $\mu_n$.  We then model the sample as having distribution $\P_{\blambda}^{(n)}\in\bP(\fY_n)$ determined by a natural parameter $\blambda$ in a measurable natural parameter space $\Lambda_n$. The natural parameter space can be infinite-dimensional, for
instance, see <ref> below, where the natural parameter space is given as the space of density functions. We assume that there exists  nonnegative function $\sq_n:\Lambda_n\times \fY_n\mapsto\R_{\ge0}$ such that, for every parameter $\blambda\in\Lambda_n$,  $\int_{\fY_n}\sq_n(\blambda,\y^{(n)})\d\mu_n(\y^{(n)})=1$ and

\begin{equation}
\label{eq:likelihood}
\P_{\blambda}^{(n)}(\d\y^{(n)})=\sq_n(\blambda, \y^{(n)}) \mu_n(\d\y^{(n)})
\end{equation}

for any $\y^{(n)}\in\fY_n$. The function $\sq_n$ is called a likelihood function or simply likelihood. We denote by $\bP(\fY_n;\sq_n,\Lambda_n)$ the set of all distributions of the form eq:likelihood.

With a countable set of model indices $\cM_n$, we consider a collection of probability models $\{\bM_{n,m}\}_{m\in\cM_n}$  for estimation based on the sample $\Y^{(n)}$, where each model is of the form

\begin{equation*}
\bM_{n,m}:=\cbr{\P_{\sT_n(\btheta)}^{(n)}\in\bP(\fY_n;\sq_n,\Lambda_n):\btheta\in\Theta_{n,m}}
\end{equation*}

with a parameter space $\Theta_{n,m}$ and a measurable map $\sT_n:\cup_{m\in\cM_n}\Theta_{n,m}\mapsto \Lambda_n$.  We call the map $\sT_n$ the natural parameterization map. We introduce additional notations which will be useful in the rest of the paper. We denote $\Theta_{n,\cM_n}:=\cup_{m\in\cM_n}\Theta_{n,m}$. For any subset $\Theta'\subset\Theta_{n,\cM_n}$, we write $\sT_n(\Theta'):=\cbr{\sT_n(\btheta'):\btheta'\in\Theta'}$.

In order to obtain computational tractability (see <ref>), we assume that the parameter spaces $\{\Theta_{n,m}\}_{m\in\cM_n}$ are disjoint for every $n\in\bN$. One can easily construct disjoint parameter spaces in many cases. For example, a model selection problem for a nested model such that the corresponding parameter spaces $\{\Theta_{n,\le m}\}_{m\in\bN}$ satisfy $\Theta_{n,\le m}\subset\Theta_{n,\le m+1}$ for every $m\in\bN$, can be formulated using the disjoint parameter spaces $\Theta_{n,m}:=\Theta_{n,\le m+1}\setminus\Theta_{n,\le m}$ for $m\in\bN$ with $\Theta_{n,\le0}:=\emptyset$. Furthermore, the natural parametrization map can mitigate any restriction of the disjointness assumption. One may construct a collection of disjoint parameter spaces, which possibly yields non-disjoint natural parameter spaces such as $\sT_n(\Theta_{n,m})\cap\sT_n(\Theta_{n,m'})\neq \emptyset$ for some $m\neq m'$.

We provide a few  examples of a collection of probability models with disjoint parameter spaces.

[Factor model]

Estimating the high-dimensional covariance matrix given a $p_n$-dimensional iid normal sample $\Y^{(n)}:=(\Y_1,\dots, \Y_n)\in\fY_n:=(\R^{p_n})^{\otimes n}$ has attracted a lot of attention, which has a distribution given as

\begin{equation*}
(\Y_1,\dots, \Y_n)\sim \P_{\bSigma}^{(n)}:=\bigotimes_{i=1}^n\N\del{\zero_{p_n}, \bSigma}
\end{equation*}

for $\bSigma\in\Lambda_n:=\textup{SPD}_{p_n}$, where $\textup{SPD}_{p_n}$ denotes the set of $p_n \times p_n$ symmetric positive definite matrices. One of the popular models for modeling the high-dimensional covariance matrix is the factor model. Let $\Theta_{n,m}:=\R^{p_n\times m}$ for $m\in\bN$ and let $\sT_n:\cup_{m\in\bN}\Theta_{n,m}\mapsto\textup{SPD}_{p_n}$ be a map defined as

\begin{equation*}
\sT_n(\bPhi)=\bPhi\bPhi^\top+\I_{p_n}
\end{equation*}

for $\bPhi\in\cup_{m\in\bN}\Theta_{n,m}$. The covariance matrix $\sT_n(\bPhi)$ has a spiked eigen spectrum, i.e., low rank structure, when the rank of the “factor loading matrix” $\bPhi$ is less than $p_n.$ The factor model is the collection of the models  $\{\bM_{n,m}\}_{m\in[p_n-1]}$ with

\begin{equation*}
\bM_{n,m}:=\cbr{\P_{\sT_n(\bPhi)}^{(n)}:=\bigotimes_{i=1}^n\N\del{\zero_{p_n}, \sT_n(\bPhi)}:\bPhi\in\Theta_{n,m}},
\end{equation*}

where the parameter spaces $\{\Theta_{n,m}\}_{m\in[p_n-1]}$ are disjoint.

[Gaussian location mixture model]

Consider a statistical experiment in which a $d$-dimensional real-valued sample $\Y^{(n)}:=(\Y_1,\dots, \Y_n)\in\fY_n:=(\R^d)^{\otimes n}$ is observed. Let $\cP^d$ denote the set of all probability density functions on $\R^d$. Let $\sq_n:\cP^d\times(\R^d)^{\otimes n}\mapsto\R_{\ge0}$ be a likelihood function such that

\begin{equation*}
\sq_n(\sp, \Y^{(n)})=\prod_{i=1}^n\sp(\Y_i)
\end{equation*}

for $\sp\in\cP^d$ and $\Y^{(n)}\in(\R^d)^{\otimes n}$. For a density $\sp\in\cP^d$, we denote by $\P_{\sp}^{(n)}$ the distribution of the sample $\Y^{(n)}$ with the probability density function $\sq_n(\sp,\cdot)$. The Gaussian location mixture model for modeling the sample $\Y^{(n)}$ is given as follows. Given the number of components $m\in\bN$, a vector of mixing weights $\bvarpi:=(\varpi_1,\dots, \varpi_m)\in\Delta_m$ and  a vector of location parameters $\bvartheta_{1:m}:=(\bvartheta_1,\dots, \bvartheta_m)\in (\R^d)^{\otimes m}$, we let $\sT:\cup_{m\in\bN}(\Delta_m\times(\R^d)^{\otimes m})\mapsto \cP^d$ be a map defined as

\begin{equation*}
\sT(\bvarpi, \bvartheta_{1:m})=\sum_{k=1}^m\varpi_k\sp_{\N(\bvartheta_k,\I_d)},
\end{equation*}

where $\sp_{\N(\bvartheta,\I_d)}$ stands for the probability density function of the multivariate normal distribution with mean $\bvartheta$ and covariance matrix $\I_d$. The Gaussian location mixture model is the collection of the models  $\{\bM_{n,m}\}_{m\in\bN}$ with

\begin{equation*}
\bM_{n,m}:=\cbr{\P_{\sT(\bvarpi,\bvartheta_{1:m})}^{(n)}:(\bvarpi,\bvartheta_{1:m})\in\Theta_{n,m}:=\Delta_m\times(\R^d)^{\otimes m}}.
\end{equation*}

The parameter spaces $\{\Theta_{n,m}\}_{m\in\bN}$ are disjoint.

[Stochastic block model]

Suppose that we observe a sample $(Y_{i,j})_{(i,j)\in[n]^2:i>j}$ from a graph of $n$ nodes, where $Y_{i,j}\in\{0,1\}$ indicates whether there is a connection between the $i$ and $j$-th nodes. For $m\in[n]$, let

\begin{align}
\cS^m&:=\cS^m([0,1]):=\cbr{\U\in [0,1]^{m\times m}:\U=\U^\top}, \label{eq:def_sm}\\
\cZ_{n, m}&:=\cbr{\Z=(\z_1,\dots,\z_n)^\top\in\{0,1\}^{n\times m}: |\z_i|_1=1}.  \label{eq:def_znm}
\end{align}

The stochastic block model assumes the connectivities between $n$ nodes follow the distribution

\begin{equation*}
(Y_{i,j})_{(i,j)\in[n]^2:i>j}\sim \P_{\sT_n(\U, \Z)}^{(n)}:=\bigotimes_{i=2}^{n}\bigotimes_{j=1}^{i-1}\Ber(\Omega_{i,j})
\mbox{ with }\Omega_{i,j}:=\z_i^\top\U\z_j
\end{equation*}

for $\U:=(U_{k,h})_{k\in[m], h\in[m]}\in\cS^m$ and $\Z:=(\z_1,\dots, \z_n)^\top\in\cZ_{n,m}$ for some $m\in[n]$, where $\sT_n$ is a map defined by $\sT_n(\U,\Z)=(\z_i^\top\U\z_j)_{(i,j)\in[n]^2:i>j}$. Here $\Ber(\Omega)$ denotes the Bernoulli distribution with probability $\Omega\in[0,1]$. Let $\Theta_{n,m}:=\cS^m\times\cZ_{n,m}$, then the parameter spaces $\{\Theta_{n,m}\}_{m\in[n]}$ are disjoint.

[Sparse linear regression model]
Consider a real-valued sample $\Y^{(n)}:=(Y_1,\dots, Y_n)\in\fY_n:=\R^n$ whose distribution is thought to be

\begin{align*}
(Y_1,\dots, Y_n)\sim\P_{ \bvartheta}^{(n)}:=\bigotimes_{i=1}^n\N\del{\vartheta_i, 1}
\end{align*}

for some “mean vector” $\bvartheta:=(\vartheta_1,\dots,\vartheta_n)\in\Lambda_n:=\R^n$. Suppose that  we have a fixed $p_n$-dimensional covariate $\x_i$  for each observation $Y_i$ and let $\X:=(\x_1,\dots, \x_n)^\top\in\R^{n\times p_n}$. We model the mean vector as $\bvartheta=\X\btheta$ for some (possibly sparse) “regression coefficeint” $\btheta\in\R^{p_n}$. Given a set of indices $S\subset[p_n]$, we denote by $\Theta_{n,S}\subset\R^{p_n}$ the set of sparse regression coefficients with support $S$, that is,

\begin{equation*}
\Theta_{n,S}:=\cbr{(\theta_j)_{j\in[p_n]}\in\R^{p_n}:\theta_j\neq0 \mbox{ if $j\in S$ and $\theta_j=0$ otherwise}}.
\end{equation*}

Then the sets $\{\Theta_{n,S}\}_{S\subset[p_n]}$ are disjoint. The sparse linear regression model is the collection of models $\{\bM_{n,S}\}_{S\subset [p_n]}$ with $\bM_{n,S}:=\cbr[0]{\P_{\X\btheta}^{(n)}:\btheta\in\Theta_{n,S}}$.

§.§ Prior and variational posterior

For each $n\in\bN$, we impose a prior distribution $\Pi_n$ on  the parameter space $\Theta_{n,\cM_n}:=\cup_{m\in\cM_n}\Theta_{n,m}$, which is of the form

\begin{equation}
\label{eq:prior}
\Pi_n=\sum_{m\in\cM_n}\alpha_{n,m}\Pi_{n,m},
\end{equation}

where $\balpha_n:=(\alpha_{n,m})_{m\in\cM_n}\in \Delta_{|\cM_n|}$ and $\Pi_{n,m}\in\bP(\Theta_{n,m})$ for each $m\in\cM_n$. Note that $\Pi_n(\btheta\in\Theta_{n,m})=\sum_{m\in\cM_n}\alpha_{n,m}\Pi_{n,m}(\btheta\in\Theta_{n,m})=\alpha_{n,m}$, that is, $\alpha_{n,m}$ can be interpreted as the (prior) probability that the parameter $\btheta$ belongs to the space $\Theta_{n,m}$. Then, with the likelihood function $\sq_n:\Lambda_n\times \fY_n\mapsto \R_{\ge0}$ and the natural parameterization map $\sT_n:\Theta_{n,\cM_n}\mapsto \Lambda_n$, the posterior distribution given the sample $\Y^{(n)}$ is defined as

\begin{equation}
\Pi_n(\btheta\in\cT|\Y^{(n)})=\frac{\int_{\cT} \sq_n(\sT_n(\btheta), \Y^{(n)})\d\Pi_n(\btheta)}{\int \sq_n(\sT_n(\btheta), \Y^{(n)})\d\Pi_n(\btheta) }
\end{equation}

for any measurable subset $\cT\subset \Theta_{n,\cM_n}.$

Throughout this paper, we focus on the variational family given by

\begin{equation}
\label{eq:vfamily}
\bV_n:=\cbr{\sum_{m\in\cM_n}\gamma_{n,m}\Xi_{n,m}:(\gamma_{n,m})_{m\in\cM_n}\in \Delta_{|\cM_n|}, \Xi_{n,m}\in\bV_{n,m}}
\end{equation}

for some $\bV_{n,m}\subset\bP(\Theta_{n,m})$ for $m\in\cM_n$.
Then we aim to find the closest distribution, denoted by $\hXi_n$, to the posterior distribution $\Pi_n(\cdot|\Y^{(n)})$ among the variational family $\bV_n$ in a sense that

We will show  in the next section that under appropriate conditions, the  distribution $\hXi_n$ can obtain optimal contraction rates adaptively, we therefore call it the adaptive variational posterior distribution.

We close this subsection by introducing some additional notations and well-known facts. Let $\sq_{n,\Pi}(\Y^{(n)}):=\int\sq_n(\sT_n(\btheta), \Y^{(n)})\d\Pi(\btheta)$  be the marginal likelihood of the sample $\Y^{(n)}$ with respect to the prior distribution $\Pi$. We define

\begin{equation}
\scE(\Xi,\Pi,\sq_n,\Y^{(n)}):=-\int\log\sq_n(\sT_n(\btheta), \Y^{(n)})\d\Xi(\btheta)+\kl(\Xi, \Pi)
\end{equation}

for any distributions $\Xi,\Pi\in\bP(\Theta_{n,\cM_n})$, where we suppress the dependence on the natural parametrization map $\sT_n$ for simplicity. Then by simple algebra, it can be shown that the minimization problem in eq:vbayes is equivalent to

\begin{equation}
\hXi_n\in\argmin_{\Xi\in\bV_n}\scE(\Xi,\Pi_n,\sq_n,\Y^{(n)}).
\end{equation}

Thus we call the function $\scE(\cdot,\cdot,\cdot,\cdot)$ the variational objective function. The negative of the variational objective function is called an evidence lower bound (ELBO) in the related literature, which is named after the property that the ELBO is a lower bound of the log marginal likelihood, i.e., $\log\sq_{n,\Pi}(\Y^{(n)})\ge-\scE(\Xi,\Pi,\sq_n,\Y^{(n)})$.

§.§ A general  and efficient scheme for computation

The minimization program in eq:vbayes  needs to be done over distributions on the whole parameter space $\Theta_{n,\cM_n}$ which is usually quite `big'  and even has varying dimensions and/or structures. Therefore this seems to be computationally intractable. However, as stated in the following theorem, if the parameter spaces are disjoint as we have assumed, the minimization program can be done by combining the results of the minimization over individual models.

Suppose that the parameter spaces $\{\Theta_{n,m}\}_{m\in\cM_n}$ are disjoint. Then the variational posterior distribution in eq:vbayes satisfies

\begin{equation}
\label{eq:vposterior}
\hXi_n=\sum_{m\in\cM_n}\hgamma_{n,m}\hXi_{n,m},
\end{equation}

\begin{equation}
\label{eq:vpost_comp}
\hXi_{n,m}\in\argmin_{\Xi\in\bV_{n,m}}\scE\del[1]{\Xi,\Pi_{n,m},\sq_n,\Y^{(n)}}
\end{equation}

\begin{equation}
\label{eq:vpost_mprob}
\hgamma_{n,m}=\frac{1}{\hat{Z}_{\gamma,n}} \alpha_{n,m}\exp\del[1]{-\scE\del[1]{\hXi_{n,m},\Pi_{n,m},\sq_n,\Y^{(n)}}}
\end{equation}

with $\hat{Z}_{\gamma,n}:=\sum_{m'\in\cM_n}\alpha_{n,m'}\exp\del[1]{-\scE(\hXi_{n,m'},\Pi_{n,m'},\sq_n,\Y^{(n)})}$ the normalizing constant.

Due to the above theorem, if a computation algorithm for computing the variational posterior for each individual model is tractable, which is the case for many situations, the adaptive variational posterior distribution is also easily computable. We summarize a general procedure for computation in <ref>.

Adaptive variational Bayes

Sample $\Y^{(n)}$, prior distribution $\Pi_n=\sum_{m\in\cM_n}\alpha_{n,m}\Pi_{n,m}$, variational family $\bV_n$.
The variational posterior distribution $\hXi_n$.
Find  $\hXi_{n,m}\in\argmin_{\Xi\in\bV_{n,m}}\scE(\Xi,\Pi_{n,m},\sq_n,\Y^{(n)})$.

Compute $\tgamma_{n,m}=\alpha_{n,m}\exp\del[0]{-\scE(\hXi_{n,m},\Pi_{n,m},\sq_n,\Y^{(n)})}$.

Compute $ \hgamma_{n,m}=\tgamma_{n,m}/(\sum_{m'\in\cM_n}\tgamma_{n,m'})$ for $m\in\cM_n$.

§.§ Posterior model probabilities

In this short subsection, we compare the quantified posterior model probabilities by the adaptive variational posterior and the original posterior $\Pi_n(\cdot|\Y^{(n)})$.
Under the prior distribution of eq:prior, the (original)  posterior distribution on $\Theta_{n,\cM_n}$ can be written as

\begin{equation}
\Pi_n(\cdot|\Y^{(n)})=\sum_{m\in\cM_n}\halpha_{n,m}\Pi_{n,m}(\cdot|\Y^{(n)}),
\end{equation}

where $\halpha_{n,m}:=\Pi_n(\btheta\in\Theta_{n,m}|\Y^{(n)})$, which is the (original) posterior probability of the model $m$, and $\Pi_{n,m}(\cdot|\Y^{(n)}):=\Pi_n(\cdot|\btheta\in\Theta_{n,m}, \Y^{(n)})\in\bP(\Theta_{n,m})$ for $m\in\cM_n$.  Note that

=1/Ẑ_α,n α_n,m∫_Θ_n,m_n(_n(), ^(n))Π̣_n,m()

=1/Ẑ_α,n α_n,mexp[1]log_n,Π_n,m(^(n)),

where $\hat{Z}_{\alpha,n}:=\sum_{m'\in\cM_n}\alpha_{n,m'}\int_{\Theta_{n,m'}}\sq_n(\sT_n(\btheta), \Y^{(n)})\d\Pi_{n,m'}(\btheta)$ is the normalizing constant.

The next proposition states that if the variational approximation error is small then the resulting variational posterior model probabilities are  close to the original posterior model probabilities.

Under the same assumptions of <ref>,

\begin{align}
\del{\sum_{m\in\cM_n}\abs{\hgamma_{n,m}-\halpha_{n,m}}}^2
&\le  2\kl\del[1]{\hXi_n, \Pi_n(\cdot|\Y^{(n)})}.
\end{align}

§ CONCENTRATION PROPERTIES OF ADAPTIVE VARIATIONAL POSTERIOR

In this section, we provide adaptive posterior contraction rates and model selection consistency of the adaptive variational posterior distribution  eq:vposterior in general situations. We start with describing some  assumptions on the prior, model and variational family, most of which are very general and can be easily verified. Absolute constants denoted by e.g., $\fc_1,\fc_2,\dots$ may vary from place to place.

§.§ Assumptions

In this subsection, we provide a series of assumptions and several sufficient conditions for some parts of the assumptions, which are helpful in applying our general theoretical results to specific examples. We begin with a set of assumptions on the statistical experiment.

A sequence of statistical experiments $((\fY_n,\bP_n))_{n\in\bN}$, where each $\fY_n$ is a measurable space equipped with a reference $\sigma$-finite measure $\mu_n$ and $\bP_n$ is a subset of $\bP(\fY_n)$, satisfies the followings:

* (Parametrization and likelihood)
For any sufficiently large $n\in\bN$, $\bP_n=\bP(\fY_n;\sq_n,\Lambda_n)$ for a natural parameter space $\Lambda_n$ equipped with a metric $\scd_n:\Lambda_n\times\Lambda_n\mapsto \R_{\ge0}$ and a likelihood function $\sq_n:\Lambda_n\times \fY_n\mapsto\R_{+}$.

(Identifiability and tail behavior) There exist absolute constants $\rho_\circ\in(0,1)$, $\rho_\blacklozenge>1$,  $\fa_2\ge\fa_1\in\R$  and $\fc_2\ge\fc_1>0$ such that

\begin{align}
\fa_1+\fc_1n\scd_n^2(\blambda_0, \blambda_1)
&\le \sD_{\rho_\circ}\del[1]{\P^{(n)}_{\blambda_0}, \P^{(n)}_{\blambda_1}}, \label{eq:tail1}\\
\sD_{\rho_\blacklozenge}\del[1]{\P^{(n)}_{\blambda_0}, \P^{(n)}_{\blambda_1}}
&\le \fa_2+\fc_2n\scd_n^2(\blambda_0, \blambda_1)\label{eq:tail2}
\end{align}

for  any $\blambda_0, \blambda_1\in\Lambda_n$ and any sufficiently large $n\in\bN$.

The existence of the likelihood function assumed in <ref> is ubiquitous in Bayesian literature. Note that we do not assume the observations in the sample are iid or even independent to each other. That is, our framework can deal with a general dependence structure between observations, including times series data and network data. In the rest of the paper, we let $\scd_n:\Lambda_n\times\Lambda_n\mapsto \R_{\ge0}$  denote the metric that have appeared in <ref>.

<ref> in <ref>  ensures the identifiability of the parameter. Moreover, <ref> restricts tail behavior of the distribution of the sample. Concretely, this assumption implies that the two divergence measures $\sD_{\rho_\circ}(\cdot, \cdot)$ and $\sD_{\rho_\blacklozenge}(\cdot, \cdot)$ are equivalent, while in  general, the latter is always larger than or equal to the former. This equivalence can be obtained for distributions with light tails. One very important example is given as follows: <ref> is met if for every pair of two parameters, the corresponding log likelihood ratio function  is a sub-Gamma random variable, which was assumed in [Han, 2021] under the name of “local Gaussianity condition”.

Suppose that there exist absolute constants $\fc_1'>1$, $\fc_2'>0$, $\fc_3'>0$ and $\fc_4'>0$ such that

\begin{align*}
\frac{1}{\fc_1'}\scd_n^2(\blambda_0,\blambda_1)
\le \frac{1}{n}\kl\del[1]{\P^{(n)}_{\blambda_0}, \P^{(n)}_{\blambda_1}}\le \fc_1'\scd_n^2(\blambda_0,\blambda_1)
\end{align*}

and that

\begin{align*}
\P_{\blambda_0}^{(n)}\e^{t\del{\log\del{\frac{\sq_n(\blambda_0,\Y^{(n)})}{\sq_n(\blambda_1,\Y^{(n)})}} - \kl\del{\P^{(n)}_{\blambda_0}, \P^{(n)}_{\blambda_1}}}}
\le \fc_2'\e^{\fc_3't^2n\scd_n^2(\blambda_0,\blambda_1)}
\end{align*}

for any $t\in[-\fc_4',\fc_4'],$ any $\blambda_0, \blambda_1\in\Lambda_n$ and any sufficiently large $n\in\bN$. Then <ref> holds.

We next impose conditions on expressibility and complexity of the models.

Let $(\Lambda_n^\star)_{n\in\bN}$ be a sequence of natural parameter spaces such that $\Lambda_n^\star\subset\Lambda_n$ for every $n\in\bN$. For each $n\in\bN$, consider a collection of probability models $\{\bM_{n,m}\}_{m\in\cM_n}$ for some model index set $\cM_n$, where each model $\bM_{n,m}:=\cbr[0]{\P_{\sT_n(\btheta)}^{(n)}\in\bP_n:\btheta\in\Theta_{n,m}}$ is given with the parameter space $\Theta_{n,m}$ and the natural parametrization map $\sT_n:\Theta_{n,\cM_n}\mapsto \Lambda_n$. There exist collections of positive numbers $\{\eta_{n,m}:m\in\cM_n,n\in\bN\}$ and $\{\zeta_{n,m}:m\in\cM_n,n\in\bN\}$ which satisfy the followings:

(Expressibility)  For any sufficiently large $n\in\bN$, $m\in\cM_n$ and $\blambda^\star\in\Lambda_n^\star$, there exists a parameter   $\btheta_{n,m}^*(\blambda^\star)\in\Theta_{n,m}$  such that

\begin{equation}
\scd_n\del[2]{\sT_n\del[1]{\btheta_{n,m}^*(\blambda^\star)},\blambda^\star}\le \eta_{n,m}.
\end{equation}

(Complexity) For any $u>0$, there exists an absolute constant $\fc_1(u)>0$ depending on $u$ such that

\begin{align}
\sup_{\zeta>\zeta_{n,m}} \log\scN\del{u\zeta,\cbr{\btheta\in\Theta_{n,m}:\scd_n(\sT_n(\btheta),\blambda^\star)\le 2\zeta}, \scd_n}\le \fc_1(u) n\zeta_{n,m}^2
\end{align}

for any $\blambda^\star\in\Lambda_n^\star$, $m\in\cM_n$ and sufficiently large $n\in\bN$.

In the rest of the paper, we interpret $(\Lambda_n^\star)_{n\in\bN}$ appearing in <ref> as a sequence of true natural parameter spaces, that is, the true  distribution of the sample $\Y^{(n)}$ is assumed to belong to the true model $\bP^\star_n:=\cbr[0]{\P_{\blambda^\star}^{(n)}:\blambda^\star\in\Lambda_n^\star}$. Note that the true model does not need to be contained in the model $\cup_{m\in\cM_n}\bM_{n,m}$ used for estimation. In <ref>, we assume that the “distance”  from the model $\bM_{n,m}$ to the true model $\bP^\star_n$ is at most  $\eta_{n,m}$ that we will call  the approximation error of the model $\bM_{n,m}$ in the rest of the paper. The parameter $\btheta_{n,m}^*(\blambda^\star)\in\Theta_{n,m}$ can be viewed as the best approximator of $\blambda^\star$ in $\Theta_{n.m}$. On the other hand, in <ref>, it is assumed that the local complexity of the model $\bM_{n,m}$ is bounded above by  $\zeta_{n,m}$ which we call the  estimation error of the model $\bM_{n,m}$ since the local complexity determines the estimation variance.

Note that, for the theoretical purpose, we do not need to assume that the parameter spaces are disjoint.

We need to impose some technical assumptions on behavior of the errors.

[Error rates]

The collections of approximation errors $\{\eta_{n,m}:m\in\cM_n,n\in\bN\}$ and estimation errors $\{\zeta_{n,m}:m\in\cM_n,n\in\bN\}$ satisfy the followings:

*   The oracle rate $(\epsilon_n)_{n\in\bN}$ defined as

\begin{equation}
\label{eq:oracle_rate}
\epsilon_n:=\inf_{m\in\cM_n}(\eta_{n,m}+\zeta_{n,m}),
\end{equation}

satisfies  that $n\epsilon_n^2\to\infty$ as $n\to\infty$.

There exist absolute constants $\fc_1>0$, $\fc_2>0$ and $\fc_3>0$ such that

\begin{align}
\sum_{m\in\cM_n:\zeta_{n,m}\ge H'\epsilon_n}\exp\del[1]{-L'n\zeta_{n,m}^2}&\le \exp\del[1]{(\fc_1-L'(H')^2)n\epsilon_n^2}, \\
\sum_{m\in\cM_n:\zeta_{n,m}< H'\epsilon_n}\exp\del[1]{L''n\zeta_{n,m}^2}&\le \exp\del[1]{(\fc_2+\fc_3L''(H')^2)n\epsilon_n^2}
\end{align}

for any sufficiently large $n$ and any $H'>1$, $L'>0$ and $L''>0$.

The sequence $(\epsilon_n)_{n\in\bN}$ is called an oracle rate since this is the contraction rate that is achievable when the best model is known. We assume that the oracle rate is at least $b_n/\sqrt{n}$ with $b_n\to\infty$ for technical reasons. Most of the Bayesian literature assumed the same type of conditions. The second condition is introduced to control complexities of some unions of the parameter spaces. We provide two sufficient conditions of <ref> given in the following two lemmas, which are useful in applications of our general results.

Suppose that  there exist absolute constants $\fc'_1>0$ and $\fc'_2>0$ such that $|\cM_n|\le \fc'_1 \exp\del[1]{\fc_2' n\epsilon_n^2}$ for any $n\in\bN$. Then <ref> holds.

Suppose that $\cM_n\subset\bN$ and $\zeta_{n,m}=m^{r/2}a_n$ for any $n,m\in\bN$ for some constant $r\ge1$ and some sequence $(a_n)_{n\in\bN}$ such that $na_n^2\to\infty$. Then  <ref> holds.

We consider a prior distribution that satisfies the following assumption.

[Prior distribution]

A sequence of prior distributions $(\Pi_n)_{n\in\bN}\in\bigotimes_{n\in\bN}\bP(\Theta_{n,\cM_n})$, where each $\Pi_n$ is of the form $\Pi_n=\sum_{m\in\cM_n}\alpha_{n,m}\Pi_{n,m}$ with $\balpha_n:=(\alpha_{n,m})_{m\in\cM_n}\in\Delta_{|\cM_n|}$ and $\Pi_{n,m}\in\bP(\Theta_{n,m})$ for every $m\in\bM_n$, satisfies the followings:

(Model probability) For each $n\in\bN$, the prior model probability $\balpha_n$ is given by

\begin{align}
\alpha_{n,m} =\frac{1}{Z_{\alpha,n}}\e^{-b_0L_nn\zeta_{n,m}^2}
\end{align}

for  $L_n\ge1$ and $b_0>0$, where $Z_{\alpha,n}:=\sum_{m\in\cM_n}\e^{-b_0L_nn\zeta_{n,m}^2}$ is the normalizing constant. We call the sequence $(L_n)_{n\in\bN}$ the prior exponent sequence.

(Prior mass near truth). There exists an absolute constant $\fc_1>0$  such that

\begin{align}
\Pi_{n,m}\del{\btheta\in\Theta_{n,m}:\scd_n(\sT_n(\btheta),\blambda^\star)\le \eta_{n,m} +  \zeta_{n,m}}\ge  \e^{-\fc_1n \zeta_{n,m}^2}
\end{align}

for any $\blambda^\star\in\Lambda_n^\star$,  $m\in\cM_n$ and sufficiently large $n\in\bN$.

* (Data-independence) The prior distribution $\Pi_n$ is independent to the  sample $\Y^{(n)}$.

In our specification of the prior distribution, we need to know the estimation errors $\{\zeta_{n,m}:m\in\cM_n,n\in\bN\}$ but not the approximation errors  $\{\eta_{n,m}:m\in\cM_n,n\in\bN\}$. Although <ref> involves the true natural parameter space $\Lambda_n^\star$, in most cases, we can find estimation errors which are independent of $\Lambda_n^\star$ but sufficiently tight (in a sense that they do not harm the contraction rate).  For example, see our applications (<ref>  as well as <ref>). Thus, our prior is fully adaptive to the true parameter.

We introduce the prior exponent sequence $(L_n)_{n\in\bN}$ in order to obtain model selection consistency. If  parameter estimation is of only interest, a bounded $L_n$, e.g., $L_n=1$, can be used and this does not hurt the corresponding contraction rate. However, as we will show later in <ref>, a strong regularization is sometimes needed to prevent selecting any overly complex model, which is attained via some diverging   $(L_n)_{n\in\bN}$. In that case, the contraction rate becomes slower $\sqrt{L_n}$ times.

We need a sufficiently large variational family that approximates the original posterior distribution well.

[Variational family]

There exists an absolute constant $\fc_1>0$ such that

\begin{equation}
\sup_{\blambda^\star\in\Lambda_n^\star}\P_{\blambda^\star}^{(n)}\sbr{\inf_{\Xi\in\bV_n}\kl(\Xi,\Pi_n(\cdot|\Y^{(n)}))}\le \fc_1L_nn\epsilon_n^2
\end{equation}

for any sufficiently large $n\in\bN$, where  $(L_n)_{n\in\bN}$ is the prior exponent sequence and $(\epsilon_n)_{n\in\bN}$ is the oracle rate defined in eq:oracle_rate.

In applications of our general theorems presented in the next subsections to specific problems, the most challenging task is to check <ref>. Unlike the related variational Bayes literature [Pati et al., 2018, Alquier and Ridgway, 2020, Yang et al., 2020, Zhang and Gao, 2020] that bounds the variational approximation error for the “properly complex” model only, our condition seemingly requires users to bound those errors for all the models including very large models. However, with the help of the following lemma, we can avoid this problem.

For any $\blambda^\star\in\Lambda_n$ and $\Xi\in\bP(\Theta_{n,\cM_n})$,

\begin{align}
\label{eq:vapprox_bound1}
\P_{\blambda^\star}^{(n)}\sbr{ \kl\del[1]{\Xi, \Pi_n(\cdot|\Y^{(n)})}}
&\le \Xi\sbr{\kl\del[1]{\P_{\sT_n(\btheta)}^{(n)}, \P_{\blambda^\star}^{(n)}}} +\kl(\Xi,\Pi_n)
\end{align}

for any sufficiently large $n\in\bN$.
In particular, if $\Xi\in\bP(\Theta_{n,m})$ for some $m\in\cM_n$, then

\begin{align}
\label{eq:vapprox_bound2}
\P_{\blambda^\star}^{(n)} \sbr{\kl\del[1]{\Xi, \Pi_n(\cdot|\Y^{(n)})}}
&\le \Xi\sbr{\kl\del[1]{\P_{\sT_n(\btheta)}^{(n)}, \P_{\blambda^\star}^{(n)}}} +\kl(\Xi,\Pi_{n,m})+ \log(\frac{1}{\alpha_{n,m}}).
\end{align}

Moreover, under <ref>,

\begin{equation}
\label{eq:vapprox_bound3}
\Xi\sbr{\kl\del[1]{\P_{\sT_n(\btheta)}^{(n)}, \P_{\blambda^\star}^{(n)}}}\le \fa_2+\fc_2\Xi\sbr{n\scd_n^2(\sT_n(\btheta),\blambda^\star)},
\end{equation}

where $\fa_2\in\R$ and $\fc_2>0$ are the absolute constants appearing in <ref>.

§.§ Adaptive contraction rates

Under the assumptions of the previous section, we can obtain the adaptive contraction rate of the variational posterior. Before the statement of the result, we provide contraction rates of each variational posterior $\hXi_{n,m}$ as a motivation.
The proof is almost similar to that of <ref> which appears right after, so we omit it.

Suppose that assume:exprmtassume:prior hold. Moreover, assume that there exists an absolute constant $\fc_1>0$ such that

\begin{equation*}
\P_{\blambda^\star}^{(n)}\sbr{\kl\del[1]{\hXi_{n,m},\Pi_{n}(\cdot|\Y^{(n)})}}\le \fc_1 L_n n(\eta_{n,m}+\zeta_{n,m})^2
\end{equation*}

for any $\blambda^\star\in\Lambda_n^\star$, $m\in\cM_n$ and sufficiently large $n\in\bN$. Then

\begin{equation}
\sup_{\blambda^\star\in\Lambda_n^\star}\sup_{m\in\cM_n}\P_{\blambda^\star}^{(n)}\sbr{\hXi_{n,m}\del{\scd_n(\sT_n(\btheta),\blambda^\star)\ge A_n\sqrt{L_n}(\eta_{n,m}+\zeta_{n,m})}}=\sco(1)
\end{equation}

for any diverging sequence $(A_n)_{n\in\bN}\to+\infty$.

<ref> states that the contraction rate of the individual variational posterior $\hXi_{n,m}$ is given by $\eta_{n,m}+\zeta_{n,m}$, which is the sum of the approximation and estimation errors of the model $\bM_{n,m}$. Clearly, if the best model $\bM_{n,m^*_n}$, where $m^*_n\in\argmin_{m\in\cM_n}(\eta_{n,m}+\zeta_{n,m})$, is known, one can obtain the oracle rate $\epsilon_n:=\inf_{m\in\cM_n}(\eta_{n,m}+\zeta_{n,m})$ by using the variational posterior $\hXi_{n,m^*_n}$. But this is in general not the case in practice. Our adaptive variational posterior can attain the oracle rate without any information about the best model, thus leading to the “adaptive” property.

Suppose that assume:exprmtassume:variational hold. Then

\begin{equation}
\label{eq:conv}
\sup_{\blambda^\star\in\Lambda_n^\star}\P_{\blambda^\star}^{(n)}\sbr{\hXi_n\del{\scd_n(\sT_n(\btheta),\blambda^\star)\ge A_n\sqrt{L_n}\epsilon_{n}}}=\sco(1) %M\sqrt{m^\star}\epsilon_n
\end{equation}

for any diverging sequence $(A_n)_{n\in\bN}\to\infty$.

If we choose $L_n=1$ for every $n\in\bN$, the fastest contraction rate can be achieved. This can be always done if our interest lies in parameter estimation rather than model selection. We will discuss this in more detail in the next subsection.

The original posterior distribution $\Pi_n(\cdot|\Y^{(n)})$ can be viewed as the variational posterior distribution obtained with the variational family $\bV_n=\bP(\Theta_{n,\cM_n})$. Since <ref> is trivially met in this case, all the results in this and next subsections hold for the original posterior distribution also.

§.§ Model selection consistency

In this subsection, we provide our main result for model selection. In the model selection context, it is natural to assume that the true distribution of the sample belongs to one of the models.

[True model]

For each $n\in\bN$, there exists a model index $m_n^\star\in\cM_n$ such that $\Lambda_n^\star\subset\sT_n(\Theta_{n,m_n ^\star})$. We call the model index $m_n^\star$ the true model index.

We first consider an overestimation problem of the true model index. We define

\begin{align}
\cM_{n}^{\textup{over}}(H, m_n^\star)&:=\cbr{m\in\cM_n:\zeta_{n,m}\ge H\zeta_{n,m_n^\star}}
\end{align}

for $H>1$, which can be viewed as a set of indices of overly complex models. The adaptive variational posterior distribution does not give a mass on the models in $\cM_{n}^{\textup{over}}(H, m_n^\star)$.

Suppose that <ref> as well as  assume:rateassume:true hold. Let $(H_n)_{n\in\bN}$ be a sequence such that $H_n\ge 1$ for any $n\in\bN$. Then if the prior exponent sequence $(L_n)_{n\in\bN}$ satisfies $L_n(H_n^2-1)\to+\infty$,

\begin{equation}
\label{eq:over}
\sup_{\blambda^\star\in\Lambda_n^\star}\P_{\blambda^\star}^{(n)}\sbr{\hXi_n\del{\cup_{m\in \cM_{n}^{\textup{over}}(H_n, m_n^\star)}\Theta_{n,m}}}
\end{equation}

We now consider an underestimation problem. Define

\begin{align}
\cM_{n}^{\textup{under}}(\eta^*):=\cM_{n}^{\textup{under}}(\eta^*, \Lambda_n^\star):=\cbr{m\in\cM_n:\inf_{\blambda^\star\in\Lambda_n^\star}\inf_{\btheta\in\Theta_{n,m}}\scd_n(\sT_n(\btheta),\blambda^\star)\ge\eta^*}
\end{align}

for $\eta^*>0$, which is a set of indices of less expressive models. Note that $m_n^\star\notin \cM_{n}^{\textup{under}}(\eta^*)$ for any $\eta^*>0$. The  variational posterior probability of the models in $\cM_{n}^{\textup{under}}$ vanishes asymptotically if the  expressibilty gap $\eta^*$ is large enough to detect.

Suppose that assume:exprmtassume:true hold. Then

\begin{equation}
\sup_{\blambda^\star\in\Lambda_n^\star}\P_{\blambda^\star}^{(n)}\sbr{\hXi_n\del[2]{\cup_{m\in \cM_{n}^{\textup{under}}(\uA_n\zeta_{n,m_n^\star})}\Theta_{n,m}}}
\end{equation}

for any diverging sequence $(\uA_n)_{n\in\bN}$ such that $\uA_n/\sqrt{L_n}\to+\infty$.

By <ref>, we can obtain the exact model selection consistency with suitably chosen the prior exponent sequence $(L_n)_{n\in\bN}$.

Suppose that  assume:exprmtassume:true hold. Moreover, assume that there exists a diverging sequence $(\uA_n)_{n\in\bN}\to+\infty$ such that

\begin{equation*}
\cbr{m\in\cM_n\setminus\{m_n^\star\}:\zeta_{n,m}\le \zeta_{n,m_n^\star}}\subset\cM_{n}^{\textup{under}}( \uA_n\zeta_{n,m_n^\star}).
\end{equation*}

\begin{align*}
\end{align*}

If the prior exponent sequence $(L_n)_{n\in\bN}$ satisfies $L_n((H_n^*)^2-1)\to+\infty$ and $\uA_n/\sqrt{L_n}\to+\infty$,

\begin{equation}
\inf_{\blambda^\star\in\Lambda_n^\star}\P_{\blambda^\star}^{(n)}\sbr{\hXi_n\del{\Theta_{n,m_n^\star}}}
\end{equation}

The condition of <ref> can be met when we have information about an upper bound of the true model index $m_n^\star$. The next two examples illustrates this.

[Factor model with bounded factor dimensionality]

Consider the factor model in <ref>. For a matrix $\B$, let $\sigma_{\max}(\B)$ and $\sigma_{\min}(\B)$ denote the largest and smallest singular values of $\B$. With the parameter space $\Lambda_n:=\{\bSigma\in\text{SPD}_{p_n}:\fc_0^{-1}\le \sigma_{\min}(\bSigma)\le \sigma_{\max}(\bSigma)\le \fc_0\}$ for $\fc_0>1$, and the Frobenius norm metric $\scd_n(\bSigma_1,\bSigma_2)=\fnorm{\bSigma_1-\bSigma_2}$ for $\bSigma_1,\bSigma_2\in\text{SPD}_{p_n}$, it is known that  <ref> is satisfied with $\zeta_{n,m}=\sqrt{m(p_n\log p_n)/n}$ (see e.g., Section 3.5.1 of [Han, 2021]). Suppose that $m_n^\star\le \bar{m}_n$ for some upper bound $\bar{m}_n>0$, which is known to us. Then since $H_n^*:=\inf_{m\in[p_n-1]:m>m_n^\star}(\zeta_{n,m}/\zeta_{n,m_n^\star})=\sqrt{(m_n^\star+1)/(m_n^\star)}\ge\sqrt{1+(\bar{m}_n)^{-1}}$, the choice $L_n=\uL_n\bar{m}_n$ for a sequence $(\uL_n)_{n\in\bN}$ diverging to $+\infty$ arbitrarily slowly satisfies $L_n((H_n^*)^2-1)\ge \uL_n\to+\infty$. Thus the expected probability of the event $\{m>m_n^\star\}$ goes to 0 under our variational posterior. But using such choice of the prior exponent sequence, we only obtain a (slightly) slower contraction rate $\sqrt{\uL_n\bar{m}_n}\sqrt{m_n^\star(p_n\log p_n)/n}$ of the  covariance matrix with respect to the Frobenius norm.  Ohn and Kim, 2021 observed a similar phenomenon in sparse factor models.

[Gaussian location mixture with the bounded number of components]

Consider Gaussian location mixture model in <ref> but further assume that the absolute value of every location parameter is bounded by some constant $B>0$, i.e., $\Theta_{n,m}=\Delta_m\times([-B,B]^d)^{\otimes m}$. Then by Proposition 2 of [Shen et al., 2013],  <ref> is satisfied with the Hellinger metric and the estimation error $\zeta_{n,m}=\sqrt{m\log n/n}$. Assume that the sample with size $n$ is generated from the Gaussian location mixture with the number of components $m_n^\star$ with $m_n^\star\le \log n$. Assuming the $\log n$ upper bound of the true number of components has been made also in some related literature [Wu and Yang, 2020, Ohn and Lin, 2020]. Then $H_n^*:=\inf_{m\in\bN:m>m_n^\star}(\zeta_{n,m}/\zeta_{n,m_n^\star})=\sqrt{(m_n^\star+1)/(m_n^\star)}\ge\sqrt{1+(\log n)^{-1}}$. Thus, the choice $L_n=(\log n)^{1+\iota}$ for arbitrary small $\iota>0$ satisfies $L_n((H_n^*)^2-1)\ge (\log n)^{\iota}\to+\infty$, and thus the expected probability of the event $\{m>m_n^\star\}$ goes to 0 under our variational posterior. But with such a choice of the prior exponent sequence, we only obtain a slower contraction rate $\sqrt{m_n^\star \log n/n}(\log n)^{(1+\iota)/2}$ for estimation of the true density function with respect to the Hellinger  metric.

§ APPLICATION: ADAPTIVE VARIATIONAL DEEP LEARNING

In this section, we propose a novel Bayesian deep learning procedure, called adaptive variational deep learning, based on the proposed adaptive variational Bayes. Applying our general result in <ref>, we show that the adaptive variationl deep learning adaptively attains minimax optimality for a number of  statistical problems including regression function estimation, conditional class probability function estimation for binary classification and intensity estimation for Poisson point processes.

Deep learning is a statistical inference method in which a target function we wish to estimate, e.g., a regression function, is modeled by (deep) neural networks.  In this section, we focus on the case the target function is a real-valued function supported on the $d$-dimensional unit cube $[0,1]^d$. But our analysis can be easily generalized to any compact input domain and multivariate output. We introduce some notations used throughout this section. Let $\cF^d$ be the set of all measurable real-valued functions supported on $[0,1]^d$. For $f\in\cF^d$ and $p\in\bN$, let $\|f\|_p:=(\int_{[0,1]^d}|f(\x)|^p\d\x)^{1/p}$ denote the usual $\scL^p$ norm and $\|f\|_\infty:=\sup_{\x\in[0,1]^d}|f(\x)|$ the $\scL^\infty$ norm.

§.§ Adaptive variational posterior over neural networks

§.§.§ Mathematical definition of neural networks

We introduce neural network models for function estimation. We first  introduce a notion of network parameters and  related notations. For $M_1,M_2\in\bN$, let $\hTheta_{M_1,M_2}:=\R^{M_2\times M_1}\times \R^{M_2}$. For $(\W,\b)\in\hTheta_{M_1,M_2}$, the matrix $\W$ and the vector $\b$  are called a \textit{weight matrix} and  a \textit{bias vector}, respectively, of a (neural network) layer which takes a $M_1$-dimensional input (form the previous layer) to produce a $M_2$-dimensional output. Moreover, for a positive integer $K∈_≥2$ larger than 2 and a $(K+1)$-dimensional vector of positive integers $_1:(K+1):=(M_1,…, M_K+1)∈^K+1$, we denote $__1:(K+1):=⊗_k=1^K_M_k,M_k+1$. We can understand $∈__1:(K+1)$ as a network parameter of a neural network with the number of layers $K$ and the number of nodes in the $k$-th layer $M_k+1$ for $k∈[K]$. So the space $Θ:=∪_K=2^∞∪__1:(K+1)∈^K+1_$ can be viewed as the set of all network parameters. Since we focus on the estimation of real-valued functions on supported on $[0,1]^d$, we only consider network parameters in the subset $Θ^d:=∪_K=2^∞∪__1:(K+1)∈^K+1:M_1=d, M_K+1=1_⊂Θ$.

We next define a  \textit{neural network function} or simply \textit{neural network} induced by a network parameter. For a network parameter $=((_k, _̱k))_k∈[K]$, let $()$ denote a  neural network function such that
\begin{align}
\net(\btheta):\x\mapsto[\W_{K},\b_{K}]\circ\text{ReLU}\circ[\W_{K-1},\b_{K-1}]\circ\cdots\circ\text{ReLU}\circ[\W_{1},\b_{1}]\x,
\end{align}
\begin{align*}
[\W_{k},\b_{k}]:\x'\mapsto \W_k\x'+\b_k
\mbox{ and }\text{ReLU}:\x'\mapsto \x'\vee\zero,
\end{align*}
that is, $[_k,_̱k]$ denotes the affine transformation represented as a multiplication by the matrix $_k$ and an addition of the vector $_̱k$ and $ReLU$ does the elementwise ReLU (rectified linear unit) activation function.

\subsubsection{Parameter spaces, prior and variational posterior}

In this application of adaptive variational Bayes to deep learning, we consider a multiple disjoint subsets of the parameter space $Θ^d$ which are indexed by \textit{depth} (i.e., the number of hidden layers) and \textit{width} (i.e., the (maximum) number of hidden nodes of each layer). For given (network) depth $K∈_≥2$ and (network) width $M∈$, we define
\begin{equation}
\Theta_{(K,M)}:=\Theta^d_{(K,M)}:=\hTheta_{(d,M_2,\dots,M_{K},1)} \mbox{ with } M_2=\cdots=M_{K}=M.
\end{equation}
We let  $p_(K,M)$ be the dimension of a space $Θ_(K,M)$. Note that
\begin{align*}
p_{(K,M)}:=(d+1)M+(K-2)(M^2+M)+(M+1)\le (d+1)KM^2.
%&\le M\{d+1+(K-1)(M+1) \} \\
%& \le KM \max\{M+1, d+1 \}\\
%& \le KM (Md+1).
\end{align*}
Since $Θ_(K,M)≅^p_(K,M)$, we can deal with each member in the parameter space $Θ_(K,M)$ as a $p_(K,M)$-dimensional real vector. To ensure  prior concentration probabilities are large enough, we restrict network parameters to be bounded. For a given \textit{magnitude bound} $B>0$, define
\begin{align}
\Theta^{\le B}_{(K,M)}:=\cbr{\btheta\in\Theta_{(K,M)}:|\btheta|_\infty\le B}.
\end{align}
In our collection of neural network models, we let each parameter space be of the form $Θ^≤B_n_(K,M)$ for each $n∈$, where the sequence of magnitude bounds ${B_n}_n∈⊂_≥0$ is a positive sequence that diverges at a suitable rate.

We describe the prior and variational posterior. Let $_n⊂_≥2×$ be a set of some model indices (i.e., pairs of depth and width) to be considered in consequent inference tasks. We impose the prior distribution
\begin{equations}
\label{eq:dnn_prior}
\Pi_n=&\sum_{(K,M)\in\cM_n}\alpha_{n,(K,M)}\Pi_{n,(K,M)},\\
&\mbox{ where }\Pi_{n,(K,M)}:=\bigotimes_{j=1}^{p_{(K,M)}}\Unif(-B_n, B_n)\\
&\mbox{ and }\alpha_{n,(K,M)}:= \frac{\e^{-b_0(KM)^2\log n}}{\sum_{(K',M')\in\cM_n}\e^{-b_0(K'M')^2\log n}}
\end{equations}
for $b_0>0$. For each model index $(K,M)∈_n$, we consider the variational family
\begin{equation}
\label{eq:dnn_variational_family}
\bV_{n,(K,M)}=\cbr{\bigotimes_{j=1}^{p_{(K,M)}}\Unif(-\psi_{1,j}, \psi_{2,j}):-B_n\le \psi_{1,j}<\psi_{2,j}\le B_n}.
\end{equation}
Since the parameter spaces $[0]Θ^≤B_n_(K,M)_(K,M)∈_n$ are disjoint, due to \cref{thm:compute}, the variational posterior is given by
\begin{align}
\hXi_n=&\sum_{(K,M)\in\cM_n}\hgamma_{n,(K,M)}\hXi_{n,(K,M)},  \label{eq:dnn_estimate}\\
& \mbox{ where } \hXi_{n,(K,M)}\in\argmin_{\Xi\in\bV_{n,(K,M)}} \scE\del[1]{\Xi, \Pi_{n,(K,M)},\sq_n,\Y^{(n)}} \label{eq:dnn_indiv}\\
&\mbox{ and }  \hgamma_{n,(K,M)}=\frac{\alpha_{n,(K,M)} \e^{-\scE\del[1]{\hXi_{n,(K,M)}, \Pi_{n,(K,M)},\sq_n,\Y^{(n)}}}}
{\sum_{(K',M')\in\cM_n}\alpha_{n,(K',M')} \e^{-\scE\del[1]{\hXi_{n,(K',M')}, \Pi_{n,(K',M')},\sq_n,\Y^{(n)}}}}.
\end{align}

Although solving the minimization problem in \labelcref{eq:dnn_indiv} is in general not a  hard task for deep learning practitioners nowadays, we provide a detailed computation algorithm for the sake of completeness in \cref{appen:deep_compute}.

\subsection{Nonparametric regression with Gaussian errors using neural networks}
\label{sec:dnn:reg}

In this subsection, we consider a nonparametric regression experiment with standard Gaussian errors. In this case, we observe $n$ independent outputs $^(n):=(Y_1,…, Y_n)∈_n:=^n$, where each $Y_i$ is associated with a fixed $d$-dimensional input $_i$ that has been rescaled so that $_i∈[0,1]^d$.  We model the sample $^(n)$ using the probability model $(^n;_n^Ga,^d)$  with the likelihood function $_n^Ga:^d×^n↦_≥0$ defined as
\begin{equation}
\sq_n^{\text{Ga}}(f, \Y^{(n)})=\sq_n^{\text{Ga}}(f, \Y^{(n)};(\x_i)_{i\in[n]})
\end{equation}
for $f∈^d,$ where $_(0,1)$ stands for the density function of the standard normal distribution. Here, we assume that the variance of the output distribution is fixed to be 1 for technical simplicity, however, we can easily extend this to an unknown variance case. We denote by $_f^(n)$ the distribution that has the density function $_n(f, ·)$.  Assuming that $^(n)∼_f^⋆^(n)$ for some \textit{true regression function} $f^⋆∈^d$, our aim is to accurately estimate $f^⋆$ in terms of the empirical $^2$ distance which is defined as
\begin{equation*}
\nnorm{f_0-f_1}=\del[3]{\frac{1}{n}\sum_{i=1}^n(f_0(\x_i)-f_1(\x_i))^2}^{1/2}
\end{equation*}
for $f_0,f_1∈^d$. The next lemma shows that  \cref{assume:exprmt} is satisfied for this regression experiment.

\begin{lemma}
\label{lem:gaussreg}
A sequence of the Gaussian regression experiments $\{(\R^n,\bP(\R^n;\sq_n^{\textup{Ga}},\cF^d))\}_{n\in\bN}$ satisfies \cref{assume:exprmt} with the metric $\scd_n:\cF^d\times \cF^d\mapsto\R_{\ge0}$ defined as $\scd_n(f_0,f_1)=\nnorm{f_0-f_1}$ for $f_0,f_1\in\cF^{d}$.
\end{lemma}

\subsubsection{Oracle contraction rate for regression function estimation}

We investigate optimality of the proposed adaptive variational deep learning for regression. We first provide the oracle type contraction property, in which the contraction rate is given by the minimum of the sum of the approximation and estimation errors over neural network architectures. With this oracle result, we can obtain adaptive optimal results for a number of regression problems without difficulty.

\begin{theorem}[Oracle contraction rate, regression function]
\label{thm:dnn:oracle}
Consider a sequence of the Gaussian regression experiments $\{(\R^n,\bP(\R^n;\sq_n^{\textup{Ga}},\cF^d))\}_{n\in\bN}$.
Let $\cF^\star\subset\cF^d$. For each $n\in\bN$, let $\cM_n\subset\bN_{\ge2}\times \bN$ be a set of some depth-width pairs such that $|\cM_n|\lesssim n^\iota$ for some absolute constant $\iota>0$ and $\max_{(K,M)\in\cM_n}(K\vee M)\lesssim n$. Assume that $1\le B_n\lesssim n^\nu$ for some $\nu>0.$ Then the variational posterior  \labelcref{eq:dnn_estimate} satisfies
\begin{equation}
\sup_{f^\star\in\cF^\star}\P_{f^\star}^{(n)}\sbr{\hXi_n\del{\nnorm{\net(\btheta)- f^\star}\ge A_n\epsilon_n(\cF^\star)}}=\sco(1)
\end{equation}
for any diverging sequence $(A_n)_{n\in\bN}\to+\infty$ , where
\begin{equation}
\label{eq:dnn_oracle}
\epsilon_n(\cF^\star):=\inf_{(K,M)\in\cM_n}\cbr{\sup_{f^*\in\cF^\star}\inf_{\btheta\in\Theta^{\le B_n}_{(K,M)}}\norm{\net(\btheta)-f^*}_\infty+KM\sqrt{\frac{\log n}{n}}}.
\end{equation}
\end{theorem}

Our result is the first adaptive optimal one for the fully connected (i.e., non-sparse) neural network models. Many of the related studies showing minimax optimality of neural network estimators assumed sparsity of the neural networks (in both adaptive and nonadaptive fashions) to control estimation variance while maintain high approximation ability, for example, Schmidt-Hieber, 2020, Suzuki, 2019, Imaizumi and Fukumizu, 2020, Ohn and Kim, 2020 and Kim et~al., 2021. However, to find an optimal sparse neural network is computationally expensive since we need to explore all possible zero-one configurations of the network parameter. Notable exceptions are Kohler and Langer, 2021 and Jiao et~al., 2021. The authors of these papers showed that optimality of non-sparse neural network estimators for regression but their results are nonadaptive since the network width and/or depth should be chosen based on the smoothness of the true regression function.

On the Bayesian side, Polson and Ro{\v{c}}kov{\'a}, 2018 proposed to use a spike-and-slab prior consisting of the Dirac-delta distribution at zero and an uniform distribution for the network parameter. They proved that the resulting posterior distribution over sparse neural networks gain the adaptive optimality for estimation of a H\"older smooth regression function. However, computing the posterior with the spike-and-slab prior is also computationally  expensive due to the exactly same reason as in finding the optimal sparse neural network. Ch{\'e}rief-Abdellatif, 2020 showed that under the similar spike-and-slab prior, the variational approximation of the posterior distribution can also have the optimal contraction rate. Moreover, by using the similar approach of the model selection variational Bayes of [Zhang and Gao, 2020], they achieved the adaptivity. However, even for the variational posterior, its computation under the spike-and-slab prior is still almost intractable. The proposed Bayesian deep learning approach overcomes such a computational drawback of  the spike-and-slab prior while not losing the adaptive optimality.

We now provide some explicit examples that the adaptive variational deep learning regression is optimal in the following two subsections.

\subsubsection{Estimation of H\"older smooth regression functions}

In this example, we assume the true regression function is H\"older $β$-smooth. This assumption is the most standard and classical assumption in the related literature.  Let $^p,d$ denote the class of $p$-times differentiable functions supported on $[0,1]^d$ for $p∈$ and $^0,d$ denote the class of continuous functions on $[0,1]^d$. The \textit{H\"older space} of smoothness $β>0$ with domain $[0,1]^d$ and radius $F_0>0$ is defined as
\begin{equation*}
\cH^{\beta, d, F_0}:=\cbr{f\in\cC^{\ceil{\beta-1},d}:\|f\|_{\cH^{\beta,d}}\le F_0},
\end{equation*}
where $·_^β,d$ denotes the \textit{H\"older norm} defined by
\begin{align*}
\|f\|_{\cH^{\beta,d}}
,\max_{\a\in\bN_0^d:|\a|_1= \ceil{\beta-1}}\sup_{\x_1,\x_2\in [0,1]^d, \x_1\neq\x_2 }\frac{|\partial^{\a}f(\x_1)-\partial^{\a}f(\x_2)|}{|\x_1-\x_2|_\infty^{\beta- \ceil{\beta-1}}}}.
\end{align*}
Here, $∂^f$ denotes the partial derivative of $f$ of order $:=(a_j)_j∈[d]$, that is, $∂^f:=∂^||_1f/∂x_1^a_1⋯∂x_d^a_d$.

The following corollary shows that adaptive varational deep learning attains the minimax optimal contraction rate $n^-β/(2β+d)$  [Stone, 1982], up to a logarithmic factor.

\begin{corollary}[H\"older smooth regression function]
\label{col:dnn:holder}
Let $\beta>0$, $d\in\bN$ and $F_0>0$. For each $n\in\bN$, consider the set of depth-width pairs given as
\begin{align}
\label{eq:architecture_set}
\cM_n=\cbr{\del{\floor{k_1\log n}\vee 2,\floor{n^{k_2/\log n}}}:k_1\in\sbr{\ceil{\log\log n}\vee1},k_2\in\sbr{0: \ceil{\log n}}}
\end{align}
and set $B_n= n^{\nu_0}$ with $\nu_0\ge1/2$. Then the variational posterior \labelcref{eq:dnn_estimate} satisfies
\begin{equation}
\sup_{f^\star\in\cH^{\beta, d, F_0}}\P_{f^\star}^{(n)}\sbr{\hXi_n\del{\nnorm{\net(\btheta)- f^\star}\ge n^{-\frac{\beta}{2\beta+d}}\log^2n}}=\sco(1).
\end{equation}
\end{corollary}

\subsubsection{Estimation of composition structured regression functions}
\label{sec:dnn:reg:dag}

To avoid the curse of dimensionality that exists in the convergence rate of estimation of H\"older smooth  regression function, one can make a certain structural assumption on the regression function. Schmidt-Hieber, 2020 and Kohler and Langer, 2021 considered so-called composition structured regression functions which include a single-index model [Ga{i}ffas and Lecu{\'e}, 2007], an additive model [Stone, 1985, Buja et~al., 1989] and a generalized additive model with an unknown link function [Horowitz and Mammen, 2007] as particular instances. Under this structure assumption, they derived the improved convergence rates of sparse and non-sparse neural networks respectively, but note again, their estimators are nonadaptive.

We here consider slightly more general notion of composition structure for functions called directed acyclic graph (DAG) structure and derives a contraction rate of our adaptive variational deep learning. To ease the description, we define ``select'' and ``append'' operators as follows. For an input dimension $s∈$ and a set of ``active'' elements  $⊂[s]$, we define a function $Sel_:^s↦^||$ as
\begin{equation*}
\text{Sel}_{\cJ}(\x)=(x_j)_{j\in\cJ}.
\end{equation*}
Moreover,  for a function  $f∈^||$, we define a function $App_f,:^s↦^s+1$ as
\begin{equation*}
\text{App}_{f,\cJ}(\x)=\del{\x,f(\text{Sel}_{\cJ}(\x))}.
\end{equation*}
For given  $(_ℓ)_ℓ∈[r]$ and $(f_ℓ)_ℓ∈[r]$ with $_ℓ⊂[d+ℓ-1]$ and $f_ℓ∈^|_ℓ|$,  we denote by $ [0]·; (f_ℓ)_ℓ∈[r],(_ℓ)_ℓ∈[r]$ the function such that
\begin{equation*}
\sF\del[1]{\x; (f_\ell)_{\ell\in[r]},(\cJ_\ell)_{\ell\in[r]}}=f_r\del[2]{\text{Sel}_{\cJ_r}\del[2]{\text{App}_{f_{r-1},\cJ_{r-1}}\del[1]{\cdots\text{App}_{f_{2},\cJ_{2}}\del[1]{\text{App}_{f_{1},\cJ_{1}}(\x)}\cdots}}}.
\end{equation*}
We can say that the function defined above has a DAG structure for the following reasons. If we define the directed graph that consists of the set of vertices $=[d+r]$ and the set of edges $$ containing a direct edge $j→ℓ$ if and only if $j∈_ℓ$, then it can be easily checked that the  graph $(, )$ is a DAG. The DAG structure for regression functions is more general than those of  [Schmidt-Hieber, 2020, Kohler and Langer, 2021] since this allows direct edges from the input to a function of any stages, for example, the function $f_r$ of the last stage may depend ``directly'' on a part of the input vector. \cref{fig:dag_example} presents an example of  a DAG structured function.

We define the class of DAG structured functions with the number of components $r∈$, smoothness $(β_ℓ)_ℓ∈[r]∈_+^r$ and \textit{intrinsic dimensions} $(s_ℓ)_ℓ∈[r]∈^r$ as
\begin{align*}
&\cF^{\textup{DAG}}\del{r,(\beta_\ell)_{\ell\in[r]},(s_\ell)_{\ell\in[r]}, F_0 }\\
&:=\cbr{\sF\del[1]{\cdot; (f_\ell)_{\ell\in[r]},(\cJ_\ell)_{\ell\in[r]}}:\cJ_\ell\subset[d+k-1], |\cJ_\ell|\le s_\ell, f_\ell\in\cH^{\beta_\ell,|\cJ_\ell|,F_0},0\le f_1,\dots,f_{r-1}\le 1}.
\end{align*}
We here assume that the outputs of the  functions $f_1,…, f_r-1$ of the intermediate stages lie in the interval $[0,1]$ since they can be inputs of the functions of the following stages. This can be easily generalized to an arbitrary compact subset of the real line, but we do not consider this generalization for technical simplicity.

The next corollary presents the contraction rate of the adaptive variational posterior for the estimation of the DAG structured regression functions.

\begin{corollary}[DAG structured regression function]
\label{col:dnn:composite}
Let $r\in\bN$, $\bbeta:=(\beta_\ell)_{\ell\in[r]}\in\R_+^r$,  $\s:=(s_\ell)_{\ell\in[r]}\in\bN^r$ and $F_0>0$. Then with the set of depth-width pair \labelcref{eq:architecture_set} and the magnitude bound $B_n= n^{\nu_0}$ with $\nu_0\ge1/2$, the variational posterior \labelcref{eq:dnn_estimate} satisfies
\begin{equation}
\sup_{f^\star\in\cF^{\textup{DAG}}\del{r,\bbeta,\s, F_0}} \P_{f^\star}^{(n)}\sbr{\hXi_n\del{\nnorm{\net(\btheta)- f^\star}\ge \max_{\ell\in[r]}n^{-\frac{\beta_{\ge \ell}}{2\beta_{\ge \ell}+s_\ell}}\log^2n}}=\sco(1),
\end{equation}
where we let $\beta_{\ge \ell}:=\beta_\ell\prod_{h=\ell+1}^r(\beta_h\wedge 1)$ for $\ell\in[r-1]$ and $\beta_{\ge r}:=\beta_r$.
\end{corollary}

If the intrinsic dimensions $s_1, …, s_r$ are much smaller than the input dimension $d$, the convergence rate improves substantially compared to the rate for the H\"older smooth functions.

\begin{figure}
\begin{center}
\begin{tikzpicture}[circ/.style={draw,circle,minimum size=28pt, thick},scale=0.8]
\tikzset{my/.style={->,line width=1pt}}
\node[circ] at (1,0) (x1) {$x_1$};
\node[circ] at (3,0) (x2) {$x_2$};
\node[circ] at (5,0) (x3) {$x_3$};
\node[circ] at (7,0) (x4) {$x_4$};
\node[circ] at (9,0) (x5) {$x_5$};
\node[circ] at (11,0) (x6) {$x_6$};
\node[circ] at (5.3,1.7) (f1) {$f_1$};
\node[circ] at (2,2.8) (f2) {$f_2$};
\node[circ] at (9,2) (f3) {$f_3$};
\node[circ] at (6,4) (f4) {$f_4$};
\draw[my] (x1) edge (f2);
\draw[my] (x2) edge (f2);
\draw[my] (x3) edge (f1);
\draw[my] (x4) edge (f1);
\draw[my] (x4) edge (f4);
\draw[my] (x5) edge (f3);
\draw[my] (x6) edge (f3);
\draw[my] (f1) edge (f2);
\draw[my] (f1) edge (f3);
\draw[my] (f2) edge (f4);
\draw[my] (f3) edge (f4);
\end{tikzpicture}
\caption{An example of a DAG structured function $\sF\del[0]{\cdot; (f_\ell)_{\ell\in[4]},(\cJ_\ell)_{\ell\in[4]}}$ with $\cJ_1=\{3,4\}$, $\cJ_2=\{1,2,7\}$, $\cJ_3=\{5,6,7\}$ and $\cJ_4=\{4,8,9\}$. In this case, the dimension of the input $\x$ is $d=6$ but the intrinsic dimensions of the function are $s_1=2$ and $s_2=s_3=s_4=3$, all of which are less than 6.}
\label{fig:dag_example}
\end{center}
\end{figure}

\subsection{Binary classification with neural networks}

In this subsection, we consider a statistical experiment of binary classification. Let $:=(,Y)∈:=[0,1]^d×{0,1}$ be a pair of input-label variables such that $∼([0,1]^d)$ and $Y|=∼(ω())$ for some \textit{conditional class probability function} $ω∈^d_[ϰ,1-ϰ]:=[0]ω∈^d:ϰ≤ω≤1-ϰ$ for $ϰ∈(0,1/2)$. Here, we assume that the conditional class probability is bounded away from 0 and 1. This is in order to to avoid diverging behavior of the likelihood near 0 or 1. Such an assumption has been commonly made in analyzing logistic models \cite[e.g.,][]{van2014asymptotically,ohn2020nonconvex}. For each $n∈$, we observe $n$ iid input-label pairs $^(n):=((_i,Y_i))_i∈[n]$ which are independently generated from the above probability model, i.e., $^(n)∼_ω^(n)$ for $ω∈^d_[ϰ,1-ϰ]$, where $_ω^(n)$ is a distribution with the likelihood $_n^(ω,·)$ defined as
\begin{equation}
\sq_n^{\textup{Ber}}(\omega,\Z^{(n)})=\prod_{i=1}^n\sbr{\omega(\X_i)^{Y_i}\del{1-\omega(\X_i)}^{1-Y_i}}.
\end{equation}

\begin{lemma}
\label{lem:logistic}
Let $\varkappa\in(0,1/2)$ be an arbitrary small positive number. A sequence of the binary classification experiments  $\del[1]{\del[0]{\fZ^{\otimes n},\bP(\fZ^{\otimes n};\sq_n^{\textup{Ber}},\cF^d_{[\varkappa,1-\varkappa]})}}_{n\in\bN}$ satisfies \cref{assume:exprmt:tail} with the metric $\scd_n:\cF^d_{[\varkappa,1-\varkappa]}\times\cF^d_{[\varkappa,1-\varkappa]}\mapsto \R_{\ge0}$ defined as $\scd_n(\omega_0,\omega_1)=\|\omega_0-\omega_1\|_2$ for $\omega_0,\omega_1\in\cF^d_{[\varkappa,1-\varkappa]}$.
\end{lemma}

We  model $ω$ itself by a neural network, i.e., the likelihood of a network parameter $$ is given by $∏_i=1^n(()(_i))^Y_i(1-()(_i))^1-Y_i$. To ensure that the neural network estimate is bounded away from 0 and 1, we apply a truncation operator to our neural network estimate as
\begin{equation*}
\net_{[\varkappa,1-\varkappa]}(\btheta):=(\net(\btheta)\vee\varkappa)\wedge(1-\varkappa)
\end{equation*}
for a network parameter $$.

The following corollary shows that if the true conditional probability function is H\"older $β$-smooth, our adaptive variational deep learning is able to estimate it with the optimal rate. The proof is almost similar to the proof of \cref{col:dnn:holder}, so we omit it.

\begin{corollary}[H\"older  smooth conditional class probability]
\label{col:dnn:logistic}
Consider a sequence of the binary classification experiments $\del[1]{\del[0]{\fZ^{\otimes n},\bP(\fZ^{\otimes n};\sq_n^{\textup{Ber}},\cF^d_{[\varkappa,1-\varkappa]})}}_{n\in\bN}$  for $d\in\bN$ and an arbitrary small $\varkappa\in(0,1/2)$. Let $\beta>0$ and $F_0>0$.
Then with the set of depth-width pair \labelcref{eq:architecture_set} and the magnitude bound $B_n= n^{\nu_0}$ with $\nu_0\ge1/2$, the variational posterior \labelcref{eq:dnn_estimate} satisfies
\begin{equation}
\sup_{\omega^\star\in\cH^{\beta, d, F_0}_{[\varkappa,1-\varkappa]}} \P_{\omega^\star}^{(n)}\sbr{\hXi_n\del{\Norm{\net_{[\varkappa,1-\varkappa]}(\btheta)- \omega^\star}_2\ge n^{-\frac{\beta}{2\beta+d}}\log^2n}}=\sco(1),
\end{equation}
where  $\cH^{\beta, d, F_0}_{[\varkappa,1-\varkappa]}:=\cbr[1]{\omega\in\cH^{\beta, d, F_0}:\varkappa\le \omega\le 1-\varkappa}.$
\end{corollary}

\subsection{Intensity estimation for Poisson point processes using neural networks}
\label{sec:dnn:ppp}

Let  $$ be the Borel $σ$-field of $[0,1]^d$. For each $n∈$, we consider a Poisson point process $Y$ on $[0,1]^d$ with intensity $λ:[0,1]^d↦_≥0$, that is, $Y$ is an integer-valued random measure such that
\begin{enumerate}
\item  for any $k\in\bN$ and any disjoint $X_1,\dots, X_k\subset\cX $, the random variables $Y(X_1),\dots ,Y(X_k)$ are independent;
\item for any $X\in\cX$, $Y(X)\sim \Pois\del{\int_X\lambda(\x)\d\x}$.
\end{enumerate}
Let $$ be a set of realizations of the Poisson point process $Y$, i.e.,
\begin{equation*}
\bar\fY:=\cbr{\sum_{j=1}^{N}\delta(\cdot;\x_{j}'):N\in\bN,\x_{j}'\in[0,1]^d}.
\end{equation*}
The aim is to recover the intensity function $λ$ of the Poisson point process $Y$ based on $n$ independent realizations $^(n):=(Y_i)_i∈[n]∈^⊗n$ of $Y$. With a reference measure being the Poisson point process with constant intensity $λ≡1$, we can consider the likelihood function $_n^:_≥0^d×^⊗n↦_≥0$ given by
\begin{equation}
\sq_n^{\textup{PPP}}(\lambda, \Y^{(n)})=\prod_{j=1}^n\exp\del{\int\log \lambda(\x)\d Y_i(\x)-\int_{[0,1]^d}(\lambda(\x)-1)\d\x},
\end{equation}
where $_≥0^d$ denotes the set of measurable nonnegative functions supported on $[0,1]^d$. Let $_λ^(n)$ be a distribution of the Poisson point process $Y$ with the likelihood function $_n(λ,·)$.

We first check \cref{assume:exprmt:tail}. For this, we assume that every intensity function is bounded away from zero and infinity for technical simplicity. That is, we assume that every intensity belongs to the class
\begin{equation*}
\cF^d_{[\varkappa_{\min},\varkappa_{\max}]}:=\cbr[0]{\lambda\in\cF^d:\varkappa_{\min}\le \lambda\le \varkappa_{\max}}
\end{equation*}
for given $ϰ_max> ϰ_min>0$.

\begin{lemma}
\label{lem:ppp}
For any $\varkappa_{\max}> \varkappa_{\min}>0$, a sequence of the Poisson point process experiments $\del[1]{\del[0]{\bar\fY^{\otimes n},\bP(\bar\fY^{\otimes n};\sq_n^{\textup{PPP}},\cF^d_{[\varkappa_{\min},\varkappa_{\max}]})}}_{n\in\bN}$ satisfies \cref{assume:exprmt:tail}  with the metric $\scd_n:\cF^d_{[\varkappa_{\min},\varkappa_{\max}]}\times\cF^d_{[\varkappa_{\min},\varkappa_{\max}]}\mapsto\R_{\ge0}$ defined as $\scd_n(\lambda_0,\lambda_1)=\|\lambda_0-\lambda_1\|_2$ for $\lambda_0,\lambda_1\in\cF^d_{[\varkappa_{\min},\varkappa_{\max}]}$.
\end{lemma}

We directly model $λ$ by a neural network in a sense that the likelihood of a network parameter $$ is given by $_n^((),^(n))$. To ensure boundeness of  neural network intensity estimates, we apply a truncation operator as
\begin{equation*}
\net_{[\varkappa_{\min},\varkappa_{\max}]}(\btheta):=(\net(\btheta)\vee\varkappa_{\min})\wedge\varkappa_{\max}
\end{equation*}
for a network parameter $$.

In the next corollary, we show that the adaptive variational deep learning attains the optimal rate for estimation of H\"older $β$-smooth intensity functions [Kutoyants, 2012]. The proof is almost similar to the proof of \cref{col:dnn:holder}, thus we omit it.

\begin{corollary}[H\"older smooth intensity]
\label{col:dnn:ppp}
Consider a sequence of the Poisson point process experiments $\del[1]{\del[0]{\bar\fY^{\otimes n},\bP(\bar\fY^{\otimes n};\sq_n^{\textup{PPP}},\cF^d_{[\varkappa_{\min},\varkappa_{\max}]})}}_{n\in\bN}$ for $d\in\bN$ and $\varkappa_{\max}>\varkappa_{\min}>0$. Let $\beta>0$ and $F_0>0$.
Then with the set of depth-width pair \labelcref{eq:architecture_set} and the magnitude bound $B_n= n^{\nu_0}$ with $\nu_0\ge1/2$, the variational posterior \labelcref{eq:dnn_estimate} satisfies
\begin{equation}
\sup_{\lambda^\star\in\cH^{\beta, d, F_0}_{[\varkappa_{\min},\varkappa_{\max}]}} \P_{\lambda^\star}^{(n)}\sbr{\hXi_n\del{\Norm{\net_{[\varkappa_{\min},\varkappa_{\max}]}(\btheta)- \lambda^\star}_2\ge n^{-\frac{\beta}{2\beta+d}}\log^2n}}=\sco(1),
\end{equation}
where $\cH^{\beta, d, F_0}_{[\varkappa_{\min},\varkappa_{\max}]}:=\cbr[1]{\lambda\in\cH^{\beta, d, F_0}:\varkappa_{\min}\le \lambda\le \varkappa_{\max}}.$
\end{corollary}

\section{Comparison with model selection variational Bayes}
\label{sec:compare}

In this section, we compare the adaptive variational posterior to the variational posterior accompanied with a model selection procedure, which is considered in   [Zhang and Gao, 2020]. Note that the adaptive variational posterior that we propose is the variational distribution supported on \textit{the entire parameter space} $Θ_n,_n$.  Whereas, Zhang and Gao, 2020 focus on the variational posterior distribution supported on  \textit{the single best parameter space}, which is obtained via the minimization program
\begin{equation}
\min_{m\in\cM_n}\min_{\Xi\in\bV_{n,m}}
\sbr{\scE(\Xi,\Pi_{n,m},\sq_n,\Y^{(n)})-\log\alpha_{n,m}}.
\end{equation}
Note that the distribution $_n,m̂_n$ with $m̂_n∈_m∈_n_n,m$ is a solution to the above minimization problem since
\begin{align*}
-\log\hgamma_{n,m}\propto -\log &\alpha_{n,m}+ \scE(\hXi_{n,m},\Pi_{n,m},\sq_n,\Y^{(n)})\\
&=\min_{\Xi_{n,m}\in\bV_{n,m}} \sbr{\scE(\Xi_{n,m},\Pi_{n,m},\sq_n,\Y^{(n)})-\log\alpha_{n,m}}.
\end{align*}
We call $_n,m̂_n$ the \textit{model selection variational posterior}.   Zhang and Gao, 2020 derived the contraction rate for the model selection variational posterior.

\begin{remark}
Note that every individual variational posterior needs to be obtained to optimize the selection criterion for the model selection variational Bayes. Therefore, computational costs of the two variational Bayes approaches are  exactly the same.
\end{remark}

The adaptive variational Bayes can be more useful in practical situations since it fully approximates the original posterior distributions over the entire model thus provides  fruitful information. A notable example is that it provides the (variational) posterior model probabilities and thus quantifies the uncertainty around the model structure. One can conduct model comparison using such posterior model probabilities. Also, this can avoid ``overconfident'' predictions and inferences, which arises due to that the model uncertainty is ignored [Draper, 1995]. To add to these practical advantages, it turns out   that the adaptive variational Bayes is a theoretically superior to the model selection variational Bayes as we show in the in the next subsection.

\subsection{Theoretical comparison:  the risk bound perspective}

We here argue that the adaptive variational posterior is theoretically superior over the model selection variational posterior in view of a nonasymptotic \textit{risk bound}. We first show that the model selection variational posterior $_n,m̂_n$ is a suboptimal approximation of the original posterior compared with the adaptive varational posterior $_n$, which is trivially obtained.

\begin{proposition}
\label{prop:inferior_msvb}
Let $\Y^{(n)}$ be a sample generated from the distribution $\P_\star^{(n)}$. Then
\begin{equation}
\kl\del{\hXi_n, \Pi_n(\cdot|\Y^{(n)})}
\le \kl\del{\hXi_{n,\hat{m}_n}, \Pi_n(\cdot|\Y^{(n)})}
\end{equation}
with $\P_\star^{(n)}$-probability 1, where $\hat{m}_n\in\argmax_{m\in\cM_n}\hgamma_{n,m}$.
\end{proposition}

\begin{proof}
The result trivially follows from the fact that $\cup_{m\in\cM_n}\bV_{n,m}\subset\bV_n$, which in particular implies that $\hXi_{n,\hat{m}_n}\in\bV_n.$
\end{proof}

We now reveal that a suboptimal variational approximation error to the original posterior leads to a suboptimal nonasymptotic risk bound. Indeed, a risk of a variational distribution can be upper bounded by the sum of the variational approximation error measured by the KL divergence and the term independent to the variational distribution as shown in the next proposition.

\begin{proposition}[Nonasymptotic risk bound]
\label{prop:compare:risk}
Consider the statistical setup in \cref{sec:method:exprmt}. Let  $\Pi_n\in\bP(\Theta_{n,\cM_n})$ be a prior distribution and $\scd_n:\Lambda_n\times\Lambda_n\mapsto\R_{\ge0}$ be a metric on the natural parameter space $\Lambda_n$. % and  denote the true parameter.Let $\Y^{(n)}$ be a sample generated from the  distribution $\P_{\blambda^\star}^{(n)}$ indexed by the  parameter $\blambda^\star\in\Lambda_n$ and
Then for any parameter $\blambda^\star\in\Lambda_n$  and any distribution $\Xi\in\bP(\Theta_{n,\cM_n})$ that may depend on the sample $\Y^{(n)}$,
\begin{equation}
\label{eq:risk_bound0}
\P_{\blambda^\star}^{(n)}\sbr{\Xi\sbr{n\scd_n^2(\sT_n(\btheta),\blambda^\star)}}
\le \Upsilon_{n,\blambda^\star}(\Xi),
\end{equation}
where we define
\begin{equation}
\label{eq:risk_upper}
\Upsilon_{n,\blambda^\star}(\Xi):= \P_{\blambda^\star}^{(n)}\sbr{ \inf_{\upsilon>0}\del{\frac{1}{\upsilon}\kl(\Xi, \Pi_n(\cdot|\Y^{(n)}))
+\log \int\e^{\upsilon [n\scd_n^2(\sT_n(\btheta),\blambda^\star)]} \Pi_n(\d\btheta|\Y^{(n)})}}.
\end{equation}
\end{proposition}

From \cref{prop:inferior_msvb}, we immediately obtain the following corollary which states that the adaptive variational posterior $_n$ has a tighter risk bound than the model selection variational posterior $_n,m̂_n$.

\begin{corollary}[Risk bound comparison]
Under the same assumption of \cref{prop:compare:risk},
\begin{equation}
\Upsilon_{n,\blambda^\star}(\hXi_n)\le  \Upsilon_{n,\blambda^\star}(\hXi_{n,\hat{m}_n}).
\end{equation}
\end{corollary}

We have shown that the adaptive variational Bayes always leads to a better nonasymptotic risk bound than the model selection variational Bayes. But in some situations, they are equivalent asymptotically. An example is the case where a variational posterior model probability of a specific model dominates the other variational posterior model probabilities.

\begin{proposition}
\label{prop:compare:equi}
Suppose that there exists a sequence of model indices $(m_n^*)_{n\in\bN}\in\bigotimes_{n\in\bN}\cM_n$   such that
\begin{equation*}
\inf_{\blambda^\star\in\Lambda_n^\star}\P_{\blambda^\star}^{(n)}\hgamma_{n,m_n^*}=1+\sco(1).
\end{equation*}
for some $\Lambda_n^\star\subset\Lambda_n$. Then
\begin{equation}
\sup_{\blambda^\star\in\Lambda_n^\star}\P_{\blambda^\star}^{(n)}\norm[1]{\hXi_n-\hXi_{n,\hat{m}_n}}_{\textsc{\textup{TV}}}=\sco(1).
\end{equation}
\end{proposition}

The condition of \cref{prop:compare:equi} is satisfied when model selection consistency holds for the adaptive variational posterior. Thus, under the same assumptions of \cref{col:selection}, we have the same conclusion of the theorem.

\subsection{Numerical examples}

In this subsection, we provide  some numerical examples to illustrate superiority of the adaptive variational Bayes over the model selection variational Bayes.

\subsubsection{Gaussian mixtures}

Consider the Gaussian location-scale mixture model for a bivariate sample:
\begin{equation*}
\Y_1,\dots, \Y_n\iidsim \sum_{k=1}^m\varpi_k\N(\bvartheta_k, \bSigma_k),
\end{equation*}
where $(ϖ_1,…, ϖ_m)∈Δ_m$, $_k∈^2$ and $_k∈_d$ for $k∈[m]$. We  generate a sample with size  $n=200$ from the true model
\begin{align*}
\sum_{k=1}^4&\varpi_k^\star \N(\bvartheta_k^\star, \bSigma_k^\star),
\mbox{ where } (\varpi_1^\star,\dots, \varpi_4^\star)=\frac{1}{10}(3,3,2,2), \\
& \vartheta_1^\star=(0,0), \vartheta_2^\star=(-4,-4), \vartheta_3^\star=(4,4), \vartheta_4^\star=(0,4)
\mbox{ and }  \bSigma_1^\star=\cdots=\bSigma_4^\star=\I_2.
\end{align*}
We applied the two variational Bayes methods to estimate the Gaussian location-scale mixture model based on the generated sample. We considered at most $6$ components model, that is, we set $=[6]$. We impose the following prior distribution, denoted by $Π$. For the number of components $m$, following the way of [Zhang and Gao, 2020], we considered the distribution $Π(m=m')∝exp(-m'logm')$ for  $m'∈[6]$.  Given the number of components $m,$ we impose the normal prior $((0,0)^⊤, 100)$ distribution on each $_k$ independently and the Wishart prior $(10, 0.1)$ on each $_k^-1$ independently.

\begin{figure}
\centering
\includegraphics[scale=0.1]{Figures/mixture_example.png}
\caption{Contour plots for the true and predictive densities by variational Bayes. ``AVB'' and ``MSVB'' are short for the adaptive variational Bayes and  model selection variatonal Bayes, respectively.}
\label{fig:mixture}
\end{figure}

The result is presented in \cref{fig:mixture}. The model selection variation Bayes chooses the 3 components mixture model, which cannot capture the true model structure. But the adaptive variational Bayes smooths the 3 and 4 components mixture models, so this provides better estimate. This experiment is an example illustrating  the  type of overconfident inference that model selection methods can lead to.

\subsubsection{Neural network regression}

We here compare prediction abilities of the two variational posteriors over neural networks: the adaptive variational posterior  \labelcref{eq:dnn_estimate}  and the model selection variational posterior $_n,(K̂_, M̂)$ with $(K̂, M̂):=(K̂_n, M̂_n)∈_(K,M)∈_n_n,(K,M)$. We consider the set of depth-width pairs given as $_n={2,3}×{25, 50, 100}$. We let the prior probability of the model with depth $K$ and width $M$ be proportional to $^-b_0K^2M^2logn$ with $b_0=10^-5$ and set $B_n=√(n)$. We train each of individual variational posteriors for 500 epochs, using the ADAM optimizer [Kingma and Ba, 2014] with a learning rate of $10^-3$.

The performances of the two variational posteriors are evaluated on four UCI data sets: `power', `energy', `boston' and `concrete'. For each data set, we split the data set into training (90\%) and test (10\%) sets randomly 50 times. We compute the posterior mean of the neural network based on 100 samples from the variational posterior and predictions are made based on the computed posterior mean. We calculate the root mean square error (RMSE) on the test set for each of the 50 splits. \cref{fig:deep_regression} illustrates these RMSE values of the two variational posteriors  and we see that the proposed adaptive variational posterior outperforms the model selection variational posterior on all the four UCI data sets.

\begin{figure}
\centering
\includegraphics[scale=0.55]{Figures/deep_reg_rmse.png}
\caption{Boxplots of RMSE values  obtained on UCI datasets across all splits. ``AVB'' and ``MSVB'' are short for the adaptive variational Bayes and  model selection variatonal Bayes, respectively.}
\label{fig:deep_regression}
\end{figure}

\section{Particle-based construction of prior and variational family}
\label{sec:particle}

In this section, we propose a general approach in construction of a prior distribution and variational family, called \textit{particle-based adaptive variational Bayes}. A theoretical advantage of the particle-based approach is that it automatically satisfies a part of the theoretical conditions imposed in our general results. The idea of using a prior distribution on the dicretized parameter space is not new, and has been incorporated in some Bayesian papers [Ghosal et~al., 1997, Ghosal et~al., 2000, Lember and van~der Vaart, 2007]. But their use is only for the original posterior, not for its variational approximation.

%\subsection{Method and theoretical properties}

We first introduce an assumption on the parameter spaces, which is a slightly stronger than \cref{assume:model}.

\begin{assumption}[Models: A stronger version]
\label{assume:particle}
A collection of parameter spaces $\{\Theta_{n,m}^\bullet:m\in\cM_n,n\in\bN\}$ satisfies \cref{assume:model:express} as well as the followings:
\begin{enumerate}[label=(G\arabic*)]
\item \label{assume:particle:metric}
For  any sufficiently large $n\in\bN$ and any $m\in\cM_n$, there exists a metric $\sce_n:\Theta_{n,\cM_n}\times\Theta_{n,\cM_n}\mapsto \R_{\ge0}$ such that
\begin{align}
\scd_n(\sT_n(\btheta_0),\sT_n(\btheta_1))
\le  \sce_{n,m}(\btheta_0,\btheta_1)
\end{align}
for any $\btheta_0,\btheta_1\in\Theta_{n,m}$.
\item  \label{assume:particle:complex}
There exists an absolute constant $\fc_1>0$ such that
\begin{equation}
N_{n,m}:=\scN(\zeta_{n,m}, \Theta_{n,m}, \sce_{n,m})\le \e^{\fc_1n\zeta_{n,m}^2}
\end{equation}
for any $m\in\cM_n$ and any sufficiently large $n\in\bN$.
\end{enumerate}
\end{assumption}

The particle-based approach first suitably discretizes a parameter space  and then design a prior distribution and variational family on the discretized parameter space. For each $n∈$ and $m∈_n$, define
\begin{equation}
\Theta_{n,m}^\bullet:=\cbr[1]{\btheta^\bullet_{n,m,k}:k\in[N_{n,m}]}
\in\argmin_{\mbox{$\Theta^\bullet$: a $\zeta_{n,m}$-net of $\Theta_{n,m}$ w.r.t. $\sce_{n,m}$}}|\Theta^\bullet|,
\end{equation}
that is, $Θ_n,m^∙$ is a minimal set such that $Θ_n,m⊂∪_k=1^N_n,m_Θ_n,m,_n,m(^∙_n,m,k,ζ_n,m)$. We propose to use the collection of the \textit{discretized} parameter spaces ${Θ_n,m^∙}_m∈_n$  instead of the original one ${Θ_n,m}_m∈_n$. To do this, we first show that this collection of discretized parameter spaces satisfies  \cref{assume:model}.

\begin{proposition}
\label{prop:particle:model}
Suppose that \cref{assume:particle} holds. Then the collection of parameter spaces $\{\Theta_{n,m}^\bullet:m\in\cM_n,n\in\bN\}$ satisfies \cref{assume:model:express} with  the approximation errors $\{\eta_{n,m}+\zeta_{n,m}:m\in\cM_n,n\in\bN\}$ and satisfy \cref{assume:model:complex} with  the  estimation errors  $\{\zeta_{n,m}:m\in\cM_n,n\in\bN\}$.
\end{proposition}

On the discretized parameter space $Θ_n,_n^∙:=∪_m∈_nΘ_n,m^∙$, we impose the \textit{discretized} prior distribution given as
\begin{equations}
\label{eq:discrete_prior}
\Pi_{n}&=\sum_{m\in\cM_n}\alpha_m \Pi_{n,m}, \\
\mbox{ where }& \Pi_{n,m}(\cdot)=\sum_{k=1}^{N_{n,m}}\frac{1}{N_{n,m}}\delta(\cdot;\btheta^\bullet_{n,m,k})
\mbox{ and }\alpha_{n,m}=\frac{\e^{-b_0n\zeta_{n,m}^2}}{\sum_{m'\in\cM_n}\e^{-b_0n\zeta_{n,m'}^2}}
\end{equations}
for some $b_0>0$. That is, we consider the discretized prior distribution which imposes equal probability mass on the minimal $ζ_n,m$-net of the parameter space $Θ_n,m$. For the discretized prior, the prior concentration condition \labelcref{assume:prior:truth} is met.

\begin{proposition}
\label{prop:particle:prior}
Under \cref{assume:particle}, the discretized prior distribution  \labelcref{eq:discrete_prior} satisfies that
\begin{align}
\Pi_{n,m}\del{\btheta^\bullet\in\Theta_{n,m}^\bullet:\scd_n(\sT_n(\btheta^\bullet),\blambda^\star)\le \eta_{n,m}+ \zeta_{n,m}}\ge \e^{-\fc_1n \zeta_{n,m}^2},
\end{align}
for any $m\in\cM_n$ and sufficiently large $n\in\bN$ for some absolute constant $\fc_1>0$.
\end{proposition}

We consider the variational family given by
\begin{equation}
\label{eq:particle_family}
\bV_n=\cbr{\sum_{m\in\cM_n}\gamma_m\Xi_m:(\gamma_m)_{m\in\cM_n}\in\Delta_{|\cM_n|},\Xi_m\in\bD_{Q_{n,m}}(\Theta_{n,m}^\bullet)},
\end{equation}
where for each $m∈_n$, we denote
\begin{equation*}
\bD_{Q_{n,m}}(\Theta_{n,m}^\bullet):=\cbr{\sum_{q=1}^{Q_{n,m}}\omega_q\delta(\cdot;\bpsi_q):(\omega_q)_{q\in[Q]}\in\Delta_{Q_{n,m}},\bpsi_q\in\Theta_{n,m}^\bullet}
\end{equation*}
with the number of particles $Q_n,m∈[N_n,m]$.
In words, each variational distribution in the variational family $_Q(Θ_n,m^∙)$ consists of \textit{particles} which are  Dirac-delta measures. In this sense, we call the variational family  \labelcref{eq:particle_family} a \textit{particle variational family}. A main theoretical advantage of the particle variational family  \labelcref{eq:particle_family} is that it automatically satisfies \cref{assume:variational}.

\begin{proposition}
\label{prop:particle:verror}
Suppose that \cref{assume:particle} holds. Then with the discretized prior distribution \labelcref{eq:discrete_prior}, the particle variational family  \labelcref{eq:particle_family} with arbitrary $(Q_{n,m})_{m\in\cM_n}\in\bigotimes_{m\in\cM_n}[N_{n,m}]$ satisfies that
\begin{equation}
\sup_{\blambda^\star\in\Lambda_n^\star}\P_{\blambda^\star}^{(n)}\sbr{\inf_{\Xi\in \bV_n}\kl(\Xi,\Pi_n(\cdot|\Y^{(n)}))}\le \fc_1n\epsilon_n^2
\end{equation}
for any sufficiently large $n\in\bN$ for some absolute constant $\fc_1>0$.
\end{proposition}

\begin{remark}
We choose a Dirac-delta measure as a distribution for each particle. In practice, continuous distributions such as normal distribution may be preferred because it leads to a smooth variational posterior distribution which can provides more flexible inference tools. For example, one can use multivariate normal distributions with  isotropic covariances [Gershman et~al., 2012] or with low rank structured covariances [Miller et~al., 2017]. However, theoretical analysis on the variational approximation error for this continuous relaxation of the particle variational family is not trivial as \cref{prop:particle:verror} for the Dirac-delta particles. We leave it as future work.
\end{remark}

Since \cref{assume:particle} implies  \cref{assume:model,assume:prior,assume:variational}, the next corollary follows immediately from \cref{thm:conv}.

\begin{corollary}[Adaptive contraction rate of particle-based variational posterior]
\label{col:particle_contraction}
Suppose that  \cref{assume:exprmt,assume:rate,assume:particle}  hold. Then with the discretized prior distribution \labelcref{eq:discrete_prior} and the particle variational family \labelcref{eq:particle_family} with arbitrary $(Q_{n,m})_{m\in\cM_n}\subset\bigotimes_{m\in\cM_n}[N_{n,m}]$, we have
\begin{equation}
\label{eq:particle_conv}
\sup_{\blambda^\star\in\Lambda_n^\star}\P_{\blambda^\star}^{(n)}\sbr{\hXi_n\del{\scd_n(\sT_n(\btheta),\blambda^\star)\ge A_n\epsilon_{n}}}=\sco(1)
\end{equation}
for any diverging sequence $(A_n)_{n\in\bN}\to\infty$.
\end{corollary}

The variational posterior obtained from the particle-based adaptive variational Bayes includes the following two estimators as particular instances. If $Q_n,m=1$, then the variational posterior $_n,m$ collapses to the Dirac-delta measure at the maximum likelihood estimator (MLE) over the discretized parameter space $Θ_n,m^∙$. On the other hand, if  $Q_n,m=N_n,m$, then  $_n,m(^∙)∝_n(^∙,^(n))$ for any $^∙∈Θ_n,m^∙$, that is, the variationl posterior is the same as the original posterior induced by the discretized prior distribution. Since our theory does not impose any restriction on the number of particles, we have provided the theoretical guarantees for the MLE and the original posterior also.

\begin{example}[Adaptive variational deep learning]
Recall the adaptive variational deep learning procedure  proposed in \cref{sec:dnn}. In this example, we provide a discretizing method for the neural network parameter spaces, which satisfies \cref{assume:particle}. By \cref{lem:dnn_lip}, for any $(K,M)\in\cM_n$ and any $\btheta_1,\btheta_2\in \Theta^{\le B_n}_{n,(K,M)}$,
\begin{equation*}
\norm{\net(\btheta_1)-\net(\btheta_2)}_{\infty} \le C_{n,(K,M)}\abs{\btheta_1-\btheta_2}_\infty
\mbox{ with } C_{n,(K,M)}:=K(B_n(M+1))^{K},
\end{equation*}
provided that $B_n\ge 1$.
From the above display, \cref{assume:particle:metric} is met with the metric $\sce_{n,(K,M)}(\btheta_1,\btheta_2)=C_{n,(K,M)}\abs{\btheta_1-\btheta_2}_\infty$. We construct discretized parameter spaces as follows. For each $(K,M)\in\cM_n$, we let
\begin{equation*}
\Theta_{n,(K,M)}^{\bullet}:=\cbr{\btheta^\bullet\in[-B_n,B_n]^{p_{(K,M)}}:\frac{C_{n,(K,M)}}{\zeta_{n,(K,M)}}\btheta^\bullet\in\bZ^{p_{(K,M)}}},
%\Theta_{(K,M)}^{\bullet,\le B_n}
%\Theta_{n,(K,M)}^{\bullet}:=\cbr{\del{\frac{\zeta_{n,(K,M)}}{C_{n,(K,M)}}Z_j}_{j\in[p_{(K,M)}]}\in\R^{p_{(K,M)}}:Z_j\in\bZ, |Z_j|\le \frac{C_{n,(K,M)}}{\zeta_{n,(K,M)}}B_n}.%\Theta_{(K,M)}^{\bullet,\le B_n}
\end{equation*}
where $\zeta_{n,(K,M)}:=KM\sqrt{\log n/n}$.
Then, it is clear that $ \Theta_{n,(K,M)}^{\bullet}$ is a $\zeta_{n,(K,M)}$-net of the parameter space $\Theta_{(K,M)}^{\le B_n}$ with respect to the metric $\sce_{n,(K,M)}$. Moreover, if we assume that $\max_{(K,M)\in\cM_n}(K\vee M)\lesssim n$ and $1\le B_n\lesssim n^{\nu}$ for some $\nu>0$ as we did in \cref{thm:dnn:oracle}, we have
\begin{align*}
&\le \del{2B_n \frac{C_{n,(K,M)}}{\zeta_{n,(K,M)}}}^{p_{(K,M)}}
\le \del{\frac{4\sqrt{n}B_n^{K+1}(M+1)^{K-1}}{\sqrt{\log n}}}^{(d+1)KM^2}\\
\le\e^{\fc_2K^2M^2\log n}=\e^{\fc_2n\zeta_{n,(K,M)}^2}
\end{align*}
for any $(K,M)\in\cM_n$ for some absolute constants $\fc_1>0$ and $\fc_2>0$, which verifies \cref{assume:particle:complex}. Therefore, all the theoretical properties in \cref{sec:dnn} can also be attained by the particle-based adaptive variational posterior over neural networks based on the above construction of parameter spaces.
\end{example}

In \cref{appen:particle_compute}, we provide a general computation strategy for the particle-based adaptive variational Bayes.

\section{Variational quasi-posteriors}
\label{sec:quasi}

In this section, we analyze variational approximation to \textit{quasi-posterior} that is defined as follows. For a given \textit{quasi-likelihood function} $_n^♮:Λ_n×_n↦_≥0$, we define the quasi-posterior given the sample $^(n)$ as
\begin{equation}
\Pi_n^\natural(\btheta\in\cT|\Y^{(n)}):=\frac{\int_{\cT} \sq_n^\natural(\sT_n(\btheta), \Y^{(n)})\d\Pi_n(\btheta)}{\int \sq_n^\natural(\sT_n(\btheta), \Y^{(n)})\d\Pi_n(\btheta) }
\end{equation}
for any measurable subset $⊂Θ_n,_n.$ Then \textit{variational quasi-posterior}, denoted by $_n^♮$, is the variational approximation of the quasi posterior, that is,
\begin{equation}
\hXi_n^\natural\in\argmin_{\Xi\in\bV_n}\kl(\Xi,\Pi_n^\natural(\cdot|\Y^{(n)})).
\end{equation}

The quasi-posterior has been used in a number of applications for various purposes. We provide some examples in order.
\begin{itemize}
\item \textit{Fractional posteriors}.  The fractional posterior is a quasi-posterior where the fractional likelihood $(\sq_n)^{\kappa}$ for  $\kappa\in(0,1)$ is used as the quasi-likelihood. This replacement gives a theoretical advantage that its concentration properties can be established without a testing construction [Bhattacharya et~al., 2019, Martin et~al., 2017, Lee et~al., 2019].  Moreover, there are empirical and theoretical evidences that the fractional posterior improves robustness to model misspecification [Gr{\"u}nwald, 2011, Gr{\"u}nwald and Van~Ommen, 2017, Miller and Dunson, 2018, Medina et~al., 2021].

\item \textit{Gibbs posteriors}. A loss function based quasi-likelihood, which measures the ``difference'' between the parameter of interest and the sample, is used. This is particularly useful when specification of the likelihood function is not straightforward and/or is too restrictive for defining the true model. A quantile regression is a representative example. Several quasi-likelihood functions including  check loss based quasi-likelihood (also known as the asymmetric Laplace likelihood) [Yang et~al., 2016] and empirical likelihood [Yang and He, 2012] have been proposed. General theories for contraction properties of Gibbs posteriors were established by Atchad{\'e}, 2017, Atchad{\'e} and Bhattacharyya, 2018 and Syring and Martin, 2020.

\item \textit{Robust posteriors}. It is known that the original posterior is not robust  to model misspecification. Some robust surrogates for the original posterior have been proposed.   Baraud and Birg{\'e}, 2020 proposed a robust quasi-posterior where the logarithmic function in the log-likelihood ratios is replaced by some bounded function. Ch{\'e}rief-Abdellatif and Alquier, 2020 and Matsubara et~al., 2021 considered a quasi-likelihood function based on the maximum mean discrepancy and Stein discrepancy measures, respectively, which are more robust to model misspecification.
\end{itemize}

\subsection{Adaptive contraction rates of variational quasi-posteriors}
\label{sec:quasi:contraction}

The main question of this subsection is what conditions on the statistical experiment and quasi-likelihood under which adaptive optimal contraction properties of the variational quasi-posteriors are achieved.

\begin{assumption}[Quasi-likelihood]
\label{assume:quasi}
A sequence of quasi-likelihood functions $(\sq_n^\natural)_{n\in\bN}$ with $\sq_n^\natural:\Lambda_n\times \fY_n\mapsto\R_{\ge0}$ satisfy the followings:
\begin{enumerate}[label=(H\arabic*)]
\item  \label{assume:quasi:learning}
(Learning ability) There exist absolute constants $\fc_1>0$ and $\fc_2>$ such that
\begin{equation}
\P_{\blambda^\star}^{(n)}\del{\frac{\sq_n^\natural(\blambda,\Y^{(n)})}{\sq_n^\natural(\blambda^\star,\Y^{(n)})}}\le \fc_1\e^{-\fc_2n\scd_n^2(\blambda,\blambda^\star)}
\end{equation}
for any $\blambda^\star\in\Lambda_n^\star$,  $\blambda\in\Lambda_n$ and sufficiently large $n\in\bN$.

\item  \label{assume:quasi:bound}
(Boundedness) There exist absolute constants $\barho>0$, $\fc_3>0$ and $\fc_4>0$ such that
\begin{equation}
\P_{\blambda^\star}^{(n)}\del{\frac{\sq_n^\natural(\blambda^\star,\Y^{(n)})}{\sq_n^\natural(\blambda,\Y^{(n)})}}^{\barho}\le \fc_3\e^{\fc_4n\scd_n^2(\blambda,\blambda^\star)}
\end{equation}
for any $\blambda^\star\in\Lambda_n^\star$,  $\blambda\in\Lambda_n$ and sufficiently large  $n\in\bN$.
\item  \label{assume:quasi:approximation}
(Variational approximation)
There exists an absolute constant $\fc_5>0$ such that
\begin{equation}
\sup_{\blambda^\star\in\Lambda_n^\star}\P_{\blambda^\star}^{(n)}\sbr{\inf_{\Xi\in\bV_n}\kl(\Xi,\Pi_n^\natural(\cdot|\Y^{(n)}))}\le \fc_5L_nn\epsilon_n^2
\end{equation}
for any sufficiently large $n\in\bN$.
\end{enumerate}
\end{assumption}

We discuss implications of the conditions  of  \cref{assume:quasi}. \cref{assume:quasi:learning} requires that for the quasi-likelihood of the true parameter is larger than those for the other parameters with high probability. This implies that the quasi-likelihood should suitably measure the  goodness-of-fit of the parameter to the sample. Contrastively, this assumption also requires the Bayesian update through the quasi-likelihood do not  reflect the information of the sample too much, that is, we demand some ``robustness'' to the sample. The original likelihood function $_n$, which learns all information of the sample, satisfy $_^⋆^(n)[_n(,^(n))/_n(^⋆,^(n))]=1$ for any $^⋆,∈Λ_n$. This robust property makes the proof quite simpler. \cref{assume:quasi:bound} plays a similar role of \cref{eq:tail2} in \cref{assume:exprmt:tail}, which requires that a distribution of each quasi-likelihood ratio is not heavy-tailed.

\cref{assume:quasi:approximation} is a variational quasi-posterior version of \cref{assume:variational}. The next lemma, which corresponds to \cref{lem:verror}, is useful to check \cref{assume:quasi:approximation}.

\begin{lemma}
\label{lem:quasi:verror}
Under  \cref{assume:quasi:learning,assume:quasi:bound},
\begin{align*}
\P_{\blambda^\star}^{(n)}\sbr{ \kl(\Xi, \Pi_n^\natural(\cdot|\Y^{(n)}))}
&\le \kl(\Xi,\Pi_n) +\Xi\sbr{\P_{\blambda^\star}^{(n)}\log(\frac{\sq_{n}^\natural(\blambda^\star,\Y^{(n)})}{\sq_n^\natural(\sT_{n}(\btheta),\Y^{(n)})})} + \log(\fc_1)\\
&\le  \kl(\Xi,\Pi_n) +\frac{\fc_4}{\barho}\Xi\sbr{n\scd_n^2(\sT_n(\btheta),\blambda^\star)}+\log\del[1]{\fc_1\fc_3^{1/\barho}}
\end{align*}
for any $\blambda^\star\in\Lambda_n$ and $\Xi\in\bP(\Theta_{n,\cM_n})$, where $\fc_1>0$, $\fc_3>0$, $\fc_4>0$ and $\barho>0$ are the absolute constants appearing in  \cref{assume:quasi:learning,assume:quasi:bound}.
\end{lemma}

\begin{remark}
The fractional likelihood is encompassed in our quasi-likelihood framework. Indeed,  under \cref{assume:exprmt}, the fractional likelihood $(\sq_n)^{\kappa}$ for  $\kappa\in(0,1-\rho_\circ)$ satisfies \cref{assume:quasi:learning,assume:quasi:bound}. This claim is verified as follows. Note first that by the assumption $1-\kappa>\rho_\circ$,
\begin{align*}
\P_{\blambda^\star}^{(n)}\del{\frac{(\sq_n(\blambda,\Y^{(n)}))^{\kappa}}{(\sq_n(\blambda^\star,\Y^{(n)}))^{\kappa}}}
&=\exp(-\kappa\sD_{1-\kappa}(\P^{(n)}_{\blambda^\star}, \P^{(n)}_{\blambda}))
\le \exp(-\kappa\sD_{\rho_\circ}(\P^{(n)}_{\blambda^\star}, \P^{(n)}_{\blambda})).
\end{align*}
Hence, the first inequality of \cref{assume:exprmt:tail} implies \cref{assume:quasi:learning}. On the other hand, by the second inequality of \cref{assume:exprmt:tail}, \cref{assume:quasi:bound} is  met with $\barho=(\rho_\blacklozenge-1)/\kappa>0$.
\end{remark}

We have the following contraction properties of the variational quasi-posteriors.

\begin{theorem}[Adaptive contraction rate of variational quasi-posterior]
\label{thm:quasi:conv}
Suppose that \cref{assume:prior,assume:quasi} hold. Then
\begin{equation}
\label{eq:conv_quasi}
\sup_{\blambda^\star\in\Lambda_n^\star}\P_{\blambda^\star}^{(n)}\sbr{\hXi_n^\natural\del{\scd_n(\sT_n(\btheta),\blambda^\star)\ge A_n\sqrt{L_n}\epsilon_{n}}}=\sco(1)
\end{equation}
for any diverging sequence $(A_n)_{n\in\bN}\to\infty$.
\end{theorem}

\subsection{Application to stochastic block models}
\label{sec:quasi:sbm}

In this subsection, we consider a stochastic block model example to illustrate a situation where the use of quasi-likelihood provides some theoretical advantages.

\subsubsection{Random graph experiment}

We first introduce some notation useful throughout this section. For $n∈$, let $_n:={(i,j)∈[n]^2:i>j}$ be the set of 2-dimensional indices with the first larger than the second. For a set $I$, let $I^_n:=(b_i,j)_(i,j)∈_n:b_i,j∈I$  be the set of lower triangular array of $I$-valued elements. Recall the definitions of the spaces $^m$ and $_n,m$ given in \Cref{eq:def_sm,eq:def_znm}, respectively.

Suppose that we observe a triangular array $^():=(Y_i,j)_(i,j)∈_n∈_:={0,1}^_n$ which represents  existence of connections between $n$-nodes. Here and afterward, we index all the notations by the number of connections $:=|_n|=n(n-1)/2$ instead of the number of nodes $n$, to make the notation consistently.

The $(i,j)$-th element $Y_i,j$ is encoded to be 1 if the $i$-th and $j$-th nodes are connected and to be 0 otherwise. We model the sample using the set of distributions
\begin{equation*}
\bP_{\ban}:=\cbr{\P^{(\ban)}_{\bOmega}:=\prod_{(i,j)\in\cL_n}\Ber(\Omega_{i,j}):\bOmega\in[0,1]^{\cL_n}},
\end{equation*}
i.e., $_:=({0,1}^_n;_,[0,1]^_n)$ where $_:{0,1}^_n×[0,1]^_n↦_≥0$ is the Bernoulli likelihood function  defined as  $_(, ^(n))=∏_(i,j)∈_nΩ_i,j^Y_i,j(1-Ω_i,j)^1-Y_i,j$. We call the parameter $∈[0,1]^_n$ the \textit{connectivity probability array}. A general goal is to provide a good estimate of the true connectivity probability array in terms of the (scaled) $^2$ error defines as
\begin{equation*}
\scd_{\ban,2}(\bOmega_0,\bOmega_1)=\del[3]{\frac{2}{n(n-1)}\sum_{(i,j)\in\cL_n}(\Omega_{0,i,j}-\Omega_{1,i,j})^2}^{1/2}
\end{equation*}
for $_0,_1∈[0,1]^_n$.

\subsubsection{Gaussian  quasi-likelihood}

Here we propose to use the \textit{Gaussian quasi-likelihood}
\begin{equation}
\sq_{\ban}^{\natural}(\bOmega, \Y^{(\ban)})=\prod_{(i,j)\in\cL_n}\e^{-(Y_{i,j}-\Omega_{i,j})^2}
\end{equation}
instead of the Bernoulli likelihood $_$. This is because we want to remove the unpleasant feature of the Bernoulli likelihood that the KL divergence $(^()__0,^()__1)$ can not be bounded by the $^2$ error between $_1$ and $_0$ due to its boundary behavior. To be specific, suppose that there are only $2$ nodes, i.e., 1 connectivity. Then the KL divergence
\begin{align*}
\kl(\P^{(1)}_{\Omega_0},\P^{(1)}_{\Omega_1})
\end{align*}
is  larger than $C(Ω_0-Ω_1)^2$ for any large $C>0$ when $Ω_1$ is very close to 0 or 1. Thus related Bayesian literature restricts the parameter space to be $[ϰ,1-ϰ]^_n$ for some sufficiently small $ϰ∈(0,1/2)$ to avoid such a problem [Ghosh et~al., 2020, Jiang and Tokdar, 2021]. Here, we show that the use of the Gaussian  quasi-likelihood detours this issue without the parameter restriction. Similarly to ours, Gao et~al., 2020 showed theoretical optimality of the quasi-posterior distribution with the Gaussian quasi-likelihood but by using a completely different proof technique. Also, they did not study its variational approximation.

In the next lemma, we show that the Gaussian  quasi-likelihood satisfies \cref{assume:quasi:learning,assume:quasi:bound}.

\begin{lemma}
\label{lem:sbm:quasi}
Let $\barho>0$. Then for any $\bOmega_0\in[0,1]^{\cL_n}$ and  $\bOmega_1\in[0,1]^{\cL_n}$,
\begin{align}
\P_{\bOmega_0}^{(\ban)}\del{\frac{\sq_{\ban}^{\natural}(\bOmega_1,\Y^{(\ban)})}{\sq_{\ban}^{\natural}(\bOmega_0,\Y^{(\ban)})}}
& \le \e^{-\frac{1}{2}\ban\scd_{\ban,2}^2(\bOmega_0,\bOmega_1)}
\label{eq:sbm_quasi_learning}\\
\P_{\bOmega_0}^{(\ban)}\del{\frac{\sq_{\ban}^{\natural}(\bOmega_0,\Y^{(\ban)})}{\sq_{\ban}^{\natural}(\bOmega_1,\Y^{(\ban)})}}^{\barho}
&\le \e^{\frac{\barho(\barho+2)}{2}\ban\scd_{\ban,2}^2(\bOmega_0,\bOmega_1)}.      \label{eq:sbm_quasi_bound}
\end{align}
\end{lemma}

\subsubsection{Model, prior and variational posterior}

The stochastic block model assumes that the nodes are grouped into several communities and the connectivity probabilities between nodes depend exclusively on their community membership. That is, the connectivity probability  array $$ is assumed to be an output of the map $_:∪_m=1^n(^m×_n,m)→[0,1]^_n$ defined by
\begin{equation*}
\sT_{\ban}(\U, \Z)=(\z_i^\top\U\z_j)_{(i,j)\in\cL_n}%=(\Z\U\Z^\top)_{\cL_n}:
\end{equation*}
for $∈^m$ and $∈_n,m$ for some $m∈[n]$.
The stochastic block model is the collection of the models
\begin{equation*}
\bM_{\ban,m}:=\cbr{\P^{(\ban)}_{\sT_{\ban}(\U, \Z)}\in\bP_{\ban}:(\U, \Z)\in\cS^m\times \cZ_{n,m}}
\end{equation*}
for $m∈[n]$, with the parameter spaces $Θ_,m:=^m×_n,m_m∈[n]$.

Let $_⊂[n] $ be a set of some model indices, that is, a set of the possible number of communities. We impose the prior distribution given by
\begin{equations}
\Pi_{\ban}=&\sum_{m\in\cM_{\ban}}\alpha_{\ban,m}\Pi_{\ban,m},  \\
&\mbox{ where }\Pi_{\ban,m}:=\bigotimes_{(k,h)\in[m]^2:k\le h}\Unif(0,1)\bigotimes_{i=1}^{n}\Cat(m^{-1}\one_m)\\
&\mbox{ and }    \alpha_{\ban,m}:=\frac{\e^{-b_0(m^2\log n+n\log m)}}{\sum_{m'\in\cM_n}\e^{-b_0((m')^2\log n+n\log m')}}
\end{equations}
for some $b_0>0$.
For each model index $m∈_n$, we consider the variational family such as
\begin{equation*}
\bV_{\ban,m}:=\cbr{\bigotimes_{(k,h)\in[m]^2:k\le h}\Unif(\psi_{1,k,h},\psi_{2,k,h})\bigotimes_{i=1}^{n}\Cat(\bnu_i):0\le\psi_{1,k,h}\le\psi_{2,k,h}\le1 ,\bnu_i\in\Delta_m}.
\end{equation*}
Since the parameter spaces  $[0]Θ_,m_m∈_n$  are disjoint, due to \cref{thm:compute}, the variational quasi posterior is given by
\begin{equations}
\label{eq:sbm_variational_posterior}
\hXi_{\ban}^\natural=&\sum_{m\in\cM_{\ban}}\hgamma_{\ban,m}\hXi_{\ban,m},\\
&\mbox{ where }  \hXi_{\ban,m}\in\argmin_{\Xi\in\bV_{\ban,m}} \scE\del[1]{\Xi, \Pi_{\ban, m}, \sq_{\ban}^{\natural}, \Y^{(\ban)}}\\
&\mbox{ and }  \hgamma_{\ban,m}=\frac{\alpha_{\ban,m}\e^{-\scE\del[1]{\hXi_{\ban,m}, \Pi_{\ban,m}, \sq_{\ban}^{\natural}, \Y^{(\ban)}}}}
{\sum_{m'\in\cM_n}\alpha_{\ban,m'} \e^{-\scE\del[1]{\hXi_{\ban,m'}, \Pi_{\ban,m'}, \sq_{\ban}^{\natural}, \Y^{(\ban)}}}}
\end{equations}

\subsubsection{Estimation of connectivity probability array}

In this subsection, we present contraction properties of the adaptive variational quasi-posterior with the Gaussian quasi-likelihood. The following oracle result provides a convenient theoretical tool to derive contraction rates in the related problems.

\begin{theorem}[Oracle contraction rate, connectivity probability array]
\label{thm:sbm:oracle}
Let $(\Lambda_{\ban}^\star)_{n\in\bN}$ be a sequence of sets such that $\Lambda_{\ban}^\star\subset[0,1]^{\cL_n}$.
The variational quasi-posterior \labelcref{eq:sbm_variational_posterior} satisfies
\begin{equation}
\sup_{\bOmega^\star\in\Omega_{\ban}^\star}\P_{\bOmega^\star}^{(\ban)}\sbr{\hXi_{\ban}\del{\scd_{\ban,2}(\sT_{\ban}(\U,\Z),\bOmega^\star)\ge A_n\epsilon_n(\Lambda_{\ban}^\star)}}=\sco(1)
\end{equation}
for any diverging sequence $(A_n)_{n\in\bN}\to+\infty$,
\begin{equation}
\epsilon_n(\Lambda_{\ban}^\star):=\inf_{m\in\cM_{\ban}}\cbr{\sup_{\bOmega^*\in\Lambda_{\ban}^\star}\inf_{(\U,\Z)\in\Theta_{\ban,m}}\scd_{\ban,2}(\sT_{\ban}(\U,\Z),\bOmega^*)+\sqrt{\frac{m^2}{n^2}\log n+\frac{\log m}{n}}}.
\end{equation}
\end{theorem}

The estimation error part of our oracle contraction rate in \cref{thm:sbm:oracle} is $log^1/2n$ times slower than the optimal rate $√(m^2/n^2 + logm/n)$ [Gao et~al., 2015]. This is because we use the standard ``prior mass and testing'' technique commonly used in the Bayesian nonparametric literature, which frequently leads to $logn$ suboptimality. Similarly, in Ghosh et~al., 2020, the same additional $log^1/2n$ term appeared in the contraction rate of their original posterior distribution over stochastic block models. But if we assume that the infimum is attained for $m≲√(n)$, then since $m^2logn≲n^2logm$, our rate becomes the same as the optimal rate. Jiang and Tokdar, 2021 used such an assumption and could avoid the  $logn$ suboptimality problem.

Our first corollary is about the situation where the sample is generated by the stochastic block model itself. This is a trivial consequence of \cref{thm:sbm:oracle}, and so we omit the proof.

\begin{corollary}[Stochastic block model]
Let $(m^\star_n)_{n\in\bN}$ be an arbitrary sequence such that $m^\star_n\in[n]$. Then   the variational quasi-posterior \labelcref{eq:sbm_variational_posterior} with the model index set $\cM_{\ban}=[n]$ satisfies
\begin{equation}
\sup_{(\U^\star,\Z^\star)\in\Theta^\star_{\ban,m^\star_n}}\P_{\sT_{\ban}(\U^\star,\Z^\star)}^{(\ban)}\sbr{\hXi_{\ban}\del{\scd_{\ban,2}(\sT_{\ban}(\U,\Z),\sT_{\ban}(\U^\star,\Z^\star))\ge A_n  \epsilon_n^{\textup{SBM}}(m_n^\star)}}=\sco(1)
\end{equation}
for any diverging sequence $(A_n)_{n\in\bN}\to+\infty$, where
\begin{equation*}
\epsilon_n^{\textup{SBM}}(m_n^\star):=\sqrt{\frac{(m^\star_n)^2}{n^2}\log n+\frac{\log m^\star_n}{n}}.
\end{equation*}
\end{corollary}

The second example is on  estimation of smooth \textit{graphton}s. We first introduce the related definitions and notations. For a function $f:[0,1]^2↦[0,1]$ and $n$-dimensional $[0,1]$-valued vector $:=(x_i)_i∈[n]∈[0,1]^n$, we define the corresponding graphon as
\begin{equation*}
\textup{Graphon}_n(f,\x):=\del[1]{f(x_i,x_j)}_{(i,j)\in\cL_n}\in[0,1]^{\cL_n}.
\end{equation*}
For a given class $$ of some $[0,1]$-valued functions supported on $[0,1]^2$, we denote
\begin{equation}
\Lambda_{\ban}^{\textsc{G}}(\widetilde\cF):=\cbr{\textup{Graphon}_n(f,\x)\in[0,1]^{\cL_n}:f\in\widetilde\cF, \x\in[0,1]^n}.
\end{equation}
Recall that we have defined $^β,2, F_0_[0,1]$ as the set of $[0,1]$-valued H\"older $β$-smooth functions supported on $[0,1]^2$. The next corollary provides the optimality of our adaptive variational posterior.

\begin{corollary}[Smooth graphon]
\label{col:graphon}
Let $\beta>0$ and $F_0>0$. Then the variational quasi-posterior \labelcref{eq:sbm_variational_posterior} with the model index set given by
\begin{equation*}
\cM_{\ban}=\cbr{\floor{n^{k/\log n}}:k\in\sbr{0:\ceil{\log n/2}} }
\end{equation*}
\begin{equation}
\sup_{\bOmega^\star\in  \Lambda_{\ban}^{\textup{\textsc{G}}}\del{\cH^{\beta,2, F_0}_{[0,1]}}}\P_{\bOmega^\star}^{(\ban)}\sbr{\hXi_{\ban}\del{\scd_{\ban,2}(\sT_{\ban}(\U,\Z),\bOmega^\star)\ge A_nn^{-\del{\frac{\beta}{\beta+1}}\wedge\frac{1}{2}}\log^{\frac{1}{2}} n}}=\sco(1)
\end{equation}
for any diverging sequence $(A_n)_{n\in\bN}\to+\infty$.
\end{corollary}

The minimax optimal rate for the smooth graphon estimation is $n^-β/(β+1)+√(logn/n)$ [Gao et~al., 2015]. Our contraction rate in \cref{col:graphon} is the exactly same as the optimal rate for sufficiently smooth graphons with $β≥1$ but is slightly slower  with the extra $log^1/2 n$ term than the optimal rate when $β∈(0,1)$.

\subsection{Application to regression with sub-Gaussian errors}
\label{sec:quasi:subgauss}

A real-valued random variable $Y∼$ or its distribution $$ is said to be \textit{sub-Gaussian} with mean $ϑ∈$ and variance proxy $ς^2>0$, and  denoted by $∈subG(ϑ,ς^2)$, if $Y=ϑ$ and $^t(Y-ϑ)≤^ς^2t^2/2$ for any $t∈$.

We here consider a nonparametric regression experiment with sub-Gaussian errors, where each real-valued output $Y_i$ is a sub-Gaussian random variable with mean $f(_i)$ and variance proxy $ς^2>0$ for some regression function $f∈^d$ and given covariates $(_i)_i∈[n]∈([0,1]^d)^⊗n$. In other words, we assume that a distribution of $(Y_i)_i∈[n]$ belongs to the set
\begin{equation*}
\text{subG}_n(f^\star,\varsigma^2):=\text{subG}_n(f^\star,\varsigma^2;(\x_i)_{i\in[n]})
\end{equation*}
This is a more general situation than the regression experiment with Gaussian errors considered in \cref{sec:dnn:reg}.

A difficulty of Bayesian inference for the regression experiment with sub-Gaussian errors is that we cannot construct a likelihood function of an explicit form. Fortunately, we can obtain optimal contraction of the adaptive variational posterior when we use the Gaussian quasi-likelihood, which is given as
\begin{equation}
\sq_{n}^{\natural,\kappa}(f, \Y^{(n)})=\prod_{i=1}^n\exp(-\frac{\kappa}{2}(Y_i-f(\x_i))^2)
\end{equation}
for $κ>0$.
Following some of the related literature [Gr{\"u}nwald, 2011, Syring and Martin, 2020, Bhattacharya and Martin, 2020], we call the constant $κ>0$ a \textit{learning rate}, because it determines how much the quasi-posterior learn from the current sample. We show that if the learning rate is sufficiently small, we obtain the desired contraction of the variational quasi-posterior.

The next lemma  proves that if we use a sufficiently small learning rate such that $κ∈(0,1/ς^2)$, \cref{assume:quasi} will be met.

\begin{lemma}
\label{lem:subgauss}
Let $\varsigma>0$, $\kappa>0$ and $\barho>0$. Then for any $f,f^\star\in\cF^d$ and $\P_{\star}^{(n)}\in\textup{subG}_n(f^\star,\varsigma)$,
\begin{align}
\P_{\star}^{(n)}\del{\frac{\sq_{n}^{\natural,\kappa}(f, \Y^{(n)})}{\sq_{n}^{\natural,\kappa}(f^\star, \Y^{(n)})}}
&\le \e^{-\frac{1}{2}\kappa(1-\kappa\varsigma^2)n\nnorm{f-f^\star}^2}, \label{eq:subg_quasi_learning}\\
\P_{\star}^{(n)}\del{\frac{\sq_{n}^{\natural,\kappa}(f^\star,\Y^{(n)})}{\sq_{n}^{\natural,\kappa}(f,\Y^{(n)})}}^{\barho}
&\le \e^{\frac{1}{2}\barho\kappa(1+\barho\kappa\varsigma^2)n\nnorm{f-f^\star}^2}.  \label{eq:subg_quasi_bound}
\end{align}
\end{lemma}

By \cref{thm:quasi:conv} together with \cref{lem:subgauss}, an adaptive optimality result can be obtained. As an example, we consider the adaptive variational quasi-posterior over neural networks with the Gaussian quasi-likelihood, which is given as
\begin{equations}
\label{eq:dnn_variational_quasi_posterior}
\hXi_n^{\natural,\kappa}=&\sum_{(K,M)\in\cM_n}\hgamma_{n,(K,M)}^{\natural,\kappa}\hXi_{n,(K,M)}^{\natural,\kappa},\\
&\mbox{ where }\hXi_{n,(K,M)}^{\natural,\kappa}\in\argmin_{\Xi\in\bV_{n,(K,M)}}  \scE\del[1]{\Xi, \Pi_{n,(K,M)}, \sq_{n}^{\natural,\kappa},\Y^{(n)}}\\
&\mbox { and } \hgamma_{n,(K,M)}^{\natural,\kappa}=\frac{ \alpha_{n,(K,M)} \e^{-\scE\del[1]{\hXi_{n,(K,M)}^{\natural,\kappa}, \Pi_{n,(K,M)},\sq_{n}^{\natural,\kappa},\Y^{(n)}}}}
{\sum_{(K',M')\in\cM_n}\alpha_{n,(K',M')} \e^{-\scE\del[1]{\hXi_{n,(K',M')}^{\natural,\kappa},\Pi_{n,(K',M')},\sq_{n}^{\natural,\kappa},\Y^{(n)}}}}.
\end{equations}
Here, the prior distribution and the variational family are given as \labelcref{eq:dnn_prior} and \labelcref{eq:dnn_variational_family}, respectively. The following corollary presents the contraction rate of the variational quasi-posterior in the above display.

\begin{corollary}[Oracle contraction rate, sub-Gaussian regression]
Under the same assumptions of \cref{thm:dnn:oracle}, the variational quasi-posterior \labelcref{eq:dnn_variational_quasi_posterior} over neural networks with $\kappa\in(0,1/\varsigma^2)$ satisfies
\begin{equation}
\sup_{f^\star\in\cF^\star}\sup_{\P_{\star}^{(n)}\in\textup{subG}_n(f^\star,\varsigma^2)}\P_{\star}^{(n)}\sbr{\hXi_n^{\natural,\kappa}\del{\nnorm{\net(\btheta)- f^\star}\ge A_n\epsilon_n(\cF^\star)}}=\sco(1)
\end{equation}
for any diverging sequence $(A_n)_{n\in\bN}\to+\infty$, where $\epsilon_n(\cF^\star)$ is defined in \labelcref{eq:dnn_oracle}.
\end{corollary}

\section*{Acknowledgement}

We would like to thank Minwoo Chae, Kyoungjae Lee and Ryan Martin for their helpful comments and suggestions.  We acknowledge the generous support of NSF grants DMS CAREER 1654579 and DMS 2113642. The first author was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2021R1A6A3A03044911).

\bibliographystyle{plainnat}
\begin{thebibliography}{60}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\expandafter\ifx\csname urlstyle\endcsname\relax
\providecommand{\doi}[1]{doi: #1}\else
\providecommand{\doi}{doi: \begingroup \urlstyle{rm}\Url}\fi

[Alquier and Ridgway, 2020]
Pierre Alquier and James Ridgway.
\newblock Concentration of tempered posteriors and of their variational
\newblock \emph{The Annals of Statistics}, 48\penalty0 (3):\penalty0
1475--1497, 2020.

[Atchad{\'e}, 2017]
Yves~A Atchad{\'e}.
\newblock On the contraction properties of some high-dimensional
quasi-posterior distributions.
\newblock \emph{The Annals of Statistics}, 45\penalty0 (5):\penalty0
2248--2273, 2017.

[Atchad{\'e} and Bhattacharyya, 2018]
Yves~A Atchad{\'e} and Anwesha Bhattacharyya.
\newblock An approach to large-scale quasi-{Bayesian} inference with
spike-and-slab priors.
\newblock \emph{arXiv preprint arXiv:1803.10282}, 2018.

[Baraud and Birg{\'e}, 2020]
Yannick Baraud and Lucien Birg{\'e}.
\newblock Robust {Bayes}-like estimation: Rho-{Bayes} estimation.
\newblock \emph{The Annals of Statistics}, 48\penalty0 (6):\penalty0
3699--3720, 2020.

[Belitser and Ghosal, 2003]
Eduard Belitser and Subhashis Ghosal.
\newblock Adaptive {Bayesian} inference on the mean of an infinite-dimensional
normal distribution.
\newblock \emph{The Annals of Statistics}, 31\penalty0 (2):\penalty0 536--559,

[Bhattacharya et~al., 2019]
Anirban Bhattacharya, Debdeep Pati, and Yun Yang.
\newblock Bayesian fractional posteriors.
\newblock \emph{The Annals of Statistics}, 47\penalty0 (1):\penalty0 39--66,

[Bhattacharya and Martin, 2020]
Indrabati Bhattacharya and Ryan Martin.
\newblock Gibbs posterior inference on multivariate quantiles.
\newblock \emph{arXiv preprint arXiv:2002.01052}, 2020.

[Buja et~al., 1989]
Andreas Buja, Trevor Hastie, and Robert Tibshirani.
\newblock Linear smoothers and additive models.
\newblock \emph{The Annals of Statistics}, 17\penalty0 (2):\penalty0 453--510,

[Ch{\'e}rief-Abdellatif, 2020]
Badr-Eddine Ch{\'e}rief-Abdellatif.
\newblock Convergence rates of variational inference in sparse deep learning.
\newblock In \emph{Proceedings of the 37th International Conference on Machine
Learning}, pages 1831--1842. PMLR, 2020.

[Ch{\'e}rief-Abdellatif and Alquier, 2018]
Badr-Eddine Ch{\'e}rief-Abdellatif and Pierre Alquier.
\newblock Consistency of variational {Bayes} inference for estimation and model
selection in mixtures.
\newblock \emph{Electronic Journal of Statistics}, 12\penalty0 (2):\penalty0
2995--3035, 2018.

[Ch{\'e}rief-Abdellatif and Alquier, 2020]
Badr-Eddine Ch{\'e}rief-Abdellatif and Pierre Alquier.
\newblock {MMD-Bayes}: Robust {Bayesian} estimation via maximum mean
\newblock In \emph{Symposium on Advances in Approximate Bayesian Inference},
pages 1--21. PMLR, 2020.

[Draper, 1995]
David Draper.
\newblock Assessment and propagation of model uncertainty.
\newblock \emph{Journal of the Royal Statistical Society: Series B
(Methodological)}, 57\penalty0 (1):\penalty0 45--70, 1995.

[Ga{i}ffas and Lecu{\'e}, 2007]
St{\'e}phane Ga{i}ffas and Guillaume Lecu{\'e}.
\newblock Optimal rates and adaptation in the single-index model using
\newblock \emph{Electronic Journal of Statistics}, 1:\penalty0 538--573, 2007.

[Gao et~al., 2015]
Chao Gao, Yu~Lu, and Harrison~H Zhou.
\newblock Rate-optimal graphon estimation.
\newblock \emph{The Annals of Statistics}, 43\penalty0 (6):\penalty0
2624--2652, 2015.

[Gao et~al., 2020]
Chao Gao, Aad van~der Vaart, and Harrison~H Zhou.
\newblock A general framework for {Bayes} structured linear models.
\newblock \emph{The Annals of Statistics}, 48\penalty0 (5):\penalty0
2848--2878, 2020.

[Gershman et~al., 2012]
Samuel Gershman, Matthew~D Hoffman, and David~M Blei.
\newblock Nonparametric variational inference.
\newblock In \emph{Proceedings of the 29th International Conference on Machine
Learning}. PMLR, 2012.

[Ghosal et~al., 1997]
Subhashis Ghosal, Jayanta~K Ghosh, and RV~Ramamoorthi.
\newblock Non-informative priors via sieves and packing numbers.
\newblock In \emph{Advances in statistical decision theory and applications},
pages 119--132. Springer, 1997.

[Ghosal et~al., 2000]
Subhashis Ghosal, Jayanta~K Ghosh, and Aad van~der Vaart.
\newblock Convergence rates of posterior distributions.
\newblock \emph{The Annals of Statistics}, 28\penalty0 (2):\penalty0 500--531,

[Ghosal et~al., 2008]