A Targeted Accuracy Diagnostic for Variational Approximations

TL;DR

This paper introduces TADDAA, a diagnostic tool that uses multiple short MCMC chains to accurately assess the error in specific posterior functionals of variational approximations, improving evaluation in complex models.

Contribution

The paper presents TADDAA, a novel targeted diagnostic method that provides lower bounds on the error of posterior functionals using short MCMC chains, addressing limitations of existing diagnostics.

Findings

TADDAA effectively estimates errors in variational posterior functionals.

The method is computationally efficient and practical for complex models.

Numerical experiments demonstrate TADDAA's utility on real and synthetic data.

Abstract

Variational Inference (VI) is an attractive alternative to Markov Chain Monte Carlo (MCMC) due to its computational efficiency in the case of large datasets and/or complex models with high-dimensional parameters. However, evaluating the accuracy of variational approximations remains a challenge. Existing methods characterize the quality of the whole variational distribution, which is almost always poor in realistic applications, even if specific posterior functionals such as the component-wise means or variances are accurate. Hence, these diagnostics are of practical value only in limited circumstances. To address this issue, we propose the TArgeted Diagnostic for Distribution Approximation Accuracy (TADDAA), which uses many short parallel MCMC chains to obtain lower bounds on the error of each posterior functional of interest. We also develop a reliability check for TADDAA to determine when the lower bounds should not be trusted. Numerical experiments validate the practical utility and computational efficiency of our approach on a range of synthetic distributions and real-data examples, including sparse logistic regression and Bayesian neural network models.

Equations

$$\alpha (x,Y)=\min \left\{1,\frac{\pi (Y)q_{\psi }(Y,x)}{\pi (x)q_{\psi }(x,Y)}\right\}.$$

$$P_{\psi }(x,\mathrm{d}y)$$

$$\hat{\pi }=\underset{\xi \in \mathcal{Q}}{\mathrm{argmin}}{\mathcal{D}}_{\pi }(\xi ).$$

$${\mathcal{D}}_{\pi }(\xi )=\mathrm{KL}(\xi \mid \pi ):=\int \log \left(\frac{\mathrm{d}\xi }{\mathrm{d}\pi }\right)\mathrm{d}\xi$$

$$\begin{array}{rl}{\epsilon }^{(0)}(\mathcal{F})& \ge |\mathcal{F}({\hat{\pi }}^{(0)})-\mathcal{F}({\hat{\pi }}^{(T)})|\\ & \ge \{\begin{array}{ll}0& \text{if ℓF≤0≤uF}\\ \min (|{\ell }_{\mathcal{F}}|,|u_{\mathcal{F}}|)& \text{otherwise}\end{array}\\ & =\text{mathds}1\left\{0\notin ({\ell }_{\mathcal{F}},u_{\mathcal{F}})\right\}\times \min (|{\ell }_{\mathcal{F}}|,|u_{\mathcal{F}}|)\\ & =:B_{\mathcal{F}}.\end{array}$$

$${\overline{x}}_i^{(T)}-{\mu }_i^{(0)}\pm \frac{s_i^{(T)}}{\sqrt{N}\text{lower}0.4pt\text{vrule}height=6.83331pt,depth=-5.46667pt}{t}_{N-1}(1-\alpha /2),$$

$$\left[\log \left(\frac{(N-1)s_i^{(T)}}{({\sigma }_i^{(0)}{)}^2{\chi }_{N-1}^2(1-\alpha /2)}\right),\log \left(\frac{(N-1)s_i^{(T)}}{({\sigma }_i^{(0)}{)}^2{\chi }_{N-1}^2(\alpha /2)}\right)\right],$$

$${\mathcal{F}}_p^i(\pi )=Q_p^i(\pi ):=\inf \{q\in \mathbb{R}:P_{X\sim \pi }(X_i\le q)\ge p\},$$

$$\left[X_{(l),i}^{(T)}-Q_p^i({\hat{\pi }}^{(0)}),X_{(u),i}^{(T)}-Q_p^i({\hat{\pi }}^{(0)})\right],$$

$$H({\psi }^{(1:t)},X_{1:N}^{(1:t+1)},$$

$$={\psi }^{(t)}+\frac{1}{\sqrt{t+1}\text{lower}0.4pt\text{vrule}height=6.44444pt,depth=-5.15558pt}({\bar{\alpha }}^{(t)}-{\bar{\alpha }}_{\ast }).$$

$$\begin{array}{rl}\underset{N\to \infty }{\lim }{\mathbb{E}}_{X_{N,1:N}}\left\{\prod _{i=1}^r{g}_i\left(X_{N,i}\right)\right\}& =\prod _{i=1}^r\int g_i(x)\bar{\nu }(\mathrm{d}x).\end{array}$$

$$(X_j^{(t+1)},Y_j^{(t+1)})\sim T(X_j^{(t)},Y_j^{(t)},h_N^{(t)},\cdot ,\cdot )$$

$$\begin{array}{rl}N& =\max \left(N_{\text{mean}},N_{\text{variance}}\right),\end{array}$$

$$N_{\text{mean}}$$

$$N_{\text{variance}}$$

$$\displaystyle T=\bigl{\lfloor}{{{{cd^{1/3}}}}}\bigr{\rfloor},$$

$$\displaystyle T=\bigl{\lfloor}{{{{cd^{1/4}/L}}}}\bigr{\rfloor},$$

$$\begin{array}{r}{\rho }_i^2(T):=\frac{{\sum }_{j=1}^N(X_{j,i}^{(0)}-{\hat{\mu }}_i^{(0)})(X_{j,i}^{(T)}-{\hat{\mu }}_i^{(T)})}{\sqrt{{\sum }_{j=1}^N(X_{j,i}^{(0)}-{\hat{\mu }}_i^{(0)}{)}^2}\text{lower}0.4pt\text{vrule}height=8.27223pt,depth=-6.61781pt\sqrt{({\sum }_{j=1}^N{X}_{j,i}^{(T)}-{\hat{\mu }}_i^{(T)}{)}^2}\text{lower}0.4pt\text{vrule}height=8.27223pt,depth=-6.61781pt}.\end{array}$$

$$\begin{array}{rl}y& \mid \beta \sim \mathsf{Bern}({\mathrm{logit}}^{-1}(X\beta )),\\ {\beta }_j& \mid \tau ,\lambda ,c\sim \mathcal{N}(0,{\tau }^2{\tilde{\lambda }}_j^2),\\ {\lambda }_j& \sim {\mathrm{C}}^{+}(0,1),\tau \sim {\mathrm{C}}^{+}\left(0,{\tau }_0\right),\\ c^2& \sim \mathsf{InvGam}(2,8),\end{array}$$

$$y_n\sim \mathsf{Bern}\left({\mathrm{logit}}^{-1}\left(\tanh \left(\tanh ({\alpha }^{\top }{z}_n)\cdot \beta \right)\cdot \gamma \right)\right),$$

$$Q^{RW}(x,\mathrm{d}y)=\mathcal{N}(\mathrm{d}y\mid x,hG).$$

$$\displaystyle Q^{L}(x,\mathrm{d}y)=\mathcal{N}\left(\mathrm{d}y\biggm{|}x+\frac{h}{2}G\nabla\log\pi(x),hG\right).$$

$$Q_i^B(x,\mathrm{d}{y}_i)=2\frac{{\mu }_{{\tau }_i}(y_i-x_i)}{1+e^{-z_i{c}_i(x)}}\mathrm{d}{y}_i,$$

$$Q^B(x,\mathrm{d}y)=\prod _{i=1}^d{Q}_i^B\left(x,\mathrm{d}{y}_i\right).$$

$$\bar{\pi }(x,\eta )\propto \exp \left\{\log (\pi (x))-\frac{1}{2}{\eta }^{\top }G\eta \right\}.$$

$$\eta (t+h/2)$$

$$x(t+h)$$

$$\eta (t+h)$$

$$X_j^{(t+1)}\sim K\left(X_j^{(t)},h_N^{(t)},\cdot \right).$$

$$X_j^{(t+1)}=\{\begin{array}{ll}{Y}_j^{(t)},& \text{ with probability }{\alpha }_j^{(t)}\\ X_j^{(t)},& \text{ with probability }1-{\alpha }_j^{(t)},\end{array}$$

$$\begin{array}{rl}{\alpha }_j^{(t)}& :=g(X_j^{(t)},Y_j^{(t)},h_N^{(t)})\\ & :=\min \left\{1,\frac{\pi \left(Y_j^{(t)}\right)q_{h_N^{(t)}}\left(X_j^{(t)},Y_j^{(t)}\right)}{\pi \left(X_j^{(t)}\right)q_{h_N^{(t)}}\left(Y_j^{(t)},X_j^{(t)}\right)}\right\},\end{array}$$

Full text

1 INTRODUCTION

Bayesian inference is widely used due to its flexibility to handle arbitrary models, the desirable decision-theoretic properties of the provided point estimates, and its coherent uncertainty quantification (Robert et al.,, 2007). However, posterior point estimates and uncertainties usually cannot be calculated analytically, so in practice users of Bayesian methods rely on approximation algorithms. Variational inference (VI) has become an increasingly popular alternative to Markov chain Monte Carlo (MCMC) due to its computational efficiency and the availability of black-box versions that can be applied to nearly arbitrary models (Jordan et al.,, 1999; Wainwright et al.,, 2008; Hoffman et al.,, 2013; Kingma and Welling,, 2013; Blei et al.,, 2017; Kucukelbir et al.,, 2017; Agrawal et al.,, 2020).

The computational efficiency of variational inference arises from optimizing a divergence between the exact posterior and a distribution constrained to a tractable approximating family. The trade-off is that even the optimal distribution within this family may be a poor approximation to the posterior, so in general there is no guarantee that variational point estimates and uncertainties will be accurate. Therefore, it is important to develop diagnostic methods that quantify the error introduced by using a variational approximation. However, existing diagnostic methods either lack interpretability (Gorham and Mackey,, 2017, 2015; Wada and Fujisaki,, 2015; Chee and Toulis,, 2018; Burda et al.,, 2015; Domke and Sheldon,, 2018), or do not provide useful information when the parameter space is high-dimensional (Yao et al.,, 2018; Xing et al.,, 2020; Huggins et al.,, 2020). For example, Yao et al., (2018) propose using Pareto smoothed importance sampling (PSIS) $\hat{k}$, which essentially quantifies how effective the approximation is as a proposal distribution for importance sampling – not whether it will provide good approximations to specific quantities of interest. The kernel Stein discrepancy (KSD; Gorham and Mackey,, 2017, 2015) can be used to evaluate the discrepancy between the approximate and the true distribution but its value is not readily interpretable and depends on a distant dissipativity condition that is not always easy to verify in practice. VI quality can also be accessed by monitoring the running average of changes in the estimated objective, which is usually the evidence lower bound (ELBO) (Wada and Fujisaki,, 2015; Chee and Toulis,, 2018; Burda et al.,, 2015; Domke and Sheldon,, 2018). However, like the KSD, the ELBO lacks interpretability because it includes an unknown constant offset (the marginal likelihood), which is hard to estimate well unless the approximation is very accurate. Moreover, all these methods attempt to characterize the accuracy of the full variational distribution, which in practice is usually poor even if specific summaries such as posterior means, variances, or credible intervals are accurate. Therefore, in practice, they provide limited useful diagnostic information in the types of complex models where variational methods offer the greatest computational savings.

In this paper, we propose a diagnostic that is applicable to high-dimensional parameter spaces and can provide lower bounds on the error of specific posterior summaries, including marginal means and covariances. Our approach is to run many short MCMC chains with stationary distribution equal to the posterior $\pi$. These chains are initialized with independent samples from a posterior approximation $\hat{\pi}^{(0)}$, which, after $T$ iterations, results in an improved empirical distribution $\hat{\pi}^{(T)}$ that is closer to $\pi$. If $\hat{\pi}^{(0)}$ is a poor approximation to $\pi$, then the improved distribution $\hat{\pi}^{(T)}$ obtained using MCMC can be significantly different from $\hat{\pi}^{(0)}$; otherwise, the distribution will not shift. Thus, we can obtain lower bounds on the error of a posterior summary by computing confidence intervals for how much that summary differs under $\hat{\pi}^{(0)}$ and $\hat{\pi}^{(T)}$. Our approach might fail to provide meaningful results, however, if the Markov kernel used mixes poorly. To ameliorate this issue, we jointly adapt the Markov chains at each iteration. We show that if the number of chains is large, they remain essentially independent despite the joint adaptation procedure; hence, our confidence intervals remain valid. We also propose a simple correlation-based check that can detect poorly mixing chains, and hence warn the user against trusting the diagnostic. Numerical experiments validate the practical utility and computational efficiency of our approach on a range of synthetic distributions and real-data examples.111A Python implementation of TADDAA and code to reproduce all of our experiments is available at https://github.com/TARPS-group/TADDAA.

2 PRELIMINARIES

Modern Bayesian statistics relies heavily on models with difficult-to-compute posterior distributions, and hence requires general-purpose algorithms to approximate expectations with respect to those distributions. Throughout we abuse notation and use $\pi$ to denote the posterior distribution and the corresponding density.

Markov Chain Monte Carlo.

Markov chain Monte Carlo (MCMC) sampling methods provide a general-purpose framework for obtaining samples that are asymptotically exact (in runtime). To construct a Markov kernel with the desired stationary distribution, a standard approach is to correct an arbitrary proposal kernel using the Metropolis–Hastings (MH) accept-reject procedure (Metropolis et al.,, 1953). For $x\in\mathbb{R}^{d}$, let $Q_{\psi}(x,\mathrm{d}{y})$ denote the proposal kernel parameterized by $\psi$ with current state $x$ and corresponding density $q_{\psi}(x,y)$. A proposed state $Y\sim Q_{\psi}(x,\cdot)$ is accepted with probability

[TABLE]

The Metropolis–Hastings Markov transition kernel is therefore given by

[TABLE]

where $r_{\psi}(x):=1-\int\alpha(x,y)Q_{\psi}(x,\mathrm{d}{y})$ is the rejection probability and $\delta_{x}$ is a Dirac measure at $x$.

Variational Inference.

Variational inference (VI) provides a potentially faster alternative to MCMC when models are complex and/or the dataset size is large (Blei et al.,, 2017). Variational inference aims to minimize some measure of discrepancy $\mathcal{D}_{\pi}(\cdot)$ over a tractable family $\mathcal{Q}$ of potential approximating distributions (Wainwright et al.,, 2008), resulting in a posterior approximation given by

[TABLE]

We will focus on the Kullback–Leibler (KL) divergence

[TABLE]

as the discrepancy since it is the most widely used: the unknown marginal likelihood does not affect the optimization, computing gradients requires estimating expectations only with respect to distribution within a tractable family $\mathcal{Q}$, and optimization is relatively stable (Blei et al.,, 2017).

Alternative Non-exact Approximation Methods.

While our experiments focus on variational inference, numerous other non-exact posterior approximation methods are widely used in practice such as expectation-propagation (Minka,, 2001; Bishop,, 2006), the Laplace approximation (Gelman et al.,, 1995; Bishop,, 2006), and the iterated nested Laplace approximation (Rue et al.,, 2009; Gomez-Rubio,, 2020). Our diagnostic applies equally well to all of these methods.

3 AN MCMC-BASED ACCURACY DIAGNOSTIC

Suppose we have an approximation $\hat{\pi}^{(0)}$ to the posterior density $\pi$. Given a sample $X^{(0)}\sim\hat{\pi}^{(0)}$, we could improve it by applying a Markov kernel $P_{\psi}(x,\mathrm{d}{y})$ with invariant distribution $\pi$, resulting in a sample $X^{(1)}\sim P_{\psi}(X^{(0)},\cdot)$ with a distribution $\hat{\pi}^{(1)}$ that is typically closer to $\pi$ than $\hat{\pi}^{(0)}$ if $\hat{\pi}^{(0)}\neq\pi$. Continuing to run the Markov chain for $T$ total steps results in a sample $X^{(T)}$ with distribution $\hat{\pi}^{(T)}$ that, when $\hat{\pi}^{(0)}$ is a poor approximation to $\pi$, can be substantially closer to $\pi$; see Fig. 1 for an illustration. More specifically, consider some posterior functional of interest such as a mean or variance, which we generically denote by $\mathcal{F}$. For example, $\mathcal{F}(\pi)=\int x_{i}\pi(\mathrm{d}x)$ would be the mean of the $i$th component. If the initial approximation error $\varepsilon^{(0)}(\mathcal{F}):=\mathcal{F}(\hat{\pi}^{(0)})-\mathcal{F}(\pi)$ is large, then a well-chosen Markov kernel will result in the improved approximation error $\varepsilon^{(T)}(\mathcal{F}):=\mathcal{F}(\hat{\pi}^{(T)})-\mathcal{F}(\pi)$ being substantially smaller. In this case, $|\mathcal{F}(\hat{\pi}^{(0)})-\mathcal{F}(\hat{\pi}^{(T)})|$ measures the decrease in the initial approximation error, and hence provides a lower bound on $|\varepsilon^{(0)}(\mathcal{F})|$. If, on the other hand, the initial approximation error is small, then $|\mathcal{F}(\hat{\pi}^{(0)})-\mathcal{F}(\hat{\pi}^{(T)})|$ should remain small as well.

While in principle the just-described approach is appealing, there are two immediate problems. First, we do not have direct access to $\hat{\pi}^{(T)}$. Second, $P_{\psi}(x,\mathrm{d}{y})$ may not be particularly efficient unless the proposal parameter $\psi$ is well-chosen – which is nearly impossible to do a priori. To address both of these problems, we propose to run $N$ Markov chains in parallel and, following standard MCMC practice, adapt the proposal parameter at each iteration (Andrieu and Thoms,, 2008). Denote the state of the $j$th chain after iteration $t$ by $X_{j}^{(t)}$ and the proposal parameter used at iteration $t$ by $\psi^{(t-1)}$. We write $X_{1:N}^{(t)}=(X_{j}^{(t)})_{j=1}^{N}$ and similarly for $X_{j}^{(1:t)}$, $X_{1:N}^{(1:t)}$, and other quantities. In the case of a Metropolis–Hastings kernel, at iteration $t$, for each $j\in\{1,\dots,N\}$, we generate proposals $Y_{j}^{(t+1)}\sim Q_{\psi^{(t)}}(X_{j}^{(t)},\cdot)$, then accept or reject the proposal to obtain $X_{j}^{(t+1)}$. Next, we can use inter-chain adaptation (INCA; Craiu et al.,, 2009), which uses all the samples available up to the current time to update the proposal parameter to $\psi^{(t+1)}=H(\psi^{(1:t)},X_{1:N}^{(1:t+1)},Y_{1:N}^{(1:t+1)})$. We call $H$ the adaptation function.

Finally, we can use the $N$ samples from iteration $T$ to construct a $(1-\alpha)$-confidence interval $(\ell_{\mathcal{F}},u_{\mathcal{F}})$ for $\mathcal{F}(\hat{\pi}^{(0)})-\mathcal{F}(\hat{\pi}^{(T)})$, and hence a high-confidence lower bound on $|\mathcal{F}(\hat{\pi}^{(0)})-\mathcal{F}(\pi)|$, if the following condition holds:

Condition 1**.**

The map $t\mapsto|\varepsilon^{(t)}(\mathcal{F})|$ is non-increasing and the sign of $\varepsilon^{(t)}(\mathcal{F})$ does not depend on $t$.

Under 1, with probability $1-\alpha$,

[TABLE]

Algorithm 1 summarizes the key steps of our proposed TArgeted Diagnostic for Distribution Approximation Accuracy (TADDAA) for computing the error lower bound $B_{\mathcal{F}}$. References to specific implementation recommendations (all described and justified in detail below) are in magenta.

Remark 3.1**.**

While in general we expect 1 to hold, it is possible that it could be violated for target distributions with complex geometries. Although we only consider reversible Markov kernels, it is plausible that non-reversible Markov kernels could also lead to violations.

We will construct confidence intervals under the following condition.

Condition 2**.**

The elements of $X_{1:N}^{(T)}$ are identically distributed, have finite second moments, and are pairwise independent.

We will focus on diagnosing the accuracy of means, variances, and quantiles. But TADDAA can be straightforwardly applied to other functionals of interest. For $t\in\{0,\dots,T\}$, define the mean vector $\mu^{(t)}=(\mu_{1}^{(t)},\ldots,\mu_{d}^{(t)})^{\top}:=\mathbb{E}_{X\sim\hat{\pi}^{(t)}}(X)$, and marginal standard deviations $\sigma^{(t)}_{i}:=\operatorname{stdev}_{X\sim\hat{\pi}^{(t)}}(X_{i})$ for $i=1,\dots,d$.

Marginal mean.

For the mean functional $\mathcal{F}_{\text{mean}}^{i}(\pi):=\mu_{i}(\pi):=\int x_{i}\,\pi(\mathrm{d}x)$, under 2, a $(1-\alpha)$ confidence interval for $\mu_{i}^{(T)}-\mu_{i}^{(0)}$ is given by (Mendenhall et al.,, 2012)

[TABLE]

where $\overline{x}_{i}^{(T)}=\frac{1}{N}\sum_{j=1}^{N}x_{j,i}^{(T)}$ is the sample mean, s_{i}^{(T)}=\frac{1}{N-1}\big{\{}\sum_{j=1}^{N}(x_{j,i}^{(T)}-\overline{x}_{i}^{(T)})^{2}\big{\}}^{1/2} is the sample standard deviation, and $t_{N-1}(q)$ is the $q$ quantile of a $t$ distribution with $N-1$ degrees of freedom.

Marginal variance.

For the variance functional $\mathcal{F}_{\text{var}}^{i}(\pi):=\sigma_{i}^{2}:=\operatorname{Var}_{X\sim\pi}(X_{i})$, under 2, a $(1-\alpha)$ confidence interval for $2\log(\sigma_{i}^{(T)}/\sigma_{i}^{(0)})$ is given by (Mendenhall et al.,, 2012)

[TABLE]

where $\chi^{2}_{N-1}(q)$ is the $q$ quantile of a $\chi^{2}$ distribution with $N-1$ degrees of freedom.

Marginal quantile.

For the $p$-quantile functional

[TABLE]

under 2, a $(1-\alpha)$ confidence interval for $Q^{i}_{p}(\hat{\pi}^{(T)})-Q^{i}_{p}(\hat{\pi}^{(0)})$ is given by (Hahn and Meeker,, 2011)

[TABLE]

where $X^{(T)}_{(1),i}\leq X^{(T)}_{(2),i}\leq\cdots\leq X^{(T)}_{(N),i}$ are the order statistics of the samples $X_{1:N}^{(T)}$, $l=B(\alpha/2,N,p)$, $u=B(1-\alpha/2,N,p)+1$, and $B(q,N,p)$ is the $q$ quantile of a binomial distribution with parameters $N$ and $p$.

While we have now described the key components of TADDAA, a practical and reliable implementation requires addressing a number of issues:

Proposal kernel and adaptation scheme. The choice of MCMC proposal kernel and adaptation function $H$ can have a dramatic effect on sampling efficiency and thus the accuracy of TADDAA. We propose to use gradient-based proposals that scale well to high-dimensional problems and focus on adapting the step size parameter of these algorithms. (Section 3.1) 2. 2.

Inter-chain dependence. INCA introduces dependence between the Markov chains, which could invalidate the confidence intervals. However, we show that the samples from the Markov chains at any fixed iteration satisfy a propagation-of-chaos property, which implies that the chains are asymptotically independent (in the limit $N\to\infty$), so an asymptotic version of 2 still holds. (Section 3.2) 3. 3.

Choice of $T$ and $N$. The length of each Markov chain $T$ and the number of Markov chains $N$ must be chosen such that the diagnostic is computationally efficient while still (i) having $T$ large enough that the improved approximation $\hat{\pi}^{(T)}$ becomes sufficiently different from the initial approximation $\hat{\pi}^{(0)}$ and (ii) having $N$ large enough that the confidence intervals are sufficiently small to provide meaningful diagnostic information. (Section 3.3) 4. 4.

Checking the reliability of the diagnostic. Since it is possible the proposal kernel will be inefficient despite adaptation, we propose a correlation-based reliability check for our diagnostic. (Section 3.4)

3.1 Markov Kernels and Adaptation

We consider four possible Markov kernels: random walk Metropolis–Hastings (RWMH), Metropolis-adjusted Langevin algorithm (MALA), the Barker algorithm (Barker), and Hamiltonian Monte Carlo (HMC), which we describe in detail in Appendix A. All four kernels rely on a step-size parameter $h$. The latter three algorithms exploit gradient information about $\log\pi$, and thus have superior efficiency when $d$ is large. High-dimensional sampling efficiency can be formalized using the theory of optimal scaling (Roberts and Rosenthal,, 2001). To guarantee sampling efficiency and obtain well-behaved limiting behavior, optimal scaling results show that step size $h$ needs to scale as $\Theta(d^{-\gamma})$, which results in a mixing time (i.e., number of iterations to obtain an effectively independent sample) of $\Theta(d^{\gamma})$. Hence, the smaller $\gamma$ is the more efficient the kernel. The values of $\gamma$ for the kernels we consider are $1$ for RWMH (Gelman et al.,, 1997), $1/3$ for MALA (Roberts and Rosenthal,, 1998), $1/3$ for Barker (Livingstone and Zanella,, 2022), and $1/4$ for HMC (Beskos et al.,, 2013). Optimal scaling theory also provides the optimal asymptotic acceptance rates ($\bar{\alpha}_{*}$) for each kernel, which are 0.234 for RWMH, 0.4 for Barker, 0.574 for MALA, and 0.651 for HMC. While HMC has the best scaling with $\gamma=1/4$, the efficiency of Barker is more robust to the precise step size (and hence precise acceptance rate) (Livingstone and Zanella,, 2022), which makes it attractive for our use-case where it may not be possible to obtain the optimal step size since adaptation occurs over a relatively small number of iterations $T$.

We adapt the step size so that the acceptance rate approaches the optimal asymptotic acceptance rate $\bar{\alpha}_{*}$ (Andrieu and Thoms,, 2008). Let $\bar{\alpha}^{(t)}:=N^{-1}\sum_{j=1}^{N}\alpha^{(t)}_{j}$ denote the average acceptance probability at iteration $t$. Letting $\psi^{(t)}=\log h^{(t)}$ be the log of the step size, the adaptation function is given by

[TABLE]

Compared with running many parallel chains independent of each other, using a joint adaptation scheme improves adaption speed and ensures the samples at each iteration follow the same distribution (Rosenthal,, 2000; Solonen et al.,, 2012). The sequence of proposal parameters $(\psi^{(t)})_{t=0}^{\infty}$ converge under the standard Robbins–Monro conditions (see Andrieu and Thoms, (2008, Section 4) for a related discussion). In order to prove convergence to the stationary distribution of the resulting Markov chain, one needs to verify an additional assumption, such as uniform ergodicity Roberts and Rosenthal, (2007, Theorem 1) or one of the weaker sufficient conditions proposed, for instance in Atchadé and Rosenthal, (2005, Assumption 3.1) or Green et al., (2015, Section 2.4).

Following Roberts and Rosenthal, (1998), we initialize the step sizes to $h^{(0)}=2.4^{2}/d$ for RWMH, and $h^{(0)}={2.4^{2}}/{d^{1/3}}$ for MALA and Barker, and $h^{(0)}={2.4^{2}}/{d^{1/4}}$ for HMC. All of the Markov kernels we consider can make use of a pre-conditioning matrix, which is beneficial when components have different scales (Roy and Zhang,, 2022; Stramer and Roberts,, 2007). While it is possible to adapt the pre-conditioning matrix in addition to the step size in practice we found it sufficed to keep the pre-conditioning matrix fixed and equal to the covariance of $\hat{\pi}^{(0)}$. If the covariance is not available in closed form we use the sample covariance of the samples $X_{1:N}^{(0)}$.

3.2 Asymptotic Independence of Adapted Markov Chains

Due to the use of adaptation, the Markov chains $X_{1:N}^{(t)}$ are not independent once $t>1$, so the final samples $X_{1:N}^{(T)}$ violate the independence requirement of 2. However, we can instead ensure that the Markov chains are asymptotically independent in the following sense:

Definition 3.2**.**

Let $X_{N,1:N}=(X_{N,1},\dots,X_{N,N})$ denote a random vector. The sequence of random vectors $\{X_{N,1:N}\}_{N=1}^{\infty}$ is $\bar{\nu}$-chaotic if, for any $r\in\mathbb{N}$ and any bounded continuous real-valued functions $g_{1},g_{2},\ldots,g_{r}$,

[TABLE]

In particular, we show that the Markov chains are chaotic after any fixed number of iterations.

Theorem 3.3**.**

Under Assumption 3.4 (given below), for any $t\in\mathbb{N}$, there exists a probability distribution $\bar{\nu}^{(t)}$ such that the sequence $\{X_{1:N}^{(t)}\}_{N=1}^{\infty}$ is $\bar{\nu}^{(t)}$-chaotic.

The proof can be found in Appendix B. Crucially, the asymptotics we require for Theorem 3.3 are in the number of Markov chains $N$, not the number of iterations $T$. Thus, for any fixed $T$, 2 holds asymptotically in $N$.

To write our assumptions, note that if we treat the proposed states $Y_{1:N}^{(t)}$ as additional random variables, we can write the MCMC update as

[TABLE]

for an appropriate Markov kernel $T(x,y,h,\mathrm{d}x^{\prime},\mathrm{d}y^{\prime})$, where $h$ is the step size.

Assumption 3.4**.**

1. (a)

The proposal probability density $q_{h}(y,x)$ is continuous with respect to $(x,y,h)$. 2. (b)

The target distribution has a continuous probability density function $\pi(\cdot)$. 3. (c)

Samples generated from the Markov transition kernel $T(x,y,h,\cdot,\cdot)$ satisfy $\mathbb{E}\|X_{j}^{(t)}\|^{2}<\infty$ and $\mathbb{E}\|Y_{j}^{(t)}\|^{2}<\infty$ for any $t\in\mathbb{N}$.

These conditions are all quite mild. All of the proposals used in our experiments satisfy Assumption 3.4.(a). Assumption 3.4.(b) is required to use all three gradient-based proposals (MALA, Barker, and HMC). Assumption 3.4.(c) just requires that for any fixed time $t\in\mathbb{N}$, the generated samples have finite second moment.

Remark 3.5**.**

Theorem 3.3* could be generalized to targets with discrete components if the proposal distribution is absolutely continuous with respect to an appropriate reference measure and the topology of the space on which the target is defined respects the discrete and continuous components of the target. (For example, if the target distribution is defined on $\mathbb{R}$ and has an atom at zero, then $\{0\}$ should be an open set.) In that case, Definition 3.2 would need to be adapted in order to account for the change in the underlying topology, and so will the details of the proof of Theorem 3.3.*

3.3 Length and Number of Markov Chains

Having established the validity of Algorithm 1, we now turn to the question of how to ensure the computational cost is not prohibitive and the diagnostic is accurate. Given a target distribution and Markov kernel, the computational cost is determined by the number of Markov chains $N$ and iterations $T$. We address each in turn.

Number of Markov Chains.

The number of Markov chains must be sufficiently large that the confidence intervals are small enough to detect meaningful errors. Our approach is to choose $N$ to be the smallest value that satisfies the user’s tolerance level for the margin of error. In the case of mean estimation, it follows from Eq. 7 that the error for estimating $\mu_{i}^{(T)}$ is given by t_{N-1}(\alpha/2)\,{s}_{i}^{(T)}/\mathchoice{{\hbox{\displaystyle\sqrt{N,}}\lower 0.4pt\hbox{\vrule height=6.83331pt,depth=-5.46667pt}}}{{\hbox{\textstyle\sqrt{N,}}\lower 0.4pt\hbox{\vrule height=6.83331pt,depth=-5.46667pt}}}{{\hbox{\scriptstyle\sqrt{N,}}\lower 0.4pt\hbox{\vrule height=4.78333pt,depth=-3.82668pt}}}{{\hbox{\scriptscriptstyle\sqrt{N,}}\lower 0.4pt\hbox{\vrule height=3.41666pt,depth=-2.73334pt}}}. Similarly, the log variance error can also be determined by Eq. 8. For the mean we use the relative error normalized by the standard deviation, as this accounts for the relevant scale of the error. Hence, if the user’s tolerance is $\delta_{\text{mean}}$ for relative mean error and $\delta_{\text{var}}$ for log variance error222The width of the quantile confidence interval is sample dependent, so we do not set a tolerance level for the quantile error., then the required sample size (number of Markov chains) is

[TABLE]

where

[TABLE]

Number of Iterations.

The number of iterations has to be large enough that the improved distribution $\hat{\pi}^{(T)}$ is sufficiently different from a poor approximating distribution $\hat{\pi}^{(0)}$. As discussed earlier, the theory of optimal scaling shows that the Markov chain requires $\Theta(d^{\gamma})$ iterations to mix. Thus, we propose to scale $T$ with respect to dimension $d$ in the same manner. For RWMH333Based on optimal scaling behavior of RWMH, Eq. 18 may not be large enough to guarantee significant change in the distribution. So, we do not recommend RWMH., MALA, and Barker, we take

[TABLE]

while for HMC we take

[TABLE]

where $L$ is the number of leapfrog steps in HMC. Scaling by $L$ is necessary for HMC to ensure $T$ accounts for the $L$ gradient evaluations required per iteration. All other methods require 0 or 1 gradient evaluation per iteration. Based on numerical results we suggest $c=50$. In Section C.1, we also run ablation studies on how the diagnostic result (several types of lower bound) would change with $T$, which shows that Eq. 18 and Eq. 19 are reasonable.

The results of Bhatia et al., (2022) suggest our choice of $T$ will ensure the cost of TADDAA will not be too large compared to the cost of using black-box variational inference algorithms. Specifically, Bhatia et al., (2022, Theorem 1) shows that the number of gradient evaluations of stochastic optimization needed to obtain a constant optimization error scales as $d^{1/3}$. Hence, the computational cost (as measured by number of gradient evaluations) is $\Theta(d^{1/3})$ for VI, $\Theta(d^{1/3})$ for MALA and Barker, and $\Theta(d^{1/4})$ for HMC. We verify the cost of TADDAA is not prohibitive in practice in our numerical experiments.

Remark 3.6**.**

TADDAA is essentially running MCMC only for a burn-in period. Crucially, the diagnostic remains valid and useful even if the chains remain far from stationarity as long as $\hat{\pi}^{(T)}$ is closer to stationarity (i.e., it improves the approximating distribution $\hat{\pi}^{(0)}$ somewhat). In challenging scenarios where moving closer to stationarity is hard even with best practices (e.g., INCA and preconditioning), the TADDAA error lower bounds can be close to zero, leading to uninformative diagnostic results. See Section 4.3 for a concrete example. We note also that there is a large literature studying the convergence rates of MCMC algorithms (Brooks et al.,, 2011; Roberts and Rosenthal,, 2004).

3.4 A Reliability Check for the Diagnostic

The reliability of TADDAA depends on the mixing behavior of the Markov chains:

• •

If the Markov chains are mixing well, we expect the diagnostic results can be trusted.

• •

If the Markov chains are not mixing well, diagnosis of a poor approximation can still be trusted but diagnosis of a good approximation may not be reliable. In the latter case, we should consider increasing the length of Markov chains or otherwise improving the Markov kernel.

The usual approach to checking MCMC mixing is some version of the $\widehat{R}$ diagnostic (Gelman and Rubin,, 1992; Vats and Knudson,, 2021; Vehtari et al.,, 2021). However, using $\widehat{R}$ requires a sufficiently large effective sample size for each Markov chain, which is not satisfied for the short Markov chains used in TADDAA. Instead, we propose to leverage the fact that we have many nearly independent chains to compute the correlation between the initial and final values of each chain. If the chains are mixing, then this correlation will be close to zero. Specifically we propose to use the worst-case correlation coefficient $\rho^{2}_{\max}(T):=\max_{i}\rho^{2}_{i}(T),$ where

[TABLE]

If $\rho^{2}_{\max}(T)$ is close to 0 (we use a cutoff of $0.1$), then the check passes; otherwise, the check fails. In addition to examining the value of $\rho^{2}_{\max}(T)$, we can also examine plots of $t$ versus $\rho^{2}_{\max}(t)$ and $t$ versus $\rho^{2}_{i}(t)$ to gain insight into, how quickly mixing is occurring, respectively, overall and for specific coordinates.

4 EXPERIMENTS

We now validate the efficacy and computational efficiency of TADDAA and the correlation-based reliability check. In our experiments, the number of Markov chains is $N=386$, which corresponds to the tolerance level $\delta_{\text{mean}}=0.1$ and $\delta_{\text{var}}=0.15$. While we often present results for all four kernels discussed in Section 3.1 (RWMH, MALA, Barker, and HMC), we generally recommend using Barker because it has the efficiency of gradient-based kernels but is less sensitive to parameter tuning than MALA and HMC. For variational inference, we use the default settings of the bbvi function in the VIABEL package (Welandawe et al.,, 2022), which used a mean-field Gaussian approximation family. We set the maximum number of iterations to $200{,}000$ but the optimization can terminate early if convergence is reached. We compare TADDAA to the $\hat{k}$ diagnostic proposed in Yao et al., (2018) and the Wasserstein-based mean and variance upper bounds proposed in Huggins et al., (2020).

4.1 Gaussian Model With Correlated Coordinates

We first investigate the performance of TADDAA in a high-dimensional setting where the variational mean-field approximation has correct mean but incorrect variance estimates, which is a common situation in practice. Specifically, we use a correlated Gaussian target $\pi=\mathcal{N}(\mu,\Sigma)$, where $\mu\in\mathbb{R}^{d}$, $\Sigma_{ii}=\sigma_{i}^{2}$, and $\Sigma_{ij}=\rho\sigma_{i}\sigma_{j}$ for $i\neq j$. We introduce correlation by setting $\rho=0.7$ and variance heterogeneity by setting $\sigma_{1}^{2}=10$ and $\sigma_{i}^{2}=1$ for $i=2,3,...,d$.

Fig. 2 shows that the lower bounds $B_{\text{mean}}$ and $B_{\text{variance}}$ remain valid and useful even when $d$ and $\hat{k}$ are large: TADDAA correctly captures that the mean estimates are accurate but the variance estimates are inaccurate. From Fig. 3, we can see that the Barker chains mix well in all cases. Figure 2 also provides a comparison to $\hat{k}$ and the Wasserstein-based upper bounds proposed in Huggins et al., (2020, Theorem 3.4). The $\hat{k}$ result suggests that the VI approximation is not reliable for almost any dimension (since $\hat{k}>0.7$). But if the functional of interest is the mean, VI should be considered reliable. Similarly, the Wasserstein upper bounds are too conservative by a factor of 10–1,000, incorrectly suggesting the VI approximation is extremely inaccurate. Moreover, the upper bounds are unavailable once $\hat{k}>0.7$ (Huggins et al.,, 2020). Overall, TADDAA provides much more precise and actionable information about the approximation quality than the other diagnostics.

4.2 Neal-Funnel Shape Model

Next, we consider a more challenging target with a more complex geometry: Neal’s funnel distribution (Neal,, 2003), which is similar to the geometries encountered in hierarchical models. For $x=(x_{1},\dots,x_{d})\in\mathbb{R}^{d}$, let $\sigma:=\exp(x_{1})$. The funnel distribution is given by $\pi(x)=\mathcal{N}(x_{1}\mid 0,1)\prod_{i=2}^{d-1}\mathcal{N}(x_{i}\mid 0,\sigma)$. Hence, if $X\sim\pi$, then $\operatorname{Var}(X_{i})=\exp(1/2)$ for $i\in\{2,\dots,d\}$. Since $X_{2},\dots,X_{d}$ are identically distributed, we focus on the accuracy of the mean-field VI approximations to the distributions of $X_{1}$ and $X_{2}$.

Figure 4 shows that the lower bounds constructed using Barker are quite precise and do not degrade as the dimension increases. (See Section C.2 for results for quantiles.) Moreover the reliability check passes for all $d$ (Fig. 13 in Section C.2). The $\hat{k}$ diagnostic, however, is noisy and provides little insight. The Wasserstein-based mean and covariance upper bounds are orders-of-magnitude too large and unavailable for most $d>10$ because $\hat{k}>0.7$.

4.3 Logistic Regression Model on Candy Power Ranking Dataset

Next, we validate TADDAA on a real-world dataset444https://www.kaggle.com/datasets/fivethirtyeight/the-ultimate-halloween-candy-power-ranking: each observation $(z_{n},y_{n})$ represents a different kind of candy, where $z_{n}\in\mathbb{R}^{11}$ are the features, and $y_{n}=1$ if the candy is chocolate and 0 otherwise. We use a logistic regression model with $y_{n}\sim{\sf{Bern}}(\operatorname{logit}^{-1}(\beta^{\top}z_{n}))$ and $\beta_{i}\sim\mathcal{N}(0,I)$.

Focusing on Barker, Fig. 5 shows that for $\beta_{5}$, the lower bound of the mean error is small but non-zero due to the fact that the variance is dramatically underestimated, leading to inefficient preconditioning (adaptive preconditioning did not significantly improve the lower bound). On the other hand, the variance error lower bound, while not tight, shows that the variance estimate is orders-of-magnitude too small. For $\beta_{11}$, where the variance is not so poorly estimated, both mean and variance lower bounds are very accurate. TADDAA uses about 2% as many gradient evaluations as the VI optimization ($4.3\times 10^{4}$ versus $2\times 10^{6}$), validating the computational efficiency of our approach. Figure 6 demonstrates the usefulness of the reliability checks, which correctly show that HMC chains fail to mix (due to a smaller number of iterations), while the chains for other kernels mix well.

4.4 Cancer Classification Using a Horseshoe Prior

We now consider a more challenging, higher-dimensional application to predicting leukemia using microarray data with $n=71$ observations and using 100 features chosen according to their $\chi^{2}$ scores (Ray et al.,, 2020). We use a logistic regression model with a sparsity-inducing horseshoe prior (Piironen and Vehtari,, 2017):

[TABLE]

where $y$ denotes the binary outcomes, $X\in\mathbb{R}^{71\times 100}$ is the features matrix, $\tau>0$ and $\lambda>0$ are global and local shrinkage parameters, and $\tilde{\lambda}_{j}^{2}=c^{2}\lambda_{j}^{2}/(c^{2}+\tau^{2}\lambda_{j}^{2})$. Hence, for our data the parameter dimensionality is $d=203$. We would expect the VI approximation to be poor due to multimodality of the model. Focusing on two representative parameters, Fig. 7 shows that TADDAA using Barker correctly captures that the mean estimate for $\lambda_{61}$ and variance estimate for $\beta_{0}$ are inaccurate. On the other hand, for the accurate mean estimate for $\beta_{0}$ and variance estimate for $\lambda_{61}$, the lower bounds are zero or close to zero. TADDAA uses about 28% as many gradient evaluations as the VI optimization ($1.1\times 10^{5}$ versus $4\times 10^{5}$). In addition, the reliability checks shown in Fig. 8 correctly reflect the worse performance of RWMH (due to $d$ being larger) and HMC (due to a smaller number of iterations).

4.5 Bayesian Neural Network on Candy Power Ranking Dataset

Finally, we validate TADDAA on a Bayesian neural network model, which is difficult for mean-field VI to approximate accurately (Izmailov et al.,, 2021). We use the same dataset described in Section 4.3, but now it is fitted with a two-layer Bayesian neural network with 5 units in each hidden layer:

[TABLE]

where $\alpha\in\mathbb{R}^{5\times 11}$, $\beta\in\mathbb{R}^{5\times 5}$, $\gamma\in\mathbb{R}^{5}$, $\alpha_{i,j}\sim\mathcal{N}(0,2)$, $\beta_{i,j}\sim\mathcal{N}(0,2)$ and $\gamma_{i}\sim\mathcal{N}(0,2)$. Hence, $x=(\alpha,\beta,\gamma)$ and $d=85$. To quantify the classification quality for the Bayesian neural network, we can use the log loss (a.k.a., the cross-entropy loss) $\mathcal{L}(x)=-\frac{1}{N}\sum_{n=1}^{N}\log p(y_{n}\mid x).$ Similar to the definitions for the marginal mean and median, we can also use the mean and median of the log loss as the evaluation functional $\mathcal{F}$. Figure 9 shows that TADDAA correctly captures that the inaccuracy of the VI approximation: error lower bounds show that both the mean and median of the log loss are overestimated by mean-field VI. TADDAA uses about 12$\%$ as many gradient evaluations as the VI optimization ($8.4\times 10^{4}$ versus $6.8\times 10^{5}$). The reliability check in Fig. 15 confirms the superiority of Barker in high-dimensional problems.

5 CONCLUSION AND LIMITATIONS

Overall, our theory and experimental results show that TADDAA provides an efficient evaluation tool for variational inference, while also being applicable to other inexact posterior approximation methods like (integrated nested) Laplace approximations. TADDAA provides precise information about specific functionals of interest such as means and standard deviations, rather than a check on just the overall quality of an approximation – which in practice is often quite poor. However, there are a number of important limitations that users must keep in mind when employing TADDAA. While we have developed the correlation check to guard against a poor diagnostic, the check can fail due to, e.g., multimodality. On the other hand, such failures are not unique to our diagnostic and can affect, e.g., diagnostics for MCMC too. One approach to guard against this possibility is to run the VI optimization multiple times with diverse initializations, as is usually recommended for MCMC. Another limitation is that TADDAA is not appropriate in settings where it is only feasible to make a small number of passes over the entire dataset, so a method like stochastic variational inference (Hoffman et al.,, 2013) is used. Finally, it is possible that 1 could be violated. Developing checks against this failure mode is an important direction for future work.

Acknowledgments

M. Kasprzak was support by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 101024264. J. H. Huggins was partially supported by the National Institute of General Medical Sciences of the National Institutes of Health under grant number R01GM144963 as part of the Joint NSF/NIGMS Mathematical Biology Program. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Appendix A METROPOLIS-HASTING KERNELS

We briefly summarize the kernels used in our experiments. Throughout this section let $x\in\mathbb{R}^{d}$ denote the current state, $h\in\mathbb{R}_{+}$ the step size, and $G\in\mathbb{R}^{d\times d}$ a positive semi-definition preconditioning matrix.

Random Walk Metropolis

The (pre-conditioned) random walk Metropolis (RWMH) Metropolis et al., (1953) with and step size $h$ has proposal kernel

[TABLE]

Metropolis-adjusted Langevin Algorithm

In the (pre-conditioned) Metropolis-adjusted Langevin algorithm (MALA) (Roberts and Tweedie,, 1996; Stramer and Roberts,, 2007) has proposal kernel

[TABLE]

Barker Proposal

Let $C$ denote the Cholesky factor of $G$, set $\tau_{i}^{2}=hG_{ii}$, let $\mu(\cdot)_{\tau_{i}}$ be the probability density function of $\mathcal{N}(0,\tau_{i}^{2})$, and let $c_{i}(x)=(\nabla\log\pi(x)\cdot C^{T})_{i}$. For a proposal state $y\in\mathbb{R}^{d}$, let $z_{i}=((C^{\top})^{-1}(y-x))_{i}$. The (pre-conditioned) Barker proposal Livingstone and Zanella, (2022) for the $i$th coordinate is

[TABLE]

and the full Barker kernel $Q^{B}$ is

[TABLE]

Hamiltonian Monte Carlo

Hamiltonian Monte Carlo (HMC) (Neal,, 2011) is defined on an extended state space with the random momentum vector $\eta\overset{\textrm{{indep}}}{\sim}\mathcal{N}(0,G)$, yielding the joint density

[TABLE]

HMC updates a new state using the leapfrog integrator, which proceeds according to the updates

[TABLE]

Taking $x(0)=x$, the new proposal state is $y=x(Lh)$ with $\eta(0)\sim\mathcal{N}(0,G)$.

Appendix B PROOFS

For readability, we write the proposal kernel as $K(x,h,\mathrm{d}y)=K_{\log h}(x,\mathrm{d}y)$, so

[TABLE]

Or, equivalently,

[TABLE]

with $Y_{j}^{(t)}$ being the proposal state based on the proposal distribution $Q_{h_{N}^{(t)}}(X_{j}^{(t)},\cdot)$ and

[TABLE]

where $q_{h}(x,\cdot)$ is the probability density of $Q_{h}(x,\cdot)$ and $h_{N}^{(t)}$ is the step size at time $t$. If we treat the proposed states $Y_{1:N}^{(t)}$ as additional random variables, we can write

[TABLE]

for an appropriate Markov kernel $T(x,y,h,\mathrm{d}x^{\prime},\mathrm{d}y^{\prime})$. Now define the empirical distribution of $N$ particles $X_{j}^{(t)}$ at $t$th iteration as

[TABLE]

Similarly, we can define the empirical distribution of $N$ particles $(X_{j}^{(t)},Y_{j}^{(t)})$ at the $t$th iteration as

[TABLE]

Taking the initial step size $h_{N}^{(0)}$ as fixed, we can rewrite the $t$th step size $h_{N}^{(t)}$ update equation given in Eq. 11) as

[TABLE]

Letting $\bar{h}^{(0)}=h_{N}^{(0)}$, $\bar{\nu}^{(0)}=\hat{\pi}^{(0)}$ and $\bar{\xi}^{(0)}(\mathrm{d}x,\mathrm{d}y)=Q_{h_{N}^{(0)}}(x,\mathrm{d}y)\hat{\pi}^{(0)}(\mathrm{d}x)$, then we can recursively define $\bar{h}^{(t)}$, $\bar{\nu}^{(t)}$ and $\bar{\xi}^{(t)}$ by

[TABLE]

where

[TABLE]

Let $\mathcal{P}(S)$ denote the set of probability measures on the measurable space $(S,\mathcal{B})$. Let $\Xi_{N}^{(t)}$ and $\Pi_{N}^{(t)}$ denote the probability distributions of the respective empirical measures $\xi_{N}^{(t)}$ and $\nu_{N}^{(t)}$, so $\xi_{N}^{(t)}\in\mathcal{P}(\mathcal{P}(\mathbb{R}^{d}))$ and $\Pi_{N}^{(t)}\in\mathcal{P}(\mathcal{P}(\mathbb{R}^{d}\times\mathbb{R}^{d}))$ respectively.

We will make use of the following two results.

Theorem B.1**.**

(Sznitman,, 1991)*

$\{P_{N}^{(t)}\}$ is $\bar{\nu}^{(t)}$-chaotic if and only if*

[TABLE]

in $\mathcal{P}(\mathcal{P}(S))$.

Theorem B.2**.**

If $X,X_{n}:\Omega\to\mathbb{R}^{d}$ are random variables such that $X_{n}\xrightarrow{p}X$, then $f(X_{n})\stackrel{{\scriptstyle\mathrm{L}^{1}}}{{\longrightarrow}}f(X)$ for all $f\in\mathrm{C}_{b}(\mathbb{R}^{d})$.

B.1 Proof of Theorem 3.3

Theorem 3.3 will follow from a series of propositions and lemmas, which we prove in the subsequent subsections.

Proposition B.3**.**

Under Assumption 3.4, if $(X_{1:N}^{(t-1)},Y_{1:N}^{(t-1)})$ is $\bar{\xi}^{(t-1)}$-chaotic and $h_{N}^{(t-1)}\overset{p}{\to}\bar{h}^{(t-1)}$, then $(X_{1:N}^{(t)},Y_{1:N}^{(t)})$ is $\bar{\xi}^{(t)}$-chaotic.

Proposition B.4**.**

Under Assumptions 3.4.(a) and 3.4.(b), if $(X_{1:N}^{(t-1)},Y_{1:N}^{(t-1)})$ is $\bar{\xi}^{(t-1)}$-chaotic and $h_{N}^{(t-1)}\overset{p}{\to}\bar{h}^{(t-1)}$, then $h_{N}^{(t)}\overset{p}{\to}\bar{h}^{(t)}$.

Note that \big{\{}(X_{j}^{(0)},Y_{j}^{(0)})\big{\}}_{i=1}^{N} is $\bar{\xi}^{(0)}$-chaotic since the samples are i.i.d. and $h_{N}^{(0)}=\bar{h}^{(0)}$ is fixed. By induction, it follows from Proposition B.3 and Proposition B.4 that for any $t\in\mathbb{N}$, \big{\{}(X_{j}^{(t)},Y_{j}^{(t)})\big{\}}_{i=1}^{N} is $\bar{\xi}^{(t)}$-chaotic and $h_{k}^{N}\overset{p}{\to}\bar{h}^{(t)}$.

The following lemma implies that there is subsequence $\{\Pi_{N_{l}}^{(t)}\}_{l}$ that weakly converges.

Lemma B.5**.**

Under Assumption 3.4.(c), for any fixed $t\in\mathbb{N}$, $\{\Pi_{N}^{(t)}\}_{N}$ is tight and relatively compact.

Lemma B.6**.**

Under Assumption 3.4, if $X_{1:N}^{(t-1)}$ is $\bar{\nu}_{t-1}$-chaotic, $h_{N}^{(t-1)}\overset{p}{\to}\bar{h}^{(t-1)}$, and $\{\Pi_{N_{l}}^{(t)}\}_{l}$ is a convergent subsequence of $\{\Pi_{N}^{(t)}\}_{N}$ with limit $\Pi^{(t)}$, then $\Pi^{(t)}$ is a Dirac measure concentrated on $\bar{\nu}^{(t)}$.

Lemma B.7**.**

Under assumptions and conditions of Lemma B.6, $X_{1:N}^{(t)}$ is $\bar{\nu}^{(t)}$-chaotic.

We also know that $X_{1:N}^{(0)}$ is $\bar{\nu}^{(0)}$-chaotic (due to the fact that the $X_{1:N}^{(0)}$ are i.i.d.), and $h_{k}^{N}\overset{p}{\to}\bar{h}^{(t)}$ for any $t\in\mathbb{N}$. Hence, by recursively using Lemma B.7, we conclude that $X_{1:N}^{(t)}$ is $\bar{\nu}^{(t)}$-chaotic for any $t\in\mathbb{N}$.

B.2 Proof of Proposition B.3

To prove Proposition B.3, we will first show that the probability distribution $\{\Xi_{N}^{(t)}\}_{N}$ of empirical measure $\{\xi_{N}^{(t)}\}_{N}$ is tight and relatively compact, which guarantees that there exists convergent subsequence $\{\Xi_{N_{l}}^{(t)}\}_{l}$. Then we will prove that every convergent subsequence $\{\Xi_{N_{l}}^{(t)}\}_{l}$ has the same limit $\Xi^{(t)}$ that is a Dirac measure concentrated on $\bar{\xi}^{(t)}$.

Lemma B.8**.**

Under Assumption 3.4.(c), for any fixed $t\in\mathbb{N}$, $\{\Xi_{N}^{(t)}\}_{N}$ is tight and relatively compact.

Lemma B.9**.**

Under Assumption 3.4, suppose $(X_{1:N}^{(t-1)},Y_{1:N}^{(t-1)})$ is $\bar{\xi}^{(t-1)}$-chaotic and $h_{N}^{(t-1)}\overset{p}{\to}\bar{h}^{(t-1)}$. Let $\{\Xi_{N_{l}}^{(t)}\}_{l}$ be a convergent subsequence of $\{\Xi_{N}^{(t)}\}_{N}$ with limit $\Xi^{(t)}$. Then $\Xi^{(t)}$ is a Dirac measure concentrated on $\bar{\xi}^{(t)}$.

The proofs of Lemma B.8 and B.9 are deferred to Section B.2.1 and B.2.2. Relative compactness of $\{\Xi_{N}^{(t)}\}_{N}$ implies that there is subsequence $\{\Xi_{N_{l}}^{(t)}\}_{l}$ which weakly converges. In this proof, with a slight abuse of notation, we use $\{\Xi_{N}^{(t)}\}_{N}$ to represent the convergent subsequence $\{\Xi_{N_{l}}^{(t)}\}_{l}$. By Lemma B.8 and Lemma B.9, we have

[TABLE]

Theorem B.1 guarantees that $(X_{1:N}^{(t)},Y_{1:N}^{(t)})$ is $\bar{\xi}^{(t)}$-chaotic.

B.2.1 Proof of Lemma B.8

The proof is similar to that of Lemma 2.2 in Sirignano and Spiliopoulos, (2020). For a fixed $t\in\mathbb{N}$, $\forall\epsilon>0$, we want to show that there is a compact subset $K$ of $\mathcal{P}(\mathbb{R}^{d}\times\mathbb{R}^{d})$ such that

[TABLE]

For each $L>0$, define $K_{L}=[-L,L]^{2d}$, and by Assumption 3.4.(c), there exists $C$ such that

[TABLE]

For $L>0$, define the compact set

[TABLE]

Observe that

[TABLE]

so $\lim_{L\to\infty}\mathbb{P}[\xi_{N}^{(t)}\notin\hat{K}_{L}]=0$. Thus, the sequence $(\xi_{N}^{(t)})_{N\in\mathbb{N}}$ is tight, which implies it is relatively compact.

B.2.2 Proof of Lemma B.9

For each $f\in C_{b}(\mathbb{R}^{d}\times\mathbb{R}^{d})$, define the map $F:\mathcal{P}(\mathbb{R}^{d}\times\mathbb{R}^{d})\rightarrow\mathbb{R}_{+}$ by

[TABLE]

We would like to show that $\lim_{N\to\infty}\mathbb{E}_{\xi\sim\Xi_{N}^{(t)}}[F(\xi)]=0$. Together with the fact that subsequence $\{\Xi_{N}^{(t)}\}_{N}$ converges in distribution to $\Xi^{(t)}$, we would have $\mathbb{E}_{\xi\sim\Xi^{(t)}}[F(\xi)]=0$. Then we can conclude that $\Xi^{(t)}$ is a Dirac measure concentrated on $\bar{\xi}^{(t)}$.

Letting $\langle f,\tilde{\xi}_{N}^{(t)}\rangle=\frac{1}{N}\sum_{j=1}^{N}\langle f,T(X_{j}^{(t-1)},Y_{j}^{(t-1)},\bar{h}^{(t-1)},\cdot,\cdot)\rangle$, we have

[TABLE]

We will prove that $\lim_{N\to\infty}\mathbb{E}_{\xi\sim\Xi_{N}^{(t)}}\left[F(\xi)\right]=0$ by the following two lemmas:

Lemma B.10**.**

Under the assumptions of Lemma B.9, we have

[TABLE]

Lemma B.11**.**

Under the assumptions of Lemma B.9, we have

[TABLE]

By Lemma B.10 and B.11, we have

[TABLE]

Since $F$ is uniformly bounded, we have

[TABLE]

Since this holds true for any $f\in C_{b}(\mathbb{R}^{d}\times\mathbb{R}^{d})$, using the fact that $F$ is nonnegative and by Theorem 5.1 from Billingsley, (2013), we have that $\Xi_{N}^{(t)}$ has a limit point that is a Dirac measure concentrated on $\bar{\xi}^{(t)}$.

B.2.3 Proof of Lemma B.10

Let $\mathcal{F}_{N}^{(t)}$ be the $\sigma$-algebra generated by $X_{1:N}^{(1:m)}$ and $Y_{1:N}^{(1:m)}$, for $i=1,2,...,N$ and $m=1,2,...,t$.

Note that

[TABLE]

Since $f\in C_{b}(\mathbb{R}^{d}\times\mathbb{R}^{d})$, there exists $M$ such that for any $i=1,2,\ldots N$, we have

[TABLE]

It follows that

[TABLE]

The last inequality holds true because $(X_{j}^{(t)},Y_{j}^{(t)})$ are independent given $\mathcal{F}_{N}^{(t-1)}$. Therefore, we have

[TABLE]

And \mathbb{E}\left[\bigg{|}\frac{1}{N}\sum_{i=1}^{N}f\left(X_{j}^{(t)},Y_{j}^{(t)}\right)-\frac{1}{N}\sum_{i=1}^{N}\langle f,T\left(X_{j}^{(t-1)},Y_{j}^{(t-1)},h_{N}^{(t-1)},\cdot,\cdot\right)\rangle\bigg{|}\mid\mathcal{F}_{N}^{(t-1)}\right] is bounded due to the fact that $f\in C_{b}(\mathbb{R}^{d}\times\mathbb{R}^{d})$, then by dominated convergence theorem, we have

[TABLE]

Note that

[TABLE]

Since we know $(X_{1}^{(t-1)},Y_{1}^{(t-1)},h_{N}^{(t-1)})\overset{p}{\to}(X_{1}^{(t-1)},Y_{1}^{(t-1)},\bar{h}^{(t-1)})$ and $\langle f,T(x,y,h,\cdot,\cdot)\rangle$ is bounded and continuous with respect to $(x,y,h)$, thanks to Theorem B.2, we have

[TABLE]

B.2.4 Proof of Lemma B.11

Note that

[TABLE]

The above expectation is taken with respect to $(X_{j}^{(t-1)},Y_{i,t-1})$. For any fixed $(X_{j}^{(t-1)},Y_{i,t-1})$

[TABLE]

By Assumption 3.4.(a) and 3.4.(b), it is bounded and continuous with respect to $(X_{j}^{(t-1)},Y_{j}^{(t-1)})$. Then thanks to that fact that $(X_{1:N}^{(t-1)},Y_{1:N}^{(t-1)})$ is $\bar{\xi}^{(t-1)}$-chaotic, we have

[TABLE]

Since $f\in C_{b}(\mathbb{R}^{d}\times\mathbb{R}^{d})$, we have \mathbb{E}\left[\bigg{|}\langle f,T\left(X_{j}^{(t-1)},Y_{j}^{(t-1)},\bar{h}^{(t-1)},\cdot,\cdot\right)\rangle-\left\langle f,\bar{\xi}^{(t)}\right\rangle\bigg{|}^{2}\right]\ <\infty. Therefore, we have

[TABLE]

B.3 Proof of Proposition B.4

Recall that the step size update is

[TABLE]

Note that

[TABLE]

The first term in Eq. 66 can be rewritten as

[TABLE]

Since $(X_{j}^{(t-1)},Y_{j}^{(t-1)},h_{N}^{(t-1)})\overset{p}{\to}(X_{j}^{(t-1)},Y_{j}^{(t-1)},\bar{h}^{(t-1)})$ and $g$ is bounded and continuous, it follows from Theorem B.2 that

[TABLE]

The second term of Eq. 66 can be rewritten as

[TABLE]

Since $(X_{1:N}^{(t-1)},Y_{1:N}^{(t-1)})$ is $\bar{\xi}^{(t-1)}$-chaotic, by Theorem B.1 we have

[TABLE]

By Assumptions 3.4.(a) and 3.4.(b), we have $g(x,y,\bar{h}^{(t-1)})$ is continuous with respect to $(x,y)$, thus

[TABLE]

It follows from Eq. 68 and Eq. 71 that

[TABLE]

and therefore $h_{N}^{(t)}\overset{p}{\to}\bar{h}^{(t)}$.

B.4 Proofs of Lemmas B.5, B.6 and B.7.

The proof of Lemma B.5 is similar to that of Lemma B.8; the proof of Lemma B.6 is similar to that of Lemma B.9; the proof of Lemma B.7 is similar to that of Proposition B.3.

Appendix C ADDITIONAL RESULTS AND EXPERIMENTS

C.1 Ablation Study On $T$

In Section 3.3, we discussed how to choose length of Markov chains properly: Eq. 18 and Eq. 19 guarantee that we can obtain efficient samples under different dimension $d$ without too large computation cost. Now, we would like to explore how the lower bounds would change under different number of iterations $T$.

Figures 10, 11 and 12 display the evolution of several lower bounds for $X_{1}$ of Gaussian model (Section 4.1), $X_{1}$ of Neal-funnel shape model (Section 4.2) and $\lambda_{61}$ of cancer classification using horseshoe prior (Section 4.4) respectively, which shows that when $\rho$ is large (greater than 0.1), the proposed lower bounds are quite underestimated. These results also confirm the lowers bounds become nearly stationary at our proposed number of iterations $T$.

C.2 Additional Results for Neal-Funnel Shape Model

In Section 4.2, we validate TADDAA on Neal’s funnel distribution approximated by mean-field Gaussian VI. Fig. 13 provides reliability checks for the Neal’s funnel experiment Fig. 4.

Next, we will validate the TADDAA on Neal’s funnel distribution approximated by mean-field t-distribution VI. The lower bounds $B_{\text{mean}}$ and $B_{\text{tail}}$ are computed according to Eq. 6 using functionals of interest $\mathcal{F}_{0.5}$ and $\mathcal{F}_{0.9}$ defined based on Eq. 9.

Fig. 14 shows that the median and tail quantile lower bounds constructed using Barker are quite precise and remain valid when $d$ and $hat{k}$ are large. As we can see, $\hat{k}$ is large ($\hat{k}>0.7$) for almost any dimension, indicating VI approximation is not reliable. But if we are only interested in whether median of $X_{2}$, mean field t-distribution VI should be considered as reliable.

C.3 Reliability Check for Bayesian Neural Network

In Section 4.5, we validate TADDAA on Bayesian neural network approximated by mean-field Gaussian VI. The reliability check is shown in Fig. 15: due to the fact that the variance is dramatically underestimated, only Barker chains pass the reliability check, which is also consistent with the fact that Barker kernel provides the best lower bounds as displayed in Fig. 9.