Beyond ELBOs: A Large-Scale Evaluation of Variational Methods for Sampling

Sampling methods are designed to address the challenge of generating approximate samples or estimating the intractable normalization constant $Z$ for a probability density $\pi$ on $\mathbb{R}^{d}$ of the form

where $\gamma:\mathbb{R}^{d}\rightarrow\mathbb{R}^{+}$ can be pointwise evaluated. This formulation has broad applications in fields such as Bayesian statistics and the natural sciences (Liu & Liu, 2001; Stoltz et al., 2010; Frenkel & Smit, 2023; Mittal et al., 2023).

Monte Carlo (MC) methods (Hammersley, 2013), including Annealed Importance Sampling (AIS) (Neal, 2001) and its Sequential Monte Carlo (SMC) extensions (Del Moral et al., 2006), have traditionally been considered the gold standard for addressing the sampling problem. An alternative approach is Variational Inference (VI) (Blei et al., 2017), where a tractable family of distributions is parameterized, and optimization tools are employed to maximize similarity to the intractable target distribution $\pi$.

In recent years, there has been a surge of interest in the development of sampling methods that merge MC with VI techniques to approximate complex, potentially multimodal distributions (Wu et al., 2020a; Zhang & Chen, 2021; Arbel et al., 2021; Matthews et al., 2022; Jankowiak & Phan, 2022; Midgley et al., 2022; Berner et al., 2022; Richter et al., 2023; Vargas et al., 2023a, b; Akhound-Sadegh et al., 2024).

However, the evaluation of these methods faces significant challenges, including the absence of a standardized set of tasks and diverse performance criteria. This diversity complicates meaningful comparisons between methods. Existing evaluation protocols, such as the evidence lower bound (ELBO), often rely on samples from the model, restricting their evaluation capabilities to the model’s support. This limitation becomes especially problematic when assessing the ability to mitigate mode collapse on target densities with well-separated modes. To overcome this challenge, others propose the use of integral probability metrics (IPMs), like maximum mean discrepancy (Arenz et al., 2018) or Wasserstein distance (Richter et al., 2023; Vargas et al., 2023a), leveraging samples from the target density to assess performance beyond the model’s support. However, these metrics often involve subjective design choices such as kernel selection or cost function determination, potentially leading to biased evaluation protocols.

To address these challenges, our work introduces a comprehensive set of tasks for evaluating variational methods for sampling. We explore existing evaluation criteria and propose a novel metric specifically tailored to quantify mode collapse. Through this evaluation, we aim to provide valuable insights into the strengths and weaknesses of current sampling methods, contributing to the future design of more effective techniques and the establishment of standardized evaluation protocols.

We provide an overview of Monte Carlo methods, Variational Inference, and combinations. The notation introduced in this section is used throughout the remainder of this work.

Monte Carlo Methods. A variety of Monte Carlo (MC) techniques have been developed to tackle the sampling problem and estimation of $Z$. In particular sequential importance sampling methods such as Annealed Importance Sampling (AIS) (Neal, 2001) and its Sequential Monte Carlo (SMC) extensions (Del Moral et al., 2006) are often regarded as a gold standard to compute $Z$. These approaches construct a sequence of distributions $(\pi_{t})_{t=1}^{T}$ that ‘anneal’ smoothly from a tractable proposal distribution $\pi_{0}$ to the target distribution $\pi_{T}=\pi$. One typically uses the geometric average, that is, $\gamma_{t}(\mathbf{x})=\pi_{0}(\mathbf{x})^{1-\beta_{t}}\gamma(\mathbf{x})^{\beta_{t}}$, with $\pi_{t}\propto\gamma_{t}$ for $0=\beta_{0}<\beta_{1}<...<\beta_{T}=1$. Approximate samples from $\pi$ are then obtained by starting from $\mathbf{x}_{0}\sim\pi_{0}(\cdot)$ and running a sequence of Markov chain Monte Carlo (MCMC) transitions that target $(\pi_{t})_{t=1}^{T}$.

Variational Inference. Variational inference (VI) (Blei et al., 2017) is a popular alternative to MCMC and SMC where one considers a flexible family of easy-to-sample distributions $q^{\theta}$ whose parameters $\theta$ are optimized by minimizing the reverse Kullback–Leibler (KL) divergence, i.e.,

It is well known that minimizing the reverse KL is equivalent to maximizing the evidence lower bound (ELBO) and that $\text{ELBO}\leq\log Z$ with equality if and only if $q^{\theta}=\pi$. Later, VI was extended to other variational objectives such as $\alpha$-divergences (Li & Turner, 2016; Midgley et al., 2022), log-variance loss (Richter et al., 2020), trajectory balance, (Malkin et al., 2022a) or general $f$- divergences (Wan et al., 2020). Typical choices for $q^{\theta}$ include mean-field approximations (Wainwright & Jordan, 2008), mixture models (Arenz et al., 2022) or normalizing flows (Papamakarios et al., 2021). To construct more flexible variational distributions (Agakov & Barber, 2004) modeled $q^{\theta}(\mathbf{x})$ as the marginal of a latent variable model, i.e. $q^{\theta}(\mathbf{x})=\int q^{\theta}(\mathbf{x},\mathbf{z})\mathrm{d}\mathbf{z}$ ¹¹1Agakov & Barber (2004) coined the term ‘augmentation’ for $\mathbf{z}$. We adopt the more established terminology and refer to $\mathbf{z}$ as a latent variable.. As this marginal is typically intractable, $\theta$ is then learned by minimizing a discrepancy measure between $q^{\theta}(\mathbf{x},\mathbf{z})$ and an extended target $p^{\theta}(\mathbf{x},\mathbf{z})=\pi(\mathbf{x})p^{\theta}(\mathbf{z}|\mathbf{x})$ where $p^{\theta}(\mathbf{z}|\mathbf{x})$ is an auxiliary conditional distribution. Using the chain rule for the KL-divergence (Cover, 1999) one obtains an extended version of the ELBO, that is,

Although the extended ELBO is often referred to as ELBO, latent variables $\mathbf{z}$ introduce additional looseness, i.e., $\,\overline{\!{\text{ELBO}}}\leq\text{ELBO}$ with equality if $p^{\theta}(\mathbf{z}|\mathbf{x})=q^{\theta}(\mathbf{x},\mathbf{z})/q^{\theta}(\mathbf{x})$. To compute expectations with respect to $q^{\theta}(\mathbf{x},\mathbf{z})$, one typically chooses tractable distributions $q^{\theta}(\mathbf{x}|\mathbf{z})$ and $q^{\theta}(\mathbf{z})$ and performs a Monte Carlo estimate using ancestral sampling.

Variational Monte Carlo Methods. Over recent years, the idea of using extended distributions has been further explored (Wu et al., 2020b; Geffner & Domke, 2021; Thin et al., 2021; Zhang et al., 2021; Doucet et al., 2022b; Geffner & Domke, 2022). In particular, these ideas marry Monte Carlo with variational techniques by constructing the variational distribution and extended target as Markov chains, i.e., $q^{\theta}(\mathbf{x}_{0:T})=\pi_{0}(\mathbf{x}_{0})\prod_{t=1}^{T}F^{\theta}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})$ and $p^{\theta}(\mathbf{x}_{0:T})=\pi(\mathbf{x}_{T})\prod_{t=0}^{T-1}B^{\theta}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t+1})$ with $\mathbf{x}=\mathbf{x}_{T}$, $\mathbf{z}=(\mathbf{x}_{0},\ldots,\mathbf{x}_{T-1})$ and tractable $\pi_{0}$. Common choices of transition kernels $F^{\theta}_{t},B^{\theta}_{t}$ include Gaussian distributions (Doucet et al., 2022b; Geffner & Domke, 2022) or normalizing flow maps (Wu et al., 2020a; Arbel et al., 2021; Matthews et al., 2022) and can be optimized by e.g. maximization of the extended ELBO via stochastic gradient ascent. Recently, Vargas et al. (2023a); Zhang & Chen (2021); Vargas et al. (2023b, 2024); Richter et al. (2023); Berner et al. (2022) explored the limit of $T\rightarrow\infty$ in which case the Markov chains converge to forward and backward time stochastic differential equations (SDEs) (Anderson, 1982; Song et al., 2020) inducing the path distributions $\mathbb{Q}^{\theta}$ and $\mathbb{P}^{\theta}$ which can be thought of as continuous time analogous of $q^{\theta}$ and $p^{\theta}$ respectively. Zhang & Chen (2021); Berner et al. (2022); Richter et al. (2023); Vargas et al. (2024) leveraged the continuous-time perspective to establish connections with Schrödinger bridges (Léonard, 2013) and stochastic optimal control (Dai Pra, 1991), resulting in the development of novel sampling algorithms.

Performance Criteria. Several performance criteria have been proposed for evaluating sampling methods, notably, those comparing the density ratio between the target and model density and integral probability metrics (IPMs).

Density ratio-based criteria make use of the ratio between the (unnormalized) target density $\gamma(\mathbf{x})$ and the model $q^{\theta}(\mathbf{x})$. Due to the intractability of $q^{\theta}(\mathbf{x})$ for methods that work with latent variables, the density ratio between the joint distributions of $\mathbf{x}$ and $\mathbf{z}$ is considered, i.e.,

respectively. Note that $w$ and $\,\overline{\!{w}}$ are also referred to as (unnormalized) importance weights. Using this notation, we can recover commonly used metrics such as the reverse effective sample size (ESS_r) or the ELBO, that is,

respectively. Here, ‘reverse’ is used to denote that expecations are computed with respect to $q^{\theta}$. In addition, if the true normalization constant is known, an importance-weighted reverse estimate of $\log Z$ is often employed to report the esimation bias, i.e., $\Delta\log Z_{r}=|\log Z-\log\hat{Z}_{r}|$ with

Please note that extended versions of these criteria are obtainable by replacing $w$ with the extended version $\,\overline{\!{w}}$ and taking expectations under the joint distribution $q^{\theta}(\mathbf{x},\mathbf{z})$.

Quantifying the ability to avoid mode collapse is difficult as identifying all modes of the target density $\pi$ and determining whether a model captures them accurately is inherently challenging. In particular, methods that are optimized using the reverse KL divergence are forced to assign high probability to regions with non-negligible probability in the target distribution $\pi$. This is referred to as mode-seeking behavior and can result in an overemphasis on a limited set of modes, leading to mode collapse. Consequently, performance criteria that use expectations under the model $q^{\theta}$, such as ELBO, (reverse) ESS, or $\Delta\log Z_{r}$, are influenced by the mode-seeking nature of the reverse KL divergence, making them less sensitive to mode collapse.

Here, we aim to explore criteria that are sensitive to mode collapse such as density-ratio based ‘forward’ criteria, that leverage expectations under $\pi$ and integral probability metrics (IPMs). Furthermore, we introduce entropic mode coverage, a novel criterion that leverages prior knowledge about the target to heuristically quantify mode coverage.

Forward Criteria. We discuss the ‘forward’ versions of the criteria discussed in Section 2. First, evidence upper bounds (EUBOs) are based on the forward KL divergence and have already been leveraged as learning objectives in VI (Ji & Shen, 2019). Here, we explore them as performance criteria that are sensitive to mode collapse. Formally, the EUBO is the sum of the forward KL and $\log Z$, that is,

with importance weights $w={\gamma(\mathbf{x})}/{q^{\theta}(\mathbf{x})}$. Due to the non-negativeness of the KL divergence, it is easy to see that $\text{EUBO}\geq\log Z$ with equality if and only if $q^{\theta}=\pi$. Hence, a lower EUBO means that $q^{\theta}$ is closer to $\pi$ in a $D_{\text{KL}}$ sense. The mass-covering nature of the forward KL leads to high EUBO values if the model fails to cover regions of non-negligible probability in the target distribution $\pi$ and is therefore well suited to quantify mode-collapse. This is further illustrated in Figure 1. We can again leverage the chain rule for the KL-divergence (Cover, 1999) to obtain an extended version of the EUBO, i.e., $\mathbb{E}_{\pi(\mathbf{x},\mathbf{z})}[\log\,\overline{\!{w}}]$ that satisfies $\,\overline{\!{\text{EUBO}}}\geq\text{EUBO}$, where the introduction of latent variables introduce additional looseness. The extended EUBO requires computing the importance weights $\,\overline{\!{w}}$ and expectations under $\pi(\mathbf{x},\mathbf{z})$. The former depends on the specific choice of sampling algorithm and is further discussed in Section 4 when introducing the methods considered for evaluation. The latter can be approximated by propagating target samples $\mathbf{x}$ back to $\mathbf{z}$ using $\pi(\mathbf{z}|\mathbf{x})$. Additionally, having access to samples from $\pi$ allows for computing forward versions of $Z$ and ESS which have already been used to quantify mode collapse by e.g. (Midgley et al., 2023). Formally, they are defined as

where expectations are taken with respect to the target $\pi$. For a detailed discussion see Appendix A.1.

Integral Probability Metrics. Alternatively, IPMs are often employed if samples from the target distribution $\pi$ are available (Arenz et al., 2018; Richter et al., 2023; Vargas et al., 2023a, 2024). Common IPMs for assessing sample quality are 2-Wasserstein distance ($\mathcal{W}_{2}$) (Peyré et al., 2019) or the maximum mean discrepancy (MMD) (Gretton et al., 2012). The former uses a cost function to calculate the minimum cost required to transport probability mass from one distribution to another while the latter assesses distribution dissimilarity by examining the differences in their mean embeddings within a reproducing kernel Hilbert space (Aronszajn, 1950). For further details see Appendix A.2.

Entropic Mode Coverage. Inspired by inception scores and distances from generative modelling (Salimans et al., 2016; Heusel et al., 2017) we propose a heuristic approach for detecting mode collapse by introducing the entropic mode coverage (EMC). To compute EMC, we partition $\mathbb{R}^{d}$ into disjoint subsets $\xi_{i},i\in\{1,\ldots,M\}$ that describe different modes of the target density $\pi$. Moreover, we introduce an auxiliary distribution that measures the probability of a sample $\mathbf{x}$ being element of a mode descriptor, i.e., $p(\xi_{i}|\mathbf{x})=p(\mathbf{x}\in\xi_{i})$. We then compute the expected entropy of the auxiliary distribution, that is,

where the expectation is approximated using a Monte Carlo estimate. Here, $N$ denotes the number of samples drawn from $q^{\theta}$. Please note that we employ the logarithm with a base of $M$ to ensure that $\text{EMC}\in[0,1]$. This choice of base facilitates a straightforward interpretation: A value of $0$ signifies a model that consistently produces samples that are elements of the same mode descriptor. In contrast, a value of $1$ represents a model that can produce samples from all mode descriptors with equal probability.

Clearly, EMC is limited to targets where mode descriptors are available which is further discussed in Section 5. Moreover, a suitable criterion for cases where mode descriptors are not equally probable is discussed in Appendix A.3.

In this section, we elaborate on the methods included in this benchmark, categorizing them into three distinct groups based on the computation of importance weights. Please refer to Table 1 for an overview of these methods and to Appendix B for further details.

Tractable Density Models. Tractable density models allow for computing the model likelihood $q^{\theta}(\mathbf{x})$. It is therefore straightforward to compute performance criteria associated with importance weights $w={\gamma(\mathbf{x})}/{q^{\theta}(\mathbf{x})}$. Notable works include factorized (‘mean-field’) Gaussian distributions (MFVI), Normalizing Flows (NFVI) (Rezende & Mohamed, 2015) and full rank Gaussian mixture models (GMMVI) (Arenz et al., 2022).

Sequential Importance Sampling Methods. Sequential importance sampling (SIS) methods define $\,\overline{\!{w}}$ in terms of incremental importance sampling (IS) weights, that is, $\,\overline{\!{w}}=\prod_{t=1}^{T}G_{t}(\mathbf{x}_{t-1},\mathbf{x}_{t})$ with

with annealed versions $\gamma_{t}$ of $\gamma$. For example, choosing $B^{\theta}_{t-1}(\mathbf{x}_{t-1}|\mathbf{x}_{{t}})=\pi_{t}(\mathbf{x}_{t-1})F^{\theta}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})/\pi_{t}(\mathbf{x}_{t})$ recovers AIS (Neal, 2001). Midgley et al. (2022) proposed to parameterize the proposal distribution $\pi_{0}$ with normalizing flows and, in combination with AIS, to minimize the $\alpha$-divergence, resulting in the Flow Annealed Importance Sampling Bootstrap (FAB) algorithm. Additionally, when AIS is coupled with resampling, it gives rise to Sequential Monte Carlo (SMC) as originally proposed by Del Moral et al. (2006).

Recent advancements include the development of Stochastic Normalizing Flows (Wu et al., 2020a), Annealed Flow Transport (AFT) (Arbel et al., 2021), and Continual Repeated AFT (CRAFT) (Matthews et al., 2022). These methods extend Sequential Monte Carlo by employing sets of normalizing flows that define deterministic transport maps between neighboring distributions $\gamma_{t}$. For further details on $F^{\theta}_{t},B^{\theta}_{t-1}$ and the corresponding $G_{t}$ see Table 7. For an in-depth exploration of the commonalities and distinctions among these methods, please refer to (Matthews et al., 2022).

Diffusion-Based Methods. Diffusion-based methods typically build on stochastic differential equations (SDEs) with parameterized drift terms (Tzen & Raginsky, 2019), i.e.,

with diffusion coefficient $\sigma_{t}$ and standard Brownian motion $\mathbf{w}_{t},\bar{\mathbf{w}}_{t}$. Using the Euler-Maruyama method (Särkkä & Solin, 2019), their discretized counterparts can be characterized by Gaussian forward-backward transition kernels

with discretization step size $\Delta_{t}$. The extended (unnormalized) importance weights $\,\overline{\!{w}}$ can then be constructed as

One line of work considers annealed Langevin dynamics to model Eq. (4). Works include Unadjusted Langevin Annealing (ULA) (Thin et al., 2021), Monte Carlo Diffusions (MCD) (Doucet et al., 2022a), Controlled Monte Carlo Diffusion (CMCD) (Vargas et al., 2024), Uncorrected Hamiltonian Annealing (UHA) (Geffner & Domke, 2021) and Langevin Diffusion Variational Inference (LDVI) (Geffner & Domke, 2022). A second line of work describes diffusion-based sampling from a stochastic optimal control perspective (Dai Pra, 1991). Works include methods such as Path Integral Sampler (PIS) (Zhang & Chen, 2021; Vargas et al., 2023b), Denoising Diffusion Sampler (DDS) (Vargas et al., 2023a), Time-Reversed Diffusion Sampler (DIS) (Berner et al., 2022) and Generative Flow Networks (GFN) (Lahlou et al., 2023; Malkin et al., 2022b; Zhang et al., 2023). Furthermore, Richter et al. (2023) identify several of these methods as special cases of a General Bridge Sampler (GBS) where both processes in Eq. 4 are freely parameterized. Specific choices of $\pi_{0},F^{\theta}_{t+1},B^{\theta}_{t-1}$ are detailed in Table 6. Lastly, we refer the interested reader to (Sendera et al., 2024) which concurrently benchmarked diffusion-based sampling methods.

Here, we briefly summarize the target densities $\pi$ considered in this work. The dimensionality of the problem, if we have access to the log normalizer $\log Z$, target samples, or mode descriptors for computing the entropic mode coverage is presented in Table 2. Further details and formal definitions of the target densities can be found in Appendix C.

Bayesian Logistic Regression. We consider four experiments where we perform inference over the parameters of a Bayesian logistic regression model for binary classification. The datasets Credit and Cancer were taken from Nishihara et al. (2014). The former distinguishes individuals as either good or bad credit risks, while the latter deals with the classification of recurrence events in breast cancer. The Ionosphere dataset (Sigillito et al., 1989) involves classifying radar signals passing through the ionosphere as either good or bad. Similarly, the Sonar dataset (Gorman & Sejnowski, 1988) tackles the classification of sonar signals bounced off a metal cylinder versus those bounced off a roughly cylindrical rock.

Random Effect Regression. The Seeds data was collected by (Crowder, 1978). The goal is to perform inference over the variables of a random effect regression model that models the germination proportion of seeds arranged in a factorial layout by seed and type of root.

Time Series Models. We consider the Brownian time series model obtained by discretizing a stochastic differential equation, modeling a Brownian motion with a Gaussian observation model, developed by (Sountsov et al., 2020).

Spatial Statistics. The log Gaussian Cox process (LGCP) (Møller et al., 1998) is a probabilistic model commonly used in statistics to model spatial point patterns. In this work, the log Gaussian Cox process is applied to modeling the positions of pine saplings in Finland.

Synthetic Targets. We additionally consider synthetic target densities as they commonly give access to the true normalization constant $Z$, target samples, and mode descriptors. The Funnel target was introduced by (Neal, 2003) and provides a complex ‘funnel’-shaped distribution. Moreover, we consider two different types of mixture models: a mixture of isotropic Gaussians (MoG) as proposed by Midgley et al. (2022), and Student-t distributions (MoS). To obtain mode descriptors for a mixture model with $K$ components, i.e., $\pi(\mathbf{x})=\sum_{k}\pi_{k}(\mathbf{x})/K$ we compute the density per component $\pi_{k}(\mathbf{x})$ and say that $\mathbf{x}\in\xi_{i}$ if $i=\text{argmax}_{k}\{\pi_{k}(\mathbf{x})\}_{k=1}^{K}$. Lastly, we follow Doucet et al. (2022a) and train NICE (Dinh et al., 2014) on a down-sampled $14\times 14$ variant of MNIST (Digits) (LeCun et al., 1998) and a $28\times 28$ variant of Fashion MNIST (Fashion) and use the trained model as target density. Here, we obtain the mode descriptors by training a classifier $p(\mathbf{c}|\mathbf{x})$ on samples from $\pi$ where the classes $\mathbf{c}$ are represented by ten different digits. If $i=\text{argmax}_{\mathbf{c}}\ p(\mathbf{c}|\mathbf{x})$ we conclude $\mathbf{x}\in\xi_{i}$.

In this section, we provide details on hyperparameter tuning. For further information, please refer to Appendix D.

Tractable Density Methods. For MFVI, we used a batch size of 2000 and performed 100k gradient steps, tuning the learning rate via grid search. For targets with known support, we adjusted the initial model variance accordingly. For GMMVI, we adhered to the default settings from (Arenz et al., 2022), utilizing 100 samples per mixture component. We initialized with 10 components and employed an adaptive scheme to add and remove components heuristically. The initial variance of the components was set based on the target support, and we conducted 3000 training iterations.

Sequential Importance Sampling Methods. In SIS methods, we employed 2000 particles for training. All methods except FAB used 128 annealing steps; FAB followed the original 12 steps as proposed by its authors. The choice and parameters of the MCMC transition kernel significantly impacted performance. Hamilton-Monte Carlo (Duane et al., 1987) generally outperformed Metropolis-Hastings (Chib & Greenberg, 1995) (see Appendix F.3). Step sizes for $\beta_{t}\geq 0.5$ and $\beta_{t}<0.5$ were tuned using grid search. For AFT and CRAFT, we used diagonal affine flows (Papamakarios et al., 2021), which yielded more robust results than complex flows like inverse autoregressive flows (Kingma et al., 2016) or neural spline flows (Durkan et al., 2019) (see Appendix F.6). FAB employed RealNVP (Dinh et al., 2016) for the proposal distribution $\pi_{0}$. Learning rates for these flows were also tuned via grid search. For targets with known support, the variance of $\pi_{0}(\mathbf{x})=\mathcal{N}({0},\sigma^{2}_{0}{I})$ was set accordingly, otherwise, a grid search was performed. We used multinomial resampling with a threshold of 0.3 (Douc & Cappé, 2005).

Diffusion-based Methods. Training involved a batch size of 2000 and 40k gradient steps. SDEs were discretized with 128 steps, $T=1$, and a fixed $\Delta t$. The diffusion coefficient was chosen as $\sigma_{t}=\sigma_{\text{max}}\cos^{2}({\pi(T-t)}/{2T})$, following (Vargas et al., 2023a) for better performance compared to linear or constant schedules. We used the architecture from (Zhang & Chen, 2021) with 2 layers of 64 hidden units each. For targets with known prior support, the initial model support was set accordingly. For all methods except PIS, this involved setting the variance of the prior distribution $\pi_{0}(\mathbf{x})=\mathcal{N}({0},\sigma^{2}_{0}{I})$. For PIS, $\sigma_{\text{max}}$ was carefully chosen. In MCD and LDVI, we learned the annealing schedule $\beta_{t}$ and $\sigma_{\text{max}}$ end-to-end by maximizing the ELBO.

Here, we offer an overview of the evaluation protocol. Next, we present the results obtained for synthetic target densities, followed by those for real targets. We provide further results in Appendix E and ablation studies in Appendix F.

Evaluation Protocol. We compute all performance criteria 100 times during training, applying a running average with a length of 5 over these evaluations to obtain robust results within a single run. To ensure robustness across runs, we use four different random seeds and average the best results from each run. We use 2000 samples to compute the performance criteria and tune key hyperparameters such as the learning rate and variance of the initial proposal distribution $\pi_{0}$. We report the EMC values corresponding to the method’s highest ELBO value to avoid high EMC values caused by a large initial support of the model.