Challenges and Opportunities in High-dimensional Variational Inference
Year: 2022

Abstract:
Current black-box variational inference (BBVI) methods require the user to make numerous design choices—such as the selection of variational objective and approximating family—yet there is little principled guidance on how to do so. We develop a conceptual framework and set of experimental tools to understand the effects of these choices, which we leverage to propose best practices for maximizing posterior approximation accuracy. Our approach is based on studying the pre-asymptotic tail behavior of the density ratios between the joint distribution and the variational approximation, then exploiting insights and tools from the importance sampling literature. Our framework and supporting experiments help to distinguish between the behavior of BBVI methods for approximating low-dimensional versus moderate-to-high-dimensional posteriors. In the latter case, we show that mass-covering variational objectives are difficult to optimize and do not improve accuracy, but flexible variational families can improve accuracy and the effectiveness of importance sampling—at the cost of additional optimization challenges. Therefore, for moderate-to-high-dimensional posteriors we recommend using the (mode-seeking) exclusive KL divergence since it is the easiest to optimize, and improving the variational family or using model parameter transformations to make the posterior and optimal variational approximation more similar. On the other hand, in low-dimensional settings, we show that heavy-tailed variational families and mass-covering divergences are effective and can increase the chances that the approximation can be improved by importance sampling.

heading: 1 Introduction
text: A great deal of progress has been made in black-box variational inference (BBVI) methods for Bayesian posterior approximation, but the interplay between the approximating family, divergence measure, gradient estimators and stochastic optimizer is non-trivial – and even more so for highdimensional posteriors [1, 10, 29, 31]. While the main focus in the machine learning literature has
∗ Equal contribution.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
been on improving predictive accuracy, the choice of method components becomes even more critical when the goal is to obtain accurate summaries of the posterior itself.
In this paper, we show that, while the choice of approximating family and divergence is often motivated by low-dimensional illustrations, the intuition from these examples do not necessarily carry over to higher-dimensional settings. By drawing a connection between importance sampling and the estimation of common divergences used in BBVI, we are able to develop a comprehensive framework for understanding the reliability of BBVI in terms of the pre-asymptotic behavior of the density ratio between the target and the approximate distribution. When this density ratio is heavy-tailed, even unbiased estimators exhibit a large bias with high probability, in addition to high variance. Such heavy tails occur when there is a mismatch between the typical sets of the approximating and target distributions. In higher dimensions, even over-dispersed distributions miss the typical set of the target [18, 27]. Thus, as illustrated in Fig. 1, the benefits of heavy-tailed approximate families and divergences favoring mass-covering diminish as dimensionality of the target distribution increases. Building on these insights, we make the following main contributions:
1. We develop a conceptual and experimental framework for predicting and empirically evaluating the reliability of BBVI based on the choice of variational objective, approximating family, and target distribution. Our framework also incorporates the Pareto k diagnostic [27] as a simple and practical approach for obtaining empirical and conceptual insights into the pre-asymptotic convergence rates of estimators of common divergences and their gradients.
2. We validate our framework through an extensive empirical study using simulated data and many commonly used real datasets with both Gaussian and non-Gaussian target distributions. We consider the exclusive and inclusive Kullback-Leibler (KL) divergences [4, 21, 24], tail-adaptive f -divergence [29], χ2 divergence [7], and α-divergences [12], and the resulting variational approximation for isotropic Gaussian and Student-t and normalising flow families.
3. Based on our framework and numerical results, we provide justified recommendations on design choices for different scenarios, including low- to moderate-dimensional and high-dimensional posteriors.

heading: 2 Preliminaries and Background
text: Let p(θ, Y ) be a joint distribution of a probabilistic model, where θ ∈ RD is a vector of model parameters and Y is the observed data. In Bayesian analysis, the posterior p(θ) := p(θ | Y ) = p(θ, Y )/p(Y ) (where p(Y ) := ∫ p(θ, Y )dθ ) is typically the object of interest, but most posterior summaries of interest are not accessible because the normalizing integral, in general, is intractable. Variational inference approximates the exact posterior p(θ | Y ) using a distribution q ∈ Q from a family of tractable distributionsQ. The best approximation is determined by minimizing a divergence D(p ‖ q), which measures the discrepancy between p and q:
qλ∗ = arg min qλ∈Q
D(p ‖ q), (1)
where λ ∈ RK is a vector parameterizing the variational family Q. Thus, the properties of the resulting approximation q are determined by the choice of variational family Q as well as the choice of divergence D.
The family Q is often chosen such that quantities of interest (e.g., moments of q) can be computed efficiently. For example, q can be used to compute Monte Carlo or importance sampling estimates of the quantities of interest. Let w(θ) := p(θ, Y )/q(θ) denote the density ratio between the joint and approximate distributions. For a function φ : RD → R, the biased self-normalized importance sampling estimator for the posterior expectation Eθ∼p[φ(θ)] is given by
Î(φ) := S∑ s=1 w(θs)∑S s′=1 w(θs′) φ(θs),
where θ1, . . . , θS ∼ q are independent. Using importance sampling can allow for computation of more accurate posterior summaries and to go beyond the limitations of the variational family. For example, it is possible to estimate the posterior covariance even when using a mean-field variational family.
Pareto Smoothed Importance Sampling. Since importance sampling estimates can have very high variance, Pareto smoothed importance sampling (PSIS) can be used to substantially reduce the variance with small additional bias [27]. This procedure modifies and stabilises extreme importance ratios using a generalized Pareto distribution fit to the upper tail of the distribution of the ratios.
Variational families. Let qλ(θ) be an approximating family parameterised by a K-dimensional vector λ ∈ RK for D-dimensional inputs θ ∈ RD. Typical choices of q include mean-field Gaussian and Student’s t families [3, 14], full and low rank Gaussians [15, 22], mixtures of exponential families [17, 19], and normalising flows [25]. We focus on the most popular mean-field and normalizing flow families. Mean-field families assume independence across the D dimensions: q(θ) = ∏D i=1 qi(θi), where each qi typically belongs to some exponential family or other simple class of distributions. Normalising flows [1] provide more flexible families that can capture correlation and non-linear dependencies. A normalizing flow is defined via the transformation of a probability density through a sequence of invertible mappings. By composing several maps, a simple distribution such as a mean-field Gaussian can be transformed into a more complex distribution [25].
f -divergences. The most commonly used divergences are examples of f -divergences [28]. For a convex function f satisfying f(1) = 0, the f -divergence is given by
Df (p ‖ q) := Eθ∼q [ f ( p(θ | Y ) q(θ) )] .
The exclusive Kullback-Leibler (KL) divergence corresponds to f(w) = − log(w), the inclusive KL divergence corresponds to f(w) = w log(w), the χ2 divergence corresponds to f(w) = (w − 1)2, and the general α-divergences correspond to (wα − w)/{α(α− 1)}. We also consider the adaptive f -divergence proposed by Wang et al. [29].
Loss estimation and stochastic optimization. In all the cases we consider, minimizing the f - divergence is equivalent to minimizing the loss function Lf (p ‖ q) := Eθ∼q[f(w(θ))] (although, see Wan et al. [28] for a different approach). Let L(λ) := Lf (p ‖ qλ) denote the loss as a function of the variational parameters λ. The loss and its gradient G(λ) := ∇λL(λ) can both be approximated using, respectively, the Monte Carlo estimates
L̂(λ) = 1S ∑S s=1 f(w(θs)) and Ĝ(λ) = 1 S ∑S s=1 g(θs), (2)
where θ1, . . . , θS are independent draws from qλ and g : RK → RK is an appropriate gradient-like function that depends on f and w. The two most popular gradient estimators in the literature are the score function and the reparameterization gradient estimator [20, 30]. The score function gradient corresponds to g(θ) = {f(w(θ)) − w(θ)f ′(w(θ))}∇λ log qλ(θ). It is a general-purpose estimator that applies to both discrete and continuous distributions q, but it is known to suffer from high variance. When this estimator is used for the mass-covering divergences such as the inclusive KL and general α-divergences with α > 1, the importance weights are usually replaced with selfnormalized importance weights w(θs)/ ∑S i=1 w(θi). The reparameterization gradient [20] requires expressing the distribution qλ as a deterministic transformation of a simpler base distribution r such that Tλ(z) ∼ qλ with z ∼ r. This allows writing an expectation with respect to qλ as an expectation over the simpler distribution r. The reparameterization estimator corresponds to using g(zs) = ∇λf(w(Tλ(zs))) (for zs ∼ r) in place of g(θ), where w implicitly depends on λ as well. In the case of the adaptive f -divergence, the importance weights w(θ1), . . . , w(θS) are sorted, and the gradients corresponding to each sample are then weighed by the empirical rank. The gradient estimates can be used in a stochastic gradient optimization scheme such that
λt+1 ← λt + ηtĜ(λt), (3) where ηt is the step size. In practice, more stable adaptive stochastic gradient optimisation methods such as RMSProp or Adam [9, 13], which smooth or normalize the noisy gradients, are often used.
Numerous prior work have studied some of the challenges tied to optimizing these divergence measures under the presence of noisy gradient estimates [1, 10, 29, 31]. Particularly, when dealing with mass-covering divergences, the gradient estimates can become so noisy that convergence is not possible in practice, as we will illustrate later on.

heading: 3 Assessing the Reliability of Black-box Variational Inference

heading: 3.1 Conceptual framework
text: How can we determine – both conceptually and experimentally – what is required to obtain reliable estimates of the variational divergence and optimal variational approximation? As we have seen, the most common variational divergences and their Monte Carlo gradient estimators can be expressed in terms of the density ratio w(θ). Reliable black-box variational inference ultimately depends on the behavior of w(θ) since (1) accurate optimization requires low-variance and (nearly) unbiased gradient estimates Ĝ(λ), and (2) determining convergence and validating the quality of variational approximations can require accurate estimates L̂(λ) of variational divergences [14, 15]. While asymptotically (in the number of iterations and Monte Carlo sample size S) there may be no issues with stochastic optimization or divergence estimation, in practice black-box variational inference operates in the pre-asymptotic regime. Therefore, the reliability of black-box variational inference depends on the pre-asymptotic behavior of the w(θ), and how it interacts with the choice of variational objective and gradient estimator.
Before accounting for the effects of the objective and gradient estimator, first consider the behavior of the density ratio w(θ), which can also be interpreted as an importance sampling weight with qλ(θ) as the proposal distribution [cf. 2, 16, 29]. Pickands [23] proved, under commonly satisfied conditions, that for u tending to infinity, the distribution of w(θ) | w(θ) > u is well-approximated by the three-parameter generalized Pareto distribution GPD(u, σ, k), which for k > 0 has density p(w | u, σ, k) = σ−1{1 + k(w − u)/σ}−1−1/k where w is restricted to (u,∞). Since w(θ) > 0, this implies its distribution is heavily skewed to the right with a power-law tail. Consider the idealized scenario of estimating the mean of w(θ) ∼ GPD(u, σ, k). We assume the mean is finite, which is equivalent to assuming k < 1 since b1/kc determines the number of finite moments. Because of the heavy right skew, most of the mass of w(θ) is below its mean. Therefore, even after averaging a large number of samples, most empirical estimates ∑S s=1 w(θs) will be smaller than the true mean. Figure 3a illustrates this behavior for different values of k: even with 1 million samples, the empirical mean is far below the true mean when k > 0.7. The highly variable sizes of the confidence intervals based on 10,000 replications further highlight the instability of the estimator. So, even though the empirical mean is an unbiased estimator, in the pre-asymptotic regime (before the generalized central limit theorem is applicable [5]), in practice the estimates are heavily biased downward with high probability. If w(θ) is not a generalized Pareto distribution, we can instead treat k as
the tail index k := inf{` > 0 : Eθ∼q{w(θ)1/`} < ∞}, which encodes the same tail behavior as GPD(u, σ, k). Crucially, we should expect k to be much larger than 0 when there is a significant mismatch between the target distribution and the variational family. Since selecting a variational family that can match the typical set tends to be more difficult in higher dimensions, we should expect k to be larger for higher-dimensional posteriors.
We can generalize our observations about pre-asymptotic estimation bias to the estimators L̂(λ) and Ĝ(λ). For the loss estimator, we replace w(θs) with f(w(θs)), where f(w) is polynomial in w and logw for the class of losses we consider. If the dominant term of f(w) is of order wα, the tail behavior will be similar to a generalized Pareto with kα = αk. Thus, L̂(λ) will have larger pre-asymptotic bias as α increases. For example, estimation of the mass-covering inclusive KL (where α = 1) – and, more generally, mass-covering α-divergences with α > 0 – will suffer from a large pre-asymptotic bias. On the other hand, for the mode-seeking exclusive KL, f(w) = log(w), so we can expect all moments to be finite and therefore a much smaller pre-asymptotic bias.
Similar considerations apply to the gradient estimator, with the details depending on the specific estimator used. However, when using self-normalized weights for α-divergences, we can expect a large pre-asymptotic bias whenever w(θ) has such bias since self-normalization involves estimating the mean ofw(θ). This bias will affect the accuracy of the solution found using stochastic optimization. Thus, the quality of the solutions found can only partially be improved by using a smaller step size since smaller step sizes will only reduce the effects of a large estimator variance, but not the effects from a large bias. We provide more details on the behavior of the score function and reparameterized gradients for each of the divergences in ????, following Geffner and Domke [11].
In summary, our framework makes two key predictions:
(P1) Estimates and gradients of mode-seeking divergences (in particular exclusive KL divergence with log dependence on w) have lower variance and are less biased than those of mass-covering divergences (in particular α-divergences with α > 0, with polynomial dependence on w).
(P2) The degree of polynomial dependence on w determines how rapidly the bias and variance will increase as approximation accuracy degrades – in particular, in high dimensions.
Because the adaptive f -divergence depends directly on the (ordered) weights, we expect it to behave similarly to the mass-covering divergences.

heading: 3.2 Experimental framework
text: In the light of potentially large non-asymptotic bias arising from the heavy right tail of w(θ), it is important to verify the pre-asymptotic behavior of the Monte Carlo estimators used in variational inference. We follow the approach developed by Vehtari et al. [27] for importance sampling and compute an empirical estimate k̂ of the tail index k by fitting a generalized Pareto distribution to the observed tail draws. In the importance sampling setting, Vehtari et al. [27] show that the minimal sample size to have a small error with high probability scales as S = O(exp{k/(1− k)2}). Vehtari
et al. [27] also demonstrate that k̂ provides a practical pre-asymptotic convergence rate estimate even when the variance is infinite and a generalized central limit theorem holds. While estimating k̂ in general requires larger sample size than is commonly used to estimate the stochastic gradients, we can still use it to diagnose and identify the challenges with different divergences. If k̂ > 0.7, the minimal sample size to obtain a reliable Monte Carlo estimate is so large that it is usually infeasible in practice. This cutoff is in agreement with our findings shown in Fig. 3a. Thus, together with our conceptual framework, we have a third key prediction:
(P3) The k̂ value can be used to diagnose pre-asymptotic reliability of variational objectives. In particular, the α-divergence with α > 0 will become unreliable when max(1, α) × k̂ > 0.7, even if w is bounded (by a very large constant).

heading: 3.3 Verification of Pre-asymptotic (Un)reliability
text: We first verify our three key predictions in a simple setting where we can compute most of the relevant quantities such as the loss function in closed form. Specifically, we fit a mean-field Gaussian to a Gaussian with constant 0.5 correlation factor using the inclusive KL, exclusive KL, χ2, and 1/2-divergences. We vary the dimensionality D from 1 to 50, which is a surrogate for the degree of mismatch between the optimal variational approximation and the target distribution. To find the optimal divergence-based approximation, we optimize the closed-form expression for the divergences between two Gaussians. Hence, we can consider on the best-case scenario and ignore the complexities and uncertainty due to the stochastic optimization. Due to space limitations, we focus on representative cases of the approximations from optimising the mode-seeking exclusive KL divergence and the mass-covering inclusive KL divergence. Results for the other divergences are included in the appendix.
(P1) Mode-seeking divergences are more stable and reliable than mass-covering ones. Figure 3b shows that as the approximation–target mismatch increases with dimension, the bias in and variance of the divergence estimates increases substantially for the inclusive KL and χ2 but only moderately for the exclusive KL. Similarly, Fig. 2 shows that gradient bias and variance increases with dimension for inclusive KL and χ2 but not exclusive KL.
(P2) Degree of polynomial dependence on w determines sensitivity to approximation–target mismatch. Figure 3b shows that divergence estimates resulting from optimising higher polynomials of w become more and more unstable as dimensions increases.
(P3) k̂ diagnoses pre-asymptotic reliability. Figure 3c shows that the k̂ values grow rapidly for the inclusive KL-based approximation, particularly for higher-degree dependence on w, which agrees with predicted behavior and large bias and variance of the inclusive KL and χ2. In contrast, the k̂ values remain fairly stable for the exclusive KL-based approximation, again in agreement with predicted and observed bias and variance behavior.

heading: 4 Experiments
text: In this section, we describe a series of experiments to study how our pre-asymptotic framework can be used for assessing the reliability of black-box variational approximations for practical applications and developing best-practices. For all posteriors, we fit mean-field Gaussian and Student-t families, a planar flow [25] with 6 layers and a non-volume preserving (NVP) flow [8] with 6 stacked neural networks with 2 hidden layers of 10 neurons each for both the translation and scaling operations with a standard Gaussian distribution for the latent variables. We use Stan [26] for model construction. For stochastic optimization we use RMSProp with initial step size of 10−3 run for either Tmax iterations or until convergence was detected using a modified version of the algorithm by Dhaka et al. [6]. For the exclusive KL we use 10 draws for gradient estimation per iteration, while for the other divergences we use 200 draws, and a warm start at the solution of the exclusive KL. In practice, we found the optimisation for χ2 divergence extremely challenging, with the solution failing to converge even for moderate dimensions D ≈ 10. Therefore, we only include results for the KL divergences and the adaptive f -divergence. We compare the accuracy of approximated posterior moments to ground-truth computed either analytically or using the dynamic Hamiltonian Monte Carlo algorithm in Stan [26]. Specifically, we consider the estimates µ̂ and Σ̂ for, respectively, the posterior mean µ and covariance matrix Σ. We also consider the mean and covariance estimates produced by PSIS and compute k̂. The experiments were carried on a laptop and an internal cluster with only CPU capability. The code for the experiments will be made available after acceptance using MIT license.

heading: 4.1 Heavy-tailed posteriors
text: First, we study the toy robust regression model previously used by Huggins et al. [14] given by
βd ∼ N(0, 10), yn | xn, β ∼ t10(β>xn, 1),
where yn ∈ R, xn ∈ RD are the target and predictors respectively, β denotes the unknown coefficients, and D is varied from 2 to 50. We generated data from the same model with covariates generated from a zero-mean Gaussian with constant correlation of 0.4. The Student’s t leads to the posterior having heavy tails, making it a more challenging target distribution. We use Tmax = 10,000.
Mode-seeking divergences are easier to optimize. Figure 4a shows that the estimated tail index k̂ generally increases with the dimension as expected. In particular, the k̂ values when using normalizing flows, which are more challenging to optimize, is low for D < 20 when using exclusive KL, but infinite when using either the inclusive KL or f -divergence. From Fig. 4b we can see that exclusive KL provides also more accurate and reliable posterior approximations than the inclusive KL and adaptive f -divergence, particularly for the normalizing flows. This observation is consistent with the prediction (P3) of the proposed framework. The better performance for normalizing flows corroborates the relative ease of stochastic optimization with the exclusive KL divergence compared to the inclusive KL or the adaptive f -divergence – despite the fact that we used 20 times as many Monte Carlo samples to estimate the gradients for the inclusive KL and the f -divergence compared to the exclusive KL. To further illustrate the relative difficulty of optimizing the inclusive KL divergence, Fig. 5 shows the largest dimension for which the stochastic optimization converged as a function of the step-size. For most step-sizes, the combination of normalizing flows and the inclusive KL divergence only converged for D = 2, whereas convergence is possible in higher dimensions for simpler variational families. These observations are consistent with predictions (P1)-(P2) of the proposed framework.
Adaptive f -divergence interpolates between the exclusive and inclusive KL divergence, but is difficult to optimize. In low dimensions, the adaptive f -divergence behaves somewhere between the two KL divergences as seen in ?? – as it was designed to [29]. As confirmed by Fig. 4, For higher-dimensional posteriors, we expect it to behave more like the exclusive KL, but it is less stable due to its functional dependence on the importance weights.
Normalizing flows can be effective but are challenging to optimize. Fig. 4 also shows that normalizing flows can be quite effective when used with exclusive KL to ensure stable optimization. However, as can be seen in ??, when using out-of-the-box optimization with no problem-specific tuning (as we have done for a fair comparison), the normalizing flows approximations can have pathological features – even in low dimensions.

heading: 4.2 Realistic models and datasets
text: We now study how the choice of divergence and approximating family compare across a diverse range of benchmark posteriors. We compare variational approximations for models and datasets from posteriordb ∗ in terms of accuracy of the estimated moments and predictive likelihood. We used an 80/20 training/test split on all datasets to compute the predictive likelihoods. We use Tmax = 15,000.
∗ https://github.com/stan-dev/posteriordb
Exclusive KL remains the most reliable for realistic posteriors. The results are summarized in Fig. 6, where the same pattern is seen: the exclusive KL is superior for higher-dimensional posteriors (e.g., D > 10) or when combined with normalizing flows, while inclusive KL is better for lower-dimensional posteriors. Despite the superior performance of the exclusive KL divergence, the large values for k̂ indicate that fitting approximations based on normalizing flows remains a challenge in high dimensions. The performance for the adaptive f -divergence is comparable to the inclusive KL divergence. Table 1 shows that the exclusive KL divergence consistently outperforms the inclusive KL divergence in terms of predictive accuracy, but can be significantly worse than HMC.
Importance sampling can substantially improve accuracy. Focusing on exclusive KL, Fig. 6b shows the relative errors of the first two moments for the variational approximation (dots) and after correcting the estimates using PSIS (triangles). In some cases, the PSIS correction dramatically improved the accuracy of the normalizing flows.
Reparameterization is an important tool for improving accuracy. The 8-schools model is lowdimensional (D = 10), but the funnel-shaped posterior makes inference challenging for variational approximations [14, 31]. As has been noted previously in the literature, and is clear from Figs. 6a and 6b, reparameterizing the model so that the posterior better matches the variational family can be an effective way to improve the accuracy of the approximation. See ?? for an illustration.

heading: 5 Discussion
text: Our conceptual framework based on the pre-asymptotic behavior of the density ratios / importance weights w along with our comprehensive experiments lead to a number of important takeaways for practitioners looking to obtain reasonably accurate posterior approximations using black-box variational inference:
• The instability of mass-covering divergences like inclusive KL and χ2 means that, given currently available methodology, users are better off using the exclusive KL divergence except for easy low-dimensional posteriors. The reliance of the adaptive f -divergence on importance weights leads to similar instability.
• Importance sampling appears to almost always be beneficial for improving accuracy, even when the k̂ diagnostic is large. However, a large k̂ does suggest the user should not expect even the PSIS-corrected estimates to be particularly accurate.
• Using normalizing flows – particularly NVP flows – together with exclusive KL and PSIS provides the best and most consistent performance across posteriors of varying dimensionality and difficulty. We therefore suggest this combination as a good default choice.
Our results suggest an important direction for future work is improving the stability of optimization with normalizing flows, which still tend to have some pathological behaviors unless they are very carefully tuned since such tuning significantly detracts from the benefits of using BBVI.

heading: 6 Limitations
text: While our experiments included a range of common statistical model types, our findings may not generalize to all types of posteriors or to other variational families. For example, we did not explore semi-implicit methods or applications to neural networks. We also did not investigate alternative divergences such as those used in importance-weighted autoencoders.

heading: Acknowledgments
text: We thank Sihan Liu for pointing out and fixing errors in the code. We thank Ari Hartikainen and PyStan developers for help with posteriordb, Academy of Finland (grants 298742, and 313122) and Finnish Center for Artificial Intelligence for partial support of the research. We also acknowledge the computational resources provided by the Aalto Science-IT project. AKD would like to thank Tomas Geffner and Eero Siivola for useful discussions.

Review summaries:
The authors propose a conceptual framework and design guidelines for variational inference based on the pre-asymptotic tail behaviour of the unnormalized importance weight. This framework is used to study the behaviour of commonly variational families and f-divergences in mid-to-high dimensions. The authors observe that mass-covering divergences, while often superior theoretically and for low dimensional problems, are hard to optimize in mid-to-high dimensions compared to the mode-seeking exclusive KL-divergence. These observations are consistent with the prediction made using their framework based on estimating the tail-index of the best-fit generalized pareto distribution.

The paper discusses the accuracy of variational inference and the choices that need to be me made with regard to divergence and approximating family, particularly as they relate to the dimensionality of the problem. Based on an analysis of the pre-asymptotic behavior of the density ratios involved, the authors propose a framework based on Pareto smoothed importance sampling (PSIS) to help practitioners analyze the success of their approximate inference procedure. The authors provide experiments demonstrating general findings for variational inference for a range of divergences and approximate posteriors. 

In this paper, the authors study the pre-asymptotic behavior of the density ratios between the joint distribution and the variational approximation. The authors use generalized Pareto distribution to model this density ratio in the pre-asymptotic regime and make three major predictions for how BBVI methods should behave. They verify their predictions through extensive experimentation, and finally, provide actionable advice to BBVI practitioners (use flows with exclusive KL and PSIS as a good default.) 

This paper provides a statistical tool to analyze different variational inference methods in the f-divergence family. It showed that for high dimensional and low dimensional problems, the behavior of the variational inference methods may differ, and the analysis based on the GPD statistics can suggestion different strategies for different situations. The predictions from this theory were validated by some commonly used models (e.g., linear regression). In summary, this paper can provide useful guidelines for variational inference practitioners.